Хабрахабр

[Из песочницы] Логистика акции по раздельному сбору вторсырья

Вместо вступления

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

К этому же времени туда подъезжает арендованный грузовик, который отвозит собранное и рассортированное вторсырье на площадку, откуда оно перераспределяется между различными перерабатывающими предприятиями. В Новосибирске такая деятельность формируется вокруг акции «Зеленая белка», в рамках которой раз в месяц обеспокоенные экологией горожане приносят накопленные перерабатываемые бытовые отходы в заранее определенные места в известное время. Для маршрутизации грузовиков одного лишь пристального взгляда стало не хватать, и мы начали разрабатывать оптимизационные модели для минимизации транспортных затрат. Акция существует с 2014 года, и с того времени значительно увеличилось число точек сбора вторсырья, а также его объемы. Первой из таких моделей и посвящена данная статья.

Далее, в разделе 2, задача минимизации транспортных затрат будет формализована в виде задачи маршрутизации разнородных транспортных средств с временными окнами (heterogenious fleet vehicle routing problem with time windows). В разделе 1 я подробно и с иллюстрациями опишу схему организации акции. Раздел 3 посвящен решению данной задачи с использованием свободно распространяемого пакета для решения смешанно-целочисленных линейных задач математического программирования GLPK.

1. Акция «Зеленая белка»

Начиная с 2014 года инициативная группа «Живая земля» проводит ежемесячную акцию по раздельному сбору вторсырья «Зеленая белка» в Новосибирске. На момент написания статьи к сдаче в переработку с рядом оговорок принимаются пластиковые отходы с маркировкой PET, PE, PP, PS, стекло, алюминий, железо, бытовая техника, макулатура и другое. В подготовке и проведении акции участвует более 50 волонтеров. Акция не является коммерческой, участие в ней бесплатное и не подразумевает денежного вознаграждения.

Точки

На фото одна из таких точек: слева вдоль стены стоят мешки для сбора различных фракций пластика, справа — грузовик, в который в дальнейшем все грузится, а в центре — волонтер с ушками. Акция проводится в 17 точках города, расположенных друг от друга на расстояниях, преодолеваемых автомобилем за время от 15 до 90 минут.

image

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

Также имеются данные по средним объемам вторсырья, собираемым на каждой из точек.

Маршруты

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

image

Грузовики

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

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

2. Формализация

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

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

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

  • формирование групп точек, составляющих отдельный маршрут;
  • определение схемы объезда каждой из групп;
  • удовлетворение требований на время прибытия грузовика на каждую из точек;
  • определение типа грузовика, необходимого для обслуживания каждого из маршрутов.

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

Формирование групп точек

Положим $inline$\bar{V} = V\cup \{0\}$inline$ Пусть $inline$V = \$inline$ есть множество индексов точек сбора вторсырья, и площадка, куда собранное вторсырья затем отвозится — будем называть ее депо — имеет индекс 0.

Через $inline$R$inline$ обозначим множество номеров маршрутов. Каждая группа точек формирует отдельный маршрут.

Тогда требование о том, что точка $inline$i\in V$inline$ должна войти в один из маршрутов можно записать равенством
Пусть величина $inline$z_{ir}, i\in I, r\in R$inline$ равняется единице, если точка $inline$i$inline$ включается в маршрут с номером $inline$r$inline$, и нулю в противном случае.

$$display$$\sum_{r\in R}z_{ir} = z_{i1} + z_{i2} + \dots + z_{in} = 1.$$display$$

Депо должно войти во все маршруты: $inline$z_{0r} = 1, r\in R$inline$

Тонкости

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

$$display$$1-z_{ir} \leq \sum_{j = 1}^{i-1}\left(1-\sum_{k = 1}^{r-1}z_{jk}\right) + \sum_{k = 1}^{r-1} z_{ik}, \quad i\in I, r\in R, r\leq i.$$display$$

Определение схемы объезда

