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:
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> 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…