# Unexpected(?) behavior using ModelingToolkit.jl, discrete\_events and irreducible

**URL:** <https://discourse.julialang.org/t/unexpected-behavior-using-modelingtoolkit-jl-discrete-events-and-irreducible/121142>\
**Category:** Modelling & Simulations\
**Created:** [October 10, 2024, 1:01pm UTC](https://discourse.julialang.org/t/unexpected-behavior-using-modelingtoolkit-jl-discrete-events-and-irreducible/121142 "2024-10-10T13:01:34Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![jmaedler](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmaedler/32/212218_2.png) [@jmaedler](https://discourse.julialang.org/u/jmaedler)\
**Post date:** [October 10, 2024, 1:01pm UTC](https://discourse.julialang.org/t/unexpected-behavior-using-modelingtoolkit-jl-discrete-events-and-irreducible/121142/1 "2024-10-10T13:01:34Z")

</div>

Hi everyone,

I am running into an issues using discrete events and irreducible variables in [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl). 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.

```julia
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](https://github.com/SciML/ModelingToolkit.jl/issues/1646), [#1797](https://github.com/SciML/ModelingToolkit.jl/issues/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 consistent initial conditions each time the callback is called. Is this intentional? What am I doing wrong here?  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/9/e975529366725809e78722ce6c8129fc1b926737.png)

**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](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/SampledData/)) the way to go?

Thank you for your help 🙂 !

---

<div class="post-metadata">

**Author:** ![aaliqureshi](https://avatars.discourse-cdn.com/v4/letter/a/c5a1d2/32.png) [@aaliqureshi](https://discourse.julialang.org/u/aaliqureshi)\
**Post date:** [October 10, 2024, 9:20pm UTC](https://discourse.julialang.org/t/unexpected-behavior-using-modelingtoolkit-jl-discrete-events-and-irreducible/121142/2 "2024-10-10T21:20:55Z")

</div>

> [@jmaedler](#):
>
> If I use the code above the volumes `V` and the levels `L` are set back to the consistent initial conditions each time the callback is called. Is this intentional? What am I doing wrong here?

The default behavior in DiffEq is to run a re-initialization scheme following a change in variables. To avoid that you can add modify you `affect!` function as:

```julia
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
    # force DiffEq to not re-initialize the variables
    u_modified!(integ,false)
end

```

 ![Screenshot 2024-10-10 at 5.19.21 PM](https://global.discourse-cdn.com/julialang/original/3X/2/d/2dc24add96c65486e96f990bbe46011f60715bd3.png)

I hope this is what you were looking for.

---

<div class="post-metadata">

**Author:** ![jmaedler](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmaedler/32/212218_2.png) [@jmaedler](https://discourse.julialang.org/u/jmaedler)\
**Post date:** [October 12, 2024, 1:40pm UTC](https://discourse.julialang.org/t/unexpected-behavior-using-modelingtoolkit-jl-discrete-events-and-irreducible/121142/3 "2024-10-12T13:40:24Z")

</div>

Yes, that solves my issue 🙂 Thank you very much!