Parameter estimation fails for Modelingtoolkit model specification but works fine for DifferentialEquations.jl specification

General Usage
, , , ,
4

Since you’re using the vector in the optimization instead of a vector of pairs, it might. So even from the same initial conditions, the simulation might be different since the ordering is messed up ( might be that p1 and p4 are switched. ). If you find the right permutation of the vector, I assume that this helps. Otherwise you can build the function yourself using Symbolics.jl, which could help.

I ran you’re code. I cannot use the objective function, but maybe this helps a little:

using ModelingToolkit
using OrdinaryDiffEq
using LeastSquaresOptim
using RecursiveArrayTools
using Random

Random.seed!(1234)

@parameters t, p1, p2, p3, p4
D= Differential(t)
@variables  x1(t), x2(t)

eqs = [
    D(x1) ~  p1*x1 - p2*x1*x2, 
    D(x2) ~  -(p3*x2) + p4*x1*x2
    ]

@named sys =  ODESystem(eqs)
simplified_sys=structural_simplify(sys)
    
u0=[ x1 => 1.0, x2 =>1.0]
p=[p1 =>1.5, p2 => 1.0, p3 => 3.0 , p4 => 1.0]


prob = ODEProblem(simplified_sys,u0,(0.0,10.0),p)
sol = solve(prob,Tsit5())


t = collect(range(0,stop=10,length=200))
randomized = VectorOfArray([(sol(t[i]) + 0.00randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)

cost_function = build_lsoptim_objective(prob,t,data,Tsit5())

# Ordering for me 
parameters(simplified_sys)

#4-element 
#Vector{Sym{Real, Base.ImmutableDict{DataType, Any}}}:
# p2
# p1
# p4
# p3

# Hand coded
idxs = [2,1,4,3]

x = [1.3,0.8,2.8,1.2]


res = optimize!(LeastSquaresProblem(x = x[idxs], f! = cost_function,
                output_length = length(t)*length(prob.u0)),
                LeastSquaresOptim.Dogleg(LeastSquaresOptim.LSMR()))

Checks out with

u0=[ x1 => 1.0, x2 =>1.0]
p=[p1 =>1.5, p2 => 1.0, p3 => 3.0 , p4 => 1.0]
p0 = [1.5, 1.0, 3.0, 1.0]

prob = ODEProblem(simplified_sys,u0,(0.0,10.0),p)
sol = solve(prob,Tsit5())
prob_ = ODEProblem(simplified_sys,u0,(0.0,10.0),p0[idxs])
sol_ = solve(prob_,Tsit5(), saveat = sol.t)
Array(sol) == Array(sol_)
1 Like