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
6

@baggepinnen Sorry to come back to this, but I’ve updated my ModelingToolkit to use v9(.20) and unfortunately I get weird errors.

(; A, B, C, D), simplified_sys = linearize(ode, [<inputs>], [<outputs>]; t=Inf, op=op)

where op is a Dict{Num, Float64} I get

julia> (; A, B, C, D), simplified_sys = linearize(ode, [fext], [θ]; t=10, op=op)
ERROR: BoundsError: attempt to access 13-element Vector{Vector{Int64}} at index [14]
Stacktrace:
  [1] getindex
    @ Base .\essentials.jl:13 [inlined]
  [2] đť‘‘neighbors
    @ ModelingToolkit.BipartiteGraphs <user>\.julia\packages\ModelingToolkit\WiCz5\src\bipartite_graph.jl:369 [inlined]
  [3] đť‘‘neighbors
    @ ModelingToolkit.BipartiteGraphs <user>\.julia\packages\ModelingToolkit\WiCz5\src\bipartite_graph.jl:368 [inlined]
  [4] check_consistency(state::TearingState{ODESystem}, orig_inputs::Set{Any})
    @ ModelingToolkit.StructuralTransformations <user>\.julia\packages\ModelingToolkit\WiCz5\src\structural_transformation\utils.jl:96
  [5] _structural_simplify!(state::TearingState{…}, io::Tuple{…}; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, dummy_derivative::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\systemstructure.jl:691
  [6] structural_simplify!(state::TearingState{…}, io::Tuple{…}; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\systemstructure.jl:634
  [7] __structural_simplify(sys::ODESystem, io::Tuple{Vector{Num}, Vector{Num}}; simplify::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\systems.jl:83
  [8] structural_simplify(sys::ODESystem, io::Tuple{Vector{Num}, Vector{Num}}; simplify::Bool, split::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\systems.jl:24
  [9] structural_simplify
    @ <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\systems.jl:21 [inlined]
 [10] io_preprocessing(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num}; simplify::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\abstractsystem.jl:1823
 [11] linearization_function(sys::ODESystem, inputs::Vector{…}, outputs::Vector{…}; simplify::Bool, initialize::Bool, op::Dict{…}, p::SciMLBase.NullParameters, zero_dummy_der::Bool, initialization_solver_alg::GeneralizedFirstOrderAlgorithm{…}, kwargs::@Kwargs{})
    @ ModelingToolkit <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\abstractsystem.jl:1877
 [12] linearization_function
    @ <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\abstractsystem.jl:1867 [inlined]
 [13] linearize(sys::ODESystem, inputs::Vector{…}, outputs::Vector{…}; op::Dict{…}, t::Int64, allow_input_derivatives::Bool, zero_dummy_der::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\abstractsystem.jl:2325
 [14] top-level scope

If I do the same, but without inputs or outputs specified, I get another error :

julia> (; A, B, C, D), simplified_sys = linearize(ode, [], []; t=10, op=op)
ERROR: MethodError: no method matching parameter_values(::ModelingToolkit.MTKParameters{…})

Closest candidates are:
  parameter_values(::ModelingToolkit.MTKParameters, ::ModelingToolkit.ParameterIndex)
   @ ModelingToolkit <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\parameter_buffer.jl:318
  parameter_values(::Any, ::Any)
   @ SymbolicIndexingInterface <user>\.julia\packages\SymbolicIndexingInterface\A1VUA\src\parameter_indexing.jl:5
  parameter_values(::SymbolicIndexingInterface.ProblemState)
   @ SymbolicIndexingInterface <user>\.julia\packages\SymbolicIndexingInterface\A1VUA\src\problem_state.jl:20
  ...

Stacktrace:
 [1] set_parameter!(sys::ModelingToolkit.MTKParameters{…}, val::Float64, idx::Nothing)
   @ SymbolicIndexingInterface <user>\.julia\packages\SymbolicIndexingInterface\A1VUA\src\parameter_indexing.jl:17
 [2] (::SymbolicIndexingInterface.SetParameterIndex{…})(prob::ModelingToolkit.MTKParameters{…}, val::Float64)
   @ SymbolicIndexingInterface <user>\.julia\packages\SymbolicIndexingInterface\A1VUA\src\parameter_indexing.jl:190
 [3] (::SymbolicIndexingInterface.ParameterHookWrapper{…})(prob::ModelingToolkit.MTKParameters{…}, args::Float64)
   @ SymbolicIndexingInterface <user>\.julia\packages\SymbolicIndexingInterface\A1VUA\src\parameter_indexing.jl:158
 [4] (::ModelingToolkit.var"#311#326"{…})(u::Vector{…}, p::Dict{…}, t::Int64)
   @ ModelingToolkit <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\abstractsystem.jl:1988
 [5] linearize(sys::ODESystem, lin_fun::ModelingToolkit.var"#311#326"{…}; t::Int64, op::Dict{…}, allow_input_derivatives::Bool, p::SciMLBase.NullParameters)
   @ ModelingToolkit <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\abstractsystem.jl:2280
 [6] linearize(sys::ODESystem, inputs::Vector{…}, outputs::Vector{…}; op::Dict{…}, t::Int64, allow_input_derivatives::Bool, zero_dummy_der::Bool, kwargs::@Kwargs{})
   @ ModelingToolkit <user>\.julia\packages\ModelingToolkit\WiCz5\src\systems\abstractsystem.jl:2331
 [7] top-level scope

Would you have any clues as to what happens here ?

7

No, but I’ve pinged someone who might :blush:

8

Hi! It’s hard to debug this without an MWE I can run. But I have one question that might help figure this out:

julia> (; A, B, C, D), simplified_sys = linearize(ode, [fext], [θ]; t=10, op=op)

How exactly did you obtain fext and θ? Are they the result of @variables? Or something like fext = ode.fext?

9

Yeah, it is just difficult to provide a MWE since the model setup is a bit long, sorry about that.
fext / \theta result from an @variables indeed. Just as in v8

10

@baggepinnen @cryptic.ax In addition, though it’s not a surprise, the same problems arise when using named_ss from ControlSystemsMTK…

11

I’m afraid I can’t help more without code to debug. Even if it isn’t minimal, just something that reproduces the error would be great. If it’s really big, you could put it in a GitHub Gist.

12

This reproduces the problem.

using ModelingToolkit, ModelingToolkitStandardLibrary, DifferentialEquations
import ModelingToolkit: connect

@variables t
D = Differential(t)

import ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition: Mass, Spring, Damper, Fixed, Force, Flange
import ModelingToolkitStandardLibrary.Blocks: Constant, RealOutput


@mtkmodel PositionSensor begin
    @components begin
        flange = Flange()
        output = RealOutput()
    end

    @variables begin
        s(t) = 0.0
    end

    @equations begin
        output.u ~ s
        flange.f ~ 0.0
    end
end


@named mass = Mass(m=1, s=1.0)
@named spring = Spring(k=10, l=0.5)
@named damper = Damper(d=2)
@named tfix = Fixed(s_0=10)
@named g = Force()
@named gval = Constant(k=-9.81)
@named pos = PositionSensor()

eqs = Equation[
    ModelingToolkit.connect(gval.output, :u, g.f)
    ModelingToolkit.connect(mass.flange, :y, pos.flange)
    ModelingToolkit.connect(g.flange, mass.flange)
    ModelingToolkit.connect(spring.flange_a, mass.flange)
    ModelingToolkit.connect(damper.flange_a, mass.flange)
    ModelingToolkit.connect(spring.flange_b, tfix.flange)
    ModelingToolkit.connect(damper.flange_b, tfix.flange)
]

@named model = ODESystem(eqs, t; systems=[tfix, spring, damper, mass, g, gval, pos])

sys = structural_simplify(model)

initial_conditions = [
    mass.s => 0.0,
]

prob = ODEProblem(sys, initial_conditions, (0.0, 100.0))
sol = solve(prob, Rodas5(), abstol=1e-8, reltol=1e-8, saveat=0:0.1:20)

prob_ss = SteadyStateProblem(prob)
sol_ss = solve(prob_ss, DynamicSS(Rodas5()), abstol=1e-8, reltol=1e-8)

(; A, B, C, D), simplified_sys = linearize(model, [:u], [:y]; t=Inf, op=Dict(unknowns(sys) .=> sol_ss))

Version

julia> versioninfo()
Julia Version 1.10.0
Commit 3120989f39 (2023-12-25 18:01 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 8 Ă— 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, tigerlake)
  Threads: 1 on 8 virtual cores

  [0c46a032] DifferentialEquations v7.13.0
  [961ee093] ModelingToolkit v9.20.0
  [16a59e39] ModelingToolkitStandardLibrary v2.7.2

EDIT : just updated the M(N)WE to show the bug properly

13

I just ran other examples from the documentation, and they run fine with the same environment.

14

I confirm that the MWE works with the last line replaced by

(; A, B, C, D), simplified_sys = linearize(model, [], []; t=Inf, op=Dict(states(sys) .=> sol_ss))

with v8.