Composing SDESystems does not work

[hstrey]

I am trying to compose SDESystems, and when I do, Julia goes into an infinite loop after I execute compose. Am I doing it the right way? (see below for a slightly better way of attempting to compose SDEs).

using ModelingToolkit

@variables t
D = Differential(t)

function van_der_pol(;name, θ=1.0,ϕ=0.1)
    params  = @parameters θ=θ ϕ=ϕ
    sts = @variables x(t)=1.0 y(t)=1.0 jcn(t)=0.0

    eqs = [D(x) ~ y + jcn,
           D(y) ~ θ*(1-x^2)*y - x]

    noiseeqs = [ϕ,ϕ]

    return SDESystem(eqs,noiseeqs,t,sts,params; name=name)
end

@named VP1 = van_der_pol()
@named VP2 = van_der_pol()

eqs = [VP1.jcn ~ VP2.x,
        VP2.jcn ~ VP1.x]

sys = [VP1,VP2]

@named coupledVP = compose(eqs,sys)

[isaacsas]

Don't eqs need to be in a system to compose?

[hstrey]

For ODESystems, that is the way we use compose and it works for ODESystems. @ChrisRackauckas told me that this should work also for SDEs, but maybe I am missing a detail or two.

[hstrey]

@isaacsas you are correct. For ODESystems the correct compose is: @named coupledVP = compose(ODESystem(eqs;name=:connected),sys) but there is no equivalent way for SDESystem.

[isaacsas]

Yeah, SDESystems are missing both compose and flattening when generating code for DifferentialEquations. Without the latter you might be able to compose systems but you wouldn't be able to actually generate a correct SDEProblem from the system. So there is some work needed to get this together I think.

[hstrey]

I just tried: @named coupledVP = compose(SDESystem(eqs,[],t,[],[];name=:connected),sys) that creates the correct set of equations, but structural_simplify fails

[isaacsas]

Does it actually get the correct noise equations too though?

[hstrey]

here is the updated code:

using ModelingToolkit

@variables t
D = Differential(t)

function van_der_pol(;name, θ=1.0,ϕ=0.1)
    params  = @parameters θ=θ ϕ=ϕ
    sts = @variables x(t)=1.0 y(t)=1.0 jcn(t)=0.0

    eqs = [D(x) ~ y + jcn,
           D(y) ~ θ*(1-x^2)*y - x]

    noiseeqs = [ϕ,ϕ]

    return SDESystem(eqs,noiseeqs,t,sts,params; name=name)
end

@named VP1 = van_der_pol()
@named VP2 = van_der_pol()

eqs = [VP1.jcn ~ VP2.x,
        VP2.jcn ~ VP1.x]

sys = [VP1,VP2]

@named coupledVP = compose(SDESystem(eqs,[],t,[],[];name=:connected),sys)
coupledVPs = structural_simplify(coupledVP)

structural_simplify fails with:

ERROR: BoundsError: attempt to access 10-element Vector{Vector{Int64}} at index [11]
Stacktrace:
  [1] getindex
    @ ./array.jl:924 [inlined]
  [2] 𝑠neighbors (repeats 2 times)
    @ ~/.julia/packages/ModelingToolkit/QxF9d/src/bipartite_graph.jl:345 [inlined]
  [3] find_eq_solvables!(state::TearingState{SDESystem}, ieq::Int64, to_rm::Vector{Int64}, coeffs::Vector{Int64}; may_be_zero::Bool, allow_symbolic::Bool, allow_parameter::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/QxF9d/src/structural_transformation/utils.jl:185
  [4] find_eq_solvables!(state::TearingState{SDESystem}, ieq::Int64, to_rm::Vector{Int64}, coeffs::Vector{Int64})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/QxF9d/src/structural_transformation/utils.jl:175
  [5] linear_subsys_adjmat!(state::TearingState{SDESystem}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/QxF9d/src/structural_transformation/utils.jl:257
  [6] linear_subsys_adjmat!
    @ ~/.julia/packages/ModelingToolkit/QxF9d/src/structural_transformation/utils.jl:245 [inlined]
  [7] alias_eliminate_graph!(state::TearingState{SDESystem}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/alias_elimination.jl:7
  [8] alias_eliminate_graph!
    @ ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/alias_elimination.jl:6 [inlined]
  [9] alias_elimination!(state::TearingState{SDESystem}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/alias_elimination.jl:52
 [10] alias_elimination!(state::TearingState{SDESystem})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/alias_elimination.jl:48
 [11] _structural_simplify!(state::TearingState{SDESystem}, io::Nothing; simplify::Bool, check_consistency::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit.SystemStructures ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/systemstructure.jl:538
 [12] #structural_simplify!#23
    @ ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/systemstructure.jl:523 [inlined]
 [13] structural_simplify(sys::SDESystem, io::Nothing; simplify::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/systems.jl:39
 [14] structural_simplify (repeats 2 times)
    @ ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/systems.jl:19 [inlined]
 [15] top-level scope
    @ ~/Documents/programming/Modeling-SDEs/coupled_SDEs.jl:27

[hstrey]

Does it actually get the correct noise equations too though?

no, the noise equations are missing.

julia> coupledVP.noiseeqs
Any[][1:0]

[isaacsas]

Right, so it needs a custom compose too in addition to flattening code in generating the needed DiffEqs functions.

[hstrey]

I just noticed this pull request: https://github.com/SciML/ModelingToolkit.jl/pull/1642 @YingboMa, is this still current?

[YingboMa]

This works:

using ModelingToolkit

@variables t
D = Differential(t)

function van_der_pol(;name, θ=1.0,ϕ=0.1)
    params  = @parameters θ=θ
    sts = @variables x(t)=1.0 y(t)=1.0 jcn(t)=0.0
    @brownian ϕ

    eqs = [D(x) ~ y + jcn + ϕ,
           D(y) ~ θ*(1-x^2)*y - x + ϕ]

    return System(eqs, t; name=name)
end

@named VP1 = van_der_pol()
@named VP2 = van_der_pol()

eqs = [VP1.jcn ~ VP2.x,
        VP2.jcn ~ VP1.x]

sys = [VP1,VP2]

@named coupledVP = compose(System(eqs,t;name=:connected),sys)
coupledVPs = structural_simplify(coupledVP)

[hstrey]

@YingboMa thanks so much

[arnold-pdev]

This works:

using ModelingToolkit

@variables t
D = Differential(t)

function van_der_pol(;name, θ=1.0,ϕ=0.1)
    params  = @parameters θ=θ
    sts = @variables x(t)=1.0 y(t)=1.0 jcn(t)=0.0
    @brownian ϕ

    eqs = [D(x) ~ y + jcn + ϕ,
           D(y) ~ θ*(1-x^2)*y - x + ϕ]

    return System(eqs, t; name=name)
end

@named VP1 = van_der_pol()
@named VP2 = van_der_pol()

eqs = [VP1.jcn ~ VP2.x,
        VP2.jcn ~ VP1.x]

sys = [VP1,VP2]

@named coupledVP = compose(System(eqs,t;name=:connected),sys)
coupledVPs = structural_simplify(coupledVP)

Hi Yingbo-

Could your provide an example of a construction of a (ODE? SDE?)Problem using this System? I am running into the issue that the number of equations and variables don't match.