# Putting linear ODESystem in matrix-vector form

**URL:** <https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333>\
**Category:** Modelling & Simulations\
**Tags:** modelingtoolkit, symbolics\
**Created:** [April 22, 2024, 12:49pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333 "2024-04-22T12:49:37Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [April 22, 2024, 12:49pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/1 "2024-04-22T12:49:37Z")

</div>

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](https://github.com/JuliaControl/ControlSystemsMTK.jl) but I was unsuccesful in getting it to work.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 22, 2024, 12:58pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/2 "2024-04-22T12:58:16Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [April 22, 2024, 1:14pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/3 "2024-04-22T13:14:55Z")

</div>

I tried the following [Symbolic Linearization](https://juliacontrol.github.io/ControlSystemsMTK.jl/dev/#Symbolic-linearization) and [From ModelingToolkit to ControlSystems](https://juliacontrol.github.io/ControlSystemsMTK.jl/dev/#From-ModelingToolkit-to-ControlSystems)

```julia
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)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 22, 2024, 1:16pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/4 "2024-04-22T13:16:37Z")

</div>

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](https://discourse.julialang.org/t/modelingtoolkit-and-linearization/113241/15):
>
> Do you have a fully working solution for your problem now?

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [April 22, 2024, 1:26pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/5 "2024-04-22T13:26:11Z")

</div>

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:

```julia
# 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:

```julia
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:

```julia
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:

```julia
@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:

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

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 22, 2024, 1:28pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/6 "2024-04-22T13:28:29Z")

</div>

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

> [@langestefan](#):
>
> ```julia
> 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.

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [April 22, 2024, 1:41pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/7 "2024-04-22T13:41:30Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 22, 2024, 1:42pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/8 "2024-04-22T13:42:37Z")

</div>

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

```julia
connect(Ta_in.output, Ta.V),

```

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [April 22, 2024, 1:50pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/9 "2024-04-22T13:50:08Z")

</div>

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

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

```

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

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 22, 2024, 2:02pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/10 "2024-04-22T14:02:07Z")

</div>

Pass

```julia
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.

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [April 22, 2024, 2:11pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/11 "2024-04-22T14:11:13Z")

</div>

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

```julia
# 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)]

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 22, 2024, 2:16pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/12 "2024-04-22T14:16:27Z")

</div>

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 😊

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 22, 2024, 2:22pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/13 "2024-04-22T14:22:22Z")

</div>

The following uses analysis points and linearizes correctly

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

```

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [April 22, 2024, 2:24pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/14 "2024-04-22T14:24:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 22, 2024, 2:31pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/15 "2024-04-22T14:31:19Z")

</div>

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

```julia
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)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/9/09596cc623a2deafcd1ae2d0068b5731c7cfdfaa.png)

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [April 22, 2024, 2:38pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/16 "2024-04-22T14:38:02Z")

</div>

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]

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 22, 2024, 2:46pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/17 "2024-04-22T14:46:36Z")

</div>

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

```julia
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
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
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) 

```

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [April 22, 2024, 3:59pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/18 "2024-04-22T15:59:37Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 22, 2024, 4:04pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/19 "2024-04-22T16:04:11Z")

</div>

> [@langestefan](#):
>
> If by analysis points you mean the point to linearize at

No, I mean [this](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/)

> [@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.

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [April 22, 2024, 4:23pm UTC](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333/20 "2024-04-22T16:23:19Z")

</div>

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

[Next page](https://discourse.julialang.org/t/putting-linear-odesystem-in-matrix-vector-form/113333.md?page=2)