Problem of rescaling the coefficients of PDEs

[Letmesleepp]

Hi, Doctor lulu. First of all, thanks for this brilliant library. Now, I’m trying to use deepxde to deal with a batch chromatography problem using LDF model: The isotherm equation:

The boundary conditions and initial conditions: ICs: BCs:

In this model, time domain is from 0 to 280s and axis coordinate domain is from 0 to 0.1m. I have tried to rescale the domain and coefficients but the result of training is still not good. Here is my code:

import deepxde as dde
import numpy as np
from deepxde.backend import tf
import matplotlib.pyplot as plt

L = 0.1  # Column length(m)
d = 0.0212  # Column diameter(m)
A = np.pi*d*d/4  # Column cross sectional area (m2)
v = 30.0/1000.0/1000/60.0/A  # Velocity (m/s)
e = 0.62  # porosity
Ha = 3.23  # Henry constant of A (-)
Hb = 6.5  # Henry constant of B (-)
f = 0.12  # Feed concentration (-)
ba = 0.21  # vol%
bb = 0.422  # vol%
te = 280  # final time (s)
t_inj = 4  # injected time (s)

scale_x = 10
scale_t = 1/280
L_scaled = L * scale_x
v_scaled = v * scale_x / scale_t
t_scaled = te * scale_t
t_inj_scaled = t_inj * scale_t

"""LDF model"""
DLa = 1.8e-7  # diffusion coefficient
DLb = 2.2e-7  # diffusion coefficient
ka = 6.4  # Mass transfer coefficient of A (1/s)
kb = 3.8  # Mass transfer coefficient of B (1/s)
DLa_scaled = DLa * scale_x**2 / scale_t
DLb_scaled = DLb * scale_x**2 / scale_t
ka_scaled = ka / scale_t
kb_scaled = kb / scale_t

def pde(x, y):
    ca, cb, qa, qb = y[:, 0:1], y[:, 1:2], y[:, 2:3], y[:, 3:]
    ca_x = dde.grad.jacobian(y, x, i=0, j=0)
    ca_t = dde.grad.jacobian(y, x, i=0, j=1)
    ca_xx = dde.grad.hessian(y, x, component=0, i=0, j=0)

    cb_x = dde.grad.jacobian(y, x, i=1, j=0)
    cb_t = dde.grad.jacobian(y, x, i=1, j=1)
    cb_xx = dde.grad.hessian(y, x, component=1, i=0, j=0)

    qa_t = dde.grad.jacobian(y, x, i=2, j=1)
    qb_t = dde.grad.jacobian(y, x, i=3, j=1)

    massbalance_liquid_a = (
            ca_t * scale_t + (1-e)/e * qa_t * scale_t - (DLa_scaled * ca_xx) * scale_x**2 + v_scaled * ca_x * scale_x
    )
    massbalance_liquid_b = (
            cb_t * scale_t + (1-e)/e * qb_t * scale_t - (DLb_scaled * cb_xx) * scale_x**2 + v_scaled * cb_x * scale_x
    )
    massbalance_solid_a = (
            ka_scaled * ((Ha * ca) / (1 + ba * ca + bb * cb) - qa) - qa_t * scale_t
    )
    massbalance_solid_b = (
            kb_scaled * ((Hb * cb) / (1 + ba * ca + bb * cb) - qb) - qb_t * scale_t
    )
    return [massbalance_liquid_a, massbalance_liquid_b, massbalance_solid_a, massbalance_solid_b]

def feed_concentration(x):
    return np.piecewise(x[:, 1:], [x[:, 1:] <= t_inj_scaled, x[:, 1:] > t_inj_scaled], [lambda x: f, lambda x: 0])

def boundary_beg(x, on_boundary):
    return on_boundary and np.isclose(x[0], 0)

def boundary_end(x, on_boundary):
    return on_boundary and np.isclose(x[0], L_scaled)

geom = dde.geometry.Interval(0, L_scaled)
timedomain = dde.geometry.TimeDomain(0, t_scaled)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

bc_beg_ca = dde.DirichletBC(
    geomtime, feed_concentration, boundary_beg, component=0
)
bc_beg_cb = dde.DirichletBC(
    geomtime, feed_concentration, boundary_beg, component=1
)
bc_end_ca = dde.NeumannBC(
    geomtime, lambda x: 0, boundary_end, component=0
)
bc_end_cb = dde.NeumannBC(
    geomtime, lambda x: 0, boundary_end, component=1
)

