Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Параллельная реализация следящего алгоритма для решения нестационарных задач линейного программирования И.М. Соколинская, Л.Б. Соколинский Южно-Уральский государственный университет В работе описывается параллельный алгоритм решения нестационарных задач ли- нейного программирования большой размерности, ориентированный на кластерные вычислительные системы. В основе алгоритма, получившего название «следящий», лежат фейеровские отображения. Алгоритм отслеживает изменения исходных дан- ных и вносит корректировки в вычислительный процесс. При этом задача разбивает- ся на большое количество подзадач, которые могут решаться независимо без обменов данными. Приводятся диаграммы деятельности UML, описывающие реализацию следящего алгоритма. Ключевые слова: нестационарная задача линейного программирования, фейеровские отображения, следящий алгоритм, диаграммы деятельности UML, массовый парал- лелизм, кластерные вычислительные системы. 1. Введение Основы современной теории нестационарных процессов математического программирова- ния были заложены в классической монографии [1], где было предложено использовать для решения нестационарных задач линейного программирования итерационные процессы фейе- ровского типа. Указанный подход обобщает релаксационный метод Моцкина-Агмона [2, 3]. Это направление исследований получило продолжение в дальнейших работах И.И. Еремина, В.В. Васина, Л.Д. Попова, Е.А. Бердниковой, С.В. Пацко и других ученых [4]. Нестационарные задачи линейного программирования большой размерности с быстро ме- няющимися входными данными достаточно часто встречаются в практике современного эко- номико-математического моделирования. Одним из примеров таких задач является задача управления портфелем ценных бумаг с использованием методов алгоритмической торговли [5, 6]. В подобных задачах количество переменных и неравенств в системе ограничений может составлять десятки и даже сотни тысяч, а период изменения исходных данных находиться в пределах сотых долей секунды. На нестационарность в таких задачах может накладываться плохая формализуемость части ограничений. В работе авторов [7] был описан параллельный алгоритм для решения задач линейного программирования с плохо формализуемыми ограни- чениями. Суть предложенного подхода заключается в синтезе методов линейного программи- рования и дискриминантного анализа. Для выполнения эффективного дискриминантного ана- лиза необходимы два набора образцов M и N , первый из которых удовлетворяет неформали- зованным ограничениям, а второй – нет. Для получения качественных наборов образцов могут применяться методы интеллектуального анализа [8] данных и анализа временных рядов [9]. Для преодоления проблемы нестационарности входных данных в работах [10, 11] был предложен «следящий» алгоритм решения задачи линейного программирования с использова- нием фейеровских отображений, ориентированный на кластерные вычислительные системы с многоядерными ускорителями. В данной работе дается полная параллельная реализация сле- дящего алгоритма с использованием диаграммы деятельности UML. Статья организована сле- дующим образом. В разделе 2 приводится формальная постановка задачи линейного програм- мирования, даются определения фейеровского процесса и операции псевдопроектирования на многогранник. В разделе 3 приводится математическое описание следящей области. В разделе 4 приводятся формулы для построения пересечения многогранника, задаваемого системой ограничений, с произвольной ячейкой следящей области. В разделе 4 дается полное описание параллельной реализации следящего алгоритма с помощью диаграмм деятельности UML. В заключении суммируются полученные результаты и определяются направления дальнейших исследований. 685 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 2. Постановка задачи Пусть задана задача линейного программирования max { c, x | Ax ≤ b, x ≥ 0 } (1) Определим фейеровское отображение ϕ : ℝn → ℝn следующим образом: m max { ai , x − bi , 0 } ϕ ( x ) = x − ∑ αiλi ⋅ ai (2) 2 i =1 ai Пусть M – многогранник, задаваемый ограничениями задачи линейного программирования (1). Такой многогранник всегда является выпуклым. Известно [12], что ϕ будет однозначным не- прерывным M-фейеровским отображением для любой системы положительных коэффициентов m { αi > 0 } , i = 1, …, m , таких, что ∑ αi = 1 , и коэффициентов релаксации 0 < λi < 2 . Пола- i =1 гая в формуле (2) λi = λ и αi = 1 m (i = 1, …, m ) , получаем формулу λ max { ai , x − bi , 0 } m ϕ (x ) = x − ∑ ⋅ ai , (3) m i =1 2 a i которая будет использоваться в следящем алгоритме. Обозначим ϕs (x ) = ϕ …ϕ(x ) . s Под фейеровским процессом, порождаемым отображением ϕ при произвольном начальном +∞ приближении x 0 ∈ ℝn , будем понимать последовательность { ϕs (x 0 )} . Известно, что ука- s =0 занный фейеровский процесс сходится к точке, принадлежащей множеству M: +∞ {ϕs (x 0 )}s =0 → x ∈ M . (4) Будем кратко обозначать это следующим образом: lim ϕs (x 0 ) = x . s →∞ Под ϕ -проектированием (псевдопроектированием) точки x ∈ ℝn на многогранник M по- ϕ нимается отображение πM (x ) = lim ϕs (x ) . s →∞ 3. Построение следящей области Без ограничения общности мы можем считать, что все процессы происходят в положитель- ной области координат. Пусть n – размерность пространства решений, K – количество ячеек в следящей области по одному измерению. Пусть P – количество MPI-процессов, используемых для распараллеливания вычислений. Будем предполагать, что всегда выполняется равенство: Kn = P , (5) то есть, количество ячеек следящей области равно количеству MPI-процессов. Зададим в про- странстве целочисленных координат u0 , …, un линейную нумерацию ячеек следящей области следующим образом. Пусть ячейка α имеет целочисленные координаты (α0 , …, αn ) . Тогда ее номер kα вычисляется по формуле: n −1 kα = ∑ αin i . (6) i =0 686 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt u1 u2 24 25 26 15 16 17 21 22 23 12 13 14 18 19 20 6 7 8 3 4 5 9 10 11 0 1 2 u0 Рис. 1. Линейная нумерация ячеек следящей области при n = 3. На рис. 1 приведен пример такой линейной нумерации при n = 3 . Например, ячейка с номе- ром 19 на рис. 1 имеет целочисленные координаты (1, 0, 2). Действительно, 19 = 1 ⋅ 30 + 0 ⋅ 31 + 2 ⋅ 32 . Выразим целочисленные координаты ячейки α через ее порядковый номер kα . Из (6) по- лучаем α0 = kα mod n ; (7) kα − α0 α1 = mod n ; (8) n k − α0 − α1n α2 = α mod n ; (9) n2 ...... Таким образом, в общем виде имеем: i −1 k α − ∑ αj n j j =0 αi = mod n . i (10) n Формула (10) содержит ресурсоемкую операцию возведения в степень. От нее можно избавить- ся следующим образом. По определению операции mod из (7) получаем α0 = kα − (kα ÷ n ) ⋅ n 1. (11) Подставив в (8) вместо α0 правую часть этого уравнения, получим kα − (kα − (k α ÷ n ) ⋅ n ) α1 = mod n n (k ÷ n ) ⋅ n = α mod n (12) n = ( k α ÷ n ) mod n. По определению операции mod отсюда следует α1 = kα ÷ n − ((kα ÷ n ) ÷ n ) ⋅ n . (13) Подставив в (9) вместо α0 правую часть уравнения (11), а вместо α1 – правую часть уравнения (13), получим 1 С помощью символа ÷ здесь обозначается целочисленное деление. 687 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt kα − α0 − α1n α2 = mod n n2 kα − (kα − (kα ÷ n ) ⋅ n ) − (kα ÷ n − ((kα ÷ n ) ÷ n ) ⋅ n ) ⋅ n = mod n (14) n2 = ( (kα ÷ n ) ÷ n ) mod n. Из (12) и (14) для i = 1, …, n − 1 по индукции получаем: αi = ( ( ( k α ÷ n )…) ÷ n ) mod n . (15) i Пусть g = (g 0 , …, gn −1 ) – нулевая вершина куба следящей области. Пусть y = (y0 , …, yn −1 ) – нулевая вершина произвольной ячейки α . Выразим координаты точки y r через координаты точки g. Обозначим s = – шаг сетки. Тогда K yi = gi + s αi (i = 0, …, n − 1) . (16) Определим в качестве центральной ячейки куба ячейку γ с целочисленными координата- ми (γ 0 , …, γn −1 ) , где γ 0 = … = γn −1 =  K 2  . (17)   Пусть q = (q 0 , …, qn −1 ) – нулевая вершина центральной ячейки γ . Выразим координаты точки q через координаты точки g, используя формулу (16): qi = gi + s γi (i = 0, …, n − 1) . (18) 4. Пересечение многогранника M с ячейкой Пусть y – нулевая вершина ячейки . Тогда область внутри ячейки (включая границы) задается системой из 2n неравенств:  −x 0 ≤ −y0   −x1 ≤ −y1   ⋯ ⋯ ⋯   −xn −1 ≤ −yn −1  (19)  x 0 ≤ y0 + s   x1 ≤ y1 + s  ⋯ ⋯ ⋯   x n −1 ≤ yn −1 + s  Эта же система в матричной форме: Aαx ≤ bα , (20) где (для n = 3 ) 688 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt  −1 0 0   −y0       0 −1 0   −y      1   0 0 −1   −y2  Aα =   , b =    1 0 0  α  y + s  . (21)    0   0 1 0   y + s   1   0 0 1   y + s      2 Положим A  b  A′ =   , b ′ =   . (22)  Aα   bα  Тогда пересечение многогранника M с ячейкой задается системой неравенств в матричной форме A′x ≤ b ′ , (23) где A′ – расширенная матрица размера (m + 2n ) × n , b ′ – расширенный столбец свободных членов. Расширенный столбец b ′ в соответствии с формулой (22) имеет инвариантную часть b, не зависящую от координат нулевой вершины ячейки , и вариативную часть bα , зависящую от координат нулевой вершины ячейки . Элементы расширенной матрицы A′ не зависят от координат нулевой вершины ячейки . 5. Реализация следящего алгоритма В данном разделе описывается полная реализация следящего алгоритма в виде диаграмм деятельности UML. 5.1 Схема головной подпрограммы Общая схема головной подпрограммы следящего алгоритма приведена на рис. 2. На шаге 1 выполняется подпрограмма init (см. раздел 5.2), выполняющая инициализацию переменных. Затем в цикле until с меткой 2 выполняется корректировка следящей области в соответствии с описанием идеи алгоритма в [11]. Одна итерация соответствует одной корректировке. Головная подпрограмма следящего алгоритма оформляется в виде независимого процесса, который вы- полняется до тех пор, пока переменная stop не примет значение true (истина). Начальную установку переменной stop в значение false (ложь) осуществляет головной процесс, соответ- ствующий основной программе. Он же присваивает переменной stop значение true , когда вы- числительный процесс нужно остановить. В качестве текущего приближения решения задачи (1) головная программа использует текущее значение нулевой вершины центральной ячейки q , координаты которой вычисляются по формуле (18). В теле цикла until выполняются следующие действия. На шаге 3 организуется K парал- лельных потоков управления (нитей), которые независимо друг от друга вычисляют псевдопро- екции из целевой точки z на пересечение i-той ячейки с многогранником M (i = 0, …, P − 1) . Напомним, что P равно количеству MPI-процессов, и в соответствии с формулой (5) равно количеству ячеек в кубической следящей области. Схема подпрограммы вычисления псевдо- проекции детально описана в разделе 5.3. В цикле for с меткой 5 для полученных на шаге 3 точек псевдопроекций x 0 , …, x P −1 вы- числяется номер kα ячейки, на которой достигается максимум C целевой функции. Для кор- ректной работы цикла 5 переменным kα и C на шаге 4 присваиваются начальные значения 689 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt MinInt и MinFloat соответственно. Значение MinInt соответствует минимальному машин- ному значению целого типа, MinFloat – минимальному машинному значению вещественного типа. Подпрограмма, вычисляющая точку x k = (x k 0 , …, x k ,n −1 ) псевдопроекции точки z на пе- ресечение многогранника M с ячейкой с номером kα , присваивает координате xk 0 значение −1 в том случае, когда точка x k псевдопроекции не принадлежит многограннику M . Эта си- туация возникает в случае, когда пересечение многогранника M с ячейкой с номером kα явля- ется пустым. 1 2 0..n-1 z[0..n-1], ⋅⋅⋅ 0..n-1 z[0..n-1],P 3 k =MinInt; C=MinFloat 4 5 6 σ = xk[0]=-1 σ C C = σ; k =k 7 k =MinInt 8 0..n-1 k 9 ‖q-q'‖>¾r ‖q-q'‖<¼r s= s s= s 10 ¼r≤‖q-q'‖≤¾ r 11 z[0..n-1] z[0..n-1] z[0..n-1] z[0..n-1] 0..n-1 0..n-1 0..n-1 0..n-1 12 0..n-1 0..n-1 13 z[0..n-1] z[0..n-1] Рис. 2. Головная подпрограмма следящего алгоритма. 690 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 0..n-1 k τ=k ÷ n [0] = [0] + (k − τ∗n)∗s [i]= [i]+(τ mod n)∗s τ=τ ÷ n Рис. 3. Схема подпрограммы zero вычисления нулевой вершины ячейки с порядковым номером k . Если же x k принадлежит многограннику, то в силу предположения о том, что все процессы находятся в положительной области координат (см. раздел 3), значение xk 0 не может быть от- рицательным. Указанное условие проверяется на шаге 6. Случаи xk 0 = −1 из рассмотрения исключаются. Если при выполнении цикла 5 оказывается, что ни одна из ячеек следящей обла- сти не имеет непустого пересечения с многогранником M , то в переменной kα сохраняется значение MinInt . Этот факт проверяется на шаге 7. В этом случае шаг сетки s , длина r ребра следящей области и координаты целевой точки z увеличиваются в w раз, где w – положи- тельная константа, являющаяся параметром алгоритма (шаг 13 на рис. 2). Если на шаге 7 выясняется, что kα ≠ MinInt , значит найдена ячейка с номером kα , име- ющая непустое пересечение с многогранником, на которой достигается максимум целевой функции. В этом случае на шаге 8 вычисляется вектор q ′ , представляющий нулевую вершину новой центральной ячейки следящей области. Схема подпрограммы вычисления нулевой вер- шины ячейки с порядковым номером k приведена на рис. 3. Вычисления осуществляются с использованием формул (11), (15) и (16). На шаге 9 (рис. 2) анализируется, насколько новая центральная ячейка далеко отстоит от 3 предыдущей. Если расстояние между новой и старой центральными ячейками превышает r 4 (где r – длина ребра кубической следящей области), то длина ребра кубической следящей об- ласти r , шаг сетки s и координаты целевой точки z на шаге 10 увеличиваются в 1.5 раза. Если 1 расстояние между новой и старой центральными ячейками меньше r , то длина ребра кубиче- 4 ской следящей области r , шаг сетки s и координаты целевой точки z на шаге 11 уменьшаются в 2 раза. Величина 1 2 , используемая в шагах 10 и 11, в общем случае является параметром 1 3 алгоритма. Если же отклонение лежит в пределах от r до r , то корректировка значений r , 4 4 s и z не производится. Величины 1 4 и 3 4 также являются в общем случае параметрами ал- горитма. На шаге 12 следящая область сдвигается по вектору (q ′ − q ) , и в качестве текущей нулевой вершины q центральной ячейки следящей области берется точка q ′ . 691 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt input n, m, R, p, K, u, L, T 1 2 n mod u = 0 3 input A[0..n-1,0..m-1], b[0..m-1], c[0..n-1], G[0..n-1] =rank 4 5 Kn=P [0..m-1]= [0..m-1] 6 Aα[0..2n-1,0..n-1]=0 7 8 Aα[i,i]=-1 Aα[i+n,i]=1 [0..m-1,0..n-1]=A[0..m-1,0..n-1] 9 [m..m+2n-1,0..n-1]=Aα[0..2n-1,0..n-1] 10 r=R; s=r/K 11 g[0..n-1]=G[0..n-1] 12 z[0..n-1]=c[0..n-1]∗T 13 0..n-1 [0..n-1]+s∗γ[0..n-1] 14 γ[0..n-1]=(6K/27,...,6K/27) 15 Рис. 4. Схема подпрограммы инициализации переменных init. 5.2 Схема подпрограммы инициализации переменных init Подпрограмма инициализации переменных init выполняет ввод исходных данных и ини- циализацию переменных. Схема подпрограммы init приведена на рис. 4. На шаге 1 вводятся значения переменных: n – размерность пространства решений; m – число неравенств в системе ограничений; R – начальное значение длины ребра следящей области, обеспечивающее покры- 692 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt тие многогранника M; p – количество итераций при построении псевдопроекции, выполняемое между обновлениями входных данных (этот параметр используется в подпрограмме dataChange обновления исходных данных); K – количество ячеек в следящей области по од- ному измерению; u – размерность подвектора; L – число независимых фейеровских итераций на подвекторах в подпрограмме вычисления псевдопроекции (см. рис. 5); T – масштабирующий коэффициент для вычисления координат целевой точки z. На шаге 2 проверяется условие (assert) n mod u = 0 , означающее, что размерность пространства решений n кратно размерно- сти подвектора u . На шаге 3 осуществляется ввод исходных данных задачи линейного про- граммирования: A – матрица коэффициентов неравенств; b – столбец свободных членов; c – вектор коэффициентов целевой функции. Также здесь вводится вектор G, содержащий началь- ные координаты нулевой вершины кубической следящей области. Шаг 4 присваивает перемен- ной P значение, равное количеству доступных MPI-процессов и вычисляемое с помощью си- стемной функции rank . На шаге 5 проверяется условие (assert) K n = P , означающее, что об- щее количество ячеек следящей области равно количеству MPI-процессов. На шаге 6 формиру- ется инвариантная часть расширенного столбца b ′ свободных членов, определяемого по фор- муле (22). На шаге 7 происходит инициализация матрицы Aα путем присваивания всем ее эле- ментам нулевых значений. На шаге 8 определяются ненулевые элементы матрицы Aα в соот- ветствии с системой неравенств (19). На шагах 9 и 10 строится расширенная матрица A′ , опре- деляемая формулой (22). На шаге 11 в качестве начального значения длины r ребра следящей области определяется значение R, обеспечивающее покрытие многогранника M, и вычисляется значение s длины ребра ячейки. На шаге 12 в качестве начальной нулевой вершины g кубиче- ской следящей области берется точка G, определяющая такое положение следящей области, при котором она полностью покрывает многогранник M. На шаге 13 вычисляются координаты целевой точки z по формуле z = Tc . На шаге 14 вычисляется нулевая вершина q центральной ячейки следящей области по формуле (18). На последнем шаге вычисляется вектор γ началь- ных целочисленных координат центральной ячейки по формуле (17). 5.3 Схема подпрограммы вычисления псевдопроекции На рис. 5 приведена схема подпрограммы вычисления псевдопроекции x = π(z, kα ) из це- левой точки z на пересечение многогранника M с ячейкой следящей области, имеющей по- рядковый номер kα , где kα вычисляется по формуле (6). Псевдопроекция вычисляется путем организации фейеровского процесса (4) (см. раздел 2). На шаге 1 выполняется инициализация переменных, необходимых для организации итерационного процесса. В качестве начального значения x k берется точка z ; с помощью подпрограммы zero (см. рис. 3) вычисляется нулевая вершина y ячейки с номером kα ; по формуле (21) определяется вариативная часть bα расши- ренного столбца b ′ системы ограничений (23), получаемой при пересечении многогранника M с ячейкой α . В цикле 2 вычисляется normsq – вектор квадратов норм строк расширенной мат- 2 рицы A′ : normsqi = ai′ . На шаге 3 организуется итерационный процесс вычисления псевдопроекции. Для обеспе- чения высокой масштабируемости процедуры вычисления псевдопроекции используется метод разбиения вектора x на h подвекторов размерности u , предложенный в работе [13]. Мы здесь предполагаем, что n = h ⋅ u . На каждом v-том подвекторе делается L независимых итераций вида (x vu , …, x(v +1)u −1 ) := (x vu , …, x(v +1)u −1 ) − λ m max { (ai,vu , …, ai,(v +1)u −1 ),(xvu , …, x(v +1)u −1 ) − bi , 0 } − ∑ m i =1 2 ⋅ (ai ,vu , …, ai ,(v +1)u −1 ). a i 693 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 0..n-1 z[0..n-1],k 1 x[0..n-1]=z[0..n-1]; [0..n-1] = zero(k ); [0..n-1]; [0..n-1]+s 2 = [i,j]* [i,j] 3 x'[0..u-1]=x[0..u-1] x'[ ]=x[ ] scal_ai_x=0 scal_ai_x=0 scal_ai_x+=A'[i,j]∗x[j] ... scal_ai_x+=A'[i,j]∗x[j] factor=max(scal_ai_x-b[i],0)/ factor=max(scal_ai_x-b[i],0)/ S[j]+=factor∗A'[i,j] S[j]+=factor∗A'[i,j] x[0..u-1]=x'[0..u-1]-(λ/m)∗ x[ ..n]=x'[ ..n]-(λ/m)∗ |x'[0..n-1]-x[0..n-1]|>ε 4 [0..n-1] k ) Рис. 5. Схема подпрограммы π вычисления псевдопроекции. 694 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 0..n-1 k q[i]-ε