How to represent Modelica style If-Equations using MTK

[JKRT]

I suppose this issue may very well be a documentation issue. Version of MTK is the latest (As of this writing 8.6.0) and the version of Julia was 1.7.2 Related PR I suppose is #1337 and one related issue is #1180

I am currently attempting to write a mockup for the following Modelica model (With the goal of generating MTK models for some of the MSL with high fidelity to the original Modelica implementation):

model IfEquationDer
  parameter Real u = 4;
  parameter Real uMax = 10;
  parameter Real uMin = 2;
  Real y;
equation
  if uMax < time then
    der(y) = uMax;
  elseif uMin < time then
    der(y) = uMin;
  else
    der(y) = u;
  end if;
end IfEquationDer;

The resulting simulation using OMEdit when simulating this system for 12 seconds is:

Using DifferentialEquations.jl I can represent this model in the following way: (Note that I do not necessarily need to use globals I suppose I can have them in the parameter array and update them in that way instead)

#=
 An attempt to model Modelica if-equations using DifferentialEquations.jl
=#
using DifferentialEquations
global ifCond1 = false
global ifCond2 = false
condition1(u, t, integrator) = integrator.p[1] - t
condition2(u, t, integrator) =  integrator.p[2] - t
function affect1!(integrator)
  @info "Affect 1 at:" integrator.t
  global ifCond1 = true
  global ifCond2 = false
end

function affect2!(integrator)
  @info "Affect 2 at:" integrator.t
  global ifCond1 = false
  global ifCond2 = true
end

cb1 = ContinuousCallback(
  condition1,
  affect1!,
  rootfind = true,
  save_positions = (true, true),
  affect_neg! = affect1!,
)

cb2 = ContinuousCallback(
  condition2,
  affect2!,
  rootfind = true,
  save_positions = (true, true),
  affect_neg! = affect2!,
)

function test(u, p, t)
  val = ifelse(ifCond1, p[1], ifelse(ifCond2, p[2], p[3]))
  @info val ifCond1 ifCond2
  [val]
end

tspan = (0.0,12.0)
prob = ODEProblem(test,
                  [0.0],
                  tspan,
                  [10.0, 2.0, 4.0])
global ifCond1 = false
global ifCond2 = false
sol = solve(prob,Tsit5(), callback = CallbackSet(cb1, cb2))
plot(sol)

Now to the issue: How could I in some way represent the same (contrived) system in MTK?

• Attempt 1.

#=
 An attempt to model Modelica if-equations
=#
#=
 An attempt to model Modelica if-equations
=#
using DifferentialEquations
using ModelingToolkit
@variables t
vars = @variables y(t)
pars = @parameters uMax uMin u
D = Differential(t)
global ifCond1 = false
global ifCond2 = false

function condition1(u, t, integrator)
  return integrator.p[1] - t
end

function condition2(u, t, integrator)
  return integrator.p[2] - t
end

function affect1!(u, t, integrator)
  global ifCond1 = true
  global ifCond2 = false
end

function affect2!(u, t, integrator)
  global ifCond1 = false
  global ifCond2 = true
end

cb1 = ContinuousCallback(
  condition1,
  affect1!,
  rootfind = true,
  save_positions = (true, true),
  affect_neg! = affect1!,
)

cb2 = ContinuousCallback(
  condition2,
  affect2!,
  rootfind = true,
  save_positions = (true, true),
  affect_neg! = affect1!,
)

eqs = [
  D(y) ~ ifelse(ifCond1, uMax, ifelse(ifCond2, uMin, u)),
]

@named ifequationDer = ODESystem(eqs, t, vars, pars)
ifequationDer = structural_simplify(ifequationDer)
tspan = (0.0,12.0)
prob = ODEProblem(ifequationDer,
                  tspan,
                  [y => 0.0],
                  [uMax => 10.0,
                   uMin => 2.0,
                   u => 4.0],
                  callbacks = CallbackSet(cb1, cb2))
sol = solve(prob,Tsit5())
plot(sol)

One issue with this code is that MTK automatically simplifies

