The nn can't converge while solving my odes

[Edward-CY]

Hi, lu Thanks for your great work! I am learning how to solve my odes with initial conditions. It is a complex odes that describe a chemical process. I reference to the example named Lotka_Volterra_Colab. It seems that the nn can't convergence while solving my odes. My code is as follows:

from future import absolute_import from future import division from future import print_function import math from deepxde.backend import tf import numpy as np import matplotlib.pyplot as plt import deepxde as dde

def ode_system(x, y):

Temperature in CSTR (K)

T = y[:, 0:1]
#Concentration of A in CSTR (mol/m^3)
Ca = y[:, 1:2]
dT_t = dde.grad.jacobian(y, x, i=0)
dCa_t = dde.grad.jacobian(y, x, i=1)

Tc = 280  # nominal = 300
# Volumetric Flowrate (m^3/sec)
q = 100   # nominal = 100
# Feed Concentration (mol/m^3)
Caf = 1 # nominal = 1
# Feed Temperature (K)
Tf = 350  # nominal = 350

#Ca_ss = 0.87725294608097;
#T_ss = 324.475443431599;

# Volume of CSTR (m^3)
V = 100
# Density of A-B Mixture (kg/m^3)
rho = 1000
# Heat capacity of A-B Mixture (J/kg-K)
Cp = 0.239
# Heat of reaction for A->B (J/mol)
mdelH = 5e4
# E - Activation energy in the Arrhenius Equation (J/mol)
# R - Universal Gas Constant = 8.31451 J/mol-K
EoverR = 8750
# Pre-exponential factor (1/sec)
k0 = 7.2e10    
# U - Overall Heat Transfer Coefficient (W/m^2-K)
# A - Area - this value is specific for the U calculation (m^2)
UA = 5e4

return [
    dT_t - q/V*(Tf - T) + mdelH/(rho*Cp)*k0*tf.math.exp(-EoverR/T)*Ca + UA/V/rho/Cp*(Tc-T),
    dCa_t - q/V*(Caf - Ca) - k0*tf.math.exp(-EoverR/T)*Ca]

def boundary(_, on_initial): return on_initial

geom = dde.geometry.TimeDomain(0.0, 10.0) ic1 = dde.IC(geom, lambda X: 324.475443431599, boundary, component=0) ic2 = dde.IC(geom, lambda X: 0.87725294608097, boundary, component=1) data = dde.data.PDE(geom, ode_system, [ic1, ic2], 3000, 2, num_test=3000)

layer_size = [1] + [64] * 3 + [2] activation = "tanh" initializer = "Glorot normal" net = dde.maps.FNN(layer_size, activation, initializer)

model = dde.Model(data, net) model.compile("adam", lr=0.001) losshistory, train_state = model.train(epochs=30000)

model.compile("L-BFGS") losshistory, train_state = model.train()

dde.saveplot(losshistory, train_state, issave=True, isplot=True)

After it runs. It shows:

Step Train loss Test loss Test metric 0 [inf, inf, 1.05e+05, 7.70e-01] [inf, inf, 1.05e+05, 7.70e-01] []
1000 [nan, nan, nan, nan] [nan, nan, nan, nan] []
2000 [nan, nan, nan, nan] [nan, nan, nan, nan] []
3000 [nan, nan, nan, nan] [nan, nan, nan, nan] []

I have tried many methods like changing lr or activation function. I wonder what the problem is. Thank you!

[sumantkrsoni]

I think you have some bug in jacobian definition:

dT_t = dde.grad.jacobian(y, x, i=0)---------------------------> dde.grad.jacobian(y, x, i=0, j= 0) dCa_t = dde.grad.jacobian(y, x, i=1)------------------------->dde.grad.jacobian(y, x, i = 1, j=0)

try it.

[Edward-CY]

Thanks for your advice! I tried your method, unfortunately, it fails again:

Initializing variables... Training model...

Step Train loss Test loss Test metric 0 [5.54e+04, 2.20e+00, 1.05e+05, 7.70e-01] [5.54e+04, 2.20e+00, 1.05e+05, 7.70e-01] []
1000 [nan, nan, nan, nan] [nan, nan, nan, nan] []
2000 [nan, nan, nan, nan] [nan, nan, nan, nan] []
3000 [nan, nan, nan, nan] [nan, nan, nan, nan] []

After that, I wonder if I can make the nn converge by feeding the exact solution to the training process. I use the scipy integrate.solve_ivp module to solve the odes, and get the exact solution data. And then I set the data as a kind of BC. This is my code:

from future import absolute_import from future import division from future import print_function import math from deepxde.backend import tf import numpy as np import matplotlib.pyplot as plt import deepxde as dde from scipy import integrate

def ode_system(x, y): # odes used in tf

Temperature in CSTR (K)

T = y[:, 0:1]
#Concentration of A in CSTR (mol/m^3)
Ca = y[:, 1:2]
dT_t = dde.grad.jacobian(y, x, i=0, j=0)
dCa_t = dde.grad.jacobian(y, x, i=1, j=0)

Tc = 280  # nominal = 300
# Volumetric Flowrate (m^3/sec)
q = 100   # nominal = 100
# Feed Concentration (mol/m^3)
Caf = 1 # nominal = 1
# Feed Temperature (K)
Tf = 350  # nominal = 350

