HW 1¶
By: Chengyi (Jeff) Chen
%load_ext autotime
%load_ext nb_black
%matplotlib inline
import os
from collections import defaultdict
import torch
import numpy as np
import pandas as pd
import cvxpy as cp
import scipy.stats
from scipy import optimize
# from sympy import *
# from torch.distributions import constraints
import matplotlib.pyplot as plt
from matplotlib import rc
# import daft
plt.rcParams["figure.dpi"] = 300
plt.rcParams["figure.figsize"] = (16, 12)
# import pyro
# import pyro.distributions as dist
# from pyro import poutine
# from pyro.infer.autoguide import AutoDelta
# from pyro.optim import Adam
# from pyro.infer import SVI, TraceEnum_ELBO, config_enumerate, infer_discrete
# smoke_test = "CI" in os.environ
# assert pyro.__version__.startswith("1.5.1")
# pyro.enable_validation(True)
time: 1.8 s (started: 2021-02-02 11:19:25 +08:00)
ohio_population = pd.read_excel(
"./data/Ohio_pop_5y.xlsx", engine="openpyxl", usecols=[0, 1, 2, 3, 4], nrows=88
)
ohio_population.rename({"Adjacent Matrix": "Adjacency List"}, inplace=True, axis=1)
ohio_population["Current Population"] = ohio_population["Current Population"].astype(
int
)
ohio_population["Population in 5 years later"] = ohio_population[
"Population in 5 years later"
].astype(int)
ohio_population["County Label"] = ohio_population["County Label"].astype(int)
ohio_population["Adjacency List"] = ohio_population["Adjacency List"].apply(
lambda x: np.array(list(map(int, x.split(","))))
)
ohio_population.head()
| County | Current Population | County Label | Adjacency List | Population in 5 years later | |
|---|---|---|---|---|---|
| 0 | Adams County | 27724 | 30 | [29, 31, 48, 49] | 27268 |
| 1 | Allen County | 102663 | 16 | [4, 15, 17, 22, 23] | 100495 |
| 2 | Ashland County | 53745 | 59 | [41, 56, 57, 58, 61, 62, 63] | 53836 |
| 3 | Ashtabula County | 97493 | 80 | [79, 78, 81] | 95249 |
| 4 | Athens County | 65818 | 68 | [52, 53, 66, 67, 69, 71] | 65904 |
time: 294 ms (started: 2021-02-02 11:19:27 +08:00)
Let’s convert the Adjacency List representing the counties that share a border into an Adjacency Matrix.
01 - Set Cover: Minimum Number of Principle places of business to cover all 88 counties in Ohio¶
Prior to 1979, banks in Ohio were only allowed to place branches in a county where the bank had a principal place of business. A new law in 1979 allowed banks to put branches in any county where the bank has a principal place of business and in any county adjacent to one in which it has a principal place of business. The question posed is, What is the minimum number of principal places of business, and in which counties should they be located to enable branches in all eighty-eight counties of Ohio?
The model for this problem is a classic set cover problem where the constraints ensure that for each of the 88 counties, the county has a principal place of business or a county adjacent to it has a principal place of business. Let \(x_j = 1\) if county \(j\) has a principal place of business and \(0\) if not; \(j = 1, 2, \cdots, 88\).
Define adjacency matrix \(A = \left[a_{ij} = 1\text{ if county }i\text{ and county }j\text{ share a border and }0\right]\) if not (note that \(a_{ii} = 1\)). We need to solve the following binary program:
A = np.array(
[
[
1 if (i + 1 in adjacent_counties) or (county_label - 1 == i) else 0
for i in range(len(ohio_population))
]
for county_label, adjacent_counties in ohio_population[
["County Label", "Adjacency List"]
].to_numpy()
]
)
assert np.alltrue(
[
set(list((np.argwhere(row) + 1).flatten()))
- set(list(ohio_population["Adjacency List"].iloc[idx]))
== set([ohio_population["County Label"].iloc[idx]])
for idx, row in enumerate(A)
]
), print("Adjacency Matrix is incorrect.")
A
array([[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
...,
[1, 1, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0]])
time: 35.8 ms (started: 2021-02-02 11:19:27 +08:00)
\begin{aligned} \underset{x}{\text{minimize }} &{\bf 1}^\top x \ \text{subject to } &Ax \geq {\bf 1} \ &x_j \in {0, 1 } ,, j = 1, 2, \cdots, 88. \ \end{aligned}
# Create one vector optimization variable.
x = cp.Variable((88,), integer=True)
# Create constraints.
constraints = [
A.T @ x >= 1,
0 <= x,
x <= 1,
]
# Form objective.
obj = cp.Minimize(cp.sum(x))
# Form and solve problem.
prob = cp.Problem(obj, constraints)
prob.solve()
print("Mixed Integer Programming Solution")
print("=" * 40)
print(f"Status: {prob.status}")
print(f"The optimal value is: {prob.value}")
print(f"The optimal solution is: x = {(np.argwhere(x.value) + 1).flatten()}")
Mixed Integer Programming Solution
========================================
Status: optimal
The optimal value is: 15.0
The optimal solution is: x = [ 3 5 7 8 9 28 35 40 45 49 51 72 75 76 81]
time: 18.8 ms (started: 2021-02-02 11:19:27 +08:00)
Let’s save this as the first solution.
first_solution = x.value
time: 276 µs (started: 2021-02-02 11:19:27 +08:00)
To enumerate other solutions, we impose another constraint such that the new optimization problem must find a new county to include in its set of principal place of business.
# Create one vector optimization variable.
x = cp.Variable((88,), integer=True)
# Create constraints.
constraints = [
A.T @ x >= 1,
0 <= x,
x <= 1,
first_solution @ x <= sum(first_solution) - 1,
cp.sum(x) == sum(first_solution),
]
# Form objective.
obj = cp.Minimize(cp.sum(x))
# Form and solve problem.
prob = cp.Problem(obj, constraints)
prob.solve()
print("Mixed Integer Programming Solution")
print("=" * 40)
print(f"Status: {prob.status}")
print(f"The optimal value is: {prob.value}")
print(f"The optimal solution is: x = {(np.argwhere(x.value) + 1).flatten()}")
Mixed Integer Programming Solution
========================================
Status: infeasible
The optimal value is: inf
The optimal solution is: x = []
time: 19.1 ms (started: 2021-02-02 11:19:27 +08:00)
02 - Maximal Cover: Maximize the population reached, given a limit \(k\) on the number of principal places of business opened¶
If the client can implement 5 principal places of business in each of the next three years, then we should provide a plan for which of the 15 locations should be opened in each of the next three years. To provide a plan for implementation, we need a metric of the quality or contribution of each principal place of business location. The population of each county is available from census data, and we may use that as a proxy for banking business. We may use an extension of the set cover model to solve a related problem known as the maximal coverage problem (MCP). We can use an MCP model to maximize the population reached, given a limit k on the number of principal places of business opened. The MCP is defined as follows:
Let \(x_j = 1\) if we place a principal place of business in county \(j\) and \(0\) if not; \(j = 1, 2, \cdots, 88\).
Let \(y_i = 1\) if county \(i\) is covered and \(0\) if not; \(i = 1, 2, \cdots, 88\)
Define \({pop}_{i}\) to be the population of county \(i\).
\begin{aligned} \underset{x, y}{\text{maximize }} &{pop}^\top y \ \text{subject to } &Ax \geq y \ &{\bf 1}^\top x \geq k \ &x_j \in {0, 1 } ,, j = 1, 2, \cdots, 88, \ &y_j \in {0, 1 } ,, j = 1, 2, \cdots, 88. \ \end{aligned}
# Create one vector optimization variable.
x = cp.Variable((88,), integer=True)
y = cp.Variable((88,), integer=True)
k = 5
# Create constraints.
constraints = [
A.T @ x >= y,
cp.sum(x) <= k,
0 <= x,
x <= 1,
0 <= y,
y <= 1,
]
# Form objective.
obj = cp.Maximize(ohio_population["Current Population"].to_numpy() @ y)
# Form and solve problem.
prob = cp.Problem(obj, constraints)
prob.solve()
print("Mixed Integer Programming Solution")
print("=" * 40)
print(f"Status: {prob.status}\n")
print(f"Population Covered: {int(prob.value)}\n")
print(
f"Principal Places of Business selected in Year 1 = {(np.argwhere(x.value) + 1).flatten() if sum(x.value) < 88 else 'All'}\n"
)
print(f"Counties covered in Year 1 = {int(sum(y.value))}")
Mixed Integer Programming Solution
========================================
Status: optimal
Population Covered: 7989928
Principal Places of Business selected in Year 1 = [28 42 55 69 71]
Counties covered in Year 1 = 37
time: 25.3 ms (started: 2021-02-02 11:19:27 +08:00)
principals_year_1 = x.value
counties_covered_year_1 = y.value
time: 394 µs (started: 2021-02-02 11:19:27 +08:00)
# Create one vector optimization variable.
x = cp.Variable((88,), integer=True)
y = cp.Variable((88,), integer=True)
k = 10
# Create constraints.
constraints = [
A.T @ x >= y,
cp.sum(x) <= k,
0 <= x,
x <= 1,
0 <= y,
y <= 1,
principals_year_1 @ x >= sum(principals_year_1),
counties_covered_year_1 @ y >= sum(counties_covered_year_1),
]
# Form objective.
obj = cp.Maximize(ohio_population["Current Population"].to_numpy() @ y)
# Form and solve problem.
prob = cp.Problem(obj, constraints)
prob.solve()
print("Mixed Integer Programming Solution")
print("=" * 40)
print(f"Status: {prob.status}\n")
print(f"Population Covered: {int(prob.value)}\n")
print(
f"Principal Places of Business selected in Year 2 = {(np.argwhere(x.value) + 1).flatten() if sum(x.value) < 88 else 'All'}\n"
)
print(f"Number of Counties covered in Year 2 = {int(sum(y.value))}")
Mixed Integer Programming Solution
========================================
Status: optimal
Population Covered: 10663087
Principal Places of Business selected in Year 2 = [10 13 28 42 44 55 61 69 71 80]
Number of Counties covered in Year 2 = 64
time: 23.4 ms (started: 2021-02-02 11:19:27 +08:00)
principals_year_2 = x.value
counties_covered_year_2 = y.value
time: 373 µs (started: 2021-02-02 11:19:27 +08:00)
# Create one vector optimization variable.
x = cp.Variable((88,), integer=True)
y = cp.Variable((88,), integer=True)
k = 15
# Create constraints.
constraints = [
A.T @ x >= y,
cp.sum(x) <= k,
0 <= x,
x <= 1,
0 <= y,
y <= 1,
principals_year_1 @ x >= sum(principals_year_1),
counties_covered_year_1 @ y >= sum(counties_covered_year_1),
]
# Form objective.
obj = cp.Maximize(ohio_population["Current Population"].to_numpy() @ y)
# Form and solve problem.
prob = cp.Problem(obj, constraints)
prob.solve()
print("Mixed Integer Programming Solution")
print("=" * 40)
print(f"Status: {prob.status}\n")
print(f"Population Covered: {int(prob.value)}\n")
print(
f"Principal Places of Business selected in Year 3 = {(np.argwhere(x.value) + 1).flatten() if sum(x.value) < 88 else 'All'}\n"
)
print(f"Number of Counties covered in Year 3 = {int(sum(y.value))}")
Mixed Integer Programming Solution
========================================
Status: optimal
Population Covered: 11563347
Principal Places of Business selected in Year 3 = [ 1 7 19 27 28 39 42 48 55 58 69 71 76 80 83]
Number of Counties covered in Year 3 = 84
time: 20.2 ms (started: 2021-02-02 11:19:27 +08:00)
a. Build a model in Pyomo based on equation (25)-(29) in ‘How to Influence and Improve Decisions Through Optimization Models’ paper. Set \(k\) =14, and use the current population (‘Current Population’ column in the ‘Ohio_pop_5y.csv’). Solve the model, report your solution and objective value.¶
# Create one vector optimization variable.
x = cp.Variable((88,), integer=True)
y = cp.Variable((88,), integer=True)
k = 14
# Create constraints.
constraints = [
A.T @ x >= y,
cp.sum(x) <= k,
0 <= x,
x <= 1,
0 <= y,
y <= 1,
]
# Form objective.
obj = cp.Maximize(ohio_population["Current Population"].to_numpy() @ y)
# Form and solve problem.
prob = cp.Problem(obj, constraints)
prob.solve()
print("Mixed Integer Programming Solution")
print("=" * 40)
print(f"Status: {prob.status}\n")
print(f"Population Covered: {int(prob.value)}\n")
print(
f"Principal Places of Business selected with k = 14: {(np.argwhere(x.value) + 1).flatten() if sum(x.value) < 88 else 'All'}\n"
)
print(f"Counties covered: {int(sum(y.value))}")
Mixed Integer Programming Solution
========================================
Status: optimal
Population Covered: 11618698
Principal Places of Business selected with k = 14: [ 5 7 19 28 31 36 39 44 45 48 49 51 69 76]
Counties covered: 85
time: 44.8 ms (started: 2021-02-02 11:19:27 +08:00)
hw_a_population_covered = int(prob.value)
hw_a_principal_places_selected = (
(np.argwhere(x.value) + 1).flatten() if sum(x.value) < 88 else "All"
)
hw_a_counties_covered = int(sum(y.value))
hw_a_y_val = y.value
time: 1.32 ms (started: 2021-02-02 11:19:27 +08:00)
b. Build the same model in (a), but use the population in 5 years later(‘Population in 5 years later’ column in the dataset). Solve the model, report your solution and objective value.¶
# Create one vector optimization variable.
x = cp.Variable((88,), integer=True)
y = cp.Variable((88,), integer=True)
k = 14
# Create constraints.
constraints = [
A.T @ x >= y,
cp.sum(x) <= k,
0 <= x,
x <= 1,
0 <= y,
y <= 1,
]
# Form objective.
obj = cp.Maximize(ohio_population["Population in 5 years later"].to_numpy() @ y)
# Form and solve problem.
prob = cp.Problem(obj, constraints)
prob.solve()
print("Mixed Integer Programming Solution")
print("=" * 40)
print(f"Status: {prob.status}\n")
print(f"Population Covered 5 years later: {int(prob.value)}\n")
print(
f"Principal Places of Business selected with k = 14, 5 years later: {(np.argwhere(x.value) + 1).flatten() if sum(x.value) < 88 else 'All'}\n"
)
print(f"Counties covered 5 years later: {int(sum(y.value))}")
Mixed Integer Programming Solution
========================================
Status: optimal
Population Covered 5 years later: 11717667
Principal Places of Business selected with k = 14, 5 years later: [ 5 7 8 9 19 28 39 40 45 48 49 51 69 76]
Counties covered 5 years later: 85
time: 50.8 ms (started: 2021-02-02 11:19:27 +08:00)
hw_b_population_covered = int(prob.value)
hw_b_principal_places_selected = (
(np.argwhere(x.value) + 1).flatten() if sum(x.value) < 88 else "All"
)
hw_b_counties_covered = int(sum(y.value))
hw_b_y_val = y.value
time: 641 µs (started: 2021-02-02 11:19:28 +08:00)
c. Compare the solution in (a) and (b). Discuss your findings.¶
hw_a_population_covered
11618698
time: 2.72 ms (started: 2021-02-02 11:19:28 +08:00)
hw_b_population_covered
11717667
time: 1.64 ms (started: 2021-02-02 11:19:28 +08:00)
hw_a_principal_places_selected
array([ 5, 7, 19, 28, 31, 36, 39, 44, 45, 48, 49, 51, 69, 76])
time: 1.57 ms (started: 2021-02-02 11:19:28 +08:00)
hw_b_principal_places_selected
array([ 5, 7, 8, 9, 19, 28, 39, 40, 45, 48, 49, 51, 69, 76])
time: 2.3 ms (started: 2021-02-02 11:19:28 +08:00)
hw_a_counties_covered
85
time: 1.69 ms (started: 2021-02-02 11:19:28 +08:00)
hw_b_counties_covered
85
time: 2.82 ms (started: 2021-02-02 11:19:28 +08:00)
print(
f"Did (a) and (b) cover the same 85 counties? {np.alltrue(hw_a_y_val == hw_b_y_val)}"
)
Did (a) and (b) cover the same 85 counties? True
time: 587 µs (started: 2021-02-02 11:19:28 +08:00)
Seeing that both (a) and (b) cover the same counties, the population that they cover will be the same. Hence, fortunately for us, using the current population data was sufficient to choose a set of principal places of business that will cover the same counties and subsequently, the same population as if we had used the population data 5 years later. If the population data 5 years later was unavailable to us, we would have needed to model the population 5 years later and create confidence intervals on the population size in each county, and run the optimization with different estimates of the predicted population sizes.