julia> eqs = [
         D(y) ~ ifelse(ifCond1, uMax, ifelse(ifCond2, uMin, u)),
       ]
1-element Vector{Equation}:
 Differential(t)(y(t)) ~ u

even though ifCond1,ifCond2 are not constants (Maybe possible that some constant folding optimization is too greedy somewhere (point of discussion) ). Furthermore, I get the following error message:

julia> prob = ODEProblem(ifequationDer,
                         tspan,
                         [y => 0.0],
                         [uMax => 10.0,
                          uMin => 2.0,
                          u => 4.0],
                         callbacks = CallbackSet(cb1, cb2))
ERROR: ArgumentError: Equations (1), states (1), and initial conditions (2) are of different lengths. To allow a different number of equations than states use kwarg check_length=false.

I suppose the last one might be that I wrote the code above wrongly in some way as well. However, I do not get why it claims that I have two initial conditions?

However, the main functionality I do not grasp how to get is how to use ifelse and update the condition in an ifelse via some continuous callback. I can see how I can use the new affect/continuous condition functionality to change the value of some state variable, but how (If possible) can I use that to change some condition in an ifelse expression?

Sorry for a long post as always..

Cheers, John

[baggepinnen]

I think you might need to use IfElse.ifelse from the IfElese.jl package. The Core.ifelse function is not a generic function in julia <=v1.7 so it can not be extended to symbolic expressions.

[JKRT]

@baggepinnen

Hi Fredrik, thanks!

Yea, I should clarify which version I used. Latest MTK and Julia 1.7.2. That is Julia >= v1.7 but I also tried now using the if-else package to see what happens:

#=
 An attempt to model Modelica if-equations
=#
using DifferentialEquations
using ModelingToolkit
using IfElse

@variables t
vars = @variables y(t)
pars = @parameters uMax uMin u
D = Differential(t)
global ifCond1 = false
global ifCond2 = false

function condition1(u, t, integrator)
  return integrator.p[1] - t
end

function condition2(u, t, integrator)
  return integrator.p[2] - t
end

function affect1!(u, t, integrator)
  global ifCond1 = true
  global ifCond2 = false
end

function affect2!(u, t, integrator)
  global ifCond1 = false
  global ifCond2 = true
end

cb1 = ContinuousCallback(
  condition1,
  affect1!,
  rootfind = true,
  save_positions = (true, true),
  affect_neg! = affect1!,
)

cb2 = ContinuousCallback(
  condition2,
  affect2!,
  rootfind = true,
  save_positions = (true, true),
  affect_neg! = affect1!,
)

eqs = [
  D(y) ~ ifelse(ifCond1, uMax, ifelse(ifCond2, uMin, u)),
]

@named ifequationDer = ODESystem(eqs, t, vars, pars)
ifequationDer = structural_simplify(ifequationDer)
tspan = (0.0,12.0)
prob = ODEProblem(ifequationDer,
                  tspan,
                  [y => 0.0],
                  [uMax => 10.0,
                   uMin => 2.0,
                   u => 4.0],
                  callbacks = CallbackSet(cb1, cb2))
sol = solve(prob,Tsit5())
plot(sol)

Same results and it also seems to remove my ifelse when I try to create the problem

julia> eqs = [
         D(y) ~ ifelse(ifCond1, uMax, ifelse(ifCond2, uMin, u)),
       ]
1-element Vector{Equation}:
 Differential(t)(y(t)) ~ u

[baggepinnen]

You're not using IfElse.ifelse

[JKRT]

@baggepinnen correct 💯 , thanks!

Like so then:

#=
 An attempt to model Modelica if-equations
=#
using DifferentialEquations
using ModelingToolkit
using Plots
import IfElse

@variables t
vars = @variables y(t)
pars = @parameters uMax uMin u
D = Differential(t)
global ifCond1 = false
global ifCond2 = false

function condition1(u, t, integrator)
  return integrator.p[1] - t
end

function condition2(u, t, integrator)
  return integrator.p[2] - t
end

