Initialization of Systems

While for simple numerical ODEs choosing an initial condition can be an easy affair, with ModelingToolkit's more general differential-algebraic equation (DAE) system there is more care needed due to the flexibility of the solver state. In this tutorial we will walk through the functionality involved in initialization of System and the diagnostics to better understand and debug the initialization problem.

Primer on Initialization of Differential-Algebraic Equations

Before getting started, let's do a brief walkthrough of the mathematical principles of initialization of DAE systems. Take a DAE written in semi-explicit form:

\[\begin{aligned} x^\prime &= f(x,y,t) \\ 0 &= g(x,y,t) \end{aligned}\]

where $x$ are the differential variables and $y$ are the algebraic variables. An initial condition $u0 = [x(t_0) y(t_0)]$ is said to be consistent if $g(x(t_0),y(t_0),t_0) = 0$.

For ODEs, this is trivially satisfied. However, for more complicated systems it may not be easy to know how to choose the variables such that all of the conditions are satisfied. This is even more complicated when taking into account ModelingToolkit's simplification engine, given that variables can be eliminated and equations can be changed. If this happens, how do you know how to initialize the system?

Bindings and Initial Conditions

ModelingToolkit represents initialization constraints in two different ways: bindings and initial conditions. Both of these are stored as stored as key-value pairs, similar to defaults in prior versions of ModelingToolkit. Bindings represent immutable relations between variables (and parameters) of a system. They cannot be changed after constructing a system. Initial conditions are simply a convenience to avoid repeatedly passing the same values to problem constructors.

using ModelingToolkit, Test
using ModelingToolkit: t_nounits as t, D_nounits as D

@variables x(t) y(t)
@parameters p q
@mtkcompile sys = System([D(x) ~ p * x + y, y ~ q * x], t; bindings = [y => x, q => 2p],
                         initial_conditions = [x => 1, p => 1])

\[ \begin{align} \frac{\mathrm{d} ~ x\left( t \right)}{\mathrm{d}t} &= y\left( t \right) + p ~ x\left( t \right) \end{align} \]

julia> bindings(sys)ReadOnlyDicts.ReadOnlyDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, ModelingToolkitBase.AtomicArrayDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}}} with 2 entries:
  q    => 2p
  y(t) => x(t)
julia> @test_throws Exception bindings(sys)[y] = 2x # throwsTest Passed
      Thrown: MethodError
julia> initial_conditions(sys)ModelingToolkitBase.AtomicArrayDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}} with 2 entries:
  x(t) => 1
  p    => 1
julia> initial_conditions(sys)[x] = 22
julia> initial_conditions(sys)ModelingToolkitBase.AtomicArrayDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}} with 2 entries:
  x(t) => 2
  p    => 1

Variables/parameters which have an entry in bindings are referred to as "bound" variables/ parameters. The term "bound symbolics" is used to refer to such symbolic variables, regardless of whether they are variables or parameters.

Bindings for variables can be constants or functions of other variables/parameters. Bindings for parameters can only be constants or functions of other parameters. It is often useful to bind a parameter to the initial value of a variable. In this case, the Initial operator can be used. In the above example, q => Initial(y) is a valid binding, whereas q => y is not.

Bound symbolics cannot be given initial conditions. No symbolic can have an entry in both bindings and initial_conditions. Since initial_conditions is a convenience for passing values to problem constructors, bound symbolics cannot be given values in the problem constructor either. The exception is solvable parameters, which have a binding of missing. Solving for parameter values is covered in more detail in a later section.

julia> ODEProblem(sys, [], (0.0, 1.0)) # fine┌ Warning: Initialization system is overdetermined. 2 equations for 1 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here.
│ 
│ 
│ To suppress this warning, pass `warn_initialize_determined = false`. To turn this warning into an error, pass `fully_determined = true`.
└ @ ModelingToolkit ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/initialization.jl:22
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: OVERDETERMINED
Non-trivial mass matrix: false
timespan: (0.0, 1.0)
u0: 1-element Vector{Float64}:
 2.0
julia> ODEProblem(sys, [x => 2, p => 2.5]) # fineERROR: MethodError: no method matching ODEProblem(::System, ::Vector{Pair{Num, Float64}})

Closest candidates are:
  ODEProblem(::System, ::AbstractArray{<:Pair}, ::Any, ::AbstractDict; kw...)
   @ ModelingToolkitBase ~/work/ModelingToolkit.jl/ModelingToolkit.jl/lib/ModelingToolkitBase/src/deprecations.jl:49
  ODEProblem(::System, ::AbstractArray, ::Any, ::AbstractDict; kw...)
   @ ModelingToolkitBase ~/work/ModelingToolkit.jl/ModelingToolkit.jl/lib/ModelingToolkitBase/src/deprecations.jl:49
  ODEProblem(::System, ::AbstractArray{<:Pair}, ::Any, ::AbstractArray{<:Pair}; kw...)
   @ ModelingToolkitBase ~/work/ModelingToolkit.jl/ModelingToolkit.jl/lib/ModelingToolkitBase/src/deprecations.jl:49
  ...
julia> @test_throws Exception ODEProblem(sys, [y => 1]) # errorsTest Passed
      Thrown: MethodError
julia> @test_throws Exception ODEProblem(sys, [q => 2]) # errorsTest Passed
      Thrown: MethodError

Bindings and initial conditions should not be cyclic. In other words, the binding or initial condition for a symbolic should not (directly or indirectly) be a function of itself.

Bindings for variables are treated equivalently to initial conditions when building problems. They form constraints in the initialization system. Bindings for floating-point-valued discrete variables (created via @discretes and used in events) are also treated in the same way. Bindings for non-floating-point discrete variables (such as ones with integer types) do not turn into constraints in the initialization system, since this would lead to mixed-integer problems. Bindings for these variables are still used to calculate initial conditions as if they were parameters. Bindings for parameters participate in initialization. The bound parameters are treated as unknowns of the initialization system.

