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 p-Median Problem
Minimize
$\sum_{i∊1}^n\sum_{j∊1}^m$ aicijxij
Subject to
$∑_{j∊m}$xij = 1, ∀ i ∊ n
$∑_{i∊n}$yj = p
xij - yj ≥ 0, ∀ i ∊ n, j ∊ m
xij, yj ∊(0,1) ∀ i ∊ n, j ∊ m
where
− i = a specific origin
− j = a specific destination
− n = the set of origins
− m = the set of destinations
− ai = weight at each node
− cij = travel costs between nodes
− xij = the decision variable at each node in the matrix
− yj = nodes chosen as service facilities
− p = the number of facilities to be sited
Adapted from:
Horner, M. W. and M. J. Widener. 2010. How do socioeconomic characteristics interact with equity and efficiency considerations? An analysis of hurricane disaster relief goods provision. Geospatial Analysis and Modelling of Urban Structure and Dynamics 99:393–414.
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_median_g = nx.Graph()
# Cplex p-Median
CPLEX_median_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], 600)
lon_client = np.random.uniform(shp_W.bbox[1], shp_W.bbox[3], 600)
lat_service = np.random.uniform(shp_W.bbox[0], shp_W.bbox[2], 300)
lon_service = np.random.uniform(shp_W.bbox[1], shp_W.bbox[3], 300)
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 weights at nodes and sum
In [10]:
# Client Weights
Ai = np.random.randint(5, 50, len(rand_coords_client))
Ai = Ai.reshape(len(Ai),1)
# Sum
AiSum = np.sum(Ai)
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), Ai[i], 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 p-Median test [p = 2]
In [15]:
t1 = time.time()
# 1. Data
# Demand
Ai = Ai
# Demand Sum
AiSum = AiSum
# Travel Costs
Cij = All_Dist_MILES
# Weighted Costs
Sij = Ai * Cij
client_nodes = range(len(Sij))
service_nodes = range(len(Sij[0]))
mPMP_GUROBI = gbp.Model(' -- p-Median -- ')
gbp.setParam('MIPFocus', 2)
# Client IxJ
client_var = []
for orig in client_nodes:
client_var.append([])
for dest in service_nodes:
client_var[orig].append(mPMP_GUROBI.addVar(vtype=gbp.GRB.BINARY,
obj=Sij[orig][dest],
name='x'+str(orig+1)+'_'+str(dest+1)))
#J
serv_var = []
for dest in service_nodes:
serv_var.append([])
serv_var[dest].append(mPMP_GUROBI.addVar(vtype=gbp.GRB.BINARY,
name='y'+str(dest+1)))
mPMP_GUROBI.update()
mPMP_GUROBI.setObjective(gbp.quicksum(Sij[orig][dest]*client_var[orig][dest]
for orig in client_nodes for dest in service_nodes),
gbp.GRB.MINIMIZE)
for orig in client_nodes:
mPMP_GUROBI.addConstr(gbp.quicksum(client_var[orig][dest]
for dest in service_nodes) == 1)
for orig in service_nodes:
for dest in client_nodes:
mPMP_GUROBI.addConstr((serv_var[orig] - client_var[dest][orig] >= 0))
mPMP_GUROBI.addConstr(gbp.quicksum(serv_var[dest][0] for dest in service_nodes) == 2)
mPMP_GUROBI.optimize()
t2P = time.time()-t1
mPMP_GUROBI.write('LP_Files/WaverlyPMP_GUROBI.lp')
print '\n*************************************************************************'
selected = []
dbf1 = ps.open('shapefiles/RandomPoints_SERVICE.dbf')
NEW_Records_PMP_GUROBI = []
for v in mPMP_GUROBI.getVars():
if 'x' in v.VarName:
pass
elif v.x > 0:
var = '%s' % v.VarName
selected.append(var)
for i in range(dbf1.n_records):
if var in dbf1.read_record(i):
x = dbf1.read_record(i)
NEW_Records_PMP_GUROBI.append(x)
else:
pass
print ' | ', var
print ' | Selected Facility Locations -------------- ^^^^ '
print ' | Candidate Facilities [p] ----------------- ', len(selected)
val = mPMP_GUROBI.objVal
print ' | Objective Value (miles) ------------------ ', val
avg = float(mPMP_GUROBI.objVal)/float(AiSum)
print ' | Avg. Value / Client (miles) -------------- ', avg
print ' | Real Time to Optimize (sec.) ------------- ', t2P
print ' | Date/Time -------------------------------- ', dt.datetime.now()
print '*************************************************************************'
print '-- Gurobi p-Median Problem --'
print '\nJames Gaboardi, 2015'
2.1b Instantiate Selected Gurobi p-median shapefile
In [16]:
SHP_Median_GUROBI = shp.Writer(shp.POINT)
# Add Points
for idy,idx,x,y in NEW_Records_PMP_GUROBI:
SHP_Median_GUROBI.point(float(x), float(y))
# Add Fields
SHP_Median_GUROBI.field('y_ID')
SHP_Median_GUROBI.field('x_ID')
SHP_Median_GUROBI.field('LAT')
SHP_Median_GUROBI.field('LON')
# Add Records
for idy,idx,x,y in NEW_Records_PMP_GUROBI:
SHP_Median_GUROBI.record(idy,idx,x,y)
# Save Shapefile
SHP_Median_GUROBI.save('shapefiles/Selected_Locations_Pmedian_GUROBI')
2.2a Cplex p-Median test [p = 2]
In [17]:
t1 = time.time()
# Demand
Ai = Ai
# Demand Sum
AiSum = AiSum
# Travel Costs
Cij = All_Dist_MILES
# Weighted Costs
Sij = Ai * Cij
client_nodes = range(len(Sij))
service_nodes = range(len(Sij[0]))
all_nodes = len(Sij) * len(Sij[0])
ALL_nodes = range(all_nodes)
x = 'x'
cli_var = []
for i in client_nodes:
for j in service_nodes:
temp = x + str(i+1) + '_' + str(j+1)
cli_var.append(temp)
client_var = np.array(cli_var)
client_var = client_var.reshape(len(Sij),len(Sij[0]))
y = 'y'
fac_var = []
for i in service_nodes:
temp = y + str(i+1)
fac_var.append(temp)
facility_var = np.array(fac_var)
facility_var = facility_var.reshape(1,len(Sij[0]))
# 2. Create Model and Add Variables
# Create Model
mPMP_CPLEX = cp.Cplex()
# Problem Name
mPMP_CPLEX.set_problem_name('\n -- P-Median -- ')
print mPMP_CPLEX.get_problem_name()
# Problem Type ==> Linear Programming
mPMP_CPLEX.set_problem_type(mPMP_CPLEX.problem_type.LP)
# Set MIP Emphasis to '2' --> Optimal
mPMP_CPLEX.parameters.emphasis.mip.set(2)
print mPMP_CPLEX.parameters.get_changed()
print '\nProblem Type\n ' + str(mPMP_CPLEX.problem_type[mPMP_CPLEX.get_problem_type()])
# Objective Function Sense ==> Minimize
mPMP_CPLEX.objective.set_sense(mPMP_CPLEX.objective.sense.minimize)
print 'Objective Sense\n ' + str(mPMP_CPLEX.objective.sense[mPMP_CPLEX.objective.get_sense()])
# Add Client Decision Variables
mPMP_CPLEX.variables.add(names = [cli_var[i] for i in ALL_nodes],
obj = [Sij[i][j] for i in client_nodes for j in service_nodes],
lb = [0] * all_nodes,
ub = [1] * all_nodes,
types = ['B'] * all_nodes)
# Add Service Decision Variable
mPMP_CPLEX.variables.add(names = [fac_var[j] for j in service_nodes],
lb = [0] * len(Sij[0]),
ub = [1] * len(Sij[0]),
types = ['B'] * len(Sij[0]))
# 3. Add Constraints
# Add Assignment Constraints
for orig in client_nodes:
assignment_constraints = cp.SparsePair(ind = [client_var[orig][dest]
for dest in service_nodes],
val = [1] * len(Sij[0]))
mPMP_CPLEX.linear_constraints.add(lin_expr = [assignment_constraints],
senses = ['E'],
rhs = [1]);
# Add Facility Constraint
facility_constraint = cp.SparsePair(ind = fac_var,
val = [1.0] * len(Sij[0]))
mPMP_CPLEX.linear_constraints.add(lin_expr = [facility_constraint],
senses = ['E'],
rhs = [2])
# Add Opening Constraint
cli_var_open = []
for i in client_nodes:
for j in service_nodes:
temp = x + str(i+1) + '_' + str(j+1)
cli_var_open.append(temp)
fac_var_open = []
for i in client_nodes:
for j in service_nodes:
temp = y + str(j+1)
fac_var_open.append(temp)
l = []
for i in ALL_nodes:
l.append([cli_var_open[i]]+[fac_var_open[i]])
for i in l:
opening_constraint = cp.SparsePair(ind = i, val = [-1.0, 1.0])
mPMP_CPLEX.linear_constraints.add(lin_expr = [opening_constraint],
senses = ['G'],
rhs = [0])
# 4. Optimize and Print Results
mPMP_CPLEX.solve()
t2C = time.time()-t1
mPMP_CPLEX.write('LP_Files/WaverlyPMP_CPLEX.lp')
solution = mPMP_CPLEX.solution
selected = []
dbf1 = ps.open('shapefiles/RandomPoints_SERVICE.dbf')
NEW_Records_PMP_CPLEX = []
for v in mPMP_CPLEX.variables.get_names():
if 'x' in v:
pass
elif (solution.get_values(v) >
mPMP_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_PMP_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 'Total cost: ' , solution.get_objective_value()
print 'Determination Time to Build Model in Python and Optimize (sec.): ', mPMP_CPLEX.get_dettime(), 'ticks'
print 'Real Time to Build Model in Python and Optimize (sec.): ', t2C
print '****************************'
for f in fac_var:
if (solution.get_values(f) >
mPMP_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 p-Median Problem -----'
print '\n-----\nJames Gaboardi, 2015'
2.2b Instantiate Selected Cplex p-median shapefile
In [18]:
SHP_Median_CPLEX = shp.Writer(shp.POINT)
# Add Points
for idy,idx,x,y in NEW_Records_PMP_CPLEX:
SHP_Median_CPLEX.point(float(x), float(y))
# Add Fields
SHP_Median_CPLEX.field('y_ID')
SHP_Median_CPLEX.field('x_ID')
SHP_Median_CPLEX.field('LAT')
SHP_Median_CPLEX.field('LON')
# Add Records
for idy,idx,x,y in NEW_Records_PMP_CPLEX:
SHP_Median_CPLEX.record(idy,idx,x,y)
# Save Shapefile
SHP_Median_CPLEX.save('shapefiles/Selected_Locations_Pmedian_CPLEX')
3. Selected locations
3.1 Gurobi & Cplex p-Median 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 p-Median
P_Med_GUROBI = ps.open('shapefiles/Selected_Locations_Pmedian_GUROBI.shp')
points_median_GUROBI = {}
for idx, coords in enumerate(P_Med_GUROBI):
GUROBI_median_g.add_node(idx)
points_median_GUROBI[idx] = coords
GUROBI_median_g.node[idx] = coords
nx.draw(GUROBI_median_g, points_median_GUROBI,
node_size=600, alpha=1, node_color='g')
# Cplex p-Median
P_Med_CPLEX = ps.open('shapefiles/Selected_Locations_Pmedian_CPLEX.shp')
points_median_CPLEX = {}
for idx, coords in enumerate(P_Med_CPLEX):
CPLEX_median_g.add_node(idx)
points_median_CPLEX[idx] = coords
CPLEX_median_g.node[idx] = coords
nx.draw(CPLEX_median_g, points_median_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 p-Median (p=2)']=GUROBI_median_g
LEGEND['Cplex Optimal p-Median (p=2)']=CPLEX_median_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 Cost: Objective Function Value (miles) '
print ' | | Gurobi ------------------ ', val, ' '
print ' | | CPLEX ------------------- ', solution.get_objective_value(), ' '
print '-------------------------------------------------------'
print ' | Total Time to Build Model and Optimize: (seconds) '
print ' | | Gurobi ------------------ ', t2P, ' '
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 '********************************************************'