function affect1!(u, t, integrator)
  global ifCond1 = true
  global ifCond2 = false
end

function affect2!(u, t, integrator)
  global ifCond1 = false
  global ifCond2 = true
end

cb1 = ContinuousCallback(
  condition1,
  affect1!,
  rootfind = true,
  save_positions = (true, true),
  affect_neg! = affect1!,
)

cb2 = ContinuousCallback(
  condition2,
  affect2!,
  rootfind = true,
  save_positions = (true, true),
  affect_neg! = affect1!,
)

eqs = [
  D(y) ~ IfElse.ifelse(ifCond1, uMax, IfElse.ifelse(ifCond2, uMin, u)),
]

@named ifequationDer = ODESystem(eqs, t, vars, pars)
ifequationDer = structural_simplify(ifequationDer)
tspan = (0.0,12.0)
prob = ODEProblem(ifequationDer,
                  tspan,
                  [y => 0.0],
                  [uMax => 10.0,
                   uMin => 2.0,
                   u => 4.0],
                  callbacks = CallbackSet(cb1, cb2))
sol = solve(prob,Tsit5())
plot(sol)

Still the same issue

[baggepinnen]

ifCond1 these parameters must be symbolic, try something like

@parameters ifCond1=false

[JKRT]

Thanks,

Still no luck:(

#Same preamble 
@variables t
vars = @variables y(t)
pars = @parameters uMax uMin u
@parameters ifCond1 = false
@parameters ifCond2 = false
D = Differential(t)

function condition1(u, t, integrator)
  return integrator.p[1] - t
end

function condition2(u, t, integrator)
  return integrator.p[2] - t
end

function affect1!(u, t, integrator)
  integrator.p[4] = true #ifCond1 = true
  itnegrator.p[5] = false #ifCond2 = false
end

function affect2!(u, t, integrator)
  integrator.p[4] = false  #ifCond1 = false
  integrator.p[5] = true   #ifCond2 = true
end

cb1 = ContinuousCallback(
  condition1,
  affect1!,
  rootfind = true,
  save_positions = (true, true),
  affect_neg! = affect1!,
)

cb2 = ContinuousCallback(
  condition2,
  affect2!,
  rootfind = true,
  save_positions = (true, true),
  affect_neg! = affect1!,
)

eqs = [
  D(y) ~ IfElse.ifelse(ifCond1, uMax, IfElse.ifelse(ifCond2, uMin, u)),
]

@named ifequationDer = ODESystem(eqs, t, vars, pars)
ifequationDer = structural_simplify(ifequationDer)
tspan = (0.0,12.0)
prob = ODEProblem(ifequationDer,
                  tspan,
                  [y => 0.0],
                  [uMax => 10.0,
                   uMin => 2.0,
                   u => 4.0],
                  callbacks = CallbackSet(cb1, cb2))
sol = solve(prob,Tsit5())
plot(sol)

This results in:

julia> @time include("IfEqModelingToolkit.jl")
ERROR: LoadError: TypeError: non-boolean (Sym{Real, Base.ImmutableDict{DataType, Any}}) used in boolean context

If I do it like this, how can I access the variables in the callbacks, now I just guessed they get idx 4 and 5 similar to my DifferentialEquations.jl example? Or can I refer to ifCond1 and 2 directly?

[baggepinnen]

Hmm, it looks like IfElse.ifelse requires an expression

eqs = [
  D(y) ~ IfElse.ifelse(ifCond1==true, uMax, IfElse.ifelse(ifCond2==true, uMin, u)),
]

Your tspan and initial state are in the wrong order btw

[JKRT]

That works to create the system, but for some reason the problem errors with the same message as before

julia> @time include("IfEqModelingToolkit.jl")
ERROR: LoadError: ArgumentError: Equations (1), states (1), and initial conditions (2) are of different lengths. To allow a different number of equations than states use kwarg check_length=false.

[baggepinnen]

Your tspan and initial state are in the wrong order

[baggepinnen]

If I do it like this, how can I access the variables in the callbacks, now I just guessed they get idx 4 and 5 similar to my DifferentialEquations.jl example? Or can I refer to ifCond1 and 2 directly?

This is currently a limitation of using manually constructed callbacks with MTK. ref: https://github.com/SciML/ModelingToolkit.jl/pull/1438#issuecomment-1023087109

[baggepinnen]

Maybe you could make use of the symbolic callback functionality in MTK? https://mtk.sciml.ai/dev/basics/Composition/#Components-with-discontinuous-dynamics

The example that looks like this

root_eqs = [x ~ 0]  # the event happens at the ground x(t) = 0
affect   = [v ~ -v] # the effect is that the velocity changes sign

@named ball = ODESystem([
    D(x) ~ v
    D(v) ~ -9.8
], t; continuous_events = root_eqs => affect) # equation => affect

can perhaps be adapted to your situation?

[JKRT]

Thanks again @baggepinnen 😅 I started out with that, but then I got into the ifelsetrouble, so I switched to how I would try to represent it using DifferentialEquations.jl I suppose what I missed when I dabbled with it was the ifCond==true part (Also missing the initial conditions which I for some reason had removed..), big thanks.

Another question, while I understand that the symbolic callbacks work for continuous-time conditions I did not get how I could configure it for discrete conditions/clocked time, is that currently possible?

My attempt with the bouncing ball suggestion:

using DifferentialEquations
using ModelingToolkit
using Plots
import IfElse

@variables t
vars = @variables y(t)
pars = @parameters uMax uMin u
@parameters ifCond1 = false
@parameters ifCond2 = false
D = Differential(t)

root_eqs = [uMax - t ~ 0,
            uMin - t ~ 0]
affect   = [(ifCond1, ifCond2 ) ~ (true, false),
            (ifCond1, ifCond2 ) ~ (false, true)
            ]
eqs = [
  D(y) ~ IfElse.ifelse(ifCond1 == true, uMax, IfElse.ifelse(ifCond2 == true, uMin, u)),
]

@named ifequationDer = ODESystem(eqs, t, vars, pars; continuous_events = root_eqs => affect)
ifequationDer = structural_simplify(ifequationDer)
tspan = (0.0,12.0)
prob = ODEProblem(ifequationDer,
                  [y => 0.0],
                  tspan,
                  [uMax => 10.0,
                   uMin => 2.0,
                   u => 4.0])
sol = solve(prob,Tsit5())
plot(sol)

julia> prob = ODEProblem(ifequationDer,
                         [y => 0.0],
                         tspan,
                         [uMax => 10.0,
                          uMin => 2.0,
                          u => 4.0],
                         callbacks = CallbackSet(cb1, cb2))
ERROR: affected variables not unique, each state can only be affected by one equation for a single `root_eqs => affects` pair.

Of course, this is since I would like to do two assignments on one effect, we have already discussed this I think (and to represent Modelica if-equations I need to activate and deactivate different branches)

So the error message I get is:

ERROR: affected variables not unique, each state can only be affected by one equation for a single `root_eqs => affects` pair.
Stacktrace:

[baggepinnen]

The following runs, but please check that it does what you expect :)

