# Component-based modelling in ModelingToolkit.jl with nonsingular mass matrix

**URL:** <https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140>\
**Category:** Modelling & Simulations\
**Created:** [October 27, 2020, 6:50pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140 "2020-10-27T18:50:05Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![Dhruva2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dhruva2/32/18475_2.png) [@Dhruva2](https://discourse.julialang.org/u/Dhruva2)\
**Post date:** [October 27, 2020, 6:50pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/1 "2020-10-27T18:50:05Z")

</div>

Hi, I essentially want to solve an **easier** version of the copy-paste example provided in the docs [here](https://mtk.sciml.ai/stable/tutorials/ode_modeling/). Instead of connecting the two lorenz systems via an algebraic equality, I want to add the state of one system as a time-varying input to another system. I tried to code this below:

```using

@parameters t σ ρ β
@parameters inp(t) #extra input to the first lorenz system I've added
@variables x(t) y(t) z(t) 
@derivatives D'~t

eqs1 = [D(x) ~ σ*(y-x) + inp,
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

eqs2 = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

lorenz1 = ODESystem(eqs1,name=:lorenz1)
lorenz2 = ODESystem(eqs2,name=:lorenz2)

connections = [lorenz1.inp ~ lorenz2.x]

connected = ODESystem(connections,t,[],[],systems=[lorenz1,lorenz2])

u0 = [lorenz1.x => 1.0,
      lorenz1.y => 0.0,
      lorenz1.z => 0.0,
      lorenz2.x => 0.0,
      lorenz2.y => 1.0,
      lorenz2.z => 0.0]

p = [lorenz1.σ => 10.0,
      lorenz1.ρ => 28.0,
      lorenz1.β => 8/3,
      lorenz2.σ => 10.0,
      lorenz2.ρ => 28.0,
      lorenz2.β => 8/3,
      lorenz1.inp => 0.]

tspan = (0.0,100.0)
prob = ODEProblem(connected,u0,tspan,p)
sol = solve(prob,Rodas5())

```

However, I can’t get it to work. For the above code, I get the error.

> ERROR: Only semi-explicit constant mass matrices are currently supported. Faulty equation: Equation(lorenz1₊inp, lorenz2₊x(t)).

Any idea how to get this to work? Thank you very much in advance!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 27, 2020, 8:47pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/2 "2020-10-27T20:47:50Z")

</div>

You need to define `inp` as a `pins`. Then it’ll be able to perform the connection. See this test as an example:

[https://github.com/SciML/ModelingToolkit.jl/blob/master/test/inputoutput.jl](https://github.com/SciML/ModelingToolkit.jl/blob/master/test/inputoutput.jl)

Note that we haven’t added documentation for pins and observables yet, so it’s not quite released, but it should work.

---

<div class="post-metadata">

**Author:** ![Dhruva2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dhruva2/32/18475_2.png) [@Dhruva2](https://discourse.julialang.org/u/Dhruva2)\
**Post date:** [October 27, 2020, 8:54pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/3 "2020-10-27T20:54:58Z")

</div>

awesome, just what I was looking for! Thanks a bunch

---

<div class="post-metadata">

**Author:** ![Dhruva2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dhruva2/32/18475_2.png) [@Dhruva2](https://discourse.julialang.org/u/Dhruva2)\
**Post date:** [October 27, 2020, 9:51pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/4 "2020-10-27T21:51:10Z")

</div>

However, the test code you gave me fails, and it doesn’t seem like the connection is taking effect in the system. Specifically, the last line of code:  
`@test isequal(simplifyeqs(equations(connected)), simplifyeqs(collapsed_eqs))`  
fails, as the `equations(connected)` have not substituted in the lorenz2 variables in the lorenz1 equations (and vice versa).

Is this just because the functionality is not coded up yet, or does this highlight a bug? If the former, do you have a guesstimate of when the functionality will be usable?

Thanks again!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 27, 2020, 9:59pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/5 "2020-10-27T21:59:25Z")

</div>

It’s a bug. @shashi we forgot to add this test to [https://github.com/SciML/ModelingToolkit.jl/blob/master/test/runtests.jl](https://github.com/SciML/ModelingToolkit.jl/blob/master/test/runtests.jl) 🤦‍♂️. We’ll get this fixed up. Sorry about that. Good thing the feature isn’t documented yet haha.

---

<div class="post-metadata">

**Author:** ![Dhruva2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dhruva2/32/18475_2.png) [@Dhruva2](https://discourse.julialang.org/u/Dhruva2)\
**Post date:** [October 27, 2020, 10:20pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/6 "2020-10-27T22:20:56Z")

</div>

No problem! You shouldn’t apologise for the bug, I should thank you for making such cool functionality!

One final point of confusion that might be better answered after the bugfix:

```julia
connected = ODESystem(Equation[],t,[],[],observed=connections,systems=[lorenz1,lorenz2])

```

Suppose my connections do not involve any observed variables (i.e. the lorenz.u in thie example), but directly connect to a state variables in the system, e.g.  
` connections = [lorenz1.F ~ lorenz2.x]`  
In this case do I still need the `observed = connections` argument, or should I put the connections as the first argument of ODESystem(). I tried doing the latter but I got the same mass matrix error I highlighted in the original post.

Cheers!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 27, 2020, 10:31pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/7 "2020-10-27T22:31:25Z")

</div>

> [@Dhruva2](#):
>
> In this case do I still need the `observed = connections` argument, or should I put the connections as the first argument of ODESystem(). I tried doing the latter but I got the same mass matrix error I highlighted in the original post.

In that case it would be a state equation.

---

<div class="post-metadata">

**Author:** ![Dhruva2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dhruva2/32/18475_2.png) [@Dhruva2](https://discourse.julialang.org/u/Dhruva2)\
**Post date:** [October 27, 2020, 11:05pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/8 "2020-10-27T23:05:03Z")

</div>

You mean it’s a standard ODE that can be written without components? If so, yes that’s true!

However, I was hoping to use components anyway, so I can separate the different subsets of equations. In my case I want to model a network of three, conductance-based neurons, each of which has quite complicated dynamics. Each neuron has similarly named state variables (V(t), Ca(t), etc…). Having the components and connections allows me to easily change the pattern of synaptic connections, and to distinguish between the different neurons in an easy manner (since I get neuron1.V(t), neuron2.V(t), etc.).

Cheers

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 28, 2020, 3:47am UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/9 "2020-10-28T03:47:22Z")

</div>

> [@Dhruva2](#):
>
> You mean it’s a standard ODE that can be written without components? If so, yes that’s true!

I mean it’s a state equation of the connector system, instead of an observed variable.

---

<div class="post-metadata">

**Author:** ![Dhruva2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dhruva2/32/18475_2.png) [@Dhruva2](https://discourse.julialang.org/u/Dhruva2)\
**Post date:** [October 29, 2020, 11:44am UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/10 "2020-10-29T11:44:26Z")

</div>

Ah I see, thanks. I tried to implement a modified form of the link you sent me as you just suggested, but it still returns an error. Code:

```nohighlight
@parameters t σ ρ β
@variables x(t) y(t) z(t) F(t)
@derivatives D'~t
eqs = [D(x) ~ σ*(y-x) + F,
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
lorenz1 = ODESystem(eqs,pins=[F],name=:lorenz1)
lorenz2 = ODESystem(eqs,pins=[F],name=:lorenz2)

connections = [lorenz1.F ~ lorenz2.x,
               lorenz2.F ~ lorenz1.x]
connected = ODESystem(connections,t,[],[],systems=[lorenz1,lorenz2])
sys = connected

@variables lorenz1₊F lorenz2₊F
u0 = [lorenz1.x => 1.0,
      lorenz1.y => 0.0,
      lorenz1.z => 0.0,
      lorenz2.x => 0.0,
      lorenz2.y => 1.0,
      lorenz2.z => 0.0]

p = [lorenz1.σ => 10.0,
      lorenz1.ρ => 28.0,
      lorenz1.β => 8/3,
      lorenz2.σ => 10.0,
      lorenz2.ρ => 28.0,
      lorenz2.β => 8/3]
prob = ODEProblem(connected, u0, (0.,10.), p)
sol = solve(prob, Tsit5()) 

```

Still returns the same error:  
ERROR: Only semi-explicit constant mass matrices are currently supported. Faulty equation: Equation(lorenz1₊F(t), lorenz2₊x(t)).

Am I doing something wrong?

Thanks

---

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [November 6, 2020, 4:51pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/11 "2020-11-06T16:51:10Z")

</div>

It should be

```julia
prob = ODEProblem(alias_elimination(connected), u0, (0.,10.), p)

```

since we need to do alias elimination to reduce the connecting equations away.

---

<div class="post-metadata">

**Author:** ![shashi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashi/32/1824_2.png) [@shashi](https://discourse.julialang.org/u/shashi)\
**Post date:** [November 6, 2020, 4:52pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/12 "2020-11-06T16:52:29Z")

</div>

> [@ChrisRackauckas](#):
>
> `sys = connected`

You also need to flatten and reduce the system:

```julia
julia> sys = ModelingToolkit.alias_elimination(ModelingToolkit.flatten(connected));

julia> ModelingToolkit.equations(sys)
6-element Array{Equation,1}:
 Equation(Differential(lorenz1₊x(t)), (lorenz1₊σ * (lorenz1₊y(t) - lorenz1₊x(t))) + (lorenz2₊x(t) + lorenz2₊y(t) + (-1 * lorenz2₊z(t))))
 Equation(Differential(lorenz1₊y(t)), (lorenz1₊x(t) * (lorenz1₊ρ - lorenz1₊z(t))) - lorenz1₊y(t))
 Equation(Differential(lorenz1₊z(t)), (lorenz1₊x(t) * lorenz1₊y(t)) - (lorenz1₊β * lorenz1₊z(t)))
 Equation(Differential(lorenz2₊x(t)), (lorenz2₊σ * (lorenz2₊y(t) - lorenz2₊x(t))) + (lorenz1₊x(t) + lorenz1₊y(t) + (-1 * lorenz1₊z(t))))
 Equation(Differential(lorenz2₊y(t)), (lorenz2₊x(t) * (lorenz2₊ρ - lorenz2₊z(t))) - lorenz2₊y(t))
 Equation(Differential(lorenz2₊z(t)), (lorenz2₊x(t) * lorenz2₊y(t)) - (lorenz2₊β * lorenz2₊z(t)))

```

So Chris’s code with this line should work.

---

<div class="post-metadata">

**Author:** ![Dhruva2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dhruva2/32/18475_2.png) [@Dhruva2](https://discourse.julialang.org/u/Dhruva2)\
**Post date:** [November 6, 2020, 7:24pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/13 "2020-11-06T19:24:37Z")

</div>

Thanks, Yingbo and Shashi!

I wasn’t aware of the flattening and alias elimination functionality. It works with them!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 6, 2020, 8:51pm UTC](https://discourse.julialang.org/t/component-based-modelling-in-modelingtoolkit-jl-with-nonsingular-mass-matrix/49140/14 "2020-11-06T20:51:04Z")

</div>

We need to document it, so not your fault :). We’re still in the process of getting the full functionality of the structural transformations together, but we’re fairly close now.