Проблема распределения объектов модели смешанного целочисленного программирования

Проблема распределения объектов МСЦП

Подробное руководство по Python для моделей p-рассеивания и максисума

Фото от Z на Unsplash

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

Kuby (1987) предложил две различные смешанные целочисленные программные формулировки для этих задач: задачу p-рассеивания и задачу максисума. В этой статье оба будут реализованы с использованием библиотеки Python Pyomo и решателя HiGHS.

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

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

Как обычно, полный код можно найти в этом репозитории git.

Выберите одну из этих локаций для размещения 5 объектов:

Возможные местоположения в задаче размещения объектов. (Изображение от автора).

Произведение бинарных переменных

При определении основных элементов этой задачи можно использовать бинарные переменные, которые принимают значение 1, если местоположение выбрано, и 0 в противном случае. Обозначим эти переменные xᵢ. Предположим, что расстояния (или некоторая другая метрика рассеивания) между двумя местоположениями вычисляются заранее, и обозначим их как dᵢⱼ. Как мы можем вычислить рассеивание для любой выбранной пары объектов?

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

Рассмотрим ориентированный граф G(V, A) и zᵢⱼ новую бинарную переменную, которая указывает, что и узлы i, и j выбраны. Если i или j не выбраны, zᵢⱼ должна быть равна 0. Это приводит к первой группе ограничений:

Первая группа ограничений в линеаризованной форме произведения бинарных переменных. (Изображение от автора).

До сих пор, даже если и i, и j выбраны, zᵢⱼ может принимать значение 0. Поэтому мы должны добавить дополнительное ограничение так, чтобы при выборе i и j, zᵢⱼ становилась равной 1.

Вторая группа ограничений в линеаризованной форме произведения бинарных переменных. (Изображение от автора).

При максимизации чего-то, пропорционального zᵢⱼ, как в задаче максисума, вторая группа ограничений не нужна, так как не имеет смысла не вычислять прирост, пропорциональный zᵢⱼ, если это возможно. Однако это может быть полезно при формулировке других MIP-моделей.

В следующем разделе давайте применим эти концепции к задаче максимума.

Модель максимума

Дискретная задача максимума должна определить расположение p объектов среди N дискретных узлов для максимизации суммы расстояний (или среднего расстояния), вычисленных между всеми парами выбранных объектов. Поэтому предположим, что объекты – это узлы, расположенные в направленном графе G(V, A). Вес каждой дуги от i до j – это известная заранее метрика расстояния (разброса) dᵢⱼ. Также рассмотрим бинарные переменные xᵢ, чтобы указать, выбрано ли местоположение i, и zᵢⱼ, чтобы указать, выбраны ли оба i и j.

Цель может быть сформулирована следующим образом:

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

Таким образом, мы можем сформулировать ограничения задачи следующим образом:

Давайте преобразуем это в код с использованием Python. Для этого мы должны начать с импорта используемых библиотек. Библиотека numpy будет использоваться для вычислений линейной алгебры и создания случайных координатных точек; функции squareform и pdist из scipy будут использоваться для вычисления расстояний на основе матрицы координат; matplotlib будет нашим инструментом визуализации; а pyomo будет использоваться в качестве языка алгебраического моделирования (AML) (с решателем HiGHS).

import numpy as npfrom scipy.spatial.distance import squareform, pdistimport matplotlib.pyplot as pltimport pyomo.environ as pyofrom pyomo.contrib.appsi.solvers.highs import Highs

Мы должны создать N точек и определить, сколько из них должны быть выбраны в качестве местоположений объектов. Зафиксировав случайное зерно (12) при каждом выполнении кода, должны быть получены точки, представленные во введении.

# Фиксация случайного зернаnp.random.seed(12)# Создание случайных точекN = 25p = 5coordinates = np.random.random((N, 2))# Вычисление попарных расстоянийdistances = squareform(pdist(coordinates))

Теперь у нас есть необходимые элементы для начала нашей модели pyomo.

