p-median

p-median 알고리즘은 위치 할당을 위한 최적화 식입니다.

수요지 $i$ 와 공급 후보지 $j$ 가 있다고 할 때 다음과 같은 구조를 가집니다.

Inputs:

$h_{i}$ : $i$ 지역에 대한 수요량

$d_{i j}$ : 수요지 $i$ 와 공급지 $j$ 사이의 거리

$p$ : 선정할 최종 공급지의 개수

Decision Variables:

$x_{j}$
- 1 : $j$ 지점에 최종 공급지를 설치한다면
- 0 : 그렇지 않을 경우
$y_{i j}$
- 1 : $i$ 의 수요량이 가장 가까운 거리에 있는 $j$ 에 할당된다면
- 0 : 그렇지 않을 경우

$Min \sum_{i} \sum_{j} h_{i} d_{i j} y_{i j}$

Constraints:

$\begin{equation} \sum_j y_{ij} = 1 \quad (\forall i) \end{equation} \begin{equation} \sum_j x_j = p \end{equation} \begin{equation} y_{ij} - x_j \leq 0 \quad (\forall i, j) \end{equation} \begin{equation} x_j \in \{0, 1\} \quad (\forall j) \end{equation} \begin{equation} y_{ij} \in \{0, 1\} \quad (\forall j) \end{equation}$

(1) 모든 수요지점에서 반드시 한 개의 가장 가까운 공급지를 이용해야 한다.
(2) 선정한 최종 공급지의 개수는 $p$ 개이다.
(3) $j$ 지점에 공급지가 설치된 경우에만 수요지점에서 가장 가까운 공급지 $y_{ij}$ 가 활성화된다.
(4) $x_{j}$는 0 또는 1의 값이다.
(5) $y_{ij}$는 0 또는 1의 값이다.

코드

import math
import pulp

# Haversine formula to calculate distance between two points on the Earth
def haversine(coord1, coord2):
    R = 6371.0  # Earth radius in kilometers

    lat1, lon1 = coord1
    lat2, lon2 = coord2

    dlat = math.radians(lat2 - lat1)
    dlon = math.radians(lon2 - lon1)

    a = math.sin(dlat / 2) ** 2 + math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) * math.sin(dlon / 2) ** 2
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))

    distance = R * c
    return distance

# Define the problem
prob = pulp.LpProblem("p-median", pulp.LpMinimize)

# Example data
# Coordinates of demand points and candidate facility locations
demand_points = [(37.4907391344, 127.0154235196), (37.4917391344, 127.0164235196)]  # Example coordinates
candidate_facilities = [(37.4927391344, 127.0174235196), (37.4937391344, 127.0184235196)]

# Demand quantities for each demand point
demand_quantities = [100, 150]  # Example demand quantities

p = 1  # Number of facilities to select

# Define decision variables
x = pulp.LpVariable.dicts("facility", range(len(candidate_facilities)), cat='Binary')
y = pulp.LpVariable.dicts("assign", [(i, j) for i in range(len(demand_points)) for j in range(len(candidate_facilities))], cat='Binary')

# Define the objective function
prob += pulp.lpSum(
    pulp.lpSum(
        y[(i, j)] * demand_quantities[i] * haversine(demand_points[i], candidate_facilities[j]) 
        for j in range(len(candidate_facilities))
    ) 
    for i in range(len(demand_points))
)

# Define constraints
for i in range(len(demand_points)):
    prob += pulp.lpSum(y[(i, j)] for j in range(len(candidate_facilities))) == 1
prob += pulp.lpSum(x[j] for j in range(len(candidate_facilities))) == p
for i in range(len(demand_points)):
    for j in range(len(candidate_facilities)):
        prob += y[(i, j)] <= x[j]

# Solve the problem
prob.solve()

# Output the results
for v in prob.variables():
    print(v.name, "=", v.varValue)

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /opt/anaconda3/lib/python3.12/site-packages/pulp/solverdir/cbc/osx/64/cbc /var/folders/5_/ym1r3my96mb6yq3bdn1242bw0000gp/T/ae306e2e5033456f8aa6cccddd1af83b-pulp.mps -timeMode elapsed -branch -printingOptions all -solution /var/folders/5_/ym1r3my96mb6yq3bdn1242bw0000gp/T/ae306e2e5033456f8aa6cccddd1af83b-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 12 COLUMNS
At line 43 RHS
At line 51 BOUNDS
At line 58 ENDATA
Problem MODEL has 7 rows, 6 columns and 14 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 49.6805 - 0.00 seconds
Cgl0003I 0 fixed, 0 tightened bounds, 0 strengthened rows, 7 substitutions
Cgl0004I processed model has 0 rows, 0 columns (0 integer (0 of which binary)) and 0 elements
Cbc3007W No integer variables - nothing to do
Cuts at root node changed objective from 49.6805 to -1.79769e+308
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                49.68052983assign_(0,_0) = 1.0
assign_(0,_1) = 0.0
assign_(1,_0) = 1.0
assign_(1,_1) = 0.0
facility_0 = 1.0
facility_1 = 0.0