notebook.community

Edit and run 





In [ ]:



    


#    GNU LESSER GENERAL PUBLIC LICENSE

#    Version 3, 29 June 2007

#    Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
#    Everyone is permitted to copy and distribute verbatim copies
#    of this license document, but changing it is not allowed.

Network-Based Model Building for Discrete Location Allocation in Python: Integrating PySAL and Gurobi

James D. Gaboardi, David C. Folch

Florida State University, Department of Geography

NARSC 2015

1. Optimization in GIS & the Current State

       • Imports

       • Reproducibility

       • Open Source

       • Proprietary

       • Selected Software

2. Demonstration of PySAL + Gurobi

       • Imports & Specs

       • Data Preparation & Creation

       • Mathematical Optimization

             → Set Cover Location Problem

             → p-Median Problem

             → p-Center Problem

             → p-Center Problem (Manual LP Creation)

       • Visualization of Selected Locations

3. Next Steps

Imports





In [3]:



    


import datetime as dt
import gurobipy as gbp
import IPython.display as IPd
import os
import platform
from prettytable import from_csv
import pysal as ps
import sys
import time

# Local path on user's machine
path = '/Users/jgaboardi/NARSC__15/'

1.     Optimization & GIS: Current State

Is it Reproducible?

Later in this section there is an example of the workflow for another paper I am currently working on using an accepted, traditional method. This involves manual data input and processing at every step in the schema which tends to decrease efficiency and reproducibility. This is very similar to an Operations Research approach to model building and Operations Research Methodology Diagram (Bertrand and Fransoo 2002, pulp documentation 2009).

Open Source Software Options

Many of the software options below don't implement the algorithms we need and this list is not exhaustive.

How do we choose and what do we choose?





In [5]:



    


# Table of Open-Source software for optimization
OpenOpt = open(path+'NARSC__Optimizers__Open.csv', "r")
OpenTable = from_csv(OpenOpt)
print OpenTable
OpenOpt.close()



    


















    






+------------------+------------------------+--------------+--------------------------------+
|   Open Source    |          Type          | last updated |             Notes              |
+------------------+------------------------+--------------+--------------------------------+
|     COIN-OR      | Open Source Initiative |              |   many resources associated    |
| CoinGraphClasses |        Library         |              |   network algorithms (COIN)    |
|       CMPL       |        Language        |              |             (COIN)             |
|      CoinMP      |   Interface Library    |              |      supports CLP-CBC-CGL      |
|       CLP        |       LP Solver        |              |         Coin LP (COIN)         |
|       CBC        |      MILP Solver       |              |   Coin Branch-and-Cut (COIN)   |
|       CGL        |        Library         |              | Cut Generation Library (COIN)  |
|     Symphony     |  MILP Solver/Library   |   10/10/15   |             (COIN)             |
|       dlib       |        Library         |              |              C++               |
|       GLPK       |         Solver         |   10/1/15    |      complicated install       |
|    JOptimizer    |        Library         |   12/20/14   |      minimization in Java      |
|      LEMON       |        Library         |              |           C++ (COIN)           |
|      Liger       | Integrated Environment |    9/7/13    |              C++               |
|      Midaco      |         Solver         |    5/1/13    |         many languages         |
|      NOMAD       |         Solver         |    3/1/15    |              C++               |
|     OpenMDAO     |   Library of Solvers   |   4/28/15    |             Python             |
|     OpenOpt      |        Library         |    8/8/15    |         Python Package         |
|       PPL        |         Solver         |    8/1/13    |             C/C++              |
|      Scilab      |        Language        |   3/31/15    |                                |
|  SciPy.Optimize  |        Library         |   7/23/15    |         Python Package         |
|       PuLP       |       LP Modeler       |    6/9/15    | Python to call solvers  (COIN) |
+------------------+------------------------+--------------+--------------------------------+

Proprietary Software Options

This list is not exhaustive.





In [6]:



    


# Table of Proprietary software for optimization
PropOpt = open(path+'NARSC__Optimizers__Proprietary.csv', "r")
PropTable = from_csv(PropOpt)
print PropTable
PropOpt.close()



    


















    






+-------------+-------------------+--------------------------------+
| Proprietary |        Type       |             Notes              |
+-------------+-------------------+--------------------------------+
|    Gurobi   |       Solver      |     free academic license      |
|    CPLEX    |       Solver      |     free academic license      |
|   TransCAD  |        GIS        | efficient cost matrix creation |
|    ArcGIS   |        GIS        |      Network Analyst Tool      |
|    MATLAB   |       Solver      |      Optimization Toolbox      |
|     AMPL    |   Solver/Modeler  |      30-day student trial      |
|    Kimeme   |       Solver      |        multi-objective         |
|    Lingo    |       Solver      |                                |
|    Maple    |       Solver      |                                |
|    MOSEK    |       Solver      |                                |
|     NAG     |       Solver      |                                |
|    SAS/OR   |       Solver      |                                |
|    TOMLAB   | Platform/Language |                                |
|    XPRESS   |       Solver      |                                |
+-------------+-------------------+--------------------------------+

PySAL 1.11.0 dev

Python Spatial Analysis Library

