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?

Initialization By Example: The Cartesian Pendulum

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

using ModelingToolkit, OrdinaryDiffEq, Plots
using ModelingToolkit: t_nounits as t, D_nounits as D

@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( x\left( t \right) \right)^{2} - \left( y\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 \mathtt{xˍtt}\left( t \right) x\left( t \right) - 2 \left( \mathtt{xˍt}\left( t \right) \right)^{2} - 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.BasicSymbolic{Real}}:
 -2xˍt(t)*x(t) - 2y(t)*yˍt(t)
 -g + y(t)*λ(t)
 -2xˍtt(t)*x(t) - 2(xˍt(t)^2) - 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}:
 -0.0
 -1.0
 -2.220446049250313e-16
  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]

1.1102230246251565e-16

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.BasicSymbolic{Real}}:
 -2xˍt(t)*x(t) - 2y(t)*yˍt(t)
 -g + y(t)*λ(t)
 -2xˍtt(t)*x(t) - 2(xˍt(t)^2) - 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 and thus 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 4rd order "free" stiffness-aware interpolation
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Vector{Float64}}:
 [0.0, 0.0, 1.0, 2.0, 1.0]

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.

Constant constraints in initialization

Consider the pendulum system again:

julia> equations(pend)5-element Vector{Equation}:
 0 ~ 1 - (x(t)^2) - (y(t)^2)
 Differential(t)(y(t)) ~ yˍt(t)
 0 ~ -2xˍt(t)*x(t) - 2y(t)*yˍt(t)
 Differential(t)(yˍt(t)) ~ -g + y(t)*λ(t)
 0 ~ -2xˍtt(t)*x(t) - 2(xˍt(t)^2) - 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 4rd 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.045681552850823336
 0.0944235554294302
 0.15548707001985734
 0.229318268820542
 0.31703803193266283
 0.41933953211430064
 0.5364992864870883
 0.6683292127136459
 0.8130876436721504
 0.9642240959164611
 1.1153605481607718
 1.2664970004050826
 1.3990265143903713
 1.5
