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: 622-element Vector{Float64}:
 0.0
 1.0e-6
 3.059507904487559e-5
 0.00010493013533892183
 0.00023059148752658964
 0.00051150251599221
 0.00117462969063216
 0.002471919421168622
 0.0038470675125998333
 0.00526081718946869
 ⋮
 4.892351055174196
 4.909461452154816
 4.926587225645396
 4.9431954338308985
 4.958649098354077
 4.972589520532182
 4.983574014858123
 4.992417461965787
 5.0
u: 622-element Vector{Vector{Float64}}:
 [1.0, 0.0, 0.0, 0.0, 0.0]
 [1.0, -2.9929944814927005e-18, -4.899999999907105e-12, -9.800000000000003e-6, -5.706217857357636e-16]
 [1.0, 5.27874353047921e-13, -4.58668844216687e-9, -0.00029983177463981217, -1.0032893840539915e-10]
 [0.9999999999999986, -1.1002571697151832e-11, -5.395063295315326e-8, -0.0010283153263207926, -6.916581938886145e-9]
 [0.999999999999966, 1.4396814835659737e-11, -2.605449271805573e-7, -0.002259796577761725, -8.908046718541218e-8]
 [0.9999999999991783, -2.1028415337145945e-10, -1.282010636128357e-6, -0.005012724656719096, -1.0782454291482401e-6]
 [0.9999999999771458, -3.803393446868107e-9, -6.760799062088173e-6, -0.011511370968192219, -1.3743965341272258e-5]
 [0.9999999995517715, -8.15099012139997e-8, -2.9940889534893144e-5, -0.024224810325656064, -0.00013084621883102933]
 [0.9999999973704503, -4.721003658109905e-7, -7.251964931606626e-5, -0.03770126160100114, -0.000496673977913053]
 [0.9999999908045084, -1.6735627425532897e-6, -0.00013561336631286557, -0.05155600830590682, -0.001274409586924759]
 ⋮
 [0.9932496873686112, -0.17497914402984224, -0.11599594342680954, -1.498726354436502, -3.4124883031445647]
 [0.9897216816612306, -0.23950279538625663, -0.14300696938905583, -1.6579652088134131, -4.207665691945125]
 [0.9849679680196319, -0.31795541278587364, -0.17273709142499097, -1.8132347497381565, -5.084823029200119]
 [0.9789568820642802, -0.40834689056711326, -0.2040672070676705, -1.9586319388961435, -6.0123467654684175]
 [0.9719120687463936, -0.5059122597311336, -0.23534429190938105, -2.088082899998507, -6.942656516406527]
 [0.9641838004793751, -0.6055826492653628, -0.2652350105336083, -2.1989126843510167, -7.836901017233438]
 [0.9570721684571752, -0.6921617778389031, -0.28984975083423736, -2.281537466286295, -8.577933213146204]
 [0.9506349430881555, -0.7671203374232561, -0.31031149339058156, -2.3446162406373703, -9.197640315608416]
 [0.9445764375225798, -0.8351672872606196, -0.32829160698944565, -2.395992636889389, -9.745289637586431]

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