# Modeling a crowbar circuit with ModelingToolkit.jl

**URL:** <https://discourse.julialang.org/t/modeling-a-crowbar-circuit-with-modelingtoolkit-jl/81968>\
**Category:** Modelling & Simulations\
**Tags:** modelingtoolkit\
**Created:** [May 31, 2022, 8:11am UTC](https://discourse.julialang.org/t/modeling-a-crowbar-circuit-with-modelingtoolkit-jl/81968 "2022-05-31T08:11:44Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Irene\_Ndindabahizi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/irene_ndindabahizi/32/36782_2.png) [@Irene\_Ndindabahizi](https://discourse.julialang.org/u/Irene_Ndindabahizi)\
**Post date:** [May 31, 2022, 8:11am UTC](https://discourse.julialang.org/t/modeling-a-crowbar-circuit-with-modelingtoolkit-jl/81968/1 "2022-05-31T08:11:44Z")

</div>

Hi All,  
I am trying to model an electromagnetic launcher using ModelingToolkit.jl. The electric circuit is basically an RLC with a crowbar circuit that cut off the capacitor when the voltage is 0V. I managed to do the simulation using callbacks (terminating the first simulation-sol1- at 0V and starting another RL circuit simulation-sol2).  
The current generated by this circuits is used for accelerating an object so I need to add the currents from the two solutions (sol1.i + sol2.i).  
The end goal of my simulations is to optimize the launcher by varying the number of capacitor banks. I would like to do it automatically by generating the code based one the number of capacitors banks.  
I tried to model the crowbar circuit as “idealCommutingSwitch” from Modelica but I didn’t manage to make it work (I have an unballanced system). Basically now I am trying to figure out how can I monitor the voltage and once the voltage hits the threshold how can I modify the circuits. Is this something possible in ModelingToolkit.jl?

Thank you for your help.

Here is a my code with a single capacitor bank:

```julia
using ModelingToolkit, Plots, DifferentialEquations

@variables t
@connector function Pin(;name)
    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
    ODESystem(Equation[], t, sts, []; name=name)
end

function Ground(;name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name=name), g)
end

function OnePort(;name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

function Resistor(;name, R = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R=R
    eqs = [
           v ~ i * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters C=C
    D = Differential(t)
    eqs = [
           D(v) ~ i / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Inductor(; name, L=1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters L=L
    D = Differential(t)
    eqs = [
        D(i) ~ v / L
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function ConstantVoltage(;name, V = 1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V=V
    eqs = [
           V ~ v
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

L1 = 0.6*1e-6; L2 = 4*1e-6
L3 = 1.1*1e-6; R1= 4*1e-3; R2 = 0.28*1e-3; R3 = 3.6*1e-3
RD = 0.1*1e-3; Lpr = 4.2*1e-9
C = 3.08*1e-3 ; Uc = 10*1e3

L = L1+L2+L3

R = R1+R2+R3

@named resistor_1 = Resistor(R=R1)
@named resistor_2 = Resistor(R=R2)
@named resistor_3 = Resistor(R=R3)
@named capacitor = Capacitor(C=C)
@named inductor_1 = Inductor(L=L1)
@named inductor_2 = Inductor(L=L2)
@named inductor_3 = Inductor(L=L3)
@named resistor_D = Resistor(R=RD)
#@named diode = Gate(false)

@named ground = Ground()

rc_eqs = [
        #connect_pins(V1.p, R1.p)
        connect_pins(resistor_1.n, inductor_1.p)
        #connect_pins(inductor_1.n, capacitor.p)
        connect_pins(capacitor.n, resistor_1.p, ground.g)
        connect_pins(inductor_1.n, resistor_2.p)
        #connect_pins(resistor_2.n, capacitor.p)
        connect_pins(resistor_2.n, inductor_2.p)
        connect_pins(inductor_2.n, resistor_3.p)
        connect_pins(resistor_3.n, inductor_3.p)
        connect_pins(inductor_3.n, capacitor.p)
        #connect_pins(capacitor.n, resistor_1.p)
        #connect_pins(V1.n, L2.n, GND.g)
         ]

rc_eqs2 = [
        
        
        #connect_pins(inductor_3.n, resistor_D.p, ground.g)
        connect_pins(resistor_D.n, resistor_2.p)
        #connect_pins(resistor_2.n, capacitor.p)
        connect_pins(resistor_2.n, inductor_2.p)
        connect_pins(inductor_2.n, resistor_3.p)
        connect_pins(resistor_3.n, inductor_3.p)
        connect_pins(inductor_3.n, resistor_D.p, ground.g)
        
]
@named _rc_model = ODESystem(rc_eqs, t)
@named _rc_model2 = ODESystem(rc_eqs2, t)

@named rc_model = compose(_rc_model,
[resistor_1,resistor_2,resistor_3, capacitor, inductor_1, inductor_2, inductor_3,ground])

@named rc_model2 = compose(_rc_model2,
[resistor_2,resistor_3, inductor_2, inductor_3,ground,resistor_D])

sys = structural_simplify(rc_model)

sys2 = structural_simplify(rc_model2)

u0 = [
        capacitor.v => -Uc
        inductor_1.i => 0.0
        inductor_2.i => 0.0
        inductor_3.i => 0.0
        inductor_2.v => 0.0
        inductor_3.v => 0.0
       ]

prob = ODEProblem(sys, u0, (0, 0.01))
function condition(u,t,integrator)
    -(integrator.sol[capacitor.v][end]+integrator.sol[resistor_1.v][end]+integrator.sol[inductor_1.v][end])<0.0
end

affect!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition,affect!, save_positions=(false,false))

sol = solve(prob, Rodas4(), reltol=1e-12,abstol=1e-12,callback=cb)

u02 = [
           inductor_2.i => sol[inductor_2.i][end]
           inductor_3.i => sol[inductor_3.i][end]
           inductor_3.v => sol[inductor_3.v][end]
          ]

prob2 = ODEProblem(sys2, u02, (sol.t[end], 0.01))
sol2 = solve(prob2, Rodas4(),reltol=1e-12,abstol=1e-12)

#plot(sol2)

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 31, 2022, 2:54pm UTC](https://discourse.julialang.org/t/modeling-a-crowbar-circuit-with-modelingtoolkit-jl/81968/2 "2022-05-31T14:54:29Z")

</div>

Hello and welcome to the community!

I’m struggling to grasp exactly what you mean by

> [@Irene\_Ndindabahizi](#):
>
> how can I monitor the voltage and once the voltage hits the threshold how can I modify the circuits

Generating `N` components in a for-loop should be straightforward, and you can also push connect statements to the vector of equations in a loop, like so:  
`push!(rc_eqs, connect_pins(resistor_n.n, inductor_n.p))`

Does this help?

---

<div class="post-metadata">

**Author:** ![Irene\_Ndindabahizi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/irene_ndindabahizi/32/36782_2.png) [@Irene\_Ndindabahizi](https://discourse.julialang.org/u/Irene_Ndindabahizi)\
**Post date:** [June 1, 2022, 7:13am UTC](https://discourse.julialang.org/t/modeling-a-crowbar-circuit-with-modelingtoolkit-jl/81968/3 "2022-06-01T07:13:54Z")

</div>

Thanks!  
This will certainly help when I start the optimization.

I am trying to simulate a current source generated by capacitor discharge. From t0, the voltage decreases and at a certain point it hits 0. At this point the capacitor needs to be short circuited from the load (which in my case is the launcher). In practice it is done through a crowbar circuit (see the attached picture for a schematic of the circuit).  
 ![circuit](https://global.discourse-cdn.com/julialang/original/3X/8/8/8864b38d8fb2898ad525124dad6559a024d3866c.jpeg)

Yesterday, I tried something new: I am now simulating the 2 complete circuit the resistor D (RD) begins with a high value (so very little current there) - this works fine - and as soon as the voltage there hits 0 I use callback to decrease the RD value to the correct value - for the moment the solve operation get stuck (or it gets really slow)  
Here I am trying to approximate switches using resistors … The best would really be to have 2 switches. I am still struggling to create a switch component in the circuit.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [June 1, 2022, 7:31am UTC](https://discourse.julialang.org/t/modeling-a-crowbar-circuit-with-modelingtoolkit-jl/81968/4 "2022-06-01T07:31:12Z")

</div>

If you could post the latest code you have that runs slowly, I could have a look at why. I’m not familiar with the electronics you are describing so I unfortunately can’t help much based on the textual description.

---

<div class="post-metadata">

**Author:** ![Irene\_Ndindabahizi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/irene_ndindabahizi/32/36782_2.png) [@Irene\_Ndindabahizi](https://discourse.julialang.org/u/Irene_Ndindabahizi)\
**Post date:** [June 1, 2022, 9:26am UTC](https://discourse.julialang.org/t/modeling-a-crowbar-circuit-with-modelingtoolkit-jl/81968/5 "2022-06-01T09:26:51Z")

</div>

Here is the code with the callbacks that changes the resistors parameter:

```julia
using ModelingToolkit, Plots, DifferentialEquations

@variables t
@connector function Pin(;name)
    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
    ODESystem(Equation[], t, sts, []; name=name)
end

function Ground(;name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name=name), g)
end

function OnePort(;name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

function Resistor(;name, R = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R=R
    eqs = [
           v ~ i * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters C=C
    D = Differential(t)
    eqs = [
           D(v) ~ i / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Inductor(; name, L=1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters L=L
    D = Differential(t)
    eqs = [
        D(i) ~ v / L
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function ConstantVoltage(;name, V = 1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V=V
    eqs = [
           V ~ v
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function VoltageSensor(; name)
    @named p = Pin()
    @named n = Pin()
    @variables v(t)=1.0
    eqs = [
        p.i ~ 0
        n.i ~ 0
        v ~ p.v - n.v
    ]
    ODESystem(eqs, t, [v], []; systems=[p, n], name=name)
end

L1 = 0.6*1e-6; L2 = 4*1e-6
L3 = 1.1*1e-6; R1= 4*1e-3; R2 = 0.28*1e-3; R3 = 3.6*1e-3
RD = 10.0#0.1*1e-3 #100.0 
Lpr = 4.2*1e-9
C = 3.08*1e-3 ; Uc = 10*1e3

L = L1+L2+L3

R = R1+R2+R3

@named resistor_1 = Resistor(R=R1)
@named resistor_2 = Resistor(R=R2)
@named resistor_3 = Resistor(R=R3)
@named capacitor = Capacitor(C=C)
@named inductor_1 = Inductor(L=L1)
@named inductor_2 = Inductor(L=L2)
@named inductor_3 = Inductor(L=L3)
@named resistor_D = Resistor(R=RD)
@named level = VoltageSensor()
#@named diode = Gate(false)

@named ground = Ground()

    rc_eqs =[
        #connect_pins(V1.p, R1.p)
        connect(resistor_1.n, inductor_1.p)
        #connect_pins(inductor_1.n, capacitor.p)
        connect(capacitor.n, resistor_1.p, ground.g)
        connect(inductor_1.n, resistor_2.p,resistor_D.n)
        #connect_pins(resistor_2.n, capacitor.p)
        connect(resistor_2.n, inductor_2.p)
        connect(inductor_2.n, resistor_3.p)
        connect(resistor_3.n, inductor_3.p)
        connect(inductor_3.n, capacitor.p,resistor_D.p)
        connect(inductor_1.n,level.p)
        connect(level.n,capacitor.p)

        #connect_pins(capacitor.n, resistor_1.p)
        #connect_pins(V1.n, L2.n, GND.g)
        #connect_pins(resistor_D.n, resistor_2.p)
        #connect_pins(inductor_3.n, resistor_D.p, ground.g)
         ]

@named _rc_model = ODESystem(rc_eqs, t)
#@named _rc_model2 = ODESystem(rc_eqs2, t)

@named rc_model = compose(_rc_model,
[resistor_1,resistor_2,resistor_3, capacitor, inductor_1, inductor_2, inductor_3,ground,resistor_D,level])

sys = structural_simplify(rc_model)

#sys2 = structural_simplify(rc_model2)

u0 = [
        capacitor.v => -Uc
        inductor_1.i => 0.0
        inductor_2.i => 0.0
        inductor_3.i => 0.0
        inductor_2.v => 0.0
        inductor_3.v => 0.0
       ]

prob = ODEProblem(sys, u0, (0, 0.01))
function condition(u,t,integrator)
    #-(integrator.sol[capacitor.v][end]+integrator.sol[resistor_1.v][end]+integrator.sol[inductor_1.v][end])<0.0
    integrator.sol[level.v][end]<0.0
end

#affect!(integrator) = terminate!(integrator)
function affect!(integrator)
    integrator.p[8] = 0.1*1e-3 
end 
cb = DiscreteCallback(condition,affect!, save_positions=(false,false))

sol = solve(prob, Rodas4(), reltol=1e-12,abstol=1e-12,callback=cb)

```

I did try again the switch and it really seems the problem is for defining a switch I have to change the equation depending in which state is the switch (open or close). I have tested the two equations separately and they work. The problem is to put them together in a component (I get an unbalanced problem because I have an extra equation).  
I will need somehow to be able to change the equation depending on the input: switch open (i~0), switch close (v~0). The parameter to change the state of the switch could be changed using a callback (this is my idea).  
Do you have an idea on how I could implement this?

---

<div class="post-metadata">

**Author:** ![ellocco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ellocco/32/31331_2.png) [@ellocco](https://discourse.julialang.org/u/ellocco)\
**Post date:** [November 4, 2022, 3:35pm UTC](https://discourse.julialang.org/t/modeling-a-crowbar-circuit-with-modelingtoolkit-jl/81968/6 "2022-11-04T15:35:10Z")

</div>

Hi, have you managed to solve your problem?  
I am still at the very beginning, the first thing was to test your script,  
unfortunately, it throws an error:

```julia
ERROR: LoadError: ArgumentError: SymbolicUtils.Symbolic{Real}[inductor_2₊iˍt(t)] are missing from the variable map.

```