Моделирование размещения хабов в pyomo06.11.2023 12:00

Задача размещения хабов (Hub Location Problem) относится к стратегическому уровню планирования сети. Это накладывает ограничения на возможность оперативной реализации и валидации решения. Одним из способов моделирования и анализа такого рода решений без рисков для текущей сети является математическое моделирование.

О задаче

Транспортные, телекоммуникационные и компьютерные сети часто используют Hub-and-Spoke архитектуру для эффективной маршрутизации потоков между множеством отправителей и получателей. Особенность такой топологии заключается в использовании специального объекта сети — хаба. Хабом называется объект сети, который обеспечивает распределение, соединение, переключение, консолидацию, сортировку или перевалку в распределенных системах много-ко-многим. Кроме того, хабы позволяют соединить большой набор пар отправитель/получатель с использованием небольшого кол-ва соединений.

Сети без использования хабов и с хабами

Задача размещения хабов состоит в размещении хабов и построении сети хабов (построение соединений/ребер сети) с целью оптимизации затрат/уровня сервиса.

Выделяют следующие преимущества использования сетей с хабами:

Экономия на масштабе (экономика масштаба). Консолидация трафика на хабах позволяет сформировать большие потоки между хабами, что приводит к меньшим издержкам на трансфер (например, более плотное заполнение транспорта в перевозке сборных грузов);
Сокращение затрат на установку соединений между «непопулярными» или труднодоступными направлениями;
Повышение уровня сервиса за счет более регулярных отправок;
Сокращение затрат на управление сетью.

Задача является актуальной в транспортных, компьютерных и телекоммуникационных сетях. В частности, авиа перевозки пассажиров или грузов, почтовые отправки, перевозки сборных грузов, службы экстренного реагирования, сети маркетплэйсов и другие. Пример хабов у российских авиакомпаний можно посмотреть здесь.

Предположения

В 1986 году Мортон О'Келли впервые сформулировал математическую постановку задачи в виде задачи квадратичного программирования. Академическая среда подхватила ее и нагенерировала множество вариаций задачи, подходов к ее решению, областей применения, и структурировала условия задачи. Основными предположениями задачи являются:

Удовлетворение всего объема спроса;
Запрет на прямые соединениями между не хабовыми узлами сети;
Матрица расстояний (затрат) удовлетворяет неравенству треугольника;
Стоимость организации ребра между узлами не учитывается;
Эффект масштаба моделируется в виде фиксированной скидки на прохождение потока между хабами (часто обозначается как $\alpha$ ).

Следствием предположения (3) является то, что поток будет проходит через не более чем два хаба.

Здесь рассмотрим одну из основных версий задачи, где каждый узел сети может быть связан только с одним хабом (single allocation) и целевая функция учитывает фиксированные затраты на установку хабов в узлах сети.

Математическая модель

Проведение практических экспериментов по размещению хабов в сети на реальных объектах безумно дорого и долго апробируемо. Поэтому используются различного рода модели, которые позволяют оценить перспективы принимаемого решения. Математическая модель является одним из вариантов такого моделирования.

В задаче размещения хабов (Hub Location Problem) существует три типа моделей в виде задачи смешанного целочисленного линейного программирования. Предлагаю рассмотреть две, наиболее жизнеспособных и популярных. Модели отличаются решающими переменными. В одном случае переменная имеет четыре индекса (Campbell 1994, Skorin-Kapov et al., 1996), задача получается размерности O(n^4) переменных и O(n^3) ограничений, где — кол-во узлов в сети. В другом случае, решающая переменная имеет 3 индекса (Ernst and Krishnamoorthy 1996), задача получается размерности O(n^3) переменных и O(n^2) ограничений.

4-индексная модель уступает по размерности и скорости поиска оптимального решения 3-индексной. Но 4-индексная модель позволяет использовать различного рода разложения, в частности, разложение Бендерса. В некоторых случаях разложение задачи значительно изменяет баланс чаши весов в пользу 4-индексной модели.

Перейдем к математическим постановкам задачи. Обе модели имеют пару схожих ограничений и переменных, поэтому сформулируем общую часть отдельно.

Общие обозначения

Индексы

$N=\{1, \dots, n\}$ — множество узлов в сети;

$K \subseteq N$ — множество узлов в сети, которые могут быть выбраны в качестве хаба;

Константы

$w_{ij}$ — поток из узла в узел ;

f_i — установка хаба в узле ;

$d_{ij}$ — расстояние между узлами и ;

$\chi$ — стоимость консолидации потоков за единицу расстояния и единицу потока;

$\alpha$ — стоимость трансфера потоков за единицу расстояния и единицу потока (между хабами);

$\delta$ — стоимость распределения потоков за единицу расстояния и единицу потока;

$O_i = \sum_{j \in N} w_{ij}$ — общий исток из узла ;

$D_i = \sum_{j \in N} w_{ji}$ — общий сток в узел ;

$F_{ijkl} = w_{ij} (\chi d_{ik} + \alpha d_{kl} + \delta d_{lj})$ — затраты на транспортировку потока по маршруту из в через хабы и ;

Размещение хабов в сети является целесообразным, когда $\alpha < \chi$ и $\alpha < \delta$ . Здесь $\alpha$ моделирует эффект масштаба.

Переменные

$z_{kk}$ — бинарная переменная, принимает значение 1, если выбран в качестве хаба;

$z_{ik}$ — бинарная переменная, привязка узла сети к хабу ;

Ограничения

Каждый узел сети связан только с одним хабом:

$\sum_{k \in K} z_{ik} = 1, \quad \forall i \in N$

Связь узловой точки с хабом возможна, только если узловая точка выбрана в качестве хаба:

$z_{ik} \le z_{kk}, \quad \forall i \in N, k \in K$

Модель 4-индекса

Постановка задачи с использованием переменных с четырьмя индексами (Skorin-Kapov и др. , 1996).

Переменные

$x_{ijkl}$ -вещественная переменная, доля всего потока из в , проходящая через хабы и ;

Ограничения

$\sum_{l \in K} x_{ijkl} = z_{ik}, \quad \forall i, j \in N, k \in K$

$\sum_{k \in K} x_{ijkl} = z_{jl}, \quad \forall i, j \in N, l \in K$

Целевая функция

$\min \sum_k f_{k} z_{kk} + \sum_{i, j \in N} \sum_{k,l \in K} F_{ijkl}x_{ijkl}$

Модель 3-индекса

Постановка задачи с использованием переменных с тремя индексами (Ernst A. T., Krishnamoorthy M., 1996)

Переменные

$y_{ikl}$ — вещественная переменная, поток из узла через хабы и ;

Ограничения

$O_ix_{ik} + \sum_{l \in K} y_{ilk} = \sum_{j \in N} w_{ij}x_{jk} + \sum_{l \in K} y_{ikl}, \quad \forall i \in N, k \in K$

Целевая функция

$\min \sum_k f_{k} z_{kk} + \sum_{i \in N, k \in K} (\chi O_i + \delta D_i)d_{ik}z_{ik} + \sum_{k,l \in K} \alpha d_{kl}y_{ikl}$

Данные

Задачу размещения хабов можно решать методами линейного программирования, квадратичного программирования, эвристическими (например, GRASP) или мета-эвристическими методами. Оценку производительности и проведение сравнительного анализа подходов к решению между собой необходимо проводить на фиксированном наборе данных (тестовый набор). Кроме того, новые концепции задачи также сравнивают с существующими, а теоретические выкладки подтверждают практическими экспериментами на этих данных. В научной среде в этих целях часто используются следующие наборы:

CAB (Civil Aeronautic Board). Данные совета гражданской авиации США: пассажиропоток авиакомпаний между 25 городами США, 1970 год.
AP (Australia Post). Почтовые отправки между 200 почтовых округов в Сиднее, Австралия.
TR (Turkish Network). Потоки между 81 городом Турции; сформированы на основе численности населения.

Наборы данных можно найти в OR Library by J.E. Beasley.

Программная реализация моделей

Компьютерную реализацию математических моделей осуществим в python посредством библиотеки Pyomo. Pyomo — это программный пакет с открытым исходным кодом на базе Python, который поддерживает разнообразный набор возможностей оптимизации для формулирования, решения и анализа оптимизационных моделей. Решение оптимизационной задачи отдадим на откуп солверу CBC, также с открытым исходным кодом.

Рассмотрим реализацию моделей для 4-индексной и 3-индексной версий модели размещения хабов в сети с фиксированными затратами на установку хабов и одинарной привязкой узла сети к хабу (single allocation). Моделирование проведу на данных CAB. Так как в выбранном наборе данных нет информации по стоимости оборудования узла сети до уровня хаба, то эти данные сгенерируем искусственно (десятичный логарифм исходящего потока из узла). Дополнительно нормируем потоки сети так, чтобы сумма всех потоков равнялась единице.

# Установка pyomo и солвера cbc
%pip install -q pyomo
%apt-get install -y -qq coinor-cbc

Код подготовки данных

# Обработка и расширение данных
import pandas as pd
import numpy as np

# Загрузка данных Civil Aeronautic Boards
df_cab = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/04_hlp_cab.csv", sep=";", encoding="cp1251")

# Стоимость установки хаба
df_f = df_cab.groupby("i", as_index=False)["w"].sum()
df_f = df_f.set_index(["i"])
df_f["f"] = 15 * np.log10(df_f["w"])
dct_f = df_f.to_dict()["f"]

# Нормируем размер пассажиропотоков
df_cab["w"] = df_cab["w"] / df_cab["w"].sum()

# Коэффициент эффекта масштаба за трансфер между хабами
alpha = 0.4

Код 4-индексной версии модели

# 4-индексная модель размещения хабов в сети
from pyomo.environ import *
from pyomo.opt import SolverStatus, TerminationCondition

# Объединяем затраты на единицу потока по маршруту i-k-l-j
df_tmp = df_cab[["i", "j", "c"]].copy()
df_F = df_tmp.set_axis(["i", "k", "c_ik"], axis=1).merge(df_tmp.set_axis(["l", "j", "c_lj"], axis=1), how="cross")
df_F = df_F.merge(df_tmp.set_axis(["k", "l", "c_kl"], axis=1), how="inner", on=["k", "l"])
df_F = df_F.merge(df_cab, how="inner", on=["i", "j"])
df_F["F"] = df_F["w"] * (df_F["c_ik"] + alpha * df_F["c_kl"] + df_F["c_lj"])
df_F = df_F.set_index(["i", "j", "k", "l"])
dct_F = df_F.to_dict()["F"]

# Инициализация модели
model = ConcreteModel()

# Инициализация узлов сети
model.nodes = RangeSet(1, df_cab["i"].max())

# Инициализация переменных
model.z = Var(model.nodes, model.nodes, initialize=0, within=Binary)
model.x = Var(model.nodes, model.nodes, model.nodes, model.nodes, initialize=0, bounds=(0,1), within=Reals)

# Ограничение: Каждый узел сети связан только с одним хабом
def rule_single_hub(model, i):
  return sum(model.z[i, k] for k in model.nodes) == 1
model.constr_sh = Constraint(model.nodes, rule=rule_single_hub)

# Ограничение: Связь узловой точки с хабом возможна, только если узловая точка выбрана в качестве хаба
def rule_if_hub(model, i, k):
  return model.z[i, k] <= model.z[k, k]
model.constr_ih = Constraint(model.nodes, model.nodes, rule=rule_if_hub)

# Ограничение: Маршрут с узла i должен проходить через связанный хаб
def rule_route_hub_first(model, i, j, k):
  return sum(model.x[i, j, k, l] for l in model.nodes) == model.z[i, k]
model.constr_rhf = Constraint(model.nodes, model.nodes, model.nodes, rule=rule_route_hub_first)

# Ограничение:  Маршрут в узел j должен проходить только через связанный с этим узлом хаб
def rule_route_hub_second(model, i, j, l):
  return sum(model.x[i, j, k, l] for k in model.nodes) == model.z[j, l]
model.constr_rhs = Constraint(model.nodes, model.nodes, model.nodes, rule=rule_route_hub_second)

# Целевая функция
def rule_obj(model):
  sum_route_costs = sum(dct_F[i, k, l, j] * model.x[i, k, l, j] for i, j, k, l in dct_F)
  sum_hub_install_costs = sum(dct_f[k] * model.z[k, k] for k in dct_f)
  return sum_route_costs + sum_hub_install_costs
model.obj = Objective(rule=rule_obj, sense=minimize)

# Инициализация солвера и решение задачи
solver = SolverFactory('cbc', executable='/usr/bin/cbc')
results = solver.solve(model, tee=True)

if (results.solver.status == SolverStatus.ok) and (results.solver.termination_condition == TerminationCondition.optimal):
  # Транспортные затраты
  trans_costs = value(sum(dct_F[i, k, l, j] * model.x[i, k, l, j] for i, j, k, l in dct_F))
  print("Найдено оптимальное решение задачи.")
  print(f"Значение целевой функции: {round(value(model.obj), 2)}")
  print(f"Транспортные затраты: {round(trans_costs, 2)}")
elif results.solver.termination_condition == TerminationCondition.infeasible:
  print("Нет решения задачи: infeasible")
else:
  print(str(results.solver))

# Извлечение результата
df_cab["allocation"] = df_cab.apply(lambda x: value(model.z[x.i, x.j]), axis=1)

# Привязка узлов к хабам
df_allocation = df_cab[df_cab["allocation"] == 1]

# Выбранные хабы
hubs = df_allocation["j"].unique()  # 450 sec

Код 3-индексной версии модели

# 3-индексная модель размещения хабов в сети
from pyomo.environ import *
from pyomo.opt import SolverStatus, TerminationCondition

# Подготовка данных O_i, D_j и преобразование спроса и затрат в словари
df_o = df_cab.groupby("i")["w"].sum()
dct_o = df_o.to_dict()
df_d = df_cab.groupby("j")["w"].sum()
dct_d = df_d.to_dict()
df_w = df_cab.set_index(["i", "j"])
dct_w = df_w.to_dict()["w"]
dct_c = df_w.to_dict()["c"]

# Инициализация модели
model = ConcreteModel()

# Инициализация узлов сети
model.nodes = RangeSet(1, df_cab["i"].max())

# Инициализация переменных
model.z = Var(model.nodes, model.nodes, initialize=0, within=Binary)
model.y = Var(model.nodes, model.nodes, model.nodes, initialize=0, bounds=(0,1), within=Reals)

# Ограничение: Каждый узел сети связан только с одним хабом
def rule_single_hub(model, i):
  return sum(model.z[i, k] for k in model.nodes) == 1
model.constr_sh = Constraint(model.nodes, rule=rule_single_hub)

# Ограничение: Связь узловой точки с хабом возможна, только если узловая точка выбрана в качестве хаба
def rule_if_hub(model, i, k):
  return model.z[i, k] <= model.z[k, k]
model.constr_ih = Constraint(model.nodes, model.nodes, rule=rule_if_hub)

# Ограничение: Баланс потоков
def rule_balance(model, i, k):
  lhs = sum(model.y[i, k, l] for l in model.nodes) - sum(model.y[i, l, k] for l in model.nodes)
  rhs = dct_o[i] * model.z[i, k] - sum(dct_w[i, j] * model.z[j, k] for j in model.nodes)
  return lhs == rhs
model.constr_b = Constraint(model.nodes, model.nodes, rule=rule_balance)

# Целевая функция
def rule_obj(model):
  sum_to_hub_costs = sum((dct_o[i] + dct_d[i]) * dct_c[i, k] * model.z[i, k] for i, k in dct_w)
  sum_between_hubs = sum(alpha * dct_c[k, l] * model.y[i, k, l] for i, k in dct_w for l in model.nodes)
  sum_hub_install_costs = sum(dct_f[k] * model.z[k, k] for k in dct_f)
  return sum_hub_install_costs + sum_to_hub_costs + sum_between_hubs
model.obj = Objective(rule=rule_obj, sense=minimize)

# Инициализация солвера и решение задачи
solver = SolverFactory('cbc', executable='/usr/bin/cbc')
results = solver.solve(model, tee=True)

if (results.solver.status == SolverStatus.ok) and (results.solver.termination_condition == TerminationCondition.optimal):
  # Транспортные затраты
  sum_to_hub_costs = value(sum((dct_o[i] + dct_d[i]) * dct_c[i, k] * model.z[i, k] for i, k in dct_w))
  sum_between_hubs = value(sum(alpha * dct_c[k, l] * model.y[i, k, l] for i, k in dct_w for l in model.nodes))
  trans_costs = sum_to_hub_costs + sum_between_hubs
  print("Найдено оптимальное решение задачи.")
  print(f"Значение целевой функции: {round(value(model.obj), 2)}")
  print(f"Транспортные затраты: {round(trans_costs, 2)}")
elif results.solver.termination_condition == TerminationCondition.infeasible:
  print("Нет решения задачи: infeasible")
else:
  print(str(results.solver))

# Извлечение результата
df_cab["allocation"] = df_cab.apply(lambda x: value(model.z[x.i, x.j]), axis=1)

# Привязка узлов к хабам
df_allocation = df_cab[df_cab["allocation"] == 1]

# Выбранные хабы
hubs = df_allocation["j"].unique()  # 18 sec

Результат

Оптимальное решение для 4-индексной и 3-индексной модели получилось одинаковым, но за разное время: ~450 и ~18 секунд соответственно. Значения целевых функций 1133.54, а транспортные затраты с учетом дисконта $\alpha$ равны 794.56, что соответствует значениям других исследователей. Оптимальный набор хабов при $\alpha = 4$ — это узлы $\{1, 18, 4, 12\}$ .

Размерность задачи для 4-индексной версии — 31875 ограничений и 391250 переменных. Размерность задачи для 3-индексной версии — 1875 ограничений и 15625 переменных. При таком явном преимуществе 3-индексной модели 4-индексная версия по прежнему активно используется в теории и на практике, благодаря возможности применять эвристики и различного рода разложения.

Решение задачи размещения хабов на наборе данных CAB

Эффект масштаба

Одним из предположений классической постановки задачи размещения хабов является постоянная скидка на трансфер потоков между хабами. Это сильное предположение, которое не сочетается с реальностью. Одним из направлений в этой области является уход от классической постановки учета скидки (постоянный дисконт) к другим моделям:

Нелинейная функция или кусочно-линейная функция скидки;
Использование пороговых значений (кусочно-постоянная функция);
Учет в модели дискретных единиц транспорта;
Разложение функции затрат на постоянную и переменную часть.

Вариации постановок

Классическая постановка задачи сформулирована в 1990-ых гг., с того момента область значительно расширилась. Здесь можно выделить два направления развития: добавление правил в рамках классической постановки и концептуальный подход к моделированию. К первому относятся такие модификации как:

Ограничение на узлы, которые потенциально могут стать хабами. Учет пропускной способности хабов;
Ограничения на связи: связь с одним, несколькими или r-связей с хабами (single-, multiple-, r-allocation). Допуск прямых соединений, другие ограничения на маршрутизацию потоков;
Сеть хабов является полным/неполным графом, формирует топологию типа звезда, кольца, древовидную или др. Учет различных типов соединений (разный тип транспорта, кабеля).

К концептуальным вариациям относятся сами подходы моделирования:

Стохастические модели (спрос, затраты, время);
Робастное размещение хабов в сети;
Динамическое/многопериодное размещение хабов;
Размещение хабов с учетом заторов;
Размещение хабов в условиях конкуренции;
Надежное размещение хабов в сети;
Размещение хабов с целевой функцией максимизации прибыли;
Размещение хабов и маршрутизация.

Заключение

Задача размещения хабов относится к стратегическому уровню планирования. Построение решения и его проверка на реальной сети затратно по времени и финансам. Поэтому для проработки решений используют различные методы моделирования. В статье рассмотрели классическую постановку задачи (MIP) и ее реализацию в среде Pyomo на наборе данных Civil Aeronautic Board (CAB), оптимизационную задачу решили с помощью солвера cbc.

Моделирование размещения хабов в pyomo06.11.2023 12:00

О задаче

Предположения

Математическая модель

Общие обозначения

Индексы

Константы

Переменные

Ограничения

Модель 4-индекса

Переменные

Ограничения

Целевая функция

Модель 3-индекса

Переменные

Ограничения

Целевая функция

Данные

Программная реализация моделей

Результат

Эффект масштаба

Вариации постановок

Заключение

Ссылки