Modeling a crowbar circuit with ModelingToolkit.jl

Specific Domains Modelling & Simulations
1

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:

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)
2

Hello and welcome to the community!

I’m struggling to grasp exactly what you mean by

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?

3

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

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.

4

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.

5

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

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?

6

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:

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