Spatial Data Analysis（六）：空间优化问题
Spatial Data Analysis（六）：空间优化问题
使用pulp库解决空间优化问题：
pulp是一个用于优化问题的Python库。它包含了多种优化算法和工具，可以用于线性规划、混合整数线性规划、非线性规划等问题。Pulp提供了一个简单的方式来定义优化问题，包括变量、约束和目标函数，并且可以使用多种求解器进行求解。Pulp也提供了可视化工具来展示优化问题的结果。Pulp是一个开源项目，可以在GitHub上获取它的源代码。
空间优化（一）：p-中值问题
这个问题需要p设施的位置，同时最小化服务所有需求的总加权距离。
每个节点都有一个关联的权重，表示该节点的需求量。
目标函数： 最小化所有设施和需求节点的需求加权总和。
决策变量： 将设施放置在何处以及哪个设施位置为哪些需求节点提供服务
限制：
- 每个节点由 1 个设施提供服务
- 仅当某个位置存在设施时，节点才可以由该设施提供服务。
- 我们必须放置p设施
- 每个节点要么是一个设施，要么不是。
pip install -q pulp
from pulp import *
import numpy as np
import geopandas as gp
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
#read a sample shapefile
georgia_shp = gp.read_file("https://raw.githubusercontent.com/Ziqi-Li/GEO4162C/main/data/georgia/G_utm.shp")
georgia_shp.shape
(172, 18)
创建一个需求和一个设施变量，表示每个需求和设施的索引。
需求节点：所有县
facility：设施将建在一些选定的县之上
#create a demand and a facilities variable, indicating the indices of each demand and facility.
#demand node: all counties
#facility: Facilities will be built on top of some chosen countiesdemand = np.arange(0,172,1)
facilities = np.arange(0,172,1)
计算距离矩阵d_ij(n×n)
#Calculate a distance matrix d_ij (n by n)
coords = list(zip(georgia_shp.centroid.x,georgia_shp.centroid.y))d = cdist(coords,coords)
每个县(hi)的需求是总人口
#the demand for each county (h_i) is the total populatoion
h = georgia_shp.TotPop90.values
声明设施变量；生成的变量名称为：X_1,X_2,…
# declare facilities variables;the resulting variable names are: X_1,X_2,...
X = LpVariable.dicts('X_%s',(facilities),cat='Binary')# declare demand-facility pair variables; the resulting variable names are Y_0_1, Y_0_2,...
Y = LpVariable.dicts('Y_%s_%s', (demand,facilities),cat='Binary')
要放置的设施数量
#Number of facilities to place
p = 3 #change this and re-run the code.#Create a new problem
prob = LpProblem('P_Median', LpMinimize)
目标函数：最小化所有设施和需求节点的加权需求距离总和
(h_i: i 处的需求；d_ij: i 和 j 之间的距离)
“for”循环用于迭代序列
# Objective function: Minimizing weighted demand-distance summed over all facilities and demand nodes
# (h_i: demand at i; d_ij: distance between i and j)
# A "for" loop is used for iterating over a sequenceprob += sum(sum(h[i] * d[i][j] * Y[i][j] for j in facilities) for i in demand)
这个约束表明我们必须精确放置 p 个设施
# This constraint indicates we must place exactly p facilitiesprob += sum([X[j] for j in facilities]) == p
这一约束意味着需求节点 i 只能由一个设施提供服务
# This constraint implies that a demand node i can only be serviced by one facilityfor i in demand:prob += sum(Y[i][j] for j in facilities) == 1
这个约束意味着需求节点 i
仅当 j 处有设施时才能由 j 处的设施提供服务
它隐式地消除了 X[j] = 0 但 Y[i][j] = 1 时的情况
（节点 i 由 j 提供服务，但 j 处没有设施）
# This constraint implies that that demand node i
# can be serviced by a facility at j only if there is a facility at j
# It implicitly removes situation when X[j] = 0 but Y[i][j] = 1
# (node i is served by j but there is no facility at j)for i in demand:for j in facilities:prob += Y[i][j] <= X[j]
%%time# Solve the above problem
prob.solve()print("Status:", LpStatus[prob.status])
Status: Optimal
CPU times: user 1.35 s, sys: 64 ms, total: 1.42 s
Wall time: 11.5 s
# The minimized total demand-distance. The unit is person * meter (total distance travelled)
print("Objective: ",value(prob.objective))
Objective: 469538765110.4489
# Print the facility node.
rslt=[]
for v in prob.variables():subV = v.name.split('_')if subV[0] == "X" and v.varValue == 1:rslt.append(int(subV[1]))print('Facility Node: ', subV[1])
Facility Node: 126
Facility Node: 30
Facility Node: 82
# Get the geomerty of the facility nodes.
fac_loc = georgia_shp.iloc[rslt,:]
fac_loc
| AREA | PERIMETER | G_UTM_ | G_UTM_ID | AREANAME | Latitude | Longitud | TotPop90 | PctRural | PctBach | PctEld | PctFB | PctPov | PctBlack | X | Y | AreaKey | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 126 | 7.315030e+08 | 117190.0 | 130 | 128 | GA, Crisp County | 31.92540 | -83.77159 | 20011 | 48.4 | 10.0 | 12.47 | 0.30 | 29.0 | 40.66 | 805648.4 | 3537103 | 13081 | POLYGON ((787012.250 3547615.750, 820243.312 3... |
| 30 | 1.385270e+09 | 274218.0 | 32 | 31 | GA, Fulton County | 33.78940 | -84.46716 | 648951 | 4.2 | 31.6 | 9.63 | 4.13 | 18.4 | 49.92 | 733728.4 | 3733248 | 13121 | POLYGON ((752606.688 3785970.500, 752835.062 3... |
| 82 | 9.179670e+08 | 121744.0 | 84 | 84 | GA, Jenkins County | 32.78866 | -81.96042 | 8247 | 53.8 | 7.7 | 13.10 | 0.21 | 27.8 | 41.51 | 970465.7 | 3640263 | 13165 | POLYGON ((989566.750 3653155.750, 981378.062 3... |
#Plot the faclities (stars) on top of the demand map.
fig, ax = plt.subplots(figsize=(5,5))georgia_shp.centroid.plot(ax=ax,markersize=georgia_shp.TotPop90/1000)#markersize is proportional to the population
fac_loc.centroid.plot(ax=ax,color="red",markersize=300,marker="*")
空间优化（二）：集合覆盖问题
在此模型中，设施可以为距设施给定覆盖距离 Dc 内的所有需求节点提供服务。 问题在于放置最少数量的设施，以确保所有需求节点都能得到服务。 我们假设设施没有容量限制。
pip install -q pulp
from pulp import *
import numpy as np
import geopandas as gp
from scipy.spatial.distance import cdistimport matplotlib.pyplot as plt
#read a sample shapefile
georgia_shp = gp.read_file("https://raw.githubusercontent.com/Ziqi-Li/GEO4162C/main/data/georgia/G_utm.shp")
georgia_shp.shape
(172, 18)
创建一个需求和一个设施变量，表示每个需求和设施的索引。
需求节点：所有县
facility：我们可以在一些县建造设施
#create a demand and a facilities variable, indicating the indices of each demand and facility.
#demand node: all counties
#facility: we could build facilities in some countiesdemand = np.arange(0,172,1)
facilities = np.arange(0,172,1)
计算距离矩阵d_ij(n×n)
#Calculate a distance matrix d_ij (n by n)
coords = list(zip(georgia_shp.centroid.x,georgia_shp.centroid.y))
d = cdist(coords,coords)
阈值覆盖距离
# Threshold coverage distance
Dc = 100000 #100km coverage, change this and re run the code.
创建一个变量，指示节点 i 是否可以被设施 j 覆盖。
#Creata a variable (alpha in the lecture slide pg.28), indicating whether a node i can be covered by facility j.
a = np.zeros(d.shape)
a[d <= Dc] = 1
a[d > Dc] = 0
声明设施变量 Xj
# declare facilities variables Xj
X = LpVariable.dicts('X_%s',(facilities),cat='Binary')
创建一个最小化问题
#Create an minimization problem
prob = LpProblem('Set_Covering', LpMinimize)
目标函数：我们要最小化放置设施的数量
# Objective function: we want to minimize the number of placed facilities
prob += sum([X[j] for j in facilities])
该约束意味着每个需求节点 i 需要至少由设施服务
# This constraint implies every demand node i needs to be served by at least facility
for i in demand:prob += sum(a[i][j]*X[j] for j in facilities) >= 1%%time
# Solve the above problem
prob.solve()print("Status:", LpStatus[prob.status])
Status: Optimal
CPU times: user 22.5 ms, sys: 1.05 ms, total: 23.6 ms
Wall time: 66.4 ms
# The minimal number of facilities with the defiened coverage.
print("Objective: ",value(prob.objective))
Objective: 8.0
# Print the facility nodes.
rslt = []
for v in prob.variables():subV = v.name.split('_')if subV[0] == "X" and v.varValue == 1:rslt.append(int(subV[1]))print('Facility Node: ', subV[1])
Facility Node: 102
Facility Node: 120
Facility Node: 145
Facility Node: 150
Facility Node: 30
Facility Node: 38
Facility Node: 9
Facility Node: 97
# Get the geomerty of the facility nodes.
fac_loc = georgia_shp.iloc[rslt,:]
#Plot the faclities (stars) on top of the demand map.
fig, ax = plt.subplots(figsize=(5,5))georgia_shp.centroid.plot(ax=ax)
fac_loc.centroid.plot(ax=ax,color="red",markersize=300,marker="*")
<Axes: >
