Географическая кластеризация

Это руководство было изначально предоставлено Мэтью Хелмом (Matthew Helm) и Матье Танно (Mathieu Tanneau).

Цель этого упражнения — объединить городов в групп, минимизировав общее попарное расстояние между городами и обеспечив относительно небольшую дисперсию общей численности населения в каждой группе.

В этом руководстве используются следующие пакеты.

using JuMP
import DataFrames
import HiGHS
import LinearAlgebra

Для этого примера мы возьмем 20 самых густонаселенных городов США.

cities = DataFrames.DataFrame(
    Union{String,Float64}[
        "New York, NY" 8.405 40.7127 -74.0059
        "Los Angeles, CA" 3.884 34.0522 -118.2436
        "Chicago, IL" 2.718 41.8781 -87.6297
        "Houston, TX" 2.195 29.7604 -95.3698
        "Philadelphia, PA" 1.553 39.9525 -75.1652
        "Phoenix, AZ" 1.513 33.4483 -112.0740
        "San Antonio, TX" 1.409 29.4241 -98.4936
        "San Diego, CA" 1.355 32.7157 -117.1610
        "Dallas, TX" 1.257 32.7766 -96.7969
        "San Jose, CA" 0.998 37.3382 -121.8863
        "Austin, TX" 0.885 30.2671 -97.7430
        "Indianapolis, IN" 0.843 39.7684 -86.1580
        "Jacksonville, FL" 0.842 30.3321 -81.6556
        "San Francisco, CA" 0.837 37.7749 -122.4194
        "Columbus, OH" 0.822 39.9611 -82.9987
        "Charlotte, NC" 0.792 35.2270 -80.8431
        "Fort Worth, TX" 0.792 32.7554 -97.3307
        "Detroit, MI" 0.688 42.3314 -83.0457
        "El Paso, TX" 0.674 31.7775 -106.4424
        "Memphis, TN" 0.653 35.1495 -90.0489
    ],
    ["city", "population", "lat", "lon"],
)

20×4 DataFrame

Row	city	population	lat	lon
	Union…	Union…	Union…	Union…
1	New York, NY	8.405	40.7127	-74.0059
2	Los Angeles, CA	3.884	34.0522	-118.244
3	Chicago, IL	2.718	41.8781	-87.6297
4	Houston, TX	2.195	29.7604	-95.3698
5	Philadelphia, PA	1.553	39.9525	-75.1652
6	Phoenix, AZ	1.513	33.4483	-112.074
7	San Antonio, TX	1.409	29.4241	-98.4936
8	San Diego, CA	1.355	32.7157	-117.161
9	Dallas, TX	1.257	32.7766	-96.7969
10	San Jose, CA	0.998	37.3382	-121.886
11	Austin, TX	0.885	30.2671	-97.743
12	Indianapolis, IN	0.843	39.7684	-86.158
13	Jacksonville, FL	0.842	30.3321	-81.6556
14	San Francisco, CA	0.837	37.7749	-122.419
15	Columbus, OH	0.822	39.9611	-82.9987
16	Charlotte, NC	0.792	35.227	-80.8431
17	Fort Worth, TX	0.792	32.7554	-97.3307
18	Detroit, MI	0.688	42.3314	-83.0457
19	El Paso, TX	0.674	31.7775	-106.442
20	Memphis, TN	0.653	35.1495	-90.0489

Особенности модели

Мы разделим эти 20 городов на 3 группы и предположим, что идеальная или целевая численность населения $P$ для группы — это просто общая численность населения 20 городов, деленная на 3:

n = size(cities, 1)
k = 3
P = sum(cities.population) / k

11.038333333333334

Получение расстояний между всеми городами

Вычислим попарные расстояния по формуле гаверсинусов между всеми городами в наборе данных и сохраним результат в переменной dm:

"""
    haversine(lat1, long1, lat2, long2, r = 6372.8)


Compute the haversine distance between two points on a sphere of radius `r`,
where the points are given by the latitude/longitude pairs `lat1/long1` and
`lat2/long2` (in degrees).
"""

function haversine(lat1, long1, lat2, long2, r = 6372.8)
    lat1, long1 = deg2rad(lat1), deg2rad(long1)
    lat2, long2 = deg2rad(lat2), deg2rad(long2)
    hav(a, b) = sin((b - a) / 2)^2
    inner_term = hav(lat1, lat2) + cos(lat1) * cos(lat2) * hav(long1, long2)
    d = 2 * r * asin(sqrt(inner_term))
    # Округляем расстояние до ближайшего километра.
    return round(Int, d)
end

haversine (generic function with 2 methods)

Матрица расстояний симметрична, поэтому мы преобразуем ее в матрицу LowerTriangular, чтобы лучше интерпретировать значение целевой функции модели:

