Skip to content

Commit 40073b6

Browse files
Fix hierarchical modeling with DDEs
Fixes #2232 The issue before was that the namespace handling was not recursing into variables, meaning that when namespace was applied, states and parameters used in the delay definition were not namespaced. This was the output before: ```julia using ModelingToolkit, DifferentialEquations @variables t D = Differential(t) @parameters x(..) a function oscillator(;name, k=1.0, τ=0.01) @parameters k=k τ=τ @variables x(..)=0.1 y(t)=0.1 jcn(t)=0.0 delx(t) eqs = [D(x(t)) ~ y, D(y) ~ -k*x(t-τ)+jcn, delx ~ x(t-τ)] return System(eqs; name=name) end @nAmed osc1 = oscillator(k=1.0, τ=0.01) @nAmed osc2 = oscillator(k=2.0, τ=0.04) eqs = [osc1.jcn ~ osc2.delx, osc2.jcn ~ osc1.delx] @nAmed coupledOsc = compose(System(eqs, t; name=:connected), [osc1, osc2]) julia> equations(coupledOsc) 8-element Vector{Equation}: osc1₊jcn(t) ~ osc2₊delx(t) osc2₊jcn(t) ~ osc1₊delx(t) Differential(t)(osc1₊x(t)) ~ osc1₊y(t) Differential(t)(osc1₊y(t)) ~ osc1₊jcn(t) - osc1₊k*osc1₊x(t - τ) osc1₊delx(t) ~ osc1₊x(t - τ) Differential(t)(osc2₊x(t)) ~ osc2₊y(t) Differential(t)(osc2₊y(t)) ~ osc2₊jcn(t) - osc2₊k*osc2₊x(t - τ) osc2₊delx(t) ~ osc2₊x(t - τ) ``` You can see the issue is that it's using `τ` instead of a namespaced `tau`. With this PR: ```julia julia> equations(coupledOsc) 8-element Vector{Equation}: osc1₊jcn(t) ~ osc2₊delx(t) osc2₊jcn(t) ~ osc1₊delx(t) Differential(t)(osc1₊x(t)) ~ osc1₊y(t) Differential(t)(osc1₊y(t)) ~ osc1₊jcn(t) - osc1₊k*osc1₊x(t - osc1₊τ) osc1₊delx(t) ~ osc1₊x(t - osc1₊τ) Differential(t)(osc2₊x(t)) ~ osc2₊y(t) Differential(t)(osc2₊y(t)) ~ osc2₊jcn(t) - osc2₊k*osc2₊x(t - osc2₊τ) osc2₊delx(t) ~ osc2₊x(t - osc2₊τ) ``` However, something is going on with structural_simplify here: ```julia simpsys = structural_simplify(coupledOsc) equations(simpsys) 6-element Vector{Equation}: Differential(t)(osc1₊x(t)) ~ osc1₊y(t) Differential(t)(osc1₊y(t)) ~ osc1₊jcn(t) - osc1₊k*osc1₊x(t - osc1₊τ) Differential(t)(osc2₊x(t)) ~ osc2₊y(t) Differential(t)(osc2₊y(t)) ~ osc2₊jcn(t) - osc2₊k*osc2₊x(t - osc2₊τ) 0 ~ osc1₊x(t - osc1₊τ) - osc1₊delx(t) 0 ~ osc2₊x(t - osc2₊τ) - osc2₊delx(t) ``` This errors: ```julia julia> prob = DDEProblem(structural_simplify(coupledOsc), [0.1,0.1,0.1,0.1], (0, 50), constant_lags=[osc1.τ, osc2.τ]) ERROR: ArgumentError: Equations (6), states (4), and initial conditions (4) are of different lengths. To allow a different number of equations than states use kwarg check_length=false. Stacktrace: [1] check_eqs_u0(eqs::Vector{Equation}, dvs::Vector{Any}, u0::Vector{Float64}; check_length::Bool, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:t, :has_difference, :constant_lags), Tuple{Int64, Bool, Vector{Num}}}}) @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\ua9Sp\src\systems\abstractsystem.jl:1714 [2] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Float64}, parammap::SciMLBase.NullParameters; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, kwargs::Base.Pairs{Symbol, Any, NTuple{4, Symbol}, NamedTuple{(:t, :has_difference, :check_length, :constant_lags), Tuple{Int64, Bool, Bool, Vector{Num}}}}) @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\abstractodesystem.jl:736 [3] (DDEProblem{true})(sys::ODESystem, u0map::Vector{Float64}, tspan::Tuple{Int64, Int64}, parammap::SciMLBase.NullParameters; callback::Nothing, check_length::Bool, kwargs::Base.Pairs{Symbol, Vector{Num}, Tuple{Symbol}, NamedTuple{(:constant_lags,), Tuple{Vector{Num}}}}) @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\abstractodesystem.jl:930 [4] DDEProblem @ C:\Users\accou\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\abstractodesystem.jl:923 [inlined] [5] DDEProblem(::ODESystem, ::Vector{Float64}, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Vector{Num}, Tuple{Symbol}, NamedTuple{(:constant_lags,), Tuple{Vector{Num}}}}) @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\abstractodesystem.jl:921 [6] top-level scope @ REPL[2]:1 ``` The reason is because the last two equations are not deleted. ```julia julia> states(simpsys) 4-element Vector{Any}: osc1₊x(t) osc1₊y(t) osc2₊x(t) osc2₊y(t) ``` ```julia julia> states(simpsys) 4-element Vector{Any}: osc1₊x(t) osc1₊y(t) osc2₊x(t) osc2₊y(t) julia> observed(simpsys) 4-element Vector{Equation}: osc2₊delx(t) ~ -0.0 osc1₊delx(t) ~ -0.0 osc1₊jcn(t) ~ osc2₊delx(t) osc2₊jcn(t) ~ osc1₊delx(t) ``` These last two equations ``` 0 ~ osc1₊x(t - osc1₊τ) - osc1₊delx(t) 0 ~ osc2₊x(t - osc2₊τ) - osc2₊delx(t) ``` should be observed equations in terms of the constant value (via delay): ``` osc1₊delx(t) ~ osc1₊x(t - osc1₊τ) osc2₊delx(t) ~ osc2₊x(t - osc2₊τ) ``` and removed from the equations, but instead they are still in there so it things there's 6 equations and errors. Is some value matching thing missing the deletion of these equations?
1 parent 74fb734 commit 40073b6

