How to use Linear Programming with ArcGIS Pro to optimise Covid vaccinating

ArcGIS Pro has excellent support for Python, and Python has excellent support for solving optimisation models. This demo shows how a linear programming model was created in Python using PuLP to solve the model. Then the script was modified to render the results in ArcGIS Pro. The Python script can easily be added to an ArcGIS Pro toolbox so that end users do not need to see any code. For this example, I have used a very basic example to show how the model can be solved with PuLP and added to ArcGIS Pro.

To make it easier to understand what is being optimised, I have used hospitals and suburbs, with the objective function for the linear program being to minimise the total travelled distance from each suburb to the nearest hospital for Covid vaccinations. For a real model for Covid vaccinations, we would not use suburb centroids and Euclidean (straight line) distances from the suburbs to the hospitals. They are only being used for this model to show how to use Python, PuLP, and ArcGIS Pro together. Likewise, for a real model, we would use a more suitable solving engine like Gurobi. All of the code for the demo is shown in the YouTube video.

Send me a message on LinkedIn if you need help running the code, or if you need help building a real world model, or if you want to me to post the QGIS (open source) alternative if you do not have access to ArcGIS Pro.

The YouTube video is here: https://www.youtube.com/watch?v=qbeN59otEr8

The Python code in the video is here:

from pulp import *
import arcpy
import sys


wksp = r"C:\gisdata\demos\lp\lp.gdb"
arcpy.env.workspace = wksp
arcpy.env.overwriteOutput = True


# create the results feature class
results_fc = "results"


out_path = wksp
out_name = results_fc
geometry_type = "POLYLINE"
template = ""
has_m = "DISABLED"
has_z = "DISABLED"
spatial_reference = arcpy.SpatialReference(4326)


if arcpy.Exists(results_fc):
    arcpy.Delete_management(results_fc)


arcpy.CreateFeatureclass_management(
    out_path, out_name, geometry_type, template, has_m, has_z, spatial_reference)


# create the fields for the results feature class
results_fields = ['HOSPITAL', 'TO', 'DISTANCE']
arcpy.AddField_management(in_table=results_fc, field_name=results_fields[0],
                          field_type="TEXT", field_length="50")
arcpy.AddField_management(in_table=results_fc, field_name=results_fields[1],
                          field_type="TEXT", field_length="50")
arcpy.AddField_management(in_table=results_fc, field_name=results_fields[2],
                          field_type="LONG")


results_fields.append('SHAPE@')


cost = {}
hospitals = []
suburbs = []


with arcpy.da.SearchCursor("hospitals", ['Name', 'Capacity', 'SHAPE@XY']) as cursor1:
    for row1 in cursor1:
        hospitals.append(
            {
                'name': row1[0],
                'capacity': row1[1]
            }
        )


        cost[row1[0]] = {}


        with arcpy.da.SearchCursor("suburbs_centroids", ["LOCALITY", 'SHAPE@XY']) as cursor2:
            for row2 in cursor2:
                if(not row2[0] in suburbs):
                    suburbs.append(row2[0])


                origin_x = row1[2][0]
                origin_y = row1[2][1]
                destination_x = row2[1][0]
                destination_y = row2[1][1]


                array = arcpy.Array(
                    [arcpy.Point(origin_x, origin_y), arcpy.Point(destination_x, destination_y)])


                spatial_reference = arcpy.SpatialReference(4326)
                line = arcpy.Polyline(array, spatial_reference)


                distance = int(line.getLength('GEODESIC', 'METERS'))


                cost[row1[0]][row2[0]] = {}
                cost[row1[0]][row2[0]]['distance'] = distance
                cost[row1[0]][row2[0]]['line'] = line


hospital_names = [h['name'] for h in hospitals]


# define the decision variables
x = LpVariable.dicts("x", [(hospital, suburb)
                           for hospital in hospital_names for suburb in suburbs], cat=LpBinary)


prob = LpProblem("ilp", LpMinimize)


# create the objective function
obj_function = ''


for hospital, suburb in x:
    obj_function += int(cost[hospital][suburb]
                        ['distance']) * x[(hospital, suburb)]


prob += obj_function, "Total cost travelled"


# each suburb is assigned to one hospital
for suburb in suburbs:
    prob += pulp.lpSum(x[(hospital, suburb)]
                       for hospital in hospital_names) == 1


# each hospital can only have a set number of points
for hospital in hospitals:
    prob += pulp.lpSum(x[(hospital['name'], suburb)]
                       for suburb in suburbs) == int(hospital['capacity'])


# solve the model
prob.solve()


print('Status: %s' % (LpStatus[prob.status]))
print("Total travel cost = %s" % (value(prob.objective)))


results = {}


for v in prob.variables():
    if (v.varValue == 1):
        r = v.name.replace("x_('", '').replace(
            "',_'", ',').replace("')", '').split(',')
        hospital, suburb = r[0], r[1].replace('_', ' ')


        if not hospital in results:
            results[suburb] = {}
            results[suburb]['hospital'] = hospital
            results[suburb]['distance'] = cost[hospital][suburb]['distance']
            results[suburb]['line'] = cost[hospital][suburb]['line']


insert_cursor = arcpy.da.InsertCursor(results_fc, results_fields)


for k, v in results.items():
    origin, to, distance, line = v['hospital'], k, v['distance'], v['line']
    insert_cursor.insertRow((origin, to, distance, line))


del insert_cursor


with arcpy.da.UpdateCursor('suburbs', ['LOCALITY', 'HOSPITAL']) as cursor:
    for row in cursor:
        row[1] = results[row[0]]['hospital']
        cursor.updateRow(row)