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

General Usage
, , , ,
1

Hi, I have confronted with something amazing !

I was trying parameter estimation on Lotka-volterra example in DiffeqParamEstim which uses DifferentialEquations.jl model specification and parameter estimation using LeastSquaresOptim.jl. That went pretty smooth. But the same example using ModelingToolkit specification and LeastSquaresOptim.jl fails throwing Bounds Error and asking for large maxiters. Increasing maxiters didn’t worked yet.

Another fact is that both produced similar simulation output. And without adding noise, same data is used for parameter estimation in both cases.

I would like to integrate aforementioned methods to my model and thats why I keep exploring all the many ways! Can someone help me out? Is that an issue internally with calling optimization in LeastSquaresOptim?

Here goes the package versions and code I have tried…

DiffEqParamEstim v1.23.1
DifferentialEquations v7.1.0
LeastSquaresOptim v0.8.3
ModelingToolkit v8.14.0
RecursiveArrayTools v2.29.2

using DifferentialEquations.jl :

using Random
Random.seed!(1234)

function f2(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end

u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(f2,u0,tspan,p)
sol = solve(prob,Tsit5())

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

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

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

using ModelingToolkit.jl :

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

x = [1.3,0.8,2.8,1.2]
res = optimize!(LeastSquaresProblem(x = x, f! = cost_function,
                output_length = length(t)*length(prob.u0)),
                LeastSquaresOptim.Dogleg(LeastSquaresOptim.LSMR()))
2

Have you checked, if the ordering of the MTK systems parameters hold without using the Pairs to initialize them? AFAIK p is not guaranteed to be in order.

3

Yes. I have checked it. In that case the simulation output is different from earlier, and I assumed that ordering might have taken differently.

But since I specify them in pairs no issues with ordering, right ? And does using pairs in model create some issue with parameter estimation? I hope it won’t because that worked with some other example I have already tried.

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
5

So that means using vector of pairs in model specification is the bottleneck for parameter estimation? The model I work upon has 11 parameters and how could I know the right permutation of vector in case if I specify without pairs?

Have to explore more on using Symbolics.jl. Thanks for the suggestion.

6

Not the bottleneck. It just gives you an unstable system due to the “misspecified” initial parameters. You simply cannot assume that the ordering is given in what you expect :).
But you could just use parameters(system) to get MTKs ordering and build up a permutation which gives you the right ordering, similar to above. So something along:

p_system = parameters(system)
my_p = [p1;p2;p3;p4]
map(my_p) do ps
     map(p_system) do pss
         isequal(ps, pss)
     end
end

Which would give you a matrix, where each parameter went.

7

Now it seems sensible. Thank you so much for your time and help…

1 Like