# How to linearize a model

**URL:** <https://discourse.julialang.org/t/how-to-linearize-a-model/99886>\
**Category:** Modelling & Simulations\
**Tags:** modelingtoolkit\
**Created:** [June 5, 2023, 12:04pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886 "2023-06-05T12:04:25Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [June 5, 2023, 12:04pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/1 "2023-06-05T12:04:25Z")

</div>

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:

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

---

<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:** [June 5, 2023, 12:13pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/2 "2023-06-05T12:13:46Z")

</div>

> [@ufechner7](#):
>
> ```julia
> eqs = [J * D(ω) * ω ~ Pr - Pg,
> λ ~ ω * R/U,
> Pr ~ 0.5 * ρ * A * U^3 * cp(λ)]
> 
> ```

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

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

```

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [June 5, 2023, 12:33pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/3 "2023-06-05T12:33:48Z")

</div>

> [@ChrisRackauckas](#):
>
> ```julia
> eqs = [D(ω) ~ (Pr - Pg) / J * ω,
> 0 ~ ω * R/U - λ,
> 0 ~ 0.5 * ρ * A * U^3 * cp(λ) - Pr]
> 
> ```

This helps indeed. The following code works:

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

---

<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:** [June 5, 2023, 12:40pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/4 "2023-06-05T12:40:08Z")

</div>

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.

---

<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 5, 2023, 1:24pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/5 "2023-06-05T13:24:44Z")

</div>

The problem

> [@ufechner7](#):
>
> `ArgumentError: SymbolicUtils.BasicSymbolic{Real}[ωˍt(t)] are missing from the variable map.`

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](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/)

[![](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:** ![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 5, 2023, 1:26pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/6 "2023-06-05T13:26:50Z")

</div>

> [@ChrisRackauckas](#):
>
> I also wonder if you structurally simplify before calling linearize if it will similarly do better.

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.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [June 6, 2023, 5:30pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/7 "2023-06-06T17:30:50Z")

</div>

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

---

<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 6, 2023, 5:54pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/8 "2023-06-06T17:54:46Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [June 6, 2023, 6:14pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/9 "2023-06-06T18:14:08Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [June 6, 2023, 6:19pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/10 "2023-06-06T18:19:23Z")

</div>

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

---

<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 6, 2023, 6:30pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/11 "2023-06-06T18:30:54Z")

</div>

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](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](https://docs.sciml.ai/ModelingToolkit/dev/basics/Linearization/#Input-derivatives)

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [June 6, 2023, 6:42pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/12 "2023-06-06T18:42:35Z")

</div>

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

Code:

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

```

---

<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 6, 2023, 6:47pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/13 "2023-06-06T18:47:03Z")

</div>

> 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

> [@ufechner7](#):
>
> Please, before making the next suggestion how to fix this issue try it yourself.

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.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [June 6, 2023, 7:34pm UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/14 "2023-06-06T19:34:25Z")

</div>

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

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!

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 9, 2023, 8:07am UTC](https://discourse.julialang.org/t/how-to-linearize-a-model/99886/15 "2023-07-09T08:07:54Z")

</div>

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:

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

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

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