I failed to try other models with glycolysis code

[chenyv118]

I am trying to solve the ODE equation based on a model that modifies the glycolysis reaction, but the loss function keeps reporting nan and I can't get the result, even if I set the loss_weights of the ODE to 0, it doesn't help, the first 15 of the 17 state variables have a loss of nan. I hope someone can help me find the cause of the problem if possible.

The ODE equation of the system is as follows.

The code for data generation is as follows： Biodiesel_ODE.py

import numpy as np
from scipy.integrate import odeint
from data_config import *
from k_config import *

K = [Vreactf, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k_1, k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, k_10, Fa, CoAlco, af, Ae]

def enzymatic_biodiesel_model(t, k=None):
    if k is None:
        k = K
    CoAlco = k[22]
    af = k[23]
    Ae = k[24]
    Vreactf = k[0]

    def ode_func(x, t):
        T, D, M, BD, FFA, G, W, CH, E, EX, ET, ED, EM, ECH, Ef, Vp, V = x
        if t <= 120:
            Fa = k[21][0]
        else:
            Fa = k[21][1]
        if x[16] >= Vreactf:
            Fa = 0

        aT = af * (Vp / V)
        Af = (aT - Ae * (E + EX + ET + ED + EM + ECH)) / Ae

        r1 = k[1] * Ef * Af - k[11] * E
        r2 = k[2] * T * E - k[12] * ET
        r3 = k[3] * ET - k[13] * EX * D
        r4 = k[4] * D * E - k[14] * ED
        r5 = k[5] * ED - k[15] * EX * M
        r6 = k[6] * M * E - k[16] * EM
        r7 = k[7] * EM - k[17] * EX * G
        r8 = k[8] * EX * W - k[18] * FFA * E
        r9 = k[9] * EX * CH - k[19] * BD * E
        r10 = k[10] * CH * E - k[20] * ECH

        rg = (r7 * V * 92 / 1261)
        rw = (-r8 * V * 18 / 1000)

        return [
            -T * Fa / V - r2,
            -D * Fa / V + (r3 - r4),
            -M * Fa / V + (r5 - r6),
            -BD * Fa / V + r9,
            -FFA * Fa / V + r8,
            -G * Fa / V + r7,
            -W * Fa / V - r8,
            (CoAlco - CH) * Fa / V - (r9 + r10),
            -E * Fa / V + (r1 - r2 - r4 - r6 + r8 + r9 - r10),
            -EX * Fa / V + (r3 + r5 + r7 - r8 - r9),
            -ET * Fa / V + (r2 - r3),
            -ED * Fa / V + (r4 - r5),
            -EM * Fa / V + (r6 - r7),
            -ECH * Fa / V + r10,
            -Ef * Fa / V - r1,
            rg + rw,
            Fa
        ]

    # Experimental data in [mole/L]
    st = 0
    tmp = data[32][st:]
    tmp[0] = 1e-4

    FAME = tmp
    FFA = data[33][st:]
    TAG = data[34][st:]
    DAG = data[35][st:]
    MAG = data[36][st:]

    # Initial Condition
    T = TAG[0]
    D = DAG[0]
    M = MAG[0]
    B = FAME[0]
    FA = FFA[0]
    G = 1e-6
    W = mH2O
    CH = Alco
    E = 0.0
    EX = 0.0
    ET = 0.0
    ED = 0.0
    EM = 0.0
    ECH = 0.0
    Ef = Enzyme
    Vp = Vpo
    V = Vo
    x0 = [T, D, M, B, FA, G, W, CH, E, EX, ET, ED, EM, ECH, Ef, Vp, V]

    return odeint(ode_func, x0, t)

def main():
    t = np.arange(0, 1500, 1)[:, None]
    np.ravel(t)

    y, info = enzymatic_biodiesel_model(np.ravel(t))
    # y = enzymatic_biodiesel_model(np.ravel(t))
    np.savetxt("enzymatic_biodiesel.dat", np.hstack((t, y)))

if __name__ == "__main__":
    main()

The neural network code is as follows： Biodisel_NN.py

import torch
import numpy as np
import deepxde as dde
from Biodiesel_ODE import enzymatic_biodiesel_model
from data_config import *