File tree

1 file changed

+9
-9
lines changed

1 file changed

+9
-9
lines changed

src/systems/abstractsystem.jl

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -425,7 +425,7 @@ function GlobalScope(sym::Union{Num, Symbolic})
425425
end
426426
end
427427

428-
renamespace(sys, eq::Equation) = namespace_equation(eq, sys)
428+
renamespace(sys, eq::Equation) = @show namespace_equation(eq, sys)
429429

430430
renamespace(names::AbstractVector, x) = foldr(renamespace, names, init = x)
431431
function renamespace(sys, x)
@@ -434,7 +434,7 @@ function renamespace(sys, x)
434434
if x isa Symbolic
435435
T = typeof(x)
436436
if istree(x) && operation(x) isa Operator
437-
return similarterm(x, operation(x),
437+
return similarterm(x, renamespace.(sys,operation(x)),
438438
Any[renamespace(sys, only(arguments(x)))])::T
439439
end
440440
let scope = getmetadata(x, SymScope, LocalScope())
@@ -495,18 +495,18 @@ function namespace_expr(O, sys, n = nameof(sys))
495495
O = unwrap(O)
496496
if any(isequal(O), ivs)
497497
return O
498-
elseif isvariable(O)
499-
renamespace(n, O)
500498
elseif istree(O)
501499
T = typeof(O)
502-
if symtype(operation(O)) <: FnType
503-
renamespace(n, O)::T
500+
renamed = let sys = sys, n = n, T = T
501+
map(a -> namespace_expr(a, sys, n)::Any, arguments(O))
502+
end
503+
if isvariable(O)
504+
similarterm(O, renamespace(n, operation(O)), renamed)::T
504505
else
505-
renamed = let sys = sys, n = n, T = T
506-
map(a -> namespace_expr(a, sys, n)::Any, arguments(O))
507-
end
508506
similarterm(O, operation(O), renamed)::T
509507
end
508+
elseif isvariable(O)
509+
renamespace(n, O)
510510
elseif O isa Array
511511
let sys = sys, n = n
512512
map(o -> namespace_expr(o, sys, n), O)

0 commit comments

Comments
 (0)