Initialization By Example: The Cartesian Pendulum

To illustrate how to perform the initialization, let's take a look at the Cartesian pendulum:

using OrdinaryDiffEq, Plots

@parameters g
@variables x(t) y(t) [state_priority = 10] λ(t)
eqs = [D(D(x)) ~ λ * x
       D(D(y)) ~ λ * y - g
       x^2 + y^2 ~ 1]
@mtkcompile pend = System(eqs, t)

\[ \begin{align} 0 &= 1 - \left( y\left( t \right) \right)^{2} - \left( x\left( t \right) \right)^{2} \\ \frac{\mathrm{d} ~ y\left( t \right)}{\mathrm{d}t} &= \mathtt{yˍt}\left( t \right) \\ 0 &= - 2 ~ \mathtt{xˍt}\left( t \right) ~ x\left( t \right) - 2 ~ y\left( t \right) ~ \mathtt{yˍt}\left( t \right) \\ \frac{\mathrm{d} ~ \mathtt{yˍt}\left( t \right)}{\mathrm{d}t} &= - g + y\left( t \right) ~ \lambda\left( t \right) \\ 0 &= - 2 ~ \left( \mathtt{xˍt}\left( t \right) \right)^{2} - 2 ~ x\left( t \right) ~ \mathtt{xˍtt}\left( t \right) - 2 ~ \left( \mathtt{yˍt}\left( t \right) \right)^{2} - 2 ~ y\left( t \right) ~ \left( - g + y\left( t \right) ~ \lambda\left( t \right) \right) \end{align} \]

While we defined the system using second derivatives and a length constraint, the structural simplification system improved the numerics of the system to be solvable using the dummy derivative technique, which results in 3 algebraic equations and 2 differential equations. In this case, the differential equations with respect to y and D(y), though it could have just as easily have been x and D(x). How do you initialize such a system if you don't know in advance what variables may defined the equation's state?

To see how the system works, let's start the pendulum in the far right position, i.e. x(0) = 1 and y(0) = 0. We can do this by:

prob = ODEProblem(pend, [x => 1, y => 0, g => 1], (0.0, 1.5), guesses = [λ => 1])

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: true
timespan: (0.0, 1.5)
u0: 5-element Vector{Float64}:
 0.0
 0.0
 1.0
 0.0
 1.0

This solves via:

sol = solve(prob, Rodas5P())
plot(sol, idxs = (x, y))

and we can check it satisfies our conditions via:

conditions = getfield.(equations(pend)[3:end], :rhs)

3-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}:
 -2xˍt(t)*x(t) - 2y(t)*yˍt(t)
 -g + y(t)*λ(t)
 -2(xˍt(t)^2) - 2x(t)*xˍtt(t) - 2(yˍt(t)^2) - 2y(t)*(-g + y(t)*λ(t))

[sol[conditions][1]; sol[x][1] - 1; sol[y][1]]

5-element Vector{Float64}:
  3.013752158750194e-21
 -1.0
 -4.462500713999945e-8
  0.0
  0.0

Notice that we set [x => 1, y => 0] as our initial conditions and [λ => 1] as our guess. The difference is that the initial conditions are required to be satisfied, while the guesses are simply a guess for what the initial value might be. Every variable must have either an initial condition or a guess, and thus since we did not know what λ would be we set it to 1 and let the initialization scheme find the correct value for λ. Indeed, the value for λ at the initial time is not 1:

sol[λ][1]

2.2312503569999724e-8

We can similarly choose λ = 0 and solve for y to start the system:

prob = ODEProblem(pend, [x => 1, λ => 0, g => 1], (0.0, 1.5); guesses = [y => 1])
sol = solve(prob, Rodas5P())
plot(sol, idxs = (x, y))

or choose to satisfy derivative conditions:

prob = ODEProblem(
    pend, [x => 1, D(y) => 0, g => 1], (0.0, 1.5); guesses = [λ => 0, y => 1])
sol = solve(prob, Rodas5P())
plot(sol, idxs = (x, y))

Notice that since a derivative condition is given, we are required to give a guess for y.

We can also directly give equations to be satisfied at the initial point by using the initialization_eqs keyword argument, for example:

prob = ODEProblem(pend, [x => 1, g => 1], (0.0, 1.5); guesses = [λ => 0, y => 1],
    initialization_eqs = [y ~ 0])
sol = solve(prob, Rodas5P())
plot(sol, idxs = (x, y))

Additionally, note that the initial conditions are allowed to be functions of other variables and parameters:

prob = ODEProblem(
    pend, [x => 1, D(y) => g, g => 1], (0.0, 3.0); guesses = [λ => 0, y => 1])
sol = solve(prob, Rodas5P())
plot(sol, idxs = (x, y))

Determinability: Underdetermined and Overdetermined Systems

For this system we have 3 conditions to satisfy:

conditions = getfield.(equations(pend)[3:end], :rhs)

3-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}:
 -2xˍt(t)*x(t) - 2y(t)*yˍt(t)
 -g + y(t)*λ(t)
 -2(xˍt(t)^2) - 2x(t)*xˍtt(t) - 2(yˍt(t)^2) - 2y(t)*(-g + y(t)*λ(t))

when we initialize with

prob = ODEProblem(pend, [x => 1, y => 0, g => 1], (0.0, 1.5); guesses = [y => 0, λ => 1])

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: true
timespan: (0.0, 1.5)
u0: 5-element Vector{Float64}:
 0.0
 0.0
 1.0
 0.0
 1.0

