Nonlinear optimization with ODE solution in nonlinear constraint

[afossa]

Hello,

I am trying to model and solve a fairly complex constrained nonlinear optimization problem. The main issue is that to compute some equality constraints I need to solve an ODE numerically.

So far, I have been able to model a “minimal” version of the problem using Optimization and OrdinaryDiffEq, and solve it with the IPNewton() method implemented in Optim.

However, I want to swap this solver with a more performant one, i.e. Ipopt, which better handles large problems. So far I tried the followings solutions, but unsuccessfully:

1. Use Ipopt through OptimizationMOI and AmplNLWriter: issues with Jacobian/Hessian computation since Optimization.AutoForwardDiff(), which works with IPNewton(), does not provide the “expression graph” required by AmplNLWriter. Switching to Optimization.AutoModelingToolkit() does not work on my model. Providing explicit expressions for the Jacobian/Hessian, other than being excessively cumbersome, does not solve the “expression graph” problem either.
2. Switch to JuMP to model my problem: my variables and constraints are vector-valued, but JuMP lacks streamlined support for these. Plus, I tried to add the nonlinear constraint that solves the ODE in different ways, by registering an @operator and providing “explicit” Jacobian/Hessian obtained with ForwardDiff, but without success.

So my main questions are:

1. Which framework (Optimization, JuMP, ModelingToolkit, …) would you suggest to model my problem?
2. Would it be feasible to model my NLP with ModelingToolkit if also the ODE is modeled with the same toolkit? If yes, how?
3. Is there a way to interface Optimization to Ipopt directly, avoiding the AmplNLWriter that is causing issues?

Finally, my NLP is very large, and I have several constraints that solve the same ODE with different initial conditions. Thus, having a model that can automatically detect and take advantage of the sparsity structure of the derivatives will be highly beneficial.

I attach a small (working) example that should be representative of my problem. It is a Hohmann transfer which can be solved analytically, but it has the same kind of constraints I have in my full NLP.

optim_ex.jl (1.5 KB)

P.S.: I use Julia 1.10.4 and the latest version of each package.

[gdalle]

I don’t have a full answer to your problem but here are some remarks.

@afossa

Providing explicit expressions for the Jacobian/Hessian, other than being excessively cumbersome, does not solve the “expression graph” problem either.

You definitely don’t want to provide full Jacobians and Hessians, that’s what autodiff is here for.

@afossa

Switch to JuMP to model my problem: my variables and constraints are vector-valued, but JuMP lacks streamlined support for these.

That doesn’t seem true to me, JuMP lets you easily specify vector variables and constraints. What have you tried that failed.

@afossa

Plus, I tried to add the nonlinear constraint that solves the ODE in different ways, by registering an @operator and providing “explicit” Jacobian/Hessian obtained with ForwardDiff, but without success.

Here too, it would be nice to have an idea of what you tried unsuccessfully, ideally an MWE.

@afossa

Finally, my NLP is very large, and I have several constraints that solve the same ODE with different initial conditions. Thus, having a model that can automatically detect and take advantage of the sparsity structure of the derivatives will be highly beneficial.

You can check out SparseConnectivityTracer.jl and DifferentiationInterface.jl (especially the sparsity handling tutorial) to compute sparse exact Jacobians and Hessians with very little work on your end. Ping @hill, the co-developer of these packages.

[afossa]

Thank you @gdalle

The best “almost working” example I got with JuMP is:

using JuMP, Ipopt

mag(Δv) = sqrt(Δv[1]^2+Δv[2]^2);
mdl = Model(Ipopt.Optimizer);

@variable(mdl, u[1:nv]);
@objective(mdl, Min, mag(u[9:10]) + mag(u[11:12]));
@constraint(mdl, x10+B*u[9:10] .== u[1:4]);
@constraint(mdl, u[5:8]+B*u[11:12] .== x2f);

@operator(mdl, op_dyna_1, 4, (u...) -> dyna(collect(u))[1]);
@operator(mdl, op_dyna_2, 4, (u...) -> dyna(collect(u))[2]);
@operator(mdl, op_dyna_3, 4, (u...) -> dyna(collect(u))[3]);
@operator(mdl, op_dyna_4, 4, (u...) -> dyna(collect(u))[4]);

