MTK... linearization of zero dynamics

[BLI]

[What is the trick of hiding/folding/collapsing code blocks?]

Consider the well-known 4-tank model by K.H. Johansson:

which has been designed to illustrate multivariable zero dynamics. The system has 2 control inputs (v_1, v_2), and two outputs (y_1, y_2); the outputs are (scaled) measures of the lower levels (h_1, h_2), while the upper levels (h_3, h_4) can be estimated. The actuators (v_1, v_2) affect all four tanks via some split-range valves.

The following code implements the model with inputs, initial states, and parameters according to Johansson’s P- case (i.e., with stable zero dynamics):

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as Dt
using DifferentialEquations
using ControlSystems, ControlSystemsMTK
using Plots, LaTeXStrings
#
# Input functions
# -- registering inputs
@register_symbolic u_u1(t)
@register_symbolic u_u2(t)
#
u_v_1(t) = 3.0
u_v_2(t) = 3.0
@register_symbolic u_r1(t)
@register_symbolic u_r2(t)
#
u_r_1(t) = 12.4*0.5
u_r_2(t) = 12.7*0.5
#
# -- Model
@mtkmodel FourTank begin
    # Ref: Johansson, K.H. (2000): "The Quadruple-Tank Process: A Multivariable 
    #               Laboratory Process with an Adjustable Zero", 
    #               IEEE Trans. Ctrl. Tech., Vol. 8, No. 3, May 2000
    # Model parameters
    @parameters begin
        # Constant 
        g = 981,            [description = "Acceleration of gravity, cm/s^2"]
        # Parameters
        # -- tank parameters
        A_1 = 28,           [description = "Tank 1, cross-section area, cm^2"]
        A_2 = 32,           [description = "Tank 2, cross-section area, cm^2"]
        A_3 = 28,           [description = "Tank 3, cross-section area, cm^2"]
        A_4 = 32,           [description = "Tank 4, cross-section area, cm^2"]
        a_1 = 0.071,        [description = "Tank 1, cross-section area of outlet hole, cm^2"]
        a_2 = 0.057,        [description = "Tank 2, cross-section area of outlet hole, cm^2"]
        a_3 = 0.071,        [description = "Tank 3, cross-section area of outlet hole, cm^2"]
        a_4 = 0.057,        [description = "Tank 4, cross-section area of outlet hole, cm^2"]
        # -- valve parameters, LHP zeros
        γ_1 = 0.70,         [description = "Valve 1 split setting, -"]
        γ_2 = 0.60,         [description = "Valve 2 split setting, -"]
        k_1 = 3.33,         [description = "Valve 1 gain, cm3/V.s"]
        k_2 = 3.35,         [description = "Valve 2 gain, cm3/V.s"]
        # -- sensor parameter
        k_c = 0.5,          [description = "Sensor gain, V/cm"]
    end
    #
    # Model variables, with initial values needed
    @variables begin
        # Variables
        # -- states
        h_1(t)=12.4,        [description = "Tank 1 level, cm"]
        h_2(t)=12.7,        [description = "Tank 2 level, cm"]
        h_3(t)=1.8,         [description = "Tank 3 level, cm"]
        h_4(t)=1.4,         [description = "Tank 4 level, cm"]
        # -- sensor signals
        y_1(t),             [description = "Output 1: level in Tank 1, cm"]
        y_2(t),             [description = "Output 2: level in Tank 2, cm"]
        # -- actuator signals
        v_1(t),             [description = "Valve 1 signal, V"]
        v_2(t),             [description = "Valve 2 signal, V"]
    end
    #
    # Model equations
    @equations begin
        # ODEs
        Dt(h_1) ~ -a_1/A_1*sqrt(2g*h_1) + a_3/A_1*sqrt(2g*h_3) + γ_1*k_1/A_1*v_1
        Dt(h_2) ~ -a_2/A_2*sqrt(2g*h_2) + a_4/A_2*sqrt(2g*h_4) + γ_2*k_2/A_2*v_2
        Dt(h_3) ~ -a_3/A_3*sqrt(2g*h_3) + (1-γ_2)*k_2/A_3*v_2
        Dt(h_4) ~ -a_4/A_4*sqrt(2g*h_4) + (1-γ_1)*k_1/A_4*v_1
        # Outputs
        y_1 ~ k_c*h_1
        y_2 ~ k_c*h_2
        # Inputs
        v_1 ~ u_v_1(t)
        v_2 ~ u_v_2(t)
    end
