import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

m = GEKKO(remote=True)

# Add 0.01 as first step
# 0,0.01,0.1,0.2,0.3,...11.9,12.0)
m.time = np.insert(np.linspace(0,12,121),1,0.01)

# change solver options
m.solver_options = ['minlp_gap_tol 0.001',\
                    'minlp_maximum_iterations 10000',\
                    'minlp_max_iter_with_int_sol 30',\
                    'minlp_branch_method 3',\
                    'minlp_integer_tol 0.001',\
                    'minlp_integer_leaves 0',\
                    'minlp_maximum_iterations 200']

c0 = 0.4; c1 = 0.2

last = m.Param(np.zeros(122))
last.value[-1] = 1

x0 = m.CV(value=0.5,lb=0)
x0.STATUS = 1; x0.SPHI=1.05; x0.SPLO=0.95
x1 = m.CV(value=0.7,lb=0)
x1.STATUS = 1; x1.SPHI=1.05; x1.SPLO=0.95

w = m.MV(value=0,lb=0,ub=1,integer=True)
w.STATUS = 1

m.Equations([x0.dt() == x0 - x0*x1 - c0*x0*w,\
             x1.dt() == - x1 + x0*x1 - c1*x1*w])

m.options.IMODE = 6
m.options.NODES = 2
m.options.SOLVER = 1
m.options.MV_TYPE = 0
m.solve()

plt.figure(figsize=(6,4))
plt.subplot(2,1,1)
plt.step(m.time,w.value,'r-',label='w')
plt.legend(loc=1); plt.grid(); plt.ylabel('Fleet Deployed')
plt.subplot(2,1,2)
plt.plot([0,12],[1.05,1.05],'g:',label=r'$SP_{HI}$')
plt.plot([0,12],[0.95,0.95],'g:',label=r'$SP_{LO}$')
plt.plot(m.time,x0.value,'b-',label=r'Prey $(x_0)$')
plt.plot(m.time,x1.value,'k-',label=r'Predator $(x_1)$')
plt.xlabel('Time'); plt.ylabel('Fish')
plt.legend(loc=1); plt.grid(); plt.tight_layout()
plt.savefig('lotka2.png',dpi=300); plt.show()