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

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

Делиться

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

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

Этот пост выполняет данное обещание.

Мы начнём с того места, где закончилась предыдущая статья, покажем, почему очевидный подход терпит крах в больших масштабах, выявим единственное наблюдение, которое нас спасает, а затем медленно разберём сам алгоритм. Математики будет немного больше, чем в прошлый раз, но ничего экзотического; самым продвинутым инструментом, который мы будем использовать, является двойственность линейного программирования, и даже она представлена в виде «мы можем записать ту же задачу по-другому».

Краткое резюме: модель, которую мы пытаемся решить.

В предыдущем посте нашим героем была немецкая компания, занимающаяся производством зимней одежды и закупающая ее в Бангладеш. Спрос ξ — случайная величина; решение о производстве x должно быть принято до того, как станет известно ξ. После того, как спрос станет известен, решение y, предусматривающее устранение дефицита, принимается с большими затратами (например, экстренная партия из Румынии).

Если предположить дискретное распределение с S сценариями — значениями ξ¹, ξ², …, ξS^S с вероятностями p₁, …, pS — двухэтапная модель компенсации сводится к одной большой задаче линейного программирования, её детерминированному эквиваленту :

Один экземпляр второго этапа на каждый сценарий, все приклеены к одному и тому же первому этапу xx. Передайте его в HiGHS или Gurobi, подождите и получите решение. Конец истории.

Но иногда ожидание становится всей историей.

Когда детерминированный эквивалент упирается в стену

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

• Для аппроксимации непрерывного распределения спроса с помощью метода выборочного среднего может потребоваться несколько тысяч выборок, прежде чем оптимальное значение xx стабилизируется.
• Многоступенчатая модель разрастается экспоненциально с увеличением числа ступеней. Три ступени, каждая из которых содержит десять ветвей, уже дают тысячу листьев; пять ступеней — и число листьев превысит сотни тысяч.
• В реальных приложениях, таких как планирование работы гидроэлектростанций, планирование цепочки поставок в условиях резких колебаний спроса или диспетчеризация энергии, десятки тысяч вариантов моделируемых сценариев — не редкость. Специалисты по гидроэнергетике регулярно создают деревья сценариев, которые слишком велики, чтобы их даже перечислить, не говоря уже о том, чтобы передать их решателю.

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

Это называется блочно-угловой структурой . Каждый сценарий вносит свой собственный диагональный блок WsW_s, плюс столбец TsT_s слева, который связывает его с xx. Общая строка внизу — это ограничение осуществимости первого этапа Ax≥bAx geq b.

Форма графика информативна, но для решателя это плохая новость. Время выполнения линейного программирования растет сверхлинейно с размером задачи. Быстрый поиск в Google покажет, что в худшем случае оно составляет примерно O(n².5)O(n².5) для симплексных методов и O(n³.5)O(n³.5) для методов внутренней точки, что означает, что удвоение количества сценариев не удваивает время выполнения — оно увеличивает его в шесть раз. В какой-то момент задача перестает помещаться в памяти, и ваш дружелюбный решатель начинает выгружать данные на диск, что является вежливым способом сказать: «Это слишком много».

Поэтому нам нужно что-то поумнее, чем «запихнуть всё в программу для решения задач и надеяться на лучшее».

Почему наивные уловки не помогают

Прежде чем перейти к самой интересной идее, давайте отбросим очевидные варианты.

Идея 1: Просто используйте меньше сценариев. Иногда это допустимо; иногда ответ существенно меняется по мере добавления сценариев, и ваш «оптимум» оказывается случайным совпадением с полученной вами выборкой. В литературе, посвященной аппроксимации среднего значения выборки, как раз и говорится о том, когда это допустимо, а когда нет.

