Typically, you want to at least be able to express the lengths of lines, and the values of angles. I don't know how you could do this without using at least quadratic terms. So you have to do at least quadratic programming, and there's some things that can add higher degree terms.

Moreover, extremely similar problems are studied in automated geometric theorem proving, and they also get pulled into using polynomials.

Unfortunately, these problems are much much harder than linear programming.

There's some really powerful tools from algebraic geometry, mainly Grobner Bases, which you would ideally want to use. They're ridiculously powerful. Not only can you use them to solve systems of polynomial equations, but you can also test whether a set of equations implies another polynomial, or conflicts with another polynomial -- super useful for dealing with over-constrained models.

Tragically, the relevant algorithm is O(exp(exp(n))) -- well, or at least, that's the best upper bound anyone can probe, in practice it isn't that bad. Still, I quickly ran into computational problems. Grobner bases, also, while great are tricky to work with if you only want real valued solutions. They like to give you complex valued solutions to problems, and complex valued solutions to your CAD constraints are usually not what you are looking for.

> How does gradient descent work with over-constrainted models ? (as in: impossible to satisfy all constraints). You need some way to back out of those or relax one of the constraints.

In my implementation, each constraint was a separate cost function, which I added together. Normally, when you run the solver, it pushes all the constraint costs to zero, and you know you've solved everything.

Sometimes it would get stuck in a local minimum and require user intervention to help it get to the right solution. That was surprisingly rare.

If you have conflicting constraints it simply fails to reach a cost of zero at all. You can handle that in your interface however you like. The obvious thing is to try and solve a constraint immediately after adding it and then reject the constraint if you can't.

(I really didn't expect a local optimization method like gradient descent to work so well, but it did. That said, if I attacked the problem again I would try to use some sort of second order optimization method, I think.)

My my view is probably way too simplistic, but I'm really wondering how you would end up with quadratic terms if you need to define stuff like:

* this line needs to be 1/3 the length of that line

* these set of lines should start 50mm left of this line

I can see how trigonometric expressions could be useful to specify angles (although a constant there covers 99% of the user-cases?).

Remember that you're working with points. So, to express the length of a line, you need something like (x1-x2)^2 + (y1-y2)^2 .

If you just wanted to do things like "this point is +(0,10) of this other one", you're doing a really simple constructive system, which is easy to evaluate, but not anywhere near as flexible or powerful. ImplicitCAD/OpenSCAD trivially do that, just without a nice interface.