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

[JaviJR]

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)

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)

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:

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:

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!

[langestefan]

@JaviJR

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

The problem may be initialization then. Can you try a non-zero initial phase like π/2 as the initial value for src A?

[JaviJR]

I already tried using a non-zero initial phase (e.g. π/2 for srcA and shifted accordingly for srcB), but unfortunately the solver still gets stuck at t = 0.0 and never advances .

[langestefan]

You can also try using guesses, see Initialization of Systems · ModelingToolkit.jl

And what may also help is a small capacitor from DC to ground. And as a solver you may also want to try Rodas5P.

[JaviJR]

Hi @langestefan, thanks for the suggestions. We tried both, adding a small capacitor from the DC bus to ground and switching to Rodas5P, and also passed initial guesses for the diode and capacitor voltages, but unfortunately the solver still gets stuck at t = 0 and we can’t get it to advance.

Honestly, we are not sure at this point whether this is a fundamental limitation of how MTK handles simultaneous switching events, or whether we are simply missing the right way to set things up. Tools like PSCAD or PSpice have dedicated teams and decades of development behind their switching solvers, so it is understandable that they handle this kind of circuit more robustly, but it is still disappointing to hit this wall in MTK, especially given how well everything else has worked for us so far.

function Converter2(; name,
                      RMSVoltage = 9.0,
                      Freq       = 50.0,
                      R_load     = 10.0,
                      C_snub     = 1e-6,
                      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 cdc   = Capacitor(C = C_snub)   
    @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),
        connect(rload.p, cdc.p),  
        connect(cdc.n,   gnd.g),
    ]
    System(eqs, t; systems = [srcA, srcB, D1, D2, D3, D4, rload, cdc, gnd], name)
end

@named ctr2 = Converter2(RMSVoltage = 9.0, Freq = 50.0, R_load = 10.0, C_snub = 1e-6)
sys4  = mtkcompile(ctr2)


guesses = [
    ctr2.D1.v => 0.0,
    ctr2.D2.v => 0.0,
    ctr2.D3.v => 0.0,
    ctr2.D4.v => 0.0,
    ctr2.cdc.v => 0.0,
]

prob4 = ODEProblem(sys4, [], (0.0, 0.04), guesses = guesses)
sol4  = solve(prob4, Rodas5P(), dt = 1e-5, adaptive = false)

[langestefan]

@JaviJR

Tools like PSCAD or PSpice have dedicated teams and decades of development behind their switching solvers, so it is understandable that they handle this kind of circuit more robustly

Well MTK is not a tool purpose built for doing circuit simulations. MTK can handle pretty much any kind of physical system you can think of. So you’re not making a fair comparison here.

@JaviJR

but it is still disappointing to hit this wall in MTK, especially given how well everything else has worked for us so far.

It’s more likely a problem with your model specifiation than MTK itself.

In this case, I think the diode implementation itself is the problem. My implementation was quite naive directly from the equations, so you probably want to look at the spice implementation where they must do something that is numerically stable.

You can see why this blows up:

With T = 300.15 K:

V_T = kT/q = 0.02587 V

So the equation is effectively i = Is·(exp(v / 0.02587) − 1).

At t = 0 with RMSVoltage = 9 (V̂ ≈ 12.73 V) and srcB phase = −2π/3:

• V_srcB(0) = 12.73 · sin(−2π/3) ≈ −11.02 V
• D4 sees v ≈ +11.02 V forward bias
• argument to exp = 11.02 / 0.02587 ≈ 426
• exp(426) ≈ 10^185
• i_D4 = 1e-6 · 10^185 ≈ 10^179 A

[ChrisRackauckas]

Is this a higher index DAE at the initial point?

[ChrisRackauckas]

@langestefan

i = Is·(exp(v / 0.02587) − 1)

SPICE simulators will usually truncate / temper this in order to make it more numerically stable in the tails.

[JaviJR]

