Zack - 2 months ago 3x
Python Question

Update initial condition in ODE solver each time step

I am wanting to solve a system of ODEs where for the first 30,000 seconds, I want one of my state variables to start from the same initial value. After those 30,000 seconds, I want to change the initial value of that state variable to something different and simulate the system for the rest of time. Here is my code:

``````def ode_rhs(y, t):
ydot[0] = -p[7]*y[0]*y[1] + p[8]*y[8] + p[9]*y[8]
ydot[1] = -p[7]*y[0]*y[1] + p[8]*y[8]
ydot[2] = -p[10]*y[2]*y[3] + p[11]*y[9] + p[12]*y[9]
ydot[3] = -p[13]*y[3]*y[6] + p[14]*y[10] + p[15]*y[10] - p[10]*y[2]*y[3] + p[11]*y[9] + p[9]*y[8] - p[21]*y[3]
ydot[4] = -p[19]*y[4]*y[5] - p[16]*y[4]*y[5] + p[17]*y[11] - p[23]*y[4] + y[7]*p[20]
ydot[5] = -p[19]*y[4]*y[5] + p[15]*y[10] - p[16]*y[4]*y[5] + p[17]*y[11] + p[18]*y[11] + p[12]*y[9] - p[22]*y[5]
ydot[6] = -p[13]*y[3]*y[6] + p[14]*y[10] - p[22]*y[6] - p[25]*y[6] - p[23]*y[6]
ydot[7] = 0
ydot[8] = p[7]*y[0]*y[1] - p[8]*y[8] - p[9]*y[8]
ydot[9] = p[10]*y[2]*y[3] - p[11]*y[9] - p[12]*y[9] - p[21]*y[9]
ydot[10] = p[13]*y[3]*y[6] - p[14]*y[10] - p[15]*y[10] - p[22]*y[10] - p[21]*y[10] - p[23]*y[10]
ydot[11] = p[19]*y[4]*y[5] + p[16]*y[4]*y[5] - p[17]*y[11] - p[18]*y[11] - p[22]*y[11] - p[23]*y[11]
ydot[12] = p[22]*y[10] + p[22]*y[11] + p[22]*y[5] + p[22]*y[6] + p[21]*y[10] + p[21]*y[3] + p[21]*y[9] + p[24]*y[13] + p[25]*y[6] + p[23]*y[10] + p[23]*y[11] + p[23]*y[4] + p[23]*y[6]
ydot[13] = p[15]*y[10] + p[18]*y[11] - p[24]*y[13]
return ydot

pysb.bng.generate_equations(model)
alias_model_components()
p = np.array([k.value for k in model.parameters])
ydot = np.zeros(len(model.odes))
y0 = np.zeros(len(model.odes))
y0[0:7] = p[0:7]
t = np.linspace(0.0,1000000.0,100000)
r = odeint(ode_rhs,y0,t)
``````

So, in other words, I want to set y0[1] to the same value (100) each time
`odeint`
is called for the first 30,000 seconds. I'm effectively trying to let the system equilibrate for an amount of time before inputing a signal into the system. I thought about doing something like
`if t < 30000: y0[1] = 100`
as the first line of my
`ode_rhs()`
function, but I'm not quite sure that works.

It sounds like you want y1(t) to remain constant (with the value 100) for the equilibration phase. You can do this by ensuring that dy1(t)/dt = 0 during this phase. There are (at least) two ways you can accomplish that. The first is to modify the calculation of `ydot[1]` in `ode_rhs` as follows:

``````if t < 30000:
ydot[1] = 0.0
else:
ydot[1] = -p[7]*y[0]*y[1] + p[8]*y[8]
``````

and use the intitial condition 100 for `y[1]`.

Note that this introduces a discontinuity in the right-hand side of the system, but the adaptive solver used by `odeint` (the Fortran code LSODA) is usually robust enough to handle it.

Here's a self-contained example. I've made `p` and `t1` arguments to `ode_rhs`. `t1` is the duration of the equilibration phase.

``````import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def ode_rhs(y, t, p, t1):
ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2]
if t < t1:
ydot[1] = 0.0
else:
ydot[1] = -p[0]*y[0]*y[1] + p[1]*y[2]
ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2]
return ydot

ydot = np.zeros(3)
p = np.array([0.01, 0.25, 0.1])
y0 = [20.0, 100.0, 0.0]
t = np.linspace(0, 200, 2001)
t1 = 20.0

sol = odeint(ode_rhs, y0, t, args=(p, t1))

plt.figure(1)
plt.clf()

plt.subplot(3, 1, 1)
plt.plot(t, sol[:, 0])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[0]')

plt.subplot(3, 1, 2)
plt.plot(t, sol[:, 1])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[1]')
plt.ylim(0, 110)

plt.subplot(3, 1, 3)
plt.plot(t, sol[:, 2])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[2]')
plt.xlabel('t')

plt.show()
``````

A slight variation of the above method is to modify the system by adding a parameter that is either 0 or 1. When the parameter is 0, the equlibration system is solved, and when the parameter is 1, the full system is solved. The code for `ydot[1]` (in my smaller example) is then

``````ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2])
``````

where `full` is the parameter.

To handle the equilibration phase, the system is solved once on 0 <= t < t1 with `full=0`. Then the final value of the equilibration solution is used as the initial condition to the second solution, run with `full=1`. The advantage of this method is that you are not forcing the solver to deal with the discontinuity.

Here's how it looks in code.

``````import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def ode_rhs(y, t, p, full):
ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2]
ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2])
ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2]
return ydot

ydot = np.zeros(3)
p = np.array([0.01, 0.25, 0.1])
y0 = [20.0, 100.0, 0.0]
t1 = 20.0  # Equilibration time
tf = 200.0  # Final time

# Solve the equilibration phase.
teq = np.linspace(0, t1, 100)
full = 0
soleq = odeint(ode_rhs, y0, teq, args=(p, full))

# Solve the full system, using the final point of the
# equilibration phase as the initial condition.
y0 = soleq[-1]
# Note: the system is autonomous, so we could just as well start
# at t0=0.  But starting at t1 makes the plots (below) align without
# any additional shifting of the time arrays.
t = np.linspace(t1, tf, 2000)
full = 1
sol = odeint(ode_rhs, y0, t, args=(p, full))

plt.figure(2)
plt.clf()
plt.subplot(3, 1, 1)
plt.plot(teq, soleq[:, 0], t, sol[:, 0])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[0]')

plt.subplot(3, 1, 2)
plt.plot(teq, soleq[:, 1], t, sol[:, 1])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[1]')
plt.ylim(0, 110)

plt.subplot(3, 1, 3)
plt.plot(teq, soleq[:, 2], t, sol[:, 2])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[2]')
plt.xlabel('t')

plt.show()
``````

And here's the plot that it generates (the plot from the first example is almost exactly the same):