u: 20-element Vector{Vector{Float64}}:
 [0.0, 0.0007886435331230285, 1.0, 0.0, 0.000788643042647863]
 [7.886432873711667e-10, 0.0007886435326230288, 0.9999996890206408, -9.99999378041376e-7, 0.0007886435316230296]
 [8.675075501083682e-9, 0.0007886434726230662, 0.9999996890206881, -1.0999993158455952e-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.0006034621608964476]
 [-1.1637759953727926e-5, -0.00025475829465943824, 0.9999999675489266, -0.04568154471913162, -0.00234156201837766]
 [-0.0003464634132308097, -0.003669250299838975, 0.9999932682695065, -0.09442281343936325, -0.012585038507219523]
 [-0.0017568578272789631, -0.011299186829808253, 0.9999361620387821, -0.15547531786542068, -0.0354748496109648]
 [-0.005847599419163959, -0.02550150203099147, 0.9996747829798569, -0.229229355413298, -0.07808179832385762]
 [-0.015671835992493246, -0.04944382990318128, 0.9987769011652969, -0.3165743809458416, -0.14990877676108802]
 [-0.03645622909215795, -0.08700238513267104, 0.9962080823980212, -0.417436178924844, -0.2625843273346279]
 [-0.07631915536658758, -0.14254336466355616, 0.9897884840435016, -0.5299425363191282, -0.4292064266557376]
 [-0.14654672643858765, -0.22035617464769985, 0.9754192731561115, -0.6487000338705867, -0.6626405133928107]
 [-0.25956238582822094, -0.3227038201335274, 0.9464996019603908, -0.7613160739284774, -0.9696662121407758]
 [-0.4196794727446767, -0.44465870144405056, 0.8956996825709975, -0.8454137454889569, -1.3354829191620285]
 [-0.6180429404913603, -0.575659836607492, 0.8176894175646954, -0.8779513913305704, -1.7283683476895555]
 [-0.8402321902732645, -0.706575830416946, 0.7076416992265765, -0.8416319086659056, -2.1207160155035387]
 [-1.0355622254188266, -0.8123393364027227, 0.5831983963104594, -0.7436449412352276, -2.4374171964346516]
 [-1.1711852892441768, -0.8818183190779684, 0.4716031686853881, -0.6265051407608526, -2.645770821188183]

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 4rd 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.886435331228014e-5, 0.0007886435331230285, 1.0, -0.1, -0.009211356957380068]
 [7.885516632767382e-5, 0.000788543532623025, 0.9999996890995002, -0.10000100000726413, -0.009211662687967071]
 [7.87630347287516e-5, 0.0007875434726225892, 0.9999996898875912, -0.10001100007986786, -0.009214662867968381]
 [7.784006805955089e-5, 0.0007775373725784321, 0.9999996977177714, -0.10011100080212508, -0.00924468116810085]
 [6.8444722873874e-5, 0.0006769263683004033, 0.9999997708853197, -0.10111100763297366, -0.009546514180934945]
 [-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.15093609309640452]
 ⋮
 [-0.12414778702003189, -0.19561177128055143, 0.9806812997023608, -0.6224044077882067, -0.5984100908092208]
 [-0.21974564219515536, -0.2871116064702765, 0.9578968653386936, -0.7331489040112288, -0.8729003121723828]
 [-0.36017628495441634, -0.3998794409635339, 0.9165673529054944, -0.8255832224794488, -1.2111747663614187]
 [-0.5398472940549877, -0.5243237832731898, 0.8515187098642965, -0.8767679660863378, -1.58444080570272]
 [-0.7480118779350776, -0.652159501606796, 0.7580827888456892, -0.8695809417578748, -1.9677513923114498]
 [-0.9660265597229094, -0.7738710224363448, 0.6333547327628859, -0.7908157668314592, -2.3321450051510846]
 [-1.134008280458917, -0.8614551129283096, 0.5078555528504797, -0.6687593111582794, -2.594293926979205]
 [-1.2480067428966068, -0.9184524219835762, 0.39556080452929626, -0.537715045385103, -2.764788819826563]
 [-1.284368287303742, -0.9359394736870904, 0.3521609714335672, -0.4832659563443791, -2.8190834123605817]

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.886435331228014e-5
julia> prob.ps[Initial(D(x))] = 1.01.0
julia> sol3 = solve(prob, Rodas5P())retcode: Success
Interpolation: specialized 4rd 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.886435331228014e-5, 0.0007886435331230285, 1.0, -0.1, -0.009211356957380068]
 [7.885516632767382e-5, 0.000788543532623025, 0.9999996890995002, -0.10000100000726413, -0.009211662687967071]
 [7.87630347287516e-5, 0.0007875434726225892, 0.9999996898875912, -0.10001100007986786, -0.009214662867968381]
 [7.784006805955089e-5, 0.0007775373725784321, 0.9999996977177714, -0.10011100080212508, -0.00924468116810085]
 [6.8444722873874e-5, 0.0006769263683004033, 0.9999997708853197, -0.10111100763297366, -0.009546514180934945]
 [-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.15093609309640452]
 ⋮
 [-0.12414778702003189, -0.19561177128055143, 0.9806812997023608, -0.6224044077882067, -0.5984100908092208]
 [-0.21974564219515536, -0.2871116064702765, 0.9578968653386936, -0.7331489040112288, -0.8729003121723828]
 [-0.36017628495441634, -0.3998794409635339, 0.9165673529054944, -0.8255832224794488, -1.2111747663614187]
 [-0.5398472940549877, -0.5243237832731898, 0.8515187098642965, -0.8767679660863378, -1.58444080570272]
 [-0.7480118779350776, -0.652159501606796, 0.7580827888456892, -0.8695809417578748, -1.9677513923114498]
 [-0.9660265597229094, -0.7738710224363448, 0.6333547327628859, -0.7908157668314592, -2.3321450051510846]
 [-1.134008280458917, -0.8614551129283096, 0.5078555528504797, -0.6687593111582794, -2.594293926979205]
 [-1.2480067428966068, -0.9184524219835762, 0.39556080452929626, -0.537715045385103, -2.764788819826563]
 [-1.284368287303742, -0.9359394736870904, 0.3521609714335672, -0.4832659563443791, -2.8190834123605817]
