How to linearize a model

Specific Domains Modelling & Simulations
1

What I did so far:

  • I can define a model using ModelingToolkit.jl
  • I can simplify the system of equations
  • I can simulate the step response of the model

Now I would like to linearize the model. In theory very easy, just one line of code. But it doesn’t work, and I don’t understand how to fix the error.

My code:

using ModelingToolkit, OrdinaryDiffEq, DataInterpolations

Cp=[0.024384377390186646,0.03811829845979822,0.05731175742000688,0.0803120051369546,0.10526323380995135,0.1312349486430945,
    0.1588590278697544,0.1882995422205645,0.21937359147848456,0.25269290703920616,0.2882162631147021,0.3240355349501659,
    0.3607290350346774,0.39295261085207084,0.4141821251393675,0.4274331262027119,0.43636422894857274,0.4414767530303278,
    0.443412806515125,0.44365385445693295,0.44317448744425714,0.44218641629727234,0.4407405317795181,0.43888039054963845,
    0.4362875461540461,0.4325702155989231,0.4278606704093836,0.4224263790651815,0.4164272616252002,0.40987027258773945,
    0.4027832291503407,0.39518785578723936,0.38719847687832065]
TSR = 2.0:0.25:10.0

cp = CubicSpline(Cp, TSR)

# initial state
U0  = 9.0         # m/s
Pg_0 = 2.302e6    # W
w0 = 1.302

# model
@variables t ω(t)=w0 Pr(t) λ(t) Pg(t)=Pg_0
@parameters J R ρ A U=U0
D = Differential(t)

eqs = [J * D(ω) * ω ~ Pr - Pg,
       λ            ~ ω * R/U,
       Pr           ~ 0.5 * ρ * A * U^3 * cp(λ)]

@named sys = ODESystem(eqs, t)

# linearisation
matrices, simplified_sys = ModelingToolkit.linearize(sys, [Pg], [ω])

If I run this code I get the error message:

julia> include("src/mwe.jl")
ERROR: LoadError: ArgumentError: SymbolicUtils.BasicSymbolic{Real}[ωˍt(t)] are missing from the variable map.
Stacktrace:
  [1] throw_missingvars(vars::Vector{SymbolicUtils.BasicSymbolic{Real}})
    @ ModelingToolkit C:\Users\uwefechner\.julia\packages\ModelingToolkit\XgH0T\src\variables.jl:121
  [2] _varmap_to_vars(varmap::Dict{Any, Any}, varlist::Vector{SymbolicUtils.BasicSymbolic{Real}}; defaults::Dict{Any, Any}, check::Bool, toterm::typeof(ModelingToolkit.default_toterm))
    @ ModelingToolkit C:\Users\uwefechner\.julia\packages\ModelingToolkit\XgH0T\src\variables.jl:115
  [3] varmap_to_vars(varmap::Dict{Any, Any}, varlist::Vector{Any}; defaults::Dict{Any, Any}, check::Bool, toterm::Function, promotetoconcrete::Nothing, tofloat::Bool, use_union::Bool)
    @ ModelingToolkit C:\Users\uwefechner\.julia\packages\ModelingToolkit\XgH0T\src\variables.jl:84
  [4] get_u0_p(sys::ODESystem, u0map::Dict{Any, Any}, parammap::SciMLBase.NullParameters; use_union::Bool, tofloat::Bool)
    @ ModelingToolkit C:\Users\uwefechner\.julia\packages\ModelingToolkit\XgH0T\src\systems\diffeqs\abstractodesystem.jl:589
  [5] get_u0_p
    @ C:\Users\uwefechner\.julia\packages\ModelingToolkit\XgH0T\src\systems\diffeqs\abstractodesystem.jl:580 [inlined]
  [6] linearize(sys::ODESystem, lin_fun::ModelingToolkit.var"#130#133"{Vector{Num}, UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}, Vector{Any}, ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#531"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x85b51431, 0x1dd0388c, 0x97d6c70c, 0x979552d7, 0x94638b6f), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x33297964, 0x424aff99, 0x47f62688, 0x0209746d, 0x7cebcf48), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#583#generated_observed#539"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#13004503558829430822"), Symbol("##arg#8953363045742149305"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa9fbe3ad, 0x73fdf9d6, 0xa00f1e3c, 0x018b49fc, 0xc99f5b11), Nothing}, ForwardDiff.Chunk{1}}; t::Float64, op::Dict{Any, Any}, allow_input_derivatives::Bool, p::SciMLBase.NullParameters)
    @ ModelingToolkit C:\Users\uwefechner\.julia\packages\ModelingToolkit\XgH0T\src\systems\abstractsystem.jl:1511
  [7] linearize(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num}; op::Dict{Any, Any}, t::Float64, allow_input_derivatives::Bool, zero_dummy_der::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit C:\Users\uwefechner\.julia\packages\ModelingToolkit\XgH0T\src\systems\abstractsystem.jl:1563
  [8] linearize(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num})
    @ ModelingToolkit C:\Users\uwefechner\.julia\packages\ModelingToolkit\XgH0T\src\systems\abstractsystem.jl:1554
  [9] top-level scope
    @ C:\Users\uwefechner\repos\WindTurbines\src\mwe.jl:30
 [10] include(fname::String)
    @ Base.MainInclude .\client.jl:478
 [11] top-level scope
    @ REPL[1]:1
