# Modeling Toolkit, linearize nonlinear ODE system without input and simplification

**URL:** <https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036>\
**Category:** General Usage\
**Tags:** modelingtoolkit, nonlinearsolve\
**Created:** [April 28, 2023, 11:32am UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036 "2023-04-28T11:32:31Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [April 28, 2023, 11:32am UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/1 "2023-04-28T11:32:31Z")

</div>

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

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

```

which I manage to solve using

```julia
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`)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 28, 2023, 12:43pm UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/2 "2023-04-28T12:43:56Z")

</div>

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

[![](https://global.discourse-cdn.com/julialang/original/3X/e/f/efc5408321efbc8e0ff1414297073a7ec6bac349.jpeg "How to Linearize a ModelingToolkit model") ](https://www.youtube.com/watch?v=-XOux-2XDGI)

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [April 28, 2023, 12:57pm UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/3 "2023-04-28T12:57:43Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 28, 2023, 1:06pm UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/4 "2023-04-28T13:06:36Z")

</div>

> [@BambOoxX](#):
>
> That may not currently be implemented, but wouldn’t it be useful to avoid a second `structural_simplify` call ?

`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

> [@BambOoxX](#):
>
> I guess this is because I have no inputs isn’t it ?

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

```julia
op = Dict(states(simple_ode) .=> last(sol.u))

```

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [April 28, 2023, 1:24pm UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/5 "2023-04-28T13:24:18Z")

</div>

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 !

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 27, 2024, 4:48pm UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/6 "2024-06-27T16:48:50Z")

</div>

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

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

```

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

```julia
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
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 ?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [June 27, 2024, 5:51pm UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/7 "2024-06-27T17:51:29Z")

</div>

> [@BambOoxX](#):
>
> Would you have any clues as to what happens here ?

No, but I’ve pinged someone who might 😊

---

<div class="post-metadata">

**Author:** ![cryptic.ax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cryptic.ax/32/220540_2.png) [@cryptic.ax](https://discourse.julialang.org/u/cryptic.ax)\
**Post date:** [June 27, 2024, 6:28pm UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/8 "2024-06-27T18:28:26Z")

</div>

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
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`?

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 27, 2024, 7:33pm UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/9 "2024-06-27T19:33:18Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 28, 2024, 11:21am UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/10 "2024-06-28T11:21:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![cryptic.ax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cryptic.ax/32/220540_2.png) [@cryptic.ax](https://discourse.julialang.org/u/cryptic.ax)\
**Post date:** [June 28, 2024, 11:37am UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/11 "2024-06-28T11:37:32Z")

</div>

> [@BambOoxX](#):
>
> 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

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.

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 28, 2024, 12:01pm UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/12 "2024-06-28T12:01:37Z")

</div>

This reproduces the problem.

```julia
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
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

```

```julia
  [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

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 28, 2024, 2:30pm UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/13 "2024-06-28T14:30:11Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 28, 2024, 3:41pm UTC](https://discourse.julialang.org/t/modeling-toolkit-linearize-nonlinear-ode-system-without-input-and-simplification/98036/14 "2024-06-28T15:41:19Z")

</div>

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

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

```

with v8.