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: 624-element Vector{Float64}:
 0.0
 1.0e-6
 3.059507904487559e-5
 0.00010493045200655489
 0.00023059678317840222
 0.0005115254907802978
 0.001174694409152639
 0.0024719939762102086
 0.0038471397651965236
 0.005260880989641164
 ⋮
 4.890863932231844
 4.907929741507869
 4.925098271831778
 4.941832049260109
 4.957493787237935
 4.971684178661713
 4.982979239992738
 4.992096239271686
 5.0
u: 624-element Vector{Vector{Float64}}:
 [1.0, 0.0, 0.0, 0.0, 0.0]
 [1.0, -2.99272169122073e-18, -4.899999999907118e-12, -9.800000000000003e-6, -5.706331226400594e-16]
 [1.0, 5.27809507593611e-13, -4.5866884421638045e-9, -0.0002998317746398121, -1.003294826634931e-10]
 [0.9999999999999986, -1.0999713035493748e-11, -5.39509585878918e-8, -0.0010283184296635964, -6.916598893733934e-9]
 [0.999999999999966, 1.4392125308898031e-11, -2.6055689441382537e-7, -0.0022598484751494875, -8.908536209054332e-8]
 [0.999999999999178, -2.1029113238084374e-10, -1.2821258050016937e-6, -0.0050129498096423575, -1.0783669302428895e-6]
 [0.9999999999771407, -3.804235958610987e-9, -6.7615440815092435e-6, -0.011512005209692906, -1.3745848457455767e-5]
 [0.9999999995517175, -8.15163470832287e-8, -2.9942695643873477e-5, -0.02422554096506339, -0.00013085389212311434]
 [0.9999999973702528, -4.7212180396378426e-7, -7.25223733556092e-5, -0.03770196967644684, -0.0004966856576259517]
 [0.9999999908040624, -1.673586706276006e-6, -0.00013561665561498618, -0.051556633547591116, -0.001274413563898522]
 ⋮
 [0.9935035021243901, -0.17002816849890137, -0.11380154462489456, -1.48477321197939, -3.3475521946432196]
 [0.9900802595596732, -0.23322623501093956, -0.14050295370893892, -1.6439008772545352, -4.133433377366214]
 [0.9854315446238384, -0.31060522651676, -0.1700725477966511, -1.799968874658277, -5.0055617613515135]
 [0.9795027585564285, -0.4004256432787717, -0.2014307514664329, -1.946957206503656, -5.933511478915601]
 [0.9724854280387512, -0.49820859066293754, -0.23296372403908652, -2.078661692434194, -6.870957885547281]
 [0.9647211038546726, -0.5988198658879121, -0.26327400643794197, -2.191937575232685, -7.777289593102764]
 [0.9574736087841169, -0.6873431238434594, -0.2885208857882799, -2.2772075015159965, -8.537018179737075]
 [0.950871548768589, -0.7643680677173623, -0.3095857153945488, -2.342407571114353, -9.174806464041236]
 [0.9445669296022601, -0.8352201602780062, -0.3283189638647778, -2.3960152209267056, -9.745349295144326]

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