Есть два подхода к моделированию проблемы в pyomo: абстрактная и конкретная модели. В первом подходе алгебраические выражения проблемы определяются до предоставления некоторых значений данных, тогда как во втором экземпляр модели создается немедленно при определении ее элементов. Подробнее об этих подходах можно узнать в документации библиотеки или в книге Байнума и др. (2021). В этой статье мы будем использовать формулировку конкретной модели.

model = pyo.ConcreteModel()

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

# Наборы узлов и дугmodel.V = pyo.Set(initialize=range(N))model.A = pyo.Set(initialize=[(i, j) for i in model.V for j in model.V if i != j])

Предварительно предоставленные параметры – это количество узлов, которые должны быть выбраны, и расстояния между парами узлов.

# Параметрыmodel.d = pyo.Param(model.A, initialize={(i, j): distances[i, j] for (i, j) in model.A})model.p = pyo.Param(initialize=p)

Затем давайте включим переменные принятия решений.

# Переменные принятия решенияmodel.x = pyo.Var(model.V, within=pyo.Binary)model.z = pyo.Var(model.A, within=pyo.Binary)

И ограничения…

# Выбрано p узловdef p_selection(model):    return sum(model.x[:]) == model.p# Если начальный узел не выбран, то дуга равна 0def dispersion_c1(model, i, j):    return model.z[i, j] <= model.x[i]# Если конечный узел не выбран, то дуга равна 0def dispersion_c2(model, i, j):    return model.z[i, j] <= model.x[j]# Включение ограничений в модельmodel.p_selection = pyo.Constraint(rule=p_selection)model.dispersion_c1 = pyo.Constraint(model.A, rule=dispersion_c1)model.dispersion_c2 = pyo.Constraint(model.A, rule=dispersion_c2)

Наконец, нам необходимо создать целевую функцию.

def disp_obj(model):    return sum(model.z[i, j] * model.d[i, j] for (i, j) in model.A)model.obj = pyo.Objective(rule=disp_obj, sense=pyo.maximize)

Теперь модель maxisum готова к решению. Затем мы должны создать совместимый с pyomo решатель для ее обработки. Решатель Highs доступен в Pyomo (проверьте импорты) с версии 6.4.3, если установлен пакет highspy. Убедитесь, что вы запустили pip install highspy.

solver = Highs()solver.solve(model)

И, после примерно 120 секунд вычислительного времени, мы получаем следующие результаты:

Результаты модели maxisum. (Изображение автора).

Обратите внимание, что даже если общая дисперсия максимизируется, два объекта в верхнем левом углу все еще находятся достаточно близко друг к другу, что может быть нежелательным. Таким образом, в качестве альтернативы формулировке maxisum у нас есть модель p-дисперсии, в которой мы максимизируем минимальное расстояние между любой парой выбранных узлов.

Модель p-дисперсии

В модели p-дисперсии нам нужна дополнительная решающая переменная для вычисления минимального расстояния между всеми парами выбранных узлов, которое является нашей целью для максимизации. Обозначим эту переменную D. Мы должны создать большое ограничение М, которое гарантирует, что если и i, и j выбраны, то D не превышает dᵢⱼ; в противном случае, мы должны гарантировать, что D неограничено. Если мы сохраняем формулировку с zᵢⱼ в виде произведения бинарных переменных, мы можем сформулировать это дополнительное ограничение следующим образом.

Максимальное минимальное ограничение с произведением бинарных переменных. (Изображение автора).

В этой формулировке M – это фиксированный параметр, произвольно большой, используемый для формулировки дизъюнктивного правила. Хороший выбор M должен быть достаточно большим, чтобы гарантировать выполнимость ограничения (i, j) для любого значения D, когда zᵢⱼ равно нулю, но настолько мал, насколько это возможно, чтобы линейная релаксация модели была похожа на целочисленную версию, что облегчает сходимость. В этой проблеме наибольшее значение dᵢⱼ может быть интересной альтернативой.

Хотя эта формулировка будет работать хорошо для этой модели, можно использовать более краткий подход, в котором переменные zᵢⱼ даже не включаются в модель.

