In this minimal example I want to solve the basic integrator ODE dT/dt = k[F(t) - T(t)] where F(t) is a forcing vector, which is selected to be a square wave. I have implemented two procedures to solve the ODE: Using Scipy's solve_ivp and using Gekko. The former takes 6 milliseconds to run, the latter takes 8 seconds to run. Am I doing something wrong with Gekko? Are there some additional parameters that I can tune to improve performance?
I have omitted the testing part of the code, but I have compared both solutions to the analytic solution, and they are both accurate within a 3-4 significant digits.
import numpy as np
from scipy.integrate import solve_ivp
from gekko import GEKKO
def step_problem_forcing(time_arr, param):
T_low = param['T_low']
T_high = param['T_high']
t_break_l = param['t_break_l']
t_break_r = param['t_break_r']
forcing = np.full(len(time_arr), T_low)
idxs_pulse = np.logical_and(time_arr >= t_break_l,
time_arr <= t_break_r)
forcing[idxs_pulse] = T_high
return forcing
def step_problem_generate(tmin, tmax):
tAvg = (tmin+tmax) / 2
tDur = tmax - tmin
return {
'k' : 10 ** np.random.uniform(-3, 3),
'T_low' : np.random.uniform(20, 50),
'T_high' : np.random.uniform(20, 50),
'T0' : np.random.uniform(20, 50),
't_break_l' : np.random.uniform(tmin, tAvg - 0.1*tDur),
't_break_r' : np.random.uniform(tAvg + 0.1*tDur, tmax)
}
def ode_scipy_step_solver(time_arr, param: dict):
f_this = lambda t: np.interp(t, time_arr, param['forcing'])
def _dxdt(t, x, param: tuple):
# if (t < tmin) or (t > tmax):
# raise ValueError(f"{t} is out of bounds [{tmin}, {tmax}]")
k, = param
return k*(f_this(t) - x)
k = param['k']
sol = solve_ivp(_dxdt, (time_arr[0], time_arr[-1]),
[param['T0'], ], t_eval=time_arr, args=([k, ],),
rtol=1.0E-6, method='RK45')
return sol['y'].T
def ode_gekko_step_solver(time_arr, param: dict) -> np.ndarray:
m = GEKKO() # create GEKKO model
# Define variables
m.time = time_arr
T = m.Var(value=param['T0'])
F = m.Param(value=param['forcing'])
k = m.Const(value=param['k'])
# t = m.Param(value=m.time)
# equations
m.Equation(T.dt() == k * (F - T))
m.options.IMODE = 7 # dynamic simulation
m.solve(disp=False) # solve locally (remote=False)
return T.value
tmin = 0
tmax = 10
time_arr = np.linspace(tmin, tmax, 1000)
param = step_problem_generate(tmin, tmax)
param['forcing'] = step_problem_forcing(time_arr, param)
# Takes 6.8ms to run
rezScipy = ode_scipy_step_solver(time_arr, param)
# Takes 8.5s to run
rezGekko = ode_gekko_step_solver(time_arr, param)
Switch to remote=False to solve locally and avoid the overhead due to network communication to the public server. Switching to IMODE=4 solves all 1000 time steps simultaneously instead of sequentially and can improve solve time. Here are 10 trials with the 1 original and 9 modified versions. I'm on Ubuntu Linux. Timing will likely be different based on the computer.
# IMODE=7, Original Script
ODE Scipy Time: 0.014239072799682617
Gekko Time: 8.92684006690979
# IMODE=4, Remote=False
ODE Scipy Time: 0.2553427219390869
Gekko Time: 0.12177252769470215
ODE Scipy Time: 0.0017511844635009766
Gekko Time: 0.11760997772216797
ODE Scipy Time: 0.0032286643981933594
Gekko Time: 0.11528849601745605
ODE Scipy Time: 0.06327152252197266
Gekko Time: 0.12140154838562012
ODE Scipy Time: 0.01037144660949707
Gekko Time: 0.12077593803405762
ODE Scipy Time: 0.0036330223083496094
Gekko Time: 0.11592268943786621
ODE Scipy Time: 0.1100163459777832
Gekko Time: 0.11649894714355469
ODE Scipy Time: 0.0032815933227539062
Gekko Time: 0.1154778003692627
ODE Scipy Time: 0.0015087127685546875
Gekko Time: 0.1158149242401123
Here's the modified script:
import numpy as np
from scipy.integrate import solve_ivp
from gekko import GEKKO
import time
def step_problem_forcing(time_arr, param):
T_low = param['T_low']
T_high = param['T_high']
t_break_l = param['t_break_l']
t_break_r = param['t_break_r']
forcing = np.full(len(time_arr), T_low)
idxs_pulse = np.logical_and(time_arr >= t_break_l,
time_arr <= t_break_r)
forcing[idxs_pulse] = T_high
return forcing
def step_problem_generate(tmin, tmax):
tAvg = (tmin+tmax) / 2
tDur = tmax - tmin
return {
'k' : 10 ** np.random.uniform(-3, 3),
'T_low' : np.random.uniform(20, 50),
'T_high' : np.random.uniform(20, 50),
'T0' : np.random.uniform(20, 50),
't_break_l' : np.random.uniform(tmin, tAvg - 0.1*tDur),
't_break_r' : np.random.uniform(tAvg + 0.1*tDur, tmax)
}
def ode_scipy_step_solver(time_arr, param: dict):
f_this = lambda t: np.interp(t, time_arr, param['forcing'])
def _dxdt(t, x, param: tuple):
# if (t < tmin) or (t > tmax):
# raise ValueError(f"{t} is out of bounds [{tmin}, {tmax}]")
k, = param
return k*(f_this(t) - x)
k = param['k']
sol = solve_ivp(_dxdt, (time_arr[0], time_arr[-1]),
[param['T0'], ], t_eval=time_arr, args=([k, ],),
rtol=1.0E-6, method='RK45')
return sol['y'].T
def ode_gekko_step_solver(time_arr, param: dict) -> np.ndarray:
m = GEKKO(remote=False) # create GEKKO model
# Define variables
m.time = time_arr
T = m.Var(value=param['T0'])
F = m.Param(value=param['forcing'])
k = m.Const(value=param['k'])
# t = m.Param(value=m.time)
# equations
m.Equation(T.dt() == k * (F - T))
m.options.IMODE = 4 # dynamic simulation
m.solve(disp=False) # solve locally (remote=False)
return T.value
tmin = 0
tmax = 10
time_arr = np.linspace(tmin, tmax, 1000)
param = step_problem_generate(tmin, tmax)
param['forcing'] = step_problem_forcing(time_arr, param)
st = time.time()
rezScipy = ode_scipy_step_solver(time_arr, param)
print(f'ODE Scipy Time: {time.time()-st}')
st = time.time()
rezGekko = ode_gekko_step_solver(time_arr, param)
print(f'Gekko Time: {time.time()-st}')
Scipy solve_ivp is a specialized ODE solver while Gekko solves higher index DAEs, mixed integer optimization, regression, and machine learning models. It uses a simultaneous method and only has basic adaptive step lengths to control error. I recommend solve_ivp if it is an ODE problem.