[https://www.pysal.readthedocs.org]

Sergio Rey at Arizona State University leads the PySAL project. [https://geoplan.asu.edu/people/sergio-j-rey]

"PySAL is an open source library of spatial analysis functions written in Python intended to support the development of high level applications. PySAL is open source under the BSD License." [https://pysal.readthedocs.org/en/latest/]

I will be only be demonstrating a portion of the functionality in PySAL.Network, but there are many other classes and functions for statistical spatial analysis within PySAL.





In [7]:



    


print dir(ps)



    


















    






['Box_Plot', 'DistanceBand', 'Equal_Interval', 'Fisher_Jenks', 'G', 'G_Local', 'Gamma', 'Geary', 'IOHandlers', 'Jenks_Caspall', 'Jenks_Caspall_Forced', 'Jenks_Caspall_Sampled', 'Join_Counts', 'K_classifiers', 'Kernel', 'LISA_Markov', 'MISSINGVALUE', 'Markov', 'Max_P_Classifier', 'Maximum_Breaks', 'Maxp', 'Maxp_LISA', 'Moran', 'Moran_BV', 'Moran_BV_matrix', 'Moran_Local', 'Natural_Breaks', 'Network', 'NetworkF', 'NetworkG', 'NetworkK', 'Percentiles', 'Quantiles', 'SpatialTau', 'Spatial_Markov', 'Std_Mean', 'Tau', 'Theil', 'TheilD', 'TheilDSim', 'Theta', 'User_Defined', 'W', 'WSP', '__builtins__', '__doc__', '__file__', '__name__', '__package__', '__path__', 'adaptive_kernelW', 'adaptive_kernelW_from_shapefile', 'base_path', 'bin', 'bin1d', 'binC', 'block_weights', 'buildContiguity', 'build_lattice_shapefile', 'cg', 'check_version', 'comb', 'common', 'config', 'config_path', 'core', 'datetime', 'directional', 'ergodic', 'esda', 'examples', 'full', 'gadf', 'hexLat2W', 'higher_order', 'higher_order_sp', 'inequality', 'json', 'kernelW', 'kernelW_from_shapefile', 'knnW', 'knnW_from_array', 'knnW_from_shapefile', 'lag_spatial', 'lat2W', 'min_threshold_dist_from_shapefile', 'network', 'open', 'order', 'os', 'pysal', 'quantile', 'queen_from_shapefile', 'query_yes_no', 'regime_weights', 'region', 'remap_ids', 'rook_from_shapefile', 'shimbel', 'spatial_dynamics', 'spreg', 'stable_release_date', 'sys', 'threshold_binaryW_from_array', 'threshold_binaryW_from_shapefile', 'threshold_continuousW_from_array', 'threshold_continuousW_from_shapefile', 'turn_off_check', 'urllib2', 'version', 'w_difference', 'w_intersection', 'w_subset', 'w_symmetric_difference', 'w_union', 'weight_convert', 'weights']

PySAL.Network

PySAL.Network was principally developed by Jay Laura at Arizona State Universty and the United States Geological Suvery. [https://geoplan.asu.edu/people/jay-laura]





In [8]:



    


print dir(ps.Network)



    


















    






['NetworkF', 'NetworkG', 'NetworkK', '__doc__', '__init__', '__module__', '_extractnetwork', '_newpoint_coords', '_round_sig', '_snap_to_edge', '_yieldneighbor', 'allneighbordistances', 'compute_distance_to_nodes', 'contiguityweights', 'count_per_edge', 'distancebandweights', 'enum_links_node', 'extractgraph', 'loadnetwork', 'nearestneighbordistances', 'node_distance_matrix', 'savenetwork', 'segment_edges', 'simulate_observations', 'snapobservations']

Gurobi 6.5

Relatively new company founded by optimization experts formerly at key positions with CPLEX.

[http://www.gurobi.com] [http://www.gurobi.com/company/about-gurobi]

Gurobipy

Python wrapper for Gurobi





In [9]:



    


print dir(gbp)



    


















    






['AttrConstClass', 'CallbackClass', 'Column', 'Constr', 'Env', 'ErrorConstClass', 'GRB', 'GurobiError', 'Iterable', 'LinExpr', 'Model', 'ParamClass', 'ParamConstClass', 'QConstr', 'QuadExpr', 'SOS', 'StatusConstClass', 'TempConstr', 'Var', '__builtins__', '__doc__', '__file__', '__name__', '__package__', '__path__', 'disposeDefaultEnv', 'fnmatch', 'gc', 'glob', 'gurobi', 'gurobipy', 'help', 'inspect', 'models', 'multidict', 'operator', 'os', 'paramHelp', 'quicksum', 'read', 'readParams', 'resetParams', 'setParam', 'sys', 'system', 'tuplelist', 'types', 'writeParams']

















In [10]:



    


# Manual Conceptual Model
My_Manual = IPd.Image(path+'/Manual_Diagram.png')
#print dir(My_Manual)
My_Manual



    


















    
Out[10]:

So why do things differently?

When solving spatial problems in a GIS, there are typically lots of files to keep track of and a lot of button clicking. This can lead to compounded errors at every stage in the process.

The following flowchart is a conceptual model of my current workflow. This workflow mimics the one described above, but the processes within are largely automated to maximize efficiency and reproducibility by minimizing the potential for human error in data handling, etc. This workflow will be demonstrated next.





In [11]:



    


# Automated Conceptual Model
My_Auto = IPd.Image(path+'/Auto_Diagram.png')
My_Auto



    


















    
Out[11]:

2.      Demonstration

Imports & Specs





In [12]:



    


from collections import OrderedDict
import networkx as nx
import numpy as np
import shapefile as shp
%pylab inline

print '********************************************************'
print ' | Platform Specs:                                    |'
print ' |  | OS X v', platform.mac_ver()[0],'                                 |'
print ' |  | Processor: ', platform.processor(), '                               |'
print ' |  | Machine: ', platform.machine(), '                               |'
print ' |  | Python: ', platform.python_version(), '                                 |'
print ' |  | PySAL: ', ps.version, '                              |'
print ' |  | Gurobi: ', gbp.gurobi.version(),'                             |'
print '********************************************************'
print ' |  | Date/Time --------- ', dt.datetime.now(), '|'
print '********************************************************'



    


















    






Populating the interactive namespace from numpy and matplotlib
********************************************************
 | Platform Specs:                                    |
 |  | OS X v 10.11.1                                  |
 |  | Processor:  i386                                |
 |  | Machine:  x86_64                                |
 |  | Python:  2.7.9                                  |
 |  | PySAL:  1.11.0dev                               |
 |  | Gurobi:  (6, 5, 0)                              |
********************************************************
 |  | Date/Time ---------  2015-11-14 16:12:17.181539 |
********************************************************

2.1      Data preparation and creation

2.1.1      Instantiate a network

First, I start out by instantiating a network by reading a shapefile in as a network object with PySAL. I use Waverly Hills, a smallish neighborhood in Tallahassee, FL where quite a few professors at Florida State University own homes. This shapefile was clipped from a US Census TIGER/Line file (US Census Bureau 2015). Additionally, there are some good topological challenges in Waverly Hills that make for good test cases, as seen below.

Note: Steps 3.1.3 through 3.1.7 are needed to create a toy dataset of aggregated population demand points (i.e. households) and service facilities (i.e. school bus stops). These steps are eliminated when reading in actual datasets.

Loop Road





In [13]:



    


# Avon Circle
Avon_Cir = IPd.Image(path+'/Avon.jpg')
Avon_Cir



    


















    
Out[13]:

Self-intersecting Road





In [14]:



    


# Millstream Road
Millstream_Rd = IPd.Image(path+'/Millstream.jpg')
Millstream_Rd



    


















    
Out[14]:
























In [15]:



    


# Instanitate network of Waverly Hills
ntw = ps.Network(path+'Waverly/Waverly.shp')

print 'Nodes in network:       ', len(ntw.nodes)
print 'Edges in network:       ', len(ntw.edges)
print 'Graph edges in network: ', len(ntw.graphedges)



    


















    






Nodes in network:        1086
Edges in network:        1122
Graph edges in network:  152

2.1.2      Instantiate all graphs to be drawn

Next, I instantiate the graphs to be drawn for visualization. This includes a representation of edges comprising the road network and nodes for the demand and service nodes. NetworkX is utilized for this step.





In [16]:



    


# 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()

## Optimized Locations
# Set Cover
setcover_g = nx.Graph()

# p-Median
median_g = nx.Graph()

# p-Center
center_g = nx.Graph()

# p-Center Manual
center_g_man = nx.Graph()

2.1.3      Create Bounding Box from Waverly.shp





In [17]:



    


# Instantiate the Waverly Hills neighborhood shapefile
shp_W = ps.open(path+'Waverly/Waverly.shp')

# Create a bounding box of the shapefile
shp_W.bbox



    


















    
Out[17]:








[-84.27732299999998, 30.470055856478194, -84.25284173452727, 30.49997682517085]

In order to create simulated demand and service nodes for this toy problem I create a bounding latitude and longitude box from the extremes of the 'Waverly.shp' of the Waverly Hills neighborhood. The total area is roughly 1.5 square miles.

2.1.4      Create Numpy arrays of random floats within a bounding box

Within the bounding box I create lists of 100 random numbers for the latitude and longitude of demand points (clients), then 15 random numbers for the latitude and longitude of the service points with Numpy.





In [18]:



    


# Client latitude
np.random.seed(850)
lat_client = np.random.uniform(shp_W.bbox[0], shp_W.bbox[2], 100)

np.random.seed(352)
# Client longitude
lon_client = np.random.uniform(shp_W.bbox[1], shp_W.bbox[3], 100)

np.random.seed(904)
# Service latitude
lat_service = np.random.uniform(shp_W.bbox[0], shp_W.bbox[2], 15)

np.random.seed(407)
# Service longitude
lon_service = np.random.uniform(shp_W.bbox[1], shp_W.bbox[3], 15)

2.1.5      Zip the latitude and longitude lists together

Next, these four lists are zipped together to form two: a list of client lat/lons and a list and service lat/lons





In [19]:



    


# Client lat/lon coordinates
rand_coords_client = map(list, zip(lat_client, lon_client))

# Service lat/lon coordinates
rand_coords_service = map(list, zip(lat_service, lon_service))

2.1.6      Create empty random points dictionaries

Dictionaries are instantiated to house the point location data for the simulated client and service facilities.





In [20]:



    


# Empty Clients dictionary
points_client = {}

# Empty Service dictionary
points_service = {}

2.1.7      Fill dictionaries of random roints

The dictionaries that have just been created are now filled with the ids and lat/lons of the points.





In [21]:



    


# CLIENT with {id: [lon, lat], } format
for idx, coords in enumerate(rand_coords_client):
    GRAPH_client.add_node(idx)
    points_client[idx] = coords
    GRAPH_client.node[idx] = coords

# SERVICE with {id: [lon, lat], } format  
for idx, coords in enumerate(rand_coords_service):
    GRAPH_service.add_node(idx)
    points_service[idx] = coords
    GRAPH_service.node[idx] = coords

2.1.8      Draw roads, simplified network, and random client & service nodes

The following is a visualization of the data I have just created and the roads in Waverly Hills using matplotlib. The roads and nodes are derived for network object created in [3.1.1].

Nodes in pink represent the actual intersections and the roads are represented in light red. The simplified network (connecting graph nodes) is gray and the remaining graph nodes (where a choice in forward destination is possible) are red. These are necessary for determination of a shortest path and may greatly would diminish the computational time. The client nodes are represented in blue and the service nodes are represented in cyan.





In [22]:



    


print dt.datetime.now()

#Instantiate Figure
figsize(10,11)

#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=20, 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=75, alpha=1, node_color='c')

# Legend (Ordered Dictionary) from collections
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,
       scatterpoints=1)

# Title
title('Waverly Hills\nTallahassee, 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.49, 0.0, 0.005, width=.0003, head_width=0.001, 
          head_length=0.0015, fc='k', ec='k',alpha=0.75,)
annotate('N', xy=(-84.2815, 30.498), fontstyle='italic', fontsize='xx-large',
            fontweight='heavy', alpha=0.75)



    


















    






2015-11-14 16:15:54.657268










    
Out[22]:








<matplotlib.text.Annotation at 0x1105cb990>

2.1.9      Create weights at nodes and sum

In order to solve the population weighted problem in this demonstration, the p-median problem, I will Numpy to generate a list of 100 random integers [1-4] to represent population at each client node. In this case we can think of the population weight as number of children in the household.





In [23]:



    


# Client Weights for demand
np.random.seed(850)
Ai = np.random.randint(1, 5, len(rand_coords_client))
Ai = Ai.reshape(len(Ai),1)

# Sum of Weights (Total Demand) 
AiSum = np.sum(Ai)

2.1.10      Instantiate client and service shapefiles

Using the python package shapefile, I instantiate .shp files of the client and service location data that was created above.





In [24]:



    


# Client Shapefile
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(path+'Simulated/RandomPoints_CLIENT')

#Service Shapefile
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(path+'Simulated/RandomPoints_SERVICE')

2.1.11      Snap observations to ntw

Next, I snap the client and service observations to the network with the ntw.snapobservations method within the instantiated network object. This allows for the creation of a distance matrix later.





In [25]:



    


t1 = time.time()
Snap_C = ntw.snapobservations(path+'Simulated/RandomPoints_CLIENT.shp', 
                     'Rand_Points_CLIENT', attribute=True)
Snap_S = ntw.snapobservations(path+'Simulated/RandomPoints_SERVICE.shp', 
                     'Rand_Points_SERVICE', attribute=True)
print round(time.time()-t1, 4), 'seconds'



    


















    






2.9429 seconds

2.1.12      Draw ntw, snapped coords, & random coords

When the graph is redrawn we can the client and service location in their original locations represented by the smaller points, and the network-snapped location represented with larger points.





In [26]:



    


# Instantiate Figure
figsize(10,11)

# 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=75, 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=75, 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 Service 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, scatterpoints=1)

# 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.49, 0.0, 0.005, width=.0003, head_width=0.001, 
          head_length=0.0015, fc='k', ec='k',alpha=0.75,)
annotate('N', xy=(-84.2815, 30.498), fontstyle='italic', fontsize='xx-large',
            fontweight='heavy', alpha=0.75)



    


















    
Out[26]:








<matplotlib.text.Annotation at 0x11156be10>

2.1.13      Instantiate snapped service shapefile

This code chunk is key for representing the results of the optimization problems accurately. I create a shapefile of the snapped service variable names and their coordinates to be used for comparison later.





In [27]:



    


# Create Lat & Lon dictionaries of the snapped service locations
lat_snapped = []
lon_snapped = []
for i,j in ntw.pointpatterns['Rand_Points_SERVICE'].snapped_coordinates.iteritems():
    lat_snapped.append(j[0]) 
    lon_snapped.append(j[1])

# Snapped Service Shapefile
service_SNAP = shp.Writer(shp.POINT)
# Add Points
for i,j in ntw.pointpatterns['Rand_Points_SERVICE'].snapped_coordinates.iteritems():
    service_SNAP.point(j[0],j[1])
# Add Fields
service_SNAP.field('y_ID')
service_SNAP.field('x_ID')
service_SNAP.field('LAT')
service_SNAP.field('LON')
counter = 0
for i in range(len(ntw.pointpatterns['Rand_Points_SERVICE'].snapped_coordinates)):
    counter = counter + 1
    service_SNAP.record('y' + str(counter), 'x' + str(counter), lat_snapped[i], lon_snapped[i])
# Save Shapefile    
service_SNAP.save(path+'Snapped/SERVICE_Snapped')

2.1.14      Create distance matrices

Finally, I create a client to service distance matrix with ntw.allneighbordistances() from the instantiated network object of Waverly Hills [Waverly.shp]. The matrix is 100x15 and the values are converted from decimal degrees to miles.

A current limitation of this is run time which is hindered by the geographic size of the network. When tested with much larger client and service point sets the Dijkstra algorithm ran in similar time. This demonstrates the limitation of calculating a shortest-path with the actual road network. As a potential improvement, the Dijkstra may be refactored and calculated on the simplified graph to see decreases in run time. The run time should decrease because there are far fewer nodes and edges in the simplified graph.





In [28]:



    


t1 = time.time()

# Define Client to Service Matrix Function
def c_s_matrix():
    global All_Dist_MILES
    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 

# Call Client to Service Matrix Function
c_s_matrix()    

seconds = round(time.time()-t1, 4)    
print seconds, 'seconds'
print 'Client to Service Matrix Shape --> ', All_Dist_MILES.shape



    


















    






6.7371 seconds
Client to Service Matrix Shape -->  (100, 15)

2.2      Mathematical Optimization

2.2.1a      Set Cover test [S = 1.25 miles]

The set cover location problem, also known as the SCLP, is a fundamental facility location problem (Beasley 1987; Daskin 1995; Beasley and Chu 1996). The objective of the SCLP is to minimize the number of service facilities to site while achieving complete (100%) coverage of client nodes.

Add Revelle Current etc...

Minimize

         $\displaystyle\sum_{j=1}^n d_j x_j$

Subject to

         $\displaystyle\sum_{j=1}^n a_{ij} x_j \geq 1,$            $i = 1 \in n$

         $\displaystyle x_j \in (0,1)$                $j = 1 \in m$

where

         − $i$ = a specific origin

         − $j$ = a specific destination

         − $n$ = the set of origins

         − $m$ = the set of destinations

         − $x_i$ = the decision variable at each node in the matrix

         − $d_{ij}$ = distance from $i$th origin to $j$th destination

         − $a_{ij}$ = binary matrix describing the coverage of each node

          where

                     − $a_{ij}$ = $1 \forall i, j \ni d_{ij} \geq S$ ($S$ is user defined)

                     − $a_{ij} $ = $0$ otherwise

Adapted from:

  • Daskin, M. S. 1995. Network and Discrete Location: Models, Algorithms, and Applications. Hoboken, NJ, USA: John Wiley & Sons, Inc.




In [29]:



    


# Define the Set Cover function
def gbpSCLP():    
    
    t1 = time.time()
    
    # Define Global Records Variable
    global NEW_Records_SCLP
    
    #     1. Read In Data
    # Cost Matrix
    Cij = All_Dist_MILES
    # Create Aij: Determine Aij (nodes within S)
    # S --> 1 = served; 0 = unserved
    S = 1.25
    # Aij
    Aij = []
    for i in np.nditer(Cij):
        if i <= S:
            outtext = 1
        else:
            outtext = 0
        Aij.append(outtext)
    rows, cols = Cij.shape
    
    # Transform Aij into an array and resphape to match with Cij
    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 = gbp.Model(" -- SCLP -- ")
    
    # Set MIP Focus to 2 for optimality
    gbp.setParam('MIPFocus', 2)

    # Add Service Decision Variables (j)
    serv_var = []
    for dest in service_nodes:
        serv_var.append(mSCLP.addVar(vtype=gbp.GRB.BINARY,
                                    ub = 1,
                                    name='x'+str(dest+1)))
    # Update Model Variables
    mSCLP.update()       

    #     3. Set Objective Function
    mSCLP.setObjective(gbp.quicksum(serv_var[dest] 
                                for dest in service_nodes), 
                                gbp.GRB.MINIMIZE)

    #     4. Add Constraints 
    # Add Coverage Constraints  
    for orig in client_nodes:
        mSCLP.addConstr(gbp.quicksum(Aij[orig][dest]*serv_var[dest] 
                                for dest in service_nodes) >= 1)        

    #     5. Optimize and Print Results
    # Solve
    try:
        mSCLP.optimize()
    except Exception as e:
        exc_type, exc_obj, exc_tb = sys.exc_info()
        fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
        print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
        print '   ################################################################'
        print ' < ISSUE : ', e, ' >'
        print '   ################################################################'
    
    # Write LP
    mSCLP.write(path+'LP_Files/WaverlySCLP.lp')
    t2 = time.time()-t1
    
    # Record and Display Results
    print '\n*****************************************************************************************'
    selected = []
    dbf1 = ps.open(path+'Snapped/SERVICE_Snapped.dbf')
    NEW_Records_SCLP = []
    for v in mSCLP.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.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 100% coverage of clients --- ', len(selected)
    print '    | Real Time to Optimize (sec.) ----------------- ', t2
    print '    | Date/Time ------------------------------------ ', dt.datetime.now()
    print '*****************************************************************************************'
    print ' -- Set Cover Location Problem -- '
    
# Call SCLP Function
try:
    gbpSCLP()
    print '\nJames Gaboardi, 2015'
except Exception as e:
    exc_type, exc_obj, exc_tb = sys.exc_info()
    fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
    print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
    print '   ################################################################'
    print ' < ISSUE : ', e, ' >'
    print '   ################################################################'



    


















    






Changed value of parameter MIPFocus to 2
   Prev: 0  Min: 0  Max: 3  Default: 0
Optimize a model with 100 rows, 15 columns and 651 nonzeros
Coefficient statistics:
  Matrix range    [1e+00, 1e+00]
  Objective range [1e+00, 1e+00]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 1e+00]
Found heuristic solution: objective 3
Presolve removed 100 rows and 15 columns
Presolve time: 0.01s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.01 seconds
Thread count was 1 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective 3.000000000000e+00, best bound 3.000000000000e+00, gap 0.0%

*****************************************************************************************
    |                                                 x12
    |                                                 x13
    |                                                 x15
    | Selected Facility Locations ------------------  ^^^^ 
    | Coverage (S) in miles ------------------------  1.25
    | Client Nodes ---------------------------------  100
    | Facilities needed 100% coverage of clients ---  3
    | Real Time to Optimize (sec.) -----------------  0.100801944733
    | Date/Time ------------------------------------  2015-11-14 16:20:00.075939
*****************************************************************************************
 -- Set Cover Location Problem -- 

James Gaboardi, 2015

2.2.1b     Instantiate selected Set Cover location shapefile

Following optimization, I created a new shapefile of the selected candidate service locations.





In [30]:



    


# Define SCLP shapefile function
def Create_SCLP():
    try:
        # Instantiate SCLP shapefile
        SHP_SetCover = shp.Writer(shp.POINT)
        # Add Points
        for idy,idx,x,y in NEW_Records_SCLP:
            SHP_SetCover.point(float(x), float(y))
        # Add Fields
        SHP_SetCover.field('y_ID')
        SHP_SetCover.field('x_ID')
        SHP_SetCover.field('LAT')
        SHP_SetCover.field('LON')
        # Add Records
        for idy,idx,x,y in NEW_Records_SCLP:
            SHP_SetCover.record(idy,idx,x,y)
        # Save Shapefile    
        SHP_SetCover.save(path+'Results/Selected_Locations_SetCover')
    except Exception as e:
        exc_type, exc_obj, exc_tb = sys.exc_info()
        fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
        print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
        print '   ################################################################'
        print ' < ISSUE : ', e, ' >'
        print '   ################################################################'

# Call SCLP shapefile function
try:
    Create_SCLP()
except Exception as e:
    exc_type, exc_obj, exc_tb = sys.exc_info()
    fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
    print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
    print '   ################################################################'
    print ' < ISSUE : ', e, ' >'
    print '   ################################################################'

2.2.2a     p-Median test [p = 2]

The p-median problem (PMP), also known as the minisum problem, was first proposed by Hakimi in the 1960s and has been utilized in wide-ranging research topics such as Mathematics and Emergency Management since then (Hakimi 1964; Teitz and Bart 1968; Tansel, Francis, and Lowe 1983; Daskin 1995; Horner and Widener 2010; Kalcsics et al. 2014). This model sites p facilities with an objective of minimizing the total travel cost by siting facilities nearest to the greatest amount of demand. Weighting decision variables with population, for example, does this. Efficiency (lowest travel cost) is the key to the PMP.

Minimize

         $\displaystyle\sum_{i \in 1}^n\sum_{j\in 1}^m a_i c_{ij} x_{ij}$

Subject to

         $\displaystyle\sum_{j\in m} x_{ij} = 1 ,$         $\forall i \in n$

         $\displaystyle\sum_{i \in n} y_j = p$

         $x_{ij} - y_j \geq 0,$        $\forall i \in n, j \in m$

         $x_{ij}, y_j \in (0,1)$      $\forall i \in n , j \in m$

where

         − $i$ = a specific origin

         − $j$ = a specific destination

         − $n$ = the set of origins

         − $m$ = the set of destinations

         − $a_i$ = weight at each node

         − $c_{ij}$ = travel costs between nodes

         − $x_{ij}$ = the decision variable at each node in the matrix

         − $y_j$ = nodes chosen as service facilities

         − $p$ = the number of facilities to be sited

Adapted from:

  • Daskin, M. S. 1995. Network and Discrete Location: Models, Algorithms, and Applications. Hoboken, NJ, USA: John Wiley & Sons, Inc.




In [31]:



    


# define p-Median function
def gbpPMP():
    
    t1 = time.time()
    
    # Define Global Variables
    global Ai
    global AiSum
    global NEW_Records_PMP 
    
    #     1. Data
    # Demand
    Ai = Ai
    # Demand Sum
    AiSum = AiSum
    # Travel Costs
    Cij = All_Dist_MILES
    # Weighted Costs
    Sij = Ai * Cij
    # Total Client and Service nodes
    client_nodes = range(len(Sij))
    service_nodes = range(len(Sij[0]))
    
    #     2. Create Model, Set MIP Focus, Add Variables, & Update Model
    # Instantiate Model
    mPMP = gbp.Model(' -- p-Median -- ')
    
    # Set MIP focus to 'Optimal'
    gbp.setParam('MIPFocus', 2)

    # Add Client Decision Variables (iXj)
    client_var = []
    for orig in client_nodes:
        client_var.append([])
        for dest in service_nodes:
            client_var[orig].append(mPMP.addVar(vtype=gbp.GRB.BINARY, 
                                                obj=Sij[orig][dest], 
                                                name='x'+str(orig+1)+'_'+str(dest+1)))
    
    # Add Service Decision Variables (j)
    serv_var = []
    for dest in service_nodes:
        serv_var.append([])
        serv_var[dest].append(mPMP.addVar(vtype=gbp.GRB.BINARY, 
                                        name='y'+str(dest+1)))
    
    # Update the model
    mPMP.update()
    
    #     3. Set Objective Function
    mPMP.setObjective(gbp.quicksum(Sij[orig][dest]*client_var[orig][dest] 
                            for orig in client_nodes for dest in service_nodes), 
                            gbp.GRB.MINIMIZE)
    
    
    #     4. Add Constraints
    # Assignment Constraints
    for orig in client_nodes:
        mPMP.addConstr(gbp.quicksum(client_var[orig][dest] 
                            for dest in service_nodes) == 1)
    # Opening Constraints
    for orig in service_nodes:
        for dest in client_nodes:
            mPMP.addConstr((serv_var[orig][0] - client_var[dest][orig] >= 0))
    
    # Facility Constraint
    mPMP.addConstr(gbp.quicksum(serv_var[dest][0] for dest in service_nodes) == 2)
    
    #     5. Optimize and Print Results
    # Solve
    try:
        mPMP.optimize()
    except Exception as e:
        exc_type, exc_obj, exc_tb = sys.exc_info()
        fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
        print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
        print '   ################################################################'
        print ' < ISSUE : ', e, ' >'
        print '   ################################################################'
    
    # Write LP
    mPMP.write(path+'LP_Files/WaverlyPMP.lp')
    t2 = time.time()-t1
    
    # Record and Display Results
    print '\n*************************************************************************'
    selected = []
    dbf1 = ps.open(path+'Snapped/SERVICE_Snapped.dbf')
    NEW_Records_PMP = []
    for v in mPMP.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.append(x)
                else:
                    pass
            print '    |                                            ', var
    print '    | Selected Facility Locations --------------  ^^^^ '
    print '    | Candidate Facilities [p] ----------------- ', len(selected)
    val = mPMP.objVal
    print '    | Objective Value (miles) ------------------ ', val
    avg = float(mPMP.objVal)/float(AiSum)
    print '    | Avg. Value / Client (miles) -------------- ', avg
    print '    | Real Time to Optimize (sec.) ------------- ', t2
    print '    | Date/Time -------------------------------- ', dt.datetime.now()
    print '*************************************************************************'
    print ' -- The p-Median Problem -- '
    
# Call p-Median Function   
try:
    gbpPMP()
    print '\nJames Gaboardi, 2015'
except Exception as e:
    exc_type, exc_obj, exc_tb = sys.exc_info()
    fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
    print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
    print '   ################################################################'
    print ' < ISSUE : ', e, ' >'
    print '   ################################################################'



    


















    






Parameter MIPFocus unchanged
   Value: 2  Min: 0  Max: 3  Default: 0
Optimize a model with 1601 rows, 1515 columns and 4515 nonzeros
Coefficient statistics:
  Matrix range    [1e+00, 1e+00]
  Objective range [4e-03, 1e+01]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 2e+00]
Presolve time: 0.01s
Presolved: 1601 rows, 1515 columns, 4515 nonzeros
Variable types: 0 continuous, 1515 integer (1515 binary)
Found heuristic solution: objective 159.8861180
Presolved: 1601 rows, 1515 columns, 4515 nonzeros


Root relaxation: objective 1.567066e+02, 432 iterations, 0.01 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0     156.7065613  156.70656  0.00%     -    0s

Explored 0 nodes (432 simplex iterations) in 0.05 seconds
Thread count was 4 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective 1.567065613457e+02, best bound 1.567065613457e+02, gap 0.0%

*************************************************************************
    |                                             y12
    |                                             y15
    | Selected Facility Locations --------------  ^^^^ 
    | Candidate Facilities [p] -----------------  2
    | Objective Value (miles) ------------------  156.706561346
    | Avg. Value / Client (miles) --------------  0.631881295749
    | Real Time to Optimize (sec.) -------------  0.17392206192
    | Date/Time --------------------------------  2015-11-14 16:21:17.361723
*************************************************************************
 -- The p-Median Problem -- 

James Gaboardi, 2015

2.2.2b     Instantiate Selected p-median shapefile

Following optimization, I created a new shapefile of the selected candidate service locations.





In [32]:



    


# define PMP shapefile function
def Create_PMP():
    try:
        # Instantiate Shapefile
        SHP_Median = shp.Writer(shp.POINT)
        # Add Points
        for idy,idx,x,y in NEW_Records_PMP:
            SHP_Median.point(float(x), float(y))
        # Add Fields
        SHP_Median.field('y_ID')
        SHP_Median.field('x_ID')
        SHP_Median.field('LAT')
        SHP_Median.field('LON')
        # Add Records
        for idy,idx,x,y in NEW_Records_PMP:
            SHP_Median.record(idy,idx,x,y)
        # Save Shapefile    
        SHP_Median.save(path+'Results/Selected_Locations_Pmedian')
    except Exception as e:
        exc_type, exc_obj, exc_tb = sys.exc_info()
        fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
        print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
        print '   ################################################################'
        print ' < ISSUE : ', e, ' >'
        print '   ################################################################'

# Call PMP shapefile function
try:
    Create_PMP()
except Exception as e:
    exc_type, exc_obj, exc_tb = sys.exc_info()
    fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
    print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
    print '   ################################################################'
    print ' < ISSUE : ', e, ' >'
    print '   ################################################################'

2.2.3a      p-Center test [p = 2]

The p-center problem (PCP), also referred to as the minimax problem, sites facilities with a focus on equity and has been widely utilized in research as a counter-measure to the efficiency-based p-median problem. The objective of the PCP is the minimization of the maximum (worst-case) travel cost from client nodes to service facilities (Minieka 1970; Shier 1977; Tansel, Francis, and Lowe 1983; Horner and Widener 2010).

Minimize

         $W$

Subject to

         $\displaystyle\sum_{j\in m} x_{ij} = 1,$              $\forall i \in n$

         $\displaystyle\sum_{i \in n} y_j = p$

         $x_{ij} - y_j \geq 0,$             $\forall i\in n, j \in m$

         $\displaystyle W \geq \sum_{j \in m} c_{ij} x_{ij}$          $\forall i \in n$

         $x_{ij}, y_j \in (0,1)$           $\forall i \in n, j \in m$

where

         − $W$ = the maximum travel cost between client and service nodes

         − $i$ = a specific origin

         − $j$ = a specific destination

         − $n$ = the set of origins

         − $m$ = the set of destinations

         − $a_i$ = weight at each node

         − $c_{ij}$ = travel costs between nodes

         − $x_{ij}$ = the decision variable at each node in the matrix

         − $y_j$ = nodes chosen as service facilities

         − $p$ = the number of facilities to be sited

Adapted from:

  • Daskin, M. S. 1995. Network and Discrete Location: Models, Algorithms, and Applications. Hoboken, NJ, USA: John Wiley & Sons, Inc.




In [33]:



    


# define PCP shapefile function
def gbpPCP():
    
    t1 = time.time()
    
    # Define Global Variables
    global Cij
    global mPCP
    global dbf1
    global selected_PCP
    global NEW_Records_PCP
    
    #     1. Data
    Cij = All_Dist_MILES
    
    # Total Client and Service nodes
    client_nodes = range(len(Cij))
    service_nodes = range(len(Cij[0]))
    
    #    2. Create Model, Set MIP Focus, Add Variables, & Update Model
    # Instantiate Model
    mPCP = gbp.Model(' -- P-Center -- ')
    # Set MIP focus to 'Optimal'
    gbp.setParam('MIPFocus', 2)
    
    # Add Client Decision Variables (iXj)
    client_var = []
    for orig in client_nodes:
        client_var.append([])
        for dest in service_nodes:
            client_var[orig].append(mPCP.addVar(vtype=gbp.GRB.BINARY, 
                                                obj=Cij[orig][dest], 
                                                name='x'+str(orig+1)+'_'+str(dest+1)))
    # Add Service Decision Variables (j)
    serv_var = []
    for dest in service_nodes:
        serv_var.append([])
        serv_var[dest].append(mPCP.addVar(vtype=gbp.GRB.BINARY, 
                                        name='y'+str(dest+1)))
    
    # Add the Maximum travel cost variable
    W = mPCP.addVar(vtype=gbp.GRB.CONTINUOUS,
                name='W')        
    
    # Update the model
    mPCP.update()  
    
    #     3. Set the Objective function
    mPCP.setObjective(W, gbp.GRB.MINIMIZE)
    
    #     4. Add Constraints
    # Add Assignment Constraints
    for orig in client_nodes:
        mPCP.addConstr(gbp.quicksum(client_var[orig][dest] 
                            for dest in service_nodes) == 1)
    
    # Add Opening constraints
    for orig in service_nodes:
        for dest in client_nodes:
            mPCP.addConstr((serv_var[orig][0] - client_var[dest][orig] >= 0))
    
    # Add Facility Constraints
    mPCP.addConstr(gbp.quicksum(serv_var[dest][0] for dest in service_nodes) == 2)
    
    # Add Maximum travel cost constraints
    for orig in client_nodes:
        mPCP.addConstr(gbp.quicksum(Cij[orig][dest]*client_var[orig][dest]
                            for dest in service_nodes) - W <= 0)
    
    #     5. Optimize and Print Results
    # Solve
    try:
        mPCP.optimize()
    except Exception as e:
        exc_type, exc_obj, exc_tb = sys.exc_info()
        fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
        print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
        print '   ################################################################'
        print ' < ISSUE : ', e, ' >'
        print '   ################################################################'  
    
    # Write LP
    mPCP.write(path+'/LP_Files/WaverlyPCP.lp')
    t2 = time.time()-t1
    print '\n*************************************************************************'
    
    # Record and Display Results
    selected_PCP = []
    dbf1 = ps.open(path+'Snapped/SERVICE_Snapped.dbf')
    NEW_Records_PCP = []
    for v in mPCP.getVars():
        if 'x' in v.VarName:
            pass
        elif 'W' in v.VarName:
            pass
        elif v.x > 0:
            var = '%s' % v.VarName
            selected_PCP.append(var)
            for i in range(dbf1.n_records):
                if var in dbf1.read_record(i):
                    x = dbf1.read_record(i)
                    NEW_Records_PCP.append(x)
                else:
                    pass
            print '    |                                            ', var,  '         '
    print '    | Selected Facility Locations -------------- ^^^^ ',  '         '
    print '    | Candidate Facilities [p] ----------------- ', len(selected_PCP), '           '
    print '    | Objective Value (miles) ------------------ ', mPCP.objVal, '     '
    print '    | Real Time to Optimize (sec.) ------------- ', t2
    print '    | Date/Time -------------------------------- ', dt.datetime.now()
    print '*************************************************************************'
    print ' -- The p-Center Problem -- '
    
# Call p-Center Function   
try:
    gbpPCP()
    print '\nJames Gaboardi, 2015'
except Exception as e:
    exc_type, exc_obj, exc_tb = sys.exc_info()
    fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
    print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
    print '   ################################################################'
    print ' < ISSUE : ', e, ' >'
    print '   ################################################################'



    


















    






Parameter MIPFocus unchanged
   Value: 2  Min: 0  Max: 3  Default: 0
Optimize a model with 1701 rows, 1516 columns and 6111 nonzeros
Coefficient statistics:
  Matrix range    [1e-03, 3e+00]
  Objective range [1e+00, 1e+00]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 2e+00]