julia> sol3[D(x), 1]7.886435331228014e-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 4rd 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.8150000000000002e-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.08488928690714088]
 [-0.017575856123422627, -0.05365382034934039, 0.99855959063321, -0.32710659895301336, -0.16096145764241876]
 [-0.04029240863590199, -0.09328356754840549, 0.9956395561082689, -0.43005070234954135, -0.2798505450196839]
 [-0.083339720839394, -0.15143450049258636, 0.9884672084192943, -0.5439881359656944, -0.4543022888723173]
 [-0.15846448549354739, -0.2324234876915769, 0.9726144534875725, -0.663124059624317, -0.6972639369758995]
 [-0.2776999534582109, -0.33784261624181144, 0.9412022005369218, -0.7736637555227581, -1.0135011116166446]
 [-0.4434091980583462, -0.4615400003136145, 0.8871190239541645, -0.8523040762825664, -1.3845397042688328]
 [-0.6461934443048096, -0.5932846176259522, 0.8049931159076568, -0.876847669913449, -1.779635751500051]
 [-0.8701878089975688, -0.7235561564783292, 0.6902715103027571, -0.8303091501346732, -2.169966459576492]
 [-1.0599348226416343, -0.8253195019734907, 0.56468157724535, -0.7254058974588831, -2.4746548451999186]
 [-1.1714440369830121, -0.882117934381316, 0.47103299825688594, -0.6255903872788648, -2.6456670826618733]
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.

Initialization of parameters

Parameters may also be treated as unknowns in the initialization system. Doing so works almost identically to the standard case. For a parameter to be an initialization unknown (henceforth referred to as "solved parameter") it must represent a floating point number (have a symtype of Real or <:AbstractFloat) or an array of such numbers. Additionally, it must have a guess and one of the following conditions must be satisfied:

1. The value of the parameter as passed to ODEProblem is an expression involving other variables/parameters. For example, if [p => 2q + x] is passed to ODEProblem. In this case, p ~ 2q + x is used as an equation during initialization.
2. The parameter has a default (and no value for it is given to ODEProblem, since that is condition 1). The default will be used as an equation during initialization.
3. The parameter has a default of missing. If ODEProblem is given a value for this parameter, it is used as an equation during initialization (whether the value is an expression or not).
4. ODEProblem is given a value of missing for the parameter. If the parameter has a default, it will be used as an equation during initialization.

All parameter dependencies (where the dependent parameter is a floating point number or array thereof) also become equations during initialization, and the dependent parameters become unknowns.

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.

The variable on the left hand side of all parameter dependencies also has an Initial variant, which is used if a constant constraint is provided for the variable.

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;
    defaults = [total => missing], guesses = [total => 1.0])

\[ \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} - x\left( t \right) - y\left( t \right) \end{align} \]

The system is now overdetermined. In fact:

[equations(initsys); observed(initsys)]

\[ \begin{align} 0 &= \mathtt{total} - x\left( t \right) - y\left( t \right) \\ \mathtt{total} &= 4 \\ y\left( t \right) &= Initial\left( y\left( t \right) \right) \\ x\left( t \right) &= Initial\left( x\left( t \right) \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} \mathtt{total} &= 4 \\ y\left( t \right) &= Initial\left( y\left( t \right) \right) \\ x\left( t \right) &= \mathtt{total} - y\left( t \right) \end{align} \]

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} 0 &= 1 - \left( x\left( t \right) \right)^{2} - \left( y\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 \mathtt{xˍtt}\left( t \right) x\left( t \right) - 2 \left( \mathtt{xˍt}\left( t \right) \right)^{2} - 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) \\ y\left( t \right) &= Initial\left( y\left( t \right) \right) \\ x\left( t \right) &= Initial\left( x\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)

