Unexpected(?) behavior using ModelingToolkit.jl, discrete_events and irreducible

[jmaedler]

Question❓

Hi everyone,

I am running into an issues using discrete events and irreducible variables in ModelingToolkit.jl. I made a post and Discourse and had brief Slack chat with @YingboMa. He asked me to ping @BenChung.

I hope the following example is sufficent as an MWE. It is a composite model of two tanks in sequence which are connected with valves. In the future, I want to set the valve position S from a soft-controller in a discrete event. For now, I just use a sine function. Furthermore, a level sensor should send back the level in the tank L to the soft-controller, which is "simulated" by println() within the Callback.

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

# Create a simple Tank model with reading and writing events
function affect!(integ,vars,pars,_)
    integ.ps[pars.S] = 0.5 + 0.5*sin(2*3.14*integ.ps[pars.freq]*integ.t) 
    println("$(integ.u[vars.L])") # Read L using the passed identifier
end

@mtkmodel Tank begin
    @constants begin
        g::Float64 = 9.81, [description = "Gravitational acceleration."]
    end
    @parameters begin
        S(t)::Float64 = 0.0
        k_V::Float64 = 10.0
        rho::Float64 = 1000.0
        A::Float64 = 1.0
        freq::Float64 = 0.1
    end
    @variables begin
        F_in(t)::Float64, [guess = 0.0]
        F_out(t)::Float64, [guess = 0.0]
        L(t)::Float64, [guess = 0.0, irreducible = true] # Making L irreducible
        V(t)::Float64 = 1.0
        Dp(t)::Float64, [guess = 0.0]
    end
    @equations begin
        D(V) ~ F_in - F_out
        V ~ A*L
        Dp ~ rho*g*L
        F_out ~ S*k_V*sign(Dp)*sqrt(abs(Dp)/100000*1000/rho)
    end
    @discrete_events begin
        1.0 => (affect!,[L],[freq,S],[S],nothing) # Pass indentfier for L to affect! function
    end
end

# Create composite model
@mtkmodel SimulateTankSystem begin
    @components begin
        tank1 = Tank(;A=1.0)
        tank2 = Tank(;A=2.0,freq=0.2)
    end
    @equations begin
        tank1.F_in ~ 2.0
        tank2.F_in ~ tank1.F_out
    end
end

# Main function
function main()
    @mtkbuild sys = SimulateTankSystem()

    prob = ODEProblem(sys, [], (0.0, 10.0);tstops=0:0.1:10.0)

    sol = solve(prob,QNDF()) # A DAE solver is needed here, as the equation for L remains the same!

    # Plotting
    plot(sol;idxs=sys.tank1.S)
    plot!(sol;idxs=sys.tank2.S)
    plot(sol;idxs=sys.tank1.L)
    plot!(sol;idxs=sys.tank2.L)
end

main()

Now the issue is, that structural_simplify() in mtkbuild() does eliminate L, if you do not apply a workarround. I studied multiple issues on this topic, e.g. #1646, #1797. They are still open but have been arround for quite a while already. They suggest to use irreducible = true to avoid that the required variables are eliminated from the system.

Now my problem: If I use the code above the volumes V and the levels L are set back to the initial each time the callback is called. I do not change anything about the variables... so I did not expect this behavior. Is this intentional? What am I doing wrong here?

Furthermore: Have other solutions to this problem perhaps been worked out in the meantime?

• I saw one suggest like reconstructing observables without using irreducible. I could hardcode this in the callback using integ.sol[integ.f.sys.tank1.L][end], but I am not sure what this does to performance of the simulation for a growing number of timesteps saved within sol. I have the impression that the whole time history is calculated each time I call integ.sol. Furthermore, I would like to construct more complex systems in future and am not sure yet, how to get the path to the specific tank in a generic form into the callback.
• Or are discrete-time variables (SampledData) the way to go?

Thank you for your help :-) !

[jmaedler]

In the discussion with @YingboMa we already figured out that the initial condition before starting the simulation is not consistent... even if the system of equations for the initialization is actually only linear I am a bit confused...

[BenChung]

Now my problem: If I use the code above the volumes V and the levels L are set back to the initial each time the callback is called. I do not change anything about the variables... so I did not expect this behavior. Is this intentional? What am I doing wrong here?

This, unfortunately, is due to a known issue with the initialization system that's ran after the callback happens. It reruns the initialization system that it runs at the beginning of the simulation, which then usually sets everything to 0s. There's a change in PR right now that fixes this behavior (it's come up a few times, so I'll probably cherry-pick that out into its own PR and then ask for that to be merged in advance of the rest of the stuff that's holding that big PR up). Right now, the only workaround to this issue is to disable the initialization system.

I'll try to get that fix out this weekend.

[jmaedler]

Hi @BenChung, thank you for the explanation. ...And thank you for pushing for a fix :-)

Is there also development going on to get rid of the irreducible = true workarround by providing a way to use observed variables in callbacks? I am just asking because this would be the ultimate fix to my issue.

[BenChung]

Yes, though it's probably going to be at least a couple of months, unfortunately.

Ultimately, we need to generate custom InitializationSystems for each callback. Once we do this, it's easy to let you use observed values in callbacks; we just stuff the new value into the initialization system and use the nonlinear solve to restore consistency. At the moment, though, our main priority is improving our handling of modular discrete systems (e.g. clocked behavior) and then proper use of InitializationSystems to restore consistency with observables will be part of what's needed to make that work correctly.

[jmaedler]

Thanks for the explanation! That helps me a lot!

Please allow me one more question: Your explanation sounds to me as if you need the individual initialization systems if you want to write the observed variables in the callback. Would a generic read access be easier to implement as an intermediate step?

Here is why I ask: I got the impression that I am not alone with my Hardware-in-the-Loop (HIL) / Software-in-the-Loop (SIL) use case for ModelingToolkit. Time-dependent parameters offer an excellent solution for writing to actuators. What is missing is a simple, generic solution for reading sensor values, which after simplification are either in the unkown or the observed variables. Such sensor values only need to be passed on to the programmable logic controller (PLC) using shared memory or similar, but do not need to be able to be written back to the simulation. This is done indirectly via the time-dependent parameters after processing by the PLC.

P.S.: Would it be interesting for you guys if I prepared my example as a tutorial for MTK and prepared a corresponding PR?