Solving Stiff-ODE using PINN via NNODE

Specific Domains Modelling & Simulations
,
1

Hi I have been trying to solve the ODE below:

dc/dt = - Acexp(-E/(R*T(t)))

Then I scaled time and Temperature it turned into:

dc/dt_ = - Actim_maxexp(-E/(RT_(t)*T_max))

!Note that the temperature increases through time slowly for 30 minutes then increases fast till it reaches a peak then reduces.

!Note that the values of everything is in the code.

I got the Temperature and time data, therefore I used Akima Interpolator
to define T as a function of t (T(t)). The problem I am having is that my residuals never get’s under 1e4 and the results are false. I have tried some technics like transforming the equation into log to make the stiffness less but I have failed since the real c would be small Negative numbers at some points. I assume the problematic term is the exp term since it produces very small values near to 0.
How can I solve this stiff problem?

using XLSX
using DataInterpolations
using Plots
using Lux, Random, NeuralPDE
using LuxCUDA, ComponentArrays
using CUDA
using DifferentialEquations
using Plots
const gdev = CUDA.device()
#########################################
#Reading the data
xf = XLSX.readxlsx("avetempsec.xlsx")
xf
XLSX.sheetnames(xf)
sh = xf["Sheet1"]
T = sh[:,2]
#########################################
#time data generation 0 to 15000 seconds
l = []
for i in range(0,15000)
    push!(l,i)
end
println(length(l))
tim = Array(l)
reshape(tim,(15001,1))
#########################################
#constant values of ODE
A_e = 5.14e25
E_e = 2.7e5
Rc = 8.314
#########################################
#defining time and Temperature
T = vec(T)
tim = vec(tim)
T = convert(Vector{Float64}, T)
tim = convert(Vector{Float64}, tim)
##########################################
T_max = maximum(T) # is 590 K
tim_max = maximum(tim) # is 15000 seconds
T_scaled = T/T_max
tim_scaled = tim/tim_max # initial (Min) temperature is 293.15 K
##########################################
# Defining Interpolation for defining T as a function of t
A = AkimaInterpolation(T_scaled, tim_scaled)
scatter(tim_scaled, T_scaled, label = "input data")
plot!(A)
#########################################
#defining ODE
u0 = 1
tspan = (minimum(tim_scaled), maximum(tim_scaled))
f(u, p, t) =  - A_e * u * tim_max * exp(-E_e/(Rc*A(t)*T_max))
prob = ODEProblem(f, u0, tspan)
#########################################
#define the Network
rng = Random.default_rng()
Random.seed!(rng, 0)
chain = Chain(Dense(1, 10, Lux.tanh),Dense(10, 10, Lux.tanh),Dense(10, 10, Lux.tanh), Dense(10, 1))
ps, st = Lux.setup(rng, chain)
ps = ps |> ComponentArray |> gpu .|> Lux.Float64
#########################################
using OptimizationOptimisers
opt = Adam(0.001)
alg = NNODE(chain, opt, init_params = ps)
sol = solve(prob, alg, verbose = true, maxiters = 1000,saveat= 0.01)

Thanks for your time. @ChrisRackauckas

2

It seems crazy to me to solve a first-order ODE like this this with a PINN. Not only would a regular ODE solver chew right through this, but isn’t this particular ODE separable? i.e. you have

\frac{dc}{c} = -A e^{-E/RT(t)} dt \implies \log c(t) = \log c(0) - \int_0^t A e^{-E/RT(t')} dt'

so you don’t even need an ODE solver: it boils down to a simple integral that you can perform with any quadrature package, like QuadGK.jl.

If T(t) is interpolated from data, then whether you are using a quadrature routine or an ODE solver you will want to break the integral into pieces at the knots of the interpolation (i.e. the data points).

Sorry if this is not directly answering your question about PINNs, but my personal bias is that the only people who should currently be using PINNs are researchers who want to study PINNs, and who are therefore already experts in NNs and differential equations. If you just want to solve a differential equation, classic methods are currently almost always faster and easier.

4 Likes
3

Many thanks for your reply.
I am trying to solve a coupled system consists of one PDE and one ODE. The ODE is as I said. My PDE is a Transient Energy balance equation (Heat-diffusion equation with heat generation) in which the heat generation is:

Constant*-dc/dt

I tried to solve the coupled system but the MSE won’t be great and my problem gives wrong answers. So I somehow wanted to debug it by this method.
I have the data and I know the solutions. When I gave the PDE the Qgen, PINN was able to solve the problem but vice versa it’s not true because I know that PINNs are susceptible to stiffness and the ODE is the problematic term.

4

Have you considered method of lines, with a standard spatial discretization for your PDE (e.g. finite differences or finite elements)? That will almost certainly work better than a PINN.

1 Like
5

PINNs are known to be unstable with stiffness. There’s some papers that go into this:

https://pubs.acs.org/doi/10.1021/acs.jpca.1c05102

Etc. There are things we can do to further mitigate it in the library, but as @stevengj alludes, if you have a hard problem then a PINN is almost certainly not the right way to solve it unless you need a very specific surrogate (but there’s also better surrogates)

1 Like
6

Many Thanks for your reply.
Unfortunately my project is to solve this coupled system by PINN.
I know the difficulties to solve such equations by PINN.
I have already solved it using:

sol = solve(prob, AutoTsit5(Rosenbrock23()), reltol = 1e-8, saveat = 0.001)

and the solution showed a great fit.

I read the paper sent by you and it was a great help. What I understood from the paper for solving these equations are using deeper networks and implementing learning rate annealing algorithms. But I have already solved the problem by ReLoBRaLo and it wasn’t correct.
Is it logical to use the data to solve the ODEs by solvers just like the one I said and solve the PDE by PINN?

7

It’s not a bad idea at all. Pre-training the neural network to match the ODE solver data, and then finishing with the PINN train, could be a good way to stabilize the training and it’s a technique we have used before.

1 Like