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

**URL:** <https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086>\
**Category:** Modelling & Simulations\
**Tags:** modelingtoolkit\
**Created:** [May 12, 2026, 8:24am UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086 "2026-05-12T08:24:44Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![JaviJR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/javijr/32/218294_2.png) [@JaviJR](https://discourse.julialang.org/u/JaviJR)\
**Post date:** [May 12, 2026, 8:24am UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/1 "2026-05-12T08:24:44Z")

</div>

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](https://discourse.julialang.org/t/ideal-diode-and-ideal-switch-instability-problem/107769/9) and [here](https://discourse.julialang.org/t/modeling-piecewise-linear-diodes-discontinuous-behaviors-with-dyad/130058/5)), 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](https://global.discourse-cdn.com/julialang/original/3X/9/c/9cfd9bc1b602af6a0f98e360b521540041e57eb2.png)

```julia-auto
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](https://global.discourse-cdn.com/julialang/original/3X/6/a/6ad68bbca00da00f2bc790ade8cdc6b8d8018055.png)

```julia-auto
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](https://global.discourse-cdn.com/julialang/original/3X/2/c/2cfbbc1747216a179aa6df7f23c599d3d3cf0129.png)

```julia-auto
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](https://global.discourse-cdn.com/julialang/original/3X/c/8/c80489103f6b960182880545d0272502c2889f40.png)

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!

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [May 12, 2026, 10:45am UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/2 "2026-05-12T10:45:53Z")

</div>

> [@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?

---

<div class="post-metadata">

**Author:** ![JaviJR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/javijr/32/218294_2.png) [@JaviJR](https://discourse.julialang.org/u/JaviJR)\
**Post date:** [May 12, 2026, 10:52am UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/3 "2026-05-12T10:52:22Z")

</div>

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 ☹ .

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [May 12, 2026, 10:55am UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/4 "2026-05-12T10:55:09Z")

</div>

You can also try using `guesses`, see [Initialization of Systems · ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/dev/tutorials/initialization/)

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

---

<div class="post-metadata">

**Author:** ![JaviJR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/javijr/32/218294_2.png) [@JaviJR](https://discourse.julialang.org/u/JaviJR)\
**Post date:** [May 12, 2026, 11:14am UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/5 "2026-05-12T11:14:57Z")

</div>

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.

```julia-auto
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)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/b/eb801ea0a19293862d139b1c3c39e25ab6978a3c.png)

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [May 12, 2026, 11:35am UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/6 "2026-05-12T11:35:36Z")

</div>

> [@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`

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 12, 2026, 12:03pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/7 "2026-05-12T12:03:53Z")

</div>

Is this a higher index DAE at the initial point?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 12, 2026, 12:04pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/8 "2026-05-12T12:04:33Z")

</div>

> [@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.

---

<div class="post-metadata">

**Author:** ![JaviJR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/javijr/32/218294_2.png) [@JaviJR](https://discourse.julialang.org/u/JaviJR)\
**Post date:** [May 12, 2026, 12:06pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/9 "2026-05-12T12:06:31Z")

</div>

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`.

```julia-auto
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

```

---

<div class="post-metadata">

**Author:** ![JaviJR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/javijr/32/218294_2.png) [@JaviJR](https://discourse.julialang.org/u/JaviJR)\
**Post date:** [May 12, 2026, 12:10pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/10 "2026-05-12T12:10:14Z")

</div>

> [@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).

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 12, 2026, 12:17pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/11 "2026-05-12T12:17:24Z")

</div>

> [@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.

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [May 12, 2026, 12:27pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/12 "2026-05-12T12:27:46Z")

</div>

> [@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 🙂

---

<div class="post-metadata">

**Author:** ![JaviJR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/javijr/32/218294_2.png) [@JaviJR](https://discourse.julialang.org/u/JaviJR)\
**Post date:** [May 12, 2026, 12:30pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/13 "2026-05-12T12:30:45Z")

</div>

* * *

With the basic Shockley model:

```julia-auto
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`):

```julia-auto
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

```

---

<div class="post-metadata">

**Author:** ![JaviJR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/javijr/32/218294_2.png) [@JaviJR](https://discourse.julialang.org/u/JaviJR)\
**Post date:** [May 12, 2026, 12:32pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/14 "2026-05-12T12:32:16Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [May 12, 2026, 12:42pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/15 "2026-05-12T12:42:42Z")

</div>

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

```julia
@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

```

---

<div class="post-metadata">

**Author:** ![pierre-haessig](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pierre-haessig/32/217129_2.png) [@pierre-haessig](https://discourse.julialang.org/u/pierre-haessig)\
**Post date:** [May 12, 2026, 1:28pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/16 "2026-05-12T13:28:58Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![JaviJR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/javijr/32/218294_2.png) [@JaviJR](https://discourse.julialang.org/u/JaviJR)\
**Post date:** [May 12, 2026, 1:48pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/17 "2026-05-12T13:48:07Z")

</div>

> [@langestefan](#):
>
> Can you try a ramped AC source to get around the initialization problem?
> 
> ```julia-auto
> @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:

```julia-auto
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])

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/a/a/aadeaef55c5dd4abf3172ea6144fea8868edd74b.png)

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

---

<div class="post-metadata">

**Author:** ![JaviJR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/javijr/32/218294_2.png) [@JaviJR](https://discourse.julialang.org/u/JaviJR)\
**Post date:** [May 12, 2026, 1:55pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/18 "2026-05-12T13:55:12Z")

</div>

> [@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:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/b/2bc4c9d7f6a292ea3592b1c02747e84020ad99d8.png)

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.

---

<div class="post-metadata">

**Author:** ![pierre-haessig](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pierre-haessig/32/217129_2.png) [@pierre-haessig](https://discourse.julialang.org/u/pierre-haessig)\
**Post date:** [May 12, 2026, 3:58pm UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/19 "2026-05-12T15:58:27Z")

</div>

> [@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

---

<div class="post-metadata">

**Author:** ![JaviJR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/javijr/32/218294_2.png) [@JaviJR](https://discourse.julialang.org/u/JaviJR)\
**Post date:** [May 13, 2026, 10:02am UTC](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086/20 "2026-05-13T10:02:08Z")

</div>

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

 ![dibujo4](https://global.discourse-cdn.com/julialang/original/3X/e/c/ec203e8d3d626d6c30e984c57a4b92ff8f96f6a7.png)

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

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/6/0613e11a70ccd5718cedb1b65e0eaaf0b91f909a.png)

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

```julia-auto
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)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/c/7c10352a72149d51daf9a683fbebced66f020662.jpeg)

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:

```julia-auto
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.

```julia-auto
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),
]

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/a/5/a5fdd50ba825c4771712675b30aba68a85a0007e.png)  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/1/d1f4c34a20915d52790d28b3448139332cf169a2.png)

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.

[Next page](https://discourse.julialang.org/t/ideal-diode-switching-with-multiple-voltage-sources-in-modelingtoolkit-solver-gets-stuck-at-t-0/137086.md?page=2)