# # Getting Started with CVXPY
# This [Jupyter notebook](http://jupyter.org/notebook.html) demonstrates use of the [CVXPY](http://www.cvxpy.org/en/latest/) 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)](https://media.readthedocs.org/pdf/cvxpy/latest/cvxpy.pdf).
# ## Installation
#
# It is necessary to install CVXPY before use. Instructions can be found [here](http://www.cvxpy.org/en/latest/install/). 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.
# ## Solving linear equations
#
# 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.
# In[ ]:
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)
# In[ ]:
# ### Using Matrix Parameters and Vector Variables
# In[9]:
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)
# ## Regression: Evidence of Climate Change
#
# 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](https://www.dnr.state.mn.us/ice_out/index.html). 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".
# In[ ]:
get_ipython().run_line_magic('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');
# Create a regression model and fit using cvxpy.
# In[ ]:
# 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))
# Plot the resulting fit and compare to the original data.
# In[ ]:
plt.plot(year, iceout, year, pred.value)
plt.title('Ice Out Day on Rainy Lake, Minnesota')
plt.ylabel('Day of Year')
plt.xlabel('Year');
# ## Integer Program
#
# 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.
# In[ ]:
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)
# 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.
# In[ ]:
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)
# 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.
# In[ ]:
print('Integer Solution x = ', int(np.round(x.value)))
print('Integer Solution y = ', int(np.round(y.value)))
#
# < [Getting Started Guides](http://nbviewer.jupyter.org/github/jckantor/CBE40455/blob/master/notebooks/01.00-Getting-Started-Guides.ipynb) | [Contents](toc.ipynb) | [Getting Started with Gurobi](http://nbviewer.jupyter.org/github/jckantor/CBE40455/blob/master/notebooks/01.02-Getting-Started-with-Gurobi.ipynb) >