Ideal diode switching with multiple voltage sources in ModelingToolkit , solver gets stuck at t=0

Specific Domains Modelling & Simulations
1

Hi everyone,

My research team and I have been working with Julia and ModelingToolkit for a while now and have obtained quite satisfactory results overall (our focus is on wind turbine modeling). Recently, however, we started exploring whether ModelingToolkit could be used for EMT-type (Electromagnetic Transient) simulations, where switching elements (diodes, thyristors) need to be modeled as realistically as possible.

Looking through the forum, we found that this topic has been discussed before (see here and here), and the implementations described there appeared to be working. So we ran some tests ourselves.

What works

Circuit 1 — Parallel rectifier (single source, two diodes)
dibujo1

function ParallelRectifierCircuit(; name,
                                    RMSVoltage = 9.0,
                                    Freq       = 50.0,
                                    Phase      = 0.0,
                                    R1         = 10.0,
                                    R2         = 10.0,
                                    I_S        = 1e-6,
                                    n_ideal    = 1.0,
                                    T_K        = 293.15)
    @named source = SourceAC(RMSVoltage = RMSVoltage, Phase = Phase, Freq = Freq)
    @named diode1 = ShockleyDiode(I_S = I_S, n_ideal = n_ideal, T_K = T_K)
    @named diode2 = ShockleyDiode(I_S = I_S, n_ideal = n_ideal, T_K = T_K)
    @named res1   = Resistor(R = R1)
    @named res2   = Resistor(R = R2)
    @named gnd    = Ground()
    eqs = [
        connect(source.p, diode1.p, diode2.p),
        connect(diode1.n, res1.p),
        connect(diode2.n, res2.p),
        connect(res1.n,   res2.n, gnd.g),
        connect(source.n, gnd.g),
    ]
    System(eqs, t; systems = [source, diode1, diode2, res1, res2, gnd], name)
end

@named circuit1 = ParallelRectifierCircuit(R1 = 10.0, R2 = 20.0)
sys1  = mtkcompile(circuit1)
prob1 = ODEProblem(sys1, [], (0.0, 0.1))
sol1  = solve(prob1, Rodas5(), dt = 1e-5, adaptive = false)
plot(sol1, idxs = [circuit1.diode1.p.i, circuit1.diode2.p.i])

This solves correctly and produces the expected waveforms.

Circuit 2 — Center-tap rectifier (single source, two diodes)

dibujo2

function CenterTapRectifier(; name,
                               RMSVoltage = 9.0,
                               Freq       = 50.0,
                               R_load     = 10.0,
                               I_S        = 1e-6,
                               n_ideal    = 1.0,
                               T_K        = 293.15)
    @named src   = SourceAC(RMSVoltage = RMSVoltage, Freq = Freq, Phase = 0.0)
    @named D1    = ShockleyDiode(I_S = I_S, n_ideal = n_ideal, T_K = T_K)
    @named D2    = ShockleyDiode(I_S = I_S, n_ideal = n_ideal, T_K = T_K)
    @named rload = Resistor(R = R_load)
    @named gnd   = Ground()
    eqs = [
        connect(src.p,  D1.p),
        connect(D1.n,   rload.p),
        connect(src.n,  D2.p),
        connect(D2.n,   D1.p),
        connect(src.n,  rload.n, gnd.g),
    ]
    System(eqs, t; systems = [src, D1, D2, rload, gnd], name)
end

@named ctr = CenterTapRectifier(RMSVoltage = 9.0, Freq = 50.0, R_load = 10.0)
sys3  = mtkcompile(ctr)
prob3 = ODEProblem(sys3, [], (0.0, 0.04))
sol3  = solve(prob3, Rodas5(), dt = 1e-5, adaptive = false)
plot(sol3, idxs = [ctr.rload.p.i])

The problem

Circuit 3 — Two AC sources (120° phase shift), four diodes

This is essentially two center-tap rectifiers sharing the same DC bus, fed by two sources shifted by 2π/3:
dibujo3

function Converter2(; name,
                      RMSVoltage = 9.0,
                      Freq       = 50.0,
                      R_load     = 10.0,
                      I_S        = 1e-6,
                      n_ideal    = 1.0,
                      T_K        = 293.15)
    @named srcA  = SourceAC(RMSVoltage = RMSVoltage, Freq = Freq, Phase = 0.0)
    @named srcB  = SourceAC(RMSVoltage = RMSVoltage, Freq = Freq, Phase = -2π/3)
    @named D1    = ShockleyDiode(I_S = I_S, n_ideal = n_ideal, T_K = T_K)
    @named D2    = ShockleyDiode(I_S = I_S, n_ideal = n_ideal, T_K = T_K)
    @named D3    = ShockleyDiode(I_S = I_S, n_ideal = n_ideal, T_K = T_K)
    @named D4    = ShockleyDiode(I_S = I_S, n_ideal = n_ideal, T_K = T_K)
    @named rload = Resistor(R = R_load)
    @named gnd   = Ground()
    eqs = [
        connect(srcA.p, D1.p),
        connect(D1.n,   rload.p),
        connect(srcA.n, D2.p),
        connect(D2.n,   D1.p),
        connect(srcB.p, D3.p),
        connect(D3.n,   rload.p),
        connect(srcB.n, D4.p),
        connect(D4.n,   D3.p),
        connect(srcA.n, srcB.n, rload.n, gnd.g),
    ]
    System(eqs, t; systems = [srcA, srcB, D1, D2, D3, D4, rload, gnd], name)
end

@named ctr2 = Converter2(RMSVoltage = 9.0, Freq = 50.0, R_load = 10.0)
sys4  = mtkcompile(ctr2)
prob4 = ODEProblem(sys4, [], (0.0, 0.04))
sol4  = solve(prob4, Rodas5(), dt = 1e-5, adaptive = false)

The solver gets stuck at t = 0.0 and never advances:

duda
duda754×169 7.24 KB

We tried adding snubber capacitors in parallel with the diodes, which sometimes allows the solver to get past t = 0, but it makes the system extremely stiff and eventually causes other failures that prevent it from completing. We also tried forcing fixed-step integration with a small dt (adaptive = false, dt = 1e-5), but that made no difference either. Has anyone run into this with multiple interacting switching elements in MTK and found a reliable solution? Is there a recommended pattern for handling simultaneous switching events other than snubbers? Overall, this feels like a significant barrier to simulating power electronic converters and similar circuits in MTK, so any guidance would be greatly appreciated.

Thanks in advance!