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)
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
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
This reads
(; 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 ?
