对于Python GEKKO中基于ODE或PDE的生态系统模型，最合适的求解方法是什么？

Question

我已经找了一段时间了，但是在任何地方都找不到这个问题的答案，如果它是重复的，很抱歉！

我已经开始构建一个基于xarray-simlab framework的python package，目标是提供一个模块化工具箱，用于构建可重现的、灵活的海洋生态系统模型。Xarray-simlab目前只支持显式步长来求解模型函数。为了更安全高效地解决复杂模型，我开始使用GEKKO作为求解器后端，因为模型语法似乎很适合。(注意:目前，我只需要随时间推移求解模型方程的功能，但我希望在以后的阶段利用GEKKO的优化功能将模型参数拟合到现场或实验室数据。)

包的当前原型创建了一个xsimlab流程类，该流程类将GEKKO模型实例m传递给所有子流程。继承模型实例的流程类基于添加到模型和参数的流程初始化m.SV、m.Param或定义m.Intermediates。SV尺寸)在运行时提供。在下一步中，所有初始化的中间件都会累积到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)返回一个无意义的解决方案(两条不同的线到达>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类似的具有自适应步长的求解器？或者，是否有另一种方法更适合于处理基于PDE (或空间离散的PDE系统)的生态模型？

任何帮助都是非常感谢的！

Answer 1

尝试使用以下命令增加每个网段的节点数：

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。如果你想要一个具有IMODE=7和错误检查的自适应步长，请考虑使用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()