Максимальное минимальное ограничение с узловыми переменными. (Изображение автора).

В этой формулировке либо xᵢ, либо xⱼ, равное нулю, является достаточным условием для гарантии, что неравенство справедливо для любого значения D. И цель просто сводится к максимизации D.

В следующем фрагменте кода на Python предположим, что у нас есть новая модель с теми же наборами и параметрами, что и предыдущая модель, а также группа решающих переменных x.

# Максимальное минимальное ограничениеdef maxmin_rule(model, i, j):    return model.D <= model.d[i, j] + model.M * (1 - model.x[i]) + model.M * (1 - model.x[j])# Новый параметр большое Mmodel.M = max(model.d[:, :])# Новая переменнаяmodel.D = pyo.Var(within=pyo.NonNegativeReals)# Новое ограничениemodel.maxmin_rule = pyo.Constraint(model.A, rule=maxmin_rule)# Целевая функцияmodel.obj = pyo.Objective(expr=model.D, sense=pyo.maximize)

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

Результаты модели p-дисперсии. (Изображение автора).

Что кажется отличным, так это то, что местоположения хорошо распределены в пространстве.

Есть ли способ улучшить эту дистрибуцию?

Мультикритериальный подход

Помните, что целевая функция модели p-распределения зависит только от минимального расстояния между любой парой выбранных узлов. Поэтому можно получить несколько решений, используя комбинации двух точек, определяющих это расстояние, и других с расстояниями друг к другу, большими или равными самой цели. Можем ли мы улучшить наше решение, выбрав лучшую из этих альтернатив? Это приводит к двухэтапному “мультикритериальному подходу”, предложенному Куби (1987).

На первом этапе решается модель p-распределения для достижения оптимальности, и текущая цель сохраняется в параметре d_opt. Затем решается модель maxisum с дополнительным ограничением, гарантирующим, что D ≥ d_opt, чтобы получить, среди оптимальных решений модели p-распределения, то, которое имеет максимальное общее разбросание.

Чтобы сделать это в Python, предположим, что у вас есть модель p-распределения, созданная также с переменными принятия решений и ограничениями модели maxisum.

# D должно быть оптимальнымdef composed_constr(model):    return model.D >= model.d_opt# Решаем p-распределениеsolver.solve(model)# Новый параметрmodel.d_opt = pyo.Param(initialize=model.obj())# Деактивируем старую целевую функциюmodel.obj.deactivate()# Дополнительные переменныеmodel.z = pyo.Var(model.A, within=pyo.Binary)# Решение не ухудшит текущее Dmodel.composed_constr = pyo.Constraint(rule=composed_constr)# Новая целевая функцияmodel.obj_disp = pyo.Objective(rule=disp_obj, sense=pyo.maximize)solver.solve(model)

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

Мультикритериальная проблема: модель p-распределения, за которой следует максимальная цель maxisum. (Изображение автора).

Дальнейшее чтение

Когда клиенты не равномерно распределены, объекты имеют ограниченную вместимость или, возможно, адекватное количество объектов неизвестно заранее, вероятно, вы сталкиваетесь с другой проблемой размещения объектов. Вы можете найти формулировку Балинского (1965), реализованную на Python с использованием PuLP, в этой удивительной статье Николо Козимо Альбанезе.

Оптимизация: Проблема размещения объектов с ограничением вместимости на Python

Найдите оптимальное количество и местоположение складов для снижения затрат и удовлетворения спроса

towardsdatascience.com

Выводы

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

Ссылки

Balinski, M. L. 1965. Integer programming: methods, uses, computations. Management science, 12(3), 253–313.

Bynum, M. L., Hackebeil, G. A., Hart, W. E., Laird, C. D., Nicholson, B. L., Siirola, J. D., … & Woodruff, D. L. 2021. Pyomo-optimization modeling in python (Vol. 67). Berlin/Heidelberg, Germany: Springer.

Kuby, M. J. 1987. Programming models for facility dispersion: The p‐dispersion and maxisum dispersion problems. Geographical Analysis, 19(4), 315–329.