end

@mtkbuild t4 = FourTank()
tspan = (0,60*8)
prob = ODEProblem(t4, [], tspan)
sol = solve(prob, reltol=1e-4)

The results indicate that the specified operating points in the referenced paper are not 100% in steady state:

Next, I can do linearization of the model with the trick of creating a new model [FourTank_noi] where I have no inputs, achieved by commenting out the two last equations in FourTank. This works fine with the latest update of ModelingToolkit.

---

Next, I want to study the zero dynamics. To do so, I relax on specifying the inputs v_1 and v_2, and instead specify the outputs y_1 and y_2:

@mtkmodel FourTank_z begin
    #
    # Model parameters
    @parameters begin
        # Constant 
        g = 981,            [description = "Acceleration of gravity, cm/s^2"]
        # Parameters
        # -- tank parameters
        A_1 = 28,           [description = "Tank 1, cross-section area, cm^2"]
        A_2 = 32,           [description = "Tank 2, cross-section area, cm^2"]
        A_3 = 28,           [description = "Tank 3, cross-section area, cm^2"]
        A_4 = 32,           [description = "Tank 4, cross-section area, cm^2"]
        a_1 = 0.071,        [description = "Tank 1, cross-section area of outlet hole, cm^2"]
        a_2 = 0.057,        [description = "Tank 2, cross-section area of outlet hole, cm^2"]
        a_3 = 0.071,        [description = "Tank 3, cross-section area of outlet hole, cm^2"]
        a_4 = 0.057,        [description = "Tank 4, cross-section area of outlet hole, cm^2"]
        # -- valve parameters, LHP zeros
        γ_1 = 0.70,         [description = "Valve 1 split setting, -"]
        γ_2 = 0.60,         [description = "Valve 2 split setting, -"]
        k_1 = 3.33,         [description = "Valve 1 gain, cm3/V.s"]
        k_2 = 3.35,         [description = "Valve 2 gain, cm3/V.s"]
        # -- sensor parameter
        k_c = 0.5,          [description = "Sensor gain, V/cm"]
    end
    #
    # Model variables, with initial values needed
    @variables begin
        # Variables
        # -- states
        h_1(t)=12.4,        [description = "Tank 1 level, cm"]
        h_2(t)=12.7,        [description = "Tank 2 level, cm"]
        h_3(t)=1.8,         [description = "Tank 3 level, cm"]
        h_4(t)=1.4,         [description = "Tank 4 level, cm"]
        # -- sensor signals
        y_1(t),             [description = "Output 1: level in Tank 1, cm"]
        y_2(t),             [description = "Output 2: level in Tank 2, cm"]
        # -- actuator signals
        v_1(t),             [description = "Valve 1 signal, V"]
        v_2(t),             [description = "Valve 2 signal, V"]
    end
    #
    # Model equations
    @equations begin
        # ODEs
        Dt(h_1) ~ -a_1/A_1*sqrt(2g*h_1) + a_3/A_1*sqrt(2g*h_3) + γ_1*k_1/A_1*v_1
        Dt(h_2) ~ -a_2/A_2*sqrt(2g*h_2) + a_4/A_2*sqrt(2g*h_4) + γ_2*k_2/A_2*v_2
        Dt(h_3) ~ -a_3/A_3*sqrt(2g*h_3) + (1-γ_2)*k_2/A_3*v_2
        Dt(h_4) ~ -a_4/A_4*sqrt(2g*h_4) + (1-γ_1)*k_1/A_4*v_1
        # Outputs
        y_1 ~ k_c*h_1
        y_2 ~ k_c*h_2
        # Inputs
        y_1 ~ u_r_1(t)
        y_2 ~ u_r_2(t)
    end