I had to make the conditions @variables instead of parameters and treat them as dynamic states (with zero dynamics)

using DifferentialEquations
using ModelingToolkit
using Plots
import IfElse

@variables t
vars = @variables y(t)
pars = @parameters uMax uMin u
@variables ifCond1(t) = false
@variables ifCond2(t) = false
D = Differential(t)

continuous_events = [
    (uMax - t ~ 0) => [ifCond1 ~ true, ifCond2 ~ false]
    (uMin - t ~ 0) => [ifCond1 ~ false, ifCond2 ~ true]
]

eqs = [
  D(y) ~ IfElse.ifelse(ifCond1 == true, uMax, IfElse.ifelse(ifCond2 == true, uMin, u)),
  D(ifCond1) ~ 0,
  D(ifCond2) ~ 0,
]

@named ifequationDer = ODESystem(eqs, t, vars, pars; continuous_events)
ifequationDer = structural_simplify(ifequationDer)
tspan = (0.0,12.0)
prob = ODEProblem(ifequationDer,
                  [y => 0.0, ifCond1=>false, ifCond2=>false],
                  tspan,
                  [uMax => 10.0,
                   uMin => 2.0,
                   u => 4.0])
sol = solve(prob,Tsit5())

plot(sol)

[JKRT]

The following runs, but please check that it does what you expect :)

