用 Scipy 求解 ODE 集

Question

我正在尝试开发一种算法（使用 scipy.integrate.odeint()）预测细胞、底物和产物浓度的变化（即 , , ）随着时间的推移，直到系统达到稳定状态（~100 或 200 小时） .生物反应器中细胞的初始浓度为0.1/，反应器中最初没有葡萄糖或产物。我想针对介于 0.01 /ℎ 和 0.25 /ℎ 之间的一系列不同流量测试算法，并分析流量对产品生产的影响（即 ⋅ in /ℎ）。最终，我想在 x-axis 上生成一个显示产品生产率 (y-axis) 与流速的关系图。我的目标是估计导致最大（或临界）生产率的流速。到目前为止，这是我的代码：

from scipy.integrate import odeint
import numpy as np

# Constants
u_max = 0.65
K_s = 0.14
K_1 = 0.48
V = 2
X_in = 0
S_in = 4
Y_s = 0.38
Y_p = 0.2

# Variables
# Q - Flow Rate (L/h), value between 0.01 and 0.25 that produces best Q * P
# X - Cell Concentration (g/L)
# S - The glucose concentration (g/L)
# P - Product Concentration (g/L)

# Equations
def func_dX_dt(X, t, S):
    u = (u_max) / (1 + (K_s / S))
    dX_dt = (((Q * S_in) - (Q * S)) / V) + (u * X)
    return dX_dt

def func_dS_dt(S, t, X):
    u = (u_max) / (1 + (K_s / S))
    dS_dt = (((Q * S_in) - (Q * S)) / V) - (u * (X / Y_s))
    return dS_dt

def func_dP_dt(P, t, X, S):
    u = (u_max) / (1 + (K_s / S))
    dP_dt = ((-Q * P) / V) - (u * (X / Y_p))
    return dP_dt

t = np.linspace(0, 200, 200)

# Q placeholder
Q = 0.01

# Attempt to solve the Ordinary differential equations
sol_dX_dt = odeint(func_dX_dt, 0.1, t, args=(S,))
sol_dS_dt = odeint(func_dS_dt, 0.1, t, args=(X,))
sol_dP_dt = odeint(func_dP_dt, 0.1, t, args=(X,S))

在程序当前状态下，似乎没有办法为 P 生成稳态值。我试图进行此修改以获得 X 的值。

sol_dX_dt = odeint(func_dX_dt, 0.1, t, args=(odeint(func_dS_dt, 0.1, t, args=(X,)),))

它产生错误：

NameError: name 'X' is not defined

此时我不确定如何前进。

（编辑 1：添加原始方程式）

第一个方程

第二个方程和第三个方程

Answer 1

您不必将函数应用于每个部分，而是 return 导数的元组，如下所示：

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

Q = 0.01
V = 2
Ys = 0.38
Sin = 4
Yp = 0.2
Xin = 0
umax = 0.65
Ks = 0.14
K1 = 0.48

def mu(S, umax, Ks, K1):
    return umax/((1+Ks/S)*(1+S/K1))

def dxdt(x, t, *args):
    X, S, P = x
    Q, V, Xin, Ys, Sin, Yp, umax, Ks, K1 = args
    m = mu(S, umax, Ks, K1)
    dXdt = (Q*Xin - Q*X)/V + m*X
    dSdt = (Q*Sin - Q*S)/V - m*X/Ys
    dPdt = -Q*P/V - m*X/Yp
    return dXdt, dSdt, dPdt

t = np.linspace(0, 200, 200)
X0 = 0.1
S0 = 0.1
P0 = 0.1
x0 = X0, S0, P0
sol = odeint(dxdt, x0, t, args=(Q, V, Xin, Ys, Sin, Yp, umax, Ks, K1))

plt.plot(t, sol[:, 0], 'r', label='X(t)')
plt.plot(t, sol[:, 1], 'g', label='S(t)')
plt.plot(t, sol[:, 2], 'b', label='P(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()

输出：

用 Scipy 求解 ODE 集

Solving set of ODEs with Scipy

python

scipy

ode

differential-equations