Hi everyone,
I am running into an issues using discrete events and irreducible variables in 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.
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 consistent initial conditions each time the callback is called. 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 usinginteg.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 withinsol. I have the impression that the whole time history is calculated each time I callinteg.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 
!