\[ \begin{align} 0 &= 1 - \left( x\left( t \right) \right)^{2} - \left( y\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 \mathtt{xˍtt}\left( t \right) x\left( t \right) - 2 \left( \mathtt{xˍt}\left( t \right) \right)^{2} - 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) \\ y\left( t \right) &= Initial\left( y\left( t \right) \right) \\ x\left( t \right) &= Initial\left( x\left( t \right) \right) \\ \mathtt{xˍtt}\left( t \right) &= x\left( t \right) \lambda\left( t \right) \end{align} \]

unknowns(isys)

6-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 xˍt(t)
 y(t)
 x(t)
 yˍt(t)
 λ(t)
 xˍ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 &= 1 - \left( x\left( t \right) \right)^{2} - \left( y\left( t \right) \right)^{2} \\ 0 &= - \mathtt{xˍtt}\left( t \right) + x\left( t \right) \lambda\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) \\ 0 &= - 2 \mathtt{xˍtt}\left( t \right) x\left( t \right) - 2 \left( \mathtt{xˍt}\left( t \right) \right)^{2} - 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} \]

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} 0 &= 1 - \left( x\left( t \right) \right)^{2} - \left( y\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 \mathtt{xˍtt}\left( t \right) x\left( t \right) - 2 \left( \mathtt{xˍt}\left( t \right) \right)^{2} - 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) \\ y\left( t \right) &= Initial\left( y\left( t \right) \right) \\ x\left( t \right) &= Initial\left( x\left( t \right) \right) \\ \mathtt{yˍt}\left( t \right) &= Initial\left( \mathtt{yˍt}\left( t \right) \right) \\ \lambda\left( t \right) &= Initial\left( \lambda\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)

\[ \begin{align} 0 &= 1 - \left( x\left( t \right) \right)^{2} - \left( y\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 \mathtt{xˍtt}\left( t \right) x\left( t \right) - 2 \left( \mathtt{xˍt}\left( t \right) \right)^{2} - 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} \]

equations(isys)

\[ \begin{align} 0 &= 1 - \left( x\left( t \right) \right)^{2} - \left( y\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 \mathtt{xˍtt}\left( t \right) x\left( t \right) - 2 \left( \mathtt{xˍt}\left( t \right) \right)^{2} - 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} \]

unknowns(isys)

1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 xˍt(t)

observed(isys)

\[ \begin{align} y\left( t \right) &= Initial\left( y\left( t \right) \right) \\ \mathtt{yˍt}\left( t \right) &= Initial\left( \mathtt{yˍt}\left( t \right) \right) \\ \lambda\left( t \right) &= Initial\left( \lambda\left( t \right) \right) \\ x\left( t \right) &= Initial\left( x\left( t \right) \right) \\ \mathtt{xˍtt}\left( t \right) &= x\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: 1-element Vector{Float64}:
 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: 1-element Vector{Float64}:
 0.0

We can see that the default solver stalls

sol.stats

SciMLBase.NLStats
Number of function evaluations:                    276
Number of Jacobians created:                       76
Number of factorizations:                          0
Number of linear solves:                           231
Number of nonlinear solver iterations:             198

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: StalledSuccess
u: 1-element Vector{Float64}:
 0.0

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

sol.resid

3-element Vector{Float64}:
   0.0
  -0.0
 -10.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: 1-element Vector{Float64}:
 0.0

gives a NonlinearLeastSquaresProblem which can be solved:

sol = solve(iprob)

retcode: Success
u: 1-element Vector{Float64}:
 0.0

sol.resid

3-element Vector{Float64}:
  0.0
 -0.0
  0.0

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

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

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

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

sol = solve(iprob)

retcode: Success
u: 4-element Vector{Float64}:
 0.0
 0.0
 1.1593499698618233e-16
 5.180172573250178e-17

sol = solve(iprob, TrustRegion())

retcode: MaxIters
u: 4-element Vector{Float64}:
 0.0
 0.0
 0.42507071480642905
 0.1350143071877229

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

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)

