You can always re-parameterize your problem such that you can fit an unconstrained derived parameter using Lev_Mar routines.
Some ideas for constraining are given for example in
this link. (Look for "One-sided boundary constraints").
In your case, you could fit for y= a1^2 + b1^2x + c1^2x^2 in which a1, b1, and c1 are unconstrained, then calculate a=a1^2 ,b=b1^2, and c=c1^2 at the end.