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

SHRUG 2015

1. Optimization in GIS

       • Imports

       • What is Spatial Optimization?

       • Reproducibility

       • 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 [1]:



    


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

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

1.     Optimization & GIS

What is Spatial Optimization?

       • Geographic Information Systems & Science

       • Operations Research / Decision Science

       • Optimization & Mathematical Programming

Maximize

       9x + 7y

Subject To

       10x + 5y <= 50

       6x + 6y <= 36

       4.5x + 18y <= 81

Bounds

       x >= 0

       y >= 0

End

The equations above represent a basic linear program designed to maximize profit.

The Shop Problem

  • x is Product 1; y is Product 2

  • \$9 for each unit of ***x*** and \$7 for each unit of y

  • Constraint 1 is Machine 1 time (hrs.)

  • Constraint 2 is Machine 2 time (hrs.)
  • Constraint 3 is Machine 3 time (hrs.)

Solution:

  • objective value = \$50
  • x = 4 units
  • y = 2 units

Originally Published:

  • Greenwald, D. U.
  • 1957
  • Linear Programming: An Explanation of the Simplex Algorithm
  • New York, NY
  • The Ronald Press Company

       • Location Allocation / Facility Location

Linear Programming (LP) Integer Programming (IP)   Mixed-Integer Programming(MIP)  
   The Traveling Salesperson Problem    p-Median Problem p-Center Problem
The Transportation Problem   Set Cover Location Problem   Capacitated p-Center Problem

How do we solve these models and is the workflow 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 [see references].

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 [ ]:



    


print dir(ps)

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 [ ]:



    


print dir(ps.Network)

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 [ ]:



    


print dir(gbp)

"Traditional" Conceptual Model





In [ ]:



    


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

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.

Automated Model





In [ ]:



    


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

2.      Demonstration

Imports & Specs





In [ ]:



    


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 '********************************************************'

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 [see references]. 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.





In [ ]:



    


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

Loop Road





In [ ]:



    


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

Self-intersecting Road





In [ ]:



    


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

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 [ ]:



    


# 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 [ ]:



    


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

# Create a bounding box of the shapefile
shp_W.bbox

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 [ ]:



    


# 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 [ ]:



    


# 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 [ ]:



    


# 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 [ ]:



    


# CLIENT with {id: [lon, lat], } dictionary 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], } dictionary 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 [2.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 [ ]:



    


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)

Network Characteristics





In [ ]:



    


print '\nNodes in original network:   ', len(ntw.nodes)
print 'Edges in original network:   ', len(ntw.edges)
print '\nNodes in simplified network: ', len(g1.node)
print 'Edges in simplified network: ', len(ntw.graphedges)

Instantiate simplified network shapefile





In [ ]:



    


# List of coords by key
LC = []
for i,j in g1.edges():
    if i in g1.node and j in g1.node:
        x=[list(g1.node[i]), list(g1.node[j])]
        x = list(x)
        LC.append(x)
lc = [LC]  
SimpGraph = shp.Writer(shp.POLYLINE)
# Add Edges [k] in list of edges by end node coord
for k in lc:
    SimpGraph.poly(shapeType=shp.POLYLINE, parts=k)
# Add Fields
SimpGraph.field('Graph_ID')
counter = 0
for i in range(len(g1.node)):
    counter = counter + 1
    SimpGraph.record(counter)
# Save Shapefile    
SimpGraph.save(path+'Waverly/Simplified_Waverly.shp')

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 [ ]:



    


# 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 [ ]:



    


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

Instantiate Simplified Network





In [ ]:



    


SimNet = ps.Network(path+'Waverly/Simplified_Waverly.shp')

2.1.11a      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 for visualization.





In [ ]:



    


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.1.11b      Snap observations to SimNet

Also, I snap the client and service observations to the simplified network with the SimNet.snapobservations method. This allows for the creation of a distance matrix.





In [ ]:



    


t1 = time.time()
Snap_C = SimNet.snapobservations(path+'Simulated/RandomPoints_CLIENT.shp', 
                     'Rand_Points_CLIENT', attribute=True)
Snap_S = SimNet.snapobservations(path+'Simulated/RandomPoints_SERVICE.shp', 
                     'Rand_Points_SERVICE', attribute=True)
print round(time.time()-t1, 4), '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 [ ]:



    


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

2.1.13      Instantiate shapefile of service nodes snapped to ntw

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 [ ]:



    


# 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 SimNet.allneighbordistances() from the instantiated network object of the simplified network of Waverly Hills [Waverly_Simplified.shp]. The cost matrix in created on the to reduce the computational burden, although it is not correct in it's current form. The matrix is 100x15 and the values are converted from decimal degrees to miles.

A current limitation of this is the run-time to calculate accurate road segment length which is hindered by the number of nodes and edges of the network. When tested with much larger client and service point sets the shortest-path Dijkstra algorithm ran in similar time. This demonstrates the limitation of calculating a cost matrix with the actual road network. As a potential improvement, the Dijkstra may be refactored with a different data structure to see decreases in run-time. The run-time may decrease with another method of data storage.





In [ ]:



    


t1 = time.time()

