`@discrete_events` not working with acausal modeling framework (w/ reproducible example)

[joseph03choi]

Here is the full code to reproduce the error with @discrete_events not being able to see connector variables.

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

@connector Stream begin
    q(t)
    c_A(t)
    c_B(t)
    c_C(t)
end

@mtkmodel Influent begin
    @components begin
        out = Stream()
    end
    @parameters begin
        qin = 4.0
        c_Ain = 1.15
        c_Bin = 0.3
        c_Cin = 0.0
    end
    @equations begin
        out.q ~ qin
        out.c_A ~ c_Ain
        out.c_B ~ c_Bin
        out.c_C ~ c_Cin
    end
end

@mtkmodel Effluent begin
    @components begin
        in = Stream()
    end
    @variables begin
        q(t)
        c_A(t)
        c_B(t)
        c_C(t)
    end
    @equations begin
        q ~ in.q
        c_A ~ in.c_A
        c_B ~ in.c_B
        c_C ~ in.c_C
    end
end

@mtkmodel CSTR begin
    @components begin
        in1 = Stream()
        in2 = Stream()
        out = Stream()
        readingsOut = Stream()
    end
    @variables begin
        q(t)
        c_A(t) = 0
        c_B(t) = 0
        c_C(t) = 0
    end
    @parameters begin
        k = 0.01
        V = 1000.0
    end
    begin
        k_La = 500/6*atan(0.02*in2.q)
    end
    @equations begin
        D(c_A) ~ in1.q/V*(in1.c_A - c_A) - k*c_A*c_B
        D(c_B) ~ in1.q/V*(in1.c_B - c_B) - k*c_A*c_B + k_La/V*(in2.c_B - c_B)
        D(c_C) ~ in1.q/V*(in1.c_C - c_C) + k*c_A*c_B
        out.q ~ in1.q
        out.c_A ~ c_A
        out.c_B ~ c_B
        out.c_C ~ c_C 
        readingsOut.c_B ~ c_B
    end
end

@mtkmodel PI_Control_Valve begin
    @components begin
        readingsIn = Stream()
        out = Stream()
    end
    @variables begin
        I_err(t) = 0
        u(t) = 0

        # B(t) = 0 # Trying intermediate variable, doesn't work
    end
    @parameters begin
        c_BSP = 0.32
        K_b = 0.0
        K_c = 9.0
        K_I = 0.18
        eventStart = 200
        eventInterval = 30
        s = 0 # No initial signal to open valve

        τ_p = 5.25
        c_Bmax = 0.61
    end
    @equations begin
        D(I_err) ~ c_BSP - readingsIn.c_B
        D(u) ~ (s - u)/τ_p
        out.q ~ 10*u
        out.c_B ~ c_Bmax

        # B ~ readingsIn.c_B # Discrete events still can't see this variable
        # D(B) ~ 0 # Discrete events can only see B if I do this, but then I can't do `B ~ readingsIn.c_B`
    end
    @continuous_events begin
        [u ~ 1] => [s ~ u]
    end
    @discrete_events begin
        (u > 1-1e-13) => [u ~ 1-1e-13]

        # Discrete events can't see readingsIn.c_B or B
        # If you replace readingsIn.c_B here with a number (1 for example), the code runs with no errors
        ((mod(t-eventStart,eventInterval) == 0) & (t >= eventStart)) => [s ~ K_b + K_c*(c_BSP - readingsIn.c_B) + K_I*I_err]
    end
end

@mtkmodel CSTR_Model begin
    @components begin
        influent = Influent()
        effluent = Effluent()
        cstr = CSTR()
        pi_control_valve = PI_Control_Valve()
    end
    @equations begin
        connect(influent.out, cstr.in1)
        connect(cstr.out, effluent.in)
        connect(cstr.readingsOut, pi_control_valve.readingsIn)
        connect(pi_control_valve.out, cstr.in2)
    end
end

@mtkbuild cstr_Model = CSTR_Model()

prob = ODEProblem(cstr_Model, [], (0,2000), [])

# ERROR: UndefVarError: `pi_control_valve₊readingsIn₊c_B` not defined
sol = solve(prob, Tsit5(), tstops = [200:30:2000;])

plot(
    sol.t, vcat(sol[1:3,:],sol[5,:]')',
    title = "Concentration & Valve Pos. vs Time",
    yaxis = ("Concentration (mol/L)",(0,1.01)),
    xaxis = ("Time (min)",(0,2000*1.01)),
    label = ["c_A" "c_B" "c_C" "u"]
)

Everything runs properly until you get to solve(), where you get the error message ERROR: UndefVarError: 'pi_control_valve₊readingsIn₊c_B' not defined. You can replace the readingsIn.c_B under @discrete_events with a number and the problem can be solved, so we know this is the issue.

[contradict]

Try using your intermediate variable and mark it irreducible B(t) = 0, [irreducible=true] to prevent it from being eliminated during simplification.

[ChrisRackauckas]

yeah right now it's a bit manual to make sure things aren't eliminated, we need to automate that.