Thanks for the detailed breakdown. Based on that we tried a SPICE-style linearization (we also tried to simply clamp the exponential, but it did not work), where instead of clamping the exponential we replace it with its tangent line beyond a linearization threshold VLIN, keeping the model continuous and differentiable everywhere while bounding the current growth to linear for large forward bias. Unfortunately it did not work either and the solver still gets stuck at t = 0.

function ShockleyDiode(;
    name,
    I_S     = 1e-12,
    n_ideal = 1.0,
    T_K     = 300.15,
    R_s     = 1e-2,
    G_min   = 1e-8,
    V_lin   = 0.8,
)
    @parameters begin
        I_s    = I_S
        n_id   = n_ideal
        T      = T_K
        Rs     = R_s
        GMIN   = G_min
        VLIN   = V_lin
        k_B    = 1.380649e-23
        q_e    = 1.602176634e-19
    end
    @variables begin
        v(t)  = 0.0
        vd(t) = 0.0
        i(t)  = 0.0
    end
    systems = @named begin
        p = Pin()
        n = Pin()
    end

    Vt = k_B * T / q_e

    I_lin  = I_s * (exp(VLIN / (n_id * Vt)) - 1)
    dI_lin = I_s * exp(VLIN / (n_id * Vt)) / (n_id * Vt)

    diode_current(vd) = ifelse(
        vd < VLIN,
        I_s * (exp(vd / (n_id * Vt)) - 1) + GMIN * vd,
        I_lin + dI_lin * (vd - VLIN) + GMIN * vd,
    )

    eqs = [
        v   ~ p.v - n.v,
        v   ~ vd + Rs * i,
        i   ~ diode_current(vd),
        p.i ~  i,
        n.i ~ -i,
    ]
    System(eqs, t; name, systems)
end

[JaviJR]

@ChrisRackauckas

Is this a higher index DAE at the initial point?

The system comes out as index 1 after mtkcompile, with 1 differential equation (the capacitor voltage, if a snubber capacitor is added) and 4 algebraic equations (the diode models).

[ChrisRackauckas]

@JaviJR

The system comes out as index 1 after mtkcompile, with 1 differential equation (the capacitor voltage, if a snubber capacitor is added) and 4 algebraic equations (the diode models).

Structurally index one, but is it numerically index 1? Share the initial Jacobian.

[langestefan]

@ChrisRackauckas

SPICE simulators will usually truncate / temper this in order to make it more numerically stable in the tails.

Indeed. And I know that I said it many times before, but if we had good base component libraries, we wouldn’t need DIY solutions

[JaviJR]

---

With the basic Shockley model:

1×1 Matrix{Float64}:
 -1.1934200746486273e191
rank: 1
cond: 1.0
det:  -1.1934200746486273e191

With the SPICE-style linearized model (tangent line beyond VLIN = 0.8 V):

5×5 Matrix{Float64}:
  0.0    0.0   0.0     1.0    0.0
 -0.0   -1.0  -0.0    -0.0   -0.0
 100.0   0.0   0.0     0.0  100.0
  0.0    0.0   1.0     0.0    1.0
 -1.0e8  0.0  39.5956  0.0  -1.001e8

rank: 5
cond: 2.4509517798668017e9
det:  1.0007934967303366e7

[JaviJR]

Totally agree. On top of that, we have been using ModelingToolkitStandardLibrary quite a bit, but since many of its components still rely on @mtkmodel (which still works but is deprecated) there are sometimes interactions that do not behave quite right. So even when the building blocks exist, you end up hitting rough edges anyway.

[langestefan]

Can you try a ramped AC source to get around the initialization problem?

@mtkmodel SourceACRamped begin
    @extend v, i = oneport = OnePort()
    @parameters begin
        RMSVoltage = 1.0
        Freq       = 50.0
        Phase      = 0.0
        τ          = 2e-3, [description = "Ramp time constant (s)"]
    end
    @equations begin
        v ~ (1 - exp(-t / τ)) *
            RMSVoltage * sqrt(2) * sin(2π * Freq * t + Phase)
    end
end

[pierre-haessig]

Thanks @JaviJR for sharing your experiments with diode rectifiers using MTK. I’ve not yet played with it myself.