@constraint(mdl, op_dyna_1(u[1:4]...) == u[5]);
@constraint(mdl, op_dyna_2(u[1:4]...) == u[6]);
@constraint(mdl, op_dyna_3(u[1:4]...) == u[7]);
@constraint(mdl, op_dyna_4(u[1:4]...) == u[8]);

optimize!(mdl);

You can run it by prepending the code attached to the original post. The ODE solution throws a lot of warnings about the step size being almost zero, and it finally fails.

For the explicit derivatives of operators, I tried this JuMP example.
I also checked out this trick, but adding scalar operators and constraints one by one is not an option for me. And I do not see any support for operators with vector outputs.

[ChrisRackauckas]

@gdalle

@afossa

Finally, my NLP is very large, and I have several constraints that solve the same ODE with different initial conditions. Thus, having a model that can automatically detect and take advantage of the sparsity structure of the derivatives will be highly beneficial.

You can check out SparseConnectivityTracer.jl and DifferentiationInterface.jl (especially the sparsity handling tutorial) to compute sparse exact Jacobians and Hessians with very little work on your end. Ping @hill, the co-developer of these packages.

That’s not the way. You will miss a lot of SIMD optimizations which are also possible. Treating it purely as a sparse matrix misses the repeated structure optimizations. For a full treatment of that topic, especially in the context of GPU usage, see:

https://www.sciencedirect.com/science/article/abs/pii/S0045782523007156

And yes if you have a lot of initial conditions I would suggest GPUing it, though even without the GPU codegening to the specific form could be like 5x even on CPU.

@afossa

Which framework (Optimization, JuMP, ModelingToolkit, …) would you suggest to model my problem?

Definitely Optimization.jl. You’ll want the flexibility in AD choice here, since likely you want to not forward mode. Also, I’ve never been successful getting an ODE into JuMP’s nonlinear constraints, so it may just be a no-go.

@afossa

• Would it be feasible to model my NLP with ModelingToolkit if also the ODE is modeled with the same toolkit? If yes, how?

Currently you won’t see an advantage in this if you’ve already written the ODE in a performant way.

@afossa

Is there a way to interface Optimization to Ipopt directly

@Vaibhavdixit02 might have an answer.

[Vaibhavdixit02]

Yes you should use Ipopt directly, you can take a look at the tutorial here which uses it Using Equality and Inequality Constraints · Optimization.jl

The main difference is passing Ipopt.Optimizer() directly as the solver choice and not wrapping it in AmplNLWriter.Optimizer(Ipopt_jll.amplexe)

[afossa]

@Vaibhavdixit02

Yes you should use Ipopt directly, you can take a look at the tutorial here which uses it Using Equality and Inequality Constraints · Optimization.jl

The main difference is passing Ipopt.Optimizer() directly as the solver choice and not wrapping it in AmplNLWriter.Optimizer(Ipopt_jll.amplexe)

Thanks! I completely missed this example. It works better indeed.

[afossa]

@ChrisRackauckas

You’ll want the flexibility in AD choice here, since likely you want to not forward mode.

Which AD backend would you suggest?

So far, only ForwardDiff seems to work. The other choices listed in the docs fail with quite cryptic errors.
For example, AutoModelingToolkit throws a no method matching Float64(::Num) error during the initialization of the step size in the ODE solution, which I think is related to my previous question:

@afossa

Would it be feasible to model my NLP with ModelingToolkit if also the ODE is modeled with the same toolkit? If yes, how?

I can provide more detailed error stack traces if you want.

I didn’t have the time yet, but I’ll definitely look into this:

@ChrisRackauckas

For a full treatment of that topic, especially in the context of GPU usage, see:

https://www.sciencedirect.com/science/article/abs/pii/S0045782523007156

And yes if you have a lot of initial conditions I would suggest GPUing it, though even without the GPU codegening to the specific form could be like 5x even on CPU.