Putting linear ODESystem in matrix-vector form

[langestefan]

Hi,

I have an ODESystem, defined in modelingtoolkit.jl by these equations:

\begin{align} \frac{\mathrm{d} Ci_{+}v\left( t \right)}{\mathrm{d}t} =& \frac{\frac{ - Ci_{+}v\left( t \right) + Ch_{+}v\left( t \right)}{Rih_{+}R} + \frac{Ce_{+}v\left( t \right) - Ci_{+}v\left( t \right)}{Rie_{+}R}}{Ci_{+}C} \\ \frac{\mathrm{d} Ce_{+}v\left( t \right)}{\mathrm{d}t} =& \frac{\frac{Ta_{in_{+}k} - Ce_{+}v\left( t \right)}{Rea_{+}R} + \frac{ - Ce_{+}v\left( t \right) + Ci_{+}v\left( t \right)}{Rie_{+}R}}{Ce_{+}C} \\ \frac{\mathrm{d} Ch_{+}v\left( t \right)}{\mathrm{d}t} =& \frac{{\Phi}cv\left( t \right) + {\Phi}hp\left( t \right) + \frac{Ci_{+}v\left( t \right) - Ch_{+}v\left( t \right)}{Rih_{+}R}}{Ch_{+}C} \end{align}

I would like to put this into a linear state-space representation, or matrix-vector form of the shape \dot{v} = Av + Bu where v is a vector of voltages and u is a vector of \Phi's. It’s not hard to do this by hand but it would be convenient if can be done programmatically.

Is there a way to do this in julia? I found ControlSystemsMTK.jl but I was unsuccesful in getting it to work.

[baggepinnen]

What error did you get when using ControlSystemsMTK? It should work

[langestefan]

I tried the following Symbolic Linearization and From ModelingToolkit to ControlSystems

using ControlSystemsMTK, ControlSystemsBase, RobustAndOptimalControl

@named model = ODESystem(eqs, t, systems = components)
sys = structural_simplify(model, allow_parameter=true)
prob = ODEProblem(sys, Pair[], (0, 10.0))

# inputs = Φhp, Φcv
# outputs = Ci.v, Ce.v, Ch.v
RobustAndOptimalControl.named_ss(model, [Φhp_ctrl.output.u, Φcv_ctrl.output.u], [Ci.v, Ce.v, Ch.v])

This will give:
BoundsError: attempt to access 63-element Vector{Vector{Int64}} at index [64]

With stacktrace:

Summary