I had to make the conditions @variables instead of parameters and treat them as dynamic states (with zero dynamics)

using DifferentialEquations
using ModelingToolkit
using Plots
import IfElse

@variables t
vars = @variables y(t)
pars = @parameters uMax uMin u
@variables ifCond1(t) = false
@variables ifCond2(t) = false
D = Differential(t)

continuous_events = [
    (uMax - t ~ 0) => [ifCond1 ~ true, ifCond2 ~ false]
    (uMin - t ~ 0) => [ifCond1 ~ false, ifCond2 ~ true]
]

eqs = [
  D(y) ~ IfElse.ifelse(ifCond1 == true, uMax, IfElse.ifelse(ifCond2 == true, uMin, u)),
  D(ifCond1) ~ 0,
  D(ifCond2) ~ 0,
]

@named ifequationDer = ODESystem(eqs, t, vars, pars; continuous_events)
ifequationDer = structural_simplify(ifequationDer)
tspan = (0.0,12.0)
prob = ODEProblem(ifequationDer,
                  [y => 0.0, ifCond1=>false, ifCond2=>false],
                  tspan,
                  [uMax => 10.0,
                   uMin => 2.0,
                   u => 4.0])
sol = solve(prob,Tsit5())

plot(sol)

It looks correct, awesome! I suppose this issue can be closed, but maybe we should make a contribution to the documentation (I can attempt to write something up). Ideally, maybe fix the issues with having to introduce the zero dynamic states with new equations (Where would I start?). Still, using the techniques you proposed I should be able to emulate Modelica if-branches, so maybe this is more documentation issue

[baggepinnen]

Yeah, there are a few things that should be fixed here

• Accept symbol as first argument in IfElse.ifelse (currently requires Equation)
• Figure out why it didn't work with ifCond1 as parameter rather than @variable
• Make this functionality more discoverable in the docs. I myself find the structure of the docs slightly hard to navigate and am not sure what the intent is. A lot of functionality appears to be somewhat hidden under "tutorials" which can make it hard to find.

[JKRT]

Accept symbol as first argument in IfElse.ifelse (currently requires Equation)

I looked into this briefly, It seems that there might be an issue with the IfElse package, not sure Calling this the ordinary way it fails on line 37 in num.jl in Symbolics

IfElse.ifelse(x::Num,y,z) = begin
  Num(IfElse.ifelse(value(x), value(y), value(z)))
end

ERROR: LoadError: MethodError: no method matching ifelse(::Term{Real, Base.ImmutableDict{DataType, Any}}, ::Sym{Real, Base.ImmutableDict{DataType, Any}}, ::Term{Real, Nothing})
Closest candidates are:
  ifelse(::SymbolicUtils.Symbolic{Bool}, ::Any, ::Any) at C:\Users\John\.julia\packages\SymbolicUtils\v2ZkM\src\methods.jl:174
  ifelse(::Static.False, ::Any, ::Any) at C:\Users\John\.julia\packages\Static\pkxBE\src\bool.jl:127
  ifelse(::Num, ::Any, ::Any) at C:\Users\John\Projects\Programming\JuliaPackages\Symbolics.jl\src\num.jl:37
  ...
Stacktrace:
 [1] ifelse(x::Num, y::Num, z::Num)
   @ Symbolics C:\Users\John\Projects\Programming\JuliaPackages\Symbolics.jl\src\num.jl:40
 [2] top-level scope
   @ C:\Users\John\Projects\Programming\JuliaPackages\Scratch\ifEqTest.jl:20
 [3] include(fname::String)
   @ Base.MainInclude .\client.jl:451
 [4] top-level scope
   @ REPL[33]:1

