Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org Решение нестационарных задач линейного программирования большой размерности на кластерных вычислительных системах И.М. Соколинская, Л.Б. Соколинский Южно-Уральский государственный университет В работе описывается параллельный алгоритм решения нестационарных задач линейного программирования большой размерности, ориентированный на кла- стерные вычислительные системы. В основе алгоритма, получившего название «следящий», лежат фейеровские отображения. Алгоритм отслеживает измене- ния исходных данных и вносит корректировки в вычислительный процесс. При этом задача разбивается на большое количество подзадач, которые могут ре- шаться независимо без обменов данными. Приводятся диаграммы деятельности UML, описывающие параллельный следящий алгоритм. 1. Введение Одной из особенностей современного экономико-математического моделирования явля- ется появление нестационарных задач линейного программирования [1] большой размер- ности с быстро меняющимися входными данными. Кроме этого, часть ограничений мо- жет оказаться плохо формализуемой [2]. Одним из примеров таких задач является задача управления портфелем ценных бумаг с использованием методов алгоритмической торгов- ли [3, 4]. В подобных задачах количество переменных и неравенств в системе ограничений может составлять десятки и даже сотни тысяч, а период изменения исходных данных нахо- диться в пределах сотых долей секунды. В работе [5] был описан параллельный алгоритм для решения задач линейного программирования с плохо формализуемыми ограничени- ями. Для преодоления проблемы нестационарности входных данных в работах [6, 7] был предложен «следящий» алгоритм решения задачи линейного программирования с исполь- зованием фейеровских отображений [8], ориентированный на кластерные вычислительные системы с многоядерными ускорителями. В данной работе дается параллельная реализация следящего алгоритма с использованием диаграммы деятельности UML. Статья организова- на следующим образом. В разделе 2 приводится формальная постановка задачи линейного программирования, даются определения фейеровского процесса и операции псевдопроек- тирования на многогранник. В разделе 3 приводится неформальное описание следящего алгоритма. В разделе 4 описывается межузловое распараллеливание следящего алгоритма, предполагающее использование технологии параллельного программирования MPI. В за- ключении суммируются полученные результаты и определяются направления дальнейших исследований. 2. Формализация задачи Пусть задана задача линейного программирования max \{ \langle c, x\rangle | Ax \leq b, x \geq 0\} . (1) Определим отображение \varphi : \BbbR n \rightarrow \BbbR n следующим образом: m \sum max \{ \langle aj , x\rangle - bj , 0\} \varphi (x) = x - \alpha j \lambda j \cdot aj . (2) j=1 \| aj \| 2 420 Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org Пусть M – многогранник, задаваемый ограничениями задачи линейного программирова- ния (1). Такой многогранник всегда является выпуклым. Известно [8], что \varphi будет одно- значным непрерывным M -фейеровским отображением для любой системы положительных m коэффициентов \{ \alpha j > 0\} (j = 1, . . . , m), таких, что \alpha j = 1, и коэффициентов релаксации \sum j=1 0 < \lambda j < 2. Полагая в формуле (2) \lambda j = \lambda и \alpha j = 1/m (j = 1, . . . , m), получаем формулу m \lambda \sum max \{ \langle aj , x\rangle - bj , 0\} \varphi (x) = x - \cdot aj , (3) m j=1 \| aj \| 2 которая будет использоваться в следящем алгоритме. Обозначим \varphi s (x) = \varphi . . . \varphi (x). \underbrace{} \underbrace{} s Под фейеровским процессом, порождаемым отображением \varphi при произвольном начальном приближении x0 \in \BbbR n , будем понимать последовательность \{ \varphi s (x0 )\} +\infty s=0 . Известно, что ука- занный фейеровский процесс сходится к точке, принадлежащей множеству M: \{ \varphi s (x0 )\} +\infty s=0 \rightarrow x \= \in M. Будем кратко обозначать это следующим образом: lim \varphi s (x0 ) = x \=. s\rightarrow \infty Под \varphi -проектированием (псевдопроектированием ) точки x \in \BbbR n на многогранник M \varphi \varphi понимается отображение \pi M : \BbbR n \rightarrow M [9], задаваемое соотношением \pi M (x) = lim \varphi s (x). s\rightarrow \infty 3. Следящий алгоритм В общем виде следящий алгоритм может быть описан следующим образом. Все про- странство \BbbR n делится гиперкубической сеткой на ячейки с переменной длиной грани. Раз- мер ячейки может динамически меняться в зависимости от скорости изменения исходных данных: чем быстрее меняются данные, тем больше становится ячейка, чтобы успевать сле- дить. Алгоритм постоянно находит (отслеживает) ячейку наименьшего размера, на которой достигается максимум целевой функции. Будем такую ячейку называть целевой. В каждой итерации алгоритма просчитываются ячейки, входящие в некоторую следящую область вокруг целевой ячейки. В простейшем случае следящая область может быть кубической формы со следящей ячейкой в центре. В общем случае следящая область может быть про- извольной выпуклой фигурой. При этом ячейки всегда имеют кубическую форму. Целевая точка – это некоторая точка, находящаяся вне следящей области. Ее координаты могут динамически вычисляться по коэффициентам целевой функции. Например, при целевой функции Q (x) = x1 + 2x2 в качестве целевой точки можно взять точку (T, 2T ), где T – достаточно большое положительное число. Нулевой вершиной кубической области будем называть вершину куба, находящуюся ближе всего к началу координат. Нулевой вершиной кубической ячейки будем называть вершину ячейки, находящуюся ближе всего к началу координат. Неформально следящий алгоритм может быть описан следующей последовательностью действий. 1. Первоначально выбирается следящая кубическая область с длиной ребра r и нулевой вершиной в точке q , заведомо покрывающая многогранник. 2. Выбирается целевая точка z = T c, расположенная вне следящей области. 3. Следящая область разбивается на K ячеек. 421 Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org 4. В условиях динамического изменения входных данных (A, b, c), для всех ячеек внутри следящей области вычисляется псевдопроекция из точки z на пересечение ячейки с многогранником. Если пересечение пусто, то такие ячейки отбрасываются. 5. Если получено пустое множество псевдопроекций, то мы увеличиваем размер следя- щей области в w раз и переходим на шаг 2. 6. Если получено непустое множество псевдопроекций, то на нем вычисляется максимум целевой функции. 7. Если расстояние от точки максимума до центра следящей области меньше 41 r, то длина ребра следящей области r уменьшается в 2 раза. 8. Если расстояние от точки максимума до центра следящей области больше 34 r, то длина ребра следящей области r увеличивается в 2 раза. 9. Центр следящей области перемещается в центр ячейки, в которой достигнут макси- мум, и переходим на шаг 2. Псевдопроекции на шаге 4 для различных ячеек следящей области могут вычисляться параллельно без обменов данными между MPI-процессами. Фейеровский процесс вычис- ления псевдопроекции, стартовавший из целевой точки на шаге 4, вычисляется до конца, несмотря на то, что координаты целевой точки и сам многогранник в процессе счета могли измениться. Выполнение условия шага 5 означает, что входные данные изменяются быстрее, чем мы вычисляем псевдопроекции. Константы 1/2 и 3/4 , используемые при выполнении шагов 7 и 8, являются значениями параметров алгоритма. 4. Межузловое распараллеливание Межузловое распараллеливание – это распараллеливание вычислений между отдель- ными узлами кластерной вычислительной системы. В качестве технологии параллельного программирования используется MPI. Каждый MPI-процесс работает на отдельном узле кластера. Обмен данными между MPI-процессами осуществляется с помощью передачи со- общений по коммуникационной сети. Один MPI-процесс считает одну псевдопроекцию на пересечение одной ячейки следящей области с многогранником. Таким образом, количе- ство ячеек следящей области всегда равно константе K, совпадающей с количеством MPI- процессов. Опишем метод межузлового распараллеливания следящего алгоритма. Без ограниче- ния общности мы можем считать, что все процессы происходят в положительной области координат. Пусть P – количество доступных MPI-процессов (количество процессорных уз- лов в многопроцессорной системе). Определим общее количество ячеек в следящей области следующим образом: \Bigl( \Bigl\lfloor \Bigr\rfloor \Bigr) n K = P 1/n . (4) Тогда количество ячеек, примыкающих к одному ребру куба следящей области, равно h = K 1/n . Зададим в пространстве целочисленных координат u0 , . . . , un линейную нумера- цию ячеек следящей области следующим образом. Пусть ячейка \alpha имеет целочисленные координаты (\alpha 0 , . . . , \alpha n ). Тогда ее номер k\alpha вычисляется по формуле: n - 1 \sum k\alpha = \alpha i ni . (5) i=0 422 Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org Рис. 1. Линейная нумерация ячеек следящей области при n = 3. На рис. 1 приведен пример такой линейной нумерации при n = 3. Например, ячейка с номером 19 на рис. 1 имеет целочисленные координаты (1, 0, 2). Действительно, 19 = 1 \cdot 30 + 0 \cdot 31 + 2 \cdot 32 . Выразим целочисленные координаты ячейки \alpha через ее порядковый номер k\alpha . Из (5) получаем \alpha 0 = k\alpha mod n; (6) k\alpha - \alpha 0 \alpha 1 = mod n; (7) n k\alpha - \alpha 0 - \alpha 1 n \alpha 2 = mod n; (8) n2 ......... Таким образом, в общем виде имеем: i - 1 \alpha j nj \sum k\alpha - j=0 \alpha i = mod n. (9) ni Формула (9) содержит ресурсоемкую операцию возведения в степень. От нее можно изба- виться следующим образом. По определению операции mod из (6) получаем \alpha 0 = k\alpha - (k\alpha \div n) \cdot n. (10) С помощью символа \div здесь обозначается целочисленное деление. Подставив в (7) вместо \alpha 0 правую часть этого уравнения, получим \alpha 1 = k\alpha - (k\alpha - (k n \alpha \div n)\cdot n) mod n = (k\alpha \div n)\cdot n mod n (11) n = (k\alpha \div n) mod n. По определению операции mod отсюда следует \alpha 1 = k\alpha \div n - ((k\alpha \div n) \div n) \cdot n. (12) Подставив в (8) вместо \alpha 0 правую часть уравнения (10), а вместо \alpha 1 – правую часть урав- нения (12), получим \alpha 2 = k\alpha - \alpha n02 - \alpha 1 n mod n = k\alpha - (k\alpha - (k\alpha \div n)\cdot n) - (kn\alpha 2 \div n - ((k\alpha \div n)\div n)\cdot n)\cdot n mod n (13) = ((k\alpha \div n) \div n) mod n. 423 Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org Из (11) и (13) для i = 1, . . . , n - 1 по индукции получаем: \alpha i = (((k\alpha \div n) . . .) \div n) modn. (14) \underbrace{} \underbrace{} i Пусть g = (g0 , . . . , gn - 1 ) – нулевая вершина куба следящей области; y = (y0 , . . . , yn - 1 ) – нулевая вершина произвольной ячейки \alpha . Выразим координаты точки y через координаты точки g . Обозначим s = K 1/n r – шаг сетки. Тогда yi = gi + s\alpha i (15) для i = 0, . . . , n - 1. Определим в качестве центральной ячейки куба ячейку \gamma с целочисленными координа- тами (\gamma 0 , . . . , \gamma n - 1 ), где \Bigl\lfloor \Big/ \Bigr\rfloor \gamma 0 = . . . = \gamma n - 1 = K 1/n 2 (16) Пусть q = (q0 , . . . , qn - 1 ) – нулевая вершина центральной ячейки \gamma . Выразим координаты точки q через координаты точки g , используя формулу (15): qi = gi + s\gamma i , (17) для i = 0, . . . , n - 1. Пусть y – нулевая вершина ячейки \alpha .Тогда область внутри ячейки \alpha (включая границы) задается системой из 2n неравенств: \left\{ - x0 \leq - y0 \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot - xn - 1 \leq - yn - 1 (18) x0 \leq y0 + s \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot xn - 1 \leq yn - 1 + s Чтобы найти пересечение ячейки \alpha с многогранником M необходимо к системе неравенств Ax \leq b добавить неравенства из системы (18). При этом получается система из m + 2n неравенств с n неизвестными. Параллельная версия следящего алгоритма изображена на рис. 2 в виде диаграммы деятельности UML. В качестве начального значения длины ребра r кубической следящей области берется значение R, а в качестве нулевой вершины g кубической следящей области берется точка G, которые заведомо обеспечивают покрытие многогранника M следящей областью. Начальный шаг сетки s определяется по формуле s = K 1/n r , где количество ячеек K следящей области вычисляется по формуле (4). Координаты целевой точки z получаются путем умножения вектора c коэффициентов целевой функции на масштабирующий коэф- фициент T , выбираемый таким образом, чтобы целевая точка находилась вне следящей области. Целочисленные координаты центральной ячейки следящего куба определяется по формуле (16). Нулевая вершина q центральной ячейки следящей области соответственно вычисляется по формуле (17). Вектор-функция \pi (z, k) с помощью фейеровского процесса, описанного в разделе 2, вычисляет псевдопроекцию целевой точки z на пересечение многогранника с ячейкой, име- ющей порядковый номер k , и возвращает точку псевдопроекции xk или вектор ( - 1, . . .), 424 Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org Рис. 2. Параллельная версия следящего алгоритма. 425 Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org если точка псевдопроекции не принадлежит многограннику (пересечение многогранника M с ячейкой пусто). Псевдопроекции вычисляются параллельно без обменов данными в цикле until. В хо- де вычисления псевдопроекций исходные данные задачи (1) могут меняться. После того как все точки псевдопроекций x0 , . . . , xK - 1 вычислены, в цикле for происходит вычисление порядкового номера k\alpha ячейки \alpha , на которой достигается максимум целевой функции в точке псевдопроекции. При этом в переменной C вычисляется максимум целевой функции. В соответствии с этим в качестве начального значения k\alpha выбирается наименьшее целое число в машинном представлении MinInt, а в качестве начального значения C выбирается наименьшее вещественное число в машинном представлении MinFloat. После завершения цикла for проверяется значение k\alpha . Если оно равно MinInt, зна- чит пересечение всех ячеек с многогранником M оказалось пустым. Это означает, что в результате изменения исходных данных задачи (1) многогранник M «уплыл» за пределы следящей области. В этом случае в w раз увеличиваются: длина r ребра следящей области, шаг сетки s и координаты целевой точки z . Константа w является параметром алгоритма. Если значение k\alpha отлично от MinInt, значит пересечение многогранника M со следящей областью не пусто. В этом случае в качестве новой центральной ячейки рассматривается ячейка с номером k\alpha и вычисляются координаты q \prime – нулевой вершины этой ячейки. Для этого используется вектор-функция zero(k\alpha ), диаграмма деятельности которой приведена на рис. 3. Рис. 3. Вектор-функция, вычисляющая нулевую вершину ячейки с порядковым номером k\alpha . Далее рассматриваются три случая в зависимости от удаленности нулевой вершины новой центральной ячейки следящей области от старой. Если \| q - q \prime \| > 34 r, то r, s и ко- ординаты z увеличиваются в полтора раза. Если \| q - q \prime \| < 41 r, то r, s и координаты z уменьшаются в полтора раза. Если 14 r \leq \| q - q \prime \| \leq 43 r то r, s и координаты z остаются прежними. После этого происходит сдвиг нулевой вершины следящей области на вектор (q \prime - q) и копирование координат q \prime в q . Далее выполняется очередная итерация цикла until. 5. Заключение В работе была описана параллельная реализация следящего алгоритм для решения нестационарных задач большой размерности на кластерных вычислительных системах. Ал- горитм был реализован на языке Си с использованием библиотеки MPI. Начальные экспери- менты на модельной задаче показали почти линейную масштабируемость параллельной ре- ализации следящего алгоритма на 480 узлах суперкомпьютера «Торнадо ЮУрГУ». Авторы 426 Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org планируют реализовать внутриузловое распараллеливание следящего алгоритма с исполь- зованием технологии параллельного программирования OpenMP. Объектом распараллели- вания на этом уровне является отдельный фейеровский процесс, вычисляющий псевдопро- екцию. Для каждого подвектора организуется отдельная нить. Вычисления предполагается выполнять с использованием ускорителей Xeon Phi, которыми оснащен вычислительный кластер «Торнадо ЮУрГУ». Эффективное использование многоядерных ускорителей до- стигается путем распараллеливания процесса вычисления отдельной псевдопроекции. При этом предполагается использовать метод разбиения вектора на подвекторы [9]. Суть мето- да заключается в том, что вектор исходной точки делится на подвекторы по числу про- цессорных ядер многоядерного ускорителя. На каждом подвекторе независимо делается v фейеровских приближений с использованием редуцированных фейеровских отображений. Такое редуцированное отображение воздействует только на «свой» подвектор, оставляя оставшуюся часть вектора без изменений. Затем нити управления через общую память об- мениваются модифицированными подвекторами и процесс продолжается. В работе [9] была доказана сходимость описанного итерационного процесса. Литература 1. Еремин И.И., Мазуров Вл.Д. Нестационарные процессы математического программирования. М.: Наука, 1979. 291 с. 2. Мазуров Вл.Д. Распознавание образов как метод формирования плохо формализуемых ограничений в моделях планирования // Оптимизация. Вып. 10(27). Новосибирск: СО АН СССР, 1973. С. 144-158. 3. Дышаев М.М., Соколинская И.М. Представление торговых сигналов на основе адаптивной скользящей средней Кауфмана в виде системы линейных неравенств // Вестник Южно-Уральского государственного университета. Серия: Вычислительная математика и информатика. 2013. Т. 2. № 4. С. 103-108. 4. Ананченко И.В., Мусаев А.А. Торговые роботы и управление в хаотических средах: обзор и критический анализ // Труды СПИИРАН. 2014. № 3 (34). С. 178-203. 5. Sokolinskaya I.M., Sokolinskii L.B. Parallel algorithm for solving linear programming problem under conditions of incomplete data // Automation and Remote Control. 2010. Vol. 71, No. 7. P. 1452-1460. 6. Соколинская И.М., Соколинский Л.Б. О применении фейеровских отображений в задачах линейной оптимизации с быстро меняющимися входными данными // Информационный бюллетень Ассоциации программирования. № 13. ИММ УрО РАН, 2015. C. 56-58. 7. Соколинская И.М., Соколинский Л.Б. Алгоритм решения нестационарных задач линейного программирования для кластерных вычислительных систем с многоядерными ускорителями // Параллельные вычислительные технологии (ПаВТ’2015): труды международной научной конференции. Челябинск: Издательский центр ЮУрГУ, 2015. С. 477-481. 8. Еремин И.И. Фейеровские методы для задач выпуклой и линейной оптимизации. Челябинск: Изд-во ЮУрГУ, 2009. 200 с. 9. Ершова А.В., Соколинская И.М. О сходимости масштабируемого алгоритма построения псевдопроекции на выпуклое замкнутое множество // Вестник Южно-Уральского государственного университета. Серия: Математическое моделирование и программирование. 2012. № 18 (277). С. 5-12. 427 Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org Solving unstable linear programming problems of high dimension on cluster computing systems Irina Sokolinskaya and Leonid Sokolinsky Keywords: linear programming, unstable problems of high dimension, cluster computing systems, Fejer mappings Solving unstable linear programming problems of high dimension on cluster computing systems with multi-core accelerators