we have two extra conditions to satisfy, x ~ 1 and y ~ 0 at the initial point. That gives 5 equations for 5 variables. However, this is not sufficient for a well-formed system. This initialization is singular, which means that at least one of the initial conditions is redundant and provides no extra information. Here, one of x ~ 1 or y ~ 0 is redundant, since the other can be inferred using the algebraic equation x ^ 2 + y ^ 2 ~ 1. Thus, this is similar to an underdetermined system. To make the system well-formed, we need to give an initial value for a derivative. For example:

prob = ODEProblem(pend, [x => 1, D(y) => 0, g => 1], (0.0, 1.5); guesses = [y => 0, λ => 1])

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: true
timespan: (0.0, 1.5)
u0: 5-element Vector{Float64}:
 0.0
 0.0
 1.0
 0.0
 1.0

Now the system is well-formed. What happens if that's not the case?

prob = ODEProblem(pend, [x => 1, g => 1], (0.0, 1.5); guesses = [y => 0, λ => 1])

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: UNDERDETERMINED
Non-trivial mass matrix: true
timespan: (0.0, 1.5)
u0: 5-element Vector{Float64}:
 0.0
 0.0
 1.0
 0.0
 1.0

Here we have 4 equations for 5 unknowns (note: the warning is post-simplification of the nonlinear system, which solves the trivial x ~ 1 equation analytical and thus says 3 equations for 4 unknowns). This warning thus lets you know the system is underdetermined and thus the solution is not necessarily unique. It can still be solved:

sol = solve(prob, Rodas5P())
plot(sol, idxs = (x, y))

and the found initial condition satisfies all constraints which were given. In the opposite direction, we may have an overdetermined system:

prob = ODEProblem(
    pend, [x => 1, y => 0.0, D(y) => 0, g => 1], (0.0, 1.5); guesses = [λ => 1])

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: OVERDETERMINED
Non-trivial mass matrix: true
timespan: (0.0, 1.5)
u0: 5-element Vector{Float64}:
 0.0
 0.0
 1.0
 0.0
 1.0

Can that be solved?

sol = solve(prob, Rodas5P())
plot(sol, idxs = (x, y))

Indeed since we saw D(y) = 0 at the initial point above, it turns out that this solution is solvable with the chosen initial conditions. However, for overdetermined systems we often aren't that lucky. If the set of initial conditions cannot be satisfied, then you will get a SciMLBase.ReturnCode.InitialFailure:

prob = ODEProblem(
    pend, [x => 1, y => 0.0, D(y) => 2.0, λ => 1, g => 1], (0.0, 1.5); guesses = [λ => 1])
sol = solve(prob, Rodas5P())

retcode: InitialFailure
Interpolation: specialized 4th (Rodas6P = 5th) order "free" stiffness-aware interpolation
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Vector{Float64}}:
 [-4.051733408186758e-13, 0.0, 1.0, 2.0, -3.000000000000004]

What this means is that the initial condition finder failed to find an initial condition. While this can be sometimes due to numerical error (which is then helped by picking guesses closer to the correct value), most circumstances of this come from ill-formed models. Especially if your system is overdetermined and you receive an InitialFailure, the initial conditions may not be analytically satisfiable!. In our case here, if you sit down with a pen and paper long enough you will see that λ = 0 is required for this equation, but since we chose λ = 1 we end up with a set of equations that are impossible to satisfy.

Note

If you would prefer to have an error instead of a warning in the context of non-fully determined systems, pass the keyword argument fully_determined = true into the problem constructor. Additionally, any warning about not being fully determined can be suppressed via passing warn_initialize_determined = false.

Constant constraints in initialization

Consider the pendulum system again:

julia> equations(pend)5-element Vector{Equation}:
 0 ~ 1 - (y(t)^2) - (x(t)^2)
 Differential(t, 1)(y(t)) ~ yˍt(t)
 0 ~ -2xˍt(t)*x(t) - 2y(t)*yˍt(t)
 Differential(t, 1)(yˍt(t)) ~ -g + y(t)*λ(t)
 0 ~ -2(xˍt(t)^2) - 2x(t)*xˍtt(t) - 2(yˍt(t)^2) - 2y(t)*(-g + y(t)*λ(t))
julia> observed(pend)1-element Vector{Equation}:
 xˍtt(t) ~ x(t)*λ(t)

Suppose we want to solve the same system with multiple different initial y-velocities from a given position.

prob = ODEProblem(
    pend, [x => 1, D(y) => 0, g => 1], (0.0, 1.5); guesses = [λ => 0, y => 1, x => 1])
sol1 = solve(prob, Rodas5P())

retcode: Success
Interpolation: specialized 4th (Rodas6P = 5th) order "free" stiffness-aware interpolation
t: 20-element Vector{Float64}:
 0.0
 1.0e-6
 1.1e-5
 0.00011099999999999999
 0.0011109999999999998
 0.011110999999999996
 0.045681559052122296
 0.0944235656195696
 0.15548708021915836
 0.22931828686064443
 0.31703805361025933
 0.4193395580340552
 0.5364993172651313
 0.6683292489583085
 0.8130876859195416
 0.9642241302333386
 1.1153605745471356
 1.2664970188609326
 1.3990265258920607
 1.5
