Events seemingly ignored in system model

Specific Domains Modelling & Simulations
,
1

Hi guys,

I am running into an issue simulating a system model. Fundamentally it is two valves in sequence between two pressure sources. This is a test for the valve models which should allow these scenarios for a virtual commissioning simulation.

Now my issue, if I run the code below the events

    @discrete_events begin
        [30] => [binary_valve_1.S ~ 0.0]
        [60] => [binary_valve_1.S ~ 1.0, binary_valve_2.S ~ 0.0]
        [120] => [binary_valve_1.S ~ 0.0]
    end

seem to be ignored entirely and I have no idea why?!? If I remove one valve from the model, everything works fine. You can see that in the plot.

image
image733×492 4.99 KB

Below you can find the “MRE”. It is quiet large since as mentioned the error does just arise if I use both valves.

using ModelingToolkit
using ModelingToolkit: D_nounits as D, t_nounits as t
using DifferentialEquations
using Plots

@connector LiquidPort begin
    p(t)::Float64, [    description = "Set pressure in bar", 
                        guess = 1.01325]
    Vdot(t)::Float64, [ description = "Volume flow rate in L/min", 
                        guess = 0.0, 
                        connect = Flow]
end

@mtkmodel PressureSource begin
    @components begin
        port = LiquidPort()
    end
    @parameters begin
        p_set::Float64 = 1.01325, [description = "Set pressure in bar"]
    end
    @equations begin
        port.p ~ p_set
    end
end

@mtkmodel BinaryValve begin
    @constants begin
        p_ref::Float64 = 1.0, [description = "Reference pressure drop in bar"]
        ρ_ref::Float64 = 1000.0, [description = "Reference density in kg/m^3"]
    end
    @components begin
        port_in = LiquidPort()
        port_out = LiquidPort()
    end
    @parameters begin 
        k_V::Float64 = 1.0, [description = "Valve coefficient in L/min/bar"]
        k_leakage::Float64 = 1e-08, [description = "Leakage coefficient in L/min/bar"]
        ρ::Float64 = 1000.0, [description = "Density in kg/m^3"]
    end
    @variables begin
        S(t)::Float64, [description = "Valve state", guess = 0.0, irreducible = true]
        Δp(t)::Float64, [description = "Pressure difference in bar", guess = 1.0]
        Vdot(t)::Float64, [description = "Volume flow rate in L/min", guess = 1.0]
    end
    @equations begin
        # Port handling
        port_in.Vdot ~ -Vdot
        port_out.Vdot ~ Vdot
        Δp ~ port_in.p - port_out.p
        # System behavior
        D(S) ~ 0.0
        Vdot ~ S*k_V*sign(Δp)*sqrt(abs(Δp)/p_ref * ρ_ref/ρ) + k_leakage*Δp # softplus alpha function to avoid negative values under the sqrt
    end
end

# Test System
@mtkmodel TestSystem begin
    @components begin
        pressure_source_1 = PressureSource(p_set = 2.0)
        binary_valve_1 = BinaryValve(S = 1.0)
        binary_valve_2 = BinaryValve(S = 1.0)
        pressure_source_2 = PressureSource(p_set = 1.0)
    end
    @equations begin
        connect(pressure_source_1.port, binary_valve_1.port_in)
        connect(binary_valve_1.port_out, binary_valve_2.port_in)
        connect(binary_valve_2.port_out, pressure_source_2.port)
    end
    @discrete_events begin
        [30] => [binary_valve_1.S ~ 0.0]
        [60] => [binary_valve_1.S ~ 1.0, binary_valve_2.S ~ 0.0]
        [120] => [binary_valve_1.S ~ 0.0]
    end
end

# Test Simulation
@mtkbuild sys = TestSystem()

# Test Simulation
prob = ODEProblem(sys, [], (0.0, 150.0))
sol = solve(prob)

# Plot
plot(sol; idxs=[sys.binary_valve_1.S,sys.binary_valve_2.S])
2

Open an issue.