ModelingToolkit.jl and differential equations with step/pulse inputs

Specific Domains Modelling & Simulations
1

In lots of differential equations applications, we want to solve a ‘forced’ DE whose evolution depends upon some input u(t). For instance we can consider a harmonic oscillator with a forcing term. The equation below makes a function create(), that takes an input u(t), and outputs an ODEProblem of the form

\dot{x}(t) = f(x,p, u(t)),
representing a harmonic oscillator being forced by your choice of input u(t), with parameters p.

function create(input)
    
  @parameters t
  @parameters k,c,m
  @derivatives D'~t
  @variables pos(t) vel(t) inp(t)

  eqs = [D(pos) ~ vel,
        D(vel) ~ input(t) - (1/m)*(c*vel + k*pos),
  ]
  ps = [k,c,m] .=>  [1.,4.,1.] 
  ics = [pos, vel] .=> [1.,0.]
  od = ODESystem(eqs, t, first.(ics),first.(ps))
  tspan = (0.,100.)
  prob = ODEProblem(od, ics, tspan, ps)
  return prob
end

If I take ‘nice’ inputs, everything works great. For instance,
prob = create(sin)

However, in a lot of important applications (eg pharmacokinetic modelling), inputs come as nondifferentiable functions, such as pulses or steps. E.g.

nasty(t) = t>1 ? 1 : 0

Now I can understand that functions with if statements and the like, are not convertable to symbolic code by ModelingToolkit.jl. Indeed, the code:

@parameters t
nasty(t)

gives an error:

ERROR: TypeError: non-boolean (Operation) used in boolean context

But it would be really nice if ModelingToolkit.jl could generate one of its’ awesome, compiled vector fields from eqs using inputs like nasty(t), without ever forming an (impossible) symbolic representation of the nasty bit. Basically I want some way of doing

prob = create(nasty)

Since a function like nasty(t) depends only on t, we don’t need to do anything symbolic (jacobians etc) with it when converting the overall system of equations eqs into the vector field of an ODEProblem.

Is there any way of doing something like this using ModelingToolkit.jl? Is it even reasonable?

Many thanks, and also a big thumbs up to everyone working on ModelingToolkit.jl. Although I’m not enough of a programmer to directly contribute, I’m really enjoying seeing it develop: it’s very cool. Hopefully I can contribute a few implementations of ODE models from the systems biology literature in the fullness of time.

2 Likes
2

A quick-and-dirty solution to the problem might be to overload the ifelse function for Operations, or create an @ifelse macro which does a similar thing. For now, you could use min and max to “simulate” an if statement (or a ternary operator)…

1 Like
3

Thanks for the reply!

So replacing my ternary operator with a minmax-style (equivalent) operator works! Pulses would get somewhat annoying with this type of syntax I imagine.

I guess I’m not sure what you mean by overloading the ifelse function. Is there a page of documentation / a file in the source code that implements something similar? Cheers!

4

Ah, on closer inspection it looks like someone tried that:

Unfortunately, ifelse is not in Base, but Core. This means that it’s a builtin function and cannot be extended, so that idea is moot. However, it is possible to do:

julia> @inline function if_else(condition::Bool, @nospecialize(trueval), @nospecialize(falseval))
           return ifelse(condition, trueval, falseval)
       end
if_else (generic function with 1 method)

julia> ModelingToolkit.@register if_else(x, trueval, falseval)
if_else (generic function with 8 methods)

julia> if_else |> methods
# 8 methods for generic function "if_else":
[1] if_else(condition::Bool, trueval, falseval) in Main at REPL[13]:2
[2] if_else(x::Expression, trueval::Expression, falseval::Expression) in Main
[3] if_else(x::Number, trueval::Expression, falseval::Expression) in Main
[4] if_else(x::Expression, trueval::Number, falseval::Expression) in Main
[5] if_else(x::Number, trueval::Number, falseval::Expression) in Main
[6] if_else(x::Expression, trueval::Expression, falseval::Number) in Main
[7] if_else(x::Number, trueval::Expression, falseval::Number) in Main
[8] if_else(x::Expression, trueval::Number, falseval::Number) in Main