Идея 2: Решить каждый сценарий отдельно и усреднить результаты. Заманчиво, но неправильно. Весь смысл наличия одного варианта xx заключается в том, что он должен работать для всех сценариев одновременно. Решение сценария ss само по себе дает вам оптимум для конкретного сценария («выжидательное» решение из предыдущего поста), а не решение на первом этапе, которому вы можете быть уверены.

Идея 3: Решить детерминированную задачу со средним спросом 𝔼[ξ], а затем разобраться с остатками. Это решение EV. Оно может быть очень ошибочным, и именно эту ошибочность измеряет VSS.

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

Это нечто существует с 1962 года.

Единственное наблюдение, которое нас спасает

Взгляните еще раз на матрицу ограничений:

Теперь представьте, что xx фиксировано, и выберите любое значение, неважно какое. Блок первого этапа Ax≥bAx geq b либо выполняется, либо нет. И остальные ограничения разделяются: ограничение сценария s становится Wsys≥hs−TsxW_s y_s ≥ h_s − T_s x, которое включает только yₛ. Связь между сценариями проходит через xx. Зафиксируйте xx, и сценарии распадаются на SS-независимые задачи линейного программирования.

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

Рецепт разложения Бендеров основан на этом наблюдении. Примерно так:

1. Угадайте решение на первом этапе xx.
2. Решите подзадачи сценария S независимо друг от друга, имея xx.
3. Используйте информацию из подзадач, чтобы скорректировать свою оценку значения xx.
4. Повторяйте до тех пор, пока ваши предположения не перестанут улучшаться.

Самая сложная (и одновременно элегантная) часть — это шаг 3: как подзадачи могут рассказать нам что-либо полезное о xx? Вот где проявляется ценность двойственности LP.

Бендеры!

Давайте замедлим темп и будем строить алгоритм по частям.

Функция ценности v(ξˢ, x)

При фиксированном решении на первом этапе xx и фиксированном сценарии ss подзадача второго этапа представляет собой всего лишь небольшую задачу линейного программирования:

Его двойственность

Эта двойственная задача имеет особенность, на которую мы будем опираться. Ее допустимая область Λs={λ≥0:λ⊤Ws≤qs}Lambda_s = {lambda geq 0 : lambda^top W_s leq q_s} представляет собой многогранник и совершенно не зависит от xx. xx появляется только в целевой функции. В сочетании с тем фактом, что в задачах линейного программирования оптимум всегда достигается в вершине допустимой области, мы получаем полезную переформулировку:

где λkslambda^s_k — это (конечное число) вершин ΛsLambda_s.

Это очень приятное выражение для размышления. Как функция от x, v(ξs,x)v(xi^s, x) является поточечным максимумом конечного набора аффинных функций . Это делает её кусочно-линейной и выпуклой . Представьте себе верхнюю огибающую веера прямых.

Ожидаемые издержки на устранение последствий представляют собой лишь взвешенную сумму этих факторов:

Взвешенная сумма кусочно-линейных выпуклых функций сама по себе является кусочно-линейной и выпуклой. Таким образом, Q(x)Q(x), то, что мы пытаемся минимизировать на первом этапе поверх c⊤xc^top x, является кусочно-линейной выпуклой функцией от xx. У нас нет явной формулы для неё, поскольку это потребовало бы перечисления каждой вершины каждого ΛsLambda_s, но мы знаем её форму.

Эта фигура — это рычаг, который тянет Бендерс.

Аппроксимация Q с помощью сечений

Мы не знаем QQ в замкнутой форме, но кусочно-линейная выпуклая функция обладает очень снисходительным свойством: в любой точке x‾bar{x}, где нам известно значение Q(x‾)Q(bar x) и субградиент g∈∂Q(x‾)g ∈ ∂Q(x̄), линейная функция

является глобальной нижней границей для QQ. Она касается QQ в точке x‾x̄ и остается ниже во всех остальных точках.

Что касается стоимости регресса, то субградиент можно вычислить. Если мы решим двойственные подзадачи и получим оптимальные двойственные вершины λk∗slambda^s_{k*} для каждого сценария, то субградиент QQ в точке x‾x̄ будет равен