Presolve time: 0.02s
Presolved: 1701 rows, 1516 columns, 6015 nonzeros
Variable types: 1 continuous, 1515 integer (1515 binary)
Found heuristic solution: objective 1.7394879
Found heuristic solution: objective 1.7390156
Found heuristic solution: objective 1.7295779
Presolve removed 112 rows and 112 columns
Presolved: 1589 rows, 1404 columns, 5567 nonzeros


Root relaxation: objective 1.227915e+00, 958 iterations, 0.02 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    1.22791    0  315    1.72958    1.22791  29.0%     -    0s
H    0     0                       1.5277094    1.22791  19.6%     -    0s
*    0     0               0       1.3697645    1.36976  0.00%     -    0s

Cutting planes:
  Gomory: 4
  Implied bound: 13
  Zero half: 120

Explored 0 nodes (4072 simplex iterations) in 0.20 seconds
Thread count was 4 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective 1.369764530495e+00, best bound 1.369764530495e+00, gap 0.0%

*************************************************************************
    |                                             y3          
    |                                             y15          
    | Selected Facility Locations -------------- ^^^^           
    | Candidate Facilities [p] -----------------  2            
    | Objective Value (miles) ------------------  1.3697645305      
    | Real Time to Optimize (sec.) -------------  0.298266887665
    | Date/Time --------------------------------  2015-11-14 16:21:55.071883