end

@mtkbuild t4z = FourTank_z()
prob_z = ODEProblem(t4z, [], tspan)

This gives a warning that the system is over specified. The reason is (of course) that I have specified initial states for all four levels. If I remove the initial states for h_1 and h_2, the warning disappears.

Simulating the zero dynamics gives

which is as expected: outputs y_1 and y_2 are specified to be constant (and half the value of h_1 and h_2), which only leaves h_3 and h_4 to move. Because I chose the operating point with stable zero dynamics, the system should be stable.

Next, I want to linearize the zero dynamics so that I can find the eigenvalue(s) of the linearized model: these (this) should equal the (transmission) zeros of the system.

With the standard trick of commenting out the “input” specification (which in this case is y_1 and y_2), I get:

mtkmodel FourTank_z_noi begin
    #
    # Model parameters
    @parameters begin
        # Constant 
        g = 981,            [description = "Acceleration of gravity, cm/s^2"]
        # Parameters
        # -- tank parameters
        A_1 = 28,           [description = "Tank 1, cross-section area, cm^2"]
        A_2 = 32,           [description = "Tank 2, cross-section area, cm^2"]
        A_3 = 28,           [description = "Tank 3, cross-section area, cm^2"]
        A_4 = 32,           [description = "Tank 4, cross-section area, cm^2"]
        a_1 = 0.071,        [description = "Tank 1, cross-section area of outlet hole, cm^2"]
        a_2 = 0.057,        [description = "Tank 2, cross-section area of outlet hole, cm^2"]
        a_3 = 0.071,        [description = "Tank 3, cross-section area of outlet hole, cm^2"]
        a_4 = 0.057,        [description = "Tank 4, cross-section area of outlet hole, cm^2"]
        # -- valve parameters, LHP zeros
        γ_1 = 0.70,         [description = "Valve 1 split setting, -"]
        γ_2 = 0.60,         [description = "Valve 2 split setting, -"]
        k_1 = 3.33,         [description = "Valve 1 gain, cm3/V.s"]
        k_2 = 3.35,         [description = "Valve 2 gain, cm3/V.s"]
        # -- sensor parameter
        k_c = 0.5,          [description = "Sensor gain, V/cm"]
    end
    #
    # Model variables, with initial values needed
    @variables begin
        # Variables
        # -- states
        h_1(t)=12.4,        [description = "Tank 1 level, cm"]
        h_2(t)=12.7,        [description = "Tank 2 level, cm"]
        h_3(t)=1.8,         [description = "Tank 3 level, cm"]
        h_4(t)=1.4,         [description = "Tank 4 level, cm"]
        # -- sensor signals
        y_1(t),             [description = "Output 1: level in Tank 1, cm"]
        y_2(t),             [description = "Output 2: level in Tank 2, cm"]
        # -- actuator signals
        v_1(t),             [description = "Valve 1 signal, V"]
        v_2(t),             [description = "Valve 2 signal, V"]
    end
    #
    # Model equations
    @equations begin
        # ODEs
        Dt(h_1) ~ -a_1/A_1*sqrt(2g*h_1) + a_3/A_1*sqrt(2g*h_3) + γ_1*k_1/A_1*v_1
        Dt(h_2) ~ -a_2/A_2*sqrt(2g*h_2) + a_4/A_2*sqrt(2g*h_4) + γ_2*k_2/A_2*v_2
        Dt(h_3) ~ -a_3/A_3*sqrt(2g*h_3) + (1-γ_2)*k_2/A_3*v_2
        Dt(h_4) ~ -a_4/A_4*sqrt(2g*h_4) + (1-γ_1)*k_1/A_4*v_1
        # Outputs
        y_1 ~ k_c*h_1
        y_2 ~ k_c*h_2
        # Inputs
        #y_1 ~ u_r_1(t)
        #y_2 ~ u_r_2(t)
    end