in expression starting at C:\Users\uwefechner\repos\WindTurbines\src\mwe.jl:30

Any idea?

2

See if solving the system a bit makes it easier on the function:

eqs = [D(ω)  ~ (Pr - Pg) / J * ω,
       0            ~ ω * R/U - λ,
       0           ~ 0.5 * ρ * A * U^3 * cp(λ) - Pr]
3

This helps indeed. The following code works:

using ModelingToolkit, OrdinaryDiffEq, DataInterpolations

Cp=[0.024384377390186646,0.03811829845979822,0.05731175742000688,0.0803120051369546,0.10526323380995135,0.1312349486430945,
    0.1588590278697544,0.1882995422205645,0.21937359147848456,0.25269290703920616,0.2882162631147021,0.3240355349501659,
    0.3607290350346774,0.39295261085207084,0.4141821251393675,0.4274331262027119,0.43636422894857274,0.4414767530303278,
    0.443412806515125,0.44365385445693295,0.44317448744425714,0.44218641629727234,0.4407405317795181,0.43888039054963845,
    0.4362875461540461,0.4325702155989231,0.4278606704093836,0.4224263790651815,0.4164272616252002,0.40987027258773945,
    0.4027832291503407,0.39518785578723936,0.38719847687832065]
TSR = 2.0:0.25:10.0

cp = CubicSpline(Cp, TSR)

# initial state
U0  = 9.0         # m/s
Pg_0 = 2.302e6    # W
w0 = 1.302

# model
@variables t ω(t)=w0 Pr(t) λ(t) Pg(t)=Pg_0
@parameters J=4.0470e+07 R=63 ρ=1.225 A=1.2469e+04 U=U0
D = Differential(t)

eqs = [D(ω)  ~ (Pr - Pg) / J * ω,
       0     ~ ω * R/U - λ,
       0     ~ 0.5 * ρ * A * U^3 * cp(λ) - Pr]

@named sys = ODESystem(eqs, t)

# linearisation
matrices, simplified_sys = ModelingToolkit.linearize(sys, [Pg], [ω])

Is manually rewriting always required?

4

It shouldn’t be. Open an issue. It’s just something that some passes don’t handle well and we need to more uniformly handle it. I also wonder if you structurally simplify before calling linearize if it will similarly do better.

5

The problem

is due to MTK inserting a dummy derivative variable and there is no initial condition specified for this variable. You can use the function ModelingToolkit.missing_variable_defaults to automatically set all dummy derivatives to 0.

Check out this video for additional details on linearization in MTK, as well as the docs on analysis points

1 Like
6

You can’t and shouldn’t do this. When you linearize, there are unbound inputs present and structural_simplify will thus complain that the system is unbalanced. linearize will internally call most of the machinery of structural_simplify, with some additional steps to properly handle the linearization inputs.

7

I created a bug report: Linearizing of a simple model fails · Issue #2185 · SciML/ModelingToolkit.jl · GitHub
@baggepinnen Your suggestion does not work, I get the same error message.

8

Once again, the problem is that you haven’t specified ωˍt. Try passing the keyword argument op = Dict(D(ω) => 0) to linearize

9

I tried, but this doesn’t help, just leads to another error.