u: 20-element Vector{Vector{Float64}}:
 [-0.0, 0.0007886435331230285, 1.0, 0.0, 0.0007886430426196882]
 [7.886432873711593e-10, 0.0007886435326230288, 0.9999996890206408, -9.999993780413668e-7, 0.0007886435316230295]
 [8.675075501083674e-9, 0.0007886434726230662, 0.9999996890206881, -1.0999993158455941e-5, 0.0007886433516231414]
 [8.753872113897968e-8, 0.00078863737262686, 0.9999996890254989, -0.00011099993096331238, 0.0007886250516345232]
 [8.754970282908355e-7, 0.0007880263730066766, 0.9999996895071696, -0.0011109993097247108, 0.0007867920527739726]
 [8.076765836480156e-6, 0.0007269164090608092, 0.9999997357962322, -0.011110993785204729, 0.0006034621608964477]
 [-1.1637774474464156e-5, -0.00025475857794436663, 0.9999999675488545, -0.045681550920426896, -0.0023415628682325307]
 [-0.000346463541473337, -0.0036692512620206245, 0.999993268265976, -0.0944228236290318, -0.012585041393764722]
 [-0.001756858189074605, -0.011299188415547803, 0.9999361620208633, -0.1554753280606326, -0.03547485436818331]
 [-0.005847600827312833, -0.02550150616631248, 0.9996747828743651, -0.22922937341747782, -0.07808181072982302]
 [-0.015671839238185957, -0.04944383676575335, 0.9987769008255678, -0.3165744024627593, -0.1499087973487985]
 [-0.03645623587248301, -0.08700239595251781, 0.9962080814530779, -0.4174362042524294, -0.26258435979409683]
 [-0.07631916844181787, -0.1425433809741752, 0.9897884816945197, -0.5299425652140476, -0.4292064755871063]
 [-0.1465467498651394, -0.2203561981597509, 0.9754192678444635, -0.6487000648224333, -0.6626405839264237]
 [-0.25956242460064916, -0.3227038522977162, 0.9464995909941124, -0.7613161029541278, -0.9696663086227203]
 [-0.4196795137961398, -0.44465873045635007, 0.8956996681682624, -0.8454137594255228, -1.3354830061996288]
 [-0.6180429777885209, -0.5756598597733689, 0.8176894012556061, -0.8779513914621667, -1.7283684171965688]
 [-0.8402322179888526, -0.706575845949237, 0.7076416837169174, -0.8416318994642342, -2.120716062142544]
 [-1.0355622418085537, -0.8123393449546973, 0.5831983843961042, -0.7436449299603659, -2.4374172222085932]
 [-1.1711852894113337, -0.8818183190756633, 0.4716031686803357, -0.6265051407709424, -2.6457708217136915]

sol1[D(y), 1]

0.0

Repeatedly re-creating the ODEProblem with different values of D(y) and x or repeatedly calling remake is slow. Instead, for any variable => constant constraint in the ODEProblem initialization (whether provided to the ODEProblem constructor or a default value) we can update the constant value. ModelingToolkit refers to these values using the Initial operator. For example:

prob.ps[[Initial(x), Initial(D(y))]]

2-element Vector{Float64}:
 1.0
 0.0

To solve with a different starting y-velocity, we can simply do

prob.ps[Initial(D(y))] = -0.1
sol2 = solve(prob, Rodas5P())

retcode: Success
Interpolation: specialized 4th (Rodas6P = 5th) order "free" stiffness-aware interpolation
t: 21-element Vector{Float64}:
 0.0
 1.0e-6
 1.1e-5
 0.00011099999999999999
 0.0011109999999999998
 0.011110999999999996
 0.041215093898905854
 0.08993330820016343
 0.14812222179944373
 0.22329866507033358
 ⋮
 0.5370414123084308
 0.671743376923075
 0.8159140935032908
 0.961394822363531
 1.1068755512237711
 1.2523562800840113
 1.3716583013796808
 1.4657632985175937
 1.5
u: 21-element Vector{Vector{Float64}}:
 [7.886435331230286e-5, 0.0007886435331230285, 1.0, -0.1, -0.009211356957380313]
 [7.885516632767382e-5, 0.000788543532623025, 0.9999996890995002, -0.10000100000726413, -0.009211662687967071]
 [7.87630347287516e-5, 0.0007875434726225892, 0.9999996898875912, -0.10001100007986786, -0.00921466286796838]
 [7.784006805955089e-5, 0.0007775373725784321, 0.9999996977177714, -0.10011100080212508, -0.00924468116810085]
 [6.8444722873874e-5, 0.0006769263683004033, 0.9999997708853197, -0.10111100763297366, -0.009546514180934943]
 [-4.26870660757205e-5, -0.00038418387425569374, 0.999999926201371, -0.11111102201026807, -0.012729844908642285]
 [-0.0005905889703810618, -0.004182196076764525, 0.9999912545780697, -0.14121379765139605, -0.024123881544559737]
 [-0.0023263484222510893, -0.01224834994352976, 0.9999249860967343, -0.18991730127525644, -0.04832234358298631]
 [-0.006200453362596942, -0.024990709673709748, 0.9996876831826028, -0.24803279435426354, -0.08654942324772835]
 [-0.015015386284505857, -0.046452933838936464, 0.998920477739512, -0.32288960619098006, -0.1509360930964045]
 ⋮
 [-0.12414778702003171, -0.19561177128055085, 0.980681299702361, -0.6224044077882076, -0.5984100908092214]
 [-0.21974564219515424, -0.2871116064702748, 0.957896865338694, -0.7331489040112299, -0.8729003121723823]
 [-0.36017628495441484, -0.3998794409635319, 0.9165673529054953, -0.8255832224794502, -1.2111747663614179]
 [-0.5398472940549871, -0.5243237832731887, 0.8515187098642971, -0.8767679660863394, -1.584440805702721]
 [-0.7480118779350768, -0.6521595016067951, 0.7580827888456901, -0.8695809417578764, -1.96775139231145]
 [-0.9660265597229108, -0.7738710224363448, 0.633354732762886, -0.7908157668314607, -2.3321450051510895]
 [-1.1340082804589167, -0.8614551129283091, 0.50785555285048, -0.6687593111582794, -2.594293926979204]
 [-1.2480067428966022, -0.9184524219835755, 0.39556080452929804, -0.5377150453851036, -2.764788819826551]
 [-1.2843682873037399, -0.9359394736870899, 0.3521609714335687, -0.4832659563443807, -2.819083412360577]

