Linearization with ModelingToolkit gives different results

Specific Domains Modelling & Simulations
, ,
9

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:

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