Stacktrace: [inlined] [2] 𝑑neighbors @ [C:\Users\lange.julia\packages\ModelingToolkit\pmsNn\src\bipartite_graph.jl:369](file:///C:/Users/lange/.julia/packages/ModelingToolkit/pmsNn/src/bipartite_graph.jl:369) [inlined] [3] 𝑑neighbors @ [C:\Users\lange.julia\packages\ModelingToolkit\pmsNn\src\bipartite_graph.jl:368](file:///C:/Users/lange/.julia/packages/ModelingToolkit/pmsNn/src/bipartite_graph.jl:368) [inlined] [4] check_consistency(state::TearingState{ODESystem}, orig_inputs::Set{Any}) @ ModelingToolkit.StructuralTransformations [C:\Users\lange.julia\packages\ModelingToolkit\pmsNn\src\structural_transformation\utils.jl:96](file:///C:/Users/lange/.julia/packages/ModelingToolkit/pmsNn/src/structural_transformation/utils.jl:96) [5] _structural_simplify!(state::TearingState{ODESystem}, io::Tuple{Vector{Num}, Vector{Num}}; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, dummy_derivative::Bool, kwargs::@Kwargs{}) @ ModelingToolkit [C:\Users\lange.julia\packages\ModelingToolkit\pmsNn\src\systems\systemstructure.jl:689](file:///C:/Users/lange/.julia/packages/ModelingToolkit/pmsNn/src/systems/systemstructure.jl:689) [6] _structural_simplify! @ [C:\Users\lange.julia\packages\ModelingToolkit\pmsNn\src\systems\systemstructure.jl:672](file:///C:/Users/lange/.julia/packages/ModelingToolkit/pmsNn/src/systems/systemstructure.jl:672) [inlined] [7] structural_simplify!(state::TearingState{ODESystem}, io::Tuple{Vector{Num}, Vector{Num}}; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{}) @ ModelingToolkit [C:\Users\lange.julia\packages\ModelingToolkit\pmsNn\src\systems\systemstructure.jl:635](file:///C:/Users/lange/.julia/packages/ModelingToolkit/pmsNn/src/systems/systemstructure.jl:635) [8] __structural_simplify(sys::ODESystem, io::Tuple{Vector{Num}, Vector{Num}}; simplify::Bool, kwargs::@Kwargs{}) @ ModelingToolkit [C:\Users\lange.julia\packages\ModelingToolkit\pmsNn\src\systems\systems.jl:74](file:///C:/Users/lange/.julia/packages/ModelingToolkit/pmsNn/src/systems/systems.jl:74) [9] structural_simplify(sys::ODESystem, io::Tuple{Vector{Num}, Vector{Num}}; simplify::Bool, split::Bool, kwargs::@Kwargs{}) @ ModelingToolkit [C:\Users\lange.julia\packages\ModelingToolkit\pmsNn\src\systems\systems.jl:22](file:///C:/Users/lange/.julia/packages/ModelingToolkit/pmsNn/src/systems/systems.jl:22) [10] structural_simplify @ [C:\Users\lange.julia\packages\ModelingToolkit\pmsNn\src\systems\systems.jl:19](file:///C:/Users/lange/.julia/packages/ModelingToolkit/pmsNn/src/systems/systems.jl:19) [inlined] [11] io_preprocessing(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num}; simplify::Bool, kwargs::@Kwargs{}) @ ModelingToolkit [C:\Users\lange.julia\packages\ModelingToolkit\pmsNn\src\systems\abstractsystem.jl:1691](file:///C:/Users/lange/.julia/packages/ModelingToolkit/pmsNn/src/systems/abstractsystem.jl:1691)

…

@ ControlSystemsMTK [C:\Users\lange.julia\packages\ControlSystemsMTK\SqLsj\src\ode_system.jl:205](file:///C:/Users/lange/.julia/packages/ControlSystemsMTK/SqLsj/src/ode_system.jl:205) [17] named_ss(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num}) @ ControlSystemsMTK [C:\Users\lange.julia\packages\ControlSystemsMTK\SqLsj\src\ode_system.jl:164](file:///C:/Users/lange/.julia/packages/ControlSystemsMTK/SqLsj/src/ode_system.jl:164)

[baggepinnen]

I can’t see your model definition so I cannot know what’s wrong. Are you making the same mistake as in this thread?

ModelingToolkit and Linearization

Do you have a fully working solution for your problem now?

[langestefan]

Sorry, this is my first time working with modelingtoolkit.jl and controlsystemsmtk.jl. I am studying the issue you linked.

In the meantime here is all the code relevant to the model here:

Components:

# Capacitors
@named Ci = Capacitor(C=3)
@named Ce = Capacitor(C=1)
@named Ch = Capacitor(C=2)

# Resistances
@named Rih = Resistor(R=0.1)
@named Rie = Resistor(R=0.5)
@named Rea = Resistor(R=3)

# Current sources
@named Φhp_c = Current()
@named Φcv_c = Current()

@named Ta_in = Constant(k=Tambient)
@named Ta = Voltage()

@named gnd = Ground()

Some dummy input signals:

hp_in_val = 100*randn(10)
Φhp(t) = max(t >= 10 ? hp_in_val[end] : hp_in_val[Int(floor(t)) + 1], 0)
@register_symbolic Φhp(t)
@named Φhp_ctrl = TimeVaryingFunction(Φhp)

cv_in_val = 500*randn(10)
Φcv(t) = max(t >= 10 ? cv_in_val[end] : cv_in_val[Int(floor(t)) + 1], 0)
@register_symbolic Φcv(t)
@named Φcv_ctrl = TimeVaryingFunction(Φcv)

Connections:

eqs = [
    # Ci -> Rih, Rie
    connect(Ci.p, Rie.p),
    connect(Ci.p, Rih.p),

    # heat sources
    # Φh, Φcv -> Ch
    connect(Φhp_c.n, Ch.p),
    connect(Φcv_c.n, Ch.p),

    # Ch -> Rih
    connect(Ch.p, Rih.n),

    # Ce -> Rie, Rea
    connect(Ce.p, Rie.n),
    connect(Ce.p, Rea.p),

    # Ta -> Rea
    connect(Ta.p, Rea.n),
    
    # input to source connections
    connect(Ta_in.output, Ta.V),
    connect(Φcv_ctrl.output, Φcv_c.I),
    connect(Φhp_ctrl.output, Φhp_c.I),

    # ground connections
    connect(gnd.g, Ci.n),
    connect(gnd.g, Ce.n),
    connect(gnd.g, Ch.n),
    connect(gnd.g, Φhp_c.p),
    connect(gnd.g, Φcv_c.p),
    connect(gnd.g, Ta.n),
];

components = [
    Ci, Ce, Ch,
    Rih, Rie, Rea,
    Ta_in, Ta,
    Φhp_c, Φcv_c, 
    Φhp_ctrl, Φcv_ctrl, 
    gnd,
];

Problem definition:

@named model = ODESystem(eqs, t, systems = components)
print(equations(model))
sys = structural_simplify(model, allow_parameter=true)
prob = ODEProblem(sys, Pair[], (0, 10.0))

Try to find the state-space model:

using ControlSystemsMTK, ControlSystemsBase, RobustAndOptimalControl

# inputs = Φhp, Φcv
# outputs = Ci.v, Ce.v, Ch.v
RobustAndOptimalControl.named_ss(model, [Φhp_ctrl.output.u, Φcv_ctrl.output.u], [Ci.v, Ce.v, Ch.v])

[baggepinnen]

Yeah it is the same problem, try disconnecting these two connections (comment them out or something).

@langestefan

connect(Φcv_ctrl.output, Φcv_c.I),
    connect(Φhp_ctrl.output, Φhp_c.I)

and then use Φcv_c.I.u (and similar for the other input) as input to the linearization instead.

A nicer solution is to use analysis points, you’ll learn about those in the videos linked in the other thread.

[langestefan]

Thanks! Trying this. I removed the control signals from the list of connections and components, now getting this error:

ExtraVariablesSystemException: The system is unbalanced. There are 60 highest order derivative variables and 59 equations.
More variables than equations, here are the potential extra variable(s):
Ci₊v(t)
Ce₊v(t)
Ch₊v(t)
Φhp_c₊I₊u(t)

[baggepinnen]

Did you remove one connection too much? Only remove those connections associated with the linearization inputs. Keep this one

connect(Ta_in.output, Ta.V),

[langestefan]

Sorry that error came from modelingtoolkit.jl, it was my fault. I meant to copy this error:

RobustAndOptimalControl.named_ss(model, [Φhp_c.I.u, Φcv_c.I.u], [Ci.v, Ce.v, Ch.v])

Initial condition underdefined. Some are missing from the variable map.
Please provide a default (`u0`), initialization equation, or guess
for the following variables:

Any[Φhp_c₊I₊u(t), Φcv_c₊I₊u(t)]

[baggepinnen]

Pass

op = Dict(Φhp_c.I.u=>value1, Φcv_c.I.u=>value2)

as keyword argument to named_ss to indicate which point to linearize around. Replace value 1, value2 with the appropriate operating point.

[langestefan]

Thanks. I must still be doing something stupid, because this:

# inputs = Φhp, Φcv
# outputs = Ci.v, Ce.v, Ch.v
u0 = Dict(Φhp_c.I.u => 1, Φcv_c.I.u => 1)
RobustAndOptimalControl.named_ss(model, [Φhp_c.I.u, Φcv_c.I.u], [Ci.v, Ce.v, Ch.v], op=u0)

Will still give:

Initial condition underdefined. Some are missing from the variable map.
Please provide a default (u0), initialization equation, or guess
for the following variables:
Any[Φhp_c₊I₊u(t), Φcv_c₊I₊u(t)]

[baggepinnen]

Hmm, the internals of MTK related to this has changed recently so there might be a new bug introduced. Mind opening an issue on ModelingToolkit describing this, with code to reproduce? Use linearize instead of named_ss in the issue to not have to load the control packages

[baggepinnen]

The following uses analysis points and linearizes correctly

using ModelingToolkit
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkit: t_nounits as t
# Capacitors
@named Ci = Capacitor(C=3)
@named Ce = Capacitor(C=1)
@named Ch = Capacitor(C=2)

# Resistances
@named Rih = Resistor(R=0.1)
@named Rie = Resistor(R=0.5)
@named Rea = Resistor(R=3)

# Current sources
@named Φhp_c = Current()
@named Φcv_c = Current()

@named Ta_in = Constant(k=20)
@named Ta = Voltage()

@named gnd = Ground()


hp_in_val = 100*randn(10)
Φhp(t) = max(t >= 10 ? hp_in_val[end] : hp_in_val[Int(floor(t)) + 1], 0)
@register_symbolic Φhp(t)
@named Φhp_ctrl = TimeVaryingFunction(Φhp)

cv_in_val = 500*randn(10)
Φcv(t) = max(t >= 10 ? cv_in_val[end] : cv_in_val[Int(floor(t)) + 1], 0)
@register_symbolic Φcv(t)
@named Φcv_ctrl = TimeVaryingFunction(Φcv)


eqs = [
    # Ci -> Rih, Rie
    connect(Ci.p, Rie.p),
    connect(Ci.p, Rih.p),

    # heat sources
    # Φh, Φcv -> Ch
    connect(Φhp_c.n, Ch.p),
    connect(Φcv_c.n, Ch.p),

    # Ch -> Rih
    connect(Ch.p, Rih.n),

    # Ce -> Rie, Rea
    connect(Ce.p, Rie.n),
    connect(Ce.p, Rea.p),

    # Ta -> Rea
    connect(Ta.p, Rea.n),
    
    # input to source connections
    connect(Ta_in.output, Ta.V),
    connect(Φcv_ctrl.output, :u2, Φcv_c.I),
    connect(Φhp_ctrl.output, :u1, Φhp_c.I),

    # ground connections
    connect(gnd.g, Ci.n),
    connect(gnd.g, Ce.n),
    connect(gnd.g, Ch.n),
    connect(gnd.g, Φhp_c.p),
    connect(gnd.g, Φcv_c.p),
    connect(gnd.g, Ta.n),
];

components = [
    Ci, Ce, Ch,
    Rih, Rie, Rea,
    Ta_in, Ta,
    Φhp_c, Φcv_c, 
    Φhp_ctrl, Φcv_ctrl, 
    gnd,
];


@named model = ODESystem(eqs, t, systems = components)
# print(equations(model))
sys = structural_simplify(model, allow_parameter=true)
prob = ODEProblem(sys, Pair[], (0, 10.0))


# using ControlSystemsMTK, ControlSystemsBase, RobustAndOptimalControl

# inputs = Φhp, Φcv
# outputs = Ci.v, Ce.v, Ch.v
ss_matrices, ssys = linearize(model, [:u1, :u2], [Ci.v, Ce.v, Ch.v])

[langestefan]

Sure. Do you mean in ControlSystemMTK.jl or in ModelingToolkit.jl?

[baggepinnen]

This is working for me, see if it works for you as well

using ModelingToolkit
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkit: t_nounits as t
# Capacitors
@named Ci = Capacitor(C=3)
@named Ce = Capacitor(C=1)
@named Ch = Capacitor(C=2)

# Resistances
@named Rih = Resistor(R=0.1)
@named Rie = Resistor(R=0.5)
@named Rea = Resistor(R=3)

# Current sources
@named Φhp_c = Current()
@named Φcv_c = Current()

@named Ta_in = Constant(k=20)
@named Ta = Voltage()

@named gnd = Ground()


hp_in_val = 100*randn(10)
Φhp(t) = max(t >= 10 ? hp_in_val[end] : hp_in_val[Int(floor(t)) + 1], 0)
@register_symbolic Φhp(t)
@named Φhp_ctrl = TimeVaryingFunction(Φhp)

cv_in_val = 500*randn(10)
Φcv(t) = max(t >= 10 ? cv_in_val[end] : cv_in_val[Int(floor(t)) + 1], 0)
@register_symbolic Φcv(t)
@named Φcv_ctrl = TimeVaryingFunction(Φcv)


eqs = [
    # Ci -> Rih, Rie
    connect(Ci.p, Rie.p),
    connect(Ci.p, Rih.p),

    # heat sources
    # Φh, Φcv -> Ch
    connect(Φhp_c.n, Ch.p),
    connect(Φcv_c.n, Ch.p),

    # Ch -> Rih
    connect(Ch.p, Rih.n),

    # Ce -> Rie, Rea
    connect(Ce.p, Rie.n),
    connect(Ce.p, Rea.p),

    # Ta -> Rea
    connect(Ta.p, Rea.n),
    
    # input to source connections
    connect(Ta_in.output, Ta.V),
    connect(Φcv_ctrl.output, :Φcv, Φcv_c.I),
    connect(Φhp_ctrl.output, :Φhp, Φhp_c.I),

    # ground connections
    connect(gnd.g, Ci.n),
    connect(gnd.g, Ce.n),
    connect(gnd.g, Ch.n),
    connect(gnd.g, Φhp_c.p),
    connect(gnd.g, Φcv_c.p),
    connect(gnd.g, Ta.n),
];

components = [
    Ci, Ce, Ch,
    Rih, Rie, Rea,
    Ta_in, Ta,
    Φhp_c, Φcv_c, 
    Φhp_ctrl, Φcv_ctrl, 
    gnd,
];


@named model = ODESystem(eqs, t, systems = components)
# print(equations(model))
sys = structural_simplify(model, allow_parameter=true)
prob = ODEProblem(sys, Pair[], (0, 10.0))


using ControlSystemsMTK, ControlSystemsBase, RobustAndOptimalControl

# inputs = Φhp, Φcv
# outputs = Ci.v, Ce.v, Ch.v
ss_matrices, ssys = linearize(model, [:Φhp, :Φcv], [Ci.v, Ce.v, Ch.v]) # Linearization without control packages

lsys = named_ss(model, [:Φhp, :Φcv], [Ci.v, Ce.v, Ch.v]) # Linearization with control packages
using Plots
bodeplot(lsys, size=(600, 800), margin=4Plots.mm)

[langestefan]

Great, that works! Thanks.

Do you perhaps know if this will also work using ModelingToolkit.linearize_symbolic?

That’s still giving me:

Some specified inputs were not found in system. The following variables were not found Any[:Φhp, :Φcv]

[baggepinnen]

The symbolic linearization does not yet support analysis points, but this works

eqs = [
    # Ci -> Rih, Rie
    connect(Ci.p, Rie.p),
    connect(Ci.p, Rih.p),

    # heat sources
    # Φh, Φcv -> Ch
    connect(Φhp_c.n, Ch.p),
    connect(Φcv_c.n, Ch.p),

    # Ch -> Rih
    connect(Ch.p, Rih.n),

    # Ce -> Rie, Rea
    connect(Ce.p, Rie.n),
    connect(Ce.p, Rea.p),

    # Ta -> Rea
    connect(Ta.p, Rea.n),
    
    # input to source connections
    connect(Ta_in.output, Ta.V),
    # connect(Φcv_ctrl.output, :Φcv, Φcv_c.I),
    # connect(Φhp_ctrl.output, :Φhp, Φhp_c.I),

    # ground connections
    connect(gnd.g, Ci.n),
    connect(gnd.g, Ce.n),
    connect(gnd.g, Ch.n),
    connect(gnd.g, Φhp_c.p),
    connect(gnd.g, Φcv_c.p),
    connect(gnd.g, Ta.n),
];

components = [
    Ci, Ce, Ch,
    Rih, Rie, Rea,
    Ta_in, Ta,
    Φhp_c, Φcv_c, 
    Φhp_ctrl, Φcv_ctrl, 
    gnd,
];

@named model = ODESystem(eqs, t, systems = components)

ss_matrices, ssys = ModelingToolkit.linearize_symbolic(model, [Φcv_c.I.u, Φhp_c.I.u], [Ci.v, Ce.v, Ch.v])

julia> ss_matrices.A
3×3 Matrix{Num}:
 (-1 / Rih₊R + -1 / Rie₊R) / Ci₊C  1 / (Ci₊C*Rie₊R)                        1 / (Ci₊C*Rih₊R)
                1 / (Ce₊C*Rie₊R)        (-1 / Rea₊R + -1 / Rie₊R) / Ce₊C                  0
      1 / (Ch₊C*Rih₊R)

you may want to further simplify the symbolic entries slightly:

julia> simplify.(ss_matrices.A)
3×3 Matrix{Num}:
             (-Rie₊R - Rih₊R) / (Ci₊C*Rie₊R*Rih₊R)  …   1 / (Ci₊C*Rih₊R)
           1 / (Ce₊C*Rie₊R)                                            0
 1 / (Ch₊C*Rih₊R) 

[langestefan]

That’s awesome, thanks a lot. If by analysis points you mean the point to linearize at, I guess I don’t need those since the model is linear in the v's and \Phi's anyway?

[baggepinnen]

@langestefan

If by analysis points you mean the point to linearize at

No, I mean this

@langestefan

I guess I don’t need those since the model is linear in the vvv’s and \PhiΦ\Phi’s anyway?

With MTK, it turns out that you need to provide the operating point anyway since the linearization is numeric (using forward-mode automatic differentiation), and the fact that the system is linear isn’t known to the AD tool when it needs the numerical input point at which to linearize.

[langestefan]

I see, that makes sense. Seems I have a lot of reading up to do.

Anyway, thanks again, and I hope others who stumble upon this can also use it.

ps. Do you think it would be nice to have a specific option for systems for which you already know they are linear in the variables (and fail if they are not)? Which would essentially only rearrange the variables into a matrix-vector product