Python GEKKO 中基于 ODE 或 PDE 的生态系统模型最合适的求解方法是什么？
What is the most appropriate solving method for ODE or PDE based ecosystem models in Python GEKKO?
我已经找了一段时间了，但在任何地方都找不到这个特定问题的答案，抱歉，如果它是重复的！
我已经开始构建 python package based on the xarray-simlab framework，目标是提供一个模块化工具箱，用于构建可重复且灵活的海洋生态系统模型。 Xarray-simlab 目前仅支持显式步长来求解模型函数。为了更安全高效地求解复杂模型，我开始使用 GEKKO 作为求解器后端，因为模型语法似乎很合适。 （注意：目前我只需要随着时间的推移求解模型方程的功能，但我想利用 GEKKO 的优化功能在后期将模型参数拟合到现场或实验室数据。）
包的当前原型创建了一个传递GEKKO模型实例的xsimlab进程class
m 所有子流程。继承模型实例的进程 classes 根据添加到模型的进程和运行时提供的参数（包括 SV 维度）初始化 m.SV、m.Param 或定义 m.Intermediates .在下一步中，所有初始化的中间体都累积到 m.Equations 中受影响的状态变量。一旦成功求解，GEKKO 变量将被重新打包到一个 xarray 数据结构中，其中包括相关的元数据，可以进一步分析。包原型可以使用 IMODE=7 求解基本模型，但我遇到了一个与该求解器的时间步长相关的问题：
我期待类似于 scipy 的 odeint 的功能，具有自适应时间步长评估，但显然情况似乎并非如此，而是它在提供的离散时间步长上评估模型。
该软件包仍在大力开发中，还有很多我仍在努力改进的功能，因此下面是一个简单恒化器模型的最小代码示例。该模型描述了在简化的流通系统中在营养物上生长的浮游植物状态变量。养分以恒定速率流入，浮游植物以恒定速率死亡并从系统中流失：
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
m = GEKKO() # create GEKKO model
halfsat_const = m.Param(0.1)
N0 = m.Param(1.)
inflow_rate = m.Param(0.1)
mortality_rate = m.Param(0.1)
N = m.SV(1)
P = m.SV(0.1)
t = np.arange(0,10,0.01)
m.time = t
# Growth under nutrient limitation is described via Monod / Michaelis-Menten kinetics
nutlim = m.Intermediate(N/(N+halfsat_const)*P)
N_influx = m.Intermediate(N0 * inflow_rate)
mortality = m.Intermediate(P * mortality_rate)
m.Equation(N.dt()==N_influx - nutlim)
m.Equation(P.dt()==nutlim - mortality)
m.options.IMODE = 7
m.solve(disp=False)
plt.plot(m.time, N, label='N')
plt.plot(m.time, P, label='P')
plt.legend()
这对于提供的时间步长非常有效，但是例如m.time = np.arange(0,10) returns 一个荒谬的解决方案（两条不同的线达到 >1e7）。 Odeint 解决没问题:
import numpy as np
from scipy.integrate import odeint
halfsat = 0.1
N0 = 1.
inflow = 0.1
mortality_rate = 0.1
def model(y,t):
N,P = y
nutlim = N/(N+halfsat)*P
influx = N0 * inflow
mortality = P * mortality_rate
dNdt = influx - nutlim
dPdt = nutlim - mortality
return [dNdt, dPdt]
model_time = np.arange(0,10)
out = odeint(model,[1,0.1],model_time)
plt.plot(model_time,out[:,0], label='N')
plt.plot(model_time,out[:,1], label='P')
plt.legend()
我用我的包构建的模型可能会变得相对复杂，有数百个状态变量，以及更多的交互，产生高度非线性的结果。我不确定如何确定我提供的时间步长是否合适，因为较小的时间步长会显着增加计算时间。
是否有 GEKKO 包含的求解器（或与 GEKKO 模型语法兼容）提供与具有自适应步长的 odeint 类似的求解器？或者是否有另一种方法更适合处理基于 ODE（或空间离散 PDE 系统）的生态模型？
非常感谢任何帮助！
尝试增加每段的节点数：
m.options.NODES = 3
这给出了更准确的解决方案，因为使用了更高阶的搭配方法。在这种情况下，时间点 [0,1,2,...9,10] 对于准确的解决方案来说过于粗糙，但 [0,0.5,1,...9.5,10] 工作正常。
此外，通过m.SV(lb=0)将状态变量的下边界设置为零将提高求解稳定性。这是生态系统模型的一个基本假设，即由例如跟踪的组件。生物量不会是负数。
我通常建议进行网格独立性测试，您可以在其中减小步长直到解不改变或与自适应步长求解器（例如 ODEINT）进行比较。 Gekko 确实为 IMODE=7 执行自适应步长，但仅当求解器在某个步骤上失败时才这样做。由用户决定离散化。 Gekko 的优势在于优化，优化中的自适应步长需要一个可能非常慢的多级策略。但是，已经有recent progress. If you'd like to have an adaptive step size with IMODE=7 and error checking, please consider a feature request.
import numpy as np
from gekko import GEKKO
from scipy.integrate import odeint
import matplotlib.pyplot as plt
m = GEKKO(remote=False) # create GEKKO model
halfsat_const = m.Param(0.1)
N0 = m.Param(1.)
inflow_rate = m.Param(0.1)
mortality_rate = m.Param(0.1)
N = m.SV(1, lb=0)
P = m.SV(0.1, lb=0)
t = np.arange(0,10,0.2)
m.time = t
# Growth under nutrient limitation is described via Monod / Michaelis-Menten kinetics
nutlim = m.Intermediate(N/(N+halfsat_const)*P)
N_influx = m.Intermediate(N0 * inflow_rate)
mortality = m.Intermediate(P * mortality_rate)
m.Equation(N.dt()==N_influx - nutlim)
m.Equation(P.dt()==nutlim - mortality)
m.options.NODES = 3
m.options.IMODE = 7
m.solve(disp=False)
halfsat = 0.1
N0 = 1.
inflow = 0.1
mortality_rate = 0.1
def model(y,t):
N,P = y
nutlim = N/(N+halfsat)*P
influx = N0 * inflow
mortality = P * mortality_rate
dNdt = influx - nutlim
dPdt = nutlim - mortality
return [dNdt, dPdt]
model_time = np.arange(0,10)
out = odeint(model,[1,0.1],model_time)
plt.plot(model_time,out[:,0], 'ro', label='N ODEINT')
plt.plot(model_time,out[:,1], 'bx', label='P ODEINT')
plt.plot(m.time, N, 'r--', label='N Gekko')
plt.plot(m.time, P, 'b--', label='P Gekko')
plt.legend()
plt.show()
我已经找了一段时间了，但在任何地方都找不到这个特定问题的答案，抱歉，如果它是重复的！
我已经开始构建 python package based on the xarray-simlab framework，目标是提供一个模块化工具箱，用于构建可重复且灵活的海洋生态系统模型。 Xarray-simlab 目前仅支持显式步长来求解模型函数。为了更安全高效地求解复杂模型，我开始使用 GEKKO 作为求解器后端，因为模型语法似乎很合适。 （注意：目前我只需要随着时间的推移求解模型方程的功能，但我想利用 GEKKO 的优化功能在后期将模型参数拟合到现场或实验室数据。）
包的当前原型创建了一个传递GEKKO模型实例的xsimlab进程class
m 所有子流程。继承模型实例的进程 classes 根据添加到模型的进程和运行时提供的参数（包括 SV 维度）初始化 m.SV、m.Param 或定义 m.Intermediates .在下一步中，所有初始化的中间体都累积到 m.Equations 中受影响的状态变量。一旦成功求解，GEKKO 变量将被重新打包到一个 xarray 数据结构中，其中包括相关的元数据，可以进一步分析。包原型可以使用 IMODE=7 求解基本模型，但我遇到了一个与该求解器的时间步长相关的问题：
我期待类似于 scipy 的 odeint 的功能，具有自适应时间步长评估，但显然情况似乎并非如此，而是它在提供的离散时间步长上评估模型。
该软件包仍在大力开发中，还有很多我仍在努力改进的功能，因此下面是一个简单恒化器模型的最小代码示例。该模型描述了在简化的流通系统中在营养物上生长的浮游植物状态变量。养分以恒定速率流入，浮游植物以恒定速率死亡并从系统中流失：
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
m = GEKKO() # create GEKKO model
halfsat_const = m.Param(0.1)
N0 = m.Param(1.)
inflow_rate = m.Param(0.1)
mortality_rate = m.Param(0.1)
N = m.SV(1)
P = m.SV(0.1)
t = np.arange(0,10,0.01)
m.time = t
# Growth under nutrient limitation is described via Monod / Michaelis-Menten kinetics
nutlim = m.Intermediate(N/(N+halfsat_const)*P)
N_influx = m.Intermediate(N0 * inflow_rate)
mortality = m.Intermediate(P * mortality_rate)
m.Equation(N.dt()==N_influx - nutlim)
m.Equation(P.dt()==nutlim - mortality)
m.options.IMODE = 7
m.solve(disp=False)
plt.plot(m.time, N, label='N')
plt.plot(m.time, P, label='P')
plt.legend()
这对于提供的时间步长非常有效，但是例如m.time = np.arange(0,10) returns 一个荒谬的解决方案（两条不同的线达到 >1e7）。 Odeint 解决没问题:
import numpy as np
from scipy.integrate import odeint
halfsat = 0.1
N0 = 1.
inflow = 0.1
mortality_rate = 0.1
def model(y,t):
N,P = y
nutlim = N/(N+halfsat)*P
influx = N0 * inflow
mortality = P * mortality_rate
dNdt = influx - nutlim
dPdt = nutlim - mortality
return [dNdt, dPdt]
model_time = np.arange(0,10)
out = odeint(model,[1,0.1],model_time)
plt.plot(model_time,out[:,0], label='N')
plt.plot(model_time,out[:,1], label='P')
plt.legend()
我用我的包构建的模型可能会变得相对复杂，有数百个状态变量，以及更多的交互，产生高度非线性的结果。我不确定如何确定我提供的时间步长是否合适，因为较小的时间步长会显着增加计算时间。
是否有 GEKKO 包含的求解器（或与 GEKKO 模型语法兼容）提供与具有自适应步长的 odeint 类似的求解器？或者是否有另一种方法更适合处理基于 ODE（或空间离散 PDE 系统）的生态模型？
非常感谢任何帮助！
尝试增加每段的节点数：
m.options.NODES = 3
这给出了更准确的解决方案，因为使用了更高阶的搭配方法。在这种情况下，时间点 [0,1,2,...9,10] 对于准确的解决方案来说过于粗糙，但 [0,0.5,1,...9.5,10] 工作正常。
此外，通过m.SV(lb=0)将状态变量的下边界设置为零将提高求解稳定性。这是生态系统模型的一个基本假设，即由例如跟踪的组件。生物量不会是负数。
我通常建议进行网格独立性测试，您可以在其中减小步长直到解不改变或与自适应步长求解器（例如 ODEINT）进行比较。 Gekko 确实为 IMODE=7 执行自适应步长，但仅当求解器在某个步骤上失败时才这样做。由用户决定离散化。 Gekko 的优势在于优化，优化中的自适应步长需要一个可能非常慢的多级策略。但是，已经有recent progress. If you'd like to have an adaptive step size with IMODE=7 and error checking, please consider a feature request.
import numpy as np
from gekko import GEKKO
from scipy.integrate import odeint
import matplotlib.pyplot as plt
m = GEKKO(remote=False) # create GEKKO model
halfsat_const = m.Param(0.1)
N0 = m.Param(1.)
inflow_rate = m.Param(0.1)
mortality_rate = m.Param(0.1)
N = m.SV(1, lb=0)
P = m.SV(0.1, lb=0)
t = np.arange(0,10,0.2)
m.time = t
# Growth under nutrient limitation is described via Monod / Michaelis-Menten kinetics
nutlim = m.Intermediate(N/(N+halfsat_const)*P)
N_influx = m.Intermediate(N0 * inflow_rate)
mortality = m.Intermediate(P * mortality_rate)
m.Equation(N.dt()==N_influx - nutlim)
m.Equation(P.dt()==nutlim - mortality)
m.options.NODES = 3
m.options.IMODE = 7
m.solve(disp=False)
halfsat = 0.1
N0 = 1.
inflow = 0.1
mortality_rate = 0.1
def model(y,t):
N,P = y
nutlim = N/(N+halfsat)*P
influx = N0 * inflow
mortality = P * mortality_rate
dNdt = influx - nutlim
dPdt = nutlim - mortality
return [dNdt, dPdt]
model_time = np.arange(0,10)
out = odeint(model,[1,0.1],model_time)
plt.plot(model_time,out[:,0], 'ro', label='N ODEINT')
plt.plot(model_time,out[:,1], 'bx', label='P ODEINT')
plt.plot(m.time, N, 'r--', label='N Gekko')
plt.plot(m.time, P, 'b--', label='P Gekko')
plt.legend()
plt.show()