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-Maxian Problem
Maximize
$\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:
Church, R. L., and R. S. Garfinkel. Locating an Obnoxious Facility on a Network. Transportation Science 12 (2):107–118.
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-Maxian
GUROBI_maxian_g = nx.Graph()
# Cplex p-Maxian
CPLEX_maxian_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], 100)
lon_client = np.random.uniform(shp_W.bbox[1], shp_W.bbox[3], 100)
lat_service = np.random.uniform(shp_W.bbox[0], shp_W.bbox[2], 50)
lon_service = np.random.uniform(shp_W.bbox[1], shp_W.bbox[3], 50)
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-Maxian test [p = 2]
In [15]:
t1 = time.time()
# 1. Data
# Weighted Costs
# Cost Matrix
Cij = All_Dist_MILES
# Weights Matrix
Ai = Ai
# Demand Sum
AiSum = AiSum
# Weighted Cost Coefficients for Decision Variables
Sij = Ai * Cij
client_nodes = range(len(Sij))
service_nodes = range(len(Sij[0]))
mPMaxP_GUROBI = gbp.Model(' -- p-Maxian -- ')
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(mPMaxP_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(mPMaxP_GUROBI.addVar(vtype=gbp.GRB.BINARY,
name='y'+str(dest+1)))
mPMaxP_GUROBI.update()
mPMaxP_GUROBI.setObjective(gbp.quicksum(Sij[orig][dest]*client_var[orig][dest]
for orig in client_nodes for dest in service_nodes),
gbp.GRB.MAXIMIZE)
for orig in client_nodes:
mPMaxP_GUROBI.addConstr(gbp.quicksum(client_var[orig][dest]
for dest in service_nodes) == 1)
for orig in service_nodes:
for dest in client_nodes:
mPMaxP_GUROBI.addConstr((serv_var[orig] - client_var[dest][orig] >= 0))
mPMaxP_GUROBI.addConstr(gbp.quicksum(serv_var[dest][0] for dest in service_nodes) == 2)
mPMaxP_GUROBI.optimize()
mPMaxP_GUROBI.write('WaverlyPMaxP.lp')
t2G = time.time()-t1
print '\n*************************************************************************'
selected = []
dbf1 = ps.open('shapefiles/RandomPoints_SERVICE.dbf')
NEW_Records_PMaxP_GUROBI = []
for v in mPMaxP_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_PMaxP_GUROBI.append(x)
else:
pass
print ' | ', var
print ' | Selected Facility Locations -------------- ^^^^ '
print ' | Candidate Facilities [p] ----------------- ', len(selected)
val_G = mPMaxP_GUROBI.objVal
print ' | Objective Value (miles) ------------------ ', val_G
avg_G = float(mPMaxP_GUROBI.objVal)/float(AiSum)
print ' | Avg. Value / Client (miles) -------------- ', avg_G
print ' | Real Time to Optimize (sec.) ------------- ', t2G
print ' | Date/Time -------------------------------- ', dt.datetime.now()
print '*************************************************************************'
print '\nJames Gaboardi, 2015'
2.1b Instantiate Selected Gurobi p-Maxian shapefile
In [16]:
SHP_Maxian_GUROBI = shp.Writer(shp.POINT)
# Add Points
for idy,idx,x,y in NEW_Records_PMaxP_GUROBI:
SHP_Maxian_GUROBI.point(float(x), float(y))
# Add Fields
SHP_Maxian_GUROBI.field('y_ID')
SHP_Maxian_GUROBI.field('x_ID')
SHP_Maxian_GUROBI.field('LAT')
SHP_Maxian_GUROBI.field('LON')
# Add Records
for idy,idx,x,y in NEW_Records_PMaxP_GUROBI:
SHP_Maxian_GUROBI.record(idy,idx,x,y)
# Save Shapefile
SHP_Maxian_GUROBI.save('shapefiles/Selected_Locations_Pmaxian_GUROBI')
2.2a Cplex p-Maxian 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
mPMaxP_CPLEX = cp.Cplex()
# Problem Name
mPMaxP_CPLEX.set_problem_name('\n -- P-Maxian -- ')
print mPMaxP_CPLEX.get_problem_name()
# Problem Type ==> Linear Programming
mPMaxP_CPLEX.set_problem_type(mPMaxP_CPLEX.problem_type.LP)
# Set MIP Emphasis to '2' --> Optimal
mPMaxP_CPLEX.parameters.emphasis.mip.set(2)
print mPMaxP_CPLEX.parameters.get_changed()
print '\nProblem Type\n ' + str(mPMaxP_CPLEX.problem_type[mPMaxP_CPLEX.get_problem_type()])
# Objective Function Sense ==> Maximize
mPMaxP_CPLEX.objective.set_sense(mPMaxP_CPLEX.objective.sense.maximize)
print 'Objective Sense\n ' + str(mPMaxP_CPLEX.objective.sense[mPMaxP_CPLEX.objective.get_sense()])
# Add Client Decision Variables
mPMaxP_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
mPMaxP_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]))
mPMaxP_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]))
mPMaxP_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])
mPMaxP_CPLEX.linear_constraints.add(lin_expr = [opening_constraint],
senses = ['G'],
rhs = [0])
# 4. Optimize and Print Results
mPMaxP_CPLEX.solve()
t2C = time.time()-t1
mPMaxP_CPLEX.write('WaverlyPMP_CPLEX.lp')
solution = mPMaxP_CPLEX.solution
selected = []
dbf1 = ps.open('shapefiles/RandomPoints_SERVICE.dbf')
NEW_Records_PMaxP_CPLEX = []
for v in mPMaxP_CPLEX.variables.get_names():
if 'x' in v:
pass
elif (solution.get_values(v) >
mPMaxP_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_PMaxP_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.): ', mPMaxP_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) >
mPMaxP_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-Maxian Problem -----'
print '\n-----\nJames Gaboardi, 2015'
2.2b Instantiate Selected Cplex p-Maxian shapefile
In [18]:
SHP_Maxian_CPLEX = shp.Writer(shp.POINT)
# Add Points
for idy,idx,x,y in NEW_Records_PMaxP_CPLEX:
SHP_Maxian_CPLEX.point(float(x), float(y))
# Add Fields
SHP_Maxian_CPLEX.field('y_ID')
SHP_Maxian_CPLEX.field('x_ID')
SHP_Maxian_CPLEX.field('LAT')
SHP_Maxian_CPLEX.field('LON')
# Add Records
for idy,idx,x,y in NEW_Records_PMaxP_CPLEX:
SHP_Maxian_CPLEX.record(idy,idx,x,y)
# Save Shapefile
SHP_Maxian_CPLEX.save('shapefiles/Selected_Locations_Pmaxian_CPLEX')
3. Selected locations
3.1 Gurobi & Cplex p-Maxian 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-Maxian
P_Max_GUROBI = ps.open('shapefiles/Selected_Locations_Pmaxian_GUROBI.shp')
points_max_GUROBI = {}
for idx, coords in enumerate(P_Max_GUROBI):
GUROBI_maxian_g.add_node(idx)
points_max_GUROBI[idx] = coords
GUROBI_maxian_g.node[idx] = coords
nx.draw(GUROBI_maxian_g, points_max_GUROBI,
node_size=600, alpha=1, node_color='g')
# Cplex p-Maxian
P_Max_CPLEX = ps.open('shapefiles/Selected_Locations_Pmaxian_CPLEX.shp')
points_max_CPLEX = {}
for idx, coords in enumerate(P_Max_CPLEX):
CPLEX_maxian_g.add_node(idx)
points_max_CPLEX[idx] = coords
CPLEX_maxian_g.node[idx] = coords
nx.draw(CPLEX_maxian_g, points_max_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-Maxian (p=2)']=GUROBI_maxian_g
LEGEND['Cplex Optimal p-Maxian (p=2)']=CPLEX_maxian_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_G, ' '
print ' | | CPLEX ------------------- ', solution.get_objective_value(), ' '
print '-------------------------------------------------------'
NpmG = NEW_Records_PMaxP_GUROBI
NpmC = NEW_Records_PMaxP_CPLEX
G_val = mPMaxP_GUROBI.objVal
C_val = solution.get_objective_value()
if G_val == C_val:
print ' | Gurobi and CPLEX came to equal Objective Values'
else:
print ' | Gurobi and CPLEX came to different Objective Values'
if NpmG == NpmC:
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 '********************************************************'