Solving DAEs with Probabilistic Numerics

ProbNumDiffEq.jl provides probabilistic numerical solvers for differential algebraic equations (DAEs). Currently, we recommend using the semi-implicit EK1 algorithm.

Solving mass-matrix DAEs with the EK1

First, define the DAE (here the ROBER problem) as an ODE problem with singular mass matrix:

using ProbNumDiffEq, Plots, LinearAlgebra, OrdinaryDiffEq, ModelingToolkit, LinearAlgebrafunction rober(du, u, p, t)    y₁, y₂, y₃ = u    k₁, k₂, k₃ = p    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃    du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2    du[3] = y₁ + y₂ + y₃ - 1    nothingendM = [1 0 0     0 1 0     0 0 0]f = ODEFunction(rober, mass_matrix=M)prob_mm = ODEProblem(f, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: true
timespan: (0.0, 100000.0)
u0: 3-element Vector{Float64}:
 1.0
 0.0
 0.0

We can solve this problem directly with the EK1:

sol = solve(prob_mm, EK1(), reltol=1e-8, abstol=1e-8)plot(    sol,    xscale=:log10,    tspan=(1e-6, 1e5),    layout=(3, 1),    legend=false,    ylabel=["u₁(t)" "u₂(t)" "u₃(t)"],    xlabel=["" "" "t"],    denseplot=false,)

Looks good!

Solving an Index-3 DAE directly

The following is based on the "Automatic Index Reduction of DAEs" tutorial by ModelingToolkit.jl, which demonstrates how the classic Rodas4 solver fails to solve a DAE due to the fact that it is of index 3; which is why ModelingToolkit's automatic index reduction is so useful.

It turns out that our probabilistic numerical solvers can directly solve the index-3 DAE!

First, define the pendulum problem as in the tutorial:

function pendulum!(du, u, p, t)    x, dx, y, dy, T = u    g, L = p    du[1] = dx    du[2] = T * x    du[3] = dy    du[4] = T * y - g    du[5] = x^2 + y^2 - L^2endpendulum_fun! = ODEFunction(pendulum!, mass_matrix=Diagonal([1, 1, 1, 1, 0]))u0 = [1.0, 0, 0, 0, 0];p = [9.8, 1];tspan = (0, 5.0);pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)

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

We can try to solve it directly with one of the classic mass-matrix DAE solvers from OrdinaryDiffEq.jl:

solve(pendulum_prob, Rodas4())

It does not work! This is because of the index of the DAE; see for example this explanation from the tutorial.

Does this also hold for the EK1 solver? Let's find out:

sol = solve(pendulum_prob, EK1())

retcode: Success
Interpolation: ODE Filter Posterior
t: 628-element Vector{Float64}:
 0.0
 1.0e-6
 3.059507904487559e-5
 0.00010493045200655489
 0.00023059678317840222
 0.0005115254907802978
 0.001174695488853979
 0.002472014211685773
 0.00384718896152788
 0.005261036144077238
 ⋮
 4.8962217947660385
 4.913405922736698
 4.930471041386539
 4.946961849708434
 4.962239003120593
 4.976007705962295
 4.9868494820907285
 4.99559623753655
 5.0