One thing I’m missing though in your examples: are the test Circuits 1/2/3 “just” test circuits (that is fragments or larger circuits like a three-phase rectifier with 6 diodes) or are they your actual target applications?

I’m asking because I’m confused by the topology of Circuit 2 (and thus Circuit 3 as well). To me it seems that the diode D2 will put the AC voltage source to near short circuit every half period, whenever the AC source is meant to be negative. Assuming that the AC source is an ideal voltage source (0 series resistance) while the diode has a continuous characteristic (i ∝ exp(v)), this means that the D2 voltage drop will be forced to a given positive voltage value of RMSVoltage*sqrt(2) ≈ 13 V, yielding a huge forward current (which is what I called “near short circuit”). Is my analysis correct?

(now perhaps this is orthogonal to the problem you are facing, but at least let’s be sure to eliminate potentially ill-defined operating points)

[Edit]: if Circuit 2 is meant to a “one leg extract” of a three-phase (6 diodes) rectifier, I would recommand to modify it to connect the “-” side of the voltage source to the mid-point of the DC side, that is cut Rload into two half resistors to create a mid-point node.

[JaviJR]

@langestefan

Can you try a ramped AC source to get around the initialization problem?

@mtkmodel SourceACRamped begin
    @extend v, i = oneport = OnePort()
    @parameters begin
        RMSVoltage = 1.0
        Freq       = 50.0
        Phase      = 0.0
        τ          = 2e-3, [description = "Ramp time constant (s)"]
    end
    @equations begin
        v ~ (1 - exp(-t / τ)) *
            RMSVoltage * sqrt(2) * sin(2π * Freq * t + Phase)
    end
end

It seems to start simulating with this approach, though I still need to run more tests to check whether it scales to more complex topologies. For reproducibility, here is the full code and the resulting waveform:

function ShockleyDiode(; name, I_S = 1e-6, n_ideal = 1.0, T_K = 293.15)

    @parameters begin
        I_s   = I_S
        n_id  = n_ideal
        T     = T_K
        k_B   = 1.380649e-23   # Constante de Boltzmann [J/K]
        q_e   = 1.602e-19      # Carga del electrón [C]
    end
    @variables begin
        v(t) =0.0
        i(t) 
    end
    systems = @named begin
        p = Pin()
        n = Pin()
    end

    eqs = [
        v   ~ p.v - n.v,
        i   ~ I_s * (exp(v / (n_id * k_B * T / q_e)) - 1),
        p.i ~  i,
        n.i ~ -i,
    ]
    System(eqs, t; name, systems)
end

function SourceAC(; name,
                          RMSVoltage = 220.0,
                          Phase      = 0.0,
                          Freq       = 50.0,
                          τ          = 2e-3)
    @named p = Pin()
    @named n = Pin()
    @parameters RMSVoltage=RMSVoltage  Phase=Phase  Freq=Freq  τ=τ
    eqs = [
        p.v - n.v ~ (1 - exp(-t / τ)) * (RMSVoltage * √2) * sin(2π * Freq * t + Phase),
        p.i + n.i ~ 0
    ]
    System(eqs, t, [], [RMSVoltage, Phase, Freq, τ];
        systems = [p, n],
        name)
end

function Converter2(; name,
                      RMSVoltage = 9.0,
                      Freq       = 50.0,
                      R_load     = 10.0,
                      C_snub     = 1e-6,
                      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 cdc   = Capacitor(C = C_snub)   # condensador DC → gnd
    @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),
        connect(rload.p, cdc.p),   # DC bus → condensador
        connect(cdc.n,   gnd.g),
    ]
    System(eqs, t; systems = [srcA, srcB, D1, D2, D3, D4, rload, cdc, gnd], name)
end

@named ctr2 = Converter2(RMSVoltage = 9.0, Freq = 50.0, R_load = 10.0, C_snub = 1e-6)

sys4  = mtkcompile(ctr2)

prob4 = ODEProblem(sys4, [], (0.0, 0.04))

sol4  = solve(prob4, Rodas5P(), dt = 1e-5, adaptive = false)

plot(sol4, idxs = [ctr2.rload.p.i])