*************************************************************************
 -- The p-Center Problem -- 

James Gaboardi, 2015

2.2.3b      Instantiate selected p-Center shapefile

Following optimization, I created a new shapefile of the selected candidate service locations.





In [34]:



    


# define PCP shapefile function
def Create_PCP():
    try:
        # Instantiate Shapefile
        SHP_Center = shp.Writer(shp.POINT)
        # Add Points
        for idy,idx,x,y in NEW_Records_PCP:
            SHP_Center.point(float(x), float(y))
        # Add Fields
        SHP_Center.field('y_ID')
        SHP_Center.field('x_ID')
        SHP_Center.field('LAT')
        SHP_Center.field('LON')
        # Add Records
        for idy,idx,x,y in NEW_Records_PCP:
            SHP_Center.record(idy,idx,x,y)
        # Save Shapefile    
        SHP_Center.save(path+'Results/Selected_Locations_Pcenter')
    except Exception as e:
        exc_type, exc_obj, exc_tb = sys.exc_info()
        fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
        print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
        print '   ################################################################'
        print ' < ISSUE : ', e, ' >'
        print '   ################################################################'

# Call PCP shapefile function
try:
    Create_PCP()
except Exception as e:
    exc_type, exc_obj, exc_tb = sys.exc_info()
    fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
    print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
    print '   ################################################################'
    print ' < ISSUE : ', e, ' >'
    print '   ################################################################'