# Define Client to Service Matrix Function
def c_s_matrix():
    global All_Dist_MILES
    All_Neigh_Dist = SimNet.allneighbordistances(sourcepattern=SimNet.pointpatterns['Rand_Points_CLIENT'],
                                             destpattern=SimNet.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

2.2      Mathematical Optimization

Set Parameters within Gurobi





In [ ]:



    


# Set Parameters
gbp.setParam('MIPFocus', 2)                 # Set MIP focus to 'Optimal' --> 2
gbp.setParam('MIPGapAbs', 0)                # Set Absolute MIP Gap  --> 0
gbp.setParam('GomoryPasses', 0)             # Set Number of Gomory Cuts  --> 0
gbp.setParam('ZeroHalfCuts', 0)             # Set Number of Zero Half Cuts  --> 0
gbp.setParam('ImpliedCuts', 0)              # Set Number of Implied Cuts  --> 0
gbp.setParam('BarConvTol', .000000001)      # Set Barrier Convergence Tolerence
gbp.setParam('FeasibilityTol', .000000001)  # Set Feasibility Tolerence
gbp.setParam('IntFeasTol', .000000001)      # Set Integer Feasibility Tolerence
gbp.setParam('OptimalityTol', .000000001)   # Set Optimality Tolerence
gbp.setParam('Method', 4)                   # Set Algorithm to 'concurrent': Dual Simplex and Barrier
gbp.setParam('DisplayInterval', 1)          # Set Display Interval to 1

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 [see references]. The objective of the SCLP is to minimize the number of service facilities to site while achieving complete (100%) coverage of client nodes.

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 [ ]:



    


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

    # Add Service Decision Variables (j)
    serv_var = []
    for dest in service_nodes:
        serv_var.append(mSCLP.addVar(vtype=gbp.GRB.BINARY,
                                    lb=0,
                                    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 '   ################################################################'

2.2.1b     Instantiate selected Set Cover location shapefile

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





In [ ]:



    


# 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 [see references]. 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 [ ]:



    


# 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, Add Variables, & Update Model
    # Instantiate Model
    mPMP = gbp.Model(' -- p-Median -- ')

    # 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,
                                                lb=0,
                                                ub=1,
                                                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,
                                          lb=0,
                                          ub=1,
                                          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 '   ################################################################'

2.2.2b     Instantiate Selected p-median shapefile

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





In [ ]:



    


# 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 [see references].

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 [ ]:



    


# 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, Add Variables, & Update Model
    # Instantiate Model
    mPCP = gbp.Model(' -- P-Center -- ')
    
    # 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, 
                                                lb=0,
                                                ub=1,
                                                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,
                                          lb=0,
                                          ub=1,
                                          name='y'+str(dest+1)))
    
    # Add the Maximum travel cost variable
    W = mPCP.addVar(vtype=gbp.GRB.CONTINUOUS,
                    lb=0.,
                    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 '   ################################################################'

2.2.3b      Instantiate selected p-Center shapefile

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





In [ ]:



    


# 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 [ ]:



    


# 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 [ ]:



    


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

2.2.4c      Instantiate selected manual p-Center shapefile

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





In [ ]:



    


# 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.5      Automated p-Center .lp File Solved

The .lp file created as an output when the intial p-center model was created and solved is then read into Gurobi through gurobipy and solved.





In [ ]:



    


# Define manual LP read PCP Function
def A_to_M_LP_PCP():
    
    global Cij
    global A_to_MPCP
    global dbf1
    global selected_PCP_A_to_M
    global NEW_Records_PCP_A_to_M
    
    t1 = time.time()

    # Instantiate Optimization model from .lp file
    A_to_MPCP = gbp.read(path+'LP_Files/WaverlyPCP.lp')
    
    # Solve
    try:
        A_to_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 '   ################################################################'    
    
    t2 = time.time()-t1

    # Record and Display Results
    print '\n*************************************************************************'
    selected_PCP_A_to_M = []
    dbf1 = ps.open(path+'Snapped/SERVICE_Snapped.dbf')
    NEW_Records_PCP_A_to_M = []
    for v in A_to_MPCP.getVars():
        if 'x' in v.VarName:
            pass
        elif 'W' in v.VarName:
            pass
        elif v.x > 0:
            var = '%s' % v.VarName
            selected_PCP_A_to_M.append(var)
            for i in range(dbf1.n_records):
                if var in dbf1.read_record(i):
                    x = dbf1.read_record(i)
                    NEW_Records_PCP_A_to_M.append(x)
                else:
                    pass
            print '    |                                            ', var,  '         '
    print '    | Selected Facility Locations -------------- ^^^^ ',  '         '
    print '    | Candidate Facilities [p] ----------------- ', len(selected_PCP_A_to_M), '           '
    print '    | Objective Value (miles) ------------------ ', A_to_MPCP.objVal, '     '
    print '    | Real Time to Optimize (sec.) ------------- ', t2
    print '    | Date/Time -------------------------------- ', dt.datetime.now()
    print '*************************************************************************'
    print ' -- The p-Center Problem Solved by rereading the Auto-creation LP -- '
    
# Call thFunction
try:
    A_to_M_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 '   ################################################################'

2.2.4d      p-Center Results compared

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





In [ ]:



    


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) --------------- ', str(mPCP.ObjVal) == str(manualPCP.ObjVal)
print '    | Same Objective Value (float)  ------------ ', mPCP.ObjVal.real == manualPCP.ObjVal.real
print '    | Same Facilities  ------------------------- ', selected_PCP == selected_PCP_manual

print '\nIdentical Parameters? ', str(mPCP.Params) == str(manualPCP.Params)

print '\n1: ', ("%.50f" % mPCP.ObjVal)
print '2: ', ("%.50f" % manualPCP.ObjVal)
print '3: ', ("%.50f" % A_to_MPCP.ObjVal)

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.

Of particular note is the frequent discrepancy of the optimal candidate locations between the PCP solutions without a seed even when run with the exact same parameters. Although the objective values are generally equal up to 10 significant after every trial, the PCP model built and solved with gurobipy generally runs either more quickly or with more simplex iterations. The cause of the discrepancy is yet unknown.

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





In [ ]:



    


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)

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/SHRUG__2015





In [ ]:



    


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

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/SHRUG__2015

Similar notebooks:

notebook.community | gallery | about