如何使用 GEKKO 描述离散参数？
How to describe a discrete parameter using GEKKO?
我有一个包含 10 个描述补料分批生物反应器的方程式的模型。基本上，每 24 小时 "food"（葡萄糖和其他成分）就会添加到系统中。我的问题是这个喂养程序目前被建模为两个时间步长 (delta_T) 的流速 (L/H)，而不是单个离散食物添加 (delta_T = 0 ).
葡萄糖方程式如下所示：
e4 = m.Intermediate(**(Fi/V)*SG** - (Fo/V)*G + (-mu/YXG - mG)*XV)
m.Equation(G.dt() == e4)
其中 G 是生物反应器中的葡萄糖浓度 (mM)，Fi 是输入进料速率 (L/h)，V 是生物反应器体积 ( L), SG为饲料中的葡萄糖浓度(mM)。
我设法通过调用此 delta_T = 0.2 hours 使系统可行，换句话说，葡萄糖浓度连续（而不是立即) 从时间 t1 的 G1 上升到时间 t1 + 0.2h 的 G2。如果我尝试降低此 delta_T，系统会显示出非常奇怪的行为。
时间数组如下所示：[..., 19.5, 20.0, 20.5, 21.0, 21.5, 22.0, 22.5, 23.0, 23.5, 24.0, 24.2, 24.5, 25.0，...]。
它以 0.5 小时为步长变化，每 24 小时，当我向生物反应器中添加葡萄糖时，我强制下一个时间步长只有 0.2，而不是 0.5。我希望这个增量为 0。
我的进给率是这样的：
[..., 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0,...]
有没有办法让这个喂食过程离散化？完整代码如下所示。谢谢！！
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
import pandas as pd
m = GEKKO(remote=False) # create GEKKO model
#constants 3L 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)
N_HOURS = 150 #Number of hours of the experiment
TIME_STEP = 0.5
feed_small_delta_t = 0.2
#create time array. It will be from 0 to N_HOURS, with time step TIME_STEP,
#and every 24h, there will be a feed_small_delta_t
t = []
for i in range(int(N_HOURS/TIME_STEP +1)):
t_value = i*TIME_STEP
t.append(t_value)
if t_value%24 == 0:
t.append(t_value + feed_small_delta_t)
m.time = t
#Create input feed-rate array
Fi = np.zeros(len(t))
for i in range(1,len(t)):
if t[i]%(24) == 0:
Fi[i] = 0.1/feed_small_delta_t #it is divided by feed_small_delta_t so volume added is constant.
#Flow, volume and concentration
Fi = m.Param(Fi) #input feed-rate (L/h)
Fo = 0 #output 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)
XTMM = m.Var(value=2,lb=-0.0000,name='XT') #total cell density (MMcells/L)
XVMM = m.Var(value=2,lb=-0.0000, name='XV') #viable cell density (MMcells/L)
XDMM = m.Var(value=0,lb=-0.0000,name='XD') #dead cell density (MMcells/L)
G = m.Var(value = 20,lb=-0.0000,name='G') #glucose concentration (mM)
Q = m.Var(value = 3.8,lb=-0.0000, name='Q') #glutamine concentration (mM)
L = m.Var(value=0.1,lb=-0.0000,name='L') #lactate concentration (mM)
A = m.Var(value=1.8,lb=-0.0000,name='A') #ammonia concentration (mM)
Ci = m.Var(lb=-0.0000, name='Ci') #inhibitor concentration (mM)
mu = m.Var(lb=-0.0000, name='mu') #growth rate (1/h)
Kd = m.Var(lb=-0.0000, name='Kd') #death rate(1/h)
# scale back to cells/L from million cells/L
XT = m.Intermediate(XTMM*1e8)
XV = m.Intermediate(XVMM*1e8)
XD = m.Intermediate(XDMM*1e8)
e1 = m.Intermediate((mu*XV - Klysis*XD - XT*Fo/V)/1e8)
e2 = m.Intermediate(((mu - Kd)*XV - XV*Fo/V)/1e8)
e3 = m.Intermediate((Kd*XV - Klysis*XD - XV*Fo/V)/1e8)
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.plot(m.time, Ci.value,label='Ci')
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()
您将葡萄糖作为脉冲喂入的策略是获得不连续输入的好方法。葡萄糖浓度离散跳跃的问题在于，存在方程式 4 中的葡萄糖导数项：m.Equation(G.dt() == e4)。如果 dG/dt 项在很短的时间内发生变化，则导数项会变得非常大。
处理离散点处的大导数的一种策略是使用m.options.NODES=2来避免在有限元上具有正交配置的附加内部节点的问题。如果没有内部节点，您可能需要增加时间点的数量以保持集成的准确性。这允许非常短持续时间的葡萄糖脉冲输入到间歇式反应器，例如 3.6 seconds 用于添加。
feed_small_delta_t = 0.001 # 3.6 seconds
馈送输入的索引偏移了 1，因此 Fi[i+1] 应该是施加脉冲的位置，而不是 Fi[i]。
Fi = np.zeros(len(t))
for i in range(1,len(t)):
if t[i]%(24) == 0:
Fi[i+1] = 0.1/feed_small_delta_t
这会得到与您之前得到的结果类似的结果，但对于向批次中添加额外糖分的离散事件，输入脉冲更短。
修改后的完整脚本
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
import pandas as pd
m = GEKKO(remote=False) # create GEKKO model
#constants 3L 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)
N_HOURS = 150 #Number of hours of the experiment
TIME_STEP = 0.5
feed_small_delta_t = 0.001
#create time array. It will be from 0 to N_HOURS, with time step TIME_STEP,
#and every 24h, there will be a feed_small_delta_t
t = []
for i in range(int(N_HOURS/TIME_STEP +1)):
t_value = i*TIME_STEP
t.append(t_value)
if t_value%24 == 0:
t.append(t_value + feed_small_delta_t)
m.time = t
#Create input feed-rate array
Fi = np.zeros(len(t))
for i in range(1,len(t)):
if t[i]%(24) == 0:
Fi[i+1] = 0.1/feed_small_delta_t #it is divided by feed_small_delta_t so volume added is constant.
#Flow, volume and concentration
Fi = m.Param(Fi) #input feed-rate (L/h)
Fo = 0 #output 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)
XTMM = m.Var(value=2,lb=-0.0000,name='XT') #total cell density (MMcells/L)
XVMM = m.Var(value=2,lb=-0.0000, name='XV') #viable cell density (MMcells/L)
XDMM = m.Var(value=0,lb=-0.0000,name='XD') #dead cell density (MMcells/L)
G = m.Var(value = 20,lb=-0.0000,name='G') #glucose concentration (mM)
Q = m.Var(value = 3.8,lb=-0.0000, name='Q') #glutamine concentration (mM)
L = m.Var(value=0.1,lb=-0.0000,name='L') #lactate concentration (mM)
A = m.Var(value=1.8,lb=-0.0000,name='A') #ammonia concentration (mM)
Ci = m.Var(lb=-0.0000, name='Ci') #inhibitor concentration (mM)
mu = m.Var(lb=-0.0000, name='mu') #growth rate (1/h)
Kd = m.Var(lb=-0.0000, name='Kd') #death rate(1/h)
# scale back to cells/L from million cells/L
XT = m.Intermediate(XTMM*1e8)
XV = m.Intermediate(XVMM*1e8)
XD = m.Intermediate(XDMM*1e8)
e1 = m.Intermediate((mu*XV - Klysis*XD - XT*Fo/V)/1e8)
e2 = m.Intermediate(((mu - Kd)*XV - XV*Fo/V)/1e8)
e3 = m.Intermediate((Kd*XV - Klysis*XD - XV*Fo/V)/1e8)
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 = 2
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.plot(m.time, Ci.value,label='Ci')
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()
我有一个包含 10 个描述补料分批生物反应器的方程式的模型。基本上，每 24 小时 "food"（葡萄糖和其他成分）就会添加到系统中。我的问题是这个喂养程序目前被建模为两个时间步长 (delta_T) 的流速 (L/H)，而不是单个离散食物添加 (delta_T = 0 ).
葡萄糖方程式如下所示：
e4 = m.Intermediate(**(Fi/V)*SG** - (Fo/V)*G + (-mu/YXG - mG)*XV)
m.Equation(G.dt() == e4)
其中 G 是生物反应器中的葡萄糖浓度 (mM)，Fi 是输入进料速率 (L/h)，V 是生物反应器体积 ( L), SG为饲料中的葡萄糖浓度(mM)。
我设法通过调用此 delta_T = 0.2 hours 使系统可行，换句话说，葡萄糖浓度连续（而不是立即) 从时间 t1 的 G1 上升到时间 t1 + 0.2h 的 G2。如果我尝试降低此 delta_T，系统会显示出非常奇怪的行为。
时间数组如下所示：[..., 19.5, 20.0, 20.5, 21.0, 21.5, 22.0, 22.5, 23.0, 23.5, 24.0, 24.2, 24.5, 25.0，...]。
它以 0.5 小时为步长变化，每 24 小时，当我向生物反应器中添加葡萄糖时，我强制下一个时间步长只有 0.2，而不是 0.5。我希望这个增量为 0。
我的进给率是这样的：
[..., 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0,...]
有没有办法让这个喂食过程离散化？完整代码如下所示。谢谢！！
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
import pandas as pd
m = GEKKO(remote=False) # create GEKKO model
#constants 3L 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)
N_HOURS = 150 #Number of hours of the experiment
TIME_STEP = 0.5
feed_small_delta_t = 0.2
#create time array. It will be from 0 to N_HOURS, with time step TIME_STEP,
#and every 24h, there will be a feed_small_delta_t
t = []
for i in range(int(N_HOURS/TIME_STEP +1)):
t_value = i*TIME_STEP
t.append(t_value)
if t_value%24 == 0:
t.append(t_value + feed_small_delta_t)
m.time = t
#Create input feed-rate array
Fi = np.zeros(len(t))
for i in range(1,len(t)):
if t[i]%(24) == 0:
Fi[i] = 0.1/feed_small_delta_t #it is divided by feed_small_delta_t so volume added is constant.
#Flow, volume and concentration
Fi = m.Param(Fi) #input feed-rate (L/h)
Fo = 0 #output 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)
XTMM = m.Var(value=2,lb=-0.0000,name='XT') #total cell density (MMcells/L)
XVMM = m.Var(value=2,lb=-0.0000, name='XV') #viable cell density (MMcells/L)
XDMM = m.Var(value=0,lb=-0.0000,name='XD') #dead cell density (MMcells/L)
G = m.Var(value = 20,lb=-0.0000,name='G') #glucose concentration (mM)
Q = m.Var(value = 3.8,lb=-0.0000, name='Q') #glutamine concentration (mM)
L = m.Var(value=0.1,lb=-0.0000,name='L') #lactate concentration (mM)
A = m.Var(value=1.8,lb=-0.0000,name='A') #ammonia concentration (mM)
Ci = m.Var(lb=-0.0000, name='Ci') #inhibitor concentration (mM)
mu = m.Var(lb=-0.0000, name='mu') #growth rate (1/h)
Kd = m.Var(lb=-0.0000, name='Kd') #death rate(1/h)
# scale back to cells/L from million cells/L
XT = m.Intermediate(XTMM*1e8)
XV = m.Intermediate(XVMM*1e8)
XD = m.Intermediate(XDMM*1e8)
e1 = m.Intermediate((mu*XV - Klysis*XD - XT*Fo/V)/1e8)
e2 = m.Intermediate(((mu - Kd)*XV - XV*Fo/V)/1e8)
e3 = m.Intermediate((Kd*XV - Klysis*XD - XV*Fo/V)/1e8)
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.plot(m.time, Ci.value,label='Ci')
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()
您将葡萄糖作为脉冲喂入的策略是获得不连续输入的好方法。葡萄糖浓度离散跳跃的问题在于，存在方程式 4 中的葡萄糖导数项：m.Equation(G.dt() == e4)。如果 dG/dt 项在很短的时间内发生变化，则导数项会变得非常大。
处理离散点处的大导数的一种策略是使用m.options.NODES=2来避免在有限元上具有正交配置的附加内部节点的问题。如果没有内部节点，您可能需要增加时间点的数量以保持集成的准确性。这允许非常短持续时间的葡萄糖脉冲输入到间歇式反应器，例如 3.6 seconds 用于添加。
feed_small_delta_t = 0.001 # 3.6 seconds
馈送输入的索引偏移了 1，因此 Fi[i+1] 应该是施加脉冲的位置，而不是 Fi[i]。
Fi = np.zeros(len(t))
for i in range(1,len(t)):
if t[i]%(24) == 0:
Fi[i+1] = 0.1/feed_small_delta_t
这会得到与您之前得到的结果类似的结果，但对于向批次中添加额外糖分的离散事件，输入脉冲更短。
修改后的完整脚本
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
import pandas as pd
m = GEKKO(remote=False) # create GEKKO model
#constants 3L 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)
N_HOURS = 150 #Number of hours of the experiment
TIME_STEP = 0.5
feed_small_delta_t = 0.001
#create time array. It will be from 0 to N_HOURS, with time step TIME_STEP,
#and every 24h, there will be a feed_small_delta_t
t = []
for i in range(int(N_HOURS/TIME_STEP +1)):
t_value = i*TIME_STEP
t.append(t_value)
if t_value%24 == 0:
t.append(t_value + feed_small_delta_t)
m.time = t
#Create input feed-rate array
Fi = np.zeros(len(t))
for i in range(1,len(t)):
if t[i]%(24) == 0:
Fi[i+1] = 0.1/feed_small_delta_t #it is divided by feed_small_delta_t so volume added is constant.
#Flow, volume and concentration
Fi = m.Param(Fi) #input feed-rate (L/h)
Fo = 0 #output 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)
XTMM = m.Var(value=2,lb=-0.0000,name='XT') #total cell density (MMcells/L)
XVMM = m.Var(value=2,lb=-0.0000, name='XV') #viable cell density (MMcells/L)
XDMM = m.Var(value=0,lb=-0.0000,name='XD') #dead cell density (MMcells/L)
G = m.Var(value = 20,lb=-0.0000,name='G') #glucose concentration (mM)
Q = m.Var(value = 3.8,lb=-0.0000, name='Q') #glutamine concentration (mM)
L = m.Var(value=0.1,lb=-0.0000,name='L') #lactate concentration (mM)
A = m.Var(value=1.8,lb=-0.0000,name='A') #ammonia concentration (mM)
Ci = m.Var(lb=-0.0000, name='Ci') #inhibitor concentration (mM)
mu = m.Var(lb=-0.0000, name='mu') #growth rate (1/h)
Kd = m.Var(lb=-0.0000, name='Kd') #death rate(1/h)
# scale back to cells/L from million cells/L
XT = m.Intermediate(XTMM*1e8)
XV = m.Intermediate(XVMM*1e8)
XD = m.Intermediate(XDMM*1e8)
e1 = m.Intermediate((mu*XV - Klysis*XD - XT*Fo/V)/1e8)
e2 = m.Intermediate(((mu - Kd)*XV - XV*Fo/V)/1e8)
e3 = m.Intermediate((Kd*XV - Klysis*XD - XV*Fo/V)/1e8)
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 = 2
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.plot(m.time, Ci.value,label='Ci')
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()