Таким образом, каждый раз, когда мы решаем подзадачи в некоторой точке-кандидате xτx^tau, мы получаем точную аффинную нижнюю границу для QQ совершенно бесплатно. Бендерс называет каждую такую нижнюю границу оптимальным ограничением .

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

Главная проблема

Мы заменяем Q(x)Q(x) на переменную-заполнитель θtheta, которая должна лежать выше каждой нижней границы, которую мы накопили к настоящему моменту. Итерация tt имеет главную задачу :

Здесь LL — любая безопасная нижняя граница для QQ — часто L=0, если затраты на исправление неотрицательны, или что-то более отрицательное, если нет — и каждая (αt,βt)(αₜ, βₜ) представляет собой оптимальное отсечение из предыдущей итерации. Главная задача — это обычная задача линейного программирования, гораздо меньшая, чем её детерминированный эквивалент, поскольку она содержит только переменные первого этапа и постоянно увеличивающийся набор отсечений.

Затем решение-кандидат (xτ,θτ)(x^tau, θ^tau), полученное из этого основного набора данных, передается в подзадачи, которые возвращают новый результат, который добавляется, и мы повторяем весь процесс.

Алгоритм

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

1. Решите основную задачу, чтобы получить (xτ,θτ)(x^tau, θ^tau).

2. Для каждого сценария s=1,…,Ss = 1, …, S, решите двойственную подзадачу в точке xτx^tau, чтобы получить оптимальную вершину λτsλˢ_τ и оптимальное значение v(ξs,xτ)v(ξˢ, x^tau).

3. Постройте оптимальное ограничение θ≥αt+βt⊤xθ ≥ α_t + β_t^top x с

и добавить его в основной файл.

4. Остановитесь, если значение θτθ^tau у мастера уже удовлетворяет новому критерию в точке xτx^tau — то есть, если θτ≥αt+βt⊤xτθ^tau ≥ α_t + β_t^top x^tau. В противном случае вернитесь к шагу 1.

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

Конечность. Каждый разрез исходит из некоторой вершины λksλ^s_k некоторого множества ΛsLambda_s. Каждое множество ΛsLambda_s имеет конечное число вершин. Поэтому в худшем случае алгоритм перечисляет все возможные разрезы и останавливается — он не может зацикливаться бесконечно. На практике он останавливается гораздо раньше: обычно имеет значение лишь небольшое количество вершин.

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

Информация, а не перечисление. Мы никогда не записываем все сценарии в одном LP. Мы позволяем подзадачам, через их двойные переменные, рассказать нам то, что нужно знать о них главному модулю. Разделы — это язык, на котором они говорят.

Однослойная против многослойной нарезки

В приведенном выше рецепте мы объединяем информацию по каждому сценарию в один разрез на QQ. Естественная альтернатива — хранить отдельные заполнители SS θ1,…,θSθ₁, …, θ_S, по одному на каждый сценарий, и добавлять отдельные разрезы SS на каждой итерации. Это вариант с несколькими разрезами .

Многократное отсечение использует больше информации за итерацию, поэтому приближение QQ, полученное с помощью основного алгоритма, становится более точным быстрее. Недостатком является то, что за итерацию добавляется одно ограничение SS, и основной алгоритм быстрее становится громоздким. Универсального победителя нет — Ю и Гроссманн (2013) сообщают о значительном ускорении многократных сечений в моделях цепочек поставок, но это зависит от конкретного случая. Разумным вариантом по умолчанию является начало с однократного сечения и попытка многократных сечений, если сходимость медленная.

Изображение одной из итераций

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

Представьте Q(x)Q(x) как истинную (неизвестную) кусочно-линейную выпуклую функцию, построенную на графике xx. Текущее приближение мастера Q^τ(x)hat Q^tau(x) — это поточечный максимум всех сечений, которые мы добавили до сих пор, — полиэдральная нижняя граница, которая находится под Q повсюду.