u: 628-element Vector{Vector{Float64}}:
 [1.0, 0.0, 0.0, 0.0, 0.0]
 [1.0, -5.2308124708835185e-17, -4.89999999990727e-12, -9.800000000000003e-6, -1.6582319408395278e-10]
 [1.0, 3.1300721393858876e-13, -4.586688442131688e-9, -0.00029983177463981244, -5.274705680235624e-9]
 [0.9999999999999986, -1.1765605240836134e-11, -5.39509585882545e-8, -0.0010283184296635915, -2.5528866190765293e-8]
 [0.999999999999966, 9.628759432573492e-12, -2.6055689441388254e-7, -0.002259848475149491, -1.3316810259149302e-7]
 [0.9999999999991781, -2.366607199924195e-10, -1.2821258050029018e-6, -0.005012949809642313, -1.1915391164539253e-6]
 [0.9999999999771407, -3.975985200826331e-9, -6.761556511038636e-6, -0.011512015790765483, -1.4084318333704079e-5]
 [0.9999999995517028, -8.252603038829201e-8, -2.9943185861260923e-5, -0.024225739272706853, -0.00013185101689948278]
 [0.9999999973701182, -4.7513850688521257e-7, -7.252422816609226e-5, -0.037702451800374515, -0.0004986227217320749]
 [0.9999999908029774, -1.6802425589217303e-6, -0.0001356246549744592, -0.05155815406055371, -0.0012775339386260874]
 ⋮
 [0.9925445275968695, -0.18845407147785992, -0.12188257074933884, -1.5350824921989754, -3.5851653149174094]
 [0.9887414815879484, -0.25632438523009793, -0.14963382979693712, -1.6941696800235013, -4.401999301823816]
 [0.9836917339321452, -0.33780483124018806, -0.17986265158236167, -1.8477770471952495, -5.29372972594089]
 [0.9773733303623875, -0.4308846746658905, -0.21152157175644473, -1.9907895169818173, -6.23085124037384]
 [0.9700488383026391, -0.5304419743723015, -0.24290997723953475, -2.1172544563963207, -7.164335281053306]
 [0.9620667288232236, -0.6316345137110422, -0.2728142633344204, -2.225160843469619, -8.058818442411416]
 [0.9547591024729455, -0.7191903733006798, -0.2973803468546597, -2.3053562837880506, -8.798218063947152]
 [0.9481513899824944, -0.7949726451672193, -0.3178190686396982, -2.366569319139831, -9.417028559032959]
 [0.9445705007069478, -0.8348820679894055, -0.328308653012859, -2.396101456925807, -9.73606169960135]

Nope! The EK1 is able to solve the index-3 DAE directly. Pretty cool!

plot(sol)

Is index-reduction still worth it?

The point of the "Automatic Index Reduction of DAEs" tutorial is to demonstrate ModelingToolkit's utility for automatic index reduction, which enables the classic implicit Runge-Kutta solvers such as Rodas5 to solve this DAE. Let's see if that still helps in this context here.

First, modelingtoolkitize the problem:

traced_sys = modelingtoolkitize(pendulum_prob)

Model ##MTKizedODE#168:
Equations (5):
  5 standard: see equations(##MTKizedODE#168)
Unknowns (5): see unknowns(##MTKizedODE#168)
  x₁(t)
  x₂(t)
  x₃(t)
  x₄(t)
  ⋮
Parameters (2): see parameters(##MTKizedODE#168)
  α₁
  α₂

(how cool is this latex output ?!?)

Next, lower the DAE index and simplify it with ModelingToolkit's dae_index_lowering and structural_simplify:

simplified_sys = structural_simplify(dae_index_lowering(traced_sys))

Model ##MTKizedODE#168:
Equations (5):
  5 standard: see equations(##MTKizedODE#168)
Unknowns (5): see unknowns(##MTKizedODE#168)
  x₄(t)
  x₂(t)
  x₁(t)
  x₃(t)
  ⋮
Parameters (2): see parameters(##MTKizedODE#168)
  α₁
  α₂

Let's build two different ODE problems, and check how well we can solve each:

prob_index3 = ODEProblem(complete(traced_sys), Pair[], tspan)
prob_index1 = ODEProblem(simplified_sys, Pair[], tspan)

sol3 = solve(prob_index3, EK1(), abstol=1e-7)
sol1 = solve(prob_index1, EK1(), abstol=1e-7)

truesol = solve(prob_index1, Rodas4(), abstol=1e-10, reltol=1e-10)

sol1_final_error = norm(sol1.u[end] - truesol.u[end])
sol1_f_evals     = sol1.stats.nf
perm = [findfirst(x -> x === s, unknowns(traced_sys)) for s in unknowns(simplified_sys)]
sol3_final_error = norm(sol3.u[end][perm] - truesol.u[end])
sol3_f_evals     = sol3.stats.nf
@info "Results" sol1_final_error sol1_f_evals sol3_final_error sol3_f_evals

The error for the index-1 DAE solve is quite a bit lower. Thus it seems that, even if the index-3 DAE could also be solved directly, index lowering might still be beneficial when solving DAEs with the EK1!

References