sol2[D(y), 1]

-0.1

Note that this only applies for constant constraints for the current ODEProblem. For example, D(x) does not have a constant constraint - it is solved for by initialization. Thus, mutating Initial(D(x)) does not have any effect:

julia> sol2[D(x), 1]7.886435331230286e-5
julia> prob.ps[Initial(D(x))] = 1.01.0
julia> sol3 = solve(prob, Rodas5P())retcode: Success
Interpolation: specialized 4th (Rodas6P = 5th) order "free" stiffness-aware interpolation
t: 21-element Vector{Float64}:
 0.0
 1.0e-6
 1.1e-5
 0.00011099999999999999
 0.0011109999999999998
 0.011110999999999996
 0.041215093898905854
 0.08993330820016343
 0.14812222179944373
 0.22329866507033358
 ⋮
 0.5370414123084308
 0.671743376923075
 0.8159140935032908
 0.961394822363531
 1.1068755512237711
 1.2523562800840113
 1.3716583013796808
 1.4657632985175937
 1.5
u: 21-element Vector{Vector{Float64}}:
 [7.886435331230286e-5, 0.0007886435331230285, 1.0, -0.1, -0.009211356957380313]
 [7.885516632767382e-5, 0.000788543532623025, 0.9999996890995002, -0.10000100000726413, -0.009211662687967071]
 [7.87630347287516e-5, 0.0007875434726225892, 0.9999996898875912, -0.10001100007986786, -0.00921466286796838]
 [7.784006805955089e-5, 0.0007775373725784321, 0.9999996977177714, -0.10011100080212508, -0.00924468116810085]
 [6.8444722873874e-5, 0.0006769263683004033, 0.9999997708853197, -0.10111100763297366, -0.009546514180934943]
 [-4.26870660757205e-5, -0.00038418387425569374, 0.999999926201371, -0.11111102201026807, -0.012729844908642285]
 [-0.0005905889703810618, -0.004182196076764525, 0.9999912545780697, -0.14121379765139605, -0.024123881544559737]
 [-0.0023263484222510893, -0.01224834994352976, 0.9999249860967343, -0.18991730127525644, -0.04832234358298631]
 [-0.006200453362596942, -0.024990709673709748, 0.9996876831826028, -0.24803279435426354, -0.08654942324772835]
 [-0.015015386284505857, -0.046452933838936464, 0.998920477739512, -0.32288960619098006, -0.1509360930964045]
 ⋮
 [-0.12414778702003171, -0.19561177128055085, 0.980681299702361, -0.6224044077882076, -0.5984100908092214]
 [-0.21974564219515424, -0.2871116064702748, 0.957896865338694, -0.7331489040112299, -0.8729003121723823]
 [-0.36017628495441484, -0.3998794409635319, 0.9165673529054953, -0.8255832224794502, -1.2111747663614179]
 [-0.5398472940549871, -0.5243237832731887, 0.8515187098642971, -0.8767679660863394, -1.584440805702721]
 [-0.7480118779350768, -0.6521595016067951, 0.7580827888456901, -0.8695809417578764, -1.96775139231145]
 [-0.9660265597229108, -0.7738710224363448, 0.633354732762886, -0.7908157668314607, -2.3321450051510895]
 [-1.1340082804589167, -0.8614551129283091, 0.50785555285048, -0.6687593111582794, -2.594293926979204]
 [-1.2480067428966022, -0.9184524219835755, 0.39556080452929804, -0.5377150453851036, -2.764788819826551]
 [-1.2843682873037399, -0.9359394736870899, 0.3521609714335687, -0.4832659563443807, -2.819083412360577]
julia> sol3[D(x), 1]7.886435331230286e-5

To enforce this constraint, we would have to remake the problem (or construct a new one).

julia> prob2 = remake(prob; u0 = [y => 0.0, D(x) => 0.0, x => nothing, D(y) => nothing]);
julia> sol4 = solve(prob2, Rodas5P())retcode: Success
Interpolation: specialized 4th (Rodas6P = 5th) order "free" stiffness-aware interpolation
t: 20-element Vector{Float64}:
 0.0
 1.0e-6
 1.1e-5
 0.00011099999999999999
 0.0011109999999999998
 0.011110999999999996
 0.047421953721844705
 0.09899191776516113
 0.16167014460606832
 0.23791158268969675
 0.3276728673290381
 0.43231140466948814
 0.5516095548372473
 0.6855651987992282
 0.8318595500535554
 0.9833110142781696
 1.1347624785027837
 1.2862139427273978
 1.4161711905701027
 1.5
u: 20-element Vector{Vector{Float64}}:
 [0.0, 0.0, 1.0, 0.0, 0.0]
 [-4.999999999999893e-19, -4.999999999999966e-13, 1.0, -9.999999999999896e-7, -1.4999999999999716e-12]
 [-6.654999999999996e-16, -6.050000000000005e-11, 1.0, -1.1000000000000023e-5, -1.8150000000000005e-10]
 [-6.838154999999942e-13, -6.1605e-9, 1.0, -0.00011100000000000012, -1.8481499999999906e-8]
 [-6.856653155000046e-10, -6.171604999999509e-7, 0.9999999999998096, -0.0011109999999997437, -1.851481499999888e-6]
 [-6.8585048037596e-7, -6.172716045067835e-5, 0.9999999980948788, -0.011110999974598648, -0.0001851814813919577]
 [-5.3322222786193873e-5, -0.001124420557837535, 0.9999993678387505, -0.04742191774786085, -0.003373261765019013]
 [-0.0004850274022298178, -0.004899676317827757, 0.9999879965002433, -0.09899049186644025, -0.014699029711803393]
 [-0.0021126982846158527, -0.013068171259272267, 0.9999146076613947, -0.16165357828752486, -0.03920451617564782]
 [-0.006731514577681339, -0.028296426994205093, 0.9995995748496768, -0.2377972696176479, -0.08488928690714087]
 [-0.017575856123422627, -0.05365382034934039, 0.99855959063321, -0.32710659895301336, -0.1609614576424188]
 [-0.04029240863590199, -0.09328356754840549, 0.9956395561082689, -0.43005070234954135, -0.2798505450196839]
 [-0.083339720839394, -0.15143450049258636, 0.9884672084192943, -0.5439881359656944, -0.45430228887231733]
 [-0.15846448549354736, -0.2324234876915769, 0.9726144534875725, -0.6631240596243169, -0.6972639369758994]
 [-0.2776999534582109, -0.3378426162418114, 0.9412022005369218, -0.7736637555227582, -1.0135011116166446]
 [-0.44340919805834594, -0.46154000031361414, 0.8871190239541646, -0.8523040762825665, -1.3845397042688325]
 [-0.6461934443048095, -0.5932846176259519, 0.804993115907657, -0.8768476699134494, -1.7796357515000512]
 [-0.8701878089975693, -0.7235561564783292, 0.690271510302757, -0.8303091501346735, -2.169966459576493]
 [-1.0599348226416356, -0.825319501973491, 0.5646815772453493, -0.7254058974588827, -2.4746548451999217]
 [-1.1714440369830115, -0.8821179343813165, 0.4710329982568851, -0.6255903872788631, -2.64566708266187]
julia> sol4[D(x), 1]0.0

Note the need to provide x => nothing, D(y) => nothing to override the previously provided initial conditions. Since remake is a partial update, the constraints provided to it are merged with the ones already present in the problem. Existing constraints can be removed by providing a value of nothing.

Initial variables

The Initial form used above is only available for a specific set of variables. ModelingToolkit allows using the Initial form on:

Unknowns (unknowns(sys)).
Observables (observables(sys)).
First derivatives of unknowns (even if the unknown in question is itself the derivative of another unknown/observable).
First derivatives of observables (likewise).
Discrete variables (created via @discretes and updated via events).
Bound parameters (bound_parameters(sys)).

Initialization of parameters

Parameters may also be treated as unknowns in the initialization system. This is automatically the case for any floating-point-typed (a symtype of Real or <:AbstractFloat) bound parameter. The binding is used as an initialization equation. It is also often useful to solve for parameters that are not explicitly bound to functions of other parameters. For example, a model for the fluid flow in a circular pipe might use the cross-sectional area A of the pipe. For convenience, the model may want to allow specifying either the area A directly or the radius r of the pipe. Binding A => pi * r * r prevents directly specifying the value of A. In such cases, the required parameters can be marked as initialization unknowns by giving them a binding of missing. These parameters are henceforth referred to as "solvable parameters". In the previous example, this would entail passing [A => missing, r => missing] to the bindings keyword of the System constructor. The relation between them can be passed as an equation A ~ pi * r * r to the initialization_eqs keyword of the System constructor.

remake will reconstruct the initialization system and problem, given the new constraints provided to it. The new values will be combined with the original variable-value mapping provided to ODEProblem and used to construct the initialization problem.

Parameter initialization by example

Consider the following system, where the sum of two unknowns is a constant parameter total.

using ModelingToolkit, OrdinaryDiffEq # hidden
using ModelingToolkit: t_nounits as t, D_nounits as D # hidden

@variables x(t) y(t)
@parameters total
@mtkcompile sys = System([D(x) ~ -x, total ~ x + y], t; bindings = [total => missing])

\[ \begin{align} \frac{\mathrm{d} ~ x\left( t \right)}{\mathrm{d}t} &= - x\left( t \right) \end{align} \]

Given any two of x, y and total we can determine the remaining variable.

prob = ODEProblem(sys, [x => 1.0, y => 2.0], (0.0, 1.0))
integ = init(prob, Tsit5())
integ.ps[total]

3.0

Suppose we want to re-create this problem, but now solve for x given total and y:

prob2 = remake(prob; u0 = [y => 1.0], p = [total => 4.0])
initsys = prob2.f.initializeprob.f.sys

\[ \begin{align} 0 &= \mathtt{total} - y\left( t \right) - x\left( t \right) \end{align} \]

The system is now overdetermined. In fact:

[equations(initsys); observed(initsys)]

\[ \begin{align} 0 &= \mathtt{total} - y\left( t \right) - x\left( t \right) \\ y\left( t \right) &= Initial\left( y\left( t \right) \right) \\ \mathtt{total} &= Initial\left( \mathtt{total} \right) \\ x\left( t \right) &= Initial\left( x\left( t \right) \right) \\ \mathtt{xˍt}\left( t \right) &= - x\left( t \right) \end{align} \]

The system can never be satisfied and will always lead to an InitialFailure. This is due to the aforementioned behavior of retaining the original variable-value mapping provided to ODEProblem. To fix this, we pass x => nothing to remake to remove its retained value.

prob2 = remake(prob; u0 = [y => 1.0, x => nothing], p = [total => 4.0])
initsys = prob2.f.initializeprob.f.sys

\[ \begin{align} \end{align} \]

The system is fully determined, and the equations are solvable.

[equations(initsys); observed(initsys)]

\[ \begin{align} y\left( t \right) &= Initial\left( y\left( t \right) \right) \\ \mathtt{total} &= Initial\left( \mathtt{total} \right) \\ x\left( t \right) &= \mathtt{total} - y\left( t \right) \\ \mathtt{xˍt}\left( t \right) &= - x\left( t \right) \end{align} \]

What constitutes an initialization system?

The initialization system considers the following as unknowns:

All unknowns of the system (unknowns(sys)).
First derivatives of all differential variables in the system.
All observables of the system (observables(sys)).
All bound parameters of the system (bound_parameters(sys)).
All parameters with a binding of missing.

The equations of the initialization system are:

All equations of the system (equations(sys)). Differential equations are used to solve for the first derivatives of differential variables.
All observed equations of the system (observed(sys)).
All bindings in the system, excluding bindings for solvable parameters (since the value is missing).
All additional initialization equations (initialization_equations(sys)). This also includes additional equations passed to the initialization_eqs keyword of the problem constructor.
All initial conditions, formed by combinding initial_conditions(sys) with those passed to the problem. Initial conditions passed to the problem override those in the system in case of conflict.

ModelingToolkit's improved simplification and index reduction may also be able to analytically find the derivatives of some algebraic variables, typically ones corresponding to algebraic equations that arise from index reduction. In such a case, these variables (along with the corresponding equations) are also present in the initialization system.

Diving Deeper: Constructing the Initialization System

To get a better sense of the initialization system and to help debug it, you can construct the initialization system directly. The initialization system is a NonlinearSystem which requires the system-level information and the additional nonlinear equations being tagged to the system.

isys = generate_initializesystem(pend; op = [x => 1.0, y => 0.0], guesses = [λ => 1])

\[ \begin{align} \mathtt{yˍtt}\left( t \right) &= - g + y\left( t \right) ~ \lambda\left( t \right) \\ 0 &= 1 - \left( y\left( t \right) \right)^{2} - \left( x\left( t \right) \right)^{2} \\ 0 &= - 2 ~ \mathtt{xˍt}\left( t \right) ~ x\left( t \right) - 2 ~ y\left( t \right) ~ \mathtt{yˍt}\left( t \right) \\ 0 &= - 2 ~ \left( \mathtt{xˍt}\left( t \right) \right)^{2} - 2 ~ x\left( t \right) ~ \mathtt{xˍtt}\left( t \right) - 2 ~ \left( \mathtt{yˍt}\left( t \right) \right)^{2} - 2 ~ y\left( t \right) ~ \left( - g + y\left( t \right) ~ \lambda\left( t \right) \right) \\ x\left( t \right) &= Initial\left( x\left( t \right) \right) \\ y\left( t \right) &= Initial\left( y\left( t \right) \right) \\ \mathtt{xˍtt}\left( t \right) &= x\left( t \right) ~ \lambda\left( t \right) \end{align} \]

We can inspect what its equations and unknown values are:

equations(isys)

unknowns(isys)

7-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}:
 xˍt(t)
 y(t)
 x(t)
 yˍt(t)
 λ(t)
 xˍtt(t)
 yˍtt(t)

Notice that all initial conditions are treated as initial equations. Additionally, for systems with observables, those observables are too treated as initial equations. We can see the resulting simplified system via the command:

isys = mtkcompile(isys; fully_determined = false)

\[ \begin{align} 0 &= - 2 ~ \mathtt{xˍt}\left( t \right) ~ x\left( t \right) - 2 ~ y\left( t \right) ~ \mathtt{yˍt}\left( t \right) \\ 0 &= - 2 ~ \left( \mathtt{xˍt}\left( t \right) \right)^{2} - 2 ~ x\left( t \right) ~ \mathtt{xˍtt}\left( t \right) - 2 ~ \left( \mathtt{yˍt}\left( t \right) \right)^{2} - 2 ~ y\left( t \right) ~ \left( - g + y\left( t \right) ~ \lambda\left( t \right) \right) \\ 0 &= 1 - \left( y\left( t \right) \right)^{2} - \left( x\left( t \right) \right)^{2} \end{align} \]

Note fully_determined=false allows for the simplification to occur when the number of equations does not match the number of unknowns, which we can use to investigate our overdetermined system:

isys = ModelingToolkit.generate_initializesystem(
    pend; op = [x => 1, y => 0.0, D(y) => 2.0, λ => 1], guesses = [λ => 1])

\[ \begin{align} \mathtt{yˍtt}\left( t \right) &= - g + y\left( t \right) ~ \lambda\left( t \right) \\ 0 &= 1 - \left( y\left( t \right) \right)^{2} - \left( x\left( t \right) \right)^{2} \\ 0 &= - 2 ~ \mathtt{xˍt}\left( t \right) ~ x\left( t \right) - 2 ~ y\left( t \right) ~ \mathtt{yˍt}\left( t \right) \\ 0 &= - 2 ~ \left( \mathtt{xˍt}\left( t \right) \right)^{2} - 2 ~ x\left( t \right) ~ \mathtt{xˍtt}\left( t \right) - 2 ~ \left( \mathtt{yˍt}\left( t \right) \right)^{2} - 2 ~ y\left( t \right) ~ \left( - g + y\left( t \right) ~ \lambda\left( t \right) \right) \\ x\left( t \right) &= Initial\left( x\left( t \right) \right) \\ \lambda\left( t \right) &= Initial\left( \lambda\left( t \right) \right) \\ \mathtt{yˍt}\left( t \right) &= Initial\left( \mathtt{yˍt}\left( t \right) \right) \\ y\left( t \right) &= Initial\left( y\left( t \right) \right) \\ \mathtt{xˍtt}\left( t \right) &= x\left( t \right) ~ \lambda\left( t \right) \end{align} \]

isys = mtkcompile(isys; fully_determined = false)

equations(isys)

unknowns(isys)

2-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}:
 λ(t)
 xˍt(t)

observed(isys)

