CC BY
0
УДК 519.852, 004.021, 004.032.24
DOI: 10.14529/ cmse240301
ЧИСЛЕННАЯ РЕАЛИЗАЦИЯ МЕТОДА ПОВЕРХНОСТНОГО ДВИЖЕНИЯ ДЛЯ РЕШЕНИЯ ЗАДАЧ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ
© 2024 Л.Б. Соколинский, H.A. Ольховский, И.М. Соколинская
Южно-Уральский государственный университет (454080 Челябинск, пр. им. В.И. Ленина, д. 76) E-mail: [email protected], [email protected], [email protected]
Поступила в редакцию: 10.06.2024
Работа посвящена численной реализации нового метода линейного программирования, получившего название «метод поверхностного движения». В основе реализации лежит оригинальный алгоритм AlFaMove, который строит на поверхности допустимого многогранника оптимальный целевой путь от произвольной граничной точки до точки, являющейся решением задачи линейного программирования. Оптимальность пути заключается в том, что направление движения по грани многогранника соответствует максимальному увеличению значения целевой функции. Для вычисления оптимального направления движения используется метод, базирующийся на операции построения псевдопроекции на линейное многообразие. Операция псевдопроекции обобщает понятие ортогональной проекции и реализуется с помощью итерационного алгоритма проекционного типа. Доказано, что в случае линейного многообразия, образуемого путем пересечения гиперплоскостей, псевдопроекция совпадает с ортогональной проекцией. Также доказано, что в случае линейного многообразия метод на основе псевдопроектирования вычисляет вектор движения в направлении максимального увеличения целевой функции. Выполнена параллельная реализация алгоритма AlFaMove. Приведены результаты вычислительных экспериментов на кластерной вычислительной системе, демонстрирующие высокую масштабируемость предложенной численной реализации.
Ключевые слова: линейное программирование, метод поверхностного движения, численная реализация, алгоритм AlFaMove, параллельная реализация, кластерная вычислительная система, исследование масштабируемости.
ОБРАЗЕЦ ЦИТИРОВАНИЯ
Соколинский Л.Б., Ольховский H.A., Соколинская И.М. Численная реализация метода поверхностного движения для решения задач линейного программирования // Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2024. Т. 13, № 3. С. 5-31. DOI: 10.14529/cmse240301.
Введение
Эпоха больших данных и индустрия 4.0 породили задачи линейного программирования (ЛП) сверхбольших размерностей, включающих в себя миллионы переменных и миллионы ограничений [1—4]. Во многих случаях объектом линейного программирования являются задачи, связанные с оптимизацией нестационарных процессов [5]. В нестационарных задачах ЛП целевая функция и/или ограничения изменяются в течение вычислительного процесса. Также среди этого класса задач встречаются приложения, в которых необходимо выполнять оптимизацию в режиме реального времени. Для решения таких задач необходимы масштабируемые методы и параллельные алгоритмы линейного программирования.
Один из стандартных подходов к решению нестационарных задач оптимизации состоит в том, чтобы рассматривать каждое изменение как появление новой задачи оптимизации, которую необходимо решать с нуля [5]. Однако такой подход часто непрактичен, поскольку решение проблемы с нуля без повторного использования информации из прошлого может
занять слишком много времени. Таким образом, желательно иметь алгоритм оптимизации, способный непрерывно адаптировать решение к изменяющейся среде, повторно используя информацию, полученную в прошлом. Этот подход применим для процессов реального времени, если алгоритм достаточно быстро отслеживает траекторию движения оптимальной точки. В случае больших задач ЛП последнее требует разработки масштабируемых методов и параллельных алгоритмов ЛП.
До настоящего времени наиболее популярными методами решения задач ЛП являются симплекс-метод [6] и методы внутренних точек [7]. Эти методы способны решать задачи с десятками тысяч переменных и ограничений. Однако масштабируемость параллельных алгоритмов, основанных на симплекс-методе, в общем случае ограничивается 16-32 процессорными узлами [8]. Что касается алгоритмов внутренних точек, они не поддаются эффективному распараллеливанию в общем случае. Это ограничивает применение указанных методов для решения сверхбольших нестационарных задач ЛП в режиме реального времени. В соответствии с этим задача разработки масштабируемых методов и эффективных параллельных алгоритмов ЛП для кластерных вычислительных систем остается актуальной.
В недавней работе [9] было дано теоретическое описание нового метода ЛП — метода поверхностного движения, строящего на поверхности допустимого многогранника1 оптимальный целевой путь к решению задачи ЛП. Под оптимальным целевым путем понимается путь по поверхности допустимого многогранника в направлении наибольшего увеличения значений целевой функции. Однако предложенный в этой статье алгоритм 1 на шаге 15 требует нахождения на границе гипердиска точки с максимальным значением целевой функции. При этом не приводится численный алгоритм, позволяющий выполнить этот шаг. В этой статье мы приводим и исследуем алгоритм AlFaMove, устраняющий допущенный пробел.
Статья организована следующим образом. Раздел 1 содержит теоретический базис, необходимый для описания метода поверхностного движения и его численной реализации. Раздел 2 посвящен описанию операции псевдопроекции, позволяющий найти вектор движения по оптимальному целевому пути для линейного многообразия, получаемого в результате пересечением гиперплоскостей. В разделе 3 дается формализованное описание алгоритма AlFaMove, представляющего собой численную реализацию метода поверхностного движения. Раздел 4 посвящен описанию параллельной версии алгоритма AlFaMove. В разделе 5 представлены информация о программной реализации алгоритма AlFaMove и результаты экспериментов на кластерной вычислительной системе по исследованию его масштабируемости. В заключении суммируются полученные результаты и намечаются направления дальнейших исследований.
1. Теоретический базис
Данный раздел содержит необходимый теоретический базис, используемый для описания алгоритма движения по граням AlFaMove. Рассмотрим задачу ЛП в следующем виде:
х = arg max {(с, х)\Ах ^ b}, (1)
хеш™
где с G Rn, Ь £ Rm, А G Rmxn, т > 1, с ф 0. Здесь (-, ■) обозначает скалярное произведение двух векторов. Мы предполагаем, что ограничение х ^ 0 также включено в матричное
Допустимый многогранник — область допустимых решений задачи линейного программирования.
неравенство Ах ^ Ь в форме
-х ^ 0.
Линейная целевая функция задачи (1) имеет вид
/(ж) = (с, х).
Вектор с в данном случае является градиентом целевой функции /(ж).
Пусть щ € Мп обозначает вектор, представляющий г-тую строку матрицы А. Мы предполагаем, что щ ф 0 для всех г е {1 ,...,т}. Обозначим через Щ замкнутое полупространство, определяемое неравенством (щ,х) ^ Ь^, а через Щ — ограничивающую его гиперплоскость:
Я* = {жеМп|(а*,ж)О0; (2)
Н{ = {х£Жп\(сц,х) = Ь{}. (3)
Определим допустимый многогранник
т
М=р(4)
г=1
представляющий множество допустимых точек задачи ЛП (1). Заметим, что М в этом случае будет замкнутым выпуклым множеством. Мы будем предполагать, что множество М является ограниченным и М ф 0, то есть задача ЛП (1) имеет решение. Дадим определение рецессивного полупространства [10].
Определение 1. Полупространство Н{ называется рецессивным, если
Ухет,У\>0:х + \с<£Щ. (5)
Геометрический смысл этого определения состоит в том, что луч, исходящий в направлении вектора с из любой точки гиперплоскости, ограничивающей рецессивное полупространство, не имеет общих точек с этим полупространством, за исключением начальной. Известно [10], что следующее условие является необходимым и достаточным для того, чтобы полупространство Щ было рецессивным:
(сц, с) > 0.
Определим
Х = {ге{1,...,т}|(а4,с)>0}! (6)
то есть X представляет множество индексов, для которых полупространство Щ является рецессивным. Поскольку допустимый многогранник М представляет собой ограниченное множество, имеем
Положим
М=р\Щ. (7)
гбХ
Очевидно, что М является выпуклым, замкнутым, неограниченным многогранником. Будем называть его рецессивным. Из (4) и (6) следует
М СМ.
Обозначим через Г(М) множество граничных точек допустимого многогранника М, а через Г(М) — множество граничных точек рецессивного многогранника М2. Согласно утверждению 3 в [10] имеем
х е Г(М),
то есть решение задачи ЛП (1) лежит на границе рецессивного многогранника М. Следуя [11], дадим определение ортогональной проекции на гиперплоскость.
Определение 2. Пусть в пространстве Мп имеется гиперплоскость
Н = {х еГ|(а,ж) = Ъ}. Ортогональная проекция 7Г#(г>) точки V е Мп на гиперплоскость Н определяется формулой
. . (а.у) — Ь
=«- ,,,2 а■ (8) ||а||
Следующее утверждение дает способ вычисления оптимального пути на гиперплоско-
сти.
Утверждение 1. Пусть в пространстве Мп задана гиперплоскость Н с нормалью аёК", проходящая через точку и £ Мп:
Н = {х еМп|(а,ж) = (а,и)}. (9)
Пусть задана линейная функция /(ж) : Мп —> Ж с градиентом с 6 Мп:
/(х) = (с,х). (10)
Пусть векторы а и с линейно независимы (не коллинеарны, и среди них нет нулевого вектора). Положим
V = и + с. (11)
Построим ортогональную проекцию 7Г#(и) точки V на гиперплоскость Н:
го = 7г#(г;). (12)
Тогда вектор с1 = гю — и однозначно задает направление максимального увеличения линейной функции /(аз), определяемой формулой (10).
Доказательство. Предположим противное, то есть существует точка го € Н такая, что
(с, го) ^ (с, го), (13)
||го — и|| = ||го — и|| и го ф го (см. рис. 1). Здесь и далее || • || обозначает евклидову норму. Вычислим (с,го). В соответствии с определением 2 ортогональная проекция 7Г#(«) точки
2Под граничной точкой множества М С Кп понимается точка в К™, для которой любая открытая ее окрест-
ность в К™ имеет непустое пересечение как с множеством М, так и с его дополнением.
Рис. 1. Иллюстрация к доказательству предложения 1. Пунктиром обозначена гиперокружность с радиусом ||го — и\\
V на гиперплоскость Н, задаваемую формулой (9), вычисляется следующим образом:
IV = V —
(а, V — и)
а.
а
Подставив вместо V правую часть формулы (11), получим
IV = и + с
{а, с)
а.
а
Используя (14), находим
/ \ / \ I II 1|2 (а1 С)
(с,го) = {с,и) + ||с||--и 2 .
\а\
(14)
(15)
Поскольку а и с линейно независимы, в соответствии с неравенством Коши-Буняковского имеем
\2 ^ 11_112 ||„||2
Отсюда следует
(а,сУ < ||а|| • ||с|| . (а, с)2
|с|12-
> 0.
а
(16)
Теперь вычислим (с,го). Определим й = 7Гь(го) — ортогональная проекция точки чЪ на прямую Ь, проходящую через точки и и го. По построению существует число 5 е К, удовлетворяющее условию
-1 < <5 < 1, (17)
такое, что
Положим
й = и + й(го — и).
V = и + — и).
(18)
Тогда точка й является ортогональной проекцией точки V на гиперплоскость Н, определяемую формулой (9), и может быть вычислена следующим образом:
и = V —
(а, V — и)
а.
а
Подставив вместо V правую часть формулы (18), отсюда получим
. ( (а.у — и) и = и + о уи — и---——а
Используя (11), произведем замену V — и = с:
й = и + . (19)
Очевидно,
го = (го - и) + й. (20)
Заменим второе слагаемое в (20) на правую часть формулы (19):
/ ~ ~ч г( (о, с) ги = [ги — и) +и + о 1с--^-а
Отсюда получаем
(с, го) = (с, ги — й) + (с, и) + 6 ( ||с||2 — ^
а
По построению вектор го — й ортогонален вектору v — u = c. Поэтому (с, го — и) = 0. Таким образом, последняя формула преобразуется к виду
(с,гЬ) = (с,и) + б(^\с\\2-{-^У (21)
Сопоставляя (15) и (21), с учетом (16) и (17) получаем
(с, го) < (с, го),
что противоречит (13). □
Возвращаясь к задаче ЛП (1), можно сказать следующее. Пусть и £ М П Г (М) и существует единственная рецессивная гиперплоскость Щ> (г' £ X) такая, что и £ Щ'. В этом случае вектор в,, определяющий направление оптимального целевого пути из точки и, в соответствии с утверждением 1 вычисляется по формуле
6 = с-{Щ''и + С\-Ь*щ,. (22)
НаИ1
В следующем разделе мы рассмотрим общий случай, когда через точку и проходит несколько гиперплоскостей.
2. Операция псевдопроектирования на линейное многообразие
Пусть J С {1 ,...,т}, и ф и Р| Щ ф 0. В этом случае множество индексов J
геЗ
определяет в пространстве Мга линейное многообразие
Ь = р| Щ. (23)
г^а
Обозначим кь — размерность линейного многообразия Ь. При кь < п — 1 многообразие Ь не является гиперплоскостью, и для определения вектора движения (I по этому многообразию в направлении максимального увеличения целевой функции нельзя применить формулу (22), так как такое линейное многообразие невозможно задать одним линейным уравнением в пространстве Мп. Однако мы можем найти указанный вектор с£ с помощью операции псевдопроектирования [10]. Определим проекционное отображение </>(■):
Известно [12], что отображение <р(х) является непрерывным Ь-фейеровским отображением, и последовательность точек
[хк = срк(х о)}^°=1, (25)
порождаемая этим отображением, начиная с произвольной точки Хо Е Мп, сходится к точке, принадлежащей Ь:
Хк х Е Ь.
Теперь мы готовы дать определение псевдопроекции на линейное многообразие, образуемое пересечением гиперплоскостей.
Определение 3. Пусть »7 С {1,... , т}, 3 ф 0, ф 0, (р(-) — проекционное отоб-
ражение, определяемое формулой (24). Псевдопроекцией PJ(x) точки х Е Мга на линейное многообразие Ь = Пге ,7 ^ называется предельная точка последовательности (25):
Нт
к—ьоо
р^х)-ук{х)
= 0.
Важным свойством псевдопроекции на линейное многообразие является то, что в этом случае псевдопроекция совпадает с ортогональной проекцией. Для доказательства этого факта нам понадобится следующая лемма.
Лемма 1. Пусть в пространстве Мп заданы гиперплоскость Щ/ и линейное многообразие Ь, принадлежащее этой гиперплоскости:
Нг> = {жбМ"| (а#,х) = Ьг'}; I' ЕЗ.
Обозначим через Р линейное многообразие, являющееся ортогональным дополнением к Ь:
Р = Ь±. (26)
Тогда для любой точки V, принадлежащей линейному многообразию Р, ее ортогональная проекция 7Гд./ (и) на гиперплоскость также будет принадлежать линейному многообразию Р:
Уу ЕР: 7ГЯ(у) Е Р.
Доказательство. Обозначим р — ортогональная проекция точки V на гиперплоскость Щ>:
_Р = *нЛ*)-_(27)
2024, т. 13, № 3 11
Без ограничения общности мы можем считать, что р € Ь (см. рис 2). Предположим, что точка р не принадлежит линейному многообразию Р. Обозначим q — точка пересечения линейного многообразия Ь с его ортогональным дополнением Р:
Я = ЬП Р.
Поскольку Р является ортогональным дополнением к линейному многообразию Ь, точка д является ортогональной проекцией точки р на линейное многообразие Р:
Ч = 1ГР(р). (28)
Рассмотрим треугольник Д(г;,р, д). В силу (27) угол ¿.р с вершиной в точке р является прямым. В в соответствии с (28) угол /д с вершиной в точке q также является прямым. Но это возможно только в случае, когда р = д, то есть
р Е Р.
Лемма доказана. □
Следующее утверждение доказывает, что псевдопроекция на линейное многообразие совпадает с ортогональной проекцией.
Утверждение 2. Пусть выполняются следующие условия:
•7С{1,...,т}, (29)
Ь = П Ни Ь + 0; (30)
где Нг = {ж € Мп| {щ, ж) = Ьг}. Обозначим 7Г£(ж) — ортогональная проекция точки ж е Мп на линейное многообразие Ь. Тогда
рь{ ж) = тгь(ж),
то есть, псевдопроекция на линейное многообразие Ь совпадает с ортогональной проекцией. 12 Вестник ЮУрГУ. Серия «Вычислительная математика и информатика»
Доказательство. Зафиксируем произвольную точку Уа & Мп. Рассмотрим линейное многообразие Р, содержащее точку Уа, и являющееся ортогональным дополнением к линейному многообразию Ь:
у0еР = (31)
Обозначим через у точку пересечения линейного многообразия Ь с его ортогональным дополнением Р:
ЬПР = {у}
(см. рис. 3). Построим ортогональную проекцию точки Уо на гиперплоскость Н^ для произвольного ] € 3 '.
Ро = ТяД« о)-
В соответствии с леммой 1 имеем
тгя3 («о) € р-
Отсюда и из (24) следует, что
г>1 = <р(г0) € Р.
Это означает, что последовательность
сходится к точке у пересечения линейного многообразия Ь с линейным многообразием Р, то есть, рь{щ) = V. С другой стороны, в силу (31) имеем 7Гь(щ) = V. Следовательно
Ух е Мп : рь(х) = 1ГЬ(х). Утверждение доказано. □
Процедура приближенного вычисления псевдопроекции на линейное многообразие представлена в виде алгоритма 1. Дадим краткий комментарий по шагам этого алгоритма. На шаге 2 счетчик итераций к устанавливается в значение ноль. Шаг 3 задает начальное приближение жо- Шаг 4 открывает итерационный цикл вычисления псевдопроекции. Шаги 5-8 для текущего приближения хк вычисляют сумму из правой части формулы (24).
Алгоритм 1 Вычисление псевдопроекции pj{x)
Require: Щ = {х £ Мп|(а*,х) = h} ; J С {1,...,m}; J ф 0; f| ф 0
ieJ
1: function pj(x) k:= 0 Жо := ж repeat
£:=0
for г € do
5 := S + ((aj, - bj)cij/||flj||2 end for ®(fc+i) :=a!fc -
imax := 0 > Максимальная невязка
for г € J do
6 := || (ai,ajfc+i) - bi\\ if & > Us then
£max •= Сг
end if end for
until £max < e^ > Малый параметр 6f > 0
return Xk end function
Шаг 9 находит следующее приближение Xk+i■ Шаги 10-16 определяют максимальную невязку £ для нового приближения относительно всех гиперплоскостей Щ, участвующих в вычислении. На шаге 17 счетчик итераций увеличивается на единицу. Шаг 18 проверяет условие выхода из итерационного цикла. Шаг 19 возвращает в качестве результата полученное приближение Xk-
3. Алгоритм движения по граням AlFaMove
В этом разделе мы опишем новый алгоритм AlFaMove (Along Faces Movement), представляющий собой численную реализацию метода поверхностного движения [9]. Алгоритм AlFaMove строит на поверхности допустимого многогранника путь из произвольной граничной точки «о G МПГ(М) до точки х, являющейся решением задачи ЛП (1). Перемещение по граням допустимого многогранника происходит в направлении наибольшего увеличения значения целевой функции. Путь, построенный в результате такого движения, называется оптимальным целевым путем.
В основе алгоритма AlFaMove лежит процедура D(u) вычисления вектора движения d по грани допустимого многогранника М из точки и в направлении максимального увеличения значения целевой функции. Процедура D(u) представлена в виде алгоритма 2. Схематично работа этого алгоритма изображена на рис. 4. Дадим краткий комментарий по шагам алгоритма 2. Шаги 2—7 строят множество U, включающее в себя индексы всех гиперплоскостей Щ, проходящих через точку и. На шаге 8 вектору направления движения d в качестве начального значения присваивается нулевой вектор. На шаге 9 значению це-
Алгоритм 2 Вычисление вектора движения d = D{u) Require: Щ = {ж е Жп\(аих) = h}; и € Г(М)
1 function D(u)
2 U := 0 > U — множество индексов гиперплоскостей Щ, проходящих через точку и
3 for i = 1... m do
4 if (a,i,u) = bi then
5 U:=UU{i}
6 end if
7 end for
8 d:= 0
9 /:= -oo > / — значение целевой функции /(ж) = (с, х)
10 ec-.= c/\\c\\
11 v\=u + Sec > Большой параметр 8 > 0
12 for J € V{U)\% do > Т(Ы) — множество всех подмножеств множества Ы
13 w:=pj(v)
14 d:=w — u
15 ed:=d/\\d\\
16 if (u + red) E M then > Малый параметр г > 0
17 if (c, u + red) > / then
18 / :=(c, u + red)
19 d:=d
20 end if
21 end if
22 end for
23 return d
24 end function
левой функции / присваивается бесконечно малое число3. На шаге 10 строится единичный вектор ес, сонаправленный с вектором с. На шаге 11 строится точка v путем прибавления вектора 5ес к вектору и (см. рис. 4). Здесь 5 — «большой» положительный параметр: чем
3В случае 64-разрядной арифметики с плавающей точкой, бесконечно малое число равно —1 • 10308
больше 5, тем точнее будет вычислен вектор направления движения d. Однако при увеличении параметра S будет увеличиваться и время вычисления псевдопроекции pj(v) на шаге 13. Цикл for на шаге 12 перебирает все возможные комбинации индексов гиперплоскостей, проходящих через точку и. Каждой такой комбинации J соответствует линейное многообразие L = iej Щ, которое также проходит через точку и. Шаг 13 вычисляет точку w, выполняя псевдопроекцию точки v на линейное многообразие, соответствующее комбинации J. На шаге 14 вычисляется возможный вектор движения d по текущему линейному многообразию в направлении максимального увеличения значений целевой функции. Шаг 15 вычисляет единичный вектор е^, сонаправленный с вектором d. Шаг 16 проверяет, что малое движение из точки и в направлении d не выходит за границы допустимого многогранника. Шаг 17, в свою очередь, проверяет, что значение целевой функции в точке (и + еа) превышает максимальное значение, полученное на предыдущих итерациях цикла for (шаг 12). В этом случае это значение запоминается как максимальное (шаг 18), а найденное направление присваивается вектору d (шаг 19). После того, как проверены все возможные комбинации, вектор d возвращается в качестве результата (шаг 23). Если ни одна комбинация не прошла проверку на шагах 16—17, то в качестве результата будет возвращен нулевой вектор. Это означает, что любое движение из точки и по поверхности допустимого многогранника не приводит к увеличению значения целевой функции.
Теперь мы готовы описать алгоритм движения по граням AlFaMove, решающий задачу ЛП (1). За основу мы берем алгоритм 1 из работы [9]. Реализация алгоритма AlFaMove на псевдокоде представлена в виде алгоритма 3. Прокомментируем шаги этого алгоритма. На шаге 1 вводится начальное приближение и. Это может быть произвольная граничная точка рецессивного многогранника М, удовлетворяющая условию
и0 еМПГ(М).
Это условие проверяется на шаге 2. Для получения подходящего начального приближения может применяться алгоритм, реализующий стадию Quest апекс-метода [10]. Шаг 3 вычисляет начальный вектор движения do- Для этого используется функция D(-), реа-
Алгоритм 3 Алгоритм движения по граням AlFaMove
т л
Require: Щ = {х <Е Мп|(а,,ж) ^ fy} ; М = f) Щ; М
¿=1
1: input Uq
2: assert и0€МП Г(М) 3: do:=D(uo) 4: assert do ф О 5: А;:=0 6: repeat
7: ujfe+i :=p(uk,dk) 8: dfc+i :=D(uk+1) 9: к := к + 1 10: until dk = 0 11: Output Uk 12: Stop
= П Н» i6l«(aj,c) > 0
i€l
> Решение задачи ЛП (1)
Рис. 5. Действие функции /л: fi(u,d) = 71 (u,d)
лизованная в алгоритме 2. Предполагается, что do ф О4. Это условие контролируется на шаге 4. На шаге 5 счетчик итераций к устанавливается в значение ноль. Шаг 6 открывает цикл repeat, выполняющий движение по граням допустимого многогранника до тех пор, пока очередной вектор движения dk не окажется нулевым вектором. Это означает, что последнее найденное приближение и^ является решением задачи ЛП (1). Шаг 7 в теле цикла вычисляет следующее приближение iifc+i с помощью вектор-функции ц (определение см. ниже). Шаг 8 вычисляет вектор движения dfc+i для вновь найденного приближения Шаг 9 увеличивает счетчик итераций к на единицу. Если последний вектор движения равен нулевому вектору, то на шаге 10 происходит выход из цикла repeat/until, после чего на шаге 11 выводятся координаты точки и^ в качестве решения задачи ЛП (1). Шаг 12 завершает работу алгоритма AlFaMove.
Вектор-функция /а(-), используемая на шаге 7, определяется следующим образом. Обозначим
Q = {i £ {1,..., т} | (<ц,и) <bi А {щ, d) >0}. (32)
Тогда
/х(и, d) = argmin{||u — ®|| | x = 7¿(u, d)} . (33)
гей
Здесь тi(u, d) обозначает вектор-функцию, вычисляющую косоугольную проекцию точки и на гиперплоскость Н\ по направлению d:
7i{u, d) = и — , btd.
Действие функции ц проиллюстрировано на рис. 5. В соответствии с рисунком гиперплоскости Н^ и не удовлетворяют неравенству (сч,и) < Ь{ в (32). Гиперплоскости Н3 и Н& не удовлетворяют неравенству (a,i,d) > 0 в (32). Таким образом, Q = {1,2}. Так как ||«-7i(«,d)|| < \\u-~f2{u,d)\\, то n{u,d) =7i(tt,d).
4Равенство вектора <1о нулевому вектору означает, что точка ид является решением задачи ЛП (1).
Сходимость алгоритма 3 к решению задачи ЛП (1) за конечное число итераций обеспечивается следующей теоремой.
Теорема 1. (Сходимость алгоритма AlFaMove) Пусть допустимый многогранник М задачи ЛП (1) является ограниченным непустым множеством. Пусть х — решение задачи ЛП (1). Тогда последовательность приближений {uk}k=i, генерируемых алгоритмом 3, является конечной (К < +оо), и (с,ик) = (с,ж), то есть ик является решением задачи ЛП (1).
Доказательство. Обозначим через dAlFaMove вектор dfc+i, вычисляемый на шаге 8 алгоритма 3. В соответствии с шагами 13, 14 алгоритма 2 имеем
dAlFaMove = Pj{v) — U.
Согласно утверждению 2 отсюда следует
dAlFaMove = TTj(v) - U, (34)
где 7Tj(v) обозначает ортогональную проекцию точки v на линейное многообразие L = f]iejHi. В соответствии с утверждением 1 это означает, что вектор dAlFaMove однозначно задает направление максимального увеличения целевой функции f(x) задачи (1). А это, в свою очередь, означает, что алгоритм 3 является численной реализацией метода поверхностного движения [9], то есть последовательности приближений {u\iFaMove} И {измм}' генерируемые алгоритмом 3 из настоящей статьи и алгоритмом 1 из [9] соответственно, совпадают. Таким образом сходимость алгоритма AlFaMove непосредственно следует из сходимости алгоритма поверхностного движения, гарантируемой теоремой 1 из статьи [9]. □
4. Параллельная версия алгоритма AlFaMove
Наиболее ресурсоемкой операцией, используемой в алгоритме AlFaMove (алгоритм 3), является операция вычисления вектора движения, выполняемая в цикле на шаге 8. При решении больших задач ЛП она занимает более 90% процессорного времени. Это объясняется тем, что функция вычисления вектора движения D(-), реализованная в виде алгоритма 2, использует в цикле на шаге 13 операцию псевдопроектирования, заключающуюся в последовательном применении отображения (р(-), задаваемого формулой (24), к исходной точке (см. алгоритм 1, реализующий вычисление псевдопроекции). Известно, что в случае больших задач ЛП проекционный метод может потребовать значительных временных затрат [13]. Кроме того, следует отметить, что алгоритм 2 в цикле на шаге 12 перебирает все непустые подмножества множества U, включающего в себя индексы гиперплоскостей, проходящих через точку и. Если, например, через точку проходит 30 гиперплоскостей, то у нас получится 2зо _ !
= 1073 741823 непустых подмножества. Перебор такого количества подмножеств потребует суперкомпьютерных мощностей. Потому мы разработали параллельную версию алгоритма AlFaMove, представленную в виде алгоритма 4.
Параллельный алгоритм 4 построен на основе модели параллельных вычислений BSF [14], ориентированной на кластерные вычислительные системы. Модель BSF использует схему распараллеливания «мастер—рабочие» и требует представление алгоритма в виде операций над списками с использованием функций высшего порядка Map и Reduce.
Алгоритм 4 Параллельная версия алгоритма AlFaMove
мастер Z-тый рабочий (I = 0,..., L — 1)
1: input n, m, A, b, u0 2: k:=0 3: repeat
4: Broadcast uk
5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15:
16: Gather Creduce
17: (dfc, fk) ■= Reducei©, Creduce)
18: if dfc = 0 then
19: exit := true
20: else
21: uk+1:=fj,(uk,dk)
22: k:=k +1 23: exit := false
24: end if 25: Broadcast exit 26: until exit 27: output Uk,fk 28: Stop
1: input n, m, A, b, c 2:
3: repeat
4: RecvFromMaster uk
5: U:=[]
6: for i = 1... m do 7: if (a,i,uk) = bi then 8: W:=H-H-[t] 9: end if
10: end for
11: K:= 2^-1
12: L := NumberOfWorkers
13: Cmapd) - = [IK/L, . . . , (I + 1 )K/L ~ 1]
14: £>reduce(l) '■= MaP(^uk, ^map(l))
15: (di, fi) := ReduceCreduce(i)) 16: SendToMaster (dj,/i) 17: 18: 19: 20: 21: 22: 23: 24:
25: RecvFromMaster exit
26: until exit
27:
28: Stop
В качестве второго параметра функции высшего порядка Map в алгоритме 4 используется список Стар = [1, • ■ ■, -К], содержащий порядковые номера всех подмножеств множества U, за исключением пустого. Здесь К = 2^1 — 1. В качестве первого параметра фигурирует параметризованная функция
Fu : {1 ->МП х М,
определенная следующим образом:
Fu(j) = (dj, Д);
(35)
U =
—оо,
(с, и + ed), if (u + red) € M Л (с, го) > (с, it); —оо, if (u + red) ^ М V (с, го) ^ (с, г»),
где
го =
pCTt0(w + £c/||c||)
(36)
и
го — и
ed =
w — и
Семантика функции Ри(-) однозначно определяется алгоритмом 2. Функция сг(-), используемая в формуле (36), отображает натуральное число ^ € {1,..., К} в j-тoe подмножество множества элементов списка 1А следующим образом. Число у преобразуется в двоичное представление, состоящее из \1А\ разрядов. Каждому разряду соответствует индекс гиперплоскости из списка 1А в естественном порядке. Если разряд содержит единицу, то соответствующий индекс входит в подмножество сг(з), если — ноль, то не входит. Например, пусть через точку и проходят гиперплоскости Н2, Н^, Н7, Н$. В этом случае 1Л = [2,4,7,9] и К = 24 - 1 = 15, то есть из множества элементов списка Ы может быть сформировано 15 различных непустых подмножеств. Найдем, например, пятое подмножество. Функция сг(-) преобразует число 5 в двоичную запись из 4-х разрядов 0101 и возвращает в качестве результата подмножество {4,9}. Таким образом функция высшего порядка Мар Стар) преобразует список Стар номеров подмножеств в список пар
Здесь dj (j = 1,..., К) является единичным вектором движения, a fj — значение целевой функции, которое достигается в точке и + dj.
Обозначим через Creduce список пар, генерируемый функцией высшего порядка Map:
£reduce = Map (Fu, Стар) = [(dl,fl), {dK, /к)] ■
Определим бинарную ассоциативную операцию
®:Г хКчйпхК,
являющуюся первым параметром функции высшего порядка Reduce:
Функция высшего порядка Reduce редуцирует список пар 4redUce к одной паре путем последовательного применения операции © ко всем элементам списка:
Map (Fu, Стар) = [Fu(l),..., Fu (К)] = [(db Д),..., (dK, fK)\.
(d',/'), если /'>/"; (d",f"), если /'</".
(37)
V
Reduce (©, Creduce) = (db Д) © ... © (dK, fK) = {drJr)
где согласно (37)
Параллельная работа алгоритма 4 организована по схеме «мастер-рабочие» и включает в себя L + 1 процесс: один процесс-мастер и L процессов-рабочих. Процесс-мастер осуществляет общее управление вычислениями, распределяет работу между процессами-рабочими, получает от них результаты и формирует итоговый результат. Для простоты будем предполагать, что количество подмножеств К кратно количеству рабочих L. На шаге 1 мастер вводит исходные данные задачи ЛП и начальную точку «о- На шаге 2 мастер устанавливает счетчик итераций к в значение ноль. Шаги 3-26 реализуют основной цикл repeat/until, вычисляющий решение задачи ЛП (1). На шаге 4 мастер рассылает текущее приближение ик всем рабочим. На шаге 16 он получает от рабочих частичные результаты, которые на шаге 17 редуцируются в пару (dk,fk)■ Если на шаге 18 выполняется условие dk = 0, то решение найдено (мы предполагаем do Ф 0). В этом случае на шаге 19 мастер присваивает логической переменной exit значение true. Если dk ф 0, то мастер на шаге 21 вычисляет следующее приближение Uk+1, увеличивает на шаге 22 счетчик итераций к на единицу и на шаге 23 присваивает логической переменной exit значение false. На шаге 25 мастер рассылает всем рабочим значение логической переменной exit. Если логическая переменная exit принимает значение «истина», цикл repeat/until завершается на шаге 26. На шаге 27 мастер выводит в качестве результата последнее приближение и^ и оптимальное значение целевой функции Д. Шаг 28 завершает работу процесса-мастера.
Все рабочие выполняют один и тот же код, но над различными данными. На шаге 1 Z-тый рабочий (I = 1,..., L) вводит исходные данные задачи ЛП. Цикл repeat/until рабочего (шаги 3-26) соответствует циклу repeat/until мастера. На шаге 4 Z-тый рабочий получает от мастера текущее приближение Uk■ Затем он формирует подсписок £map(i) своих порядковых номеров подмножеств для последующей обработки (шаги 5—13). Подсписки различных рабочих не пересекаются:
I' ф I" тар{1') Ф £'map(l")i (38)
а их конкатенация дает полный список:
Стар = Аттр(О) -Н- • • ■ -Н- Anap(L-l) ■ (39)
На шаге 14 рабочий вызывает функцию высшего порядка Map, которая, в свою очередь, применяет параметризованную функцию FUfc, определенную формулами (35), ко всем элементам подсписка Стар{1) > формируя На выходе ПОДСПИСОК пар £reduce(l) • Этот подсписок на шаге 15 редуцируется рабочим в единственную пару (di,fi) с помощью функции высшего порядка Reduce, которая последовательно применяет бинарную операцию ©, определенную по формуле (37), ко всем элементам подсписка A-e<iuce(Z) • Результат пересылается мастеру на шаге 16. На шаге 25 рабочий получает от мастера значение логической переменной exit. Если эта переменная принимает значение true, то рабочий процесс завершается. В противном случае продолжает выполняться цикл repeat/until. Операторы обмена Broadcast, Gather, RecvFromMaster и SendToMaster обеспечивают неявную синхронизацию работы процесса-мастера и процессов-рабочих.
j' = arg max /,-.
5. Вычислительные эксперименты
Мы реализовали параллельную версию алгоритма AlFaMove на языке С++ с использованием программного BSF-каркаса [15], базирующегося на модели параллельных вычислений BSF [14]. BSF-каркас инкапсулирует все аспекты, связанные с распараллеливанием программы на основе библиотеки MPI. Исходные коды параллельной реализации алгоритма AlFaMove решения задач ЛП свободно доступны в репозитории GitHub по адресу https://github.com/leonid-sokolinsky/AlFcLMove. Разработанная программа была протестирована на большом количестве задач ЛП, взятых из различных источников. Все эти задачи в формате MTX [16] доступны по адресу https://github.com/leonid-sokolinsky/ Set-of-LP-Problems. В качестве тестов мы также использовали искусственные задачи, полученные с помощью генератора случайных задач ЛП FRaGenLP [17]. Эти задачи доступны по адресу https://github.com/leonid-sokolinsky/Set-of-LP-Problems/tree/ main/Rnd-LP. Мы не смогли протестировать реализацию AlFaMove на задачах из репозито-рия Netlib-LP [18], так как во всех этих задачах количество гиперплоскостей, проходящих через начальную точку uq, превышало число 30, что соответствует количеству возможных комбинаций, равному 1073 741824. Массивы таких размеров не допускаются доступными нам компиляторами С++.
С помощью разработанной программы мы исследовали масштабируемость алгоритма AlFaMove. В экспериментах мы использовали параметризованную задачу ЛП «гиперкуб с отсеченной вершиной», для которой размерность пространства п является параметром. Ограничения этой задачи содержат 2n + 1 неравенств следующего вида:
XI
<
Xi +
„ xi ^ О,
Градиент целевой функции
Задача предполагает нахождение максимума целевой функции и имеет единственное решение в точке (100, 200,..., 200) со значением целевой функции, равным 100(n2 + п — 1). Для произвольного п эта задача может быть получена в формате MTX с помощью генератора FRaGenLP, если в качестве количества случайных неравенств задать 0. Эти задачи ЛП для различных п доступны по адресу https://github.com/leonid-sokolinsky/ Set-of-LP-Problems/tree/main/Rnd-LP под именами lp_rnd<n>-0, где в качестве <п> указана размерность пространства. Вычислительные эксперименты проводились на суперкомпьютере «Торнадо ЮУрГУ» [19], характеристики которого представлены в табл. 1. Все вычисления выполнялись с двойной точностью, при которой число с плавающей точкой занимает в оперативной памяти 64 бита.
В первой серии экспериментов исследовалась зависимость ускорения и эффективности распараллеливания алгоритма AlFaMove от количества используемых процессорных узлов кластера при решении задач «гиперкуб с отсеченной вершиной». Результаты этих экспериментов представлены на рис. 6. В данном контексте ускорение a(L) определялось как отно-
sC 200
х2 ^ 200
! (40)
хп ^ 200
х2 ■ ■ ■ + хп ^ 200(п - 1) +100 Х2 ^ 0, • ■ ■ , хп ^ 0.
задается вектором
с=(1,2,...,п). (41)
Таблица 1. Характеристики кластера «Торнадо ЮУрГУ»
Параметр Значение
Количество процессорных узлов Процессоры Число процессоров на узел Память на узел Соединительная сеть Операционная система 480 Intel Xeon Х5680 (6 cores, 3.33 GHz) 2 24 GB DDR3 InfiniBand QDR (40 Gbit/s) Linux CentOS
шение времени Т(1) решения задачи на конфигурации с узлом-мастером и единственным узлом-рабочим ко времени Т(Ь) решения той же задачи на конфигурации с узлом-мастером и Ь узлами-рабочими:
Т( 1)
a(L) =
тщ-
Эффективность распараллеливания /?(!*) вычислялась по формуле
Т(1)
№ =
ь-т(ьу
Вычисления проводились для размерностей 16, 18 и 20. Число ограничений соответственно составило 33, 37 и 41. В качестве начальной точки всегда выбиралась вершина допустимого многогранника со следующими координатами:
xi = 0,
хп/2 = 0, хп/2+х = 200,
хп = 200.
(42)
Эксперименты продемонстрировали хорошую масштабируемость алгоритма А1РаМоуе на задачах «гиперкуб с отсеченной вершиной», начиная с размерности п = 18. В этом случае алгоритм демонстрировал ускорение, близкое к линейному. На меньших размерностях затраты на обмены и латентность начинают превалировать над вычислительными затратами, что приводит к резкому уменьшению границы масштабируемости алгоритма5. Для
Под границей масштабируемости понимается максимальное количество процессорных узлов, на котором наблюдается рост ускорения.
60 100 140 180 Количество процессорных узлов
60 100 140 180 Количество процессорныхузлов
220
Рис. 6. Ускорение и эффективность распараллеливания алгоритма А1РаМоуе
600
500
— 400
CD О
к 300
Z
ф
m 200
100
о
Рис. 7. Зависимость времени работы А1ГаМоуе от размерности задачи на различных многопроцессорных конфигурациях (Ь — количество рабочих узлов)
Размерность задачи
п = 16 эта граница оказалась равной 180 узлам. Эксперименты также показали, что с увеличением размерности задачи эффективность распараллеливания на малом количестве процессорных узлов (меньше 120) падает. Однако, при большем количестве процессорных узлов наблюдается обратная тенденция. Так, для размерности п = 16 эффективность распараллеливания на 20 процессорных узлах составила 61%, после чего снизилась до 23% на 220 узлах. При этом, для п = 20 эффективность распараллеливания оказалась равной 50% и 40% соответственно.
В следующей серии экспериментов исследовалась зависимость времени вычислений от размерности задачи «гиперкуб с отсеченной вершиной» на различных многопроцессорных конфигурациях с количеством рабочих узлов L = 60, L = 120hL = 180. Результаты этих экспериментов представлены на рис. 7. Размерность варьировалась от п = 16 до п = 24 с шагом 2. Для размерности п = 24 каждый из списков Сгтар и С>reduce включал в себя 16 777215 элементов. Это — максимальный размер, допускаемый используемым компилятором. В качестве начальной точки всегда выбиралась вершина допустимого многогранника с координатами (42). На всех исследуемых конфигурациях эксперименты показали экспоненциальный рост времени вычислений с увеличением размерности задачи. Однако, конфигурации с большим количеством рабочих узлов демонстрировали существенно меньшее время работы алгоритма AlFaMove.
В третьей серии экспериментов мы исследовали поведение алгоритма AlFaMove на кубе Кле—Минти (Klee—Minty), представляющем собой параметризованную задачу ЛП с размерностью пространства в качестве параметра. Область допустимых решений этой задачи представляет собой гиперкуб с перекошенными углами и может быть задана следующей системой неравенств:
XI < 5
4Х! + х2 25
8Х! + 4^2 + Хз < 125
2пХ! + 2п~1х2 • • • + 4жп_1 + Хп ^ 5П
Хл > 0, Х2 > 0, i хп ^ 0.
Градиент целевой функции задается вектором
с = (2п-\ 2п~2, ... ,2, 1).
Задача предполагает нахождение максимума целевой функции и имеет единственное решение в точке (0, ... ,0, 5П) со значением целевой функции, равным 5П. Кле и Минти в статье [20] показали, что в случае, когда начальной вершиной является центр координат, классический симплекс-метод при решении этой задачи посещает все 2П вершин, делая 2П — 1 итераций. Известно, что большинство оптимизационных алгоритмов демонстрируют плохую производительность применительно к кубу Кле—Минти. Мы применили алгоритм А1ЕаМоуе для решения кубов Кле—Минти размерности от 5 до 9. Результаты экспериментов, приведенные в табл. 2, показывают, что алгоритм А1РаМоуе во всех случаях находил решение за 2п — 1 итераций, в то время, как симплекс-метод совершал 2" — 1 итераций.
Таблица 2. Эксперименты с кубами Кле—Минти
Размер- AlFaMove Simplex
ность Граница Время Относит. Кол-во
масшт-ти (сек.) погреш-ть итераций Кол-во итераций
5 10 0.2 0.9 • Ю-12 9 31
6 15 2 0.2 • Ю-12 11 63
7 20 13 0.8 • Ю-11 13 127
8 25 126 0.8 • Ю-11 15 255
9 30 1445 0.2 • Ю"10 17 511
Относительная погрешность вычислялась по формуле
fexact fapprox
ô =
f exact
где fexact — точное максимальное значение целевой функции, fapprox — значение, вычисленное алгоритмом AlFaMove. Подсчет итераций симплекс-метода производился с помощью онлайн калькулятора, доступного по адресу https://www.pmcalculators.com/ simplex-method-calculator. Эксперименты также показали, что порог масштабируемости алгоритм AlFaMove на кубах Кле—Минти рос линейно с ростом размерности. При этом одновременно наблюдался экспоненциальный рост времени вычислений.
Заключение
В статье представлен алгоритм AlFaMove, являющийся численной реализацией метода поверхностного движения для решения задач линейного программирования. Ключевой особенностью этого метода является построение оптимального пути на поверхности допустимого многогранника от начальной точки к решению задачи линейного программирования. Под оптимальным путем понимается путь по поверхности допустимой области в
направлении максимального увеличения значений целевой функции. Практическая значимость предложенного метода состоит в том, что он открывает возможность применения искусственных нейронных сетей прямого распространения для решения нестационарных многомерных задач линейного программирования в режиме реального времени.
Теоретической основой алгоритма АШаМоуе является операция построения псевдопроекции на линейные многообразия, которые формируют грани допустимого многогранника. Псевдопроекция реализуется на основе фейеровского процесса и является обобщением понятий ортогональной проекции на линейное многообразие и метрической проекции на выпуклое множество. В случае грани, образуемой гиперплоскостью, псевдопроекция сводится к ортогональной проекции. Доказано, что путь по гиперплоскости, построенный с помощью градиента целевой функции и ортогональной проекции, является оптимальным. Приведен алгоритм проекционного типа для построения псевдопроекции на линейные многообразия меньших размерностей, образуемые пересечением гиперплоскостей. Доказано, что в этом случае псевдопроекция совпадает с ортогональной проекцией. Представлено формализованное описание алгоритма А1РаМоуе, строящего оптимальный путь на поверхности допустимого многогранника. В основе алгоритма А1РаМоуе лежит процедура вычисления вектора движения по грани допустимого многогранника из точки текущего приближения в направлении максимального увеличения значения целевой функции. Дано формализованное описание этой процедуры.
Алгоритмы проекционного типа характеризуются низкой скоростью сходимости, зависящей от углов между гиперплоскостями, образующими линейное многообразие. Также отмечено, что при вычислении вектора движения возникает переборная задача комбинаторного типа, имеющая высокую пространственную и временную сложность. Представлена параллельная версия алгоритма АШаМоуе, ориентированная на кластерные вычислительные системы. Параллельная версия реализована на языке С++ с использованием программного ВЭР-каркаса, основанного на модели параллельных вычислений ВЯР. Проведены эксперименты по исследованию масштабируемости алгоритма А1РаМоуе на кластерной вычислительной системе. Вычислительные эксперименты показали, что задача линейного программирования с 24 переменными и 49 ограничениями демонстрирует ускорение близкое к линейному на 320 процессорных узлах кластера. Задачи большей размерности приводили к ошибке компилятора, связанной с превышением максимального допустимого размера массивов.
В качестве направлений дальнейших исследований выделим следующие. Мы планируем разработать новый, более эффективный метод построения пути на поверхности допустимого многогранника, приводящего к решению задачи линейного программирования. Основная идея состоит в сокращении перебираемых комбинаций гиперплоскостей при определении направления движения. Это может быть достигнуто, если мы ограничимся перемещениями только по ребрам многогранника (линейным многообразиям размерности один). Проблема пространственной сложности может быть решена путем перехода к стохастическим методам выбора направления движения.
Исследование выполнено при финансовой поддержке РНФ (проект №23-21-00356).
Литература
1. Optimization in Large Scale Problems: Industry 4.0 and Society 5.0 Applications / ed. by M. Fathi, M. Khakifirooz, P.M. Pardalos. Cham, Switzerland: Springer, 2019. XI, 340 p. DOI: 10.1007/978-3-030-28565-4.
2. Kopanos G.M., Puigjaner L. Solving Large-Scale Production Scheduling and Planning in the Process Industries. Cham, Switzerland: Springer, 2019. 291 p. DOI: 10.1007/978-3-030-01183-3.
3. Schlenkrich M., Parragh S.N. Solving large scale industrial production scheduling problems with complex constraints: an overview of the state-of-the-art // 4th International Conference on Industry 4.0 and Smart Manufacturing. Procedia Computer Science. Vol. 217 / ed. by F. Longo, M. Affenzeller, A. Padovano, W. Shen. Elsevier, 2023. P. 1028-1037. DOI: 10.1016/J.PR0CS.2022.12.301.
4. Соколинская И.М., Соколинский JI.Б. О решении задачи линейного программирования в эпоху больших данных // Параллельные вычислительные технологии (ПаВТ'2017). Короткие статьи и описания плакатов. Челябинск: Издательский центр ЮУрГУ, 2017. С. 471-484. URL: http://omega.sp.susu.ru/pavt2017/short/014.pdf.
5. Branke J. Optimization in Dynamic Environments // Evolutionary Optimization in Dynamic Environments. Genetic Algorithms and Evolutionary Computation, vol. 3. Boston, MA: Springer, 2002. P. 13-29. DOI: 10.1007/978-1-4615-0911-0_2.
6. Dantzig G.B. Linear programming and extensions. Princeton, N.J.: Princeton university press, 1998. 656 p.
7. Зоркальцев В.И., Мокрый И.В. Алгоритмы внутренних точек в линейной оптимизации // Сибирский журнал индустриальной математики. 2018. Т. 21, 1 (73). С. 11-20. DOI: 10.17377/sibj im. 2018.21.102.
8. Mamalis В., Pantziou G. Advances in the Parallelization of the Simplex Method // Algorithms, Probability, Networks, and Games. Lecture Notes in Computer Science, vol. 9295 / ed. by C. Zaroliagis, G. Pantziou, S. Kontogiannis. Cham: Springer, 2015. P. 281-307. DOI: 10.1007/978-3-319-24024-4_17.
9. Ольховский H.A., Соколинский Л.Б. О новом методе линейного программирования // Вычислительные методы и программирование. 2023. Т. 24, № 4. С. 408-429. DOI: 10. 26089/NumMet.v24r428.
10. Соколинский Л.В., Соколинская И.М. О новой версии апекс-метода для решения задач линейного программирования // Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2023. Т. 12, № 2. С. 5-46. DOI: 10.14529/cmse230201.
11. Мальцев А.И. Основы линейной алгебры. Москва: Наука. Главная редакция физико-математической литературы, 1970. 402 с.
12. Васин В.В., Ерёмин И.И. Операторы и итерационные процессы фейеровского типа. Теория и приложения. Екатеринбург: УрО РАН, 2005. 211 с.
13. Gould N.I. How good are projection methods for convex feasibility problems? // Computational Optimization and Applications. 2008. Vol. 40, no. 1. P. 1-12. DOI: 10.1007/S10589-007-9073-5.
14. Sokolinsky L.B. BSF: A parallel computation model for scalability estimation of iterative numerical algorithms on cluster computing systems // Journal of Parallel and Distributed Computing. 2021. Vol. 149. P. 193-206. DOI: 10.1016/j .jpdc.2020.12.009.
15. Sokolinsky L.B. BSF-skeleton: A Template for Parallelization of Iterative Numerical Algorithms on Cluster Computing Systems // MethodsX. 2021. Vol. 8. Article number 101437. DOI: 10.1016/j .тех.2021.101437.
16. Boisvert R.F., Pozo R., Remington K.A. The Matrix Market Exchange Formats: Initial Design: tech. rep. / NISTIR 5935. National Institute of Standards; Technology. Gaithersburg, MD, 1996. P. 14. URL: https://nvlpubs.nist.gov/nistpubs/Legacy/IR/nistir5935. pdf.
17. Соколинский JI.В., Соколинская И.М. О генерации случайных задач линейного программирования на кластерных вычислительных системах // Вестник Южно-Уральского государственного университета. Серия: Вычислительная математика и информатика. 2021. Т. 10, № 2. С. 38-52. DOI: 10.14529/cmse210103.
18. Gay D.M. Electronic mail distribution of linear programming test problems // Mathematical Programming Society COAL Bulletin. 1985. Vol. 13. P. 10-12.
19. Dolganina N., Ivanova E., Bilenko R., Rekachinsky A. HPC Resources of South Ural State University // Parallel Computational Technologies. PCT 2022. Communications in Computer and Information Science, vol. 1618 / ed. by L. Sokolinsky, M. Zymbler. Cham: Springer, 2022. P. 43-55. DOI: 10.1007/978-3-031-11623-0_4.
20. Klee V., Minty G.J. How good is the simplex algorithm? // Inequalities - III. Proceedings of the Third Symposium on Inequalities Held at the University of California, Los Angeles, Sept. 1-9, 1969 / ed. by O. Shisha. New York-London: Academic Press, 1972. P. 159-175.
Соколинский Леонид Борисович, д.ф.-м.н., профессор, заведующий кафедрой системного программирования, Южно-Уральский государственный университет (национальный исследовательский университет) (Челябинск, Российская Федерация)
Ольховский Николай Александрович, аспирант, кафедра системного программирования, Южно-Уральский государственный университет (национальный исследовательский университет) (Челябинск, Российская Федерация)
Соколинская Ирина Михайловна, к.ф.-м.н., доцент, кафедра вычислительной математики и высокопроизводительных вычислений, Южно-Уральский государственный университет (национальный исследовательский университет) (Челябинск, Российская Федерация)
JT.B. Cokojihhckhh, H.A. Ojibxobckhh, H.M. Cokojihhckelh
DOI: 10.14529/ cmse240301
NUMERICAL IMPLEMENTATION OF SURFACE MOVEMENT METHOD IN LINEAR PROGRAMMING
© 2024 L.B. Sokolinsky, N.A. Olkhovsky, I.M. Sokolinskaya
South Ural State University (pr. Lenina 76, Chelyabinsk, 454080 Russia) E-mail: [email protected], [email protected], [email protected]
Received: 10.06.2024
The article is devoted to the numerical implementation of a new linear programming method called the surface movement method. The implementation is based on the AlFaMove algorithm, which builds, on the surface of the feasible polytope, the optimal objective path from an arbitrary boundary point to a point that is a solution to the linear programming problem. The optimality of the path lies in the fact that the direction of movement along the face of the polytope corresponds to the maximum increase in the value of the objective function. To calculate the optimal direction of movement, a method based on the operation of pseudoprojection onto a linear manifold is used. The pseudoprojection operation is a generalization of the notion of orthogonal projection. Pseudo projection is implemented using an iterative projection-type algorithm. The proposition is proved that the pseudoprojection coincides with the orthogonal projection in the case of a linear manifold that is the intersection of hyperplanes. It is also proved that, in the case of a linear manifold, the pseudoprojection method calculates a movement vector that coincides with the direction of maximum increase of the objective function. A parallel implementation of the AlFaMove algorithm is performed. The results of computational experiments on a cluster computing system are presented, which demonstrate the high scalability of the proposed numerical implementation.
Keywords: linear programming, surface movement method, numerical implementation, AlFaMove algorithm, parallel implementation, cluster computing system, scalability evaluation.
FOR CITATION
Sokolinsky L.B., Olkhovsky N.A., Sokolinskaya I.M. Numerical Implementation of Surface Movement Method in Linear Programming. Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering. 2024. Vol. 13, no. 3. P. 5-31. (in Russian) DOI: 10.14529/ cmse240301.
This paper is distributed under the terms of the Creative Commons Attribution-Non Commercial 4-0 License which permits non-commercial use, reproduction and distribution of the work without further permission provided the original work is properly cited.
References
1. Optimization in Large Scale Problems: Industry 4.0 and Society 5.0 Applications / ed. by M. Fathi, M. Khakifirooz, P.M. Pardalos. Cham, Switzerland: Springer, 2019. XI, 340 p. DOI: 10.1007/978-3-030-28565-4.
2. Kopanos G.M., Puigjaner L. Solving Large-Scale Production Scheduling and Planning in the Process Industries. Cham, Switzerland: Springer, 2019. 291 p. DOI: 10.1007/978-3-030-01183-3.
3. Schlenkrich M., Parragh S.N. Solving large scale industrial production scheduling problems with complex constraints: an overview of the state-of-the-art. 4th International Conference on Industry 4.0 and Smart Manufacturing. Procedia Computer Science. Vol. 217 / ed. by F. Longo, M. Affenzeller, A. Padovano, W. Shen. Elsevier, 2023. P. 1028-1037. DOI: 10.1016/J.PR0CS.2022.12.301.
4. Sokolinskaya I.M., Sokolinsky L.B. On the Solution of Linear Programming Problems in the Age of Big Data. Parallel Computational Technologies. PCT 2017. Communications in Computer and Information Science, vol. 753. / ed. by L. Sokolinsky, M. Zymbler. Cham, Switzerland: Springer, 2017. P. 86-100. DOI: 10.1007/978-3-319-67035-5_7.
5. Branke J. Optimization in Dynamic Environments. Evolutionary Optimization in Dynamic Environments. Genetic Algorithms and Evolutionary Computation, vol. 3. Boston, MA: Springer, 2002. P. 13-29. DOI: 10.1007/978-l-4615-0911-0_2.
6. Dantzig G.B. Linear programming and extensions. Princeton, N.J.: Princeton university press, 1998. 656 p.
7. Zorkaltsev V., Mokryi I. Interior point algorithms in linear optimization. Journal of applied and industrial mathematics. 2018. Vol. 12, no. 1. P. 191-199. DOI: 10. 1134/ S1990478918010179.
8. Mamalis B., Pantziou G. Advances in the Parallelization of the Simplex Method. Algorithms, Probability, Networks, and Games. Lecture Notes in Computer Science, vol. 9295 / ed. by C. Zaroliagis, G. Pantziou, S. Kontogiannis. Cham: Springer, 2015. P. 281-307. DOI: 10.1007/978-3-319-24024-4_17.
9. Olkhovsky N.A., Sokolinsky L.B. Surface Movement Method for Linear Programming. 2024. DOI: 10.48550/arXiv.2404.12640. arXiv: 2404.12640.
10. Sokolinsky L.B., Sokolinskaya I.M. Apex Method: A New Scalable Iterative Method for Linear Programming. Mathematics. 2023. Vol. 11, no. 7. P. 1654. DOI: 10.3390/MATH11071654.
11. Maltsev A. The basics of linear algebra. Moskow: Science. The main editorial office of the phys-math literature, 1970. 402 p. (in Russian).
12. Vasin V.V., Eremin I.I. Operators and Iterative Processes of Fejer Type. Theory and Applications. Berlin, New York: Walter de Gruyter, 2009. 155 p. Inverse and Ill-Posed Problems Series. DOI: 10.1515/9783110218190.
13. Gould N.I. How good are projection methods for convex feasibility problems?. Computational Optimization and Applications. 2008. Vol. 40, no. 1. P. 1-12. DOI: 10.1007/S10589-007-9073-5.
14. Sokolinsky L.B. BSF: A parallel computation model for scalability estimation of iterative numerical algorithms on cluster computing systems. Journal of Parallel and Distributed Computing. 2021. Vol. 149. P. 193-206. DOI: 10.1016/j .jpdc.2020.12.009.
15. Sokolinsky L.B. BSF-skeleton: A Template for Parallelization of Iterative Numerical Algorithms on Cluster Computing Systems. MethodsX. 2021. Vol. 8. Article number 101437. DOI: 10.1016/j.mex.2021.101437.
16. Boisvert R.F., Pozo R., Remington K.A. The Matrix Market Exchange Formats: Initial Design: tech. rep. / NISTIR 5935. National Institute of Standards; Technology. Gaithersburg, MD, 1996. P. 14. URL: https://nvlpubs.nist.gov/nistpubs/Legacy/IR/nistir5935. pdf.
17. Sokolinsky L.B., Sokolinskaya I.M. FRaGenLP: A Generator of Random Linear Programming Problems for Cluster Computing Systems. Parallel Computational Technologies. PCT
2021. Communications in Computer and Information Science, vol. 1437 / ed. by L. Sokolinsky, M. Zymbler. Cham: Springer, 2021. P. 164-177. DOI: 10.1007/978-3-030-81691-9_12.
18. Gay D.M. Electronic mail distribution of linear programming test problems. Mathematical Programming Society COAL Bulletin. 1985. Vol. 13. P. 10-12.
19. Dolganina N., Ivanova E., Bilenko R., Rekachinsky A. HPC Resources of South Ural State University. Parallel Computational Technologies. PCT 2022. Communications in Computer and Information Science, vol. 1618 / ed. by L. Sokolinsky, M. Zymbler. Cham: Springer,
2022. P. 43-55. DOI: 10.1007/978-3-031-11623-0_4.
20. Klee V., Minty G.J. How good is the simplex algorithm?. Inequalities - III. Proceedings of the Third Symposium on Inequalities Held at the University of California, Los Angeles, Sept. 1-9, 1969 / ed. by O. Shisha. New York-London: Academic Press, 1972. P. 159-175.
