Connecting components by adding term to DE in ModelingToolkit

Specific Domains Modelling & Simulations
1

I am trying to connect components by adding function of variables from each component to one of the DEs in each component. In other words, for components 1 and 2 (and their respective variables z_1 and z_2)

\frac{d z_1}{dt} = g(x_1,y_1,z_1, \beta) +f (z_1, z_2, p)

\frac{d z_2}{dt} = g(x_2,y_2, z_2, \beta) + f(z_2, z_1, p)

where \frac{d z_1}{dt} = g(x_1,y_1,z_1, \beta) is the original DE in each component and f(z_1, z_2, p) is the connection with the other component that I would like to add.

I could presumably do this using a continuous callback in the ODEProblem, but was wondering if there is a more elegant way to do this with ModelingToolkit.

I have attempted to do this using compose, but am having issues getting the equations to combine.

using ModelingToolkit, OrdinaryDiffEq

function Lorenz(; name)
    @parameters σ ρ β
    @variables t x(t) y(t) z(t)
    D = Differential(t)

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

    @named sys = ODESystem(eqs; name=name)
end

@named Lorenz1 = Lorenz()
@named Lorenz2 = Lorenz()

f(b1, b2, a) = a * (b1- b2)

@parameters c
connections = [D(Lorenz1.z) ~ + f(Lorenz1.z , Lorenz2.z, c),
                D(Lorenz2.z) ~ + f(Lorenz2.z , Lorenz1.z, c)]
connected = compose(ODESystem(connections, t,  name = :connected), Lorenz1, Lorenz2)
connected_simp = structural_simplify(connected)

u0 = [D(Lorenz1.x) => 2.0,
    Lorenz1.x => 1.0,
    Lorenz1.y => 0.0,
    Lorenz1.z => 0.0,
    D(Lorenz2.x) => 2.0,
    Lorenz2.x => 1.0,
    Lorenz2.y => 0.0,
    Lorenz2.z => 0.0]

p = [Lorenz1.σ => 28.0,
    Lorenz1.ρ => 10.0,
    Lorenz1.β => 8 / 3,
    Lorenz2.σ => 28.0,
    Lorenz2.ρ => 10.0,
    Lorenz2.β => 8 / 3,
    a => 0.005]


prob = ODEProblem(connected, u0, tspan, p)
sol = solve(prob, Tsit5())
using Plots;
plot(sol, idxs = (Lorenz1.x, Lorenz1.y))

julia> connected_simp = structural_simplify(connected)
ERROR: ExtraEquationsSystemException: The system is unbalanced. There are 6 highest order derivative variables and 8 equations.
More equations than variables, here are the potential extra equation(s):
 Differential(t)(Lorenz1₊z(t)) ~ Lorenz1₊x(t)*Lorenz1₊y(t) - Lorenz1₊β*Lorenz1₊z(t)
 Differential(t)(Lorenz2₊z(t)) ~ Lorenz2₊x(t)*Lorenz2₊y(t) - Lorenz2₊β*Lorenz2₊z(t)
Stacktrace:
 [1] error_reporting(state::TearingState{ODESystem}, bad_idxs::Vector{Int64}, n_highest_vars::Int64, iseqs::Bool, orig_inputs::Set{Any})
   @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/ZmLra/src/structural_transformation/utils.jl:44
 [2] check_consistency(state::TearingState{ODESystem}, orig_inputs::Set{Any})
   @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/ZmLra/src/structural_transformation/utils.jl:85
 [3] _structural_simplify!(state::TearingState{ODESystem}, io::Nothing; simplify::Bool, check_consistency::Bool, fully_determined::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit.SystemStructures ~/.julia/packages/ModelingToolkit/ZmLra/src/systems/systemstructure.jl:612
 [4] _structural_simplify!
   @ ~/.julia/packages/ModelingToolkit/ZmLra/src/systems/systemstructure.jl:600 [inlined]
 [5] structural_simplify!(state::TearingState{ODESystem}, io::Nothing; simplify::Bool, check_consistency::Bool, fully_determined::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit.SystemStructures ~/.julia/packages/ModelingToolkit/ZmLra/src/systems/systemstructure.jl:566
 [6] structural_simplify(sys::ODESystem, io::Nothing; simplify::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/ZmLra/src/systems/systems.jl:39
 [7] structural_simplify
   @ ~/.julia/packages/ModelingToolkit/ZmLra/src/systems/systems.jl:19 [inlined]
 [8] structural_simplify(sys::ODESystem)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/ZmLra/src/systems/systems.jl:19
 [9] top-level scope
   @ ~/connection_MWE:50
2

The short answer is you can’t modify equations by composing or extending, only add new equations. There used to be a nice explanation in another almost identical question, but by search skills are failing to find it.

The longer answer is that you have to prepare a place to connect systems by adding a variable whose value can be provided externally. The @connector macro was created to help with this and there is a nice tutorial here.

A simple way to do it in your case with minimal changes would be something like (untested).

function Lorenz(; name)
    ps = @parameters σ ρ β
    states = @variables t x(t) y(t) z(t) c(t)
    D = Differential(t)

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

    ODESystem(eqs, t, states, ps; name)
end

and then

connections = [Lorenz1.c ~ f(Lorenz1.z , Lorenz2.z, c),
               Lorenz2.c ~ f(Lorenz2.z , Lorenz1.z, c)]
3

Thank you for your response! The simple way you suggest works:

using ModelingToolkit, OrdinaryDiffEq

function Lorenz(; name)
    @parameters σ ρ β
    @variables t x(t) y(t) z(t) c(t)
    D = Differential(t)

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

    @named sys = ODESystem(eqs; name=name)
end

@named Lorenz1 = Lorenz() 
@named Lorenz2 = Lorenz()

f(b1, b2, a) = a * (b1- b2)
@variables t
@parameters C
D = Differential(t)
connections = [Lorenz1.c ~ f(Lorenz1.z , Lorenz2.z, C),
               Lorenz2.c ~ f(Lorenz2.z , Lorenz1.z, C)]
connected = compose(ODESystem(connections, t, name = :connected), Lorenz1, Lorenz2)
connected = structural_simplify(connected)

u0 = [D(Lorenz1.x) => 2.0,
    Lorenz1.x => 1.0,
    Lorenz1.y => 0.0,
    Lorenz1.z => 0.0,
    Lorenz1.c => 0.0, 
    D(Lorenz2.x) => 5.0,
    Lorenz2.x => 0.0,
    Lorenz2.y => 1.0,
    Lorenz2.z => 0.5,
    Lorenz2.c  => 0.0,]

p = [Lorenz1.σ => 28.0,
    Lorenz1.ρ => 10.0,
    Lorenz1.β => 8 / 3,
    Lorenz2.σ => 28.0,
    Lorenz2.ρ => 10.0,
    Lorenz2.β => 8 / 3,
    C => 0.005]

prob = ODEProblem(connected, u0, [0., 1000.], p)
sol = solve(prob, Tsit5())
using Plots;
plot(sol, idxs = (Lorenz1.x, Lorenz1.y))

The tutorial you linked is very helpful. Do I understand it correctly that using the @connector macro would do something similar to this?