Промышленное математическое программирование - тема раскрученная в академической среде для стандартизированных случаев, но детали реальных внедрений раскрываются редко и спустя много лет.

В настоящей статье делюсь опытом разработки и внедрения в процессы компании оптимизационного решения на базе математического программирования. Материал расширил исследовательскими элементами и локальным мини benchmark'ом.

Бизнес задача

На входе имеем крупного производителя сетевого оборудования, у которого одним из этапов контроля качества продукции является тестирование пользовательских, нагрузочных или функциональных сценариев работы на собственных тестовых стендах.

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

Различные подразделения компании работают над улучшением определенных характеристик продукции или производят проверку сценариев ее использования. Каждый такой отдел разрабатывает свои тестовые сценарии и выполняет их на тестовом стенде.

Из чего состоит тест? Тест состоит из скрипта и требований к оборудованию, на котором он может быть запущен. Содержание скрипта теста - черный ящик, и этого понимания достаточно для решения задачи. В свою очередь, требования к оборудованию распадаются на топологию и специальные свойства оборудования/соединений. На основе этих данных происходит маппинг требований теста и участка тестового стенда для его выполнения.

Накладываем задачу на данные

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

Второй фактор, который нам не удалось побороть - ожидаемая длительность выполнения теста. Эта информация помогла бы разбивать тестовый стенд согласно ожидаемой нагрузке. Но раскрытие содержимого скрипта тестового сценария находилось в тени от нас. История (или прогнозное время выполнения тестового сценария) - не собиралась и не анализировалась.

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

Не всегда удается разделить стенд согласно указанным требованиям к покрытию тестов. Возникает вопрос надеженности модели: выдавать ошибку или выдать решение с нарушениями (мягкие ограничения)? Безусловно, второй вариант наиболее практичный с точки зрения пользователя, тем более у отделов есть маркировка - приоритет.

Размер НЕ имеет значения! Выполнение теста на одной единице оборудования равносильно выполнению теста на 5 единицах оборудования. Это облегчающее обстоятельство уводит нас от многокритериальной задачи. Такое допущение принято для упрощения материала статьи, в реальности иначе.

Наша интерпретация условий и наши возможности:

• Предложение предварительно рассчитать все или подмножество вариантов назначения теста на участки тестового стенда. Т.е. с каждым тестовым кейсом связать множество допустимых подсетей тестового стенда.

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

• Ограничение: необходимо разбить тестовый стенд таким образом, чтобы каждое подразделение могло выполнить на своем "куске" стенда определенный процент от своего набора тестов (условие минимального покрытия).

Мат. модель

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

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

Наш выбор пал на программирование в ограничениях. Потому что часть работ с нас снимает готовый solver, меньше времени требуется на разработку алгоритмов поиска решения и его тестирование. Это позволяет быстрее доставлять результаты бизнесу. Действовали по принципу 20 на 80.

Перейду непосредственно к постановке задачи в виде задачи целочисленного линейного программирования.

Индексы и множества

- множество тестовых сценариев;

- множество индексов подсетей тестового стенда, где может быть выполнен тест ;

- множество оборудования, из которого состоит тестовый стенд;

- множество департаментов компании, между которыми требуется разделить тестовый стенд;

- подмножество подсетей тестовых сценариев, которые принадлежат отделу и используют оборудование ;

Постоянные величины

- кол-во тестовых сценариев отдела , которые необходимо покрыть;

- весовой приоритет выполнения объемного покрытия тестовых сценариев для отдела ;

- минимальное кол-во единиц оборудования, которое нужно передать отделу ;

Решающие переменные

- бинарная переменная, принимает значение 1, если оборудование назначено отделу ; 0 в противном случае;

- бинарная переменная, принимает значение 1, если тест отдела может быть выполнен на участке тестового стенда закрепленным за отделом; 0 в противном случае;

