# ModelingToolkit and Linearization

Continuing a former discussion in Optimization of parameters of a dynamic problem with ModelingToolkit - #4 by BambOoxX, 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

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

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

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

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

``````(; 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`

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

See

and

1 Like

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.

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.

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

• 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

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

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

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

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?

I’m working with

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

a `resolve` keeps things the same.

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

``````(; 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

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

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

when linearizing, unless you use analysis points.

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 !

1 Like

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

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

Working on the PR.

1 Like

EDIT: The PR was merged.

1 Like

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

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