Dynamic simulation of biological system with Gekko Python not working

Question

I'm trying to simulate a basic biological system: growth of activated sludge in a reactor fed by wastewater.

It has two variables and two reactions. This is a very basic model and used extensively (more elaborate) in wastewater engineering.

The solver seems very unstable: frequently doesn't find a solution of goes haywire (large overshoots). The model is very smooth, no rapid changes after an initial little jump.

I tried the IMODE = 4 (dynamic simultaneous) but fails most of the time. Currently running with IMODE = 7 (dynamic sequential). Works for time interval from 0.1 - 4 days, but fails with larger period (5 days and up). This while the system is more or less stable after 0.2 days (see below).

The solver also takes quiet a long time to solve. While I was hoping to use it for much larger work.

My questions: Is there something wrong with my code, or is it the solver? Can the code be improved to speed up the simulation?

This test ran in a Jupyter notebook (chrome) with following packages installed: 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')

Answer 1

Gekko's adaptive time-step refinement for IMODE=7 is still under development (track feature request on GitHub). The solution fails when the solver thinks that the biomass concentration should be negative. There is exponential aerobic growth at the very beginning and very slow aerobic decay. If you only need simulation, I recommend using an adaptive time-step solver such as scipy.integrate.ODEINT (tutorials). However, if you need to perform optimization such as regression or batch control then Gekko is a good option. To solve successfully with either IMODE=4 or IMODE=7 here are some tips:

1. Add constraints to prevent biomass or substrate from going negative.

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

2. Set a time-step of 1 with a final time that is adjustable by changing n. Add time points 0.001, 0.002, ... , 0.8 to have an accurate solution with the fast dynamics of aerobic growth. When the adaptive time-step refinement is completed, inserting the extra time-steps won't be needed.

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

3. Solve locally instead of using the default server.

m = GEKKO(remote=False)

Here is the complete script that solves successfully with any final time and with IMODE=4 or 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, stop=n, num=n+1)
m.time = np.insert(t,1,[0.001,0.002,0.004,0.008,0.02,0.04,0.08,\
                   0.2,0.4,0.8])

# create state variables
S = m.Var(value=100,lb=0)  # substrate 
X = m.Var(value=500,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()