使用Gekko Python无法对生物系统进行动态仿真

问题描述

我正在尝试模拟一个基本的生物系统：废水供给的反应器中活性污泥的生长。

它有两个变量和两个反应。 这是一个非常基本的模型，在废水工程中得到了广泛的使用（更详细）。

求解器似乎非常不稳定：经常找不到解决方案（大型超调）。 模型非常平滑，初始后没有快速变化 小跳。

我尝试了IMODE = 4（动态同步），但是大多数时候都失败了。 当前以IMODE = 7（动态顺序）运行。 可以在0.1到4天的时间间隔内工作，但是会在更长的时间（5天以上）内失败。这样一来，系统在0.2天后或多或少就会稳定（请参阅下文）。

求解器还需要很长时间才能解决。当我希望将其用于更大的工作时。

我的问题： 我的代码有问题吗，还是它的求解器？ 可以改进代码以加快仿真速度吗？

该测试在装有以下软件包的Jupyter笔记本（chrome）中运行：gekko，numpy，matplotlib。

# Basic setup:
# two variables: S en X
# two reaction: aerobic growth and decay
# one reactor: in and out

# ASM parameters (activated sludge model)
Y = 0.67 # g/g COD
fp = 0.08  # -
mu_h = 6.0  # /d
Ks = 20 # g/m3 COD
b_h = 0.62  # /d

# setup parameters
V = 100  # m3,reactor volume
Q = 50  # m3/d,flow rate to reactor
Sin = 1000 # g/m3
Xin = 0 # g/m3

# Build gekko model:
m = GEKKO()
m.options.IMODE = 7  # dynamic,sequential,simulation
m.time = np.linspace(start=0,stop=3,num=100)  # FAILS AT LARGER TIME FRAME

# create state variables
S = m.Var(value=100)  # substrate 
X = m.Var(value=500) # biomass

# intermediates: reaction rates
r_ag = m.Intermediate(mu_h*S/(Ks + S)*X)  # reaction rate,aerobic growth
r_ad = m.Intermediate(b_h*X)  # rate,aerobic decay

# constraint equations
m.Equation(V*S.dt() == Q*Sin - Q*S +(-1/Y*r_ag + (1-fp)*r_ad)*V)  # mass balance for S
m.Equation(V*X.dt() == Q*Xin - Q*X + (1*r_ag -1*r_ad)*V)  # mass balance for X

# solve
m.solve(disp=False)

# plot results
plt.plot(m.time,S,'g:',label='S')
plt.plot(m.time,X,'b-',label='X')
plt.legend(loc='best')

解决方法

Gekko对IMODE=7的自适应时间步细化仍在开发中（track feature request on GitHub）。当求解器认为生物量浓度应为负值时，解决方案将失败。在开始时有指数的有氧生长，有氧衰减非常缓慢。如果只需要仿真，建议您使用自适应时步求解器，例如scipy.integrate.ODEINT（tutorials）。但是，如果您需要执行优化（例如回归或批次控制），则Gekko是一个不错的选择。要使用IMODE=4或IMODE=7成功解决问题，请注意以下几点提示：

1. 添加约束以防止生物量或底物变负。

S = m.Var(value=100,lb=0)  # substrate 
X = m.Var(value=500,lb=0)  # biomass

1. 将时间步长设置为1，最终时间可以通过更改n进行调整。添加时间点0.001,0.002,...,0.8，以快速解决好氧运动的动态变化。自适应时间步细化完成后，无需插入额外的时间步长。

n = 100
t = np.linspace(start=0,stop=n,num=n+1)
m.time = np.insert(t,1,[0.001,0.004,0.008,0.02,0.04,0.08,\
                   0.2,0.4,0.8])

1. 本地解决，而不使用默认服务器。

m = GEKKO(remote=False)

这是完整的脚本，可以在任何最终时间和IMODE=4或IMODE=7的情况下成功解决。

# Basic setup:
# two variables: S en X
# two reaction: aerobic growth and decay
# one reactor: in and out

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

# ASM parameters (activated sludge model)
Y = 0.67 # g/g COD
fp = 0.08  # -
mu_h = 6.0  # /d
Ks = 20 # g/m3 COD
b_h = 0.62  # /d

# setup parameters
V = 100  # m3,reactor volume
Q = 50  # m3/d,flow rate to reactor
Sin = 1000 # g/m3
Xin = 0 # g/m3

# Build gekko model:
m = GEKKO(remote=False)
m.options.IMODE = 7  # dynamic,sequential,simulation
n = 100
t = np.linspace(start=0,0.8])

# create state variables
S = m.Var(value=100,lb=0)  # biomass

# intermediates: reaction rates
r_ag = m.Intermediate(mu_h*(S/(Ks + S))*X)  # reaction rate,aerobic growth
r_ad = m.Intermediate(b_h*X)  # rate,aerobic decay

# constraint equations
m.Equation(V*S.dt() == Q*Sin - Q*S +(-(1/Y)*r_ag + (1-fp)*r_ad)*V)  # mass balance for S
m.Equation(V*X.dt() == Q*Xin - Q*X + (1*r_ag -1*r_ad)*V)  # mass balance for X

# solve
m.options.NODES=3
m.solve(disp=False)

# plot results
plt.plot(m.time,S,'g:',label='S')
plt.plot(m.time,X,'b-',label='X')
plt.legend(loc='best')
plt.show()