zekeriyasari / Causal.jl

Causal.jl - A modeling and simulation framework adopting causal modeling approach.
https://zekeriyasari.github.io/Causal.jl/dev/
Other
115 stars 7 forks source link

usage of callbacks on computed value #72

Closed YayeIrene closed 3 years ago

YayeIrene commented 3 years ago

Hello, I am trying to simulate a power source which consists in a series capacitors banks which are discharged successively. Each capacitor bank is an RLC circuit with a switch. When the tension on the switch reaches 0 the capacitor is taken out of the circuit. My idea is to simulate with causal.jl. I have built the RLC circuit as an ode_system using @def_ode_system. This seems to work. (My next step will be to try to build it component by component!) I then created a model which has an RLC circuit node and other nodes necessary to compute the tension. I created also a callback, the condition is that the tension is zero. The tension is an output of an Adder component. I am having trouble accessing this value during simulation. using take!(writer) the simulation get stuck at the run stage! Is there a way to read the output of the Adder during simulation to be used by a callback condition function?

Thank you! Irene

zekeriyasari commented 3 years ago

Hi @YayeIrene

Since you did not provide me with your code, I could not understand the details of the system you are trying to simulate.

using take!(writer) the simulation get stuck at the run stage!

You may rarely use take!(writer) for your simulations! take! function is one of the core internal functions to make the data flow through the connections of the model. I am sure there must be something else to solve the problem, but I need to see your use case to elaborate!

Hoping that it helps, here is a minimum working example that simulates an RLC circuit driven by step source.

The block diagram model of a series RLC circuit is as below. Untitled Diagram

Here ds_rl is a dynamical system consisting of the series combination of the components resistor R and inductor L. ds_c is the dynamical system corresponding to capacitor C. The equations modeling the ds_rl` are

\dot{x}(t) = -R / L *  x(t) + 1 / L *  u(t)
y(t) = x(t)

where x(t), u(t), y(t) are the state, input and output, respectively. The equations modelling ds_c are

\dot{x}(t) = 1 / C * u(t)
y(t) = x(t) 

So, here is the example code for the simulation of an RLC circuit (with component values R=1 Ohm, L=1 Henry, C=1 Farads).

using Causal 
using Plots 

@defmodel model begin 
    @nodes begin 
        gen = StepGenerator()
        ds_rl = ODESystem(
            righthandside = (dx, x, u, t) -> (dx[1] = -x[1] + u[1](t)), 
            readout = (x, u, t) -> x, 
            state = rand(1), 
            input = Inport(1), 
            output = Outport(1)
        )
        ds_c = ODESystem(
            righthandside = (dx, x, u, t) -> (dx[1] = u[1](t)), 
            readout = (x, u, t) -> x, 
            state = rand(1), 
            input = Inport(1), 
            output = Outport(1)
        )
        adder = Adder(
            signs = (+, -),
            input = Inport(2), 
            output = Outport(1)
        )
        writer = Writer(
            input = Inport(4) 
        ) 
    end 
    @branches begin 
        gen[1] => adder[1] 
        adder[1] => ds_rl[1] 
        ds_rl[1] => ds_c[1] 
        ds_c[1] => adder[2] 
        gen[1] => writer[1] 
        ds_rl[1] => writer[2] 
        ds_c[1] => writer[3] 
        adder[1] => writer[4] 
    end 
end

simulate!(model, 0., 0.01, 50.)

t, x = read(getnode(model, :writer).component)

plt = plot(layout=(2,2))
plot!(t, x[:, 1], label="u(t)", subplot=1, ylims=(0, 1.5))
plot!(t, x[:, 2], label="i(t)", subplot=2)
plot!(t, x[:, 3], label="v0(t)", subplot=3)
plot!(t, x[:, 4], label="vi(t) - v0(t)", subplot=4)

Here are the simulation results.

rlc

YayeIrene commented 3 years ago

Hello,

Thank you for answer! I am trying to simulate a circuit with a crowbar (D) and a switch gap spark (E):

image

As long as the voltage across the crowbar is positive the current flows through the switch. When the voltage change sign the capacitor is short-circuited. The issue I am having is to detect the moment of the switch and act on it. Here is my code (I think I will use something similar to yours, it will then be easy to drop the capacitor). I didn't manage to use callbacks with the macro (@defmodel):

`using Causal using Plots

function statefunc(dx, x, u, t)

dx[1] = x[2]
dx[2] = -1/(L1+L2+L3)*((R1+R2+R3)*x[2]+1/C*x[1])

end

function outputfunc(x, u, t) return x[1] end

@def_ode_system mutable struct CapacitorBank{RH, RO, IP, OP} <: AbstractODESystem righthandside::RH = statefunc readout::RO = outputfunc state::Vector{Float64} = [0.,Uc/L] input::IP = nothing #Inport(1) output::OP = Outport(1) end

L1 = 0.61e-6; L2 = 41e-6 L3 = 1.11e-6; R1= 41e-3; R2 = 0.281e-3; R3 = 3.61e-3

C = 3.081e-3 ; Uc = 101e3

L = L1+L2+L3

R = R1+R2+R3

function condition(railgun) take!(railgun.nodes[10].component.input)<0.0 end

function action(railgun) terminate!(railgun) end

cb = Callback(condition=condition, action=action)

ti, dt, tf = 0, 0.00001, 0.01

railgun = Model(callbacks=cb)

