# JuMP optimizing parametric DifferentialEquation

**URL:** <https://discourse.julialang.org/t/jump-optimizing-parametric-differentialequation/70302>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, forwarddiff, differentialequation\
**Created:** [October 24, 2021, 7:46pm UTC](https://discourse.julialang.org/t/jump-optimizing-parametric-differentialequation/70302 "2021-10-24T19:46:41Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![cricci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cricci/32/25160_2.png) [@cricci](https://discourse.julialang.org/u/cricci)\
**Post date:** [October 24, 2021, 7:46pm UTC](https://discourse.julialang.org/t/jump-optimizing-parametric-differentialequation/70302/1 "2021-10-24T19:46:41Z")

</div>

Hello,  
I am trying to do the following: I have to optimize some parameters of an ODE, using JuMP (I have an objective function that depends on the solution ode the ODE). Now the problem is that the ODE itself has some other parameters (that I don’t want to optimize) cooked into a struct (using Parameters.jl).  
So my idea has been that of defining an objective function for JuMP via `JuMP.register` that takes the parameter that I need to optimize and the constant one and bash them into the real objective function (joining them together as a tuple).

What happens is that Julia complains that I am trying to convert Dual numbers into Float64. I guess that is an issue when I extract the parameter that I need to optimize and the constant one (stored in a struct) inside the derivative in the ODE.

I have a minimal (non-)working example which is basically the adaptation of [Using JuMP with DiffEqParamEstim · DiffEqParamEstim.jl](https://diffeqparamestim.sciml.ai/dev/tutorials/jump/)  
by taking an ODE with additional parameter as in my case (also usign Ipopt because I like it more than NLopt).

Thanks in advance for the help!

```julia
using DifferentialEquations, JuMP, Ipopt, Parameters

@with_kw struct parameters
    extra = 1.0
end
param = parameters()

function g(du,u,p_join,t)
    p,param = p_join
    σ,ρ,β = p
    x,y,z = u
    du[1] = dx = σ*(y-x)*param.extra
    du[2] = dy = x*(ρ-z) - y
    du[3] = dz = x*y - β*z

end

u0 = [1.0;0.0;0.0]
t = 0.0:0.01:1.0
tspan = (0.0,1.0)
model_ode(p_) = ODEProblem(g, u0, tspan,p_)
solve_model(mp_) = DifferentialEquations.solve(model_ode(mp_), Tsit5(),saveat=0.01)

#test that the ODE works
solve_model(([8,15,10/3],param)) 

function loss_objective(mp_)
    mp_join_ = (mp_,param)
    sol = solve_model(mp_join_)
    return sol[1,end]
end
juobj(args...) = loss_objective(args)

jumodel = Model(Ipopt.Optimizer)

JuMP.register(jumodel, :juobj, 3, juobj, autodiff=true)
JuMP.@variables jumodel begin
    σ,(start=8)
    ρ,(start=25.0)
    β,(start=10/3)
end
@NLobjective(jumodel, Min, juobj(σ, ρ, β))
JuMP.unset_silent(jumodel)

JuMP.optimize!(jumodel)

```

output (up to relatively helpful lines):

```julia
julia> JuMP.optimize!(jumodel)
This is Ipopt version 3.13.4, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 0

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#138#140"{typeof(juobj)}, Float64}, Float64, 3})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:760
  (::Type{T})(::VectorizationBase.Double{T}) where T<:Union{Float16, Float32, Float64, VectorizationBase.Vec{var"#s33", var"#s32"} where {var"#s33", var"#s32"<:Union{Float16, Float32, Float64}}, VectorizationBase.VecUnroll{var"#s31", var"#s30", var"#s29", V} where {var"#s31", var"#s30", var"#s29"<:Union{Float16, Float32, Float64}, V<:Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8, SIMDTypes.Bit, VectorizationBase.AbstractSIMD{var"#s30", var"#s29"}}}} at /home/cricci/.julia/packages/VectorizationBase/TQLhj/src/special/double.jl:84
  ...
Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#138#140"{typeof(juobj)}, Float64}, Float64, 3})
    @ Base ./number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#138#140"{typeof(juobj)}, Float64}, Float64, 3}, i1::Int64)
    @ Base ./array.jl:839
  [3] g(du::Vector{Float64}, u::Vector{Float64}, p_join::Tuple{Tuple{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#138#140"{typeof(juobj)}, Float64}, Float64, 3}, ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#138#140"{typeof(juobj)}, Float64}, Float64, 3}, ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#138#140"{typeof(juobj)}, Float64}, Float64, 3}}, parameters}, t::Float64)
    @ Main ~/actualpath/JumpDiffEq.jl:7

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 24, 2021, 9:21pm UTC](https://discourse.julialang.org/t/jump-optimizing-parametric-differentialequation/70302/2 "2021-10-24T21:21:53Z")

</div>

JuMP is not a good choice for this sort of problem. See [https://jump.dev/JuMP.jl/stable/background/should\_i\_use/#Black-box,-derivative-free,-or-unconstrained-optimization](https://jump.dev/JuMP.jl/stable/background/should_i_use/#Black-box,-derivative-free,-or-unconstrained-optimization).

Since this is unconstrained, you’re better off using Optim.jl or GalacticOptim.jl for this.

Your error arises because your ODE solve needs to support generic number types. I assume the DifferentialEquations documentation has some advice on this somewhere.

---

<div class="post-metadata">

**Author:** ![cricci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cricci/32/25160_2.png) [@cricci](https://discourse.julialang.org/u/cricci)\
**Post date:** [October 24, 2021, 9:28pm UTC](https://discourse.julialang.org/t/jump-optimizing-parametric-differentialequation/70302/3 "2021-10-24T21:28:40Z")

</div>

Thanks!  
Sorry, this is by bad. My example here is unconstrained optimization, while the problem I actually have to solve is constrained (I just wanted to provide a short working example)

> Your error arises because your ODE solve needs to support generic number types. I assume the DifferentialEquations documentation has some advice on this somewhere.

I completely agree. It seems to me that while without the extra parameters DifferentialEquations is able to detect the type of `u` or `du` when `p` is a Dual number. While when I add the extra parameter then It doesn’t and just assumes Float64 that doesn’t work then.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 25, 2021, 12:17am UTC](https://discourse.julialang.org/t/jump-optimizing-parametric-differentialequation/70302/4 "2021-10-25T00:17:28Z")

</div>

I assume you’ll need something like

```julia
function model_ode(p_::AbstractVector{T}) where {T}
    u = [one(T), zero(T), zero(T)]
    return ODEProblem(g, u, tspan, p_)
end

```

Passing in `u0` which is a `Vector{Float64}` won’t work.

---

<div class="post-metadata">

**Author:** ![cricci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cricci/32/25160_2.png) [@cricci](https://discourse.julialang.org/u/cricci)\
**Post date:** [October 25, 2021, 7:32am UTC](https://discourse.julialang.org/t/jump-optimizing-parametric-differentialequation/70302/5 "2021-10-25T07:32:30Z")

</div>

Wow, how can I have been so dumb not to think of It. Thank you very much sir for you help.  
I will post here the pieces of code that need to be changed, in the spirit of what you suggested, in order to make it work.

```julia
function model_ode(p_::Tuple{AbstractVector{T}, parameters}) where {T}
    u0 = [one(T), zero(T), zero(T)]
    return ODEProblem(g, u0, tspan, p_)
end

function loss_objective(mp_)
    mp_join_ = (collect(mp_),param)
    sol = solve_model(mp_join_)
    return sol[1,end]
end

```