На итерации τtau мастер минимизирует c⊤x+Q^τ(x)c^top x + Q̂^tau (x) и попадает на некоторое значение xτx^tau, где его приближение считает оптимумом. Затем мы задаём подзадачам вопрос: «Каково истинное значение QQ в точке xτx^tau?» Они отвечают как значением Q(xτ)Q(x^tau), так и субградиентом (двойственным решением). Этот субградиент даёт нам новую линейную функцию, которая касается QQ в точке xτx^tau — точку отсечения. Мы добавляем её; приближение мастера возрастает ровно настолько, чтобы быть точным в точке xτx^tau; и значение xτ+1x^{tau+1} на следующей итерации (обычно) будет другим.

Цикл останавливается, когда приближение мастера уже предсказывает истинное значение в точке xτx^tau. Визуально это означает, что наша полиэдральная нижняя огибающая догнала QQ точно в той точке, которая нас интересует.

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

Это не волшебный инструмент.

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

Медленная сходимость в реальных условиях. Теоретически алгоритм Бендерса завершается за конечное число итераций; на практике «конечное число» может означать довольно много итераций, и каждая из них решает задачу линейного программирования. Работа Рахманиани и др. (2017) представляет собой стандартный обзор литературы по десятку с лишним приемов, используемых для ускорения процесса: выбор отсечения, удаление отсечения, когда старые отсечения становятся доминирующими, области доверия для предотвращения резких колебаний xτx^tau и так далее.

Основная задача становится узким местом. Нередко более 90% общего времени выполнения уходит на решение (постоянно растущей) основной задачи. Распространенные решения включают в себя отказ от решения основной задачи до достижения полной оптимальности на ранних итерациях, удаление ограничений, которые больше не помогают, и использование структуры, когда основная задача сама по себе является структурированной задачей линейного программирования.

Подзадачи тоже могут доминировать. Особенно когда SS велик или каждая подзадача сама по себе является сложной задачей линейного программирования, в дело вступает второй этап. Решение в основном сводится к инженерным решениям: распараллеливание по различным сценариям, запуск каждой подзадачи с нуля на основе предыдущей итерации, объединение задач в пакеты, когда это возможно.

Скрытое предположение: относительно полная возможность решения. Все вышеизложенное неявно предполагало, что подзадача второго этапа является выполнимой для любого xx, которое может предложить мастер. Это свойство (относительно) полной возможности решения, о котором говорилось в предыдущем посте: для каждого выполнимого x первого этапа и каждого сценария ξξ существует возможность решения y≥0y ≥ 0, удовлетворяющая условию Wy≥h(ξ)−T(ξ)xWy ≥ h(ξ) − T(ξ)x.

Когда это свойство нарушается, некоторые значения xτx^tau делают подзадачу неразрешимой — на втором этапе линейного программирования нет подходящего значения y, поэтому неявная стоимость равна +∞+ infty, что означает, что xτx^tau изначально не был допустимым вариантом на первом этапе. Бендерс решает эту проблему с помощью ограничений на допустимость , которые выводятся путем решения вспомогательной задачи линейного программирования, измеряющей неразрешимость подзадачи, и использования ее двойственной задачи для построения гиперплоскости, которая отсекает xτx^tau из допустимой области основной задачи. Без алгебраических преобразований структура та же, что и у ограничений на оптимальность: решаем небольшую задачу линейного программирования, рассматриваем двойственную задачу, строим линейное ограничение, возвращаем задачу основной. Таким образом, в каждой итерации Бендерса вы либо добавляете ограничение на оптимальность (если каждая подзадача была допустимой), либо ограничение на допустимость (если какая-либо подзадача была недопустимой). Модель с гарантией полного возмещения убытков по-прежнему работает; вам просто нужны оба вида компенсации.

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

