=Paper=
{{Paper
|id=Vol-1576/185
|storemode=property
|title=Implementation of Parallel Pursuit Algorithm for Solving Unstable Linear Programming Problems
|pdfUrl=https://ceur-ws.org/Vol-1576/185.pdf
|volume=Vol-1576
|authors=Irina Sokolinskaya,Leonid Sokolinsky
}}
==Implementation of Parallel Pursuit Algorithm for Solving Unstable Linear Programming Problems==
https://ceur-ws.org/Vol-1576/185.pdf
Параллельные вычислительные технологии (ПаВТ’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]-ε