# Linearization with ModelingToolkit gives different results

**URL:** <https://discourse.julialang.org/t/linearization-with-modelingtoolkit-gives-different-results/105131>\
**Category:** Modelling & Simulations\
**Tags:** question, package, modelingtoolkit\
**Created:** [October 18, 2023, 4:52pm UTC](https://discourse.julialang.org/t/linearization-with-modelingtoolkit-gives-different-results/105131 "2023-10-18T16:52:44Z")\
**Posts on this page:** 9\
**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:** [October 18, 2023, 4:52pm UTC](https://discourse.julialang.org/t/linearization-with-modelingtoolkit-gives-different-results/105131/1 "2023-10-18T16:52:45Z")

</div>

Using the methods linearize\_symbolic and linearize gives different results, even though the results should be the same. Could be a bug in my code, or in ModelingToolkit. I share the code here, even though it is still pretty long for an MWE:

```julia
using ModelingToolkit, ControlSystemsBase

@variables t
@syms ω(t) Pg(t) Pge(t) Pgc(t) Uest(t)
@variables R a b Γ ρ U J A Ku Γest λnom s

struct Equilibrium
    Γ
    u_est_bar # estimated wind speed
    a_bar
    b_bar
    pg_bar # generator power
    ω_bar # rotor speed
end

mutable struct Settings
    ρ0::Float64
    U0::Float64
    A0::Float64
    R0::Float64
    J0::Float64
    λ_nom::Float64
    Ku0::Float64
end

se::Settings = Settings(1.225, 7.0, 12469.0, 63.0, 4.047e7, 8.0, 2.25e-7)
eq = Equilibrium(0.99, 6.97822, -0.012273163568697454, 0.5344744027433362, 1132254, 0.88612)

function model_open_loop(se, eq, Γest0)
    @variables t ω(t)=eq.ω_bar λ(t) λest(t) Pr(t)=eq.pg_bar Pr_est(t)=eq.pg_bar Pg(t) 
    @variables Uest(t)=se.U0 r_ω(t) e_ω(t) e_p(t) Pgc(t)=eq.pg_bar Pge(t)=0
    @parameters J=se.J0 U=se.U0 ρ=se.ρ0 A=se.A0 Ku=se.Ku0 λnom=se.λ_nom R=se.R0 Γ=eq.Γ Γest=Γest0
    D = Differential(t)
    eqs = [D(ω) ~ (Pr - Pg) / (J * ω),
           D(Uest) ~ Ku * e_p,
           r_ω ~ Uest * λnom / R,
           e_ω ~ r_ω - ω,
           e_p ~ Pg - Pr_est + J * D(ω) * ω,
           λest ~ ω * R / Uest,
           Pr_est ~ 0.5 * ρ * A * Uest^3 * Γest * (eq.a_bar * λest + eq.b_bar),
           λ ~ ω * R/U,
           Pr ~ 0.5 * ρ * A * U^3 * Γ * (eq.a_bar * λ + eq.b_bar),
           Pg ~ Pgc + Pge]
    @named sys = ODESystem(eqs, t)
    sys, Pge, Pgc, e_ω, e_p
end

function model_open_loop_symbolic()
    @variables t ω(t) Pr(t) Pg(t) Pge(t) Pgc(t) Pr_est(t)
    @variables e_p(t) e_ω(t) r_ω(t) λ(t) λest(t) Uest(t)
    @parameters J U Ku ρ A Γ Γest a b R λnom
    D = Differential(t)
    eqs = [D(ω) ~ (Pr - Pg) / (J * ω),
           D(Uest) ~ Ku * e_p,
           r_ω ~ Uest * λnom / R,
           e_ω ~ r_ω - ω,
           e_p ~ Pg - Pr_est + J * D(ω) * ω,
           λest ~ ω * R / Uest,
           Pr_est ~ 0.5 * ρ * A * Uest^3 * Γest * (a*λest + b),
           λ ~ ω * R/U,
           Pr ~ 0.5 * ρ * A * U^3 * Γ * (a*λ + b),
           Pg ~ Pgc + Pge]
    @named sys = ODESystem(eqs, t)
    sys, Pge, Pgc, e_p, e_ω
end

function calc_lin_numeric()
    sys, Pge, Pgc, e_ω, e_p = model_open_loop(se, eq, 1.0)
    matrices, simplified_sys = ModelingToolkit.linearize(sys, [Pgc, Pge], [e_ω, e_p])
    A_, B, C, D = matrices
    A_[2]
end

function numeric(eq::Equilibrium, expr)
    expr = simplify(expr)
    expr = substitute(expr, Dict([Ku => se.Ku0, R => se.R0, ρ => se.ρ0, A => se.A0, J => se.J0, U => se.U0, 
                    λnom => se.λ_nom]))
    expr = substitute(expr, Dict([Γest => 1, Γ => eq.Γ, Pge(t) => 0, Pgc(t) => eq.pg_bar, Uest(t) => eq.u_est_bar, 
                    ω(t) => eq.ω_bar, a => eq.a_bar, b => eq.b_bar]))
end

function calc_lin_symbolic(num=true)
    sys, Pge, Pgc, e_p, e_ω = model_open_loop_symbolic()
    matrices, simplified_sys = ModelingToolkit.linearize_symbolic(sys, [Pge, Pgc], [e_ω, e_p])
    A_, B, C, D = matrices
    A_[2] = eval(Meta.parse(replace(repr(A_[2]), "Pge(t)" => "0")))
    if num
        return numeric(eq, A_[2])
    end
    A_[2]
end

println("Symbolic linearization: ", calc_lin_symbolic())
println("Numeric linearization: ", calc_lin_numeric())

```

Output:

```julia
julia> include("mwes\\mwe14.jl")
Symbolic linearization: 0.0002465406240476484
Numeric linearization: 0.0006510485677325923

```

Expected output: Both numbers should be the same.

I am hunting this bug now for 10 days without success, perhaps someone here spots it…

---

<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:** [October 18, 2023, 5:14pm UTC](https://discourse.julialang.org/t/linearization-with-modelingtoolkit-gives-different-results/105131/2 "2023-10-18T17:14:33Z")

</div>

There’s a lot going on here, from a quick look, the two functions `model_open_loop, model_open_loop_symbolic` are not identical. To find the source of the problem, it would probably be a good idea to reduce the example as much as possible, there’s too much code right now to dig through.

---

<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:** [October 18, 2023, 5:15pm UTC](https://discourse.julialang.org/t/linearization-with-modelingtoolkit-gives-different-results/105131/3 "2023-10-18T17:15:59Z")

</div>

> [@baggepinnen](#):
>
> , it would probably be a good idea to reduce the example as much as possible,

I will try to in the coming days…

---

<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:** [October 19, 2023, 2:47am UTC](https://discourse.julialang.org/t/linearization-with-modelingtoolkit-gives-different-results/105131/4 "2023-10-19T02:47:51Z")

</div>

Are you aware of the `op` keyword to `linearize` to set the operating point in which to linearize? You shouldn’t have to build the operating point into the system equations like it seems you are doing with the `eq` struct?

---

<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:** [October 19, 2023, 4:47am UTC](https://discourse.julialang.org/t/linearization-with-modelingtoolkit-gives-different-results/105131/5 "2023-10-19T04:47:38Z")

</div>

```julia
sys, Pge, Pgc, e_ω, e_p = model_open_loop(se, eq, 1.0)
matrices, simplified_sys = ModelingToolkit.linearize(sys, [Pgc, Pge], [e_ω, e_p])
sys_numeric = ss(matrices...)

matricess, _ = ModelingToolkit.linearize_symbolic(sys, [Pgc, Pge], [e_ω, e_p])
sys_symbolic = ss(ssdata(matricess)...)

defs = ModelingToolkit.defaults(simplified_sys)
u0, pars = ModelingToolkit.get_u0_p(simplified_sys, defs, defs)

using ContorlSystemsMTK
fun = Symbolics.build_function(sys_symbolic, states(simplified_sys), ModelingToolkit.parameters(simplified_sys); expression=Val{false}, force_SA=true)
sys_symbolic_numeric = fun(u0, pars)

```

```julia
julia> sys_numeric
ControlSystemsBase.StateSpace{Continuous, Float64}
A = 
 -0.007988194954103444 -0.0
  0.0006510485677325923 -0.11852648581281108
B = 
 -2.788523165896014e-8 -2.788523165896014e-8
  0.0 0.0
C = 
   -1.0 0.12698412698412698
 2893.5491899226326 -526784.3813902715
D = 
 0.0 0.0
 0.0 0.0

Continuous-time state-space model

julia> sys_symbolic_numeric
StaticStateSpace{Continuous, 2, 2, 2, SMatrix{2, 2, Float64, 4}, SMatrix{2, 2, Float64, 4}, SMatrix{2, 2, Float64, 4}, SMatrix{2, 2, Int64, 4}}
A = 
 -0.007988194954103444 0.0
  0.0006510485677325923 -0.11852648581281106
B = 
 -2.788523165896014e-8 -2.788523165896014e-8
  0.0 0.0
C = 
   -1.0 0.12698412698412698
 2893.5491899226326 -526784.3813902714
D = 
 0 0
 0 0

Continuous-time state-space model

```

---

<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:** [October 20, 2023, 10:28am UTC](https://discourse.julialang.org/t/linearization-with-modelingtoolkit-gives-different-results/105131/6 "2023-10-20T10:28:42Z")

</div>

Well, I had the hope you would have found an issue in my code, but you didn’t.

Your code:

```julia
sys, Pge, Pgc, e_ω, e_p = model_open_loop(se, eq, 1.0)
matrices, simplified_sys = ModelingToolkit.linearize(sys, [Pgc, Pge], [e_ω, e_p])
sys_numeric = ss(matrices...)

matricess, _ = ModelingToolkit.linearize_symbolic(sys, [Pgc, Pge], [e_ω, e_p])
sys_symbolic = ss(ssdata(matricess)...)

defs = ModelingToolkit.defaults(simplified_sys)
u0, pars = ModelingToolkit.get_u0_p(simplified_sys, defs, defs)

using ControlSystemsMTK
fun = Symbolics.build_function(sys_symbolic, states(simplified_sys), ModelingToolkit.parameters(simplified_sys); expression=Val{false}, force_SA=true)
sys_symbolic_numeric = fun(u0, pars)

```

For the calculation of sys\_symbolic\_numeric you are using the simplified\_sys that was obtained using numeric linearization… But what I was doing is to replace the symbolic values in the symbolic solution by the numeric parameters and the numeric steady state values…

So I still have no idea what is wrong in my code…

---

<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:** [October 20, 2023, 2:13pm UTC](https://discourse.julialang.org/t/linearization-with-modelingtoolkit-gives-different-results/105131/7 "2023-10-20T14:13:31Z")

</div>

> [@ufechner7](#):
>
> you are using the simplified\_sys that was obtained using numeric linearization

It doesn’t matter that I get the state variables from when I used the numerical linearization, the `simplified_sys` is still completely symbolic, and I use the equations from the symbolic linearization to build the function in the end.

> [@ufechner7](#):
>
> Well, I had the hope you would have found an issue in my code, but you didn’t.

I showed you an alternative way of doing the same thing, also proving that the symbolic and the numerical linearizations are equivalent.

I wish you good luck finding out what is wrong with your code, I have no idea what’s wrong with it

---

<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:** [October 20, 2023, 2:46pm UTC](https://discourse.julialang.org/t/linearization-with-modelingtoolkit-gives-different-results/105131/8 "2023-10-20T14:46:52Z")

</div>

OK, I tried to simplify my MWE by using only one model. But this fails. Code:

```julia
using ModelingToolkit, ControlSystemsBase

@variables t
@syms ω(t) Pg(t) Pge(t) Pgc(t) Uest(t)
@variables R a b Γ ρ U J A Ku Γest λnom s

struct Equilibrium
    Γ
    u_est_bar # estimated wind speed
    a_bar
    b_bar
    pg_bar # generator power
    ω_bar # rotor speed
end

mutable struct Settings
    ρ0::Float64
    U0::Float64
    A0::Float64
    R0::Float64
    J0::Float64
    λ_nom::Float64
    Ku0::Float64
end

se::Settings = Settings(1.225, 7.0, 12469.0, 63.0, 4.047e7, 8.0, 2.25e-7)
eq = Equilibrium(0.99, 6.97822, -0.012273163568697454, 0.5344744027433362, 1132254, 0.88612)
op_ = Dict([Ku => se.Ku0, R => se.R0, ρ => se.ρ0, A => se.A0, J => se.J0, U => se.U0, 
            λnom => se.λ_nom, Γest => 1, Γ => eq.Γ, Pge(t) => 0, Pgc(t) => eq.pg_bar, Uest(t) => eq.u_est_bar, 
            ω(t) => eq.ω_bar, a => eq.a_bar, b => eq.b_bar])

function model_open_loop_symbolic()
    @variables t ω(t) Pr(t) Pg(t) Pge(t) Pgc(t) Pr_est(t)
    @variables e_p(t) e_ω(t) r_ω(t) λ(t) λest(t) Uest(t)
    @parameters J U Ku ρ A Γ Γest a b R λnom
    D = Differential(t)
    eqs = [D(ω) ~ (Pr - Pg) / (J * ω),
           D(Uest) ~ Ku * e_p,
           r_ω ~ Uest * λnom / R,
           e_ω ~ r_ω - ω,
           e_p ~ Pg - Pr_est + J * D(ω) * ω,
           λest ~ ω * R / Uest,
           Pr_est ~ 0.5 * ρ * A * Uest^3 * Γest * (a*λest + b),
           λ ~ ω * R/U,
           Pr ~ 0.5 * ρ * A * U^3 * Γ * (a*λ + b),
           Pg ~ Pgc + Pge]
    @named sys = ODESystem(eqs, t)
    sys, Pge, Pgc, e_p, e_ω
end

function calc_lin_numeric()
    sys, Pge, Pgc, e_ω, e_p = model_open_loop_symbolic()
    matrices, simplified_sys = ModelingToolkit.linearize(sys, [Pgc, Pge], [e_ω, e_p], op=op_)
    A_, B, C, D = matrices
    A_[2]
end

function calc_lin_symbolic(num=true)
    sys, Pge, Pgc, e_p, e_ω = model_open_loop_symbolic()
    matrices, simplified_sys = ModelingToolkit.linearize_symbolic(sys, [Pge, Pgc], [e_ω, e_p])
    A_, B, C, D = matrices
    A_[2] = eval(Meta.parse(replace(repr(A_[2]), "Pge(t)" => "0")))
    if num; return substitute(A_[2], op_); else return A_[2] end
end

println("Symbolic linearization: ", calc_lin_symbolic())
println("Numeric linearization: ", calc_lin_numeric())

```

Error message:

```julia
julia> include("mwes/mwe14.jl")
Symbolic linearization: 0.00024654062404763535
ERROR: LoadError: ArgumentError: SymbolicUtils.BasicSymbolic[ω(t), Uest(t)] are missing from the variable map.
Stacktrace:
  [1] throw_missingvars(vars::Vector{SymbolicUtils.BasicSymbolic})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/rrbUl/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 ~/.julia/packages/ModelingToolkit/rrbUl/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 ~/.julia/packages/ModelingToolkit/rrbUl/src/variables.jl:84
  [4] get_u0_p(sys::ODESystem, u0map::Dict{Any, Any}, parammap::SciMLBase.NullParameters; use_union::Bool, tofloat::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/rrbUl/src/systems/diffeqs/abstractodesystem.jl:591
  [5] linearize(sys::ODESystem, lin_fun::ModelingToolkit.var"#132#135"{Vector{Num}, UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}, Vector{Any}, ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#533"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6c1101f6, 0x373c10a3, 0x40367bac, 0x761f4595, 0x0a686aa7), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xc9ec9000, 0x1549f32a, 0x2492f671, 0x805d41e3, 0x939f98c0), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#583#generated_observed#541"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#14586650869325627829"), Symbol("##arg#1597596539981100354"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xe078a58e, 0x8dc33125, 0x2d3218e9, 0x029a32cb, 0x52a0ec5a), Nothing}, ForwardDiff.Chunk{2}}; t::Float64, op::Dict{Num, Float64}, allow_input_derivatives::Bool, p::SciMLBase.NullParameters)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/rrbUl/src/systems/abstractsystem.jl:1518
  [6] linearize(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num}; op::Dict{Num, Float64}, t::Float64, allow_input_derivatives::Bool, zero_dummy_der::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/rrbUl/src/systems/abstractsystem.jl:1570
  [7] calc_lin_numeric()
    @ Main ~/repos/WindTurbines/mwes/mwe14.jl:53
  [8] top-level scope
    @ ~/repos/WindTurbines/mwes/mwe14.jl:67
  [9] include(fname::String)
    @ Base.MainInclude ./client.jl:478
 [10] top-level scope
    @ REPL[1]:1
in expression starting at /home/ufechner/repos/WindTurbines/mwes/mwe14.jl:67

```

Any idea?

---

<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:** [October 23, 2023, 12:48pm UTC](https://discourse.julialang.org/t/linearization-with-modelingtoolkit-gives-different-results/105131/9 "2023-10-23T12:48:53Z")

</div>

I found and fixed the issues. There where two issues: Firstly, I could not use the same model for symbolic and numeric linearization because of the error shown in the last message. This was fixed by using `@var` instead of `@syms` and to not use the notation with `function(t)` outside of the model. You need it inside of the model.

The second issue was a different operation point for the numeric and symbolic model. I initialized `Uest(t)` for `t=0` wrongly with `se.U0` instead of `eq.u_est_bar`.

For the sake of completeness the working code:

```julia
using ModelingToolkit, ControlSystemsBase

@variables t
@variables ω(t) Pg(t) Pge(t) Pgc(t) Uest(t)
@variables R a b Γ ρ U J A Ku Γest λnom s

struct Equilibrium
    Γ
    u_est_bar # estimated wind speed
    a_bar
    b_bar
    pg_bar # generator power
    ω_bar # rotor speed
end

mutable struct Settings
    ρ0::Float64
    U0::Float64
    A0::Float64
    R0::Float64
    J0::Float64
    λ_nom::Float64
    Ku0::Float64
end

se::Settings = Settings(1.225, 7.0, 12469.0, 63.0, 4.047e7, 8.0, 2.25e-7)
eq = Equilibrium(0.99, 6.97822, -0.012273163568697454, 0.5344744027433362, 1132254, 0.88612)
op_ = Dict([Ku => se.Ku0, R => se.R0, ρ => se.ρ0, A => se.A0, J => se.J0, U => se.U0, 
            λnom => se.λ_nom, Γest => 1, Γ => eq.Γ, Pge => 0, Pgc => eq.pg_bar, Uest => eq.u_est_bar, 
            ω => eq.ω_bar, a => eq.a_bar, b => eq.b_bar])

function model_open_loop_symbolic()
    @variables t ω(t) Pr(t) Pg(t) Pge(t) Pgc(t) Pr_est(t)
    @variables e_p(t) e_ω(t) r_ω(t) λ(t) λest(t) Uest(t)
    @parameters J U Ku ρ A Γ Γest a b R λnom
    D = Differential(t)
    eqs = [D(ω) ~ (Pr - Pg) / (J * ω),
           D(Uest) ~ Ku * e_p,
           r_ω ~ Uest * λnom / R,
           e_ω ~ r_ω - ω,
           e_p ~ Pg - Pr_est + J * D(ω) * ω,
           λest ~ ω * R / Uest,
           Pr_est ~ 0.5 * ρ * A * Uest^3 * Γest * (a*λest + b),
           λ ~ ω * R/U,
           Pr ~ 0.5 * ρ * A * U^3 * Γ * (a*λ + b),
           Pg ~ Pgc + Pge]
    @named sys = ODESystem(eqs, t)
    sys, Pge, Pgc, e_p, e_ω
end

function calc_lin_numeric()
    sys, Pge, Pgc, e_ω, e_p = model_open_loop_symbolic()
    matrices, simplified_sys = ModelingToolkit.linearize(sys, [Pgc, Pge], [e_ω, e_p], op=op_)
    A_, B, C, D = matrices
    A_[2]
end

function calc_lin_symbolic(num=true)
    global matrices
    sys, Pge, Pgc, e_p, e_ω = model_open_loop_symbolic()
    matrices, simplified_sys = ModelingToolkit.linearize_symbolic(sys, [Pge, Pgc], [e_ω, e_p])
    A_, B, C, D = matrices
    if num; return substitute(A_[2], op_); else return substitute(matrices.A[1], Dict([Pge => 0])) end
end

function model_open_loop(se, eq, Γest0)
    @variables t ω(t)=eq.ω_bar λ(t) λest(t) Pr(t)=eq.pg_bar Pr_est(t)=eq.pg_bar Pg(t) 
    @variables Uest(t)=eq.u_est_bar r_ω(t) e_ω(t) e_p(t) Pgc(t)=eq.pg_bar Pge(t)=0
    @parameters J=se.J0 U=se.U0 ρ=se.ρ0 A=se.A0 Ku=se.Ku0 λnom=se.λ_nom R=se.R0 Γ=eq.Γ Γest=Γest0
    D = Differential(t)
    eqs = [D(ω) ~ (Pr - Pg) / (J * ω),
           D(Uest) ~ Ku * e_p,
           r_ω ~ Uest * λnom / R,
           e_ω ~ r_ω - ω,
           e_p ~ Pg - Pr_est + J * D(ω) * ω,
           λest ~ ω * R / Uest,
           Pr_est ~ 0.5 * ρ * A * Uest^3 * Γest * (eq.a_bar * λest + eq.b_bar),
           λ ~ ω * R/U,
           Pr ~ 0.5 * ρ * A * U^3 * Γ * (eq.a_bar * λ + eq.b_bar),
           Pg ~ Pgc + Pge]
    @named sys = ODESystem(eqs, t)
    sys, Pge, Pgc, e_ω, e_p
end

function calc_lin_numeric_old()
    sys, Pge, Pgc, e_ω, e_p = model_open_loop(se, eq, 1.0)
    matrices, simplified_sys = ModelingToolkit.linearize(sys, [Pgc, Pge], [e_ω, e_p])
    A_, B, C, D = matrices
    A_[2]
end

println("Symbolic linearization: ", calc_lin_symbolic(false))
println("Symbolic linearization, numeric result: ", calc_lin_symbolic())
println("Numeric linearization: ", calc_lin_numeric())
println("Old numeric linearization: ", calc_lin_numeric_old())

```

which has the output:

```julia
julia> include("mwes\\mwe17.jl")
Symbolic linearization: (0.5A*R*(U^2)*a*Γ*ρ) / (J*ω(t)) + (-(-Pgc(t) + 0.5A*(U^3)*(b + (R*a*ω(t)) / U)*Γ*ρ)) / (J*(ω(t)^2))
Symbolic linearization, numeric result: 0.00024654062404763535
Numeric linearization: 0.00024654062404766154
Old numeric linearization: 0.00024654062404766154

```