#Ca_ss = 0.87725294608097;
#T_ss = 324.475443431599;

# Volume of CSTR (m^3)
V = 100
# Density of A-B Mixture (kg/m^3)
rho = 1000
# Heat capacity of A-B Mixture (J/kg-K)
Cp = 0.239
# Heat of reaction for A->B (J/mol)
mdelH = 5e4
# E - Activation energy in the Arrhenius Equation (J/mol)
# R - Universal Gas Constant = 8.31451 J/mol-K
EoverR = 8750
# Pre-exponential factor (1/sec)
k0 = 7.2e10   
# U - Overall Heat Transfer Coefficient (W/m^2-K)
# A - Area - this value is specific for the U calculation (m^2)
UA = 5e4

return [
    dT_t - q/V*(Tf - T) - mdelH/(rho*Cp)*k0*tf.math.exp(-EoverR/T)*Ca - UA/V/rho/Cp*(Tc-T),
    dCa_t - q/V*(Caf - Ca) + k0*tf.math.exp(-EoverR/T)*Ca]

def func(x, y): # odes used in scipy

Temperature in CSTR (K)

#Concentration of A in CSTR (mol/m^3)
T, Ca = y

Tc = 280  # nominal = 300
# Volumetric Flowrate (m^3/sec)
q = 100   # nominal = 100
# Feed Concentration (mol/m^3)
Caf = 1 # nominal = 1
# Feed Temperature (K)
Tf = 350  # nominal = 350

#Ca_ss = 0.87725294608097;
#T_ss = 324.475443431599;

# Volume of CSTR (m^3)
V = 100
# Density of A-B Mixture (kg/m^3)
rho = 1000
# Heat capacity of A-B Mixture (J/kg-K)
Cp = 0.239
# Heat of reaction for A->B (J/mol)
mdelH = 5e4
# E - Activation energy in the Arrhenius Equation (J/mol)
# R - Universal Gas Constant = 8.31451 J/mol-K
EoverR = 8750
# Pre-exponential factor (1/sec)
k0 = 7.2e10   
# U - Overall Heat Transfer Coefficient (W/m^2-K)
# A - Area - this value is specific for the U calculation (m^2)
UA = 5e4

dT_t = q/V*(Tf - T) + mdelH/(rho*Cp)*k0*np.exp(-EoverR/T)*Ca + UA/V/rho/Cp*(Tc-T),
dCa_t = q/V*(Caf - Ca) - k0*np.exp(-EoverR/T)*Ca
return np.array(dT_t), np.array(dCa_t)

t = np.linspace(0, 3, 1000) sol = integrate.solve_ivp(func, (0, 10), (324.475443431599, 0.87725294608097), t_eval=t) T_true, Ca_true = sol.y T_true = T_true.reshape(1000, 1) Ca_true = Ca_true.reshape(1000, 1) plt.plot(t, T_true, color="black", label="T_true") plt.show() plt.plot(t, Ca_true, color="blue", label="Catrue") plt.show() t = t.reshape(1000, 1) def boundary(, on_initial): return on_initial

observe_y0 = dde.PointSetBC(t, T_true, component=0) observe_y1 = dde.PointSetBC(t, Ca_true, component=1) geom = dde.geometry.TimeDomain(0.0, 3.0) ic1 = dde.IC(geom, lambda X: 324.475443431599, boundary, component=0) ic2 = dde.IC(geom, lambda X: 0.87725294608097, boundary, component=1) data = dde.data.PDE(geom, ode_system, [ic1, ic2, observe_y0, observe_y1], num_domain=1000, num_boundary=2, anchors=t)

layer_size = [1] + [64] * 4 + [2] activation = "tanh" initializer = "Glorot normal" net = dde.maps.FNN(layer_size, activation, initializer) model = dde.Model(data, net) model.compile("adam", lr=0.001) losshistory, train_state = model.train(epochs=10000) model.compile("L-BFGS") losshistory, train_state = model.train()

dde.saveplot(losshistory, train_state, issave=True, isplot=True)

After I run it, the nn still can't converge:

Step Train loss Test loss Test metric 0 [8.75e+05, 3.25e+00, 1.05e+05, 7.70e-01] [8.75e+05, 3.25e+00, 1.05e+05, 7.70e-01] []
1000 [nan, nan, nan, nan] [nan, nan, nan, nan] []
2000 [nan, nan, nan, nan] [nan, nan, nan, nan] []

It puzzles me for a long time. Can you give me more advice? Thank you!

[lululxvi]

The loss at step 0 has 1e5 and 1e7, and then SGD blows up, and gives you nan. You should either normalize the PDE, or assign a loss weight.

[jiangxingyuan]

The loss at step 0 has 1e5 and 1e7, and then SGD blows up, and gives you nan. You should either normalize the PDE, or assign a loss weight.

I also met this problem that the loss is very large. What do you mean normalize the PDE? Hope to get your reply. Thank you!

[lululxvi]

Check FAQ "Q: I failed to train the network or get the right solution, e.g., large training loss, unbalanced losses.".