2.2.4a      Manual p-Center .lp File Creation

The following is a python script I developed to read in the data from a cost matrix and produce a linear/integer programming file.





In [35]:



    


# p-Center Facility Location Problem
# This script creates a linear programming file to be read into an optimizer.
'''
GNU LESSER GENERAL PUBLIC LICENSE
                       Version 3, 29 June 2007

 Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
 Everyone is permitted to copy and distribute verbatim copies
 of this license document, but changing it is not allowed.
'''
# Developed by:  James D. Gaboardi, MSGIS
#                03/2015
#                James Gaboardi

#   Terminology & General Background for Facility Location and Summation Notation:

#   *        The objective of the p-center Facility Location Problem is to minimize the maximum cost 
#            of travel between service facilities and clients on a network.

#   *   [i] - a specific origin
#   *   [j] - a specifc destination
#   *   [n] - the set of origins
#   *   [m] - the set of destinations
#   *   [Cij] - travel costs between nodes
#   *   [W] - the maximum travel costs between service facilities and clients 
#   *   [x#_#] - the client decision variable
#   *   [y#] - the service decision variable
#   *   [p] - the number of facilities to be sited

#    DEFINED FUNCTIONS
# Assignment Constraints
def get_assignment_constraints():
    outtext = ' '
    for i in range(1,rows+1):
        temp = ' '
        for j in range(1,cols+1):
            temp += 'x' + str(i) + '_' + str(j) + ' + '
        outtext += temp[:-2] + '= 1\n'
    return outtext