end

@named t4z_noi = FourTank_z_noi()
t4z_noi = complete(t4z_noi)

leading to

OK. Now, the problems are popping up… I start with attempting symbolic linearization [this works fine for the system from valve signal inputs v_j to outputs y_j, but not for the zero dynamics]:

julia> mats,t4z_ = ModelingToolkit.linearize_symbolic(t4z_noi, [t4z_noi.y_1, t4z_noi.y_2], [t4z_noi.h_3, t4z_noi.h_4])
BoundsError: attempt to access 1-element Vector{Any} at index [2]

Stacktrace:
 [1] getindex(A::Vector{Any}, i1::Int64)
   @ Base .\essentials.jl:13
 [2] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Real}, simplify::Bool; occurrences::SymbolicUtils.BasicSymbolic{Real}) (repeats 5 times)
   @ Symbolics C:\Users\Bernt_Lie\.julia\packages\Symbolics\6WqId\src\diff.jl:257
 [3] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Real}, simplify::Bool; occurrences::Nothing)
   @ Symbolics C:\Users\Bernt_Lie\.julia\packages\Symbolics\6WqId\src\diff.jl:257
 [4] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Real}, simplify::Bool)
   @ Symbolics C:\Users\Bernt_Lie\.julia\packages\Symbolics\6WqId\src\diff.jl:183
 [5] jacobian(ops::Vector{SymbolicUtils.BasicSymbolic{Real}}, vars::Vector{Num}; simplify::Bool, scalarize::Bool)
   @ Symbolics C:\Users\Bernt_Lie\.julia\packages\Symbolics\6WqId\src\diff.jl:470
 [6] jacobian(ops::Vector{SymbolicUtils.BasicSymbolic{Real}}, vars::Vector{Num})
   @ Symbolics C:\Users\Bernt_Lie\.julia\packages\Symbolics\6WqId\src\diff.jl:465
 [7] linearize_symbolic(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num}; simplify::Bool, allow_input_derivatives::Bool, eval_expression::Bool, eval_module::Module, kwargs::@Kwargs{})
   @ ModelingToolkit C:\Users\Bernt_Lie\.julia\packages\ModelingToolkit\UXr3S\src\systems\abstractsystem.jl:2454
 [8] linearize_symbolic(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num})
   @ ModelingToolkit C:\Users\Bernt_Lie\.julia\packages\ModelingToolkit\UXr3S\src\systems\abstractsystem.jl:2434
 [9] top-level scope
   @ c:\Users\Bernt_Lie\...\jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_X36sZmlsZQ==.jl:1

Next, I’m trying numeric/AD linearization:

julia> mats,t4z_ = ModelingToolkit.linearize(t4z_noi, [t4z_noi.y_1, t4z_noi.y_2], [t4z_noi.h_3, t4z_noi.h_4];
                op=Dict(t4z_noi.y_1 => 6.2, t4z_noi.y_2 => 6.35))
ArgumentError: Cannot convert Sym to Float64 since Sym is symbolic and Float64 is concrete. Use `substitute` to replace the symbolic unwraps.