dm = LinearAlgebra.LowerTriangular([
    haversine(cities.lat[i], cities.lon[i], cities.lat[j], cities.lon[j])
    for i in 1:n, j in 1:n
])

20×20 LinearAlgebra.LowerTriangular{Int64, Matrix{Int64}}:
    0     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅  …     ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
 3937     0     ⋅     ⋅     ⋅     ⋅     ⋅        ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
 1145  2805     0     ⋅     ⋅     ⋅     ⋅        ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
 2282  2207  1516     0     ⋅     ⋅     ⋅        ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
  130  3845  1068  2157     0     ⋅     ⋅        ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
 3445   574  2337  1633  3345     0     ⋅  …     ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
 2546  1934  1695   304  2423  1363     0        ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
 3908   179  2787  2094  3812   481  1813        ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
 2206  1993  1295   362  2089  1424   406        ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
 4103   492  2958  2588  4023   989  2336        ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
 2432  1972  1577   235  2310  1398   118  …     ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
 1036  2907   265  1394   938  2409  1609        ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
 1345  3450  1391  1321  1221  2883  1626        ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
 4130   559  2986  2644  4052  1051  2394        ⋅     ⋅     ⋅     ⋅     ⋅  ⋅
  767  3177   444  1598   668  2679  1834        0     ⋅     ⋅     ⋅     ⋅  ⋅
  855  3405   946  1490   725  2863  1777  …   560     0     ⋅     ⋅     ⋅  ⋅
 2251  1944  1327   382  2134  1375   387     1511  1543     0     ⋅     ⋅  ⋅
  774  3186   382  1780   711  2716  1994      264   813  1646     0     ⋅  ⋅
 3054  1130  2010  1081  2945   559   804     2292  2398   864  2374     0  ⋅
 1534  2576   777   780  1415  2028  1017      820   837   722  1003  1565  0

Построение модели

Разобравшись с основами, мы можем настроить модель, создать переменные решения, добавить ограничения, а затем решить ее.

Сначала настроим модель, использующую решатель Cbc. Далее настроим двоичную переменную $x_{i,k}$, которая принимает значение $1$, если город $i$ находится в группе $k$, и $0$ в противном случае. Каждый город должен входить в группу, поэтому мы добавим ограничение $\sum_k x_{i,k} = 1$ для каждого $i$.

model = Model(HiGHS.Optimizer)
set_silent(model)
@variable(model, x[1:n, 1:k], Bin)
@constraint(model, [i = 1:n], sum(x[i, :]) == 1);
# Для уменьшения симметрии мы решаем, что первый город должен принадлежать первой группе.
fix(x[1, 1], 1; force = true)

Общая численность населения в группе $k$ равна $Q_k = \sum_ix_{i,k}q_i$, где $q_i$ — это просто $i$-е значение из столбца population в DataFrame cities. Добавим ограничение, согласно которому $\alpha \leq (Q_k - P) \leq \beta$. Установим $\alpha$ равным $-3$ миллионам, а $\beta$ — равным $3$. Изменяя эти пороговые значения, можно заметить, что существует компромисс между относительно равномерной численностью населения в группах и географически близким расположением городов в каждой группе. Иными словами, чем больше абсолютные значения $\alpha$ и $\beta$, тем ближе друг к другу будут города в группе, но тем выше будет разница между численностью населения в группах.

@variable(model, -3 <= population_diff[1:k] <= 3)
@constraint(model, population_diff .== x' * cities.population .- P)

3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 -8.405 x[1,1] - 3.884 x[2,1] - 2.718 x[3,1] - 2.195 x[4,1] - 1.553 x[5,1] - 1.513 x[6,1] - 1.409 x[7,1] - 1.355 x[8,1] - 1.257 x[9,1] - 0.998 x[10,1] - 0.885 x[11,1] - 0.843 x[12,1] - 0.842 x[13,1] - 0.837 x[14,1] - 0.822 x[15,1] - 0.792 x[16,1] - 0.792 x[17,1] - 0.688 x[18,1] - 0.674 x[19,1] - 0.653 x[20,1] + population_diff[1] = -11.038333333333334
 -8.405 x[1,2] - 3.884 x[2,2] - 2.718 x[3,2] - 2.195 x[4,2] - 1.553 x[5,2] - 1.513 x[6,2] - 1.409 x[7,2] - 1.355 x[8,2] - 1.257 x[9,2] - 0.998 x[10,2] - 0.885 x[11,2] - 0.843 x[12,2] - 0.842 x[13,2] - 0.837 x[14,2] - 0.822 x[15,2] - 0.792 x[16,2] - 0.792 x[17,2] - 0.688 x[18,2] - 0.674 x[19,2] - 0.653 x[20,2] + population_diff[2] = -11.038333333333334
 -8.405 x[1,3] - 3.884 x[2,3] - 2.718 x[3,3] - 2.195 x[4,3] - 1.553 x[5,3] - 1.513 x[6,3] - 1.409 x[7,3] - 1.355 x[8,3] - 1.257 x[9,3] - 0.998 x[10,3] - 0.885 x[11,3] - 0.843 x[12,3] - 0.842 x[13,3] - 0.837 x[14,3] - 0.822 x[15,3] - 0.792 x[16,3] - 0.792 x[17,3] - 0.688 x[18,3] - 0.674 x[19,3] - 0.653 x[20,3] + population_diff[3] = -11.038333333333334

Осталось добавить в модель последнюю двоичную переменную $z_{i,j}$, которая будет использоваться для вычисления общего расстояния между городами в группах по формуле $\sum_{i,j}d_{i,j}z_{i,j}$. Переменная $z_{i,j}$ будет равна $1$, если города $i$ и $j$ находятся в одной группе, и $0$ в противном случае.

Чтобы равенство $z_{i,j} = 1$ соблюдалось только тогда, когда города $i$ и $j$ находятся в одной группе, мы добавим ограничения $z_{i,j} \geq x_{i,k} + x_{j,k} - 1$ для каждой пары $i,j$ и каждого $k$:

@variable(model, z[i = 1:n, j = 1:i], Bin)
for k in 1:k, i in 1:n, j in 1:i
    @constraint(model, z[i, j] >= x[i, k] + x[j, k] - 1)
end

Теперь можно добавить в модель цель, которая будет заключаться в минимизации скалярного произведения $z$ и матрицы расстояний dm.

@objective(model, Min, sum(dm[i, j] * z[i, j] for i in 1:n, j in 1:i));

После этого можно вызвать optimize! и просмотреть результаты.

optimize!(model)
@assert is_solved_and_feasible(model)

Просмотр результатов

Получив результаты, мы можем добавить столбец группы в DataFrame cities и выполнить перебор по переменной $x$, чтобы назначить группу каждому городу. После этого мы можем посмотреть общую численность населения и города в каждой группе, чтобы убедиться в том, что они сгруппированы по географической близости.

cities.group = zeros(n)

for i in 1:n, j in 1:k
    if round(Int, value(x[i, j])) == 1
        cities.group[i] = j
    end
end

for group in DataFrames.groupby(cities, :group)
    @show group
    println("")
    @show sum(group.population)
    println("")
end

group = 7×5 SubDataFrame
 Row │ city              population  lat      lon       group
     │ Union…            Union…      Union…   Union…    Float64
─────┼──────────────────────────────────────────────────────────
   1 │ New York, NY      8.405       40.7127  -74.0059      1.0
   2 │ Philadelphia, PA  1.553       39.9525  -75.1652      1.0
   3 │ Indianapolis, IN  0.843       39.7684  -86.158       1.0
   4 │ Jacksonville, FL  0.842       30.3321  -81.6556      1.0
   5 │ Columbus, OH      0.822       39.9611  -82.9987      1.0
   6 │ Charlotte, NC     0.792       35.227   -80.8431      1.0
   7 │ Detroit, MI       0.688       42.3314  -83.0457      1.0

sum(group.population) = 13.944999999999999

group = 7×5 SubDataFrame
 Row │ city             population  lat      lon       group
     │ Union…           Union…      Union…   Union…    Float64
─────┼─────────────────────────────────────────────────────────
   1 │ Chicago, IL      2.718       41.8781  -87.6297      2.0
   2 │ Houston, TX      2.195       29.7604  -95.3698      2.0
   3 │ San Antonio, TX  1.409       29.4241  -98.4936      2.0
   4 │ Dallas, TX       1.257       32.7766  -96.7969      2.0
   5 │ Austin, TX       0.885       30.2671  -97.743       2.0
   6 │ Fort Worth, TX   0.792       32.7554  -97.3307      2.0
   7 │ Memphis, TN      0.653       35.1495  -90.0489      2.0

sum(group.population) = 9.909

group = 6×5 SubDataFrame
 Row │ city               population  lat      lon       group
     │ Union…             Union…      Union…   Union…    Float64
─────┼───────────────────────────────────────────────────────────
   1 │ Los Angeles, CA    3.884       34.0522  -118.244      3.0
   2 │ Phoenix, AZ        1.513       33.4483  -112.074      3.0
   3 │ San Diego, CA      1.355       32.7157  -117.161      3.0
   4 │ San Jose, CA       0.998       37.3382  -121.886      3.0
   5 │ San Francisco, CA  0.837       37.7749  -122.419      3.0
   6 │ El Paso, TX        0.674       31.7775  -106.442      3.0

sum(group.population) = 9.261000000000001