# Facility Constraint  
def get_p_facilities():
    outtext = ''
    for i in range(1, cols+1):
        temp = ''
        temp += 'y' + str(i)
        outtext += temp + ' + '
    outtext = ' ' + outtext[:-2] + '= 2\n'
    return outtext

# Opening Constraints
def get_opening_constraints_p_center():
    outtext = ' '
    for i in range(1, cols+1):
        for j in range(1, rows+1):
            outtext += ' - x' + str(j) + '_' + str(i) + ' + ' + 'y' + str(i) +  ' >= 0\n'
    return outtext

# Maximum Cost Constraints
def get_max_cost():
    outtext = ''
    for i in range(rows):
        temp = '  '
        for j in range(cols):
            temp += str(Cij[i,j]) + ' x' + str(i+1) + '_' + str(j+1) + ' + '
        outtext += temp[:-2] + '- W <= 0\n'
    return outtext

# Declaration of Bounds
def get_bounds_allocation():
    outtext = ' '
    for i in range(rows):
        temp = ''
        for j in range(cols):
            temp += ' 0 <= x' + str(i+1) + '_' + str(j+1) + ' <= 1\n'
        outtext += temp    
    return outtext

def get_bounds_facility():
    outtext = ''
    for i in range(cols):
        outtext += ' 0 <= y' + str(i+1) + ' <= 1\n'
    return outtext

# Declaration of Decision Variables (form can be: Binary, Integer, etc.)
def get_decision_variables_p_center():
    outtext = ' '
    for i in range(1, rows+1):
        temp = ''
        for j in range(1, cols+1):
            temp += 'x' + str(i) + '_' + str(j) + ' '
        outtext += temp
    return outtext
    
