# Nonlinear optimization with ODE solution in nonlinear constraint

**URL:** <https://discourse.julialang.org/t/nonlinear-optimization-with-ode-solution-in-nonlinear-constraint/116730>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, optimization, nonlinear-optimizati\
**Created:** [July 7, 2024, 4:56pm UTC](https://discourse.julialang.org/t/nonlinear-optimization-with-ode-solution-in-nonlinear-constraint/116730 "2024-07-07T16:56:49Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![afossa](https://avatars.discourse-cdn.com/v4/letter/a/eb9ed0/32.png) [@afossa](https://discourse.julialang.org/u/afossa)\
**Post date:** [July 7, 2024, 4:56pm UTC](https://discourse.julialang.org/t/nonlinear-optimization-with-ode-solution-in-nonlinear-constraint/116730/1 "2024-07-07T16:56:49Z")

</div>

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](https://en.wikipedia.org/wiki/Hohmann_transfer_orbit) which can be solved analytically, but it has the same kind of constraints I have in my full NLP.

[optim\_ex.jl](https://discourse.julialang.org/uploads/short-url/siblhIl3QBOwRQji0fgPIokvR1s.jl) (1.5 KB)

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

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [July 7, 2024, 5:34pm UTC](https://discourse.julialang.org/t/nonlinear-optimization-with-ode-solution-in-nonlinear-constraint/116730/2 "2024-07-07T17:34:31Z")

</div>

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](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/dev/tutorial2/)) to compute sparse exact Jacobians and Hessians with very little work on your end. Ping @hill, the co-developer of these packages.

---

<div class="post-metadata">

**Author:** ![afossa](https://avatars.discourse-cdn.com/v4/letter/a/eb9ed0/32.png) [@afossa](https://discourse.julialang.org/u/afossa)\
**Post date:** [July 7, 2024, 6:04pm UTC](https://discourse.julialang.org/t/nonlinear-optimization-with-ode-solution-in-nonlinear-constraint/116730/3 "2024-07-07T18:04:08Z")

</div>

Thank you @gdalle

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

```julia
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](https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/operator_ad/#ForwardDiff).  
I also checked out [this trick](https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/tips_and_tricks/#User-defined-operators-with-vector-outputs), 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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 7, 2024, 11:34pm UTC](https://discourse.julialang.org/t/nonlinear-optimization-with-ode-solution-in-nonlinear-constraint/116730/4 "2024-07-07T23:34:13Z")

</div>

> [@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](https://juliahub.com/ui/Packages/General/SparseConnectivityTracer) and [DifferentiationInterface.jl](https://juliahub.com/ui/Packages/General/DifferentiationInterface) (especially the [sparsity handling tutorial](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/dev/tutorial2/)) 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](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.

---

<div class="post-metadata">

**Author:** ![Vaibhavdixit02](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhavdixit02/32/2916_2.png) [@Vaibhavdixit02](https://discourse.julialang.org/u/Vaibhavdixit02)\
**Post date:** [July 8, 2024, 12:06am UTC](https://discourse.julialang.org/t/nonlinear-optimization-with-ode-solution-in-nonlinear-constraint/116730/5 "2024-07-08T00:06:42Z")

</div>

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](https://docs.sciml.ai/Optimization/stable/tutorials/constraints/)

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

---

<div class="post-metadata">

**Author:** ![afossa](https://avatars.discourse-cdn.com/v4/letter/a/eb9ed0/32.png) [@afossa](https://discourse.julialang.org/u/afossa)\
**Post date:** [July 8, 2024, 5:34pm UTC](https://discourse.julialang.org/t/nonlinear-optimization-with-ode-solution-in-nonlinear-constraint/116730/6 "2024-07-08T17:34:48Z")

</div>

> [@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](https://docs.sciml.ai/Optimization/stable/tutorials/constraints/)
> 
> 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.

---

<div class="post-metadata">

**Author:** ![afossa](https://avatars.discourse-cdn.com/v4/letter/a/eb9ed0/32.png) [@afossa](https://discourse.julialang.org/u/afossa)\
**Post date:** [July 8, 2024, 5:51pm UTC](https://discourse.julialang.org/t/nonlinear-optimization-with-ode-solution-in-nonlinear-constraint/116730/7 "2024-07-08T17:51:54Z")

</div>

> [@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](https://docs.sciml.ai/Optimization/stable/API/ad/) 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](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.