julia> include("src/mwe.jl")
ERROR: LoadError: Input derivatives appeared in expressions (-g_z\g_u != 0), the following inputs appeared differentiated: Any[Pg(t)]. Call `linear_statespace` with keyword argument `allow_input_derivatives = true` to allow this and have the returned `B` matrix be of double width (2), where the last 1 inputs are the derivatives of the first 1 inputs.
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:35
 [2] linearize(sys::ODESystem, lin_fun::ModelingToolkit.var"#130#133"{Vector{Num}, UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}, Vector{Any}, ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#531"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x85b51431, 0x1dd0388c, 0x97d6c70c, 0x979552d7, 0x94638b6f), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x33297964, 0x424aff99, 0x47f62688, 0x0209746d, 0x7cebcf48), Nothing}}, Matrix{Float64}, Nothing, Nothing, 
Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#583#generated_observed#539"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#13004503558829430822"), Symbol("##arg#8953363045742149305"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa9fbe3ad, 0x73fdf9d6, 0xa00f1e3c, 0x018b49fc, 0xc99f5b11), Nothing}, ForwardDiff.Chunk{1}}; t::Float64, op::Dict{Num, Int64}, allow_input_derivatives::Bool, p::SciMLBase.NullParameters)
   @ ModelingToolkit C:\Users\uwefechner\.julia\packages\ModelingToolkit\XgH0T\src\systems\abstractsystem.jl:1544
 [3] linearize(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num}; op::Dict{Num, Int64}, t::Float64, allow_input_derivatives::Bool, zero_dummy_der::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit C:\Users\uwefechner\.julia\packages\ModelingToolkit\XgH0T\src\systems\abstractsystem.jl:1563
 [4] top-level scope
   @ C:\Users\uwefechner\repos\WindTurbines\src\mwe.jl:32
 [5] include(fname::String)
   @ Base.MainInclude .\client.jl:478
 [6] top-level scope
   @ REPL[2]:1
in expression starting at C:\Users\uwefechner\repos\WindTurbines\src\mwe.jl:32
10

And if I add the parameter suggested by the error message it works, but the result is not what I would expect. I get the result:

julia> matrices
(A = [0.0 1.0; 0.0 -0.019333043055042454], B = [0.0 0.0; 0.0 -1.8978234621841593e-8], C = [1.0 0.0], D = [0.0 0.0])

Expected result with the manually simplified equation system:

julia> matrices
(A = [-0.019336236721183935;;], B = [-1.8978234621841593e-8;;], C = [1.0;;], D = [0.0;;])

I mean, it is not a problem for me if I have to manually simplify an equation system
before linearization, but if this is required it should be documented.

11

The result is not wrong, it’s just not simplified all the way down to an ODE system. You can read more on the problem and the solution using linear algebra here
https://juliacontrol.github.io/ControlSystemsMTK.jl/dev/#Internals:-Transformation-of-non-proper-models-to-proper-statespace-form

It’s all present in the docs as well, on the page labeled “linearization”:
https://docs.sciml.ai/ModelingToolkit/dev/basics/Linearization/#Input-derivatives

12

Your suggestion does not work, just results in a new error:

Code:

using ModelingToolkit, OrdinaryDiffEq, DataInterpolations
import DescriptorSystems

Cp=[0.024384377390186646,0.03811829845979822,0.05731175742000688,0.0803120051369546,0.10526323380995135,0.1312349486430945,
    0.1588590278697544,0.1882995422205645,0.21937359147848456,0.25269290703920616,0.2882162631147021,0.3240355349501659,
    0.3607290350346774,0.39295261085207084,0.4141821251393675,0.4274331262027119,0.43636422894857274,0.4414767530303278,
    0.443412806515125,0.44365385445693295,0.44317448744425714,0.44218641629727234,0.4407405317795181,0.43888039054963845,
    0.4362875461540461,0.4325702155989231,0.4278606704093836,0.4224263790651815,0.4164272616252002,0.40987027258773945,
    0.4027832291503407,0.39518785578723936,0.38719847687832065]
TSR = 2.0:0.25:10.0

cp = CubicSpline(Cp, TSR)

# initial state
U0  = 9.0         # m/s
Pg_0 = 2.302e6    # W
w0 = 1.302

# model
@variables t ω(t)=w0 Pr(t) λ(t) Pg(t)=Pg_0
@parameters J=4.0470e+07 R=63 ρ=1.225 A=1.2469e+04 U=U0

D = Differential(t)

simple=false

if simple
        eqs = [D(ω)  ~ (Pr - Pg) / (J * ω),
               0     ~ ω * R/U - λ,
               0     ~ 0.5 * ρ * A * U^3 * cp(λ) - Pr]
 else
        eqs = [J * D(ω) * ω ~ Pr - Pg,
               λ            ~ ω * R/U,
               Pr           ~ 0.5 * ρ * A * U^3 * cp(λ)]
 end

@named sys = ODESystem(eqs, t)

# linearisation
matrices, simplified_sys = ModelingToolkit.linearize(sys, [Pg], [ω];  op = Dict(D(ω) => 0), allow_input_derivatives = true)

simplified_sys=DescriptorSystems.dss2ss(simplified_sys)

nothing

Error message:

julia> include("src/mwe.jl")
ERROR: LoadError: MethodError: no method matching dss2ss(::ODESystem)

Closest candidates are:
  dss2ss(::DescriptorSystems.DescriptorStateSpace{T, ET} where ET<:Union{LinearAlgebra.UniformScaling{Bool}, Matrix{T}}) where T
   @ DescriptorSystems C:\Users\uwefechner\.julia\packages\DescriptorSystems\cZGte\src\order_reduction.jl:143
  dss2ss(::DescriptorSystems.DescriptorStateSpace{T, ET} where ET<:Union{LinearAlgebra.UniformScaling{Bool}, Matrix{T}}, ::Vector; state_mapping, simple_infeigs, fast, atol, atol1, atol2, rtol) where 
l}, Matrix{T}}, ::Vector; state_mapping, simple_infeigs, fast, atol, atol1, atol2, rtol) where T
   @ DescriptorSystems C:\Users\uwefechner\.julia\packages\DescriptorSystems\cZGte\src\order_reduction.jl:143