- бинарная переменная, принимает значение 1, если все оборудование подсети теста отдела закреплено за отделом; 0 в противном случае;

- целочисленная переменная, объем нарушения уровня покрытия тестов для отдела ;

Ограничения

1. Тестовый сценарий отдела можно выполнить на подсети, если все оборудование подсети закреплено за отделом

2. Каждая единица оборудования может быть закреплена не более чем за одним отделом

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

4. Условие покрытия объема тестов для каждого отдела. Ограничение сформулировано мягким с возможностью нарушить (слаковая переменная )

5. Условие справедливого дележа согласно объему тестов

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

• минимизация взвешенного объема нарушений покрытия тестовых сценариев:

Python реализация

Данные

Обременение политикой конфиденциальности не позволяет 3D погружение в задачу. Поэтому воспользуюсь AI (случайным генератором) для дополнения картины недостающими данными.

Входной набор данных состоит из одной таблицы с четырьмя полями:

• dep_id - идентификатор отдела/подразделения;

• tc_id - идентификатор тестового сценария;

• sn_id - идентификатор подсети, где может быть выполнен тестовый сценарий;

• device_id - идентификатор оборудования тестового стенда.

Таблица с данными минималистична, но содержит достаточно информации для принятия решения по разбиению тестового стенда. Ниже приведу сводную информацию по входным данным: общее кол-во тестовых сценариев у департамента, суммарное кол-во подсетей для запуска и кол-во задействованных единиц оборудования тестового стенда в подсетях (тестовый стенд состоит из 45 единиц оборудования).

В модели фигурирует параметр кол-ва тестовых сценариев, которые требуется покрыть у каждого отдела . Поставим это значение равным 60% от общего кол-ва тестов отдела.

Разные отделы имеют разный объем тестовых сценариев для запуска. Как следствие, нагрузка на подсети отличается. Распределение в зависимости от нагрузки обеспечивает ограничение (5), а параметр указывает минимальное кол-во оборудования для отдела. Значение параметров тоже потребуется сгенерировать.

Сформируем бюджет суммы значений : вычтем из общего кол-ва оборудования 7 единиц (сдвиг). Сдвиг делаем для дополнительной вариативности задачи. Полученные 38 единиц оборудования распределим пропорционально кол-ву тестов у отдела. Введем дополнительное ограничение снизу: каждому отделу должно достаться хотя бы 2 единицы оборудования. Значения приведены в таблице выше.

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

Python код

Синтаксис описания моделей для CP и MIP в ORtools отличается, но схема задания модели похожа. Ниже приведу код решения задачи для cp-sat солвера.

Параметр is_agg подробнее разберу в разделе Численный эксперимент. Для случая is_agg=True, постановка задачи полностью соответствует математическому описанию выше.

# Загрузка и обработка входных данных
import pandas as pd
import numpy as np

df_sns = pd.read_csv("https://raw.githubusercontent.com/Lozkins/ORP/master/data/input/01_subnetworks.csv", sep=";", encoding="cp1251")
coverage_rate = 0.60  # Процент покрытия тестов отдела
device_cnt_shift = 7  # Сдвиг общего кол-ва оборудования для генерации нижней границы кол-ва оборудования для отдела
device_cnt = df_sns["device_id"].nunique()

# Генерируем значения для справедливого дележа (минимальное кол-во оборудования для отдела)
dct_tc_parts = (
    df_sns
    .pipe(lambda _df: _df.groupby(["dep_id"], as_index=False)["tc_id"].nunique())
    .assign(tc_part=lambda _df: _df["tc_id"] / _df["tc_id"].sum())
    .assign(device_min=lambda _df: round(_df["tc_part"] * (device_cnt - device_cnt_shift)).astype(int))
    .assign(device_min=lambda _df: np.where(_df["device_min"] < 2, 2, _df["device_min"]))
    .set_index("dep_id")
    .to_dict()["device_min"]
)

import pandas as pd
import numpy as np
from ortools.sat.python import cp_model


