# Some problem when optimizing the ODEProblem parameters

**URL:** <https://discourse.julialang.org/t/some-problem-when-optimizing-the-odeproblem-parameters/113746>\
**Category:** General Usage\
**Tags:** optimization, modelingtoolkit\
**Created:** [May 2, 2024, 2:40pm UTC](https://discourse.julialang.org/t/some-problem-when-optimizing-the-odeproblem-parameters/113746 "2024-05-02T14:40:16Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![chooron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chooron/32/208300_2.png) [@chooron](https://discourse.julialang.org/u/chooron)\
**Post date:** [May 2, 2024, 2:40pm UTC](https://discourse.julialang.org/t/some-problem-when-optimizing-the-odeproblem-parameters/113746/1 "2024-05-02T14:40:16Z")

</div>

Use ModelingToolkit.jl to build ODEproblem and Optimization respectively to optimize the parameters of ODEproblem with the following code:

```julia
using ModelingToolkit, Optimization, ComponentArrays
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using OptimizationOptimisers

@mtkmodel FOL begin
    @parameters begin
        a
        b # parameters
    end
    @variables begin
        x(t) # dependent variables
    end
    @equations begin
        D(x) ~ (1 - x) / a + b
    end
end

@mtkbuild fol = FOL()
prob = ODEProblem(fol, [fol.x => 0.0], (1.0, 10.0), [fol.a => 3.0, fol.b => 4.0])
sol = OrdinaryDiffEq.solve(prob, Tsit5(), saveat=1.0)

@variables begin
    vs[1:2] = [3.0, 4.0] # , [bounds = (-2.0, 2.0)]
end
x_axes = getaxes(ComponentVector(a=1, b=2))
@parameters st=0.0

v1 = 1:10
target = @.((v1 - (v1^2)/2)/3.0+4.0*v1)
loss = begin
    params = ComponentVector(vs, x_axes)
    new_prob = remake(prob, u0=[fol.x => st], p=[fol.a => params[:a], fol.b => params[:b]])
    sol = OrdinaryDiffEq.solve(new_prob, Tsit5(), saveat=1.0)
    sol_u = [v[1] for v in sol.u]
    sum(abs.(target .- sol_u))
end

@mtkbuild sys = OptimizationSystem(loss, vs, [st]) # , constraints=cons
u0 = [vs[1] => 3.0
    vs[2] => 4.0]
prob = OptimizationProblem(sys,
    u0,
    grad=true,
    hess=true,
    cons_j=true,
    cons_h=true)
solve(prob, OptimizationOptimisers.Adam(), maxiters=10)

```

Then the code gets this error

```julia
ERROR: MethodError: promote_rule(::Type{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Num}, Num, 1}}, ::Type{Num}) is ambiguous.

Candidates:
  promote_rule(::Type{ForwardDiff.Dual{T, V, N}}, ::Type{R}) where {T, V, N, R<:Real}
    @ ForwardDiff D:\software\Julia-1.10.0\pkg\packages\ForwardDiff\PcZ48\src\dual.jl:428
  promote_rule(::Type{<:Number}, ::Type{<:Num})
    @ Symbolics D:\software\Julia-1.10.0\pkg\packages\Symbolics\Eas9m\src\num.jl:97

Possible fix, define
  promote_rule(::Type{ForwardDiff.Dual{T, V, N}}, ::Type{R}) where {T, V, N, R<:Num}

Stacktrace:
 [1] promote_type(::Type{Num}, ::Type{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Num}, Num, 1}})
   @ Base .\promotion.jl:313
 [2] promote_eltype(::Type{Vector{Num}}, ::Type{ForwardDiff.Dual{ForwardDiff.Tag{…}, Num, 1}})
   @ ArrayInterface D:\software\Julia-1.10.0\pkg\packages\ArrayInterface\NIzGA\src\ArrayInterface.jl:114
 [3] wrapfun_iip(ff::Function, inputs::Tuple{Vector{…}, Vector{…}, ModelingToolkit.MTKParameters{…}, Float64})
   @ DiffEqBase D:\software\Julia-1.10.0\pkg\packages\DiffEqBase\NaUtB\src\norecompile.jl:59
 [4] promote_f(f::ODEFunction{…}, ::Val{…}, u0::Vector{…}, p::ModelingToolkit.MTKParameters{…}, t::Float64)
   @ DiffEqBase D:\software\Julia-1.10.0\pkg\packages\DiffEqBase\NaUtB\src\solve.jl:1260
 [5] get_concrete_problem(prob::ODEProblem{…}, isadapt::Bool; kwargs::@Kwargs{…})
   @ DiffEqBase D:\software\Julia-1.10.0\pkg\packages\DiffEqBase\NaUtB\src\solve.jl:1173
 [6] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::ModelingToolkit.MTKParameters{…}, args::Tsit5{…}; kwargs::@Kwargs{…})
   @ DiffEqBase D:\software\Julia-1.10.0\pkg\packages\DiffEqBase\NaUtB\src\solve.jl:1074
 [7] solve(prob::ODEProblem{…}, args::Tsit5{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
   @ DiffEqBase D:\software\Julia-1.10.0\pkg\packages\DiffEqBase\NaUtB\src\solve.jl:1003
 [8] top-level scope
   @ f:\julia\DeepFlex.jl\test\pkgs\test_tmk_opt.jl:35
Some type information was truncated. Use `show(err)` to see complete types.

```

In addition varmap\_to\_vars to construct new parameters (the rest remains the same), see below:

```julia
loss = begin
    params = ComponentVector(vs, x_axes)
    # unew = [fol.x => st]
    # pnew = [fol.a => params[:a], fol.b => params[:b]]
    unew = ModelingToolkit.varmap_to_vars([fol.x => st], unknowns(fol))
    pnew = ModelingToolkit.MTKParameters(sys, [fol.a => params[:a], fol.b => params[:b]], unew)
    new_prob = remake(prob, u0=unew, p=pnew)
    sol = OrdinaryDiffEq.solve(new_prob, Tsit5(), saveat=1.0)
    sol_u = [v[1] for v in sol.u]
    sum(abs.(target .- sol_u))
end

```

However, the following problems also arise:

```julia
ERROR: MethodError: no method matching AbstractFloat(::Type{SymbolicUtils.BasicSymbolic{Real}})

Closest candidates are:
  AbstractFloat(::Bool)
   @ Base float.jl:264
  AbstractFloat(::UInt32)
   @ Base float.jl:272
  AbstractFloat(::UInt16)
   @ Base float.jl:271
  ...

Stacktrace:
 [1] float(x::Type)
   @ Base .\float.jl:294
 [2] promote_to_concrete(vs::Vector{SymbolicUtils.BasicSymbolic{Real}}; tofloat::Bool, use_union::Bool)
   @ ModelingToolkit D:\software\Julia-1.10.0\pkg\packages\ModelingToolkit\kByuD\src\utils.jl:634
 [3] varmap_to_vars(varmap::Vector{…}, varlist::Vector{…}; defaults::Dict{…}, check::Bool, toterm::Function, promotetoconcrete::Nothing, tofloat::Bool, use_union::Bool)
   @ ModelingToolkit D:\software\Julia-1.10.0\pkg\packages\ModelingToolkit\kByuD\src\variables.jl:163
 [4] varmap_to_vars(varmap::Vector{Pair{Num, Num}}, varlist::Vector{SymbolicUtils.BasicSymbolic{Real}})
   @ ModelingToolkit D:\software\Julia-1.10.0\pkg\packages\ModelingToolkit\kByuD\src\variables.jl:125
 [5] top-level scope
   @ f:\julia\DeepFlex.jl\test\pkgs\test_tmk_opt.jl:37
Some type information was truncated. Use `show(err)` to see complete types.

```

Hopefully someone can solve my problem

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [May 2, 2024, 4:49pm UTC](https://discourse.julialang.org/t/some-problem-when-optimizing-the-odeproblem-parameters/113746/2 "2024-05-02T16:49:10Z")

</div>

Here is a version of your code that runs:

> **Code**
>
> ```julia
> using ModelingToolkit, Optimization, ComponentArrays
> using ModelingToolkit: t_nounits as t, D_nounits as D
> using OrdinaryDiffEq
> using OptimizationOptimisers
> 
> @mtkmodel FOL begin
> @parameters begin
> a
> b # parameters
> end
> @variables begin
> x(t) # dependent variables
> end
> @equations begin
> D(x) ~ (1 - x) / a + b
> end
> end
> 
> @mtkbuild fol = FOL()
> prob = ODEProblem(fol, [fol.x => 0.0], (1.0, 10.0), [fol.a => 3.0, fol.b => 4.0])
> sol = OrdinaryDiffEq.solve(prob, Tsit5(), saveat=1.0)
> 
> begin
> vs = [3.0, 4.0] # , [bounds = (-2.0, 2.0)]
> end
> x_axes = getaxes(ComponentVector(a=1, b=2))
> st=0.0
> 
> v1 = 1:10
> target = @.((v1 - (v1^2)/2)/3.0+4.0*v1)
> function loss(vs, st)
> params = ComponentVector(vs, x_axes)
> new_prob = remake(prob, u0=[fol.x => st], p=[fol.a => params[:a], fol.b => params[:b]])
> sol = OrdinaryDiffEq.solve(new_prob, Tsit5(), saveat=1.0)
> sol_u = sol[1,:]
> return sum(abs.(target .- sol_u))
> end
> 
> u0 = [3.0, 4.0]
> optfun = OptimizationFunction(loss, Optimization.AutoForwardDiff())
> optprob = OptimizationProblem(optfun,
> u0, 
> st)
> solve(optprob, OptimizationOptimisers.Adam(), maxiters=10)
> 
> ```

The changes that were needed are:

- `vs` and `st` should be `Float64`s and not abstract parameters.
- the variable `prob` was first the ODE and then the optimisation problem (but still used in `remake`).
- The loss function is too complicated to be expressed as a symbolic term, so, you cannot carry abstract parameters all the way through it. That is still fine, as you can use automatic differentiation to get the derivatives, but it doesn’t work as you planned to do it.

Hope that helps 🙂

---

<div class="post-metadata">

**Author:** ![chooron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chooron/32/208300_2.png) [@chooron](https://discourse.julialang.org/u/chooron)\
**Post date:** [May 3, 2024, 3:43am UTC](https://discourse.julialang.org/t/some-problem-when-optimizing-the-odeproblem-parameters/113746/3 "2024-05-03T03:43:16Z")

</div>

Thanks for your help, the code you gave does work and get the calculation results, but I prefer to be able to use ModelingToolkit to automatically calculate the gradient of the problem when building ODEProblem, I have tried a similar calculation method to your code, but there are some problems in my project, so I rebuilt the above code to test its feasibility. But from your suggestion, I don’t feel that this way is feasible either, so I will continue to consider other options. Finally, thank you very much for your advice.

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [May 3, 2024, 3:59am UTC](https://discourse.julialang.org/t/some-problem-when-optimizing-the-odeproblem-parameters/113746/4 "2024-05-03T03:59:07Z")

</div>

Yeah, I think there was work on deriving the sensitivity equations on a symbolic level, but as far as I know it is not done yet. See

> [@Direct Forward Sensitivity Analysis of ODEs in ModelingToolkit.jl for a Subset of Parameters](https://discourse.julialang.org/t/direct-forward-sensitivity-analysis-of-odes-in-modelingtoolkit-jl-for-a-subset-of-parameters/104903):
>
> Hello, I am working on a project that involves setting up and solving an ODE system. I’ve defined my ODE system using ModelingToolkit.jl and not as a function. My ODE system has many parameters, and I’m interested in performing a Direct Forward Sensitivity Analysis on a subset of these parameters, not all of them. Post-analysis, I intend to apply Bayesian inference to the results further. Here’s a simplified overview of what I am trying to do: using ModelingToolkit, Symbolics, DifferentialEqu…

You might already tried it, but you can get the adjoint equations with this package:

> **[Direct Adjoint Sensitivities of Differential Equations · SciMLSensitivity.jl](https://docs.sciml.ai/SciMLSensitivity/dev/manual/direct_adjoint_sensitivities/)**
>
> Documentation for SciMLSensitivity.jl.

(The adjoint ODE computes the derivative of a loss with respect to parameters. It’s mathematically kind of the ‘symbolic derivative’, which you can get to compute this gradient.)

---

<div class="post-metadata">

**Author:** ![chooron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chooron/32/208300_2.png) [@chooron](https://discourse.julialang.org/u/chooron)\
**Post date:** [May 3, 2024, 4:06am UTC](https://discourse.julialang.org/t/some-problem-when-optimizing-the-odeproblem-parameters/113746/5 "2024-05-03T04:06:03Z")

</div>

I haven’t tried this method yet, I will give it a try, thanks for the advice 😀