Формально все они начинаются в одной из точек множества $inline$V$inline$ и заканчиваются в депо, однако проще будет считать, что все маршруты циклические. Маршруты представляют из себя чередующуюся последовательность точек и переездов между ними. Это предположение не меняет сути, но делает все вершины одинаковыми: в этом случае нет концевых вершин, все они транзитные.

Для точек $inline$i, j\in \bar{V}$inline$ и маршрута $inline$r\in R$inline$ заведем переменную $inline$x_{ijr}$inline$, равную единице, если в маршруте c номером $inline$r$inline$ грузовик совершает переезд из точки $inline$i$inline$ в точку $inline$j$inline$, и нулю иначе.

Тогда потребуем, во-первых, чтобы грузовик, двигающийся по маршруту $inline$r\in R$inline$ посещал точку $inline$i\in V$inline$, если она после разделения на группы попала в группу с номером $inline$r$inline$:

$$display$$\sum_{j\in \bar{V}}x_{jir} = z_{ir}, \quad i\in \bar{V}, r \in R.$$display$$

Во-вторых, грузовик после приезда на точку $inline$i\in \bar{V}$inline$ должен с нее выехать.

$$display$$\sum_{j\in \bar{V}}x_{jir} = \sum_{j\in \bar{V}}x_{ijr}, \quad i \in \bar{V}, r\in R.$$display$$

Маршруты такого вида, безусловно, не имеют смысла, и требуется добавить ряд ограничений, чтобы этого избежать. Можно заметить, что данные ограничения позволяют величинам $inline$x_{ijr}$inline$ задавать маршруты, представляющие собой набор непересекающихся циклов.

Об устранении подциклов

Одним из способов может быть введение вспомогательных неотрицательных величин $inline$f_{ijr}, i, j\in \bar{V}, r\in R$inline$ Записанная с использованием данных величин совокупность соотношений устраняет комбинации значений $inline$x_{ijr}$inline$, задающие маршруты, которые не являются связным циклом.

$$display$$\sum_{i\in V}f_{0jr} = \sum_{i\in V}z_{ir}, \quad r \in R,$$display$$

$$display$$f_{ijr}\leq nx_{ijr}, \quad i\in \bar{V}, j \in \bar{V}, r\in R,$$display$$

$$display$$\sum_{j\in \bar{V}}f_{jir} = \sum_{j\in \bar{V}}f_{ijr} + z_{ir}, \quad i\in V, r \in R.$$display$$

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

Удовлетворение требований на время прибытия грузовика на каждую из точек

Через $inline$B_i$inline$ и $inline$E_i$inline$ обозначим, соответственно, начало и конец временного интервала, в котором куратор точки $inline$i$inline$ может на ней присутствовать. Другими словами, необходимо посещать точки только внутри временных окон, указанных кураторами.

Через $inline$L_i, i\in V$inline$ обозначим длительность погрузки на точке $inline$i$inline$, а через $inline$D_{ij}, i, j\in \bar{V}$inline$ — время переезда из точки $inline$i$inline$ в точку $inline$j$inline$ Оговоримся, что $inline$D_{v_0i} = 0$inline$ для всех $inline$i\in \bar{V}$inline$, и в целом может выполняться $inline$D_{ij} \neq D_{ji}$inline$ для любых $inline$i\neq j$inline$ Для отслеживания выполнения данных ограничений, нам потребуется информация о времени, затрачиваемом грузовиком во время стоянок и переездов на маршруте.

Для введенных величин должны выполняться следующие соотношения. Введем переменные неотрицательные величины $inline$a_i, i\in \bar{V}$inline$ и $inline$w_{ir}, i\in \bar{V}, r\in R$inline$, равные временам прибытия и ожидания на точке $inline$i$inline$ при движении по маршруту $inline$r$inline$, соответственно.

Время ожидания не может быть меньше времени, требуемого для погрузки

$$display$$w_{ir}\geq L_iz_{ir}, \quad i\in \bar{V}, r\in R,$$display$$

и равно нулю, если точка не относится к маршруту $inline$r$inline$