(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)

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 8-element Vector{Float64}:
  0.0
  0.5203838316719773
  2.5786944719549396
  5.8760572779107605
  7.485534094375067
  8.596602375120318
  9.644019757516785
 10.0
u: 8-element Vector{Vector{Float64}}:
 [2.7920470927578725e-21, 5.354980650644457e-21]
 [5.982284964269159e-22, 1.1688314292610062e-20]
 [2.3815285323636405e-20, 2.4689818424178436e-19]
 [1.977363579511111e-17, 2.689066514306813e-17]
 [1.5804438258000038e-16, 2.9730504130951645e-16]
 [1.151372944343759e-16, 1.5716107268327463e-15]
 [5.89020855341433e-17, 7.554600579626735e-15]
 [2.0265659178800287e-17, 1.2885802663228949e-14]

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())

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 40-element Vector{Float64}:
  0.0
  0.07205457770130536
  0.19555442536537837
  0.34017815218853753
  0.5197315269179786
  0.7296409614368518
  0.9859476024857492
  1.3125688506305393
  1.6001373094389693
  1.9779709033648663
  ⋮
  7.638661765365626
  7.922169811926683
  8.211731307428632
  8.566084321082675
  8.910984213988325
  9.277801668292575
  9.72291901385274
  9.995336042698368
 10.0
u: 40-element Vector{Vector{Float64}}:
 [1.5, 0.5]
 [1.252874833326255, 0.5046075911240547]
 [0.9219106728662867, 0.5315611614577552]
 [0.6476528652853633, 0.590222615065694]
 [0.42415789793088654, 0.7028731420169776]
 [0.26693504796941986, 0.8969640810827857]
 [0.16238845397704782, 1.2488273595868689]
 [0.10189942173478911, 1.955280057819758]
 [0.08599023235056862, 2.9316165842420534]
 [0.11950416087725532, 4.983888436619836]
 ⋮
 [0.8424187326340977, 0.5449087676403777]
 [0.4282819270088137, 0.7013256443322156]
 [0.22864669024338416, 0.9881483587275364]
 [0.12337903584443591, 1.584339864178958]
 [0.08841783056090323, 2.565927273418115]
 [0.10096000448919976, 4.30441614658586]
 [0.3691037566200744, 7.73597698628302]
 [1.7498047811154076, 9.214228119813201]
 [1.801213120529808, 9.202400764435039]

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

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

40-element Vector{Float64}:
  0.0
  0.1247012350613691
  0.3072898341575289
  0.5030745547950242
  0.7561805185954678
  1.1060149716134986
  1.6704458951727528
  2.7336781795082388
  4.145334485121325
  6.880237249405644
  ⋮
  0.3583217980237515
  0.751623368083024
  1.256285686398821
  2.1810354713761324
  3.6220171872344893
  6.022050346396097
  8.748587312661025
 -2.3017582386173707
 -2.7718838506213714

plot(sol)

Summary of Initialization API

Initial(x)

generate_initializesystem(
    sys::ModelingToolkit.AbstractSystem;
    time_dependent_init,
    kwargs...
) -> System

initialization_equations(
    sys::ModelingToolkit.AbstractSystem
) -> Any

guesses(sys::ModelingToolkit.AbstractSystem) -> Any

defaults(sys::ModelingToolkit.AbstractSystem) -> Any