SciML
UndefVarError of an irreducible variable
Thanks a ton to Yingbo, Chris and contributors for supporting this wondeful package!
I'm trying to build a simple circuit where an AC source is supplying two identical parallel circuits, each with an RL through a diode. The code passes when there is only one RL-diode circuit but fails with two identical ones in parallel.
The diode switching is implemented using a continuous event. The algebraic variable is marked irreducible=true.
using DifferentialEquations
using IfElse: ifelse
using ModelingToolkit
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Sine
@variables t
function Diode(; name, v_start=0.0, i_start=0.0, Rc=0.01, Rn=1e6, VFwdDrop=0.0)
@named p = Pin()
@named n = Pin()
@variables v(t) [irreducible=true]
@variables i(t) = i_start
eqs = [v ~ p.v - n.v
0 ~ p.i + n.i
i ~ p.i
i ~ ifelse(v >= VFwdDrop,
v / Rc, # on
v / Rn, # off when V <= VFwdDrop
)
]
compose(ODESystem(eqs, t, [v, i], [];
continuous_events = [v ~ VFwdDrop],
name=name), p, n)
end
"""
Single phase voltage source using sine function
V = Vm sin(ωt+ϕ)
"""
function ACVoltageSource(; name, Vrms = 1.0, ω = 2π, ϕ = 0)
@named oneport = OnePort()
@unpack v, i = oneport
ps = @parameters Vrms=Vrms ω=ω ϕ=ϕ
eqs = [
v ~ sqrt(2) * Vrms * sin(ω * t + ϕ)
]
extend(ODESystem(eqs, t, [], ps; name), oneport)
end
@named source_ac = ACVoltageSource(Vrms = 230, ω = 60*2π, ϕ = 0)
@named r_source = Resistor(R = 0.001)
@named ground = Ground()
@named resistor = Resistor(R = 100)
@named inductor = Inductor(L = 500e-3)
@named diode = Diode()
@named resistor2 = Resistor(R = 100)
@named inductor2 = Inductor(L = 500e-3)
@named diode2 = Diode()
connections = [
connect(source_ac.p, r_source.n)
connect(r_source.p, resistor.p)
connect(resistor.n, inductor.p)
connect(inductor.n, diode.p)
connect(diode.n, source_ac.n, ground.g)
connect(r_source.p, resistor2.p)
connect(resistor2.n, inductor2.p)
connect(inductor2.n, diode2.p)
connect(diode2.n, source_ac.n, ground.g)
]
@named model = ODESystem(connections, t;
systems = [
ground,
r_source,
source_ac,
resistor,
inductor,
diode,
resistor2,
inductor2,
diode2,
])
sys = structural_simplify(model)
states(sys)
gives
4-element Vector{Any}:
inductor₊i(t)
inductor2₊i(t)
diode₊v(t)
diode2₊v(t)
prob = ODAEProblem(sys, Pair[], (0.0, 0.06))
@assert isirreducible(sys.diode₊v)
@assert isirreducible(sys.diode2₊v)
integrator = init(prob,
abstol=1e-4,
reltol=1e-6,
)
@time solve!(integrator)
gives
UndefVarError: `diode₊v(t)` not defined
Stacktrace:
[1] macro expansion
@ ~/.julia/packages/SymbolicUtils/H684H/src/code.jl:352 [inlined]
[2] macro expansion
@ ~/.julia/packages/RuntimeGeneratedFunctions/b1u1T/src/RuntimeGeneratedFunctions.jl:141 [inlined]
[3] macro expansion
@ ./none:0 [inlined]
[4] generated_callfunc
@ ./none:0 [inlined]
The following code implements the AC source supplying the RL through the diode. It works fine.
connections = [
connect(source_ac.p, r_source.n)
connect(r_source.p, resistor.p)
connect(resistor.n, inductor.p)
connect(inductor.n, diode.p)
connect(diode.n, source_ac.n, ground.g)
# connect(r_source.p, resistor2.p)
# connect(resistor2.n, inductor2.p)
# connect(inductor2.n, diode2.p)
# connect(diode2.n, source_ac.n, ground.g)
]
@named model = ODESystem(connections, t;
systems = [
ground,
r_source,
source_ac,
resistor,
inductor,
diode,
# resistor2,
# inductor2,
# diode2,
])
sys = structural_simplify(model)
prob = ODAEProblem(sys, Pair[], (0.0, 0.06))
# @assert isirreducible(sys.diode₊v)
# @assert isirreducible(sys.diode2₊v)
integrator = init(prob,
abstol=1e-4,
reltol=1e-6,
)
@time solve!(integrator)
Thank you so much in advance.
I forgot to mention the version information:
Julia Version 1.9.1
Commit 147bdf428cd (2023-06-07 08:27 UTC)
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 64 × AMD Ryzen Threadripper PRO 5975WX 32-Cores
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-14.0.6 (ORCJIT, znver3)
Threads: 2 on 64 virtual cores
[6e4b80f9] BenchmarkTools v1.3.2
[0c46a032] DifferentialEquations v7.8.0
[615f187c] IfElse v0.1.1
[961ee093] ModelingToolkit v8.59.1
[16a59e39] ModelingToolkitStandardLibrary v1.17.0
[91a5bcdd] Plots v1.38.16
I was browsing other open issues and realized that this may be the same issue as https://github.com/SciML/ModelingToolkit.jl/issues/2173.
Just a quick update on this issue.
MTK 8.60.0 can properly build the ODEProblem with
prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0.0, 1.0));
But solving the problem yields InitialFailure.
Thanks again for all developers' efforts in perfecting MTK.
A tangential question is related to the irreducibility of v of Diode. Since diode.v ~ p.v - n.v, does MTK automatically mark p.v and n.v as irreducible? If not, could this cause issues?
function Diode(; name, v_start=0.0, i_start=0.0, Rc=0.01, Rn=1e6, VFwdDrop=0.0)
@named p = Pin()
@named n = Pin()
@variables v(t) [irreducible=true] # <-- here
@variables i(t) = i_start
eqs = [v ~ p.v - n.v # <-- also here
0 ~ p.i + n.i
i ~ p.i
i ~ ifelse(v >= VFwdDrop,
v / Rc, # on
v / Rn, # off when V <= VFwdDrop
)
]
compose(ODESystem(eqs, t, [v, i], [];
continuous_events = [v ~ VFwdDrop],
name=name), p, n)
end
