补料分批生物反应器的 ODE 方程的 GEKKO 不可行系统
GEKKO Infeasible system of ODE equations of a fed-batch Bioreactor
我是 GEKKO 的新手，也是生物反应器建模的新手，所以我可能遗漏了一些明显的东西。
我有一个包含 10 个 ODE 的系统，用于描述补料分批生物反应器。给出了所有常数。下图显示了该模型的预期行为（摘自论文）。然而，我发现唯一可行的解决方案是活细胞密度 (XV) = 0，并且在所有时间 t 内都保持为 0，或者时间 T 非常小 (<20)。如果下边界 >= 0 或初始值设置为 XV 且 t > 20，则系统变得不可行。
方程式和常数被多次检查。我尝试为我的变量赋予初始值，但它也没有用。我只能想到两个问题：我没有正确启动变量，或者我没有正确使用 GEKKO。有任何想法吗？谢谢！！
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
m = GEKKO(remote=False) # create GEKKO model
#constants 3L continuous fed-batch
KdQ = 0.001 #degree of degradation of glutamine (1/h)
mG = 1.1*10**-10 #glucose maintenance coefficient (mmol/cell/hour)
YAQ = 0.90 #yield of ammonia from glutamine
YLG = 2 #yield of lactate from glucose
YXG = 2.2*10**8 #yield of cells from glucose (cells/mmol)
YXQ = 1.5*10**9 #yield of cells from glutamine (cells/mmol)
KL = 150 #lactate saturation constant (mM)
KA = 40 #ammonia saturation constant (mM)
Kdmax = 0.01 #maximum death rate (1/h)
mumax = 0.044 #maximum growth rate (1/h)
KG = 1 #glucose saturation constant (mM)
KQ = 0.22 #glutamine saturation constant (mM)
mQ = 0 #glutamine maintenance coefficient (mmol/cell/hour)
kmu = 0.01 #intrinsic death rate (1/h)
Klysis = 2*10**-2 #rate of cell lysis (1/h)
Ci_star = 100 #inhibitor saturation concentration (mM)
qi = 2.5*10**-10 #specific inhibitor production rate (1/h)
#Flow, volume and concentration
Fo = 0.001 #feed-rate (L/h)
Fi = 0.001 #feed-rate (L/h)
V = 3 #volume (L)
SG = 653 #glucose concentration in the feed (mM)
SQ = 58.8 #glutamine concentration in the feed (mM)
# create GEKKO parameter
t = np.linspace(0,120,121)
m.time = t
XT = m.Var(name='XT') #total cell density (cells/L)
XV = m.Var(lb=0, name='XV') #viable cell density (cells/L)
XD = m.Var(name='XD') #dead cell density (cells/L)
G = m.Var(value = 30, name='G') #glucose concentration (mM)
Q = m.Var(value = 4.5, name='Q') #glutamine concentration (mM)
L = m.Var(name='L') #lactate concentration (mM)
A = m.Var(name='A') #ammonia concentration (mM)
Ci = m.Var(name='Ci') #inhibitor concentration (mM)
mu = m.Var(name='mu') #growth rate (1/h)
Kd = m.Var(name='Kd') #death rate(1/h)
# create GEEKO equations
m.Equation(XT.dt() == mu*XV - Klysis*XD - XT*Fo/V)
m.Equation(XV.dt() == (mu - Kd)*XV - XV*Fo/V)
m.Equation(XD.dt() == Kd*XV - Klysis*XD - XV*Fo/V)
m.Equation(G.dt() == (Fi/V)*SG - (Fo/V)*G + (-mu/YXG - mG)*XV)
m.Equation(Q.dt() == (Fi/V)*SQ - (Fo/V)*Q + (-mu/YXQ - mQ)*XV - KdQ*Q)
m.Equation(L.dt() == -YLG*(-mu/YXG -mG)*XV-(Fo/V)*L)
m.Equation(A.dt() == -YAQ*(-mu/YXQ - mQ)*XV +KdQ*Q-(Fo/V)*A)
m.Equation(Ci.dt() == qi*XV - (Fo/V)*Ci)
m.Equation(mu.dt() == (mumax*G*Q*(Ci_star-Ci)) / (Ci_star*(KG+G)*(KQ+Q)*(L/KL + 1)*(A/KA + 1)))
m.Equation(Kd.dt() == Kdmax*(kmu/(mu+kmu)))
# solve ODE
m.options.IMODE = 4
m.open_folder()
m.solve(display = False)
plt.plot(m.time, XV.value)
使用完全相同模型的文章：
1) 使用GEKKO的硕士论文"MODELING OF MAMMALIAN CELL CULTURE"
link:
2) 描述方程式的原始论文：“工艺模型比较和跨生物反应器规模的可转移性和
哺乳动物细胞生物过程的操作模式
link: https://sci-hub.tw/10.1002/btpr.1664
3) 控制系统使用该模型的论文："Glucose concentration control of a fed-batch mammalian cell bioprocess using a nonlinear model predictive controller"
link: https://sci-hub.tw/https://doi.org/10.1016/j.jprocont.2014.02.007
归根结底这不是编程问题，而是阅读方程式并正确翻译它们的问题。
mu 和 Kd 不是动态变量，它们是状态向量（只有 8 维）的普通函数。对于这样的中间变量，Gekko 具有构造函数 m.Intermediate
mu = m.Intermediate((mumax*G*Q*(Ci_star-Ci)) / (Ci_star*(KG+G)*(KQ+Q)*(L/KL + 1)*(A/KA + 1)), name='mu') #growth rate (1/h)
Kd = m.Intermediate(Kdmax*(kmu/(mu+kmu))) #death rate(1/h)
通过此更改以及 XT 和 XV 的初始值 5e8，脚本给出了下面的图表，这些图表显示了您引用的论文中可以找到的内容。
有几个问题：
- 最后两个方程是代数方程，不是微分方程。应该是
mu==... 而不是 mu.dt()==...
- 一些方程式可能被零除。
x.dt() = z/y 等式可以用 y * x.dt()==z 代替，这样如果 y 趋近于零，则等式变为 0 * x.dt() == z。
- 一些初始条件未设置，因此它们使用默认值 0。这可能会创建一个零解。
我输入了一些不同的值并使用 m.options.COLDSTART=2 来帮助它找到初始解决方案。我还使用 Intermediates 来帮助可视化任何变大的术语。我将细胞浓度以每升数百万个细胞为单位来帮助缩放。
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
m = GEKKO(remote=False) # create GEKKO model
#constants 3L continuous fed-batch
KdQ = 0.001 #degree of degradation of glutamine (1/h)
mG = 1.1e-10 #glucose maintenance coefficient (mmol/cell/hour)
YAQ = 0.90 #yield of ammonia from glutamine
YLG = 2 #yield of lactate from glucose
YXG = 2.2e8 #yield of cells from glucose (cells/mmol)
YXQ = 1.5e9 #yield of cells from glutamine (cells/mmol)
KL = 150 #lactate saturation constant (mM)
KA = 40 #ammonia saturation constant (mM)
Kdmax = 0.01 #maximum death rate (1/h)
mumax = 0.044 #maximum growth rate (1/h)
KG = 1 #glucose saturation constant (mM)
KQ = 0.22 #glutamine saturation constant (mM)
mQ = 0 #glutamine maintenance coefficient (mmol/cell/hour)
kmu = 0.01 #intrinsic death rate (1/h)
Klysis = 2e-2 #rate of cell lysis (1/h)
Ci_star = 100 #inhibitor saturation concentration (mM)
qi = 2.5e-10 #specific inhibitor production rate (1/h)
#Flow, volume and concentration
Fo = 0.001 #feed-rate (L/h)
Fi = 0.001 #feed-rate (L/h)
V = 3 #volume (L)
SG = 653 #glucose concentration in the feed (mM)
SQ = 58.8 #glutamine concentration in the feed (mM)
# create GEKKO parameter
t = np.linspace(0,50,121)
m.time = t
XTMM = m.Var(value=1,name='XT') #total cell density (MMcells/L)
XVMM = m.Var(value=1,lb=0, name='XV') #viable cell density (MMcells/L)
XDMM = m.Var(value=1.0,name='XD') #dead cell density (MMcells/L)
G = m.Var(value = 20, name='G') #glucose concentration (mM)
Q = m.Var(value = 4.5, name='Q') #glutamine concentration (mM)
L = m.Var(value=1,name='L') #lactate concentration (mM)
A = m.Var(value=1.6,name='A') #ammonia concentration (mM)
Ci = m.Var(value=0.1,name='Ci') #inhibitor concentration (mM)
mu = m.Var(value=0.1,name='mu') #growth rate (1/h)
Kd = m.Var(value=0.5,name='Kd') #death rate(1/h)
# scale back to cells/L from million cells/L
XT = m.Intermediate(XTMM*1e7)
XV = m.Intermediate(XVMM*1e7)
XD = m.Intermediate(XDMM*1e7)
e1 = m.Intermediate((mu*XV - Klysis*XD - XT*Fo/V)/1e7)
e2 = m.Intermediate(((mu - Kd)*XV - XV*Fo/V)/1e7)
e3 = m.Intermediate((Kd*XV - Klysis*XD - XV*Fo/V)/1e7)
e4 = m.Intermediate((Fi/V)*SG - (Fo/V)*G + (-mu/YXG - mG)*XV)
e5 = m.Intermediate((Fi/V)*SQ - (Fo/V)*Q + (-mu/YXQ - mQ)*XV - KdQ*Q)
e6 = m.Intermediate(-YLG*(-mu/YXG -mG)*XV-(Fo/V)*L)
e7 = m.Intermediate(-YAQ*(-mu/YXQ - mQ)*XV +KdQ*Q-(Fo/V)*A)
e8 = m.Intermediate(qi*XV - (Fo/V)*Ci)
e9a = m.Intermediate((Ci_star*(KG+G)*(KQ+Q)*(L/KL + 1)*(A/KA + 1)))
e9b = m.Intermediate((mumax*G*Q*(Ci_star-Ci)))
e10a = m.Intermediate((mu+kmu))
e10b = m.Intermediate(Kdmax*kmu)
# create GEEKO equations
m.Equation(XTMM.dt() == e1)
m.Equation(XVMM.dt() == e2)
m.Equation(XDMM.dt() == e3)
m.Equation(G.dt() == e4)
m.Equation(Q.dt() == e5)
m.Equation(L.dt() == e6)
m.Equation(A.dt() == e7)
m.Equation(Ci.dt() == e8)
m.Equation(e9a * mu == e9b)
m.Equation(e10a*Kd == e10b)
# solve ODE
m.options.IMODE = 4
m.options.SOLVER = 1
m.options.NODES = 3
m.options.COLDSTART = 2
#m.open_folder()
m.solve(display=False)
plt.figure()
plt.subplot(3,1,1)
plt.plot(m.time, XV.value,label='XV')
plt.plot(m.time, XT.value,label='XT')
plt.plot(m.time, XD.value,label='XD')
plt.legend()
plt.subplot(3,1,2)
plt.plot(m.time, G.value,label='G')
plt.plot(m.time, Q.value,label='Q')
plt.plot(m.time, L.value,label='L')
plt.plot(m.time, A.value,label='A')
plt.legend()
plt.subplot(3,1,3)
plt.plot(m.time, mu.value,label='mu')
plt.plot(m.time, Kd.value,label='Kd')
plt.legend()
plt.xlabel('Time (hr)')
plt.figure()
plt.plot(m.time, e1.value,'r-.',label='eqn1')
plt.plot(m.time, e2.value,'g:',label='eqn2')
plt.plot(m.time, e3.value,'b:',label='eqn3')
plt.plot(m.time, e4.value,'b--',label='eqn4')
plt.plot(m.time, e5.value,'y:',label='eqn5')
plt.plot(m.time, e6.value,'m--',label='eqn6')
plt.plot(m.time, e7.value,'b-.',label='eqn7')
plt.plot(m.time, e8.value,'g--',label='eqn8')
plt.plot(m.time, e9a.value,'r:',label='eqn9a')
plt.plot(m.time, e9b.value,'r--',label='eqn9b')
plt.plot(m.time, e10a.value,'k:',label='eqn10a')
plt.plot(m.time, e10b.value,'k--',label='eqn10b')
plt.legend()
plt.show()
我是 GEKKO 的新手，也是生物反应器建模的新手，所以我可能遗漏了一些明显的东西。
我有一个包含 10 个 ODE 的系统，用于描述补料分批生物反应器。给出了所有常数。下图显示了该模型的预期行为（摘自论文）。然而，我发现唯一可行的解决方案是活细胞密度 (XV) = 0，并且在所有时间 t 内都保持为 0，或者时间 T 非常小 (<20)。如果下边界 >= 0 或初始值设置为 XV 且 t > 20，则系统变得不可行。
方程式和常数被多次检查。我尝试为我的变量赋予初始值，但它也没有用。我只能想到两个问题：我没有正确启动变量，或者我没有正确使用 GEKKO。有任何想法吗？谢谢！！
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
m = GEKKO(remote=False) # create GEKKO model
#constants 3L continuous fed-batch
KdQ = 0.001 #degree of degradation of glutamine (1/h)
mG = 1.1*10**-10 #glucose maintenance coefficient (mmol/cell/hour)
YAQ = 0.90 #yield of ammonia from glutamine
YLG = 2 #yield of lactate from glucose
YXG = 2.2*10**8 #yield of cells from glucose (cells/mmol)
YXQ = 1.5*10**9 #yield of cells from glutamine (cells/mmol)
KL = 150 #lactate saturation constant (mM)
KA = 40 #ammonia saturation constant (mM)
Kdmax = 0.01 #maximum death rate (1/h)
mumax = 0.044 #maximum growth rate (1/h)
KG = 1 #glucose saturation constant (mM)
KQ = 0.22 #glutamine saturation constant (mM)
mQ = 0 #glutamine maintenance coefficient (mmol/cell/hour)
kmu = 0.01 #intrinsic death rate (1/h)
Klysis = 2*10**-2 #rate of cell lysis (1/h)
Ci_star = 100 #inhibitor saturation concentration (mM)
qi = 2.5*10**-10 #specific inhibitor production rate (1/h)
#Flow, volume and concentration
Fo = 0.001 #feed-rate (L/h)
Fi = 0.001 #feed-rate (L/h)
V = 3 #volume (L)
SG = 653 #glucose concentration in the feed (mM)
SQ = 58.8 #glutamine concentration in the feed (mM)
# create GEKKO parameter
t = np.linspace(0,120,121)
m.time = t
XT = m.Var(name='XT') #total cell density (cells/L)
XV = m.Var(lb=0, name='XV') #viable cell density (cells/L)
XD = m.Var(name='XD') #dead cell density (cells/L)
G = m.Var(value = 30, name='G') #glucose concentration (mM)
Q = m.Var(value = 4.5, name='Q') #glutamine concentration (mM)
L = m.Var(name='L') #lactate concentration (mM)
A = m.Var(name='A') #ammonia concentration (mM)
Ci = m.Var(name='Ci') #inhibitor concentration (mM)
mu = m.Var(name='mu') #growth rate (1/h)
Kd = m.Var(name='Kd') #death rate(1/h)
# create GEEKO equations
m.Equation(XT.dt() == mu*XV - Klysis*XD - XT*Fo/V)
m.Equation(XV.dt() == (mu - Kd)*XV - XV*Fo/V)
m.Equation(XD.dt() == Kd*XV - Klysis*XD - XV*Fo/V)
m.Equation(G.dt() == (Fi/V)*SG - (Fo/V)*G + (-mu/YXG - mG)*XV)
m.Equation(Q.dt() == (Fi/V)*SQ - (Fo/V)*Q + (-mu/YXQ - mQ)*XV - KdQ*Q)
m.Equation(L.dt() == -YLG*(-mu/YXG -mG)*XV-(Fo/V)*L)
m.Equation(A.dt() == -YAQ*(-mu/YXQ - mQ)*XV +KdQ*Q-(Fo/V)*A)
m.Equation(Ci.dt() == qi*XV - (Fo/V)*Ci)
m.Equation(mu.dt() == (mumax*G*Q*(Ci_star-Ci)) / (Ci_star*(KG+G)*(KQ+Q)*(L/KL + 1)*(A/KA + 1)))
m.Equation(Kd.dt() == Kdmax*(kmu/(mu+kmu)))
# solve ODE
m.options.IMODE = 4
m.open_folder()
m.solve(display = False)
plt.plot(m.time, XV.value)
使用完全相同模型的文章：
1) 使用GEKKO的硕士论文"MODELING OF MAMMALIAN CELL CULTURE" link:
2) 描述方程式的原始论文：“工艺模型比较和跨生物反应器规模的可转移性和 哺乳动物细胞生物过程的操作模式
link: https://sci-hub.tw/10.1002/btpr.1664
3) 控制系统使用该模型的论文："Glucose concentration control of a fed-batch mammalian cell bioprocess using a nonlinear model predictive controller"
link: https://sci-hub.tw/https://doi.org/10.1016/j.jprocont.2014.02.007
归根结底这不是编程问题，而是阅读方程式并正确翻译它们的问题。
mu 和 Kd 不是动态变量，它们是状态向量（只有 8 维）的普通函数。对于这样的中间变量，Gekko 具有构造函数 m.Intermediate
mu = m.Intermediate((mumax*G*Q*(Ci_star-Ci)) / (Ci_star*(KG+G)*(KQ+Q)*(L/KL + 1)*(A/KA + 1)), name='mu') #growth rate (1/h)
Kd = m.Intermediate(Kdmax*(kmu/(mu+kmu))) #death rate(1/h)
通过此更改以及 XT 和 XV 的初始值 5e8，脚本给出了下面的图表，这些图表显示了您引用的论文中可以找到的内容。
有几个问题：
- 最后两个方程是代数方程，不是微分方程。应该是
mu==...而不是mu.dt()==... - 一些方程式可能被零除。
x.dt() = z/y等式可以用y * x.dt()==z代替，这样如果y趋近于零，则等式变为0 * x.dt() == z。 - 一些初始条件未设置，因此它们使用默认值 0。这可能会创建一个零解。
我输入了一些不同的值并使用 m.options.COLDSTART=2 来帮助它找到初始解决方案。我还使用 Intermediates 来帮助可视化任何变大的术语。我将细胞浓度以每升数百万个细胞为单位来帮助缩放。
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
m = GEKKO(remote=False) # create GEKKO model
#constants 3L continuous fed-batch
KdQ = 0.001 #degree of degradation of glutamine (1/h)
mG = 1.1e-10 #glucose maintenance coefficient (mmol/cell/hour)
YAQ = 0.90 #yield of ammonia from glutamine
YLG = 2 #yield of lactate from glucose
YXG = 2.2e8 #yield of cells from glucose (cells/mmol)
YXQ = 1.5e9 #yield of cells from glutamine (cells/mmol)
KL = 150 #lactate saturation constant (mM)
KA = 40 #ammonia saturation constant (mM)
Kdmax = 0.01 #maximum death rate (1/h)
mumax = 0.044 #maximum growth rate (1/h)
KG = 1 #glucose saturation constant (mM)
KQ = 0.22 #glutamine saturation constant (mM)
mQ = 0 #glutamine maintenance coefficient (mmol/cell/hour)
kmu = 0.01 #intrinsic death rate (1/h)
Klysis = 2e-2 #rate of cell lysis (1/h)
Ci_star = 100 #inhibitor saturation concentration (mM)
qi = 2.5e-10 #specific inhibitor production rate (1/h)
#Flow, volume and concentration
Fo = 0.001 #feed-rate (L/h)
Fi = 0.001 #feed-rate (L/h)
V = 3 #volume (L)
SG = 653 #glucose concentration in the feed (mM)
SQ = 58.8 #glutamine concentration in the feed (mM)
# create GEKKO parameter
t = np.linspace(0,50,121)
m.time = t
XTMM = m.Var(value=1,name='XT') #total cell density (MMcells/L)
XVMM = m.Var(value=1,lb=0, name='XV') #viable cell density (MMcells/L)
XDMM = m.Var(value=1.0,name='XD') #dead cell density (MMcells/L)
G = m.Var(value = 20, name='G') #glucose concentration (mM)
Q = m.Var(value = 4.5, name='Q') #glutamine concentration (mM)
L = m.Var(value=1,name='L') #lactate concentration (mM)
A = m.Var(value=1.6,name='A') #ammonia concentration (mM)
Ci = m.Var(value=0.1,name='Ci') #inhibitor concentration (mM)
mu = m.Var(value=0.1,name='mu') #growth rate (1/h)
Kd = m.Var(value=0.5,name='Kd') #death rate(1/h)
# scale back to cells/L from million cells/L
XT = m.Intermediate(XTMM*1e7)
XV = m.Intermediate(XVMM*1e7)
XD = m.Intermediate(XDMM*1e7)
e1 = m.Intermediate((mu*XV - Klysis*XD - XT*Fo/V)/1e7)
e2 = m.Intermediate(((mu - Kd)*XV - XV*Fo/V)/1e7)
e3 = m.Intermediate((Kd*XV - Klysis*XD - XV*Fo/V)/1e7)
e4 = m.Intermediate((Fi/V)*SG - (Fo/V)*G + (-mu/YXG - mG)*XV)
e5 = m.Intermediate((Fi/V)*SQ - (Fo/V)*Q + (-mu/YXQ - mQ)*XV - KdQ*Q)
e6 = m.Intermediate(-YLG*(-mu/YXG -mG)*XV-(Fo/V)*L)
e7 = m.Intermediate(-YAQ*(-mu/YXQ - mQ)*XV +KdQ*Q-(Fo/V)*A)
e8 = m.Intermediate(qi*XV - (Fo/V)*Ci)
e9a = m.Intermediate((Ci_star*(KG+G)*(KQ+Q)*(L/KL + 1)*(A/KA + 1)))
e9b = m.Intermediate((mumax*G*Q*(Ci_star-Ci)))
e10a = m.Intermediate((mu+kmu))
e10b = m.Intermediate(Kdmax*kmu)
# create GEEKO equations
m.Equation(XTMM.dt() == e1)
m.Equation(XVMM.dt() == e2)
m.Equation(XDMM.dt() == e3)
m.Equation(G.dt() == e4)
m.Equation(Q.dt() == e5)
m.Equation(L.dt() == e6)
m.Equation(A.dt() == e7)
m.Equation(Ci.dt() == e8)
m.Equation(e9a * mu == e9b)
m.Equation(e10a*Kd == e10b)
# solve ODE
m.options.IMODE = 4
m.options.SOLVER = 1
m.options.NODES = 3
m.options.COLDSTART = 2
#m.open_folder()
m.solve(display=False)
plt.figure()
plt.subplot(3,1,1)
plt.plot(m.time, XV.value,label='XV')
plt.plot(m.time, XT.value,label='XT')
plt.plot(m.time, XD.value,label='XD')
plt.legend()
plt.subplot(3,1,2)
plt.plot(m.time, G.value,label='G')
plt.plot(m.time, Q.value,label='Q')
plt.plot(m.time, L.value,label='L')
plt.plot(m.time, A.value,label='A')
plt.legend()
plt.subplot(3,1,3)
plt.plot(m.time, mu.value,label='mu')
plt.plot(m.time, Kd.value,label='Kd')
plt.legend()
plt.xlabel('Time (hr)')
plt.figure()
plt.plot(m.time, e1.value,'r-.',label='eqn1')
plt.plot(m.time, e2.value,'g:',label='eqn2')
plt.plot(m.time, e3.value,'b:',label='eqn3')
plt.plot(m.time, e4.value,'b--',label='eqn4')
plt.plot(m.time, e5.value,'y:',label='eqn5')
plt.plot(m.time, e6.value,'m--',label='eqn6')
plt.plot(m.time, e7.value,'b-.',label='eqn7')
plt.plot(m.time, e8.value,'g--',label='eqn8')
plt.plot(m.time, e9a.value,'r:',label='eqn9a')
plt.plot(m.time, e9b.value,'r--',label='eqn9b')
plt.plot(m.time, e10a.value,'k:',label='eqn10a')
plt.plot(m.time, e10b.value,'k--',label='eqn10b')
plt.legend()
plt.show()