$$display$$w_{ir}\leq Mz_{ir}, \quad i\in \bar{V}, r\in R,$$display$$

Время прибытия на точку $inline$j$inline$ вычисляется через соответствующие времена для предшествующей точки $inline$i$inline$ Здесь достаточно большая константа $inline$M$inline$ используется для устранения ненужных зависимостей в случае, когда переезд между $inline$i$inline$ и $inline$j$inline$ не совершается.

$$display$$a_i + \sum_{r\in R}w_{ir} + D_{ij}\leq a_j + M(1 - \sum_{r\in R}x_{ijr}), \quad i\in I, j\in V,$$display$$

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

$$display$$a_i\geq B_i, \quad i\in V,$$display$$

$$display$$a_i + \sum_{r\in R}w_{ir}\leq E_i, \quad i\in V.$$display$$

Определение типа грузовика, необходимого для обслуживания каждого из маршрутов
Обозначим множество типов грузовиков, доступных для аренды, через $inline$T$inline$ Грузовик типа $inline$t\in T$inline$ характеризуется объемом кузова $inline$C_t$inline$ и стоимостью часовой аренды $inline$P_t$inline$ Минимальное время аренды грузовика типа $inline$t$inline$ обозначим через $inline$U_t^0$inline$ Объем вторсырья, собираемый на точке $inline$i\in V$inline$ обозначим через $inline$G_i$inline$

Введем переменные величины $inline$y_{tr}, t\in T, r\in R$inline$, равные единице, если грузовик типа $inline$t$inline$ назначается на обслуживание маршрута с номером $inline$r$inline$, и нулю иначе.

Целочисленные переменные $inline$u_{tr}, t\in T, r\in R$inline$ задают время аренды грузовика типа $inline$t$inline$, обслуживающего маршрут с номером $inline$r$inline$
Для каждого маршрута, $inline$r\in R$inline$, необходимо назначить тип грузовика:

$$display$$\sum_{t\in T}y_{tr} = 1, \quad r\in R$$display$$

В соответствии с разбиением точек между маршрутами, некоторые маршруты могут оказаться тривиальными, то есть содержать только депо. Если маршрут $inline$r$inline$ нетривиальный, то грузовик, назначенный на него, арендуется не менее чем на $inline$U_t^0$inline$ часов.

$$display$$u_{tr}\geq U_t^0 (y_{tr} - \sum_{i\in V}z_{ir}), \quad t\in T, r\in R.$$display$$

В то же время, длительность аренды должна превосходить также суммарную длительность стоянок и переездов в маршруте.

$$display$$u_{tr}\geq \sum_{i\in V}\sum_{j\in \bar{V}}D_{ij}x_{ijr}+\sum_{i\in \bar{V}}w_{ir}-M(1-y_{tr}), \quad r\in R, t\in T.$$display$$

Добавим ограничения обеспечивающие свойство: если грузовик типа $inline$t$inline$ не назначен на маршрут $inline$r$inline$, время его аренды равно нулю.

$$display$$u_{tr}\leq My_{tr}.$$display$$

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

$$display$$\sum_{i\in V}G_iz_{ir}\leq \sum_{t\in T}C_ty_{tr}, \quad r\in R.$$display$$

Наконец, наша цель состоит в минимизации затрат на аренду грузовиков, которая с помощью введенных обозначений записывается как $inline$\sum_{t\in T}\sum_{r\in R}P_tu_{tr}$inline$

Поиск решения

Несложно убедиться в том, что все выражения, участвующие в оптимизационной модели, являются линейными функциями от переменных величин, поэтому для нахождения точного и приближенных решений можно пользоваться стандартными пакетами для решения задач смешанно-целочисленного программирования таких, как Gurobi, CPLEX, Xpress, ORtools, SCIP, BLIS и т.д.
Запишем модель минимизации транспортных затрат на языке GMPL. Это позволит использовать для наших целей свободно распространяемый пакет GLPK. Для написания кода и отладки модели удобно будет скачать среду разработки GUSEK, которая уже содержит GLPK в своем составе.