На начальных этапах, связанных с целочисленными переменными, возникают сложности. Если xx должно быть целым числом (например, решения о мощности, местоположении объекта), то задача становится целочисленной, и теория выпуклых матриц Бендерса теряет часть своей эффективности. Существуют расширения, такие как обобщенный Бендерс, целочисленная L-образная модель, модель ветвления и отсечения Бендерса и т. д., но они выходят за рамки этой статьи.

Всё это не означает, что метод Бендерса плох. Это означает, что это метод, а не чудодейственное средство. При вдумчивом использовании — с правильным разделением основной задачи на подзадачи, правильным вариантом сокращения и правильными приёмами ускорения — это разница между решением проблемы и её нерешением.

Подводя итоги

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

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

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

1. Функция ценности второго этапа v(ξ,x)v(ξ, x) является кусочно-линейной и выпуклой по xx. То же самое относится и к ожидаемой стоимости компенсации Q(x)Q(x).
2. Кусочно-линейную выпуклую функцию можно аппроксимировать с произвольной точностью с помощью растущего набора аффинных нижних границ (срезов), каждая из которых получается бесплатно из двойственного решения.
3. На первом этапе не нужно точно знать QQ — достаточно точно определить его нижнюю границу в интересующих точках. Цикл «главная задача/подзадача» как раз и обеспечивает это знание.
4. Однократное и многократное отсечение — это два способа объединения информации по каждому сценарию. Отсечение по возможности обрабатывает случай, когда подзадачи могут оказаться неразрешимыми.

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

В следующем посте этой серии мы рассмотрим, что делать, когда даже двухэтапная модель Бендерса оказывается неподходящей — когда мир не просто вежливо разделяется на «решить, наблюдать, исправить», а бесконечно повторяет «решить, наблюдать, решить, наблюдать, решить…». Это территория SDDP. Если вы когда-либо задавались вопросом, как операторы гидроэлектростанций планируют использование водохранилищ на многолетний период с тысячами сценариев совместного притока воды, ответ вы найдете в следующем посте.

Благодарности и ссылки

Как и в предыдущем посте, этот основан на лекциях доктора Рубена ван Беестена (Норвежский университет науки и технологий) с его курса по стохастическому программированию, прочитанного в октябре 2023 года в Тронхейме. Педагогическая структура взята непосредственно из его слайдов; любые неточности в пересказе — мои. Также все изображения, за исключением миниатюры, в этом посте созданы мной.

Четыре источника, о которых стоит знать:

• Бендерс, Дж. Ф. (1962). Процедуры разбиения для решения задач программирования со смешанными переменными. Numerische Mathematik, 4(1), 238–252. Оригинальная статья. Удивительно легко читается.
• Ван Слайк, Р., и Ветс, Р. (1969). L-образные линейные программы с приложениями к оптимальному управлению и стохастическому программированию. Журнал SIAM по прикладной математике, 17(4), 638–663. Версия, адаптированная под стохастическое программирование; отсюда и название «L-образный метод».
• Фэнцзи Ю и Игнасио Э. Гроссманн. (2013). Алгоритм декомпозиции многосекционных бендеров для планирования цепочки поставок в условиях неопределенности. Анналы операционных исследований, 210:191–211.
• Рахманиани, Р., Крайник, Т.Г., Жандро, М., и Рей, В. (2017). Алгоритм декомпозиции Бендерса: обзор литературы. Европейский журнал операционных исследований, 259(3), 801–817. Справочник по ускорениям, вариантам и современному использованию.

А вот предыдущие посты из этой серии: Пять вопросов, которые помогут вам лучше моделировать целочисленные линейные программы, Четыре шага для повышения устойчивости вашей линейной программы и Простое введение в стохастическое программирование.

Беренд Маркхорст Посмотреть все от Беренд Маркхорст

Источник: towardsdatascience.com