Put linear ODE system in matrix-vector form

Specific Domains Modelling & Simulations
1

Hi,

I am trying to put my linear ODE system in state space form.

My model is quite simple:

function TiThTe(input_funcs::Dict)
    @parameters begin
        C_i
        C_e
        C_h
        R_ie
        R_ih
        R_ea
        A_w
        A_e
        T_i_0
        T_e_0
        T_h_0
    end
    @variables begin
        T_i(t) = T_i_0
        T_e(t) = T_e_0
        T_h(t) = T_h_0
        T_a(t)
        P_h(t)
        P_s(t)
    end

    # Define time-varying inputs from the input_funcs dictionary
    @named T_a_out = TimeVaryingFunction(input_funcs[:T_a])
    @named P_h_out = TimeVaryingFunction(input_funcs[:P_h])
    @named P_s_out = TimeVaryingFunction(input_funcs[:P_s])

    systems = [T_a_out, P_h_out, P_s_out]

    # define the system of differential equations
    eqs = [
        # dTi equation
        C_i * D(T_i) ~ (1 / R_ih) * (T_h - T_i) + (1 / R_ie) * (T_e - T_i) + A_w * P_s,

        # dTe equation
        C_e * D(T_e) ~ (1 / R_ie) * (T_i - T_e) + (1 / R_ea) * (T_a - T_e) + A_e * P_s,

        # dTh equation
        C_h * D(T_h) ~ (1 / R_ih) * (T_i - T_h) + P_h
    ]

    # Define the inputs (mapping them to their respective time-varying functions)
    inputs = [
        T_a ~ T_a_out.output.u,
        P_h ~ P_h_out.output.u,
        P_s ~ P_s_out.output.u
    ]
    append!(eqs, inputs)

    @named sys = ODESystem(eqs, t; systems=systems)
    return sys
end

My input functions are given by DataInterpolations.jl.

I try to form the matrices in this way, defining my inputs and outputs:

# inputs
@variables t T_a(t) P_h(t) P_s(t) 
inputs = [T_a, P_h, P_s]

# outputs
@variables T_i(t) T_e(t) T_h(t)
outputs= [T_i, T_e, T_h]

# linearize
ss_matrices, ssys = ModelingToolkit.linearize_symbolic(sys, inputs, outputs)

Unfortunately this method complains about missing variables which I dont want to see in my linearization:

ExtraEquationsSystemException: The system is unbalanced. There are 6 highest order derivative variables and 9 equations.
More equations than variables, here are the potential extra equation(s):
0 ~ -T_a_out₊output₊u(t) + DataInterpolations.AkimaInterpolation{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}([8.167, 8.187, 8.077, 7.988, 7.98, 7.814, 7.784, 7.531, 7.638, 8.057

But I tried adding them anyway to see what result I get:

@variables t T_a(t) P_h(t) P_s(t) T_a_out₊output₊u(t) P_h_out₊output₊u(t) P_s_out₊output₊u(t)
inputs = [T_a, P_h, P_s, T_a_out₊output₊u,  P_h_out₊output₊u, P_s_out₊output₊u]

This throws a boundserror in AbstractSystem.jl:

BoundsError: attempt to access 9×3 Matrix{Num} at index [1:3, 4:9]
Stacktrace:
[1] throw_boundserror(A::Matrix{Num}, I::Tuple{UnitRange{Int64}, UnitRange{Int64}})
@ Base .\abstractarray.jl:737
[2] checkbounds
@ .\abstractarray.jl:702 [inlined]
[3] _getindex
@ .\multidimensional.jl:888 [inlined]
[4] getindex(::Matrix{Num}, ::UnitRange{Int64}, ::UnitRange{Int64})
@ Base .\abstractarray.jl:1291
[5] 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\lange.julia\packages\ModelingToolkit\Vsl3C\src\systems\abstractsystem.jl:2496
[6] linearize_symbolic(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num})
@ ModelingToolkit C:\Users\lange.julia\packages\ModelingToolkit\Vsl3C\src\systems\abstractsystem.jl:2472

Could somebody advise me what to do?

2

The first error, extra equations, is expected. This video explains why it appears and how to change your model to make it work

Tldr is that you cannot have the DataInterpolations connected while linearizing, the inputs must be unconnected.

This other video shared some additional tools for linear analysis you might find interesting

3

Here’s a working example

using ModelingToolkit, ModelingToolkitStandardLibrary
import ModelingToolkit: t_nounits as t, D_nounits as D
function TiThTe(; name)
    @parameters begin
        C_i
        C_e
        C_h
        R_ie
        R_ih
        R_ea
        A_w
        A_e
        T_i_0
        T_e_0
        T_h_0
    end
    @variables begin
        T_i(t) = T_i_0
        T_e(t) = T_e_0
        T_h(t) = T_h_0
        T_a(t)
        P_h(t)
        P_s(t)
    end

    # define the system of differential equations
    eqs = [
        # dTi equation
        C_i * D(T_i) ~ (1 / R_ih) * (T_h - T_i) + (1 / R_ie) * (T_e - T_i) + A_w * P_s,

        # dTe equation
        C_e * D(T_e) ~ (1 / R_ie) * (T_i - T_e) + (1 / R_ea) * (T_a - T_e) + A_e * P_s,

        # dTh equation
        C_h * D(T_h) ~ (1 / R_ih) * (T_i - T_h) + P_h
    ]

    ODESystem(eqs, t; name)
end

@named sys = TiThTe()
sys = complete(sys)

# inputs
inputs = [sys.T_a, sys.P_h, sys.P_s]

# outputs
outputs= [sys.T_i, sys.T_e, sys.T_h]

# linearize
ss_matrices, ssys = ModelingToolkit.linearize_symbolic(sys, inputs, outputs)

julia> ss_matrices.A
3×3 Matrix{Num}:
 (-1 / R_ie + -1 / R_ih) / C_i                 1 / (C_i*R_ie)   1 / (C_i*R_ih)
                1 / (C_e*R_ie)  (-1 / R_ea + -1 / R_ie) / C_e                0
                1 / (C_h*R_ih)                              0  -1 / (C_h*R_ih)
4

Thanks a lot, that works great :slight_smile:

Is there perhaps a way to apply this to an already existing model, basically remove the connections? So I don’t have to keep two model definitions.

5

Yes, this is accomplished by the analysis points described in the videos. You can also define the model without connections, and a thin outer model that adds the input connections when they are needed. This is a more extensible approach anyways so I recommend following this approach

6

Yes I also like this approach. I think it can be done using the new MTK v9 syntax for datainterpolations? Like here: https://github.com/SciML/ModelingToolkit.jl/issues/2823#issuecomment-2366342737

(last time I checked there weren’t any docs on how to do this exactly and I couldn’t figure out the syntax myself for my model)