How to use Linear Programming with ArcGIS Pro to optimise Covid vaccinating
Rhys Donoghue
MBA | Technology Leader | Full Stack Developer | Esri-certified system designer and developer
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)
Creative and Strategic Leader. Woman of Steel.
3 年That’s pretty amazing!! Hope all is going well. Long time no see. :)
Kaiwhakahaere Kaupapa me te Hōtaka | Project and Programme Manager | MSocSci | PMP? | NZPSM
3 年Top marks Rhys.
Geospatial Data Engineer at INDRA | Data | GIS
3 年This is awesome Rhys, thanks for sharing!