Thank you all A LOT. As soon as i do more tests I will let you know it this solution does actually scale well.

[JaviJR]

@pierre-haessig

Thanks @JaviJR for sharing your experiments with diode rectifiers using MTK. I’ve not yet played with it myself.

One thing I’m missing though in your examples: are the test Circuits 1/2/3 “just” test circuits (that is fragments or larger circuits like a three-phase rectifier with 6 diodes) or are they your actual target applications?

I’m asking because I’m confused by the topology of Circuit 2 (and thus Circuit 3 as well). To me it seems that the diode D2 will put the AC voltage source to near short circuit every half period, whenever the AC source is meant to be negative. Assuming that the AC source is an ideal voltage source (0 series resistance) while the diode has a continuous characteristic (i ∝ exp(v)), this means that the D2 voltage drop will be forced to a given positive voltage value of RMSVoltage*sqrt(2) ≈ 13 V, yielding a huge forward current (which is what I called “near short circuit”). Is my analysis correct?

(now perhaps this is orthogonal to the problem you are facing, but at least let’s be sure to eliminate potentially ill-defined operating points)

[Edit]: if Circuit 2 is meant to a “one leg extract” of a three-phase (6 diodes) rectifier, I would recommand to modify it to connect the “-” side of the voltage source to the mid-point of the DC side, that is cut Rload into two half resistors to create a mid-point node.

That is a very valid point, thanks for catching it. Circuits 1, 2 and 3 were just toy examples to explore how MTK handles diode switching incrementally, not physically accurate designs. The actual target is a full three-phase rectifier like this one:

Eventually we also want to extend this to IGBT-based topologies to simulate more general power electronics systems, but one step at a time. The broader goal behind all of this is to explore whether EMT simulation of electrical equipment can be done efficiently using ModelingToolkit, if that turns out to be feasible at scale, it would be pretty mindblowing for the electrical engineering field. As we make more progress and run more tests we will share the results here.

[pierre-haessig]

@JaviJR

The actual target is a full three-phase rectifier like this one:

OK then my proposition to rewire Circuit 2 with the negative pin of the AC source to a newly created mid point of the DC side is a meaninful simplification.

From my experience with the simulation of rectifiers, I would suggest testing all the cross product combinations of the following variants of AC and DC side impedance:

• AC side:
  • ideal voltage source
  • voltage source + small resistance in series
  • voltage source + small inductances in series (e.g. grid or transformer leakage)
