I have been looking for a while, but could not find the answer to this specific question anywhere, sorry if it is a duplicate!
I have started to build a python package based on the xarray-simlab framework with the goal to provide a modular toolbox for building reproducible and flexible marine ecosystem models. Xarray-simlab at the moment only supports explicit step-sizes to solve the model functions. In order to solve complex models more safely & efficiently, I have instead started using GEKKO as a solver backend, as the model syntax seems well suited. (Note: At the moment I will only need functionality to solve the model equations over time, but I would like to make use of GEKKO's optimization functionality to fit model parameters to field or lab data at later stages.)
The current prototype of the package creates a xsimlab process class that passes the GEKKO model instance
m to all sub-processes. Process classes that inherit the model instance initialize m.SV, m.Param or define m.Intermediates based on the processes added to the model & parameters (incl. SV dimensions) supplied at runtime. In the next step all initialized intermediates are accumulated to the affected state variables in m.Equations. Once successfully solved, GEKKO variables are repackaged into a xarray data structure, that includes relevant metadata and can be analysed further. The package prototype can solve basic models using IMODE=7, but I have come across one issue related to the time steps of that solver:
I was expecting functionality similar to scipy's odeint, with adaptive time step evaluation, but obviously this does not seem to be the case and instead it evaluates the model at the discrete time-steps supplied.
The package is still under heavy development, and there are plenty of features that I am still trying to improve, so below is a minimal code example of a simple chemostat model. The model describes a phytoplankton state variable growing on a nutrient in a simplified flow-through system. The nutrient flows in at a constant rate, and phytoplankton dies and is lost from the system at a constant rate:
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()
This works perfectly for the supplied time-step, but e.g. m.time = np.arange(0,10) returns a nonsensical solution (two divergent lines reaching >1e7). Odeint has no problem solving it:
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()
The models I am building with my package can get relatively complex, with hundreds of state variables, and a much larger number of interactions, yielding highly non-linear results. I am not sure how I can be sure that my supplied time step is appropriate, since smaller time steps significantly increase computational time.
Is there a solver included with GEKKO (or compatible with the GEKKO model syntax) that provides a similar solver to odeint with adaptive step size? Or is there another approach that is better suited to deal with ecological models based on ODEs (or spatially-discretized PDE systems)?
Any help is very much appreciated!
