Network not learning the ODE in inverse problem

[isaingit]

Dear all,

I am trying to solve an inverse problem with a singe external variable. I do not have any data for this external variable but I would like it to be as small as possible.

When I train the model, I reaches a plateau in the loss function after some iterations.

Inputs and outputs have already been scaled and their ranges are between 0 and 1.

Could anyone help on this?

import numpy as np
import deepxde as dde
import tensorflow as tf
import matplotlib.pyplot as plt
import pandas as pd
import pickle

# Set random seed
seed = 0
np.random.seed(seed)
tf.random.set_seed(seed)
dde.backend.tf.random.set_random_seed(seed)

save_ckpts_path = r'C:\Users\Usuario\Desktop\PinnCkpts'

#region DATA
# Model known parameters
diameter = 1.5 # m
S = np.pi * diameter**2 / 4 # m^2
molFlow0 = 324.1 #kmol/h
molFrac0 = [0.341, 0.659, 0.] # propylene , benzene, cumene

R = 8.314 # kJ/kmol/K = m3·Pa/K/mol
E = 56000. # kJ/kmol
A = 3500 # m3/kmol/s
A = A * 3600 # m3/kmol/h
conversion = 0.891835592484582 # 90% of Gibbs equilibrium conversion
stochCoef = [-1., -1., 1.] # stoichiometric coefficients of the reaction

# Initial conditions
n0 = np.multiply(molFlow0,molFrac0) # kmol/h
T0 = 350. # C
P0 = 3075. # kPa
P0 = P0 * 1000 # Pa

# Boudnaries of the normalized time interval
lmin, lmax = (0.,1.)

# Training points
n_domain = 1000
n_boundary = 2

# NN hyperparameters
n_outputs = 1 # number of outputs of the NN
n_balances = 1 # number of material balances
n_optVar = 1 # number of optimization variables
n_bcs = 2
Lguess = 10. # initial guess for the length of the reactor

activation_fcn = "tanh"
epochs_adam = 2000
lr = .001
batch_size = 100 # input attribute to generate the resampler callback
n_hidden_layers = 2
n_neurons_per_layer = 64

w_balances = 1.
w_optVar = 1e-5
w_bcs = 1.
weights_array = [w_balances]*n_balances + [w_optVar]*n_optVar + [w_bcs]*n_bcs
#endregion DATA

#STEP 1: IDENTIFY THE UNKNOWN PARAMETERS
L = dde.Variable(Lguess) # Length of the reactor. 

#STEP 2: CREATE TIME DOMAIN
geom = dde.geometry.TimeDomain(0.0, 1.0) # normalized length. Time domain because length is to PFR as time is to Batch reactors

#STEP 3: DEFINE THE ODE
def ode(l, u):
    X = u[:,0:1]
    T = T0 # ºC
    P = P0 # Pa

    dX_dl = dde.grad.jacobian(u, l, i=0)

    n = n0 + np.multiply(stochCoef,n0[0]) * X
    
    ctot = (P/R/(T+273.15)) # mol/m3
    ctot /= 1000 # kmol/m3

    k = A * tf.exp(-E/(R*(T + 273.15))) # m3/kmol/h

    odeX = dX_dl - L * S  * k * (ctot**2) * (1 - X) * n[1] / (tf.reduce_sum(n)**2) # ODE for the mole balance

    # ODE TO ENSURE MINIMUM LENGTH
    ode2 = L * tf.ones_like(X)

    return [odeX, ode2]

#STEP 4: DEFINE BOUNDARY CONDITIONS
inX  = dde.icbc.DirichletBC(geom, lambda x: 0., lambda l, on_boundary: on_boundary * dde.utils.isclose(l[0],lmin), component=0)
endX = dde.icbc.DirichletBC(geom, lambda x: conversion, lambda l, on_boundary: on_boundary * dde.utils.isclose(l[0],lmax), component=0)

#STEP 5: CREATE THE DATA OBJECT (IT CONTAINS THE DISCRETIZATION OF THE PROBLEM ACROSS THE TIME DOMAIN)
data = dde.data.PDE(geom, ode, [inX,endX], num_domain=n_domain, num_boundary=n_boundary)

#STEP 6: BUILD THE NEURAL NETWORK
def transform_function(nn_inputs, nn_outputs):
    X = nn_outputs[:,0:1]

    return tf.sigmoid(X) # sigmoid rather than tanh because the output is non-negative

net = dde.nn.FNN(layer_sizes=[1]+[n_neurons_per_layer]*n_hidden_layers+[n_outputs], activation=activation_fcn, kernel_initializer="Glorot normal") # sigmoid rather than tanh because the output is non-negative
net.apply_output_transform(transform_function)

#STEP 7: COMPILE AND TRAIN THE MODEL
resampler = dde.callbacks.PDEPointResampler(period=batch_size)
variable = dde.callbacks.VariableValue(L, period=50)
prefix = save_ckpts_path +'\checkpoint'
ckpts = dde.callbacks.ModelCheckpoint(prefix, save_better_only=True) # This command only saves the best model. Here, the default value foor period = 1 iteration.

model = dde.Model(data, net)
model.compile("adam", lr=lr, loss_weights=weights_array, external_trainable_variables=L)
losshistory, train_state = model.train(display_every=10, iterations=epochs_adam, callbacks=[resampler, variable, ckpts])

# PLOT RESUTLS
# Figure 1 will display the learning curve and Figure 2, the output predictions vs inputs for the best trained result
dde.saveplot(losshistory, train_state, issave=False, isplot=True)

Thanks in advance for the attention