def cp_sat(time_limit, is_agg=True):
  """
    Решение задачи разбиения сети на подсети с помощью cp-sat solver
  """

  m = cp_model.CpModel()

  # Инициализация переменных
  # Переменные привязки оборудования к отделу
  df_x = df_sns[["dep_id", "device_id"]].drop_duplicates()
  df_x["var_id"] = np.ogrid[:df_x.shape[0]]
  df_x["var"] = df_x["var_id"].apply(lambda x: m.NewBoolVar(f"x_{x}"))

  # Переменная возможности выполнить тест
  df_y = df_sns[["dep_id", "tc_id"]].drop_duplicates()
  df_y["var_id"] = np.ogrid[:df_y.shape[0]]
  df_y["var"] = df_y["var_id"].apply(lambda y: m.NewBoolVar(f"y_{y}"))

  # Переменная возможности выполнить тест на определенной подсети
  df_z = df_sns[["dep_id", "tc_id", "sn_id"]].drop_duplicates()
  df_z["var_id"] = np.ogrid[:df_z.shape[0]]
  df_z["var"] = df_z["var_id"].apply(lambda s: m.NewBoolVar(f"s_{s}"))

  # Слаковая переменная - объем нарушения уровня покрытия тестов для отдела
  df_p = df_sns.groupby("dep_id", as_index=False)["tc_id"].nunique()
  df_p = df_p.rename({"tc_id": "w"}, axis=1)  # Добавляем веса для целевой
  df_p["var_id"] = np.ogrid[:df_p.shape[0]]
  df_p["var"] = df_p.apply(lambda p: m.NewIntVar(0, round(p.w * coverage_rate), f"p_{p.var_id}"), axis=1)

  # Добавляем ограничения в модель
  # Ограничение: Тестовый сценарий отдела можно выполнить на подсети, если все оборудование подсети закреплено за отделом
  dct_x = df_x.set_index(["dep_id", "device_id"]).to_dict()["var"]
  # Добавим оборудование в таблицу
  df_tmp = df_z.merge(df_sns[["dep_id", "tc_id", "sn_id", "device_id"]], how="left",
                      on=["dep_id", "tc_id", "sn_id"])
  df_cnstr_1 = df_tmp.groupby(["dep_id", "device_id"], as_index=False).agg({"var": "sum", "sn_id": "count"})
  if is_agg:
    for cnstr in df_cnstr_1.itertuples():
      key = cnstr.dep_id, cnstr.device_id
      m.Add(cnstr.var <= cnstr.sn_id * dct_x.get(key, 0))
  else:
    for cnstr in df_tmp.itertuples():
      key = cnstr.dep_id, cnstr.device_id
      m.Add(cnstr.var <= dct_x.get(key, 0))

  # Ограничение: Каждая единица оборудования может быть закреплена не более чем за одним отделом
  df_cnstr_2 = df_x.groupby("device_id").agg({"var": "sum"})
  for cnstr in df_cnstr_2.itertuples():
    m.Add(cnstr.var <= 1)

  # Ограничение: Тестовый сценарий отдела можно выполнить, если хотя бы одна подсеть использует оборудование, полностью закрепленное за этим отделом
  dct_y = df_y.set_index(["dep_id", "tc_id"]).to_dict()["var"]
  df_cnstr_3 = df_z.groupby(["dep_id", "tc_id"], as_index=False).agg({"var": "sum"})
  for cnstr in df_cnstr_3.itertuples():
    key = cnstr.dep_id, cnstr.tc_id
    m.Add(cnstr.var >= dct_y.get(key, 0))

  # Ограничение: Условие покрытия объема тестов для каждого отдела
  dct_p = df_p.set_index("dep_id").to_dict()["var"]
  df_cnstr_4 = df_y.groupby("dep_id", as_index=False).agg({"var": sum, "tc_id": "count"})
  for cnstr in df_cnstr_4.itertuples():
    m.Add(cnstr.var + dct_p.get(cnstr.dep_id, 0) >= round(cnstr.tc_id * coverage_rate))

  # Ограничение: Условие справедливого дележа согласно объему тестов
  df_constr_5 = df_x.groupby("dep_id", as_index=False).agg({"var": "sum"})
  for cnstr in df_constr_5.itertuples():
    m.Add(cnstr.var >= dct_tc_parts.get(cnstr.dep_id, 0))

  # Целевая функция: минимизация взвешенного объема нарушений покрытия тестовых сценариев
  m.Minimize(sum(df_p["var"]))

  # Инициализация solver
  solver = cp_model.CpSolver()

  # Устанавливаем ограничения солвера
  solver.parameters.max_time_in_seconds = time_limit
  solver.parameters.log_search_progress = True

  # Решение задачи
  status = solver.Solve(m)

  # Проверяем статус
  print("Найдено оптимальное решение: ", status == cp_model.OPTIMAL)
  print('Целевая функция = %f' % solver.ObjectiveValue())

  # Извлекаем решение
  df_x["sol"] = df_x["var"].apply(lambda x: solver.value(x))
  df_y["sol"] = df_y["var"].apply(lambda x: solver.value(x))
  df_z["sol"] = df_z["var"].apply(lambda x: solver.value(x))
  df_p["sol"] = df_p["var"].apply(lambda x: solver.value(x))

  return dict(df_x=df_x, df_y=df_y, df_z=df_z, df_p=df_p)

