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

Specific Domains Modelling & Simulations
1

Hi, I essentially want to solve an easier version of the copy-paste example provided in the docs here. 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:


@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!

2

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

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

1 Like
3

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

4

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!

5

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

6

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:

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!

7

In that case it would be a state equation.

8

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

9

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

10

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:

@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

11

It should be

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

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

1 Like
12

You also need to flatten and reduce the system:

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.

1 Like
13

Thanks, Yingbo and Shashi!

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

14

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.

1 Like