Stacktrace:
  [1] Float64(x::SymbolicUtils.BasicSymbolic{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{true, ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#815"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x68d7fc22, 0x9cd0330d, 0x51936f95, 0xda8a0357, 0xa6ac43de), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x077ceea2, 0x254bf500, 0x74c8eda1, 0x54c835c5, 0x5d3dd8d7), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Float64, MTKParameters{Vector{Float64}, Tuple{}, Tuple{}, Tuple{}}}, Float64}, Float64, 2}})
    @ Symbolics C:\Users\Bernt_Lie\.julia\packages\Symbolics\6WqId\src\Symbolics.jl:171
  [2] convert
    @ C:\Users\Bernt_Lie\.julia\packages\ForwardDiff\PcZ48\src\dual.jl:435 [inlined]
  [3] setindex!(A::Vector{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{true, ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#815"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x68d7fc22, 0x9cd0330d, 0x51936f95, 0xda8a0357, 0xa6ac43de), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x077ceea2, 0x254bf500, 0x74c8eda1, 0x54c835c5, 0x5d3dd8d7), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Float64, MTKParameters{Vector{Float64}, Tuple{}, Tuple{}, Tuple{}}}, Float64}, Float64, 2}}, x::SymbolicUtils.BasicSymbolic{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{true, ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#815"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x68d7fc22, 0x9cd0330d, 0x51936f95, 0xda8a0357, 0xa6ac43de), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x077ceea2, 0x254bf500, 0x74c8eda1, 0x54c835c5, 0x5d3dd8d7), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Float64, MTKParameters{Vector{Float64}, Tuple{}, Tuple{}, Tuple{}}}, Float64}, Float64, 2}}, i1::Int64)
    @ Base .\array.jl:1021
  [4] macro expansion
    @ C:\Users\Bernt_Lie\.julia\packages\SymbolicUtils\5UmB7\src\code.jl:430 [inlined]
  [5] macro expansion
    @ C:\Users\Bernt_Lie\.julia\packages\Symbolics\6WqId\src\build_function.jl:546 [inlined]
  [6] macro expansion
    @ C:\Users\Bernt_Lie\.julia\packages\SymbolicUtils\5UmB7\src\code.jl:387 [inlined]
  [7] macro expansion
    @ C:\Users\Bernt_Lie\.julia\packages\RuntimeGeneratedFunctions\M9ZX8\src\RuntimeGeneratedFunctions.jl:163 [inlined]
  [8] macro expansion
    @ .\none:0 [inlined]
  [9] generated_callfunc
    @ .\none:0 [inlined]
 [10] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x077ceea2, 0x254bf500, 0x74c8eda1, 0x54c835c5, 0x5d3dd8d7), Nothing})(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{true, ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#815"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x68d7fc22, 0x9cd0330d, 0x51936f95, 0xda8a0357, 0xa6ac43de), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x077ceea2, 0x254bf500, 0x74c8eda1, 0x54c835c5, 0x5d3dd8d7), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Float64, MTKParameters{Vector{Float64}, Tuple{}, Tuple{}, Tuple{}}}, Float64}, Float64, 2}}, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{true, ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#815"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x68d7fc22, 0x9cd0330d, 0x51936f95, 0xda8a0357, 0xa6ac43de), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x077ceea2, 0x254bf500, 0x74c8eda1, 0x54c835c5, 0x5d3dd8d7), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Float64, MTKParameters{Vector{Float64}, Tuple{}, Tuple{}, Tuple{}}}, Float64}, Float64, 2}}, ::Vector{Float64}, ::Float64)
    @ RuntimeGeneratedFunctions C:\Users\Bernt_Lie\.julia\packages\RuntimeGeneratedFunctions\M9ZX8\src\RuntimeGeneratedFunctions.jl:150
 [11] f
    @ C:\Users\Bernt_Lie\.julia\packages\ModelingToolkit\UXr3S\src\systems\diffeqs\abstractodesystem.jl:351 [inlined]
...
    @ ModelingToolkit C:\Users\Bernt_Lie\.julia\packages\ModelingToolkit\UXr3S\src\systems\abstractsystem.jl:2657
 [20] linearize(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num}; op::Dict{Num, Float64}, t::Float64, allow_input_derivatives::Bool, zero_dummy_der::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\Bernt_Lie\.julia\packages\ModelingToolkit\UXr3S\src\systems\abstractsystem.jl:2708
 [21] top-level scope
    @ c:\Users\Bernt_Lie\...\jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_X24sZmlsZQ==.jl:1
Output is truncated. View as a scrollable element or open in a text editor. Adjust cell output settings...

Questions:

1. Do I do something wrong?
2. Or is there a limitation in MTK?

[baggepinnen]

Looks like a bug in MTK, please open an issue (preferably with a single code block that reproduces the error), thanks

[baggepinnen]

Here’s an explanation and a workaround

MTK.linearize_symbolic and MTK.linearize fails for "simple" system - different error messages

I’m not surprised that this happens, to invert a strictly proper model you will need some order of derivatives of the inputs. When simulating, this may work out if the input function is differentiable, but when linearizing there is nothing to differentiate, so the compiler has to be augmented with support for this use case, i.e., by

• Erroring with a useful error message.
• Or better, indicate that the user must also provide the derivatives of inputs as another set of inputs

What you can do right now is to introduce new input variables yd_1 and define y_1 as D(y_1) ~ yd_1

The issue comment includes a working example

---

PS. being from Lund, I have endured many hours listening to the sound of those pumps

[BLI]

Nice!

Of course, this trick introduces two artificial integrators in the zero dynamics, but that is probably fine. Here is what I get:

Sorry for reminding you about the pumps

[baggepinnen]

@BLI

Of course, this trick introduces two artificial integrators in the zero dynamics, but that is probably fine.

That’s unavoidable if you want to represent the model on statespace form The reason is the same as inv(sys) not working for a strictly proper system sys. IIRC, dymola and possibly also open-modelica has a slightly nicer error message in situations like this, indicating that input derivatives are required.

[baggepinnen]

Those two extra integrator state-variables are are required to indicate the operating point of the input, since now the new “input” is the derivative of the input.

[BLI]

…I more or less figured out things. But some simple questions about MTK linearize vs. named_ss…

mats, hps_ = linearize(hps_LIO, [hps_LIO.u_v], [hps_LIO.w]; op=op_pt, allow_input_derivatives=true)
sys = ss(mats...)

This gives a system sys with extra columns in the B and D matrices for input derivative.

sys_c = named_ss(hps_LIO, [hps_LIO.u_v], [hps_LIO.w]; op=op_pt, allow_input_derivatives=true)

sys_c gives matrices without input derivative. So the algorithm probably uses linearize as a backbone, introduces u as an extra state, sets up the system from input \dot{u} via states (x,u) to output y, and then removes the integrator.

In my case, sys_c gives the following transfer function:

julia> zpk(sys_c)
TransferFunction{Continuous, ControlSystemsBase.SisoZpk{Float64, ComplexF64}}
                               (1.0s - 2.1160129025661592e-9)(1.0s - 0.5250534353649084)
99.58081303429545-------------------------------------------------------------------------------------
                 (1.0s + 0.4409211191696938)(1.0s + 0.14243813925501536)(1.0s - 2.2168795262579917e-9)

Continuous-time transfer function model

So – a little bit ugly with some stray numerics – two terms which are approximately s that really should cancel out.

I try to do a minimal realization to see if this fixes things, and get:

julia> zpk(minrel(sys_c))
TransferFunction{Continuous, ControlSystemsBase.SisoZpk{Float64, ComplexF64}}
                            1.0
99.58081303429545-------------------------
                 1.0s + 1.1084126936887508

Continuous-time transfer function model

which is weirdly wrong.

• Is there a way to make minrel do a better job?

---

Finally, I also played around with the original system from linearize:

A_1 = [sys.A sys.B[:,1]; zeros(1,4) 0]
B_1 = [sys.B[:,2]; 1]
C_1 = [sys.C sys.D[:,1]]
D_1 = sys.D[:,2]
sys_1 = ss(A_1,B_1,C_1,D_1)

and then

julia> zpk(sys_1)
TransferFunction{Continuous, ControlSystemsBase.SisoZpk{Float64, ComplexF64}}
                                  (1.0s - 0.5250534374809216)(1.0s - 4.19885971173042e-16)
99.58081303429545------------------------------------------------------------------------------------------
                 (1.0s + 0.44092112239911896)(1.0s + 0.142438133808711)(1.0s)(1.0s - 2.216063822164251e-16)

Continuous-time transfer function model

The same types of stray numerics as from named_ss + the artificial pole due to using \dot{u} as input.

OK – let me try to get rid of the artificial pole by multiplying with tf("s").

julia> zpk(sys_1)*tf("s")
TransferFunction{Continuous, ControlSystemsBase.SisoZpk{Float64, ComplexF64}}
                               1.0s - 0.5250534374809217
99.58081303429545------------------------------------------------------
                 (1.0s + 0.44092112239911896)(1.0s + 0.142438133808711)

Continuous-time transfer function model

This is interesting!! Now, the stray numerics is gone, and I have found the transfer function that I really searched for, without using minrel or anything…

I’m curious:

• How does named_ss get rid of the artificial pole? Probably not by multiplying with tf("s")?
• Is the idea of multiplying with tf("s") a good way to do it, or are there possible numeric problems with it?
• What sort of numeric fixing takes place when I do zpk(sys_1)*tf("s")?

[baggepinnen]

@BLI

So the algorithm probably uses linearize as a backbone

It does, which can be readily verified by reading the docstring of named_ss or inspecting the source.

@BLI

Is the idea of multiplying with tf("s") a good way to do it, or are there possible numeric problems with it?

Transfer-function representations are always a bit iffy since they use polynomials. If the system that is differentiated is strictly proper (D matrix is zero), you can achieve the multiplication by s on statespace form by postmultiplying C with A.

@BLI

What sort of numeric fixing takes place when I do zpk(sys_1)*tf("s")?

I don’t know what system you are using so I cannot easily figure out. You have not included any code to create hps_LIO.

@BLI

Is there a way to make minrel do a better job?

I don’t even know what system type you call minreal on, a statespace object, transfer function of polynomials or a ZPK model? They all have different implementations of minreal.

In general, many of your questions, including those from the deleted post above, can be answered by reading the docstring of the function/method your question is about. For example, the aforementioned docstring of named_ss will lead to to information about how the simplification is performed when input derivatives are present. The docstring for the statespace version of minreal contains a lot of customization of numeric tolerances etc. that can be used to tweak the result.

I would guess that your system is borderline uncontrollable, you can check this by checking the minimum singular values indicated by

julia> controllability((sminreal(sys_c)))
(iscontrollable = true, ranks = [3, 3, 3], sigma_min = [0.081129755005732, 0.018320012744819, 1.0238394749593724e-11])

if some of those are very small, the system has almost uncontrollable modes which may be removed by minreal.

If you have a poorly balanced statespace model, e.g.,

B = [
 -0.5963028955165542
 -4.793145772316456
 99.58081303429545
  1.0428129567742857e9]

C = [
 0.0  0.0  1.0  0.0]

you may also improve numerics by calling balance_statespace:

julia> sys_c
StateSpace{Continuous, Float64}
A = 
 -0.7340579610440452   0.03618086493757771  0.00020843434232478795  -0.0
 -4.793295206811445    0.0                  0.0                      0.0
  0.0                 26.158882437925314    0.1506987048362153       0.0
  0.0                  2.739365216029021e8  1.5781208968256433e6     0.0
B = 
 -0.5963028955165542
 -4.793145772316456
 99.58081303429545
  1.0428129567742857e9
C = 
 0.0  0.0  1.0  0.0
D = 
 0.0

Continuous-time state-space model

julia> balance_statespace(sys_c)[1]
StateSpace{Continuous, Float64}
A = 
 -0.7340579610440452  0.14472345975031084  0.006669898954393214  0.0
 -1.1983238017028612  0.0                  0.0                   0.0
  0.0                 3.2698603047406642   0.1506987048362153    0.0
  0.0                 4.081972265286775    0.1881266709358267    0.0
B = 
 -1.1926057910331085
 -2.396572886158228
  6.223800814643465
  7.7695619819632595
C = 
 0.0  0.0  16.0  0.0
D = 
 0.0

Continuous-time state-space model

minreal will perform much better on the latter system here. You can also try calling sminreal before minreal, if sminreal removes any dimensions you can also improve the numerics.