Численный эксперимент

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

• Реализация в парадигме CP-SAT;

• Реализация в парадигме MIP.

Набор солверов следующий: SCIP, коммерческий солвер (MIP) и cp-sat (CP). Два солвера open source и один коммерческий.

Размерность исходной задачи: 985 ограничений и 19702 бинарные переменные и 5 целочисленных переменных (случай is_agg=True).

Альтернативная модель

Построение модели - процесс творческий. Математическая формулировка задачи влияет на скорость поиска решения. В рамках эксперимента продемонстрирую влияние этого фактора на производительность солверов.

Рассмотрим ограничения типа (1): Тестовый сценарий отдела можно выполнить на подсети, если все оборудование подсети закреплено за отделом.

Левая часть представляет собой сумму допустимых подсетей отдела. Если хотя бы одно оборудование не находится в юрисдикции отдела, то эта сумма должна равняться 0. Предлагаю перейти к следующему набору ограничений:

Логика ограничения не изменилась, но кол-во ограничений увеличится, т.к. для каждой допустимой подсети будет создан набор ограничений равный кол-ву единиц оборудования в топологии тестового сценария.

При замене ограничений (1) на новые получим задачу следующей размерности: 101258 ограничений и 19702 бинарных переменных и 5 целочисленных переменных (случай is_agg=False). Кол-во ограничений в модели выросло в ~106 раз. Целесообразность такого перехода рассудит эксперимент.

Производительность моделей и солверов

Рассмотрел три солвера: cp, scip и com; каждый решал две задачи: agg (исходное ограничение (1)) и unzip (распакованное ограничение (1)). На графике выведены значения целевой функции (по оси y) из логов солверов с течением времени (ось x). Все солверы запускал с базовыми настройками. Значение целевой функции - размер нарушения покрытия, выраженное в кол-ве тестовых сценариев, которое требуется для удовлетворения требуемого покрытия.

Разберем аутсайдера эксперимента: scip. Первое решение находит спустя 1959 сек. и 3261 сек. для agg и unzip моделей, соответственно. Низкую производительность солвера можно связать с однопоточностью солвера на этапе решения MIP задачи.

В логах солвера SCIP при решении задачи agg можно обнаружить следующую запись:

10 deleted vars, 235 deleted constraints, 100498 added constraints, 5 tightened bounds, 0 added holes, 0 changed sides, 4 changed coefficients

