Modeling Toolkit, linearize nonlinear ODE system without input and simplification

General Usage
,
1

I have a nonlinear ode system that has a complex structure, modeled with ModelingToolkit. This system shows a damped oscillatory response and is defined as

@named ode = ODESystem(eqs, t, <variables>, <parameters>, tspan=(0.0, 10.0))
simple_ode = structural_simplify(ode)

which I manage to solve using

prob = ODEProblem(simple_ode, <u0>, (0, 10), ps)
sol = solve(prob, Rodas5(), abstol=1e-8, reltol=1e-8)

Now I know the time response of my system, but I would like to get the state matrix A from it’s linearization around the final position so that I know its stability properties around that point.

I’ve tried a naive (; A, B, C, D), simplified_sys = linearize(simple_ode,[],[],t=last(sol.t),op=last(sol.u) since I do not care about the inputs and outputs, just the states, but I get a ERROR: ExtraVariablesSystemException: The system is unbalanced.
It looks like this function tries to re-simplify the system structure, but that goes wrong… and I do not know how to retrieve the state matrix through linearization_function which has a simplify=false option (which does not seem available in linearize)

2

To linearize the model, you pass the non-simplified system model into linearize. Have a look at this video for some additional details, tips and tricks

3

Thanks for the link !

That may not currently be implemented, but wouldn’t it be useful to avoid a second structural_simplify call ?

If I apply linearize to my initial ode, with (; A, B, C, D), simplified_sys = linearize(ode,[],[];t=last(sol.t),op = last(sol.u)) now I get a ERROR: MethodError: no method matching merge(::Dict{Any, Any}, ::Vector{Float64})

I guess this is because I have no inputs isn’t it ? Is there a way to call linearization without inputs ? My system is an autonomous one, so there is no input at all. Note : calling linearize(ode,[],[var];kwargs...) gives the same error

4

linearize will, in the presence of inputs, modify the model, so you have to perform structural_simplify again. If you have no inputs, you may also call calculate_jacobian which I believe works after simplification

No, it’s because op is supposed to be a Dict. Try

op = Dict(states(simple_ode) .=> last(sol.u))
3 Likes
5

Two out of two !

  • calculate_jacobian indeed works on the simplified system
  • indeed op when passed as a Dict works well, though the simplification complains that my model parameters are missing. Providing op = merge(Dict(states(simple_ode) .=> last(sol.u)), <parameters_dict>) solves this !

Thanks a lot !

1 Like