Декомпозиционный алгоритм для задач линейного программирования, не имеющих блочной структуры, в среде Everest Текст научной статьи по специальности «Математика»

Декомпозиционный алгоритм для задач линейного программирования, не имеющих блочной структуры, в среде Everest

В.В. Волошинов

Аннотация — В статье описывается декомпозиционный метод решения задач линейного программирования общего вида, когда в матрице ограничений нельзя или затруднительно выделить блочную структуру, используемую классическими декомпозиционными алгоритмами Данцига-Вульфа или Бендерса. Вместо этого применяется произвольное разделение матрицы ограничений на некоторое количество подматриц, соответствующих группам строк. Как и в алгоритме Данцига-Вульфа, вычислительный метод соответствует поиску максимума вогнутой функции, зависящей от двойственных переменных, субградиентным методом. Решение исходной задачи проводится итеративно, включая этапы решения наборов независимых подзадач меньшей размерности, чем у исходной. Выбор числа подматриц (и подзадач) зависит от размерности исходной задачи и вычислительной мощности распределенной среды, где при возможно параллельное решение указанных независимых подзадач. Конкретно, предлагается использовать сервисы решения задач оптимизации, развернутые на платформе Everest, everest.distcomp.org. Для программирования алгоритма расчетов и обмена данными с решателями, подключенными к Everest, можно использовать систему AMPLX, использующую алгебраический язык оптимизационного моделирования AMPL, ampl.com. Вместо коммерческого транслятора AMPL возможно применять интерпретатор языка Python и свободно доступный пакет Pyomo, pyomo.org, обеспечивающий взаимодействие с AMPL-совместимыми решателями. Программная реализация алгоритма в распределенной среде сервисов оптимизация Everest позволяет получить заметное ускорение при решении масштабных задач.

Ключевые слова — линейное программирование, лагранжевая релаксация, ленточная декомпозиция, распределенные вычисления, платформа Everest.

I. Введение

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

Здесь описывается декомпозиционная схема решения задач линейного программирования. Рассматривается общий случай, когда в матрице ограничений нельзя или затруднительно выделить блочную структуру, используемую классическими декомпозиционными алгоритмами Данцига-Вульфа или Бендерса [1], [2]. Вместо этого предлагается использовать, вообще

Статья получена 20 декабря 2018. Работа выполнена за счет гранта Российского научного фонда (проект №16-11-10352). В.В. Волошинов, Институт проблем передачи информации им. А.А. Харкевича РАН (e-mail: [email protected]).

говоря, произвольную ленточную строчную («горизонтальную») или столбцовую («вертикальную») декомпозицию матрицы ограничений на некоторое количество подматриц. Количество этих подматриц зависит от размерности исходной задачи и вычислительной мощности распределенной среды, где происходит выполнение алгоритма. Конкретно предлагается использовать сервисы решения задач оптимизации, развернутые на платформе Everest, everest.distcomp.org, [3]. Для программирования алгоритма расчетов и обмена данными с решателями, подключенными к Everest, можно использовать систему AMPLX [4], использующую один из самых популярных алгебраических языков оптимизационного

моделирования AMPL, ampl.com.

Также, вместо лицензированного, коммерческого, транслятора AMPL возможно применять интерпретатор языка Python и свободно доступный в исходных кодах пакет Pyomo (PYthon Optimization Modeling Objects) [6], pyomo.org, обеспечивающий взаимодействие с AMPL-совместимыми решателями. Общая схема применения Pyomo для расчетов в среде сервисов оптимизации Everest приведена в статье [7].

Предлагаемый алгоритм был разработан автором самостоятельно и казался оригинальным. Однако, в книге [8], Гл. 2, в п. 5 «Метод расслоения ограничений» описывается похожая схема для случая строчной декомпозиции на две подматрицы. Там же указаны ссылки на более ранние работы ([9], [10]), где применялась такая декомпозиция.

II. Общая схема декомпозиции

Рассмотрим задачу линейного программирования достаточно общего вида:

(1)