This is since the Base variant of ifelse requires a Bool I suppose. I guess is it because the final PR restricted the type of the first parameter to a Boolean:

So the backport code in is no longer in synch with some of the logic of symbolic I think, it seems that the code is written s.t the first argument could be more generic.

I experimented a bit with IfElse.jl and forced it to call Core directly (That is removing the code in the PR above)

ERROR: LoadError: TypeError: non-boolean (Term{Real, Base.ImmutableDict{DataType, Any}}) used in boolean context
Stacktrace:
 [1] ifelse(::Term{Real, Base.ImmutableDict{DataType, Any}}, ::Sym{Real, Base.ImmutableDict{DataType, Any}}, ::Term{Real, Nothing})
   @ IfElse C:\Users\John\Projects\Programming\JuliaPackages\IfElse.jl\src\IfElse.jl:3
 [2] ifelse(x::Num, y::Num, z::Num)
   @ Symbolics C:\Users\John\Projects\Programming\JuliaPackages\Symbolics.jl\src\num.jl:40
 [3] top-level scope
   @ C:\Users\John\Projects\Programming\JuliaPackages\Scratch\ifEqTest.jl:20
 [4] include(fname::String)
   @ Base.MainInclude .\client.jl:451
 [5] top-level scope
   @ REPL[33]:1

However, that also does not seem to work (With reservation that this digging might be wrong I am not 100% familiar with the Symbolic packages) having a user being forced to write ifCond1 == true might be a bit verbose, but on the other hand it can probably be simplified on the user side by a macro

[JKRT]

Figure out why it didn't work with ifCond1 as a parameter rather than @variable

I think this was something we both encountered as a "false positive". I tried to change your second example to parameters. However, I am not sure that it should, because ifCond1 is not a time derivative. if I remove the dummy equations I encounter the same issue:

ERROR: LoadError: UndefVarError: ifCond1 not defined
Stacktrace:
  [1] macro expansion
    @ C:\Users\John\.julia\packages\SymbolicUtils\v2ZkM\src\code.jl:394 [inlined]
  [2] macro expansion
    @ C:\Users\John\Projects\Programming\JuliaPackages\Symbolics.jl\src\build_function.jl:452 [inlined]
  [3] macro expansion
    @ C:\Users\John\.julia\packages\SymbolicUtils\v2ZkM\src\code.jl:351 [inlined]
  [4] macro expansion
    @ C:\Users\John\.julia\packages\RuntimeGeneratedFunctions\KrkGo\src\RuntimeGeneratedFunctions.jl:129 [inlined]
  [5] macro expansion
    @ .\none:0 [inlined]
  [6] generated_callfunc
    @ .\none:0 [inlined]

Probably having equations defines the parameters in the correct scope s.t the affect equation is able to make sense of the symbol.

Code for replication:

using DifferentialEquations
using ModelingToolkit
using Plots
import IfElse

@variables t
vars = @variables y(t)
pars = @parameters uMax uMin u
@parameters ifCond1(t) = false
@parameters ifCond2(t) = false
D = Differential(t)

continuous_events = [
    (uMax - t ~ 0) => [ifCond1 ~ true, ifCond2 ~ false]
    (uMin - t ~ 0) => [ifCond1 ~ false, ifCond2 ~ true]
]

eqs = [
  D(y) ~ IfElse.ifelse(ifCond1 == true, uMax, IfElse.ifelse(ifCond2 == true, uMin, u)),
  D(ifCond1) ~ 0,
  D(ifCond2) ~ 0,
]

@named ifequationDer = ODESystem(eqs, t, vars, pars; continuous_events)
ifequationDer = structural_simplify(ifequationDer)
tspan = (0.0,12.0)
prob = ODEProblem(ifequationDer,
                  [y => 0.0, ifCond1=>false, ifCond2=>false],
                  tspan,
                  [uMax => 10.0,
                   uMin => 2.0,
                   u => 4.0])
sol = solve(prob,Tsit5())

plot(sol)

Cheers, John

[ChrisRackauckas]

No longer needs IfElse.jl, just core's ifelse works fine.