# Putting linear ODESystem in matrix-vector form

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.

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

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)

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

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


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

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.

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)

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


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


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.

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

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

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


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

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)


Great, that works! Thanks.

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

That’s still giving me:

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)


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?

No, I mean this

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.

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