• DC side :
  • only resistance
  • DC : resistance + capacitive filter in // (i.e. voltage source like behavior)
  • DC : smoothing inductive in series with load (R or R//C, then it doesn’t matter), i.e. current source like behaviour

[JaviJR]

Following the advice, we have been running more tests. We implemented the following circuit (two sources at 0° and 180°):

To evaluate the results we implemented the same system in both Simulink and ModelingToolkit. In Simulink the load current looks as expected:

However, in ModelingToolkit the solver detects instabilities and stops prematurely:

function ShockleyDiode(; name, I_S = 1e-6, n_ideal = 1.0, T_K = 293.15)
    @parameters begin
        I_s   = I_S
        n_id  = n_ideal
        T     = T_K
        k_B   = 1.380649e-23
        q_e   = 1.602e-19
    end
    @variables begin
        v(t) = 0.0
        i(t)
    end
    systems = @named begin
        p = Pin()
        n = Pin()
    end
    eqs = [
        v   ~ p.v - n.v,
        i   ~ I_s * (exp(v / (n_id * k_B * T / q_e)) - 1),
        p.i ~  i,
        n.i ~ -i,
    ]
    System(eqs, t; name, systems)
end

function SourceAC(; name, RMSVoltage = 220, Phase = 0, Freq = 50)
    @named p = Pin()
    @named n = Pin()
    @parameters RMSVoltage=RMSVoltage  Phase=Phase  Freq=Freq
    eqs = [
        p.v - n.v ~ (RMSVoltage * √2) * sin(2π * Freq * t + Phase),
        p.i + n.i ~ 0
    ]
    System(eqs, t, [], [RMSVoltage, Phase, Freq]; systems = [p, n], name)
end

function Circuit4(; 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 = π)
    @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 gnd   = Ground()
    @named cdc   = Capacitor(C = 1e-6)
    @named rload = Resistor(R = R_load)
    eqs = [
        connect(srcA.p, D1.p, D2.n),
        connect(srcB.p, D3.p, D4.n),
        connect(rload.p, D1.n, D3.n),
        connect(rload.n, D2.p, D4.p),
        connect(srcA.n, srcB.n, gnd.g),
    ]
    System(eqs, t; systems = [srcA, srcB, D1, D2, D3, D4, gnd, rload, cdc], name)
end

@named circuit4 = Circuit4(RMSVoltage = 9.0, Freq = 50.0, R_load = 100.0)
sys4  = mtkcompile(circuit4)
prob4 = ODEProblem(sys4, [], (0.0, 0.1))
sol4  = solve(prob4, Rodas5P(), dt = 1e-5, adaptive = false)

We also tried adding series resistances and inductances on the AC side, a capacitor in parallel with the load, the ramped source from earlier, and the softer diode model, but none of them made a difference:

function ShockleyDiodeSOFT(;
    name,
    I_S     = 1e-12,
    n_ideal = 1.0,
    T_K     = 300.15,
    R_s     = 1e-2,
    G_min   = 1e-8,
    V_lin   = 0.8,
)
    @parameters begin
        I_s    = I_S
        n_id   = n_ideal
        T      = T_K
        Rs     = R_s
        GMIN   = G_min
        VLIN   = V_lin
        k_B    = 1.380649e-23
        q_e    = 1.602176634e-19
    end
    @variables begin
        v(t)  = 0.0
        vd(t) = 0.0
        i(t)  = 0.0
    end
    systems = @named begin
        p = Pin()
        n = Pin()
    end
    Vt = k_B * T / q_e
    I_lin  = I_s * (exp(VLIN / (n_id * Vt)) - 1)
    dI_lin = I_s * exp(VLIN / (n_id * Vt)) / (n_id * Vt)
    diode_current(vd) = ifelse(
        vd < VLIN,
        I_s * (exp(vd / (n_id * Vt)) - 1) + GMIN * vd,
        I_lin + dI_lin * (vd - VLIN) + GMIN * vd,
    )
    eqs = [
        v   ~ p.v - n.v,
        v   ~ vd + Rs * i,
        i   ~ diode_current(vd),
        p.i ~  i,
        n.i ~ -i,
    ]
    System(eqs, t; name, systems)
end

function SourceACSOFT(; name,
                        RMSVoltage = 220.0,
                        Phase      = 0.0,
                        Freq       = 50.0,
                        τ          = 2e-3)
    @named p = Pin()
    @named n = Pin()
    @parameters RMSVoltage=RMSVoltage  Phase=Phase  Freq=Freq  τ=τ
    eqs = [
        p.v - n.v ~ (1 - exp(-t / τ)) * (RMSVoltage * √2) * sin(2π * Freq * t + Phase),
        p.i + n.i ~ 0
    ]
    System(eqs, t, [], [RMSVoltage, Phase, Freq, τ]; systems = [p, n], name)
end

There is a tempting workaround which is to tie the source neutral directly to the DC negative bus, as it makes the system compile and run correctly and the results match Simulink when the same ground is used on the load side. However, as pointed out earlier in this thread, this is physically wrong because the lower diodes end up short-circuiting the sources during the negative half cycles, producing unrealistically large instantaneous currents.

eqs = [
    connect(srcA.p, D1.p, D2.n),
    connect(srcB.p, D3.p, D4.n),
    connect(rload.p, D1.n, D3.n),
    connect(rload.n, D2.p, D4.p, gnd.g),
    connect(srcA.n, srcB.n, gnd.g),
]

At this point we are not sure whether this is a fundamental limitation of MTK, whether tools like Simulink or PSpice are doing something under the hood to handle this case, or whether there is simply a bug somewhere in our code that we have not spotted. Any ideas would be greatly appreciated.