fitting a variables in a system of ODE as a function of time

[pabloacera]

Hi,

Thanks for developing this awesome package!
I want to solve an inverse problem modelled as a system of ODE evolving through time. I want to calculate several parameters that also depend on time. My problem is similar to the Lorenz system example but some of the parameters I want to estimate (𝜎, 𝜌, and 𝛽) are a function of time, so instead of one value they will have a different value per time point, and I have no idea the shape or nature of that function.
Please, Is there a way to fo this with your package? if so, what would be the best way of doing this?
Thanks a lot for your time!
Cheers,
Pablo.

[haison19952013]

How about considering (𝜎, 𝜌, and 𝛽) as additional outputs of your network? Therefore, they will be a function of time...Plz refer to this example to get some ideas: https://deepxde.readthedocs.io/en/latest/demos/pinn_inverse/elliptic.inverse.field.html.

[pabloacera]

Hi, that is a good idea.
I have tried to do that. I tried to train the networks to estimate 'a' as an output of the network and 'b' and 'y' as static parameters. The code runs and the loss of the ODEs is very very small but all the parameters are way off compared to the ones I used to simulate the data. Here is a code snippet that can be reproducible.

# Generate data 

# function that returns dy/dt
def model(z,t,a):
    u = z[0]
    s = z[1]
    # unsplice
    b = 1 # splicing rate
    Y = 0.8 # degradation rate
    dudt = a - b*u
    dsdt = b*u - Y*s
    dzdt = [dudt,dsdt]
    return dzdt

# initial condition
z0 = [0,0]

# number of time points
n = 100

# time points
time = np.linspace(0,100,n)

a = np.zeros(n)

# change to 2.0 at time = 5.0
a[25:50] = 2.0
a[50:100] = 0

u = np.empty_like(time)
s = np.empty_like(time)

# record initial conditions
u[0] = z0[0]
s[0] = z0[1]


# solve ODE
for i in range(1,n):
    # span for next time step
    tspan = [time[i-1],time[i]]
    # solve for next step
    z = odeint(model,z0,tspan,args=(a[i],))
    # store solution for plotting
    u[i] = z[1][0]
    s[i] = z[1][1]
    # next initial condition
    z0 = z[1]



# define the static variables
b = dde.Variable(1.0)
y = dde.Variable(1.0)

def RNA_dynamics(in_, out_):
    '''
    in_ -> network input
    out_ -> network output
    a is also part of the output of the function
    '''
    u, a, s = out_[:, 0:1], out_[:, 1:2], out_[:, 2:]
    dy1_x = dde.grad.jacobian(out_, in_, i=0)
    dy2_x = dde.grad.jacobian(out_, in_, i=2)
    return [
        dy1_x - a - b*u,
        dy2_x - b*u - y*s
        ]

def boundary(_, on_initial):
    return on_initial

maxtime = 3
geom = dde.geometry.TimeDomain(0, maxtime)

# Training data
u[0] = 0
s[0] = 0

# Initial conditions
ic1 = dde.icbc.IC(geom, lambda X: 0, boundary, component=0)
ic2 = dde.icbc.IC(geom, lambda X: 0, boundary, component=2)

# Get the training data
observe_t, ob_y = time.reshape(-1, 1), np.vstack((u, s)).T
observe_u = dde.icbc.PointSetBC(observe_t, ob_y[:, 0:1], component=0)
observe_s = dde.icbc.PointSetBC(observe_t, ob_y[:, 1:], component=2)

# define data object
data = dde.data.PDE(
                    geom,
                    RNA_dynamics,
                    [ic1, ic2, observe_u, observe_s],
                    num_domain=100,
                    num_boundary=2,
                    anchors=observe_t)


net = dde.nn.PFNN([1, [20, 20, 20], [20, 20, 20], [20, 20, 20], 3], "tanh", "Glorot uniform")

model = dde.Model(data, net)
model.compile("adam", lr=0.0001, external_trainable_variables=[b, y])

# callbacks for storing results
fnamevar = "variables.dat"
variable = dde.callbacks.VariableValue([b, y], period=100, filename=fnamevar)

losshistory, train_state = model.train(iterations=20000, callbacks=[variable])

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

# view results
yhat = model.predict(observe_t)
uhat, ahat, shat = yhat[:, 0:1], yhat[:, 1:2], yhat[:, 2:]

# Plots
utrue = u
print("l2 relative error for u: " + str(dde.metrics.l2_relative_error(utrue, uhat)))
plt.figure()
plt.plot(observe_t, utrue, "-", label="u_true")
plt.plot(observe_t, uhat, "--", label="u_NN")
plt.legend()

strue = s
print("l2 relative error for q: " + str(dde.metrics.l2_relative_error(strue, shat)))
plt.figure()
plt.plot(observe_t, strue, "-", label="s_true")
plt.plot(observe_t, shat, "--", label="s_NN")
plt.legend()

atrue = a
plt.figure()
plt.plot(observe_t, atrue, "-", label="a_true")
plt.plot(observe_t, ahat, "--", label="a_NN")
plt.legend()

plt.show()


# reopen saved data using callbacks in fnamevar
lines = open(fnamevar, "r").readlines()
# read output data in fnamevar (this line is a long story...)
Chat = np.array(
    [
        np.fromstring(
            min(re.findall(re.escape("[") + "(.*?)" + re.escape("]"), line), key=len),
            sep=",",
        )
        for line in lines
    ]
)

b_true = 1 # splicing rate
y_true = 0.8 # degradation rate

l, c = Chat.shape
plt.plot(range(l), Chat[:, 0], "r-")
plt.plot(range(l), Chat[:, 1], "g-")
plt.plot(range(l), np.ones(Chat[:, 0].shape) * b_true, "r--")
plt.plot(range(l), np.ones(Chat[:, 1].shape) * y_true, "k--")
plt.legend(["b_hat", "y_hat", "True b", "True y"], loc="right")
plt.xlabel("Epoch")

Do you think it could be an error in the way I am defining the ODE or the network? In case the solution is valid, can I put an extra constraint in the parameters? for example, all of them both static and dynamic values should be >0? Any insights from you would be very helpful. Thank you so much for your time in advance!
Pablo.

[lululxvi]

You can check manually if the found solution is current or not. If it is, then yes, you can add extra constraints.