GNU LESSER GENERAL PUBLIC LICENSE
Version 3, 29 June 2007
Copyright (C) 2007 Free Software Foundation, Inc.
Everyone is permitted to copy and distribute verbatim copies of this license document, but changing it is not allowed.
</font>The Set Cover Location Problem
Minimize
$\sum_{j=1}^n$ djxj
Subject to
$∑_{j=1}^n$aijxj ≥ 1, i = 1 ∊ n
xj ∊(0,1) j = 1 ∊ m
where
− i = a specific origin
− j = a specific destination
− n = the set of origins
− m = the set of destinations
− xi = the decision variable at each node in the matrix
− aij = binary matrix describing the coverage of each node
where
− aij = 1 ∀ i, j ∋ dij ≥ S (S is user defined)
− aij = 0 otherwise
− dij = distance from ith origin to jth destination
Adapted from:
Daskin, M. S. 1995. Network and Discrete Location: Models, Algorithms, and Applications. Hoboken, NJ, USA: John Wiley & Sons, Inc.
0. Imports
In [1]:
import pysal as ps
import numpy as np
import networkx as nx
import shapefile as shp
import gurobipy as gbp
import cplex as cp
import datetime as dt
import time
from collections import OrderedDict
import IPython.display as IPd
%pylab inline
from mpl_toolkits.basemap import Basemap
1. Data preparation and creation
1.1 Instantiate a network
In [2]:
ntw = ps.Network('Waverly/Waverly.shp')
print dt.datetime.now()
print 'PySAL.Network\n'
print dir(ntw)
1.2 Instantiate all graphs to be drawn
In [3]:
# Roads and Nodes
g = nx.Graph()
# Graph of Roads and Nodes
g1 = nx.MultiGraph()
# Clients
GRAPH_client = nx.Graph()
# Snapped Clients
g_client = nx.Graph()
# Service
GRAPH_service = nx.Graph()
# Snapped Service
g_service = nx.Graph()
# Gurobi p-Median
GUROBI_setcover_g = nx.Graph()
# Cplex p-Median
CPLEX_setcover_g = nx.Graph()
1.3 Create Bounding Box from 'Waverly.shp'
In [4]:
shp_W = ps.open('Waverly/Waverly.shp')
shp_W.bbox
Out[4]:
1.4 Create numpy arrays of random floats within a bounding box
In [5]:
lat_client = np.random.uniform(shp_W.bbox[0], shp_W.bbox[2], 500)
lon_client = np.random.uniform(shp_W.bbox[1], shp_W.bbox[3], 500)
lat_service = np.random.uniform(shp_W.bbox[0], shp_W.bbox[2], 500)
lon_service = np.random.uniform(shp_W.bbox[1], shp_W.bbox[3], 500)
1.5 Zip the latitude and longitude lists together
In [6]:
rand_coords_client = map(list, zip(lat_client, lon_client))
rand_coords_service = map(list, zip(lat_service, lon_service))
1.6 Create Empty Random Points Dictionaries
In [7]:
points_client = {}
points_service = {}
1.7 Fill dictionaries of random roints
In [8]:
# CLIENT
for idx, coords in enumerate(rand_coords_client):
GRAPH_client.add_node(idx)
points_client[idx] = coords
GRAPH_client.node[idx] = coords
# SERVICE
for idx, coords in enumerate(rand_coords_service):
GRAPH_service.add_node(idx)
points_service[idx] = coords
GRAPH_service.node[idx] = coords
1.8 Draw network, simplified network, and random client & service nodes
In [9]:
print dt.datetime.now()
#Instantiate Figure
figsize(10,10)
# Draw Graph of Actual Nodes and Roads
for e in ntw.edges:
g.add_edge(*e)
nx.draw(g, ntw.node_coords, node_size=5, alpha=0.25, edge_color='r', width=2)
# Draw only unique edges in graph
for e in ntw.graphedges:
g1.add_edge(*e)
# highlights cases where start and end node are the same
if e[0]==e[1]:
g1.add_node(e[0])
for node_id in g1.node:
g1.node[node_id] = ntw.node_coords[node_id]
nx.draw(g1, ntw.node_coords, node_size=10, alpha=0.5)
# Draw Graph of Random Client Points
nx.draw(GRAPH_client, points_client,
node_size=75, alpha=1, node_color='b')
# Draw Graph of Random Service Points
nx.draw(GRAPH_service, points_service,
node_size=100, alpha=1, node_color='c')
# Legend (Ordered Dictionary)
LEGEND = OrderedDict()
LEGEND['Network Nodes']=g
LEGEND['Roads']=g
LEGEND['Graph Vertices']=g1
LEGEND['Graph Edges']=g1
LEGEND['Client Nodes']=GRAPH_client
LEGEND['Service Nodes']=GRAPH_service
legend(LEGEND, loc='lower right', fancybox=True, framealpha=0.5)
# Title
title('Waverly Hills\n Tallahassee, Florida', family='Times New Roman',
size=40, color='k', backgroundcolor='w', weight='bold')
# North Arrow and 'N' --> Must be changed for different spatial resolutions, etc.
arrow(-84.281, 30.507, 0.0, 0.01, width=.0003, head_width=0.0012,
head_length=0.002, fc='k', ec='k',alpha=0.75,)
annotate('N', xy=(-84.2815, 30.52), fontstyle='italic', fontsize='xx-large',
fontweight='heavy', alpha=0.75)
Out[9]:
1.9 Create S
In [10]:
S = .5
1.10 Instantiate client and service shapefiles
In [11]:
# Client
client = shp.Writer(shp.POINT)
# Add Random Points
for i,j in rand_coords_client:
client.point(i,j)
# Add Fields
client.field('client_ID')
#client.field('Weight')
client.field('LAT')
client.field('LON')
counter = 0
for i in range(len(rand_coords_client)):
counter = counter + 1
client.record('client_' + str(counter), lat_client[i], lon_client[i])
# Save Shapefile
client.save('shapefiles/RandomPoints_CLIENT')
#Service
service = shp.Writer(shp.POINT)
# Add Random Points
for i,j in rand_coords_service:
service.point(i,j)
# Add Fields
service.field('y_ID')
service.field('x_ID')
service.field('LAT')
service.field('LON')
counter = 0
for i in range(len(rand_coords_service)):
counter = counter + 1
service.record('y' + str(counter), 'x' + str(counter), lat_service[i], lon_service[i])
# Save Shapefile
service.save('shapefiles/RandomPoints_SERVICE')
1.11 Snap Observations to NTW
In [12]:
t1 = time.time()
ntw.snapobservations('shapefiles/RandomPoints_CLIENT.shp',
'Rand_Points_CLIENT', attribute=True)
ntw.snapobservations('shapefiles/RandomPoints_SERVICE.shp',
'Rand_Points_SERVICE', attribute=True)
print round(time.time()-t1, 4), 'seconds'
1.12 Draw NTW, snapped coords, & random coords
In [13]:
# Instantiate Figure
figsize(10,10)
# Draw Graph of Roads
for e in ntw.edges:
g.add_edge(*e)
nx.draw(g, ntw.node_coords, node_size=5, alpha=0.25, edge_color='r', width=2)
# Draw Graph of Snapped Client Nodes
g_client = nx.Graph()
for p,coords in ntw.pointpatterns['Rand_Points_CLIENT'].snapped_coordinates.iteritems():
g_client.add_node(p)
g_client.node[p] = coords
nx.draw(g_client, ntw.pointpatterns['Rand_Points_CLIENT'].snapped_coordinates,
node_size=100, alpha=1, node_color='b')
# Draw Graph of Snapped Service Nodes
g_service = nx.Graph()
for p,coords in ntw.pointpatterns['Rand_Points_SERVICE'].snapped_coordinates.iteritems():
g_service.add_node(p)
g_service.node[p] = coords
nx.draw(g_service, ntw.pointpatterns['Rand_Points_SERVICE'].snapped_coordinates,
node_size=100, alpha=1, node_color='c')
# Draw Graph of Random Client Points
nx.draw(GRAPH_client, points_client,
node_size=20, alpha=1, node_color='y')
# Draw Graph of Random Client Points
nx.draw(GRAPH_service, points_service,
node_size=20, alpha=1, node_color='w')
# Legend (Ordered Dictionary)
LEGEND = OrderedDict()
LEGEND['Network Nodes']=g
LEGEND['Roads']=g
LEGEND['Snapped Client']=g_client
LEGEND['Snapped Service']=g_service
LEGEND['Client Nodes']=GRAPH_client
LEGEND['Service Nodes']=GRAPH_service
legend(LEGEND, loc='lower right', fancybox=True, framealpha=0.5)
# Title
title('Waverly Hills\n Tallahassee, Florida', family='Times New Roman',
size=40, color='k', backgroundcolor='w', weight='bold')
# North Arrow and 'N' --> Must be changed for different spatial resolutions, etc.
arrow(-84.281, 30.507, 0.0, 0.01, width=.0003, head_width=0.0012,
head_length=0.002, fc='k', ec='k',alpha=0.75,)
annotate('N', xy=(-84.2815, 30.52), fontstyle='italic', fontsize='xx-large',
fontweight='heavy', alpha=0.75)
Out[13]:
1.13 Create distance matrix
In [14]:
t1 = time.time()
All_Neigh_Dist = ntw.allneighbordistances(sourcepattern=ntw.pointpatterns['Rand_Points_CLIENT'],
destpattern=ntw.pointpatterns['Rand_Points_SERVICE'])
All_Dist_MILES = All_Neigh_Dist * float(10000/90) * 0.6214
seconds = round(time.time()-t1, 4)
print seconds, 'seconds'
print 'Client [i] x Service [j] Matrix Shape --> ', All_Dist_MILES.shape
2. Mathematical Optimization
2.1a Gurobi Set Cover test
In [15]:
t1 = time.time()
# 1. Read In Data
# Cost Vector
Cij = All_Dist_MILES
# Create Aij: Determine Aij (nodes within S)
# S --> 1 = served; 0 = unserved
Aij = []
for i in np.nditer(Cij):
if i <= S:
outtext = 1
else:
outtext = 0
Aij.append(outtext)
rows, cols = Cij.shape
Aij = np.array(Aij)
Aij = Aij.reshape(len(Cij),len(Cij[0]))
client_nodes = range(len(Cij))
service_nodes = range(len(Cij[0]))
# 2. Create Model, Set MIP Focus, Add Variables, & Update Model
mSCLP_GUROBI = gbp.Model(" -- SCLP -- ")
# Set MIP Focus to 2 for optimality
gbp.setParam('MIPFocus', 2)
# 3. Add Service Decision Variables
serv_var = []
for dest in service_nodes:
serv_var.append(mSCLP_GUROBI.addVar(vtype=gbp.GRB.BINARY,
ub = 1,
name='x'+str(dest+1)))
# Update Model Variables
mSCLP_GUROBI.update()
# 4. Set Objective Function
mSCLP_GUROBI.setObjective(gbp.quicksum(serv_var[dest]
for dest in service_nodes),
gbp.GRB.MINIMIZE)
# 5. Add Constraints
#Add Coverage Constraints
for orig in client_nodes:
mSCLP_GUROBI.addConstr(gbp.quicksum(Aij[orig][dest]*serv_var[dest]
for dest in service_nodes) >= 1)
# 6. Optimize and Print Results
mSCLP_GUROBI.optimize()
mSCLP_GUROBI.write('GUROBI_LP.lp')
t2G = time.time()-t1
print '\n*****************************************************************************************'
selected = []
dbf1 = ps.open('shapefiles/RandomPoints_SERVICE.dbf')
NEW_Records_SCLP_GUROBI = []
for v in mSCLP_GUROBI.getVars():
if v.x > 0:
var = '%s' % v.VarName
selected.append(v.x)
for i in range(dbf1.n_records):
if var in dbf1.read_record(i):
x = dbf1.read_record(i)
NEW_Records_SCLP_GUROBI.append(x)
else:
pass
print ' | ', var
print ' | Selected Facility Locations --------------------------- ^^^^ '
print ' | Coverage (S) in miles --------------------------------- ', S
print ' | Client Nodes ------------------------------------------ ', len(client_nodes)
print ' | Facilities needed for 100% coverage of client nodes --- ', len(selected)
print ' | Real Time to Optimize (sec.) -------------------------- ', t2G
print ' | Date/Time --------------------------------------------- ', dt.datetime.now()
print '*****************************************************************************************'
print '\nJames Gaboardi, 2015'
2.1b Instantiate Selected Gurobi Set Cover shapefile
In [16]:
SHP_SetCover_GUROBI = shp.Writer(shp.POINT)
# Add Points
for idy,idx,x,y in NEW_Records_SCLP_GUROBI:
SHP_SetCover_GUROBI.point(float(x), float(y))
# Add Fields
SHP_SetCover_GUROBI.field('y_ID')
SHP_SetCover_GUROBI.field('x_ID')
SHP_SetCover_GUROBI.field('LAT')
SHP_SetCover_GUROBI.field('LON')
# Add Records
for idy,idx,x,y in NEW_Records_SCLP_GUROBI:
SHP_SetCover_GUROBI.record(idy,idx,x,y)
# Save Shapefile
SHP_SetCover_GUROBI.save('shapefiles/Selected_Locations_SCLP_GUROBI')
2.2a Cplex Set Cover test
In [17]:
t1 = time.time()
# 1. Read In Data
# Cost Vector
Cij = All_Dist_MILES
# Create Aij: Determine Aij (nodes within S)
# S --> 1 = served; 0 = unserved
Cij = All_Dist_MILES
# Create Aij: Determine Aij (nodes within S)
# S --> 1 = served; 0 = unserved
Aij = []
for i in np.nditer(Cij):
if i <= S:
outtext = 1
else:
outtext = 0
Aij.append(outtext)
Aij = np.array(Aij)
Aij = Aij.reshape(len(Cij),len(Cij[0]))
# Indices & Variable Names
nodes = len(Cij)
Nodes = range(len(Cij))
all_nodes = len(Cij) * len(Cij)
ALL_nodes = range(all_nodes)
x = 'x'
cli_var = []
for i in Nodes:
for j in Nodes:
temp = x + str(j+1)
cli_var.append(temp)
client_var = np.array(cli_var)
results_var = []
for i in Nodes:
temp = x + str(i+1)
results_var.append(temp)
# 2. Create Model and Add Variables
# Create Model
mSCLP_CPLEX = cp.Cplex()
# Problem Name
mSCLP_CPLEX.set_problem_name('\n -- Set Cover Location Problem -- ')
print mSCLP_CPLEX.get_problem_name()
# Problem Type ==> Linear Programming
mSCLP_CPLEX.set_problem_type(mSCLP_CPLEX.problem_type.LP)
# Set MIP Emphasis to '2' --> Optimal
mSCLP_CPLEX.parameters.emphasis.mip.set(2)
print mSCLP_CPLEX.parameters.get_changed()
print '\nProblem Type\n ' + str(mSCLP_CPLEX.problem_type[mSCLP_CPLEX.get_problem_type()])
# Objective Function Sense ==> Minimize
mSCLP_CPLEX.objective.set_sense(mSCLP_CPLEX.objective.sense.minimize)
print 'Objective Sense\n ' + str(mSCLP_CPLEX.objective.sense[mSCLP_CPLEX.objective.get_sense()])
# Add Client Decision Variables
mSCLP_CPLEX.variables.add(names = [cli_var[i] for i in Nodes],
obj = [1] * nodes,
lb = [0] * nodes,
ub = [1] * nodes,
types = ['B'] * nodes)
# 3. Add Constraints
#Add Coverage Constraints
for orig in Nodes:
coverage_constraints = cp.SparsePair(ind = [client_var[dest]
for dest in Nodes],
val = [Aij[orig][dest]for dest in Nodes])
mSCLP_CPLEX.linear_constraints.add(lin_expr = [coverage_constraints],
senses = ['G'],
rhs = [1]);
# 4. Optimize and Print Results
mSCLP_CPLEX.solve()
t2C = time.time()-t1
solution = mSCLP_CPLEX.solution
mSCLP_CPLEX.write('CPLEX_LP.lp')
selected = []
dbf1 = ps.open('shapefiles/RandomPoints_SERVICE.dbf')
NEW_Records_SCLP_CPLEX = []
for v in mSCLP_CPLEX.variables.get_names():
if (solution.get_values(v) >
mSCLP_CPLEX.parameters.mip.tolerances.integrality.get()):
var = '%s' % v
selected.append(var)
for i in range(dbf1.n_records):
if var in dbf1.read_record(i):
x = dbf1.read_record(i)
NEW_Records_SCLP_CPLEX.append(x)
else:
pass
# solution.get_status() returns an integer code
print 'Solution status = ' , solution.get_status(), ':',
# the following line prints the corresponding string
print solution.status[solution.get_status()]
# Display solution.
print 'Selected Facility Locations = ' , solution.get_objective_value()
print 'Determination Time to Build Model in Python and Optimize (sec.): ', mSCLP_CPLEX.get_dettime(), 'ticks'
print 'Real Time to Build Model in Python and Optimize (sec.): ', t2C
print '****************************'
for f in results_var:
if (solution.get_values(f) >
mSCLP_CPLEX.parameters.mip.tolerances.integrality.get()):
print ' Facility %s is open' % f
else:
pass #print ' Facility %s is closed' % f
print '****************************'
print '\n----- Date/Time ------------------- ', dt.datetime.now()
print '\n----- Cplex Set Cover Loation Problem -----'
print '\n-----\nJames Gaboardi, 2015'
2.2b Instantiate Selected Cplex Set Cover shapefile
In [18]:
SHP_SCLP_CPLEX = shp.Writer(shp.POINT)
# Add Points
for idy,idx,x,y in NEW_Records_SCLP_CPLEX:
SHP_SCLP_CPLEX.point(float(x), float(y))
# Add Fields
SHP_SCLP_CPLEX.field('y_ID')
SHP_SCLP_CPLEX.field('x_ID')
SHP_SCLP_CPLEX.field('LAT')
SHP_SCLP_CPLEX.field('LON')
# Add Records
for idy,idx,x,y in NEW_Records_SCLP_CPLEX:
SHP_SCLP_CPLEX.record(idy,idx,x,y)
# Save Shapefile
SHP_SCLP_CPLEX.save('shapefiles/Selected_Locations_SCLP_CPLEX')
3. Selected locations
3.1 Gurobi & Cplex Set Cover Selected locations
In [19]:
figsize(10,10)
# Draw Network Actual Roads and Nodes
nx.draw(g, ntw.node_coords, node_size=5, alpha=0.25, edge_color='r', width=2)
# Draw Graph
nx.draw(g1, ntw.node_coords, node_size=10, alpha=0.5)
# Gurobi Set Cover
SetCover_GUROBI = ps.open('shapefiles/Selected_Locations_SCLP_GUROBI.shp')
points_SetCover_GUROBI = {}
for idx, coords in enumerate(SetCover_GUROBI):
GUROBI_setcover_g.add_node(idx)
points_SetCover_GUROBI[idx] = coords
GUROBI_setcover_g.node[idx] = coords
nx.draw(GUROBI_setcover_g, points_SetCover_GUROBI,
node_size=600, alpha=1, node_color='g')
# Cplex Set Cover
SetCover_CPLEX = ps.open('shapefiles/Selected_Locations_SCLP_CPLEX.shp')
points_SetCover_CPLEX = {}
for idx, coords in enumerate(SetCover_CPLEX):
CPLEX_setcover_g.add_node(idx)
points_SetCover_CPLEX[idx] = coords
CPLEX_setcover_g.node[idx] = coords
nx.draw(CPLEX_setcover_g, points_SetCover_CPLEX,
node_size=300, alpha=1, node_color='r')
# Draw Graph of Random Service
nx.draw(GRAPH_client, points_client,
node_size=15, alpha=.5, node_color='k')
# Draw Graph of Random Service
nx.draw(GRAPH_service, points_service,
node_size=50, alpha=1, node_color='k')
# Legend (Ordered Dictionary)
LEGEND = OrderedDict()
LEGEND['Network Nodes']=g
LEGEND['Roads']=g
LEGEND['Graph Vertices']=g1
LEGEND['Graph Edges']=g1
LEGEND['Gurobi Optimal Set Cover']=GUROBI_setcover_g
LEGEND['Cplex Optimal Set Cover']=CPLEX_setcover_g
LEGEND['Client Nodes']=GRAPH_client
LEGEND['Service Nodes']=GRAPH_service
legend(LEGEND, loc='lower right', fancybox=True, framealpha=0.5)
# Title
title('Waverly Hills\n Tallahassee, Florida', family='Times New Roman',
size=40, color='k', backgroundcolor='w', weight='bold')
# North Arrow and 'N' --> Must be changed for different spatial resolutions, etc.
arrow(-84.281, 30.507, 0.0, 0.01, width=.0003, head_width=0.0012,
head_length=0.002, fc='k', ec='k',alpha=0.75,)
annotate('N', xy=(-84.2815, 30.52), fontstyle='italic', fontsize='xx-large',
fontweight='heavy', alpha=0.75)
Out[19]:
3.2 Optimized Values
In [20]:
print '********************************************************'
print ' | Total Facilities to be opened for [ S =', S, ']'
print ' | | Gurobi ------------------ ', mSCLP_GUROBI.objVal, ' '
print ' | | CPLEX ------------------- ', solution.get_objective_value(), ' '
print '-------------------------------------------------------'
NSG = NEW_Records_SCLP_GUROBI
NSC = NEW_Records_SCLP_CPLEX
G_val = mSCLP_GUROBI.objVal
C_val = solution.get_objective_value()
if G_val == C_val:
print ' | Gurobi and CPLEX chose the same number of facilities'
else:
print ' | Gurobi and CPLEX chose different numbers of facilities'
if NSG == NSC:
print ' | ***'
print ' | Gurobi and CPLEX chose the same facilities'
else:
print ' | ***'
print ' | Gurobi and CPLEX chose different facilities'
print '-------------------------------------------------------'
print ' | Total Time to Build Model and Optimize: (seconds) '
print ' | | Gurobi ------------------ ', t2G, ' '
print ' | | CPLEX ------------------- ', t2C, ' '
print '********************************************************'
3.3 System Information
In [21]:
print '********************************************************'
print ' | Platform Specs: |'
print ' | | OS X Yosemite v10.10.4 |'
print ' | | MacBook Pro (Retina, 13-inch, Early 2015) |'
print ' | | Processor: 3.1 GHz Intel Core i7 |'
print ' | | Memory: 16 GB 1867 MHz DDR3 |'
print '********************************************************'
print ' | | Date/Time --------- ', dt.datetime.now(), '|'
print '********************************************************'