\[ \begin{align} \mathtt{yˍt}\left( t \right) &= Initial\left( \mathtt{yˍt}\left( t \right) \right) \\ x\left( t \right) &= Initial\left( x\left( t \right) \right) \\ y\left( t \right) &= Initial\left( y\left( t \right) \right) \\ \mathtt{xˍtt}\left( t \right) &= x\left( t \right) ~ \lambda\left( t \right) \\ \mathtt{yˍtt}\left( t \right) &= - g + y\left( t \right) ~ \lambda\left( t \right) \end{align} \]

After simplification we see that we have 5 equatinos to solve with 3 variables, and the system that is given is not solvable.

Numerical Isolation: InitializationProblem

To inspect the numerics of the initialization problem, we can use the InitializationProblem constructor which acts just like an ODEProblem or NonlinearProblem constructor, but creates the special initialization system for a given sys. This is done as follows:

iprob = ModelingToolkit.InitializationProblem(pend, 0.0,
    [x => 1, y => 0.0, D(y) => 2.0, λ => 1, g => 1], guesses = [λ => 1])

NonlinearLeastSquaresProblem with uType Vector{Float64}. In-place: true
u0: 2-element Vector{Float64}:
 1.0
 0.0

We can see that because the system is overdetermined we receive a NonlinearLeastSquaresProblem, solvable by NonlinearSolve.jl. Using NonlinearSolve we can recreate the initialization solve directly:

using NonlinearSolve
sol = solve(iprob)

retcode: StalledSuccess
u: 2-element Vector{Float64}:
 3.056507338356202e-11
 0.0

Note

For more information on solving NonlinearProblems and NonlinearLeastSquaresProblems, check out the NonlinearSolve.jl tutorials!.

We can see that the default solver stalls

sol.stats

SciMLBase.NLStats
Number of function evaluations:                    156
Number of Jacobians created:                       68
Number of factorizations:                          0
Number of linear solves:                           188
Number of nonlinear solver iterations:             148

after doing many iterations, showing that it tried to compute but could not find a valid solution. Trying other solvers:

sol = solve(iprob, GaussNewton())

retcode: InternalLinearSolveFailed
u: 2-element Vector{Float64}:
 1.0
 0.0

gives the same issue, indicating that the chosen initialization system is unsatisfiable. We can check the residuals:

sol.resid

4-element Vector{Float64}:
 -0.0
  0.0
  1.0
 -1.0

to see the problem is not equation 2 but other equations in the system. Meanwhile, changing some of the conditions:

iprob = ModelingToolkit.InitializationProblem(pend, 0.0,
    [x => 1, y => 0.0, D(y) => 0.0, λ => 0, g => 1], guesses = [λ => 1])

NonlinearLeastSquaresProblem with uType Vector{Float64}. In-place: true
u0: 2-element Vector{Float64}:
 0.0
 0.0

gives a NonlinearLeastSquaresProblem which can be solved:

sol = solve(iprob)

retcode: InternalLinearSolveFailed
u: 2-element Vector{Float64}:
 0.0
 0.0

sol.resid

4-element Vector{Float64}:
 -0.0
  0.0
  1.0
  0.0

In comparison, if we have a well-conditioned system:

iprob = ModelingToolkit.InitializationProblem(pend, 0.0,
    [x => 1, D(x) => 0.0, g => 1], guesses = [λ => 1, y => 0])

SCCNonlinearProblem with uType Vector{Float64}. In-place: false
u0: 3-element Vector{Float64}:
 0.0
 0.0
 1.0

notice that we instead obtained a NonlinearProblem. In this case we can use different solvers which can take advantage of the fact that the Jacobian is square.

sol = solve(iprob)

retcode: Stalled
u: 3-element Vector{Float64}:
 0.0
 0.0
 0.0

sol = solve(iprob, TrustRegion())

retcode: Stalled
u: 3-element Vector{Float64}:
 0.0
 0.0
 0.0

More Features of the Initialization System: Steady-State and Observable Initialization

Example currently disabled

This section's examples are currently disabled due to a compatibility issue with the initialization system and the current ModelingToolkit stack.

Let's take a Lotka-Volterra system:

@variables x(t) y(t) z(t)
@parameters α=1.5 β=1.0 γ=3.0 δ=1.0

eqs = [D(x) ~ α * x - β * x * y
       D(y) ~ -γ * y + δ * x * y
       z ~ x + y]

@named sys = System(eqs, t)
simpsys = mtkcompile(sys)
tspan = (0.0, 10.0)

Using the derivative initializations, we can set the ODE to start at the steady state by initializing the derivatives to zero:

prob = ODEProblem(simpsys, [D(x) => 0.0, D(y) => 0.0], tspan, guesses = [x => 1, y => 1])
sol = solve(prob, Tsit5(), abstol = 1e-16)

Notice that this is a "numerical zero", not an exact zero, and thus the solution will leave the steady state in this instance because it's an unstable steady state.

Additionally, notice that in this setup we have an observable z ~ x + y. If we instead know the initial condition for the observable we can use that directly:

prob = ODEProblem(simpsys, [D(x) => 0.0, z => 2.0], tspan, guesses = [x => 1, y => 1])
sol = solve(prob, Tsit5())

We can check that indeed the solution does satisfy that D(x) = 0 at the start:

sol[α * x - β * x * y]

plot(sol)

Summary of Initialization API

Missing docstring.

Missing docstring for Initial. Check Documenter's build log for details.

Missing docstring.

Missing docstring for isinitial. Check Documenter's build log for details.

Missing docstring.

Missing docstring for generate_initializesystem. Check Documenter's build log for details.

Missing docstring.

Missing docstring for initialization_equations. Check Documenter's build log for details.

Missing docstring.

Missing docstring for guesses. Check Documenter's build log for details.

Missing docstring.

Missing docstring for bindings. Check Documenter's build log for details.

Missing docstring.

Missing docstring for initial_conditions. Check Documenter's build log for details.