This notebook contains course material from CBE40455 by Jeffrey Kantor (jeff at nd.edu); the content is available on Github. The text is released under the CC-BY-NC-ND-4.0 license, and code is released under the MIT license.
This Jupyter notebook demonstrates use of the CVXPY package for describing and solving convex optimization problems. CVXPY provides a modeling language embedded within Python that is well suited to linear and quadratic optimization problems with linear constraints. A limited set of solvers are distributed with CVXPY, and interfaces are available for others. Further documentation is available here (.pdf).
It is necessary to install CVXPY before use. Instructions can be found here. If you happen to be using the Anaconda distribution of Python, the following cell should be sufficient for installation.
!conda install -c cvxgrp -y cvxpy
CVXPY is already installed in Google Colaboratory, so no extra installation steps are required on that platform.
One of the simpliest use cases for cvxpy is solving a system of linear equations. For example, consider the pair of the equations
\begin{align*} 3\,x + 4\,y & = 26 \\ 2\,x - 3\,y & = -11 \end{align*}where $x$ and $y$ are the unknowns. Each equation is a linear equality constraint. The problem is easily solved with a few simple cvxpy constructss.
import numpy as np
import cvxpy as cp
# create the variables to solve for
x = cp.Variable()
y = cp.Variable()
# create a list of equality constraints
eqns = [
3*x + 4*y == 26,
2*x - 3*y == -11
]
# create a problem instance
prob = cp.Problem(cp.Minimize(0), eqns)
# solve
prob.solve()
print(prob.status)
# display solution
print(x.value)
print(y.value)
optimal 2.0000000000000004 4.999999999999999
import cvxpy as cp
import numpy as np
# specify problem data as numpy matrix and vectors
A = np.array([[3, 4], [2, -3]])
b = np.array([26, -11])
# create a vector variables
x = cp.Variable(2)
# specify equality constraints
constraints = [A*x == b]
# form objective.
objective = cp.Minimize(cp.norm(x))
# form and solve problem.
problem = cp.Problem(objective, constraints)
problem.solve()
# display the solution
print(" Problem Status:", problem.status)
print("Objective Value:", problem.value)
print(" Solution x[0]:", x[0].value)
print(" Solution x[1]:", x[1].value)
Problem Status: optimal Objective Value: 5.385164802505173 Solution x[0]: 2.0 Solution x[1]: 4.999999999999999
Get the historical record of ice out dates on Rainy Lake, Minnesota. Rainy Lake is a large lake, approximately 1000 sq. km., located on the Minnesota/Ontario border. Reliable records for ice-out date goes back to 1935, where for this lake ice-out is defined as the first day the lake can be navigated from one end to the other following a specified route.
The next cell downloads and extracts the historical ice-out date from Climatology office of the Minnesota Department of Natural Resources. The url needed to extract the historical information is found by clicking on a particular lake, then on the link labeled "all ice-out data for this lake".
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from urllib.request import urlopen
import json
from dateutil.parser import parse
import datetime
# Lake ID for Rainy Lake
lakeid = '69069400'
# Get Ice-out Data
url = 'https://maps1.dnr.state.mn.us/cgi-bin/climatology/ice_out_by_lake.cgi?id=' + lakeid
response = urlopen(url).read().decode('utf8')
data = json.loads(response)
print('Data Download Status: ' + data['status'])
# Calculate ice out 'day of year' for each date
dates = [parse(d['date']) for d in data['result']['values']]
year = np.array([d.year for d in dates])
iceout = np.array([1 + datetime.date.toordinal(d) - datetime.date(d.year,1,1).toordinal()
for d in dates ])
plt.plot(year, iceout)
plt.title('Ice Out Day on Rainy Lake, Minnesota')
plt.ylabel('Julian Day of Year')
plt.xlabel('Year');
Data Download Status: OK
Create a regression model and fit using cvxpy.
# variables are slope, intercept, and predicted values
slope = cp.Variable()
intercept = cp.Variable()
pred = cp.Variable(len(year))
con = [pred == (intercept + slope*(year-1930))]
obj = cp.Minimize(cp.norm(iceout-pred,2))
prob = cp.Problem(obj,con)
prob.solve()
print(" Slope: {0:8.4f}".format(slope.value))
print("Intercept: {0:8.4f}".format(intercept.value))
Slope: -0.0753 Intercept: 126.7995
Plot the resulting fit and compare to the original data.
plt.plot(year, iceout, year, pred.value)
plt.title('Ice Out Day on Rainy Lake, Minnesota')
plt.ylabel('Day of Year')
plt.xlabel('Year');
Let's consider the problem
\begin{align*} \text{Minimize}\quad & 4 x + 6y \\ \text{subject to}\quad & 2x + 2y \geq 5 \\ & x - y \leq 1 \\ & x,y \geq 0 \\ & x,y \in \text{Integer} \end{align*}Let's first solve this problem ignoring the integer constraints on $x$ and $y$. This is known as the 'relaxation' of an integer or binary programming problem.
import numpy as np
import cvxpy as cp
# create the variables to solve for
x = cp.Variable()
y = cp.Variable()
# create a list of equality constraints
constraints = [
2*x + 2*y >= 5,
x - y <= 1,
x >= 0,
y >= 0
]
obj = cp.Minimize(4*x + 6*y)
# create a problem instance
prob = cp.Problem(obj,constraints)
# solve
prob.solve()
print("status = ", prob.status)
# display solution
print("x = ", x.value)
print("y = ", y.value)
print("Obj = ", prob.value)
status = optimal x = 1.7499999999999998 y = 0.7499999999999998 Obj = 11.499999999999998
Since this is a minimization, the relaxed solution provides a lower bound on the solution of the integer problem. The next step is to define variables $x$ and $y$ as integer variables using cvx.Variable(integer=True)
function.
import numpy as np
import cvxpy as cp
# create the variables to solve for
x = cp.Variable(integer=True)
y = cp.Variable(integer=True)
# create a list of equality constraints
constraints = [
2*x + 2*y >= 5,
x - y <= 1,
x >= 0,
y >= 0
]
obj = cp.Minimize(4*x + 6*y)
# create a problem instance
prob = cp.Problem(obj, constraints)
# solve
prob.solve()
print("status = ", prob.status)
# display solution
print("x = ", x.value)
print("y = ", y.value)
print("Obj = ", prob.value)
status = optimal_inaccurate x = 1.9999999888557622 y = 1.0000000116464522 Obj = 14.000000025301762
The solvers distributed with CVXPY tyically employ interior-point algorithms augmented with a branch and bound strategy for integer or boolean variables. A consequence is that integer solutions are reported as floating point numbers that may need rounding for some applications.
print('Integer Solution x = ', int(np.round(x.value)))
print('Integer Solution y = ', int(np.round(y.value)))
Integer Solution x = 2 Integer Solution y = 1