julia> nasty(i) = if_else(i > 1, 1, 0)
nasty (generic function with 1 method)

julia> @variables t
(t,)

julia> nasty(t)
if_else(t > 1, 1, 0)

This should work in a differential equation.

Indeed,

julia> prob = create(nasty)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 100.0)
u0: [1.0, 0.0]

cc @ChrisRackauckas

Edit: Oh, I forgot about this…functions are eval’ed in ModelingToolkit’s scope, so you would get

var"##MTIIPVar#258"[2] = (ModelingToolkit).:-((ModelingToolkit).if_else((ModelingToolkit).:>(t, 1), 1, 0), (ModelingToolkit).:*((ModelingToolkit).:/(1, m), (ModelingToolkit).:+((ModelingToolkit).:*(c, vel), (ModelingToolkit).:*(k, pos))))

this IR. Note the ModelingToolkit.if_else - you may have to @eval ModelingToolkit using Main: if_else.

1 Like
5

Not the most general solution to your problem, but:

less_nasty(t) = 1 / (1 + exp(1e3 * (1 - t)))

Or generalizing:

less_nasty(t, nastiness, step_time) = 1 / (1 + exp(nastiness * (step_time - t)))

With a closure:

less_nasty(nastiness, step_time) = t -> less_nasty(t, nastiness, step_time)
1 Like
6

Yeah… this is … something that hasn’t worked out so well. We need a better solution on this problem. Maybe just make if_else and go with it, but that’s not great.

7

Nice! For steps that works perfectly. I guess with a bit of care you could also do pulses in a similar manner. Useful as a backup hack.

1 Like
8

Thanks! That clarifies things a lot.

Unfortunately, while

prob = create(nasty)

works for your construction of nasty(),

prob |> solve

returns a long error message that I’m not sure how to interpret. I’m happy to use continuous hacks for now. Thanks again!

9

post it?

It’s probably even better, because it’ll jumps actually cause order loss in ODE solvers, so relaxing them will improve simulations.

1 Like
10

By long I mean really long. Doesn’t fit on my terminal long! I’ll post two snippets? The first is below, and looks like a portion of an assembly language representation of my function. Most of the error message looks similar. The second is the end of the error message, which might be more enlightening?

Is that enough?
Cheers

 begin
            $(Expr(:inbounds, true))
            local var"#1060#val" = begin
                        #= /.julia/packages/ModelingToolkit/ErYw5/src/build_function.jl:252 =#
                        let (pos, vel, k, c, m, t) = (var"##MTKArg#412"[1], var"##MTKArg#412"[2], var"##MTKArg#413"[1], var"##MTKArg#413"[2], var"##MTKArg#413"[3], var"##MTKArg#414")
                            begin
                                var"##MTIIPVar#416"[1] = vel
                                var"##MTIIPVar#416"[2] = (ModelingToolkit).:-((ModelingToolkit).if_else((ModelingToolkit).:>(t, 1), 1, 0), (ModelingToolkit).:*((ModelingToolkit).:/(1, m), (ModelingToolkit).:+((ModelingToolkit).:*(c, vel), (ModelingToolkit).:*(k, pos))))
                            end
                        end
                    end
            $(Expr(:inbounds, :((ModelingToolkit).pop)))
            var"#1060#val"
        end
        #= /.julia/packages/ModelingToolkit/ErYw5/src/build_function.jl:69 =#
        (ModelingToolkit).nothing
    end
end),:function}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Array{Symbol,1},Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}},OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true,ModelingToolkit.var"#f#86"{GeneralizedGenerated.NGG.RuntimeFn{TypeEncoding(list(##MTKArg#412, ##MTKArg#413, ##MTKArg#414)),TypeEncoding(nil(GeneralizedGenerated.NGG.Argument)),TypeEncoding(begin

Second

 end),:function}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Array{Symbol,1},Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}) at /.julia/packages/DiffEqBase/TiYWr/src/solve.jl:66
 [18] top-level scope at REPL[30]:1
 [19] eval(::Module, ::Any) at ./boot.jl:331
 [20] eval_user_input(::Any, ::REPL.REPLBackend) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [21] run_backend(::REPL.REPLBackend) at /.julia/packages/Revise/MgvIv/src/Revise.jl:1023
 [22] top-level scope at none:0
