Some problem when optimizing the ODEProblem parameters

[chooron]

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

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

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:

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:

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

[SteffenPL]

Here is a version of your code that runs:

Code

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 Float64s 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

[chooron]

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.

[SteffenPL]

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

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

(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.)

[chooron]

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