у которой, либо заведомо нет, либо неизвестно и трудно выявить характерные свойства «блочности»: т. е. наличие выделенных групп переменных и/или ограничений (см. Рис. 1. Схематичная структура блочных матриц.

1

Пусть матрица А ограничений задачи ( )

разрезана на горизонтальные «ленты» (см. Рис. Ошибка! Неверная ссылка закладки.). Здесь матрицы, имеют размерность и

соответственно,

где

Рис. 2). Также возможно ситуация, когда блочная структура в задаче большой размерности наличествует, но она слишком «крупна», чтобы дать заметный прирост производительности. Во всех указанных случаях непосредственное применение известных

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

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

Рис. 1. Схематичная структура блочных матриц. Пусть (тхп) матрица А

1

А), и Ак (к 1: /\ ; ограничений задачи ( ) разрезана на горизонтальные «ленты» (см. Рис. Ошибка! Неверная ссылка закладки.). Здесь

матрицы, имеют

то+ T,k=l:K тк = ГП.

размерность соответственно,

(т0х-п)

где

Рис. 2. «Строчная» декомпозиция 1

Тогда задача ( ) может быть записана в следующей эквивалентной форме:

со «вспомогательными» переменными и

соответствующими ограничениями (вектора -'о и Ьк(к—1:К) _ результат разделения компонент вектора правой части ограничений Ь). Слева от вертикальных 2

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

2

двойственную задачу для ( ), укажем вид двойственной 1

задачи для ( ):

(3)

Теперь, «по аналогии», запишем двойственную задачу 2

для ( )

(4)

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

(5)

Нетрудно заметить, что оптимальное значение

4

критерия в двойственной задаче ( ) совпадает со значением максимума суммы вышеуказанных функций:

(6)

Прежде, чем перейти к описанию схемы алгоритма

и

(7)

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

5

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

По правилу строчной декомпозиции этот вектор можно записать в виде У°~ (&о> У° > ■ • ■ > Ук> ■ ■ ■ > У к), где подвектора Уо, Ук (к=1'-к) соответствуют подматрицам 2

А0, Ак (см. ( )). По определению, для этих подвекторов выполнены неравенства:

у^Ао + Y, vfAk S ст

к=1:К

Положим

4 = yfAk Из ( ) и

(8)

определении получаем:

III. Описание алгоритма

Алгоритмически, предлагаемая вычислительная схема соответствует поиску максимума вогнутой функции Ф(г) 6

( ) субградиентным методом [2]. Вогнутость этой функции следует из вогнутости ее слагаемых

7

фо(г), (pk(z) (Ä—l:iO что видно из вьфажений ( ) (минимум аффинных по г функций по множествам, независящим от г).

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

Теорема III. 1. Пусть множества допустимых решений^ = Ах=Ь, х -— 0} и

i ) У А = с }, соответственно, в исходной и двойственной к ней задачах не пусты. Тогда:

Доказательство. Если Q т : , то и все множества Q о, 7

Qk (k=l:K)^ ( ^ также не пусты. Это следует из очевидных включений Q(~Qo и Q'-~Qk (k=l:K)_ Таким образом можно считать, что п. 1) доказан. Другими словами, для произвольного вектора z, значения всех функций Фо(з)-, 4>k(ß) (k=l:K) ограничены сверху. Для доказательства п. 2) построим такой вектор г',

5

чтобы все множества --о(-г")-(k=l:R) из ( ) были бы непусты. Рассмотрим произвольный вектор у0

Теорема доказана.

Заметим, что если все коэффициенты целевой функции исходной задачи неотрицательны, ' =0 (что часто имеет место на практике), то можно положить

Таким образом, из совместности прямой и двойственной задач следует, что все функции Фо(х), Фк(.(к=1:К) являются собственными вогнутыми функциями, а их значения определяются в результате

7

решения задач линейного программирования ( ).

Пусть для некоторого I, после решения указанных подзадач найдены вектора жо И, 4 И ... Х1И ... хк И:

Как известно из выпуклого анализа, эти соотношения позволяют определить вектора из субдифференциалов функций М*)» Фк(~) (к=ЪК)вТочке г:

(-4W-.

1'0 L--J • ■

К'П

(9)

(0,

,0,4[г], 0,... ,0) €вфк{г)сЖк-п (к=1 :К)

из непустого множества Д т. е. У'

оТ

А < ,

В первой вектор-строке К раз повторяется п-мерный вектор 4 И, взятый с обратным знаком. Во второй

9

группе соотношений ( ), для к-й функции все компоненты субградиента равны нулю, кроме к-й

группы из п компонент вектора хк И. Из ( ) следует, что все элементы субдифференциала $$'(-5) могут быть представлены в виде сумм указанных векторов:

5

(10)

Проверим, что принадлежность нуля соответствует решению исходной задачи. Действительно, если для некоторого 0€£?Ф(з*)5 т0 эт0

означает, что нашелся вектор --о L-- J, для которого

Это можно записать в эквивалентной форме.

Но последнее и означает, что

- решение

останавливается, выдавая решение хо. Если это не так, то формируется вектор

(13)

9

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

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

Шаг 0. Выбор начального вектора г. На этом шаге нужно выбрать некоторый вектор г°, для которого значение функции Ф(г0) конечно. Здесь можно использовать Теорему 111.1. Например, если вектор с неотрицателен, то можно положить г0 равным нулю. Если удалось найти какое-нибудь допустимое решение двойственной задачи, то можно положить 2 к = У к

Шаг 2$+1 (s=0,1,...). Решение подзадач. На этом, «нечетном», шаге решаются независимые подзадачи (всего — (К+1) подзадач)

который, как следует из ( ) принадлежит субдифференциалу ¿^'("''Ч Естественно теперь сделать «шаг» вдоль субдифференциала в надежде получить большее значение функции Ф(г). Однако длину такого шага надо заранее ограничить, чтобы остаться в множестве конечных значений Ф(г). Используем для этого двойственное определение слагаемых этой

5

функции соотношениями ( ). Заметим, что найдя 6

решения задач ( ) любым ЛМРЬ-совместимым пакетом, помимо векторов «прямых» переменных, становятся известны соответствующие вектора двойственных переменных Уо > Ух > У-г > • • • < У к > ■ • ■ > Ук, являющиеся

5

решениями двойственных задач ( ). Исходя из этого

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

(.+1)

вектора г по правилу:

где число определяется в результате решения следующей «координирующей»-задачи (мастер-задача в англоязычной терминологии)

Шаг 2(s+1) ^=0,1,...). Проверка критерия остановки и смена вектора г..

На этом, «четном», шаге вначале проверяется выполнение с требуемой точностью равенств

Если равенства выполнены, то алгоритм

(14)

IV. О ВОЗМОЖНОСТЯХ УСКОРЕНИЯ И ПРИМЕНИМОСТИ ПРЕДЛОЖЕННОГО АЛГОРИТМА

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

Во-первых, «нечетный» шаг алгоритма ((2в+1)-й) заключается в поиске решения (К+1) задач линейного программирования с матрицами ограничений размера (т„хп) и {'411.,, п) (к=1:К )_ Несмотря на то, что число переменных в этих задачах остается прежним (п), число строк может быть, вообще говоря, выбрано произвольным. Требуется лишь, чтобы выполнялось равенство тоо+ 52к=1:К п1к = т.

Учтем теперь давно известное «эмпирическое» правило, согласно которому трудоемкость симплекс-

3

метода пропорциональна т п, где т - число

35

ограничений (без учета условий неотрицательности на переменные), а n - размерность вектора переменных (см. [11]). Это замечание позволяет надеяться, что даже если декомпозиционный алгоритм совершит большее число итераций (чередуя параллельное решение независимых подзадач с решение одной мастер-задачи), чем тот же

1

симплекс, примененный к исходной задаче ( ), то общее время поиска решения сократится. Например, если матрицу исходной задачи «разрезать» на 100 «горизонтальных лент» (K=99), то нечетный шаг алгоритма, при наличии сотни доступных пакетов решения ЛП-задач, подключенных к системе AMPLX, может выполняться, примерно в 106 (миллион!) раз быстрее, чем решение исходной, «полноразмерной», задачи. Даже если общее число итераций будет больше, чем потребовалось бы для исходной задачи, эта «фора» параллельной фазы декомпозиционного алгоритма позволяет надеяться на ускорение.

Условием практической применимости этой схемы является наличие удобных программных средств: декомпозиции исходной задачи на указанные подзадачи; формирования «промежуточной» задачи четного шага (пересчет вектора z); компоновки исходных данных новых подзадач по новому вектору z. Все эти действия удобно совершать именно в системе AMPLX, поскольку она основана на языке высокоуровневом оптимизационного моделирования AMPL. Кроме возможностей автоматического создания исходных данных всех задач в виде т. н. стаб-файлов (их часто еще называют NL-файлами), которые можно непосредственно отправлять AMPL-совместимым пакетам, результаты решения, SOL-файлы импортируются в AMPL-программу с решениями, полученными от пакетов (вместе с двойственными переменными).

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

Возможности AMPLX по реализации алгоритмов Данцига-Вулфа и Бендерса были успешно проверены на демонстрационных примерах их реализации на языке AMPL. Читатель может сравнить примеры по именам multi.mod (*.dat) и benders.mod (*.dat) на сайте ampl.com, страница https://ampl.com/resources/the-ampl-

book/example-files, и соответствующую модификацию для системы AMPLX по ссылкам http://distcomp.ru/~vladimirv/restopt/amplx/dw/ и

http://distcomp.ru/~vladimirv/restopt/amplx/benders/.

На тех же принципах возможна реализация вычислений на языке Python средствами пакета Pyomo (см. [7]).

У.Заключение

Разработан декомпозиционный алгоритм решения задачи линейного программирования общего вида, у которой, либо заведомо нет, либо трудно выделить «блочные» группы переменных и/или ограничений,

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

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

Исследование выполнено за счет гранта Российского научного фонда (проект №16-11-10352).

Библиография

[1] Лэсдон Л.С, Оптимизация больших систем. М.: Наука, 1975.

[2] Мину М. Математическое программирование. Теория и алгоритмы: Пер. с фр. А. И. Штерна. М.: Наука, 1990.

[3] Sukhoroslov O., Volkov S., Afanasiev A. "A Web-Based Platform for Publication and Distributed Execution of Computing Applications", in Parallel and Distributed Computing (ISPDC), 14th International Symposium, IEEE, 2015, pp. 175-184.

[4] Smirnov S., Voloshinov V., Sukhosroslov O. "Distributed Optimization on the Base of AMPL Modeling Language and Everest Platform" in Procedia Computer Science, vol. 101, 2016, pp. 313322.

[5] Robert Fourer, David M. Gay, Brian W. Kernighan. AMPL: A mathematical programming language. Thomson/Brooks/Cole, 2003

[6] Hart, William E., Carl D. Laird, Jean-Paul Watson, David L. Woodruff, Gabriel A. Hackebeil, Bethany L. Nicholson, John D. Siirola. Pyomo - Optimization Modeling in Python. 2nd Edition. Vol. 67. Springer, 2017

[7] Афанасьев А.П., Волошинов В.В., Соколов А.В. «Обратные задачи моделирования на основе регуляризации и распределенных вычислений в среде Everest». В сборнике: Аналитика и управление данными в областях с интенсивным использованием данных Сборник научных трудов XIX Международной конференции DAMDID / RCDL'2017. Под ред. Л.А. Калиниченко, Я. Манолопулос, Н.А. Скворцова, В.А. Сухомлина. 2017. С. 132-140.

[8] Авербах И.Л., Цурков В.И. Оптимизация в блочных задачах с целочисленными ограничениями. М.: Наука, 1995

[9] Guignard M., Kim S. Lagrangean decomposition: A model yielding stronger lagrangean bounds, in Mathematical Programming, vol. 39, iss. 2, 1987, pp. 215-228

[10] Glover F., Klingman D. Layering strategies for creating exploitable structure in linear and integer programs, in Mathematical Programming, vol. 40, iss. 1-3, 1988, pp. 165-181

[11] Муртаф Б. Современное линейное программирование. М.: Мир, 1984.

Decomposition algorithm for linear programming problems without a block structure in the Everest computing environment

V.V. Voloshinov

Abstract — The article describes decomposition method for solving linear programming problems in a general case when it is impossible or difficult to reveal the block structure (of the constraint matrix) used by the classic Danzig-Wulf or Benders decomposition algorithms. Instead, an arbitrary partition of the constraint matrix into a number of submatrices corresponding to groups of rows is used. As in the Danzig-Wolfe algorithm, one searches a maximum of a concave function (on dual variables) by a subgradient method. Iterative routine of solving the original problem includes the phases of solving sets of independent subproblems of a smaller dimension than the original one. Number of submatrices (and subproblems) depends on the dimension of the original problem and the computational power of the distributed environment, where parallel solving of above independent subproblems is done. In concrete term, it is proposed to use the optimization services deployed on the Everest platform, everest.distcomp.org. Programming of the algorithm and data exchange with solvers connected to Everest, may be done by AMPLX system based on algebraic optimization modeling language AMPL, ampl.com. Also, instead of the commercial translator AMPL, it is possible to use the Python interpreter and the freely available Pyomo package, pyomo.org, which provides interaction with AMPL-compatible solvers. The software implementation of the algorithm by Everest platform, allows to hope to get a noticeable acceleration when solving large-scale problems.

Keywords — linear programming, Lagrangian relaxation, tape decomposition, distributed computing, Everest platform.

References

[1] Ljesdon L.S, Optimizacija bol'shih sistem. M.: Nauka, 1975.

[2] Minu M. Matematicheskoe programmirovanie. Teorija i algoritmy: Per. s fr. A. I. Shterna. M.: Nauka, 1990.

[3] Sukhoroslov O., Volkov S., Afanasiev A. "A Web-Based Platform for Publication and Distributed Execution of Computing Applications", in Parallel and Distributed Computing (ISPDC), 14th International Symposium, IEEE, 2015, pp. 175-184.

[4] Smirnov S., Voloshinov V., Sukhosroslov O. "Distributed Optimization on the Base of AMPL Modeling Language and Everest Platform" in Procedia Computer Science, vol. 101, 2016, pp. 313-322.

[5] Robert Fourer, David M. Gay, Brian W. Kernighan. AMPL: A mathematical programming language. Thomson/Brooks/Cole, 2003

[6] Hart, William E., Carl D. Laird, Jean-Paul Watson, David L. Woodruff, Gabriel A. Hackebeil, Bethany L. Nicholson, John D. Siirola. Pyomo - Optimization Modeling in Python. 2nd Edition. Vol. 67. Springer, 2017

[7] Afanas'ev A.P., Voloshinov V.V., Sokolov A.V. «Obratnye zadachi modelirovanija na osnove reguljarizacii i

raspredelennyh vychislenij v srede Everest». V sbornike: Analitika i upravlenie dannymi v oblastjah s intensivnym ispol'zovaniem dannyh Sbornik nauchnyh trudov XIX Mezhdunarodnoj konferencii DAMDID / RCDL'2017. Pod red. L.A. Kalinichenko, Ja. Manolopulos, N.A. Skvorcova, V.A. Suhomlina. 2017. S. 132-140.

[8] Averbah I.L., Curkov V.I. Optimizacija v blochnyh zadachah s celochislennymi ogranichenijami. M.: Nauka, 1995

[9] Guignard M., Kim S. Lagrangean decomposition: A model yielding stronger lagrangean bounds, in Mathematical Programming, vol. 39, iss. 2, 1987, pp. 215-228

[10] Glover F., Klingman D. Layering strategies for creating exploitable structure in linear and integer programs, in Mathematical Programming, vol. 40, iss. 1-3, 1988, pp. 165181

[11] Murtaf B. Sovremennoe linejnoe programmirovanie. M.: Mir, 1984.