def get_facility_decision_variables_p_center():  
    outtext = ''
    for i in range (1, cols+1):
        outtext += 'y' + str(i) + ' '
    return outtext    

#    DATA READS & VARIABLE DECLARATION
Cij = All_Dist_MILES
rows,cols = Cij.shape

#    START TEXT FOR .lp FILE
# Declaration of Objective Function
text = 'Minimize\n'
text += ' obj: W\n'
# Declaration of Constraints
text += 'Subject To\n'                    
text += get_assignment_constraints()
text += get_p_facilities()
text += get_opening_constraints_p_center()
text += get_max_cost()
# Declaration of Bounds
text += 'Bounds\n'
text += get_bounds_allocation()
text += get_bounds_facility()
# Declaration of Decision Variables form: Binaries
text += 'Binaries\n'
text += get_decision_variables_p_center()
text += get_facility_decision_variables_p_center()
text += '\n'
text += 'End\n'
text += "'''\n"
text += "James Gaboardi, 2015"
                
#   CREATE & WRITE .lp FILE TO DISK
# Fill path name  --  File name must not have spaces.
outfile = open(path+'LP_Files/pCenter_Manual.lp', 'w')
outfile.write(text)
outfile.close()

2.2.4b      Manual p-Center .lp File Solved

The .lp file created above is then read into Gurobi through gurobipy and solved.





In [36]:



    


# Define manual LP read PCP Function
def Manual_LP_PCP():
    
    global Cij
    global manualPCP
    global dbf1
    global selected_PCP_manual
    global NEW_Records_PCP_Man
    
    t1 = time.time()

    # Instantiate Optimization model from .lp file
    manualPCP = gbp.read(path+'LP_Files/pCenter_Manual.lp')

    # Set MIP focus to 'Optimal'
    gbp.setParam('MIPFocus', 2)

    # Solve
    try:
        manualPCP.optimize()
    except Exception as e:
        exc_type, exc_obj, exc_tb = sys.exc_info()
        fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
        print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
        print '   ################################################################'
        print ' < ISSUE : ', e, ' >'
        print '   ################################################################'    
    
    t2 = time.time()-t1

    # Record and Display Results
    print '\n*************************************************************************'
    selected_PCP_manual = []
    dbf1 = ps.open(path+'Snapped/SERVICE_Snapped.dbf')
    NEW_Records_PCP_Man = []
    for v in manualPCP.getVars():
        if 'x' in v.VarName:
            pass
        elif 'W' in v.VarName:
            pass
        elif v.x > 0:
            var = '%s' % v.VarName
            selected_PCP_manual.append(var)
            for i in range(dbf1.n_records):
                if var in dbf1.read_record(i):
                    x = dbf1.read_record(i)
                    NEW_Records_PCP_Man.append(x)
                else:
                    pass
            print '    |                                            ', var,  '         '
    print '    | Selected Facility Locations -------------- ^^^^ ',  '         '
    print '    | Candidate Facilities [p] ----------------- ', len(selected_PCP_manual), '           '
    print '    | Objective Value (miles) ------------------ ', manualPCP.objVal, '     '
    print '    | Real Time to Optimize (sec.) ------------- ', t2
    print '    | Date/Time -------------------------------- ', dt.datetime.now()
    print '*************************************************************************'
    print ' -- The p-Center Problem Manual LP Creation-- '
    
# Call Function
try:
    Manual_LP_PCP()
    print '\nJames Gaboardi, 2015'
except Exception as e:
    exc_type, exc_obj, exc_tb = sys.exc_info()
    fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
    print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
    print '   ################################################################'
    print ' < ISSUE : ', e, ' >'
    print '   ################################################################'



    


















    






Parameter MIPFocus unchanged
   Value: 2  Min: 0  Max: 3  Default: 0
Optimize a model with 1701 rows, 1516 columns and 6111 nonzeros
Coefficient statistics:
  Matrix range    [1e-03, 3e+00]
  Objective range [1e+00, 1e+00]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 2e+00]
Presolve time: 0.02s
Presolved: 1701 rows, 1516 columns, 6015 nonzeros
Variable types: 1 continuous, 1515 integer (1515 binary)
Found heuristic solution: objective 1.7394879
Found heuristic solution: objective 1.7390156
Found heuristic solution: objective 1.7295779
Presolve removed 112 rows and 112 columns
Presolved: 1589 rows, 1404 columns, 5567 nonzeros


Root relaxation: objective 1.227915e+00, 1001 iterations, 0.03 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    1.22791    0  306    1.72958    1.22791  29.0%     -    0s
H    0     0                       1.6827877    1.22791  27.0%     -    0s
H    0     0                       1.6500480    1.22791  25.6%     -    0s
H    0     0                       1.6459287    1.22791  25.4%     -    0s
*    0     0               0       1.3697645    1.36976  0.00%     -    0s

Cutting planes:
  Gomory: 6
  Implied bound: 22
  Zero half: 99

Explored 0 nodes (5622 simplex iterations) in 0.19 seconds
Thread count was 4 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective 1.369764530500e+00, best bound 1.369764530500e+00, gap 0.0%

*************************************************************************
    |                                             y3          
    |                                             y15          
    | Selected Facility Locations -------------- ^^^^           
    | Candidate Facilities [p] -----------------  2            
    | Objective Value (miles) ------------------  1.3697645305      
    | Real Time to Optimize (sec.) -------------  0.1971180439
    | Date/Time --------------------------------  2015-11-14 16:23:10.773808
*************************************************************************
 -- The p-Center Problem Manual LP Creation-- 

James Gaboardi, 2015

2.2.4c      Instantiate selected manual p-Center shapefile

Following optimization, I create another new shapefile for the selected manual candidate service locations.





In [37]:



    


# define Manual PCP shapefile function
def Create_PCP_Man():
    try:
        # Instantiate Shapefile
        SHP_Center_Man = shp.Writer(shp.POINT)
        # Add Points
        for idy,idx,x,y in NEW_Records_PCP_Man:
            SHP_Center_Man.point(float(x), float(y))
        # Add Fields
        SHP_Center_Man.field('y_ID')
        SHP_Center_Man.field('x_ID')
        SHP_Center_Man.field('LAT')
        SHP_Center_Man.field('LON')
        # Add Records
        for idy,idx,x,y in NEW_Records_PCP_Man:
            SHP_Center_Man.record(idy,idx,x,y)
        # Save Shapefile    
        SHP_Center_Man.save(path+'Results/Selected_Locations_Pcenter_Man')
    except Exception as e:
        exc_type, exc_obj, exc_tb = sys.exc_info()
        fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
        print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
        print '   ################################################################'
        print ' < ISSUE : ', e, ' >'
        print '   ################################################################'

# Call Manual PCP shapefile function
try:
    Create_PCP_Man()
except Exception as e:
    exc_type, exc_obj, exc_tb = sys.exc_info()
    fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
    print exc_type, fname, 'Line Number -- ',exc_tb.tb_lineno
    print '   ################################################################'
    print ' < ISSUE : ', e, ' >'
    print '   ################################################################'

2.2.4d      p-Center Results compared

Comaprison of the automated and manual p-Center solution results:





In [38]:



    


print ' \n Automated p-Center'
print '    | Candidate Facilities [p] ----------------- ', len(selected_PCP), '           '
print '    | Selected Facility Locations -------------- ', selected_PCP
print '    | Objective Value (miles) ------------------ ', mPCP.objVal, '     '

print '\n Manual p-Center'
print '    | Candidate Facilities [p] ----------------- ', len(selected_PCP_manual), '           '
print '    | Selected Facility Locations -------------- ', selected_PCP_manual
print '    | Objective Value (miles) ------------------ ', manualPCP.objVal, '     '

print '\n\n    | Same Objective Value  -------------------- ', str(mPCP.ObjVal) == str(manualPCP.ObjVal)
print '    | Same Facilities  ------------------------- ', selected_PCP == selected_PCP_manual



    


















    






 
 Automated p-Center
    | Candidate Facilities [p] -----------------  2            
    | Selected Facility Locations --------------  ['y3', 'y15']
    | Objective Value (miles) ------------------  1.3697645305      

 Manual p-Center
    | Candidate Facilities [p] -----------------  2            
    | Selected Facility Locations --------------  ['y3', 'y15']
    | Objective Value (miles) ------------------  1.3697645305      


    | Same Objective Value  --------------------  True
    | Same Facilities  -------------------------  True