Stacktrace:
 [1] top-level scope
   @ C:\Users\uwefechner\repos\WindTurbines\src\mwe.jl:42
 [2] include(fname::String)
   @ Base.MainInclude .\client.jl:478
 [3] top-level scope
   @ REPL[1]:1
in expression starting at C:\Users\uwefechner\repos\WindTurbines\src\mwe.jl:42
13

This model may be converted to a proper statespace model (if the system is indeed proper) using DescriptorSystems.dss2ss. All of this is handled automatically by named_ss.

Use named_ss

I’m spending my free time trying to help you use software you’re not paying for from my phone, please show some gratitude or I promise I’ll stop doing this immediately.

14

Dear Fredrik,

I am very grateful for the effort you put into the development of the Julia control ecosystem, but I would probably indeed prefer if you would not post suggestions from your phone after work that you did not have the time to try yourself.

Better wait until you have enough time for a good answer than giving a quick answer. I am not in a hurry. And let’s try to keep the good spirit of creating a free alternative to the commercial, closed source software!

15

After a month of additional experiance I can say that there were two issues in my original example:

  • on the left hand side of any equation there should either be a differential or a variable or zero, but not an expression
  • all parameters need an initial value

The following code works:

using ModelingToolkit, OrdinaryDiffEq, DataInterpolations

Cp=[0.024384377390186646,0.03811829845979822,0.05731175742000688,0.0803120051369546,0.10526323380995135,0.1312349486430945,
    0.1588590278697544,0.1882995422205645,0.21937359147848456,0.25269290703920616,0.2882162631147021,0.3240355349501659,
    0.3607290350346774,0.39295261085207084,0.4141821251393675,0.4274331262027119,0.43636422894857274,0.4414767530303278,
    0.443412806515125,0.44365385445693295,0.44317448744425714,0.44218641629727234,0.4407405317795181,0.43888039054963845,
    0.4362875461540461,0.4325702155989231,0.4278606704093836,0.4224263790651815,0.4164272616252002,0.40987027258773945,
    0.4027832291503407,0.39518785578723936,0.38719847687832065]
TSR = 2.0:0.25:10.0

cp = CubicSpline(Cp, TSR)

# initial state
U0  = 9.0         # m/s
Pg_0 = 2.302e6    # W
w0 = 1.302
J0 = 4.047e7
R0 = 63
ρ0 = 1.225
A0 = 12469.0

# model
@variables t ω(t)=w0 Pr(t) λ(t) Pg(t)=Pg_0
@parameters J=J0 R=R0 ρ=ρ0 A=A0 U=U0
D = Differential(t)

eqs = [D(ω) ~ (Pr - Pg)/(J * ω),
       λ    ~ ω * R/U,
       Pr   ~ 0.5 * ρ * A * U^3 * cp(λ)]

@named sys = ODESystem(eqs, t)

# linearisation
matrices, simplified_sys = ModelingToolkit.linearize(sys, [Pg], [ω])

Alternatively you can also make a symbolic linearization, but then you (currently) cannot use a piecewise defined function like cp() in the example above. I replaced cp() with the quadratic approximation a*λ^2 + b*λ + c, and now I can do a symbolic linearization:

using ModelingToolkit, OrdinaryDiffEq

# model
@variables t ω(t) Pr(t) λ(t) Pg(t)
@parameters J R ρ A U a b c
D = Differential(t)

eqs = [D(ω) ~ (Pr - Pg)/(J * ω),
       λ    ~ ω * R/U,
       Pr   ~ 0.5 * ρ * A * U^3 * (a*λ^2 + b*λ + c)]

@named sys = ODESystem(eqs, t)

# linearisation
matrices, simplified_sys = ModelingToolkit.linearize_symbolic(sys, [Pg], [ω])
println("A: ", matrices.A[:])
println("B: ", matrices.B[:])
println("C: ", matrices.C[:])
println("D: ", matrices.D[:])

This has the output:

A: Num[(0.5A*ρ*(U^3)*((R*b) / U + (2a*(R^2)*ω(t)) / (U^2))) / (J*ω(t)) - J*((0.5A*ρ*(c + a*(((R*ω(t)) / U)^2) + (R*b*ω(t)) / U)*(U^3) - Pg(t)) / ((J^2)*(ω(t)^2)))]
B: Num[-1 / (J*ω(t))]
C: Num[1]
D: Num[0]

Thanks @baggepinnen and @ChrisRackauckas for the good work! :slight_smile: