Моделирование размещения хабов в pyomo / forpes.ru

Главная
Моделирование размещения хабов в pyomo

Моделирование размещения хабов в pyomo +3

06.11.2023 08:50

Lozkins 7 825 Источник

Задача размещения хабов (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$

Общие ограничения (1) и (2)

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

Минимизация затрат на установку хаба и транспортные затраты:

$\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}$ - вещественная переменная, поток из узла через хабы и ;

Ограничения

Балансирование потоков (flow conservation constraint). Рассмотрим входящие и исходящие потоки хаба . Входящие потоки $O_ix_{ik}$ и $\sum_{l \in K} y_{ilk}$ , поток от отправителя напрямую и поток из первого хаба соответственно. Исходящие потоки $\sum_{j \in N} w_{ij}x_{jk}$ и $\sum_{l \in K} 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$

Общие ограничения (1) и (2)

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

Минимизация затрат на установку хаба и транспортные затраты:

$\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.

Ссылки

Различные наборы данных для задачи размещения хабов: OR Library by J.E. Beasley;
Обзор задач размещения хабов;
Предыдущая статья по хакам линейного программирования;
Ссылка на Jupyter Notebook.

Комментарии (7)

ZekaVasch
06.11.2023 11:28
#26129314
Заглавную картинку рисовала нейросеть или живой человек ?

Lozkins Автор
06.11.2023 11:28
#26129332
Живой человек )
1. rezdm
  06.11.2023 11:28
  #26129492
  Глядя на карту, сразу хочется спросить: "Крым чей"?
  1. Lozkins Автор
    06.11.2023 11:28
    #26129572
    Крымчан, очевидно

deadmoroz14
06.11.2023 11:28
#26133226
Я так и не понял, почему статья попала в хаб "машинное обучение". Нет тут его
1. Lozkins Автор
  06.11.2023 11:28
  #26133430
  К сожалению, у habra нет хаба для исследования операций/математического моделирования/математической оптимизации/индустриальной математики. Поэтому приходится выбирать смежные области, чтобы собрать свою целевую аудиторию из нескольких смежных хабов.
  1. deadmoroz14
    06.11.2023 11:28
    #26134726
    Понятно.
    
    Пора уже что-то такое вводить хабру. Далеко не первая статья по тематике

Моделирование размещения хабов в pyomo +3

О задаче

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

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

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

Индексы

Константы

Переменные

Ограничения

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

Переменные

Ограничения

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

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

Переменные

Ограничения

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

Данные

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

Результат

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

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

Заключение

Ссылки

Комментарии (7)

ZekaVasch

Lozkins Автор

rezdm

Lozkins Автор

deadmoroz14

Lozkins Автор

deadmoroz14