# ModelingToolkit and Linearization

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241>\
**Category:** Modelling & Simulations\
**Tags:** question, modelingtoolkit\
**Created:** [April 19, 2024, 2:10pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241 "2024-04-19T14:10:10Z")\
**Posts on this page:** 16\
**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 19, 2024, 2:10pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/1 "2024-04-19T14:10:10Z")

</div>

Continuing a former discussion in [Optimization of parameters of a dynamic problem with ModelingToolkit - #4 by BambOoxX](https://discourse.julialang.org/t/optimization-of-parameters-of-a-dynamic-problem-with-modelingtoolkit/103784/4), I’m now trying to model and linearize a mass-spring-damper system as an example to test MTKs capabilities.

# Objective

The idea is to be able to compute first a steady state solution, then linearize around it to retrieve a state-space model with given inputs and outputs.

# Model

The model is built with

```julia-auto
using ModelingToolkit, ModelingToolkitStandardLibrary, DifferentialEquations
import ModelingToolkit: connect
using GLMakie

@variables t

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

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

eqs = Equation[
    ModelingToolkit.connect(gval.output, g.f)
    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])

sys = structural_simplify(model)

```

# Simulation

A _long_-time transient solution, and a steady-state solution can be obtained with

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

f,ax,p = plot(sol)
hlines!(ax,sol_ss.u[1:1], linestyle=:dash)
hlines!(ax, sol_ss.u[2:2], linestyle=:dash)
axislegend(ax)

```

which gives a nice

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/9/c93d9a0fc2c7a1561754283fe1561ba5a3d4baa8.png)

# Linearization

Now the question is how to obtain the steady state model. Obviously linearization is the way to go. That gives (without inputs or outpus)

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

```

where only `A` is not empty.

# Question

Now I’d like to add some input/output handling. I figured I should just take at least one input and one output variables in `model` and see how it goes. For instance, the most straightforward way to go seemed to use `ModelingToolkit.inputs / outputs`

This reads

```julia-auto
(; A, B, C, D), simplified_sys = linearize(model, ModelingToolkit.inputs(model), ModelingToolkit.outputs(model); t=Inf, op=Dict(unknowns(sys) .=> sol_ss))

```

but that fails with

```julia-auto
julia> (; A, B, C, D), simplified_sys = linearize(model, ModelingToolkit.inputs(model), ModelingToolkit.outputs(model); t=Inf, op=Dict(unknowns(sys) .=> sol_ss))
ERROR: ExtraEquationsSystemException: The system is unbalanced. There are 23 highest order derivative variables and 24 equations.
More equations than variables, here are the potential extra equation(s):
 0 ~ damper₊flange_a₊s(t) - mass₊s(t)
 0 ~ -gval₊output₊u(t) + g₊f₊u(t)
 0 ~ -damper₊flange_b₊s(t) + tfix₊flange₊s(t)

```

what am I doing wrong ?

---

<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 19, 2024, 2:24pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/2 "2024-04-19T14:24:52Z")

</div>

See

[![](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)

and

[![](https://global.discourse-cdn.com/julialang/original/3X/6/d/6ddf069bf7cd4092f2a81b7df511a3a11feb981a.jpeg "Linear Analysis of ModelingToolkit Models | Fredrik Bagge Carlson | JuliaCon 2023") ](https://www.youtube.com/watch?v=aThEtmuuazo)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 19, 2024, 2:48pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/3 "2024-04-19T14:48:55Z")

</div>

You’re not doing anything wrong. The linearization pass is missing a pretty crucial aspect for DAEs and we’re going to be fixing that up next week. Effectively, if you do the linearization at a consistently defined state then you’re fine. But the linearization functionality is missing the ability to enforce that consistency.

This is actually causing MTK tests to fail right now and so it’s an extremely high priority that we plan to fix ASAP, though the person who was working on the PR got sick. It didn’t show because inputs and outputs used to set defaults instead of guesses, which isn’t necessarily correct but removing that exposed this issue in the linearization.

---

<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 19, 2024, 2:51pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/4 "2024-04-19T14:51:45Z")

</div>

This problem is not about initialization, the problem is explained in the first video linked above. Linearization adds one equation per input, so if the model is fully connected such that it simulates, it will not linearize. The solution is to either remove the input binding, or use the analysis-points feature.

---

<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 19, 2024, 4:05pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/5 "2024-04-19T16:05:27Z")

</div>

@baggepinnen Thanks ! I finally understood the concept of the analysis points with your video, I think I overlooked that at first.

So to be able to conform to your example, I added

- a custom `PositionSensor` that is absent from `ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition`. (the one from `ModelingToolkitStandardLibrary.Mechanical.Position` has another connector type).
- one input and one output at the proper places

However, I get a `display` error when assembling the equations with the analysis points, and the linearization fails with

```julia
ERROR: UndefVarError: `states` not defined
Stacktrace:
  [1] ap_var(sys::ODESystem)
    @ ModelingToolkitStandardLibrary.Blocks user\.julia\packages\ModelingToolkitStandardLibrary\nf8EA\src\Blocks\analysis_points.jl:26
  [2] (::ModelingToolkitStandardLibrary.Blocks.var"#198#199"{…})(ap::ModelingToolkitStandardLibrary.Blocks.AnalysisPoint, ns::Nothing)
    @ ModelingToolkitStandardLibrary.Blocks user\.julia\packages\ModelingToolkitStandardLibrary\nf8EA\src\Blocks\analysis_points.jl:426
  [3] generate_connection_set!(connectionsets::Vector{…}, domain_csets::Vector{…}, sys::ODESystem, find::ModelingToolkitStandardLibrary.Blocks.var"#177#178"{…}, replace::ModelingToolkitStandardLibrary.Blocks.var"#198#199"{…}, scalarize::Bool, namespace::Nothing)
    @ ModelingToolkit user\.julia\packages\ModelingToolkit\QP8aH\src\systems\connectors.jl:325
  [4] generate_connection_set!(connectionsets::Vector{…}, domain_csets::Vector{…}, sys::ODESystem, find::Function, replace::Function, scalarize::Bool)
    @ ModelingToolkit user\.julia\packages\ModelingToolkit\QP8aH\src\systems\connectors.jl:312
  [5] #generate_connection_set#331
    @ user\.julia\packages\ModelingToolkit\QP8aH\src\systems\connectors.jl:302 [inlined]
  [6] generate_connection_set
    @ user\.julia\packages\ModelingToolkit\QP8aH\src\systems\connectors.jl:298 [inlined]
  [7] expand_connections(sys::ODESystem, find::Function, replace::Function; debug::Bool, tol::Float64, scalarize::Bool)
    @ ModelingToolkit user\.julia\packages\ModelingToolkit\QP8aH\src\systems\connectors.jl:485
  [8] expand_connections
    @ user\.julia\packages\ModelingToolkit\QP8aH\src\systems\connectors.jl:483 [inlined]
  [9] linearization_function(sys::ODESystem, input_name::Vector{…}, output_name::Vector{…}; loop_openings::Nothing, kwargs::@Kwargs{…})
    @ ModelingToolkitStandardLibrary.Blocks user\.julia\packages\ModelingToolkitStandardLibrary\nf8EA\src\Blocks\analysis_points.jl:436
 [10] linearize(sys::ODESystem, input_name::Vector{…}, output_name::Vector{…}; loop_openings::Nothing, kwargs::@Kwargs{…})
    @ ModelingToolkitStandardLibrary.Blocks user\.julia\packages\ModelingToolkitStandardLibrary\nf8EA\src\Blocks\analysis_points.jl:486

```

Here is the complete modified example

```julia
using ModelingToolkit, ModelingToolkitStandardLibrary, DifferentialEquations
import ModelingToolkit: connect
using GLMakie

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

f, ax, p = plot(sol)
hlines!(ax, sol_ss.u[1:1], linestyle=:dash)
hlines!(ax, sol_ss.u[2:2], linestyle=:dash)
axislegend(ax)

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

```

---

<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 19, 2024, 4:16pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/6 "2024-04-19T16:16:46Z")

