Optimization on a manifold

I am 100% with you. Replication is essential.

But the other thing I have found is that by using more robust optimizers it is often easier to get a sense of the domains of attraction from initial conditions (and, frequently, problems in my models which need to be fixed or corners that I didn’t anticipate would ever bind).

In particular, you can tell these optimizers to emphasize feasibility over optimizality and then it can carefully walk along the F(a,p) = 0 manifold which might be essential for the stability of your other code. Another thing you can do is fix an initial a and tell it to find a feasible point without worrying about the objective. Once you know “the solution” you can backwards engineer conditions to guide the less robust optimizers to get to it and drop the commercial one.

But with all that, I completely understand your aversion to experimenting with commercial optimizers. And I suspect that black box derivative-free ones are never going to work with a problem of this type because feasibility on the constraint is too important.

Try setting some options in Ipopt, e.g. start_with_resto = true or max_soc = 50. You can also use second order Ipopt or AugLag by using the option first_order = false. The default in Nonconvex is first order.

Note that you can differentiate M(a,p(a)) with respect to a very efficiently and easily using an adjoint method, a.k.a. implicit functions — given the solution of F(a,p) = 0 and the Jacobian \partial F/\partial p, you only have to solve one linear system with \partial M / \partial p as a right-hand side. See e.g. section 3 of my course notes.

You could, alternatively, just use an optimization algorithm like SLSQP that supports nonlinear equality constraints directly, and then just solve \min_{a,p} M(a,p) subject to F(a,p)=0.

(For example, a better version of your “penalty” method is the augmented Lagrangian method, which adds both a penalty and a Lagrangian term and is designed so that the penalty does not need to become extremely large in order to enforce the constraint to a small tolerance — unlike your simple L2 penalty.)

1 Like

Thanks. In fact I tried using a constraint (via SLSQP and the augmented Lagrangian), see here:

but I keep getting FORCED_STOP from NLopt and I am not quite sure how to debug that.

Is your code throwing an exception?

I don’t think so. My reading of the manual suggests that NLopt would rethrow those and I would notice, is this correct?