Солвер scip добавил в модель 100498 ограничений на этапе presolve. Это наводит на мысль, что солвер сам распаковал ограничения типа (1). Если знакомы с внутренним устройством солвера scip, прошу прокомментировать догадку.

Коммерческий солвер сразу нашел допустимое решение задачи agg, для задачи unzip спустя 34 сек. Но решение в 265 непокрытых тестовых сценариев нашел быстрее в случае unzip (3483 сек.) против 3600 сек. для agg модели. Сходимость задачи agg лучше в первой половине в то время, как unzip сходит немного быстрее во второй половине отведенного времени.

Программирование в ограничениях оказалось наиболее подходящим для решения задачи даже с учетом использования open source солвера cp-sat. Помимо лучшего значения целевой функции (262 шт.) среди рассмотренных солверов, в случае unzip решение найдено через 479 сек. Поздравляем с победой в локальном benchmark'е!

Кол-во произведенных симуляций не велико, заключение не то чтобы объективно. Однако, формулировка с распакованным ограничением типа (1) - unzip, показывает лучшую сходимость и является более "удобной" для солверов (для данной задачи).

Провел две симуляции эксперимента (2х6=12шт. расчетов), поведение солверов имеет сходство. Ниже аналогичный график для второй симуляции.

Коммерческий и scip солверы работают достаточно стабильно, в то время как cp-sat переобулся и показал обратную динамику. Во втором прогоне лучшее решение 264 непокрытых сценариев на версии agg. Нестабильность - цена доступности.

Дополнительно, хочу отметить следующее: все солверы не смогли улучшить нижнюю границу решения и отклониться от уровня 0.

Результат

Разберем решения первой симуляции. Доказанного оптимального решения получить не удалось в установленные лимиты времени. Это не было камнем преткновения для использования решения в промышленной эксплуатации, заказчик доволен и пользуется моделью.

Эксперименты с ограничением по времени в 12 часов тоже не привели к доказанному оптимальному решению. Но и улучшить значение целевой функции тоже не удалось. Выглядит, что задача доставляет трудности алгоритмам оценки нижней границы решения.

В таблице представлено кол-во не покрытых тестовых сценариев в разрезе разделов. Столбец cover_req требуемое кол-во тестов для покрытия. Больше всего пострадал первый отдел, у него ни один тест не покрыт. Лучшее покрытие у самого большого отдела (dep_id = 5). Настройка весов позволяет сместить приоритет с самого "крупного" отдела, кроме того более умная настройка так же позволит повысить вариабельность задачи (сделать дележ более демократичным).

Распределение кол-ва единиц оборудования между отделами одинаково во всех решениях {1: 2, 2: 5, 3: 8, 4: 10, 5: 20} (ограничение на минимальный объем привязки : {1: 2, 2: 5, 3: 8, 4: 10, 5: 13}). Т.е. все свободное оборудование (сдвиг в 7 единиц) назначено самому "крупному" отделу.

Заключение

Разработали решение разбиения тестового стенда на непересекающиеся подсети. Инструмент условно не ограничен кол-вом разбиений и позволяет получить приемлемое решение в пределах часа. Имеет пару управляемых параметров и для настройки модели. Предоставляет возможность моделировать различные гипотезы с меньшими затратами в сравнении с ручным моделированием.

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

Что пробовали:

• У cp солверов есть условные ограничения типа EnforceOnly, которые можно применить для модели из статьи. В моих экспериментах значимого ускорения от их использования не произошло.

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

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

Пути по ускорению поиска решения. Можно заниматься настройкой солверов под задачу, которые позволят ускорить поиск решений (например, в статье описан auto-tuner для MIP, или здесь). Альтернатива погружению в детали работы солвера - варьирование или декомпозиция модели на базе экспертных предположений. Мы пошли по второму пути и преуспели в этом.

Ссылки

• Jupyter блокнот статьи

• Задача о разбиении множества