Выглядит GUSEK следующим образом:

image

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

Полное описание модели

param numOfPoints > 0, integer; #число точек
param numOfTypes > 0, integer; #число типов грузовиков
param numOfRoutes = numOfPoints;#максимальное число маршрутов set V := 1 .. numOfPoints; #множество точек
set Vbar := V union {0}; #множество точек с пунктом разгрузки (депо)
set T := 1 .. numOfTypes; #множество типов грузовиков
set R := 1 .. numOfPoints; #множество маршрутов #####Точки##########
param WDL >= 0, default 8; #длительность рабочего дня
param B{i in Vbar} >= 0; #начало временного окна
param E{i in Vbar} >= 0; #конец временного окна
param L{i in Vbar} >= 0; #минимальное время погрузки
param D{i in Vbar, j in Vbar} >= 0, <= WDL; #время переезда
param G{i in V}, >= 0; #объем вторсырья, м3
#####Грузовики#####
param C{t in T}, >= 0; #вместительность кузова
param P{t in T}, >= 0; #стоимость аренды за час
param U0{t in T}, >= 0; #минимальное время аренды, часов /**********************************************
* Формирование групп точек
**********************************************/
var z{Vbar, R} binary; #равняется единице, если точка включается в маршрут, и нулю в противном случае s.t. pointToGroup 'point to group' {i in V}: sum{r in R} z[i, r] == 1;
s.t. depotToGroup 'depot to group' {r in R}: z[0, r] == 1;
s.t. lexMinGroups 'lexicographycally minimal division' {i in V, r in R: r <= i}: 1 - z[i, r] <= sum{j in V: j <= i - 1}(1 - sum{k in R: k <= r - 1} z[j, k]) + sum{k in R: k <= r - 1}z[i, k] ; /**********************************************
* Определение схемы объезда
**********************************************/
var x{Vbar, Vbar, R} binary; #равна единице, если в маршруте c номером r грузовик совершает переезд из точки i в точку j, и нулю иначе. s.t. visitPoint 'visit point' {i in Vbar, r in R}: sum{j in Vbar} x[i, j, r] = z[i, r];
s.t. keepMoving 'keep moving' {i in Vbar, r in R}: sum{j in Vbar} x[j, i, r] = sum {j in Vbar} x[i, j, r]; var f{Vbar, Vbar, R} >= 0; #Потоки, устраняющие подциклы. s.t. flowFromDepot 'flow from depot' {r in R}: sum{i in V} f[0, i, r] == sum{i in V} z[i, r];
s.t. flowAlongActiveArcs 'flow along active arcs' {i in Vbar, j in Vbar, r in R}: f[i, j, r] <= numOfPoints * x[i, j, r];
s.t. flowConservation 'flow conservation' {i in V, r in R}: sum{j in Vbar} f[j, i, r] == sum{j in Vbar} f[i, j, r] + z[i, r]; var a{i in V} >= 0; #время прибытия грузовика на точку
var w{i in Vbar, r in R} >= 0; #время нахождения грузовика на точке s.t. wait 'wait'{i in Vbar, r in R}: w[i, r] >= L[i] * z[i, r];
s.t. dontWait 'dont wait'{i in Vbar, r in R}: w[i, r] <= (E[i] - B[i]) * z[i, r];
s.t. arrivalTime 'arrival time' {i in V, j in V}: a[i] + sum{r in R}w[i, r] + D[i,j] <= a[j] + 3 * WDL * (1 - sum{r in R} x[i, j, r]);
s.t. arriveAfter 'arrive after' {i in V}: a[i] >= B[i];
s.t. departBefore 'depart before' {i in V}: a[i] + sum{r in R}w[i, r] <= E[i]; /**********************************************
* Определение типа грузовика, необходимого для обслуживания каждого из маршрутов
**********************************************/ var y{t in T, r in R}, binary; #равные единице, если грузовик типа t назначается на обслуживание маршрута с номером r, и нулю иначе.
var u{t in T, r in R}, integer, >= 0; #время аренды грузовика типа t, обслуживающего маршрут с номером r. s.t. assignVehicle 'assign vehicle' {r in R}: sum{t in T} y[t,r] == 1;
s.t. rentTime 'rent time' {r in R, t in T}: u[t, r] >= sum{i in V, j in Vbar}D[i, j] * x[i, j, r] + sum{i in Vbar}w[i, r] - WDL * (1 - y[t, r]);
s.t. minRentTime 'minimal rent time' {r in R, t in T}: u[t, r] >= U0[t] * (y[t, r] - sum{i in V}z[i, r]);
s.t. noRent 'no rent' {t in T, r in R}: u[t, r] <= WDL * y[t, r];
s.t. fitCapacity 'fit capacity' {r in R}: sum{i in V} G[i] * z[i, r] <= sum{t in T} C[t] * y[t, r]; minimize rentCost: sum{t in T, r in R} P[t] * u[t, r]; solve; /**********************************************
* Вывод решения на экран
**********************************************/
printf{i in V, r in R} (if 0.1 < z[i,r] then "point %s to group %s\n" else ""), i, r, z[i,r];
printf{r in R, i in Vbar, j in Vbar} (if 0.1 < x[i, j, r] then "route %s: %s -> %s\n" else ""), r, i, j;
printf{i in V} "point %s arrive between %s and %s (actual = %s)\n", i, B[i], E[i], a[i];
end;

