위치 할당 최적화
이 예에서는 오픈 소스 Python을 사용하여 대규모 할당 및 위치 할당 문제에 접근하는 방법을 보여줍니다. 워크플로우는 익명화된 데이터를 사용하며 모델 입력을 준비하고, 출발지-목적지 매트릭스를 구축하고, 최적화를 해결하고, 출력을 다시 GeoJSON으로 내보내 GIS 소프트웨어에서 검토할 수 있는 방법을 보여줍니다.
계획 질문은 간단하지만 유용합니다. 각 동네나 교외의 인구, 각 시설의 수용 능력, 각 시설이 제공할 수 있는 최대 거리를 고려할 때 전체 인구를 서비스할 수 있는 시설이 충분한지, 그렇지 않은 경우 현재 제약으로 인해 할당되지 않은 지역은 어디입니까?
모델이 하는 일
이는 할당 및 위치 할당 모델입니다. 각 수요 위치는 육각형으로 표시되고, 각 시설은 사용 가능한 공급을 나타내며, 모델은 서비스 거리 및 용량 규칙을 존중하면서 모든 육각형을 단일 최적 시설에 할당합니다. 적용 범위, 서비스 영역 및 액세스 격차를 이해해야 하는 의료 시설, 학교, 은행, 소매점 또는 기타 서비스 네트워크에 동일한 모델링 패턴을 사용할 수 있습니다.
출력은 단순한 할당 테이블이 아닙니다. 또한 각 시설에서 실현 가능한 위치, 할당되지 않은 폴백으로 푸시되는 수요 영역, 용량, 거리 또는 목표가 조정될 때 결과가 어떻게 변경되는지를 보여주는 서비스 영역 및 무역 영역 레이어로 전환될 수 있습니다.
목적 함수
이 예의 목표는 모든 수요 육각형에서 할당된 시설까지의 총 거리를 최소화하는 것입니다. 이는 현재 가정 하에서 가장 효율적인 할당을 제공합니다. 비즈니스 질문이 변경되면 동일한 프레임워크를 다른 목적 함수로 다시 실행할 수 있습니다. 예를 들어, 거리 대신 운영 비용을 최소화하거나, 직선 거리 대신 이동 시간을 최소화하거나, 수요를 재조정하여 고비용 시설이 운영 프로필을 정당화할 수 있는 충분한 양을 확보할 수 있습니다.
키 제약조건
- 각 시설에는 최대 서비스 거리가 있으므로 해당 범위를 벗어나는 수요는 해당 시설에 할당할 수 없습니다.
- 각 시설에는 최소 및 최대 인구 제한이 있으므로 할당된 총 수요는 해당 용량 제한 내에서 유지되어야 합니다.
- 각 수요 육각형은 단일 시설에 할당되어야 합니다.
- 육각형을 실제 시설에 할당할 수 없는 경우, 모형
unassigned시설에 할당하여 모델의 실행 가능성을 유지하고 Coverage Gap을 확인할 수 있습니다. - 예시에서는 빠른 프로토타이핑을 위해 육각형 중심과 시설 사이의 유클리드 거리를 사용하지만, 필요한 경우 동일한 워크플로우를 실제 주행 가능한 경로로 전환할 수 있습니다.
작업흐름
워크플로에는 세 가지 주요 단계가 있습니다.
- 관심 영역을 육각형 테셀레이션으로 변환하고 각 육각형에 대한 인구를 도출합니다.
- 각 수요 육각형과 각 시설 간의 출발지-목적지 매트릭스를 생성합니다.
- 최적화 모델을 해결하고 결과를 GeoJSON으로 내보내 QGIS, ArcGIS Pro 또는 기타 매핑 소프트웨어에서 검토할 수 있습니다.
구현 스택에는 Shapely, pyproj, sqlite3, pyomo, CBC 솔버 및 GeoJSON 출력을 갖춘 오픈 소스 Python이 포함되어 있습니다. 데이터 준비 패턴은 의도적으로 모듈식으로 되어 있어 수요 계층 교체, 시설 제약 조건 변경, 목표 교체 또는 추가 비즈니스 규칙으로 모델 확장을 쉽게 할 수 있습니다.
구현 참고사항- PuLP는 원래 테스트되었지만 훨씬 더 큰 모델을 더 안정적으로 처리하기 때문에 대신 pyomo가 선택되었습니다.
- 모델은 오픈 소스 CBC 솔버로 해결되었으며, 이 접근 방식은 해당 설정을 통해 1시간 이내에 5천만 개 이상의 의사 결정 변수로 확장되었습니다.
- 더 큰 규모의 인스턴스의 경우 라이선스가 허용하는 경우 구로비를 고려할 수 있습니다.
- GeoJSON에 대규모 출력을 작성하는 것은 모델 자체를 해결하는 것보다 시간이 더 오래 걸릴 수 있으므로 대규모 생산 실행의 경우 데이터베이스에 직접 작성하는 것이 더 효율적일 수 있습니다.
- 이와 같은 모델을 구축하는 실용적인 방법은 제약 조건을 테스트하면서 큰 육각형과 빠른 유클리드 거리로 시작한 다음, 모델 동작이 검증된 후 더 미세한 테셀레이션과 더 현실적인 경로 비용으로 전환하는 것입니다.
- 추가 제약 조건은 점진적으로 추가할 수 있지만 각각의 추가 비즈니스 규칙으로 인해 모델을 실행 불가능하게 만들 위험이 증가하므로 신중하게 도입해야 합니다.
예제 소스 코드
아래 코드는 이 페이지에서 직접 엔드투엔드 워크플로를 보여줍니다.
1. generate_hexagons.py
import json
import math
import os
import pyproj
from shapely.geometry import shape
# for converting the coordinates to and from geographic and projected coordinates
TRAN_4326_TO_3857 = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3857")
TRAN_3857_TO_4326 = pyproj.Transformer.from_crs("EPSG:3857", "EPSG:4326")
# the area of interest used for generating the hexagons
input_geojson_file = "input/area_of_interest.geojson"
# load the area of interest into a JSON object
with open(input_geojson_file) as json_file:
geojson = json.load(json_file)
# the area of interest coordinates (note this is for a single-part / contiguous polygon)
geographic_coordinates = geojson["features"][0]["geometry"]["coordinates"]
# create an area of interest polygon using shapely
aoi = shape({"type": "Polygon", "coordinates": geographic_coordinates})
# get the geographic bounding box coordinates for the area of interest
(lng1, lat1, lng2, lat2) = aoi.bounds
# get the projected bounding box coordinates for the area of interest
[W, S] = TRAN_4326_TO_3857.transform(lat1, lng1)
[E, N] = TRAN_4326_TO_3857.transform(lat2, lng2)
# the area of interest height
aoi_height = N - S
# the area of interest width
aoi_width = E - W
# the length of the side of the hexagon
l = 200
# the length of the apothem of the hexagon
apo = l * math.sqrt(3) / 2
# distance from the mid-point of the hexagon side to the opposite side
d = 2 * apo
# the number of rows of hexagons
row_count = math.ceil(aoi_height / l / 1.5)
# add a row of hexagons if the hexagon tessallation does not fully cover the area of interest
if(row_count % 2 != 0 and row_count * l * 1.5 - l / 2 < aoi_height):
row_count += 1
# the number of columns of hexagons
column_count = math.ceil(aoi_width / d) + 1
# the total height and width of the hexagons
total_height_of_hexagons = row_count * l * 1.5 if row_count % 2 == 0 else 1.5 * (row_count - 1) * l + l
total_width_of_hexagons = (column_count - 1) * d
# offsets to center the hexagon tessellation over the bounding box for the area of interest
x_offset = (total_width_of_hexagons - aoi_width) / 2
y_offset = (row_count * l * 3 / 2 - l / 2 - aoi_height - l) / 2
# create an empty feature collection for the hexagons
feature_collection = { "type": "FeatureCollection", "features": [] }
oid = 1
hexagon_count = 0
for i in range(0, column_count):
for j in range(0, row_count):
if(j % 2 == 0 or i < column_count - 1):
x = W - x_offset + d * i if j % 2 == 0 else W - x_offset + apo + d * i
y = S - y_offset + l * 1.5 * j
coords = []
for [lat, lng] in [
TRAN_3857_TO_4326.transform(x, y + l),
TRAN_3857_TO_4326.transform(x + apo, y + l / 2),
TRAN_3857_TO_4326.transform(x + apo, y - l / 2),
TRAN_3857_TO_4326.transform(x, y - l),
TRAN_3857_TO_4326.transform(x - apo, y - l / 2),
TRAN_3857_TO_4326.transform(x - apo, y + l / 2),
TRAN_3857_TO_4326.transform(x, y + l)
]:
coords.append([lng, lat])
hexagon = shape({"type": "Polygon", "coordinates": [coords]})
# check if the hexagon is within the area of interest
if aoi.intersects(hexagon):
hexagon_count += 1
if(hexagon_count % 1000 == 0):
print('Generated {} hexagons'.format(hexagon_count))
population = 0
hexagon_names = []
# open the geojson file with the population data
with open("input/population_areas.geojson") as json_file:
geojson = json.load(json_file)
for feature in geojson["features"]:
polygon = shape(
{
"type": "Polygon",
"coordinates": feature["geometry"]["coordinates"]
}
)
# check if hexagon is within the polygon and derive the population for that intersected part of the hexagon
if hexagon.intersects(polygon):
if not feature["properties"]["Name"] in hexagon_names:
hexagon_names.append(feature["properties"]["Name"])
population += (
hexagon.intersection(polygon).area
/ polygon.area
* feature["properties"]["Population"]
)
hexagon_names.sort()
f = {
"type": "Feature",
"properties": {
"id": oid,
"name": ', '.join(hexagon_names),
"population": population
},
"geometry": {
"type": "Polygon",
"coordinates": [coords]
}
}
# add the hexagon to the feature collection
feature_collection['features'].append(f)
oid += 1
print('Generated {} hexagons'.format(hexagon_count))
# output the feature collection to a geojson file
with open("output/hexagons.geojson", "w") as output_file:
output_file.write(json.dumps(feature_collection))
# Play a sound when the script finishes (macOS)
for i in range(1, 2):
os.system('afplay /System/Library/Sounds/Glass.aiff')
# Play a sound when the script finishes (Windows OS)
# import time
# import winsound
# frequency = 1000
# duration = 300
# for i in range(1, 10):
# winsound.Beep(frequency, duration)
# time.sleep(0.1)2. generate_origin_destination_matrix.py
import json
import math
import os
import pyproj
import sqlite3
def getDistance(x1,y1,x2,y2):
distance = math.sqrt((x2-x1)**2+(y2-y1)**2)
return int(distance)
def getHexagonCentroid(hexagon):
coordinates = hexagon['geometry']['coordinates'][0]
# remove the last pair of coordinates in the hexagon
coordinates.pop()
lat = sum(coords[1] for coords in coordinates) / 6
lng = sum(coords[0] for coords in coordinates) / 6
return lat, lng
# for converting the coordinates to and from geographic and projected coordinates
TRAN_4326_TO_3857 = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3857")
TRAN_3857_TO_4326 = pyproj.Transformer.from_crs("EPSG:3857", "EPSG:4326")
# create a sqlite database for the results
db = 'output/results.sqlite'
# delete the database if it already exists
if(os.path.exists(db)):
os.remove(db)
# create a connection to the sqlite database
conn = sqlite3.connect(db)
# create cursors for the database connection
c1 = conn.cursor()
c2 = conn.cursor()
# create the facilities table
c1.execute('''
CREATE TABLE facilities (
facility_id INT,
facility_x REAL,
facility_y REAL,
trade_area_distance_constraint REAL,
min_population_constraint INT,
max_population_constraint INT
);
''')
c1.execute('''
CREATE TABLE od_matrix (
facility_id INT,
hexagon_id INT,
facility_x REAL,
facility_y REAL,
hexagon_x REAL,
hexagon_y REAL,
distance INT,
optimal INT
);
''')
# the geojson for the facilities
input_geojson_file = "input/facilities.geojson"
# load the area of interest into a JSON object
with open(input_geojson_file) as json_file:
geojson = json.load(json_file)
facilities = geojson['features']
for facility in facilities:
facility_id = facility['properties']['OID']
trade_area_distance_constraint = facility['properties']['trade_area_dist_constraint']
min_population_constraint = facility['properties']['min_population_constraint']
max_population_constraint = facility['properties']['max_population_constraint']
[lng, lat] = facility['geometry']['coordinates']
[x, y] = TRAN_4326_TO_3857.transform(lat, lng)
sql = '''
INSERT INTO facilities (
facility_id,
facility_x,
facility_y,
trade_area_distance_constraint,
min_population_constraint,
max_population_constraint
)
VALUES ({},{},{},{},{},{});
'''.format(facility_id, x, y, trade_area_distance_constraint, min_population_constraint, max_population_constraint)
c1.execute(sql)
# create an empty feature collection for the hexagons
feature_collection = { "type": "FeatureCollection", "features": [] }
# the geojson for the hexagons
input_geojson_file = "output/hexagons.geojson"
# load the area of interest into a JSON object
with open(input_geojson_file) as json_file:
geojson = json.load(json_file)
hexagons = geojson['features']
for i, hexagon in enumerate(hexagons):
hexagon_id = hexagon['properties']['id']
lat, lng = getHexagonCentroid(hexagon)
[hexagon_x, hexagon_y] = TRAN_4326_TO_3857.transform(lat, lng)
rows = c1.execute('''
SELECT facility_id, facility_x, facility_y, trade_area_distance_constraint
FROM facilities;
''').fetchall()
for row in rows:
facility_id, facility_x, facility_y, trade_area_distance_constraint = row
distance = getDistance(hexagon_x, hexagon_y, facility_x, facility_y)
if(distance > int(trade_area_distance_constraint)):
distance = 100000
sql = '''
INSERT INTO od_matrix(facility_id, facility_x, facility_y, hexagon_id, hexagon_x, hexagon_y, distance)
VALUES ({},{},{},{},{},{},{});
'''.format(facility_id, facility_x, facility_y, hexagon_id, hexagon_x, hexagon_y, distance)
c2.execute(sql)
(hexagon_lat, hexagon_lng) = TRAN_3857_TO_4326.transform(hexagon_x, hexagon_y)
(facility_lat, facility_lng) = TRAN_3857_TO_4326.transform(facility_x, facility_y)
coords = [[facility_lng, facility_lat],[hexagon_lng, hexagon_lat]]
if(distance <= trade_area_distance_constraint):
f = {
"type": "Feature",
"properties": {},
"geometry": {
"type": "LineString",
"coordinates": coords
}
}
# add the hexagon to the feature collection
feature_collection['features'].append(f)
if((i+1) % 1000 == 0):
print('Processed {} hexagons'.format(i+1))
conn.commit()
conn.close()
# output the feasible facility-hexagon pairs to geojson
with open("output/routes.geojson", "w") as output_file:
output_file.write(json.dumps(feature_collection))
for i in range(1,2):
os.system('afplay /System/Library/Sounds/Glass.aiff')3.solv_the_model.py
from pyomo.environ import *
from pyomo.opt import SolverFactory
import datetime
import json
import os
import pyproj
import sqlite3
# for converting the coordinates to and from geographic and projected coordinates
TRAN_4326_TO_3857 = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3857")
TRAN_3857_TO_4326 = pyproj.Transformer.from_crs("EPSG:3857", "EPSG:4326")
# the input database created from the previous script
db = 'output/results.sqlite'
# create a database connection
conn = sqlite3.connect(db)
# create a cursor for the database connection
c = conn.cursor()
# the demand, supply, and cost matrices
Demand = {}
Supply = {}
Cost = {}
'''
Supply['S0'] is for infeasible results, i.e. hexagons that do not
have any facilities when the nearest facility is too far away,
or when the population constraint for the facilities means the
hexagon cannot be assigned to that facility
'''
# the population capacity constraint of the "unassigned" facility
Supply['S0'] = {}
Supply['S0']['min_pop'] = 0
Supply['S0']['max_pop'] = 1E10
sql = '''
SELECT DISTINCT hexagon_id
FROM od_matrix
ORDER BY 1;
'''
# the assignment constraint, i.e. each hexagon can only be assigned to one facility
for row in c.execute(sql):
hexagon_id = row[0]
d = 'D{}'.format(hexagon_id)
Demand[d] = 1
# the infeasible case for each hexagon
Cost[(d,'S0')] = 1E4
sql = '''
SELECT DISTINCT facility_id, min_population_constraint, max_population_constraint
FROM facilities
ORDER BY 1;
'''
# the facility capacity constraint - cannot supply more hexagons than the facility has capacity for
for row in c.execute(sql):
[facility_id, min_population_constraint, max_population_constraint] = row
s = 'S{}'.format(facility_id)
Supply[s] = {}
Supply[s]['min_pop'] = min_population_constraint
Supply[s]['max_pop'] = max_population_constraint
sql = '''
SELECT facility_id, hexagon_id, distance
FROM od_matrix;
'''
# creating the Cost matrix
for row in c.execute(sql):
(facility_id, hexagon_id, distance) = row
d = 'D{}'.format(hexagon_id)
s = 'S{}'.format(facility_id)
Cost[(d,s)] = distance
print('Building the model')
# creating the model
model = ConcreteModel()
model.dual = Suffix(direction=Suffix.IMPORT)
# Step 1: Define index sets
dem = list(Demand.keys())
sup = list(Supply.keys())
# Step 2: Define the decision variables
model.x = Var(dem, sup, domain=NonNegativeReals)
# Step 3: Define Objective
model.Cost = Objective(
expr = sum([Cost[d,s]*model.x[d,s] for d in dem for s in sup]),
sense = minimize
)
# Step 4: Constraints
model.sup = ConstraintList()
# each facility cannot supply more than its population capacity
for s in sup:
model.sup.add(sum([model.x[d,s] for d in dem]) >= Supply[s]['min_pop'])
model.sup.add(sum([model.x[d,s] for d in dem]) <= Supply[s]['max_pop'])
model.dmd = ConstraintList()
# each hexagon can only be assigned to one facility
for d in dem:
model.dmd.add(sum([model.x[d,s] for s in sup]) == Demand[d])
'''
There is no need to add a constraint for the service/trade area distances
for the facilities. We are already handling this when we generate
the origin destination matrix. If any hexagon falls outside of all
facility trade areas, then it gets assigned to the "unassigned" facility.
'''
print('Solving the model')
# use the CBC solver and solve the model
results = SolverFactory('cbc').solve(model)
# for c in dem:
# for s in sup:
# print(c, s, model.x[c,s]())
# if the model solved correctly
if 'ok' == str(results.Solver.status):
print("Objective Function = ", model.Cost())
print('Outputting the results to GeoJSON') # note it would be faster to write the results directly to a database, e.g. Postgres / SQL Server
# print("Results:")
for s in sup:
for d in dem:
if model.x[d,s]() > 0:
# print("From ", s," to ", d, ":", model.x[d,s]())
facility_id = s.replace('S','')
hexagon_id = d.replace('D','')
c.execute('''
UPDATE od_matrix
SET optimal = 1
WHERE facility_id = {}
AND hexagon_id = {};
'''.format(facility_id, hexagon_id))
# create an empty feature collection for the results
feature_collection = { "type": "FeatureCollection", "features": [] }
rows = c.execute('''
SELECT facility_id, hexagon_id, facility_x, facility_y, hexagon_x, hexagon_y, distance
FROM od_matrix
WHERE optimal = 1;
''')
for row in rows:
(facility_id, hexagon_id, facility_x, facility_y, hexagon_x, hexagon_y, distance) = row
(hexagon_lat, hexagon_lng) = TRAN_3857_TO_4326.transform(hexagon_x, hexagon_y)
(facility_lat, facility_lng) = TRAN_3857_TO_4326.transform(facility_x, facility_y)
coords = [[facility_lng, facility_lat],[hexagon_lng, hexagon_lat]]
f = {
"type": "Feature",
"properties": {},
"geometry": {
"type": "LineString",
"coordinates": coords
}
}
# add the route to the feature collection
feature_collection['features'].append(f)
# output the optimally assigned pairs to geojson
with open("output/optimal_results.geojson", "w") as output_file:
output_file.write(json.dumps(feature_collection))
# update the hexagons
with open('output/hexagons.geojson') as json_file:
geojson = json.load(json_file)
hexagons = geojson['features']
for hexagon in hexagons:
facility_id = -1
hexagon_id = hexagon['properties']['id']
sql = '''
SELECT facility_id
FROM od_matrix
WHERE hexagon_id = {}
AND optimal = 1;
'''.format(hexagon_id)
row = c.execute(sql).fetchone()
if row:
facility_id = row[0]
hexagon['properties']['facility_id'] = facility_id
# load the area of interest into a JSON object
with open("output/hexagon_results.geojson", "w") as output_file:
output_file.write(json.dumps(geojson))
else:
print("No Feasible Solution Found")
finish_time = datetime.datetime.now()
# print('demand nodes = {} and supply nodes = {} | time = {}'.format(no_of_demand_nodes, no_of_supply_nodes, finish_time-start_time))
conn.commit()
conn.close()
for i in range(1,2):
os.system('afplay /System/Library/Sounds/Glass.aiff')