initial_condition_ca = dde.IC(
    geomtime, lambda x: 0, lambda _, on_initial: on_initial, component=0
)
initial_condition_cb = dde.IC(
    geomtime, lambda x: 0, lambda _, on_initial: on_initial, component=1
)
initial_condition_qa = dde.IC(
    geomtime, lambda x: 0, lambda _, on_initial: on_initial, component=2
)
initial_condition_qb = dde.IC(
    geomtime, lambda x: 0, lambda _, on_initial: on_initial, component=3
)

"""Choose the Training points uniformly from the domain"""
z = np.linspace(0, L_scaled, num=30)
t = np.linspace(0, t_scaled, num=30)
Z, T = np.meshgrid(z, t)
Training_Points = np.vstack((Z.flatten(), T.flatten())).T

data = dde.data.TimePDE(
    geomtime,
    pde,
    [initial_condition_ca,
     initial_condition_cb,
     initial_condition_qa,
     initial_condition_qb,
     bc_beg_ca,
     bc_beg_cb,
     bc_end_ca,
     bc_end_cb],
    num_domain=0,
    num_initial=1000,
    num_boundary=1000,
    num_test=1200,
    anchors=Training_Points
)

layer_size = [2] + [50] * 6 + [4]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)

model = dde.Model(data, net)
model.compile("adam", lr=0.001, loss_weights=[1, 1, 1e-5, 1e-5, 1000, 1000, 100, 100, 100, 100, 100, 100])
losshistory, train_state = model.train(epochs=50000)

"""Get the domain: x = L_scaled and t from 0 to t_scaled"""
X_nn = L_scaled * np.ones((100, 1))
T_nn = np.linspace(0, t_scaled, 100).reshape(100, 1)
X_pred = np.append(X_nn, T_nn, axis=1)

y_pred = model.predict(X_pred)
ca_pred, cb_pred, qa_pred, qb_pred = y_pred[:, 0:1], y_pred[:, 1:2], y_pred[:, 2:3], y_pred[:, 3:]
plt.figure()
plt.plot(T_nn / scale_t, ca_pred, color='blue', linewidth=3., label='Concentration A')
plt.plot(T_nn / scale_t, cb_pred, color='red', linewidth=3., label='Concentration B')

plt.legend()
plt.grid(True)
plt.xlabel('t')
plt.ylabel('Concentration')

plt.show()

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

And the training results: 50000 [6.83e-06, 5.00e-06, 4.81e-06, 6.09e-06, 1.16e-06, 1.47e-06, 2.13e-07, 6.40e-07, 3.02e-04, 3.12e-04, 2.06e-06, 2.16e-06] [6.00e-06, 3.19e-06, 2.45e-06, 3.46e-06, 1.16e-06, 1.47e-06, 2.13e-07, 6.40e-07, 3.02e-04, 3.12e-04, 2.06e-06, 2.16e-06] []

Best model at step 47000: train loss: 5.91e-04 test loss: 5.79e-04 test metric: []

'train' took 8946.198636 s Here is the accurate one:

I also tried training 100000 steps, but the result become worse. Can you give me some advices? Thank you.

[lululxvi]

Could you first try t from 0 to 150? Also, rescale y

[Letmesleepp]

I have tried to scale the output of the neural network using the command line:

scale_y = 1e-3
net.apply_output_transform(lambda x, y: y * scale_y)

Because the solution is of order 1e-3, I multiply the network output by 1e-3.

But there's other problems.

• The boundary condition when axis coordinate x=0 let concentration of component A and B as 0.12. I wonder whether I should multiply this number by 1e-3 or not.
• Should I rescale concentration in PDEs?

I hope you can give me more details about how to rescale y. Thank you for your help.

[lululxvi]

• A simpler way: Let y = 1e-3 y. Change of variables from y to y, and solve for y.
• Yes, you should. Make everything O(1)

You may see FAQ for more examples/suggestions.

[MINE0126]

hello, @Letmesleepp I also encountered a similar problem: the parameters of the models are very different, more than a hundred times different, resulting in very suboptimal training results. Did you solve the problem? Can you give some relevant adjustment suggestions? Thank you very much！

[Letmesleepp]

@MINE0126 I followed issue#345 to scale all of my parameters. You can ask the author for more details about the method. By the way, there's no need to make all of your parameters O(1).

[MINE0126]

@Letmesleepp Thank you for your reply.