11

That’s not the right spot.

12

I think the right spot is the beginning of the error message. But I can’t see it because the output saturates my terminal. Even the portion of the error message displayed in my terminal can’t be shown here, as it breaks the character limit. I could send a .txt file to you with the message? Or just the script that generated it, which is below. Sorry about that.

using ModelingToolkit
using DifferentialEquations


@inline function if_else(condition::Bool, @nospecialize(trueval), @nospecialize(falseval))
           return ifelse(condition, trueval, falseval)
end

ModelingToolkit.@register if_else(x, trueval, falseval)


nasty(i) = if_else(i > 1, 1, 0)
@variables t
nasty(t)


function create(input)
    
  @parameters t
  @parameters k,c,m
  @derivatives D'~t
  @variables pos(t) vel(t) inp(t)

  eqs = [D(pos) ~ vel,
        D(vel) ~ input(t) - (1/m)*(c*vel + k*pos),
  ]

  ps = [k,c,m] .=>  [1.,4.,1.] 
  ics = [pos, vel] .=> [1.,0.]
  od = ODESystem(eqs, t, first.(ics),first.(ps))
  tspan = (0.,100.)
  prob = ODEProblem(od, ics, tspan, ps)
  return prob
end

prob = create(nasty)
sol = solve(prob)
13

For me, this works:

using ModelingToolkit
using DifferentialEquations


@inline function if_else(condition::Bool, @nospecialize(trueval), @nospecialize(falseval))
           return ifelse(condition, trueval, falseval)
end

ModelingToolkit.@register if_else(x, trueval, falseval)
# This is needed to bring if_else into MTK's scope.
@eval ModelingToolkit using ..Main: if_else

nasty(i) = if_else(i > 1, 1, 0)
@variables t
nasty(t)


function create(input)
    
  @parameters t
  @parameters k,c,m
  @derivatives D'~t
  @variables pos(t) vel(t) inp(t)

  eqs = [D(pos) ~ vel,
        D(vel) ~ input(t) - (1/m)*(c*vel + k*pos),
  ]

  ps = [k,c,m] .=>  [1.,4.,1.] 
  ics = [pos, vel] .=> [1.,0.]
  od = ODESystem(eqs, t, first.(ics),first.(ps))
  tspan = (0.,100.)
  prob = ODEProblem(od, ics, tspan, ps)
  return prob
end

prob = create(nasty)
sol = solve(prob)

sol

1 Like
14

Ignore me, sorry. It works fine! My REPL was dirty I think. When I refreshed it, there was no error :slight_smile:

15

Indeed, it works for me too. It was a stupid bug due to (I think) a dirty REPL. Thanks a lot for all your help!

16

Hello,

I just wonder is there a better solution now?

I have a similar problem like the original post, but I need a square wave input function.

I have tried the method discussed above.

@inline function if_else(condition::Bool, @nospecialize(trueval), @nospecialize(falseval))
           return ifelse(condition, trueval, falseval)
end

ModelingToolkit.@register if_else(x, trueval, falseval)
@eval ModelingToolkit using ..Main: if_else

If the input is a step function, it works well. Like this.

nasty(i) = if_else(i > 1, 1, 0)

However, if I limit the range.

nasty(i) = if_else(1 < i < 2, 1, 0)

It throws an error.

TypeError: non-boolean (Operation) used in boolean context

I’m excited to use ModelingToolkit.jl to do some simulations, but I don’t know how to solve this. If someone would help, I’d be very grateful!

17

What about if you define

step(x) = if_else( a < x, 1, 0) * if_else(x < b, 1, 0)
18

Yes, it works!

It is really a simple solution, many thanks!