</div>

The problem appears to be that you have ended up with incompatible versions of MTK and MTKstdlib, that should not have been possible 🤔

`unknowns` used to be called `states` in a previous version of MTK. What’s the version of MTKstlib?

---

<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 19, 2024, 4:19pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/7 "2024-04-19T16:19:54Z")

</div>

I’m working with

```julia
⌃ [961ee093] ModelingToolkit v9.11.0
  [16a59e39] ModelingToolkitStandardLibrary v2.6.0

```

a `resolve` keeps things the same.

---

<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 19, 2024, 4:25pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/8 "2024-04-19T16:25:17Z")

</div>

I watched the first video you linked to, but I can’t figure what I am doing differently than what you did in your example.

In your system, everything looks connected as in mine, yet, if I do

```julia
(; A, B, C, D), simplified_sys = linearize(model, [g.f.u], [pos.output.u])

```

I keep getting the `ExtraEquationSystemException`, which is the one that I have and you don’t ? Especially since I follow the same logic than you

1. Form connections
2. Simulate
3. 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 19, 2024, 6:32pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/9 "2024-04-19T18:32:00Z")

</div>

In the video, I indicate that you have to remove the connection

```julia
ModelingToolkit.connect(gval.output, g.f)

```

when linearizing, _unless_ you use analysis points.

---

<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 19, 2024, 7:44pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/10 "2024-04-19T19:44:15Z")

</div>

Ok, then again, my apologies, I should have paid more attention to your video. The analysis points are surely useful to avoid changing the model’s equations between simulation and linearization !

---

<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 20, 2024, 3:01pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/11 "2024-04-20T15:01:53Z")

</div>

For the record, updating `ModelingToolkit` to `v9.12.0` and `ModelingToolkitStandardLibrary` to `v2.7.0` did not change the `states` problem. This line is still present on #main [ModelingToolkitStandardLibrary.jl/src/Blocks/analysis\_points.jl at fb8a224840fce60c25adf8aa48ae89cf7ff13a9e · SciML/ModelingToolkitStandardLibrary.jl · GitHub](https://github.com/SciML/ModelingToolkitStandardLibrary.jl/blob/fb8a224840fce60c25adf8aa48ae89cf7ff13a9e/src/Blocks/analysis_points.jl#L26)

---

<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 20, 2024, 3:07pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/12 "2024-04-20T15:07:44Z")

</div>

Would you mind opening an issue describing the problem, or a PR fixing it?

---

<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 20, 2024, 3:41pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/13 "2024-04-20T15:41:21Z")

</div>

Working on the PR.

---

<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 20, 2024, 3:52pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/14 "2024-04-20T15:52:09Z")

</div>

See [Fix `states` -\> `unknowns` by BambOoxX · Pull Request #289 · SciML/ModelingToolkitStandardLibrary.jl · GitHub](https://github.com/SciML/ModelingToolkitStandardLibrary.jl/pull/289)

EDIT: The PR was merged.

---

<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 22, 2024, 7:45am UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/15 "2024-04-22T07:45:10Z")

</div>

Do you have a fully working solution for your problem now?

---

<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 22, 2024, 7:44pm UTC](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/16 "2024-04-22T19:44:18Z")

</div>

Yes, for this example everything works fine now (using #main). I have another question, but that is for another thread ! Thanks.