GEKKO Optimization Suite¶
Miftahur Rahman, Ph.D.
March 31, 2022
Overview¶
GEKKO is a Python package for machine learning and optimization of mixed-integer and differential algebraic equations. It is coupled with large-scale solvers for linear, quadratic, nonlinear, and mixed integer programming (LP, QP, NLP, MILP, MINLP). Modes of operation include parameter regression, data reconciliation, real-time optimization, dynamic simulation, and nonlinear predictive control. GEKKO is an object-oriented Python library to facilitate local execution of APMonitor.
More of the backend details are available at What does GEKKO do? and in the GEKKO Journal Article. Example applications are available to get started with GEKKO.
Installation¶
A pip package is available (see current download stats):
pip install gekko
In [1]:
from gekko import GEKKO
m = GEKKO() # create GEKKO model
x = m.Var() # define new variable, default=0
y = m.Var() # define new variable, default=0
m.Equations([5*x+6*y==1, 3*x+4*y==0]) # equations
m.solve(disp=False) # solve
print(x.value,y.value) # print solution
In [2]:
from gekko import GEKKO
m = GEKKO() # create GEKKO model
x = m.Var(value=0) # define new variable, initial value=0
y = m.Var(value=1) # define new variable, initial value=1
m.Equations([5*x + 6*y==0, 3*x**2+4*y**2==1]) # equations
m.solve(disp=False) # solve
print([x.value[0],y.value[0]]) # print solution
Interpolation with Cubic Spline¶
In [10]:
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
xm = np.array([0,1,2,3,4,5])
ym = np.array([0.1,0.2,0.5,0.7,1.5,0.7])
m = GEKKO() # create GEKKO model
m.options.IMODE = 2 # solution mode
x = m.Param(value=np.linspace(-1,6)) # prediction points
y = m.Var() # prediction results
m.cspline(x, y, xm, ym) # cubic spline
m.solve(disp=False) # solve
# create plot
plt.plot(xm,ym,'bo')
plt.plot(x.value,y.value,'r--',label='cubic spline')
plt.legend(loc='best')
Out[10]:
Linear and Polynomial Regression¶
In [4]:
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
xm = np.array([0,1,2,3,4,5])
ym = np.array([0.2,0.3,0.5,0.7,0.9,3.0])
#### Solution
m = GEKKO()
m.options.IMODE=2
# coefficients
c = [m.FV(value=0) for i in range(4)]
x = m.Param(value=xm)
y = m.CV(value=ym)
y.FSTATUS = 1
# polynomial model
m.Equation(y==c[0]+c[1]*x+c[2]*x**2+c[3]*x**3)
# linear regression
c[0].STATUS=1
c[1].STATUS=1
m.solve(disp=False)
p1 = [c[1].value[0],c[0].value[0]]
# quadratic
c[2].STATUS=1
m.solve(disp=False)
p2 = [c[2].value[0],c[1].value[0],c[0].value[0]]
# cubic
c[3].STATUS=1
m.solve(disp=False)
p3 = [c[3].value[0],c[2].value[0],c[1].value[0],c[0].value[0]]
# plot fit
plt.plot(xm,ym,'ko',markersize=10)
xp = np.linspace(0,5,100)
plt.plot(xp,np.polyval(p1,xp),'b--',linewidth=2)
plt.plot(xp,np.polyval(p2,xp),'r--',linewidth=3)
plt.plot(xp,np.polyval(p3,xp),'g:',linewidth=2)
plt.legend(['Data','Linear','Quadratic','Cubic'],loc='best')
plt.xlabel('x')
plt.ylabel('y')
Out[4]:
Nonlinear Regression¶
In [5]:
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
# measurements
xm = np.array([0,1,2,3,4,5])
ym = np.array([0.1,0.2,0.3,0.7,0.9,2.5])
# GEKKO model
m = GEKKO()
# parameters
x = m.Param(value=xm)
a = m.FV()
a.STATUS=1
# variables
y = m.CV(value=ym)
y.FSTATUS=1
# regression equation
m.Equation(y==0.1*m.exp(a*x))
# regression mode
m.options.IMODE = 2
# optimize
m.solve(disp=False)
# print parameters
print('Optimized, a = ' + str(a.value[0]))
plt.plot(xm,ym,'bo')
plt.plot(xm,y.value,'r-')
Out[5]:
Solve Differential Equation(s)¶
Solve the following differential equation with initial condition y(0)=5:
where k=15. The solution of y(t) should be reported from an initial time 0 to final time 20. Create a plot of the result for y(t) versus t.:
In [6]:
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO()
m.time = np.linspace(0,20,100)
k = 15
y = m.Var(value=5.0)
t = m.Param(value=m.time)
m.Equation(k*y.dt()==-t*y)
m.options.IMODE = 4
m.solve(disp=False)
plt.plot(m.time,y.value)
plt.xlabel('time')
plt.ylabel('y')
Out[6]:
Mixed Integer Nonlinear Programming¶
In [7]:
from gekko import GEKKO
m = GEKKO() # Initialize gekko
m.options.SOLVER=1 # APOPT is an MINLP solver
# optional solver settings with APOPT
m.solver_options = ['minlp_maximum_iterations 500', \
# minlp iterations with integer solution
'minlp_max_iter_with_int_sol 10', \
# treat minlp as nlp
'minlp_as_nlp 0', \
# nlp sub-problem max iterations
'nlp_maximum_iterations 50', \
# 1 = depth first, 2 = breadth first
'minlp_branch_method 1', \
# maximum deviation from whole number
'minlp_integer_tol 0.05', \
# covergence tolerance
'minlp_gap_tol 0.01']
# Initialize variables
x1 = m.Var(value=1,lb=1,ub=5)
x2 = m.Var(value=5,lb=1,ub=5)
# Integer constraints for x3 and x4
x3 = m.Var(value=5,lb=1,ub=5,integer=True)
x4 = m.Var(value=1,lb=1,ub=5,integer=True)
# Equations
m.Equation(x1*x2*x3*x4>=25)
m.Equation(x1**2+x2**2+x3**2+x4**2==40)
m.Obj(x1*x4*(x1+x2+x3)+x3) # Objective
m.solve(disp=False) # Solve
print('Results')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('x4: ' + str(x4.value))
print('Objective: ' + str(m.options.objfcnval))
Model Predictive Control¶
In [8]:
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO()
m.time = np.linspace(0,20,41)
# Parameters
mass = 500
b = m.Param(value=50)
K = m.Param(value=0.8)
# Manipulated variable
p = m.MV(value=0, lb=0, ub=100)
p.STATUS = 1 # allow optimizer to change
p.DCOST = 0.1 # smooth out gas pedal movement
p.DMAX = 20 # slow down change of gas pedal
# Controlled Variable
v = m.CV(value=0)
v.STATUS = 1 # add the SP to the objective
m.options.CV_TYPE = 2 # squared error
v.SP = 40 # set point
v.TR_INIT = 1 # set point trajectory
v.TAU = 5 # time constant of trajectory
# Process model
m.Equation(mass*v.dt() == -v*b + K*b*p)
m.options.IMODE = 6 # control
m.solve(disp=False)
# get additional solution information
import json
with open(m.path+'//results.json') as f:
results = json.load(f)
plt.figure()
plt.subplot(2,1,1)
plt.plot(m.time,p.value,'b-',label='MV Optimized')
plt.legend()
plt.ylabel('Input')
plt.subplot(2,1,2)
plt.plot(m.time,results['v1.tr'],'k-',label='Reference Trajectory')
plt.plot(m.time,v.value,'r--',label='CV Response')
plt.ylabel('Output')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()
Optimization of Multiple Linked Phases¶
In [9]:
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
# Initialize gekko model
m = GEKKO()
# Number of collocation nodes
nodes = 5
# Number of phases
n = 4
# Time horizon (for all phases)
m.time = np.linspace(0,1,100)
# Input (constant in IMODE 4)
u = [m.Var(1,lb=-2,ub=2,fixed_initial=False) for i in range(n)]
# Example of same parameter for each phase
tau = 5
# Example of different parameters for each phase
K = [2,3,5,1,4]
# Scale time of each phase
tf = [1,2,4,8,16]
# Variables (one version of x for each phase)
x = [m.Var(0) for i in range(5)]
# Equations (different for each phase)
for i in range(n):
m.Equation(tau*x[i].dt()/tf[i]==-x[i]+K[i]*u[i])
# Connect phases together at endpoints
for i in range(n-1):
m.Connection(x[i+1],x[i],1,len(m.time)-1,1,nodes)
m.Connection(x[i+1],'CALCULATED',pos1=1,node1=1)
# Objective
# Maximize final x while keeping third phase = -1
m.Obj(-x[n-1]+(x[2]+1)**2*100)
# Solver options
m.options.IMODE = 6
m.options.NODES = nodes
# Solve
m.solve()
# Calculate the start time of each phase
ts = [0]
for i in range(n-1):
ts.append(ts[i] + tf[i])
# Plot
plt.figure()
tm = np.empty(len(m.time))
for i in range(n):
tm = m.time * tf[i] + ts[i]
plt.plot(tm,x[i])