2.3      Selected locations

In the final step, I visualize the selected locations on the SCLP, PMP, and PCP by calling the shapefiles created after optimizing the problems.

The following plot shows the optimal locations as represented by snapped coordinates of the service sites.

Of particular note is the frequent discrepency of the optimal candidate locations between the PCP solutions without a seed. Although the objective value appears to be equal (or is nearly equal) after every trial, the PCP model built and solved with gurobipy generally runs either more quickly or with more simplex iterations demonstrating the benefits of building models within a wrapper.

Selected Set Cover, p-median, and p-center locations</u>





In [39]:



    


figsize(10,11)
# Draw Network Actual Roads and Nodes
nx.draw(g, ntw.node_coords, node_size=5, alpha=0.25, edge_color='r', width=2)

# Set Cover
SetCover = ps.open(path+'Results/Selected_Locations_SetCover.shp')
points_setcover = {}
for idx, coords in enumerate(SetCover):
    setcover_g.add_node(idx)
    points_setcover[idx] = coords
    setcover_g.node[idx] = coords
nx.draw(setcover_g, points_setcover, 
        node_size=1400, alpha=1, node_color='g')

# p-Median
P_Med = ps.open(path+'Results/Selected_Locations_Pmedian.shp')
points_median = {}
for idx, coords in enumerate(P_Med):
    median_g.add_node(idx)
    points_median[idx] = coords
    median_g.node[idx] = coords
nx.draw(median_g, points_median, 
        node_size=1000, alpha=1, node_color='r')

# p-Center
P_Cent = ps.open(path+'Results/Selected_Locations_Pcenter.shp')
points_center = {}
for idx, coords in enumerate(P_Cent):
    center_g.add_node(idx)
    points_center[idx] = coords
    center_g.node[idx] = coords
nx.draw(center_g, points_center, 
        node_size=700, alpha=1, node_color='b')

# p-Center Manual
P_Cent_Man = ps.open(path+'Results/Selected_Locations_Pcenter_Man.shp')
points_center_man = {}
for idx, coords in enumerate(P_Cent_Man):
    center_g_man.add_node(idx)
    points_center_man[idx] = coords
    center_g_man.node[idx] = coords
nx.draw(center_g_man, points_center_man, 
        node_size=300, alpha=1, node_color='y', node_shape='d' )

# Draw Graph of Random Client
nx.draw(GRAPH_client, points_client, 
        node_size=15, alpha=.5, node_color='k')

# Draw Graph of Snapped Service
nx.draw(g_service, ntw.pointpatterns['Rand_Points_SERVICE'].snapped_coordinates, 
        node_size=50, alpha=1, node_color='k')

# Legend (Ordered Dictionary)
LEGEND = OrderedDict()
LEGEND['Network Nodes']=g
LEGEND['Roads']=g
LEGEND['Optimal Set Cover (S=1.25)']=setcover_g
LEGEND['Optimal p-Median (p=2)']=median_g
LEGEND['Optimal p-Center (p=2)']=center_g
LEGEND['Optimal p-Center Manual(p=2)']=center_g_man
LEGEND['Client Nodes']=GRAPH_client
LEGEND['Snapped Service Nodes']=g_service
legend(LEGEND, 
       loc='lower right', 
       fancybox=True, 
       framealpha=0.5, 
       scatterpoints=1)

# 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.49, 0.0, 0.005, width=.0003, head_width=0.001, 
          head_length=0.0015, fc='k', ec='k',alpha=0.75,)
annotate('N', xy=(-84.2815, 30.498), fontstyle='italic', fontsize='xx-large',
            fontweight='heavy', alpha=0.75)



    


















    
Out[39]:








<matplotlib.text.Annotation at 0x115c66390>

3. Next Steps

  • Potential for refactoring the Dijkstra algorithm in PySAL.Network to increase efficiency

  • Improve cases for determining the simplified network to tackle to issue of network size

  • Integrate PySAL with an open source solver

GitHub

https://github.com/jGaboardi/NARSC__15





In [2]:



    


IPd.HTML('https://github.com/jGaboardi')



    


















    
Out[2]:









  
    
    
    
    
    
    
    jGaboardi (James Gaboardi) · GitHub
    
    
    
    
    
    
    

      
      
      
    
    
    
    
    

    
    
    

    
    






  
    




    
        
    

      
      

    

    


  
    Skip to content

    
    
    



      
      



    


    





    


        

          



  

    

        

      

  James Gaboardi
  jGaboardi








      

        
          0
          Followers
        
        
          5
          Starred
        
        
          13
          Following
        
      


    


    



      

          

              
    
  

          


        
      


      

          


  

    


      

        
Popular repositories

        
      

    


    

  


  

        

            Public contributions
        

    

      

          
            
          
      

      

        

            Summary of pull requests, issues opened, and commits.
          
            Learn how we count contributions.
        

        

          Less
          
    
            
  • 
            
  • 
            
  • 
            
  • 
            
  • 
          

          More
        

      


        

          Contributions in the last year
          298 total
          Nov 14, 2014 – Nov 14, 2015
        

        

          Longest streak
            10 days
            
              August 14 –
              August 23
              
        

        

          Current streak
            0 days
            Last contributed 
        

    

  


  

    


  
  




Contribution activity







    

      
        
        11 commits
      
    


    





  




      

    


  




        

        

    


      

  





    
    
    

    

      
      
      Something went wrong with that request. Please try again.

References

  • Church, R. L. 2002. Geographical information systems and location science. Computers and Operations Research 29:541–562.
  • Church, R. L., and A. T. Murray. 2009. Business Site Selections, Locational Analysis, and GIS. Hoboken, NJ, USA: John Wiley & Sons, Inc.
  • Current, J., and M. O. Kelly. 1992. Locating Emergency Warning Sirens. Decision Sciences 23:221–234.
  • Daskin, M. 1997. Network and Discrete Location: Models, Algorithms and Applications. Hoboken, NJ, USA: John Wiley & Sons, Inc.
  • Gass, S. I., and A. A. Assad. 2005. An Annotated Timeline of Operations Research: An Informal History. New York: Springer.
  • Hadley, G. 1963. Linear Programming. Reading, MA: London: Addison-Wesley Pub. Co.
  • Hakimi, S. L. 1964. Optimum Locations of Switching Centers and the Absolute Centers and Medians of a Graph. Operations Research 12 (3):450–459.
  • Hoffman, K. L., M. Padberg, and G. Rinaldi. 2011. The Traveling Salesman Problem. Encyclopedia of Operations Research and Management Science 3 (1):1573–1578.
  • Longley, P. A., M. F. Goodchild, D. J. Maguire, and D. W. Rhind. 2011. Geographic Information Systems & Science 3rd ed. Hoboken, NJ, USA: John Wiley & Sons, Inc.
  • Miller, H. J., and S.-L. Shaw. 2001. Geographic Information Systems for Transportation. New York: Oxford University Press.
  • Minieka, E. 1970. The m-Center Problem. SIAM Review 12:38–39.
  • Owen, S. H., and M. S. Daskin. 1998. Strategic facility location: A review. European Journal of Operational Research 111 (3):423–447.
  • Peet, R. 1998. Introduction: Geography, Philosophy, and Social Theory. In Modern Geographical Thought, 1–33. Wiley-Blackwell.
  • Ram, K. 2013. Git can facilitate greater reproducibility and increased transparency in science. Source code for biology and medicine 8 (1):7.
  • ReVelle, C. S. 2005. Location analysis: A synthesis and survey. European Journal of Operational Research 165:1–19.
  • ReVelle, C. S., and R. W. Swain. 1970. Central facilities location. Geographical Analysis 2 (1):30–42.
  • Teitz, M. B., and P. Bart. 1968. Heuristic Methods for Estimating the Generalized Vertex Median of a Weighted Graph. Operations Research 16 (5):955–961.

Content source: jGaboardi/NARSC__15

Similar notebooks:

notebook.community | gallery | about