addnode!(railgun, CapacitorBank(), label=:ds1) addnode!(railgun, Writer(), label=:writer) addnode!(railgun, Gain(input=Inport(1), gain=R1), label=:res) addnode!(railgun, Integrator(state=zeros(1), t=0., modelargs=(), solverargs=()), label=:inti) addnode!(railgun, ConstantGenerator(;amplitude=C), label=:cond) addnode!(railgun, Multiplier(ops=(, )), label=:mlt) addnode!(railgun, ConstantGenerator(;amplitude=L1), label=:ind) addnode!(railgun, Multiplier(ops=(*, /)), label=:mlt2) addnode!(railgun, Adder(signs=(+,+,-)), label=:adder) addnode!(railgun, Writer(), label=:writer2)

addbranch!(railgun, :ds1 => :writer, 1 => 1) addbranch!(railgun, :ds1 => :res, 1 => 1) addbranch!(railgun, :ds1 => :inti, 1 => 1) addbranch!(railgun, :ds1 => :mlt, 1 => 1) addbranch!(railgun, :ind => :mlt, 1 => 2) addbranch!(railgun, :inti => :mlt2, 1 => 1) addbranch!(railgun, :cond => :mlt2, 1 => 2) addbranch!(railgun, :res => :adder, 1 => 1) addbranch!(railgun, :mlt => :adder, 1 => 2) addbranch!(railgun, :mlt2 => :adder, 1 => 3) addbranch!(railgun, :adder => :writer2, 1 => 1)

simulate!(railgun, ti, dt, tf, withbar=false)

t, x = read(getnode(railgun, :writer).component) t, v = read(getnode(railgun, :writer2).component)

plot(t, x, xlabel="t", ylabel="x", label="") plot(t, v, xlabel="t", ylabel="v", label="") `

zekeriyasari commented 3 years ago

Hi @YayeIrene,

Sorry, but I could not follow the circuit schematic from your script. For example, there are no L1, L2, L3 in the schematics. And I could not get how you could represent capacitor bank dynamics by

dx[1] = x[2]
dx[2] = -1/(L1+L2+L3)*((R1+R2+R3)*x[2]+1/C*x[1])

Besides, you should not call take! manually in your callback. Again, you should not call terminate in the callback neither. Those are internal functions used by Causal.jl internally while simulating the model.

But I understand that you are trying to simulate a switched system. In case it helps, here is the simulation of a simple buck converter. The circuit schematics is given below. buck_circuit

The circuit has two modes: when the switch is closed and opened. The dynamics corresponding to these modes are state_space

So, the right-hand side of the differential equation system is piecewise-defined corresponding to the two modes of the circuit. So, we construct an ODESystem that models the buck converter circuit. ODESystem has two inputs: a PWM signal and a constant voltage source. The value of the PWM signal determines if the switch is closed or opened. Thus, we have a piecewise-defined right-hand-side function. The block diagram and the corresponding script are given below.

block_diagram

using Causal 
using Plots 

# Circuit parameters
Vin = 50. 
R = 20 
L = 400e-6 
C = 100e-6
fs = 20e3
ts = 1 / fs 

# System matrices 
A = [0 -1 / L; 1 / C -1 / (R * C)]
B_switch_closed = [1 / L, 0] 
B_switch_openned = [0., 0.] 
righthandside_switch_closed(dx, x, u, t)  = ( dx .= A * x + B_switch_closed * u(t) ) 
righthandside_switch_opened(dx, x, u, t) = ( dx .= A * x + B_switch_openned * u(t) ) 
function righthandside(dx, x, u, t)
    switch, input = u[1](t), u[2] 
    if switch ≥ 0.5 
        righthandside_switch_closed(dx, x, input, t)
    else 
        righthandside_switch_opened(dx, x, input, t)
    end
end 
readout(x, u, t) = x 

# Construct model 
@defmodel model begin 
    @nodes begin 
        pwm = SquarewaveGenerator(period = ts, duty=0.4)
        gen = ConstantGenerator(amplitude = Vin)
        buck = ODESystem(
            righthandside = righthandside, 
            readout = readout, 
            state = [0., 20.], 
            input = Inport(2), 
            output = Outport(2)
        )
        writer = Writer(input = Inport(4))
    end 
    @branches begin 
        pwm[1] => buck[1] 
        gen[1] => buck[2] 
        pwm[1] => writer[1]
        gen[1] => writer[2]
        buck[1:2] => writer[3:4]
    end 
end 

# Simulate model 
ti, dt, tf = 0., 0.01 * ts, 0.1
simulate!(model, ti, dt, tf)

# Read simulation data 
t, x = getcomponent(model, :writer) |> read 

# Plot solutions 
n = length(t)
nr = n - 1000 : n  # Last 1000 samples 
plot(t[nr], x[nr, 1], label="pwm") 
plot!(t[nr], x[nr, 2], label="input") 
plot!(t[nr], x[nr, 3], label="iL") 
plot!(t[nr], x[nr, 4], label="vC") 

buck

Let me reiterate that this is not the circuit that you want to simulate. Here, I gave this buck converter example as it is at least a switched system and it may help you to simulate your circuit. Hope it helps.

I think you do not need Callbacks in your model. You can just insert an if -else -end block in the right-hand-side of the dynamical system as you can see in the representation of the buck converter.

function righthandside(dx, x, u, t)
    switch, input = u[1](t), u[2] 
    if switch ≥ 0.5 
        righthandside_switch_closed(dx, x, input, t)
    else 
        righthandside_switch_opened(dx, x, input, t)
    end
end