def enzymatic_biodiesel_sbinn(data_t, data_y, noise):

    k1_ = dde.Variable(0.0)
    k2_ = dde.Variable(0.0)
    k3_ = dde.Variable(0.0)
    k4_ = dde.Variable(0.0)
    k5_ = dde.Variable(0.0)
    k6_ = dde.Variable(0.0)
    k7_ = dde.Variable(0.0)
    k8_ = dde.Variable(0.0)
    k9_ = dde.Variable(0.0)
    k10_ = dde.Variable(0.0)
    k_1_ = dde.Variable(0.0)
    k_2_ = dde.Variable(0.0)
    k_3_ = dde.Variable(0.0)
    k_4_ = dde.Variable(0.0)
    k_5_ = dde.Variable(0.0)
    k_6_ = dde.Variable(0.0)
    k_7_ = dde.Variable(0.0)
    k_8_ = dde.Variable(0.0)
    k_9_ = dde.Variable(0.0)
    k_10_ = dde.Variable(0.0)
    var_list = [k1_, k2_, k3_, k4_, k5_, k6_, k7_, k8_, k9_, k10_, k_1_, k_2_, k_3_, k_4_, k_5_, k_6_, k_7_, k_8_, k_9_, k_10_]

    def ODE(t, y):
        row, col = y.shape
        T, D, M, BD, FFA, G, W, CH, E, EX, ET, ED, EM, ECH, Ef, Vp, V = [y[:, i:i + 1] for i in range(col)]

        Fa_ = torch.tensor(Fa)[(t > 120).type(torch.long)]
        Fa_ = Fa_ * ((V < Vreactf).type(torch.long).reshape(len(V), -1))

        k1 = torch.nn.functional.softplus(k1_)
        k2 = torch.nn.functional.softplus(k2_)
        k3 = torch.nn.functional.softplus(k3_)
        k4 = torch.nn.functional.softplus(k4_)
        k5 = torch.nn.functional.softplus(k5_)
        k6 = torch.nn.functional.softplus(k6_)
        k7 = torch.nn.functional.softplus(k7_)
        k8 = torch.nn.functional.softplus(k8_)
        k9 = torch.nn.functional.softplus(k9_)
        k10 = torch.nn.functional.softplus(k10_)
        k_1 = torch.nn.functional.softplus(k_1_)
        k_2 = torch.nn.functional.softplus(k_2_)
        k_3 = torch.nn.functional.softplus(k_3_)
        k_4 = torch.nn.functional.softplus(k_4_)
        k_5 = torch.nn.functional.softplus(k_5_)
        k_6 = torch.nn.functional.softplus(k_6_)
        k_7 = torch.nn.functional.softplus(k_7_)
        k_8 = torch.nn.functional.softplus(k_8_)
        k_9 = torch.nn.functional.softplus(k_9_)
        k_10 = torch.nn.functional.softplus(k_10_)

        aT = af * (Vp / V)
        Af = (aT - Ae * (E + EX + ET + ED + EM + ECH)) / Ae

        r1 = k1 * Ef * Af - k_1 * E
        r2 = k2 * T * E - k_2 * ET
        r3 = k3 * ET - k_3 * EX * D
        r4 = k4 * D * E - k_4 * ED
        r5 = k5 * ED - k_5 * EX * M
        r6 = k6 * M * E - k_6 * EM
        r7 = k7 * EM - k_7 * EX * G
        r8 = k8 * EX * W - k_8 * FFA * E
        r9 = k9 * EX * CH - k_9 * BD * E
        r10 = k10 * CH * E - k_10 * ECH

        rg = (r7 * V * 92 / 1261)
        rw = (-r8 * V * 18 / 1000)

        return [
            dde.grad.jacobian(y, t, i=0, j=0) - (-T * Fa_ / V - r2),
            dde.grad.jacobian(y, t, i=1, j=0) - (-D * Fa_ / V + (r3 - r4)),
            dde.grad.jacobian(y, t, i=2, j=0) - (-M * Fa_ / V + (r5 - r6)),
            dde.grad.jacobian(y, t, i=3, j=0) - (-BD * Fa_ / V + r9),
            dde.grad.jacobian(y, t, i=4, j=0) - (-FFA * Fa_ / V + r8),
            dde.grad.jacobian(y, t, i=5, j=0) - (-G * Fa_ / V + r7),
            dde.grad.jacobian(y, t, i=6, j=0) - (-W * Fa_ / V - r8),
            dde.grad.jacobian(y, t, i=7, j=0) - ((CoAlco - CH) * Fa_ / V - (r9 + r10)),
            dde.grad.jacobian(y, t, i=8, j=0) - (-E * Fa_ / V + (r1 - r2 - r4 - r6 + r8 + r9 - r10)),
            dde.grad.jacobian(y, t, i=9, j=0) - (-EX * Fa_ / V + (r3 + r5 + r7 - r8 - r9)),
            dde.grad.jacobian(y, t, i=10, j=0) - (-ET * Fa_ / V + (r2 - r3)),
            dde.grad.jacobian(y, t, i=11, j=0) - (-ED * Fa_ / V + (r4 - r5)),
            dde.grad.jacobian(y, t, i=12, j=0) - (-EM * Fa_ / V + (r6 - r7)),
            dde.grad.jacobian(y, t, i=13, j=0) - (-ECH * Fa_ / V + r10),
            dde.grad.jacobian(y, t, i=14, j=0) - (-Ef * Fa_ / V - r1),
            dde.grad.jacobian(y, t, i=15, j=0) - (rg + rw),
            dde.grad.jacobian(y, t, i=16, j=0) - Fa_
        ]

    geom = dde.geometry.TimeDomain(data_t[0, 0], data_t[-1, 0])

    def boundary(x, _):  
        return np.isclose(x[0], data_t[-1, 0])  

    y1 = data_y[-1]
    bc = [dde.DirichletBC(geom, lambda X: y1[i], boundary, component=i) for i in range(0, 17)]

    n = len(data_t)
    idx = np.append(
        np.random.choice(np.arange(1, n - 1), size=n // 4, replace=False), [0, n - 1]
    )
    ic = [dde.PointSetBC(data_t[idx], data_y[idx, i:i + 1], component=0) for i in range(5)]
    np.savetxt("biodiesel_input.data", np.hstack((data_t[idx], data_y[idx, 0:1], data_y[idx, 1:2], data_y[idx, 2:3], data_y[idx, 3:4], data_y[idx, 4:5])))

    data = dde.data.PDE(geom, ODE, bc + ic, anchors=data_t)

    net = dde.maps.FNN([1] + [17] + [128] * 3 + [17], "relu", "Glorot normal")

    # def feature_transform(t):
    #     return torch.concat(
    #         (
    #             t,
    #             torch.exp(-t),
    #             torch.sin(2 * t),
    #             torch.sin(3 * t),
    #             torch.sin(4 * t),
    #             torch.sin(5 * t),
    #             torch.sin(6 * t),
    #         ),
    #         axis=1,
    #     )
    #
    # net.apply_feature_transform(feature_transform)
    # ？？？？？？TODO
    # def output_transform(t, y):
    #     return (
    #             torch.as_tensor(data_y[0]) + torch.tanh(t) * torch.tensor([1, 1, 0.1, 0.1, 0.1, 1, 0.1]) * y
    #     )
    #
    # net.apply_output_transform(output_transform)

    model = dde.Model(data, net)
    checkpointer = dde.callbacks.ModelCheckpoint(
        "./model/model.ckpt", verbose=1, save_better_only=True, period=1000
    )
    variable = dde.callbacks.VariableValue(
        var_list, period=1000, filename="variables.dat", precision=3,
    )
    callbacks = [checkpointer, variable]
    # TODO
    bc_weights = [1] * 17
    if noise >= 0.1:
        bc_weights = [w * 10 for w in bc_weights]
    # TODO
    data_weights = [1] * 5
    # Large noise requires small data_weights
    if noise >= 0.1:
        data_weights = [w / 10 for w in data_weights]
    model.compile("adam", lr=1e-3, loss_weights=[0] * 17 + bc_weights + data_weights, external_trainable_variables=var_list)
    model.train(iterations=1000, display_every=1000)
    # TODO
    ode_weights = [1] * 17
    # Large noise requires large ode_weights
    if noise > 0:
        ode_weights = [10 * w for w in ode_weights]
    model.compile("adam", lr=1e-3, loss_weights=ode_weights + bc_weights + data_weights)

    losshistory, train_state = model.train(
        epochs=900000 if noise == 0 else 2000000,
        display_every=1000,
        callbacks=callbacks,
        disregard_previous_best=True,
        # model_restore_path="./model/model.ckpt-"
    )
    dde.saveplot(losshistory, train_state, issave=True, isplot=True)
    var_list = [model.sess.run(v) for v in var_list]
    return var_list

def main():
    data_set = np.loadtxt("enzymatic_biodiesel.dat")
    t = data_set[:, 0].reshape(1500, -1)
    y = data_set[:, 1:]
    noise = 0
    # Add noise
    # if noise > 0:
    #     std = noise * y.std(0)
    #     y[1:-1, :] += np.random.normal(0, std, (y.shape[0] - 2, y.shape[1]))
    #     np.savetxt("enzymatic_biodiesel_noise.dat", np.hstack((t, y)))

    # Train
    var_list = enzymatic_biodiesel_sbinn(t, y, noise)

    # Prediction
    y = enzymatic_biodiesel_model(np.ravel(t), *var_list)
    np.savetxt("glycolysis_pred.dat", np.hstack((t, y)))

if __name__ == "__main__":
    main()

I can sort out the parameters in k_config and data_config later if needed. One thing to note is that Fa will change after 120 minutes and will equal 0 when Fa>Vreactf, and this mutation may have some impact. I'm not sure.

[lululxvi]

Please see new code at https://github.com/lu-group/sbinn

[chenyv118]

Please see new code at https://github.com/lu-group/sbinn

Yes that's right, my code didn't have a restriction on the division operation of state variables, using this new code solved the problem, thanks a lot!