Practical Optimization for Stats Nerds¶
Ryan J. O'Neil
Data Science DC
March 20, 2017
ryanjoneil@gmail.com
https://ryanjoneil.github.io
@ryanjoneil
What to expect¶
- We'll take familiar models, show how they work from an optimization perspective, then show applications.
- This talk is (mostly) about model formulation.
- It might get a little technical, so...
Take-aways¶
- Many statistical techniques are based on some sort of optimization.
- Optimization has lots of uses, such as solving decision models.
- Learning to structure problems you already know for optimization solvers is a great way to understand them!
Model 1: Least Squares¶
We observe noisy data from an unknown function. We want to infer that function.
# Generate some random data.
import numpy as np
import random
# Sort the data so they're easier to plot later.
x = [random.uniform(-10, 10) for _ in range(500)]
x.sort()
y = []
for xi in x:
eps = random.normalvariate(0, 25)
yi = 3*xi**2 - 2*xi + 10 + eps
y.append(yi)
x = np.array(x)
y = np.array(y)
from bokeh.charts import Scatter, output_notebook, show
output_notebook()
scatter = Scatter({'x': x, 'y': y}, width=750, height=400)
show(scatter)
Least Squares the way you do it...¶
...assuming you use scikit-learn like every other sane Python programmer.
We looked at a chart of our data and decided to describe it with:
- A quadratic term
- A linear term
- An offset
from sklearn.linear_model import LinearRegression
# Note: A is our feature matrix.
# We intentionally add a "1" for the offset, instead of letting
# sklearn do that for us. This will make sense soon.
X = np.array([[xi**2, xi, 1] for xi in x])
lin = LinearRegression(fit_intercept=False)
lin.fit(X, y)
print(lin.coef_)
How'd we do?
from bokeh.charts import Line
y_hat = lin.predict(X)
show(Line({'x': x, 'y': y_hat}, x='x', y='y', width=750, height=400))
Least Squares the way your grandparents did it...¶
...with chalk and a slab of slate.
Construct a function to calculate the sum of squared residuals...
$$ \begin{align} \text{min}\ f(\beta) & = \frac{1}{2} ||y - X \beta||^2 \\ & \\ & = \frac{1}{2} (y - X \beta)'(y - X \beta) \\ & \\ & = \frac{1}{2} y'y - y'X\beta + \frac{1}{2} \beta'X'X\beta \\ & \\ & = \frac{1}{2} y'y - y'X\beta + \frac{1}{2} \beta'X'X\beta \\ & \\ \end{align} $$First Order Necessary Conditions¶
...and take its derivative to find a closed-form solution.
$$\nabla f(\beta) = X'X\beta - y'X\beta = 0$$
$$\beta = (X'X)^{-1}X'y$$
from numpy.linalg import inv
# beta = (X'X)^-1 * X * y
Xt = X.transpose()
pseudo_inv = inv(np.matmul(Xt, X))
beta = np.matmul(np.matmul(pseudo_inv, Xt), y)
print(beta)
How'd grandma and grandpa do?
y_hat = [beta[0]*xi**2 + beta[1]*xi + beta[2] for xi in x]
show(Line({'x': x, 'y': y_hat}, x='x', y='y', width=750, height=400))
Least Squares the way your crazy uncle Eddie does it...¶
...'cause he used to work at NASA and code in Forth.
So we need to convert from
$$\frac{1}{2} \beta'X'X\beta - y'X\beta + \frac{1}{2} y'y $$to another form
$$\frac{1}{2} \beta'P\beta + q'\beta$$which is simply
$$P = X'X, q = -y'X$$import cvxopt as cvx
P = cvx.matrix(np.matmul(Xt, X))
q = cvx.matrix(-1 * np.matmul(y.transpose(), X))
solution = cvx.solvers.qp(P, q)
beta = solution['x'] # unrelated to our x
print(beta)
How'd Crazy Uncle Eddie do?
y_hat = [beta[0]*xi**2 + beta[1]*xi + beta[2] for xi in x]
show(Line({'x': x, 'y': y_hat}, x='x', y='y', width=750, height=400))
So what's different about Crazy Uncle Eddie?¶
Well, besides the obvious.
While all three techniques produced the same result, Crazy Uncle Eddie's is interesting because it is more general than the others.
Crazy Eddie can solve any quadratic optimization problem, of which Least Squares is just one instance.
If we change the structure of the problem slightly:
- We probably can't solve it with
scikit-learn - Grams and Gramps have to go back to their chalkboard
- Crazy Eddie can update the inputs to his problem and reoptimize
Application 1: Portfolio Optimization¶
We have a big pot of money to allocate among different investments. Lucky us!
Some investment returns are correlated. They go up and down together.
Other returns are anticorrelated. They tend to do the opposite things.
How do we allocate our money to maximize our expected return, subject to our tolerance for risk?
We'll use 100 months of exchange rate data from the Fed, circa 2014.
Since this talk isn't about data wrangling, I've already cleaned it up into the important pieces.
- Expected monthly return data for each foreign currency
- A covariance matrix for those investments
Model:¶
$x$ tells me how much of my total budget to put in each investment.
$$ \begin{align} \text{max} \ \ \ & \ \mu'x - \alpha x \Sigma x \\ \text{s.t.} \ \ \ & \ e'x = 1 \\ & \ x \ge 0 \end{align} $$But wait! This is a maximization problem! The last model used $\text{min}$.
$$ \begin{align} \text{min} \ \ \ & \ \alpha x \Sigma x - \mu'x \\ \text{s.t.} \ \ \ & \ e'x = 1 \\ & \ x \ge 0 \end{align} $$The only differences between this model and the least squares model are the constraints we've added.
This one forces the model to allocate all of my budget into investments:
$$e'x = 1$$This one disallows the model from making negative investments:
$$x \ge 0$$# Read in the returns and covariance data.
import pandas as pd
exp_returns = pd.read_csv('portfolio-optimization/currency-returns.csv')
exp_returns.head()
returns_cov = pd.read_csv('portfolio-optimization/currency-covariance.csv', header=None)
returns_cov.head()
# A model that will return an optimal portfolio for any risk aversion.
def portfolio(alpha):
P = cvx.matrix(alpha * returns_cov.as_matrix())
q = cvx.matrix(-exp_returns['mean'].as_matrix())
G = cvx.matrix(0.0, (len(q),len(q)))
G[::len(q)+1] = -1.0
h = cvx.matrix(0.0, (len(q),1))
A = cvx.matrix(1.0, (1,len(q)))
b = cvx.matrix(1.0)
solution = cvx.solvers.qp(P, q, G, h, A, b)
return exp_returns['mean'].dot(solution['x'])[0]
risk_aversion = [ra/2.0 for ra in range(41)]
returns = [portfolio(alpha) for alpha in risk_aversion]
show(Line(
{'risk aversion': risk_aversion, 'expected return': returns},
x='risk aversion',
y='expected return',
width=750,
height=400
))
Model 2: Support Vector Machines¶
We want to draw a line that separates two sets with as little misclassification as possible.
If $f(x) \le 5$, the point is of type $-1$, otherwise it is type $+1$.
$$f(x) = (x_1 + \epsilon_1) - (x_2+ \epsilon_2)$$ $$\epsilon_i \sim N\left(0, 1.25^2\right)\ \forall\ i \in {1, 2}$$x1 = [random.uniform(0, 10) for _ in range(150)]
x2 = [random.uniform(0, 10) for _ in range(150)]
y = []
for xi_1, xi_2 in zip(x1, x2):
eps_1 = random.normalvariate(0, 1.25)
eps_2 = random.normalvariate(0, 1.25)
if (xi_1 + eps_1) - (xi_2 + eps_2) <= 5:
y.append(-1)
else:
y.append(+1)
X = np.array(list(zip(x1, x2)))
y = np.array(y)
show(Scatter(
{'x1': x1, 'x2': x2, 'y': y},
x='x1', y='x2', color='y', marker='y', palette=['lightblue', 'orange'],
width=500, height=500
))
Support Vector Machines a la scikit-learn¶
from sklearn import svm
clf = svm.SVC(kernel='linear')
clf.fit(X, y)
The SVM coefficients give us a separating hyperplane with minimum classification error.
$$b + w'x = 0$$print('b =', clf.intercept_[0])
print('w =', clf.coef_[0])
We classify a point $x = (x_1, x_2)'$ based on which side of that hyperplane it is on.
$$b + w'x \ge 0\ \rightarrow\ \text{Type +1}$$
$$b + w'x < 0\ \rightarrow\ \text{Type -1}$$
from sklearn import metrics
y_hat = clf.predict(X)
print('confusion matrix:')
metrics.confusion_matrix(y_true=y, y_pred=y_hat)
# We need a function to plot classifications and the separating hyperplane.
def plot_svm_results(correct, b_val, w1_val, w2_val):
# Classifications
plot = Scatter(
{'x1': x1, 'x2': x2, 'y': y, 'correct': correct},
x='x1', y='x2',
color='correct', palette=['lightgrey', 'darkred'], marker='correct',
width=500, height=500
)
# Separating hyperplane: x2 = -(b + (w1*x1)) / w2
hyperplane_x1 = np.linspace(0, 10)
hyperplane_x2 = [
-(b_val + (w1_val * x1_i)) / w2_val for x1_i in hyperplane_x1
]
plot.line(hyperplane_x1, hyperplane_x2, color='darkblue')
show(plot)
How'd scikit-learn do?
correct = ['Correct' if yi == yh_i else 'Incorrect' for yi, yh_i in zip(y, y_hat)]
plot_svm_results(correct, clf.intercept_[0], clf.coef_[0][0], clf.coef_[0][1])
Support Vector Machines a la Linear Optimization¶
PuLP is a modeling tool that lets us solve anything of this form.
Where $c$ and $b$ are input vectors, $A$ is an input matrix, and $x$ is a vector of real-valued decision variables.
Linear Optimization is well solved. You can assume that a linear optimizer is robust and fast for any model in that form.
PuLP can also solve discrete models, which are not so computationally easy.
How do we get our SVM problem into a linear optimizer?
Inputs:¶
$$(x_{i1}, x_{i2})' = \text{column vector of inputs for observation}\ i$$ $$y_i \in \{+1, -1\} = \text{known classification of observation}\ i $$Model:¶
The decision variables $b$, $w_1$, and $w_2$, give us a separating hyperplane $b + w'x = 0$.
For each misclassified observation $i$, $\xi_i$ is its classification error. We find the best hyperplane by minimizing total error.
$$ \begin{align} \text{min} \ \ \ & \ \sum_i \xi_i \ & \ \\ \text{s.t.} \ \ \ & \ b + w_1 x_{i1} + w_2 x_{i2} \le -1 + \xi_i\ & \ \forall\ y_i = -1 \\ & \ b + w_1 x_{i1} + w_2 x_{i2} \ge +1 - \xi_i\ & \ \forall\ y_i = +1 \\ & \ \xi \ge 0 \end{align} $$We can simplify the model by multiplying each constraint by $y_i$.
$$ \begin{align} \text{min} \ \ \ & \ \sum_i \xi_i \\ \text{s.t.} \ \ \ & \ y_i (b + w_1 x_{i1} + w_2 x_{i2}) \ge 1 - \xi_i\ \forall\ i \\ & \ \xi \ge 0 \end{align} $$But wait! Constraints like $a'x \ge b$ aren't in the form you gave us!
That's OK. They get converted by the optimizer.
$$a'x - e = b$$$$e \ge 0$$So do constraints of the form $a'x \le b$.
$$a'x + s = b$$$$s \ge 0$$As do any variables that are unrestricted in sign.
$$x_i = x_i^+ - x_i^-$$$$x_i^+, x_i^- \ge 0$$As long as your objective function is linear and your constraints are linear and of the form $=$, $\le$, or $\ge$, you can solve it with a linear optimizer.
Let's get hacking. First make variables for the separating hyperplane.
Note these are unrestricted in sign.
import pulp
b = pulp.LpVariable('b')
w1 = pulp.LpVariable('w1')
w2 = pulp.LpVariable('w2')
We also need variables for the classification errors.
These are nonnegative since we only want to penalize errors, not reward the model for correct classifiations.
errors = []
for i in range(len(y)):
e = pulp.LpVariable('e_%d' % i, lowBound=0) # lowBound=0 means >=0
errors.append(e)
Build a model that records classification error based on the separating hyperplane in the appropriate variables.
m = pulp.LpProblem(sense=pulp.LpMinimize)
for x_i1, x_i2, y_i, e_i in zip(x1, x2, y, errors):
m += y_i * (b + w1 * x_i1 + w2 * x_i2) >= (1 - e_i)
So far the model represents the feasible region of all possible separating hyperplanes.
We want the best one.
m.setObjective(sum(errors))
assert m.solve() == pulp.LpStatusOptimal
print('b =', b.value())
print('w =', [w1.value(), w2.value()])
We don't have a predict function, so let's build one.
def classify(xi_1, xi_2):
if b.value() + (w1.value() * xi_1) + (w2.value() * xi_2) >= 0:
return +1
return -1
y_hat = [classify(xi_1, xi_2) for xi_1, xi_2 in zip(x1, x2)]
print('confusion matrix:')
metrics.confusion_matrix(y_true=y, y_pred=y_hat)
How'd the linear optimizer do?
correct = ['Correct' if yi == yh_i else 'Incorrect' for yi, yh_i in zip(y, y_hat)]
plot_svm_results(correct, b.value(), w1.value(), w2.value())
Interlude: Problem Shapes¶
In Linear Algebra 101 we learned to solve problems that look like this.
$$Ax = b$$$$ \begin{bmatrix} \cdot & \cdot \\ \cdot & \cdot \end{bmatrix} \begin{bmatrix} \cdot \\ \cdot \end{bmatrix} = \begin{bmatrix} \cdot \\ \cdot \end{bmatrix} $$Many statistical inference problems are overdetermined and unsolvable directly.
$$ \begin{bmatrix} \cdot & \cdot \\ \cdot & \cdot \\ \cdot & \cdot \\ \cdot & \cdot \\ \cdot & \cdot \\ \cdot & \cdot \\ \cdot & \cdot \\ \cdot & \cdot \\ \cdot & \cdot \\ \cdot & \cdot \end{bmatrix} \begin{bmatrix} \cdot \\ \cdot \end{bmatrix} = \begin{bmatrix} \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \end{bmatrix} $$Decision models tend to involve a large number of potential actions in relation to the constraint set.
$$ \begin{bmatrix} \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{bmatrix} \begin{bmatrix} \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \end{bmatrix} = \begin{bmatrix} \cdot \\ \cdot \end{bmatrix} $$In both cases, we start with a problem we can't solve easily and reformulate it as something we can.
For statistical inference, we solve a secondary problem that minimizes error against the observed data.
For decision models, we choose a set of actions to take (i.e. a basis) and solve the remaining system.
Model 3: Clustering¶
We have a set of unclassified data and we want to cluster it around logical centers.
We generate 25 points around each center $(x_{10}, x_{20}) \in \{(0.3, 0.3), (0.3, 0.7), (0.7, 0.7)\}$.
$$x_{1i} = x_{10} + \epsilon$$ $$x_{2i} = x_{20} + \epsilon$$$$\epsilon \sim N\left(0, 0.2^2\right)$$centers = [(0.3, 0.3), (0.3, 0.7), (0.7, 0.7)]
x1 = []
x2 = []
for x10, x20 in centers:
for _ in range(7):
eps_x1 = random.normalvariate(0, 0.2)
eps_x2 = random.normalvariate(0, 0.2)
x1i = min(max(0.0, x10+eps_x1), 1.0)
x2i = min(max(0.0, x20+eps_x2), 1.0)
x1.append(x1i)
x2.append(x2i)
X = np.array(list(zip(x1, x2)))
show(Scatter({'x1': x1, 'x2': x2}, x='x1', y='x2', width=500, height=500))
Clustering a la scikit-learn¶
Assume we know a prioi how many clusters there are.
from sklearn.cluster import KMeans
clust = KMeans(n_clusters=3)
clust.fit(X)
The KMeans model has two outputs:
- The center point of each cluster
- Which cluster each point is assigned to
clust.cluster_centers_
We need a function to easily plot clustering results.
def plot_clusters(clusters, cluster_centers):
plot = Scatter(
{'x1': x1, 'x2': x2, 'clusters': clusters},
x='x1', y='x2',
color='clusters', marker='clusters',
width=500, height=500
)
for x10, x20 in cluster_centers:
plot.circle(x10, x20, alpha=0.5, color='grey', size=15)
show(plot)
How'd scikit-learn do?
plot_clusters(clust.predict(X), clust.cluster_centers_)
Clustering a la Mixed Integer Linear Optimization¶
Let's extend the linar optimization model support the concept of binary variables.
$$ \begin{align} \text{min} \ \ \ & \ c'x + d'y \\ \text{s.t.} \ \ \ & \ A x + D y = b \\ & \ x \ge 0 \\ & \ y \in \{0, 1\} \end{align} $$$y$ is a vector of variables which can only take values $0$ or $1$. These represent discrete decisions and allow logical relationships among them.
Discrete optimization is very powerful, but not well solved the way linear optimization is. Model solvability is highly dependent on formulation.
What's a clustering model look like with binary variables?
Inputs:¶
$$(x_{1i}, x_{2i}) = \text{location of observed point} \in [0, 1]\ \forall\ i$$Model:¶
We need to make two major types of decisions:
- Where are the cluster centers, $\{w_1, w_2, w_3\}$?
- Which cluster is each point assigned to? $v_{ij} = 1$ if point $i$ is part of cluster $j$, $0$ if not.
We minimize total taxicab distance from each point to its assigned cluster
- $d_{x_{1}i}$ is the distance from point $i$ to its cluster center along the $x_1$ axis
- $d_{x_{2}i}$ is the distance from point $i$ to its cluster center along the $x_2$ axis
Primary decision variables = locations of cluster centers
w = []
for j in range(len(centers)):
x10 = pulp.LpVariable('x_10%d' % j, lowBound=0, upBound=1)
x20 = pulp.LpVariable('x_20%d' % j, lowBound=0, upBound=1)
w.append([x10, x20])
Secondary decision variables = do we associate point $i$ with cluster $j$?
v = []
for i in range(len(x1)):
vi = []
for j in range(len(w)):
vij = pulp.LpVariable('v_%d%d' % (i, j), cat=pulp.LpBinary)
vi.append(vij)
v.append(vi)
Each point is associated with exactly one cluster.
m = pulp.LpProblem(sense=pulp.LpMinimize)
for vi in v:
m += sum(vi) == 1
Make a list of taxicab distance variables to guide the optimizer.
distances = []
for i, (vi, x1i, x2i), in enumerate(zip(v, x1, x2)):
for j, (wj, vij) in enumerate(zip(w, vi)):
x10, x20 = wj
dx1 = pulp.LpVariable('dx_1%d%d' % (i, j), lowBound=0, upBound=1)
dx2 = pulp.LpVariable('dx_2%d%d' % (i, j), lowBound=0, upBound=1)
m += dx1 >= (x1i - x10) - (1 - vij)
m += dx1 >= (x10 - x1i) - (1 - vij)
m += dx2 >= (x2i - x20) - (1 - vij)
m += dx2 >= (x20 - x2i) - (1 - vij)
distances.extend([dx1, dx2])
WARNING WARNING WARNING
Symmetry elimination is beyond the scope of this talk.
We don't technically need it to solve this model, but it sure can help.
for j, (x01_1, x02_1) in enumerate(w):
for x01_2, x02_2 in w[j+1:]:
m += x01_1 + x02_1 <= x01_2 + x02_2
Total distance from points to their cluster centers is our objective function.
m.setObjective(sum(distances))
assert m.solve() == pulp.LpStatusOptimal
cluster_centers = []
for wj in w:
cluster_centers.append([wj[0].value(), wj[1].value()])
for c in cluster_centers:
print(c)
There's no predict method, so we yank that information out of the $v_{ij}$ variables.
clusters = []
for vi in v:
for c, vij in zip([0, 1, 2], vi):
if vij.value() > 0.5:
clusters.append(c)
break
How'd mixed integer optimization do?
plot_clusters(clusters, cluster_centers)
Hmmm. That model took a while to solve. What gives?
- When we introduce discrete choices into models, we lose the characteristics that make them easy to solve.
- This matter a lot less now that it did in the 90s.
- But effective formulation of discrete models is still (may always be?) an art.
Let's look at something more realistic.
Application 2: Pedicab Sharing¶
Peuber is the Uber of pedicabs.
We have 7 pedicabs in our fleet and 11 customer requesting service. It's a prime day for Peuber!
pedicabs = range(7)
customers = range(11)
Generate locations for each of these.
rng = lambda: random.uniform(0, 100)
pedicab_locations = [(rng(), rng()) for _ in pedicabs]
pickup_locations = [(rng(), rng()) for _ in customers]
dropoff_locations = [(rng(), rng()) for _ in customers]
This function will help us visualize the routing problem and its solutions.
from bokeh.models import Range1d
def build_pedicabs_and_customers_plot():
all_points = pedicab_locations + pickup_locations + dropoff_locations
point_types = ['Pedicab'] * len(pedicabs) + \
['Pickup'] * len(customers) + \
['Dropoff'] * len(customers)
plot = Scatter(
{
'x': [l[0] for l in all_points],
'y': [l[1] for l in all_points],
'type': point_types
},
x='x', y='y',
color='type', marker='type',
width=600, height=600
)
plot.legend.orientation = "horizontal"
plot.x_range = plot.y_range = Range1d(-15, 115)
return plot
So what's our problem look like right now?
plot = build_pedicabs_and_customers_plot()
# Connect pickups to dropoffs
for p, d in zip(pickup_locations, dropoff_locations):
plot.line(
[p[0], d[0]], [p[1], d[1]],
line_alpha=0.5,
line_color='grey',
line_dash='dashed'
)
show(plot)
Column generation is a popular technique for solving problems like this. We generate columns (i.e. binary variables) which are pre-defined routes.
- $\Omega = \{\text{set of allowable routes}\}$
- $\Omega_i = \{\text{set of routes covering customer i}\}$
- $\Omega_j = \{\text{set of routes involving pedicab j}\}$
- $c_r = \{\text{cost of route r}\}$
- $x_r = 1$ if we use route $r$, and $0$ otherwise.
Our model minimizes overall cost while covering each customer and assigning no more than one route to each pedicab.
$$ \begin{align} \text{min} \ \ \ & \ \sum_r c_r x_r \ & \ \\ \text{s.t.} \ \ \ & \ \sum_{r \in \Omega_i} x_r = 1 & \ \forall\ i \\ & \ \sum_{r \in \Omega_j} x_r \le 1 & \ \forall\ j \\ & x_r \in \{0, 1\} \end{align} $$We need to track routes, their costs, and what actors they are associated with.
from collections import defaultdict
import math
l2 = lambda pt1, pt2: math.sqrt(sum((x1-x2)**2 for x1, x2 in zip(pt1, pt2)))
routes = []
route_costs = {}
routes_by_pedicab = defaultdict(list)
routes_by_customer = defaultdict(list)
def add_route(pedicab, customers, locations):
stops = [pedicab_locations[pedicab]] + list(locations)
route_num = len(routes)
routes.append(stops)
route_costs[route_num] = sum(l2(s1, s2) for s1, s2 in zip(stops, stops[1:]))
routes_by_pedicab[pedicab].append(route_num)
for cust in customers:
routes_by_customer[cust].append(route_num)
Generate $\Omega$! Let's assume we're not interested in assigning more than 2 customers to any pedicab.
from itertools import permutations
pickups_dropoffs = list(zip(pickup_locations, dropoff_locations))
for cab in pedicabs:
# Single-customer routes
for cust in customers:
add_route(cab, [cust], pickups_dropoffs[cust])
# Doubles
for c1, c2 in permutations(customers, 2):
p1, d1 = pickups_dropoffs[c1]
p2, d2 = pickups_dropoffs[c2]
add_route(cab, [c1, c2], [p1, p2, d1, d2])
add_route(cab, [c1, c2], [p1, p2, d2, d1])
add_route(cab, [c1, c2], [p1, d1, p2, d2])
len(routes)
Binary variables indicate which routes to use.
x = []
for r in range(len(routes)):
x.append(pulp.LpVariable('x%d' % r, cat=pulp.LpBinary))
The entire constraint set to the problem:
prm = pulp.LpProblem('Pedicab Routing Model', sense=pulp.LpMinimize)
# No more than one route assigned to each cab
for cab_routes in routes_by_pedicab.values():
prm += sum(x[r] for r in cab_routes) <= 1
# Each customer is picked up by exactly one pedicab
for cust_routes in routes_by_customer.values():
prm += sum(x[r] for r in cust_routes) == 1
Objective = minimize total distance.
z = pulp.LpVariable('z', lowBound=0)
prm += z == sum(route_costs[r] * x[r] for r in range(len(routes)))
prm.setObjective(z)
assert prm.solve() == pulp.LpStatusOptimal
best_total_cost = z.value()
print('total cost:', best_total_cost)
We need a function to plot solutions.
def plot_routes():
plot = build_pedicabs_and_customers_plot()
# Show
for r, stops in enumerate(routes):
if x[r].value() < 0.5:
continue
plot.line(
[s[0] for s in stops],
[s[1] for s in stops],
line_alpha=0.5,
line_color='grey',
line_dash='dashed'
)
show(plot)
plot_routes()
Some of our pedicabs complain that they have to travel much farther than others. Let's fix that by minimizing the maximum trip length.
z2 = pulp.LpVariable('z2', lowBound=0)
for r in range(len(routes)):
prm += z2 >= route_costs[r] * x[r]
prm.setObjective(z2)
assert prm.solve() == pulp.LpStatusOptimal
print('total cost:', z.value())
print('max route length:', z2.value())
plot_routes()
That can hurt our global solution a lot. What if we are willing to accept a 10% cut to the global solution for more balance among individual routes?
# Only willing to accept 10% degradation in global optimality
prm += z <= best_total_cost * 1.1
assert prm.solve() == pulp.LpStatusOptimal
print('total cost:', z.value())
print('max route length:', z2.value())
plot_routes()
Where to find more information¶
- Introduction to Operations Research by Hillier and Lieberman
- Linear and Nonlinear Optimization by Griva, Nash, and Sofer
- Applied Integer Programming by Chen, Batson, and Dang
- Model Building in Mathematical Programming by Williams