Для быстрого старта добавлю также подготовленные для использования в модели данные, взятые из головы:

Входные данные

data; param numOfPoints := 9; #число точек
param numOfTypes := 6; #число типов грузовиков #####Точки##########
param : B E L := 0 0 8 1 1 0 8 1 2 0 8 1 3 0 8 1 4 0 8 1 5 0 8 1 6 0 8 1 7 0 8 1 8 0 8 1 9 0 8 1;
param D default 0
: 0 1 2 3 4 5 6 7 8 9 :=
0 . . . . . . . . . .
1 0.1 0.3 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1
2 0.3 0.2 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1
3 0.4 0.3 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1
4 0.4 0.4 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1
5 0.1 0.2 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1
6 0.5 0.5 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1
7 0.3 0.2 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1
8 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1
9 0.5 0.2 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1;
param G := 1 1 2 2 3 1 4 2 5 1 6 2 7 1 8 2 9 1;
#####Грузовики#####
param : C P := 1 8 500 2 10 500 3 14 600 4 18 650 5 25 700 6 35 800;
param U0 default 2; #минимальное время аренды, часов end;

После копирования кода модели в файлик с именем, к примеру, model.mod, а входные данные — в файлик data.dat, все готово к запуску. Установим ограничение на время вычисления 100 секунд (это делается с помощью ключа --tmlim [время в секундах]), передадим путь к файлу с входными данными (с помощью ключа -d [путь к файлу]),

image

В случае успеха в окне справа начнут появляться сообщения, и через сто секунд у нас на руках будет лучшее решение, которое GLPK успел найти за отведенное время. и жмем F5.

image

Как видно, оно время от времени уменьшается. В синем тексте нас интересует значение после надписи «mip = ». Следующее число — это оценка снизу для наилучшего существующего, т. В процессе работы находятся все новые решения, значение транспортных затрат на лучшем из которых выводится в этом столбце (пока там 14700). оптимального решения. е. Значения слева и справа сближаются, и относительный разрыв между ними на момент скриншота составляет 54. Изначально оценка ощутимо занижена, но со временем уточняется, то есть возрастает. Как только это число станет равным нулю, алгоритм докажет, что лучшее найденное решение является наилучшим возможным. 1%. Дожидаться этого события на практике не всегда оправданно, и не только потому, что это долго, однако важно оговориться, что как правило, хорошее решение находится относительно быстро, и основные временные затраты связаны с уточнением оценки, требуемым для доказательства неулучшаемости.

Вместо заключения

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

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

Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть