Solving DAEs with Physics-Informed Neural Networks (PINNs)

This tutorial is an introduction to using physics-informed neural networks (PINNs) for solving differential algebraic equations (DAEs).

Solving an DAE with PINNs

Let's solve a simple DAE system:

using NeuralPDE, Random, OrdinaryDiffEq, Statistics, Lux, OptimizationOptimisers

example = (du, u, p, t) -> [cos(2pi * t) - du[1], u[2] + cos(2pi * t) - du[2]]
u₀ = [1.0, -1.0]
du₀ = [0.0, 0.0]
tspan = (0.0, 1.0)

(0.0, 1.0)

Differential_vars is a logical array which declares which variables are the differential (non-algebraic) vars

differential_vars = [true, false]

2-element Vector{Bool}:
 1
 0

prob = DAEProblem(example, du₀, u₀, tspan; differential_vars = differential_vars)
chain = Lux.Chain(Lux.Dense(1, 15, cos), Lux.Dense(15, 15, sin), Lux.Dense(15, 2))
opt = OptimizationOptimisers.Adam(0.1)
alg = NNDAE(chain, opt; autodiff = false)
sol = solve(prob, alg, verbose = true, dt = 1 / 100.0, maxiters = 3000, abstol = 1e-10)

retcode: Success
Interpolation: 1st order linear
t: 0.0:0.01:1.0
u: 101-element Vector{Vector{Float64}}:
 [1.0, -1.0]
 [1.0104371745734326, -0.9809714151006407]
 [1.0206948793337574, -0.9603745847290079]
 [1.03075136985688, -0.9381357096073094]
 [1.0405848482158406, -0.9141887609005457]
 [1.0501734609532696, -0.8884761220595233]
 [1.059495306523612, -0.8609491994044184]
 [1.068528452801886, -0.8315689965540741]
 [1.0772509651692859, -0.8003066478678361]
 [1.0856409455864178, -0.7671439061922686]
 ⋮
 [0.922840989871401, -0.8758496540539453]
 [0.9317701463289463, -0.900229674606199]
 [0.9409795308331141, -0.9211569671930178]
 [0.9504300972968065, -0.9386775769693866]
 [0.9600808508541842, -0.9528610336540203]
 [0.9698890031656155, -0.9637989970420656]
 [0.9798101597366279, -0.9716036833674075]
 [0.9897985394803436, -0.9764060940783827]
 [0.9998072260914305, -0.9783540711849884]

Now let's compare the predictions from the learned network with the ground truth which we can obtain by numerically solving the DAE.

function example1(du, u, p, t)
    du[1] = cos(2pi * t)
    du[2] = u[2] + cos(2pi * t)
    nothing
end
M = [1.0 0.0; 0.0 0.0]
f = ODEFunction(example1, mass_matrix = M)
prob_mm = ODEProblem(f, u₀, tspan)
ground_sol = solve(prob_mm, Rodas5(), reltol = 1e-8, abstol = 1e-8)

retcode: Success
Interpolation: specialized 4rd order "free" stiffness-aware interpolation
t: 9-element Vector{Float64}:
 0.0
 1.0e-6
 1.1e-5
 0.00011099999999999999
 0.0011109999999999998
 0.011110999999999996
 0.11111099999999996
 0.7838612882520571
 1.0
u: 9-element Vector{Vector{Float64}}:
 [1.0, -1.0]
 [1.0000010000000001, -0.9999999999802608]
 [1.0000109999999913, -0.9999999976115558]
 [1.0001109999910016, -0.9999997567932182]
 [1.0011109909770284, -0.9999756355789903]
 [1.011101976768569, -0.9975640989587868]
 [1.1023025821686998, -0.7660448918691988]
 [0.8238691280829226, -0.2111552898963995]
 [0.9794681556857161, -1.0]

using Plots
plot(ground_sol, tspan = tspan, layout = (2, 1), label = "ground truth")
plot!(sol, tspan = tspan, layout = (2, 1), label = "dae with pinns")