Параллельный алгоритм целочисленного квадратичного программирования

Тип работы:
Реферат
Предмет:
Физико-математические науки


Узнать стоимость

Детальная информация о работе

Выдержка из работы

Вычислительные технологии
Том 9, № 1, 2004
ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ ЦЕЛОЧИСЛЕННОГО КВАДРАТИЧНОГО ПРОГРАММИРОВАНИЯ*
Г. И. Злвиняко, Е. А. Котельников Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: zabin@rav. sscc. ru
The parallel algorithm of integer and mixed-integer quadratic programming, based on the branch and bound method. The algorithm was realized in FORTRAN using the MPI system of parallel programming. The efficiency of the parallel and the sequential algorithms are compared for test problems.
Введение
Задача целочисленного и частично целочисленного квадратичного программирования (ЦКП) представляется в следующем виде:
минимизировать f (ж) = ^(х, Ох) + (с, х) (1)
при условиях
Ах = Ь, (2)
а & lt- х & lt- в, (3)
х^ - целые числа для ]? 3. (4)
Здесь, А — матрица размера т х п- О — симметричная матрица размера I х I (I & lt- п) — О & gt- 0- с, х, а, в? Яп- Ь € Кт- 3 — список целочисленных переменных.
Алгоритм решения задач ЦКП основывается на методе ветвей и границ с односторонним ветвлением, который представляет собой обобщение алгоритма целочисленного линейного программирования (ЦЛП). Решение исходной задачи сводится к решению последовательности оценочных задач квадратичного программирования.
Распараллеливание осуществляется асинхронным исполнением на каждом из процессоров алгоритма ветвей и границ с односторонним ветвлением. В результате получается
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 01−07−90 367).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.
некоторое сочетание алгоритмов с односторонним и одновременным ветвлением [1]. Использование MPI существенно упрощает программирование такого сочетания алгоритмов. Во время исполнения процессы обмениваются короткими сообщениями на фоне выполнения большого числа арифметических операций, что обеспечивает высокую эффективность распараллеливания как на вычислительных системах с общей памятью, так и на кластерах.
1. Решение оценочных задач
Оценочные задачи отличаются от исходной (1)-(4) тем, что в них отсутствует требование целочисленности (4) и в каждой г-й оценочной задаче условия (3) заменяются условиями аг & lt- х & lt- вг. Векторы аг и вг формируются по правилам метода ветвей и границ, которые рассматриваются далее.
Для решения оценочных задач используется метод приведенного градиента [2], в котором переменные х подразделяются на базисные хв, супербазисные х^ и небазисные хN. В матрице А, соответственно, выделяются В-базисная т х т, Б-супербазисная т х в и N-небазисная т х (п — (т + в)) матрицы. Предполагается, что В — невырожденная матрица, а переменные XN принимают свои граничные значения.
В текущем подпространстве супербазисных переменных минимизация проводится методом сопряженных градиентов. На к-м шаге метода последовательно определяются:
1) вектор двойственных переменных ук = (В-1)ТV/(хкв) —
2) вектор приведенного градиента Ьк = V/(х|) — БТук-
3) вектор направления в подпространстве супербазисных переменных
кк, 1И1 2 «к-1.
— hk I 11 112
ps — -h + Tnjk-uviPs
Il hk ||2 l|hk-1 Il?'-
4) вектор направления в подпространстве базисных переменных рВ = -В 1Бр|-
г ^ к (к, ря) к к к
5) оптимальный шаг вдоль р, А к = -, к ~ к,, где вектор р составлен из рВ, р? и
(р ,)
pN = 0. Полагается х к+1 = х к + Ар к, где, А = шт{А к, Атах}, а Атах = argmax{А: аг & lt- хк+
Л
Арк & lt- вг}.
Если, А к & gt- Атах, то по методу приведенного градиента меняется состав супербазисных и других переменных, процесс по методу сопряженных градиентов начинается в новом подпространстве. В противном случае процесс по методу сопряженных градиентов продолжается до достижения с заданной точностью равенства нулю приведенного градиента.
Для уменьшения влияния ошибок вычислений на процесс по методу сопряженных градиентов осуществляется контроль точности решения систем уравнений ВТу к = V/в (хкв), ВрВ = Бр| и, если необходимо, их решения уточняются.
Против зацикливания метода сопряженных градиентов в плохо обусловленных задачах используется прием из [3], который состоит в проверке через определенное число итераций выполнения неравенства ||кк||2 & lt- с (к)||кк — кк||2, где величина с (к) возрастает с увеличением к, кк — приведенный градиент с непосредственным вычислением V/(х|), а в кк используется рекуррентно вычисленный V/(х|). При выполнении неравенства процесс по методу сопряженных градиентов в данном супербазисном подпространстве завершается.
После завершения оптимизации методом сопряженных градиентов в текущем супербазисном подпространстве на основе двойственных оценок небазисных столбцов проверяется
возможность сформировать новый набор супербазисных переменных. Если кандидатов на перевод из небазисных переменных в супербазисные нет, то оценочная задача решена. В противном случае используется метод сопряженных градиентов в новом супербазисном подпространстве.
2. Алгоритм целочисленного квадратичного программирования
Оценки границ целевой функции (1) с учетом (4) получаются из решений оценочных задач квадратичного программирования. Выбор переменной ветвления в ЦЛП часто осуществляют с использованием штрафов [1, 2]. В рассматриваемом алгоритме ЦКП также делается попытка, когда это возможно, для выбора переменной ветвления использовать штрафные оценки.
Вначале рассмотрим линейный случай. Значение базисной переменной Xj для ]? 3 в оптимальном базисе очередной оценочной задачи представим в виде Xj = [х,] + Vj, где [xj] - целая часть х,. Штраф за увеличение х, на величину 1 — V, обозначим через Р+, а за уменьшение на величину Vj — через Р-.
На примере Р+ покажем алгоритм определения штрафных оценок. Номера небазисных переменных зададим в двух списках: в I- поместим номера переменных, фиксированных на нижней границе, а в 1+ - на верхней границе. Если переменная х, занимает в базисе р-ю позицию, то вычислим г = е^В-1, ер — р-й орт в Кт. Тогда
P+ = min & lt- min |(1 — Vj)(-min |(1 — Vj)(- q
I q: q€/-, apq& lt-0 ^ q: q? l+, apq& gt-0 ^ ^ a-
¦pq
где apq = (z, Aq), Aq — q-й небазисный столбец матрицы A- dq — его приведенная оценка. Если q-я небазисная переменная целочисленная, то величины (1 — Vj)(dq/ - apq) или (1 — Vj)(-dq/apq), которые меньше |dq |, заменяются на |dq |. Начальные значения Pj- и Pj+ полагаются равными
Штрафы для квадратичной функции f заменим штрафами для линейной функции g (x) = (Vf (ж), ж — ж*), где x* - решение очередной оценочной задачи квадратичного программирования. Поскольку исходная функция выпуклая, штрафные оценки для g (x) будут не больше штрафов для f (ж). В нелинейном случае при вычислении штрафов необходимо дополнительно учесть супербазисные переменные (если они присутствуют в оптимальном решении оценочной задачи).
Величины dq для супербазисных переменных равны нулю (dq — q-я компонента приведенного градиента), и штрафы при рассмотрении очередной j-й переменной не обнуляются, если apq = 0. Для любой j-й супербазисной переменной Pj+ = Pj- = 0.
Пусть Pmax = max max {Pj-, P+}. В задачах ЦЛП Pmax = 0 получается в оценочных
задачах, которые имеют не единственный оптимальный базис. В нелинейном случае возможность получить Pmax = 0 увеличивается за счет супербазисных переменных. Однако численные эксперименты показали, что использование штрафов в отдельных задачах ЦКП позволяет сократить время решения в несколько раз.
В случае Pmax = 0 для ветвления выбирается базисная или супербазисная переменная Xj, которая имеет значение, наиболее уклоняющееся от целочисленного. При этом только для определения переменной ветвления временно полагаются P+ = 1 — Vj и Pj- = Vj.
Схемы метода ветвей и границ с односторонним ветвлением позволяют использовать компактную форму списков оценочных задач. Пусть на некотором уровне к для ветвления выбрана переменная х, со штрафами Р- и Р+. Если Р- & lt- Р+, то в качестве очередной
^ и и и и
оценочной задачи выбирается задача, соответствующая Р-, а в списке оценочных задач (во вспомогательном массиве к) необходимо запомнить информацию об альтернативной задаче. Для этого производятся присвоения: к (1, к) = к (2, к) = в, где в, — верхняя граница переменной х, в г-й оценочной задаче на предшествующем уровне (к-1) — к (3, к) = 1, если /г + Р+ & gt- гг и к (3,к) = 0 в противном случае, здесь /г — оптимальное значение / в г-й оценочной задаче, а гг — текущее значение рекорда- к (4,к) = /г- к (5, к) = Р+. При Р+ & lt- Р, выбирается оценочная задача, отвечающая штрафу Р+, а в массиве к запоминаются величины: к (1, к) = -к (2, к) = а, к (3, к), равная 0 либо 1 в зависимости от выполнения неравенства /г + Р- & gt- гг, к (4, к) = /г, к (5, к) = Р-.
В дальнейшем при рассмотрении параллельного алгоритма ветвь, которой соответствует к (3,к) = 1, будем называть помеченной и не помеченной при к (3, к) = 0.
Алгоритм.
Шаг 0. Положить г = 0, к = 0, значение рекорда г0 = а0 = а и в0 = в.
Шаг 1. Решить текущую задачу квадратичного программирования. Пусть хг — ее оптимальное решение и /г = / (хг).
А. Если /г & gt- гг или система ограничений несовместна, то присвоить гг+1 = гг, г = г + 1 и перейти на шаг 3.
Б. Пусть /г & lt- гг. Если вектор хг нецелочисленный, то перейти на шаг 2. Если решение хг целочисленное, то присвоить гг+1 = /г, г = г + 1. При /г = /0 перейти на шаг 5, иначе на шаг 3.
Шаг 2. Для тех базисных и супербазисных переменных х, ]? 3, у которых х, — нецелое, вычислить штрафы Р+, Р,-. Среди них найти минимальный штраф Рт-п:
А. Если /г + Рт-п & gt- гг, то перейти к шагу 3.
Б. При /г + Рт-п & lt- гг осуществить ветвление по базисной или супербазисной переменной х, которой соответствует максимальный штраф Р+ или Р,-. Из двух величин Р+, Р- выбрать меньшую. Пусть Р- & lt- Р+, тогда выбрать очередную оценочную задачу, соответствующую Р -. Присвоить к = к + 1, записать в список задачу, отвечающую штрафу Р+. Изменить верхнюю границу переменной х, положив ее равной [х, ]. Перейти на шаг 1.
(Если Р+ & lt- Р-, то в списке к запоминается задача, отвечающая Р-, а в формируемой задаче нижняя граница переменной х, полагается равной [х,] + 1.)
Шаг 3. Если к = 0, то перейти к шагу 5.
Если к & gt- 0 и к (3, к) = 1 или к (4, к) + к (5, к) & gt- гг, то перейти к шагу 4.
Иначе, используя массив к, сформировать задачу, альтернативную образованной в шаге 2 Б. Присвоить к (3, к) = 1 и перейти на шаг 1.
Шаг 4. Установить двусторонние ограничения на переменную х, используя значения к (1, к), к (2, к). Присвоить к = к — 1 и перейти на шаг 3.
Шаг 5. Остановиться.
Схемы ветвей и границ с односторонним ветвлением позволяют с небольшими затратами осуществить переход к очередной оценочной задаче. Рассмотренная реализация схемы во многом сходна с реализацией для ЦЛП в [4].
3. Параллельный алгоритм
Параллельный алгоритм устроен так, что на каждом из процессорных элементов исполняется рассмотренный ранее алгоритм с небольшими изменениями. Среди процессорных элементов выделяется элемент с нулевым номером, который кроме исполнения алгоритма ветвей и границ выполняет некоторые диспетчерские функции. По ходу решения на нулевом процессоре накапливается справочная информация, необходимая для организации параллельных вычислений.
На начальном этапе данные о задаче считываются с дисковой памяти нулевым процессором и пересылаются всем остальным. Далее на нулевом процессоре решается нулевая оценочная задача (исходная задача без учета условий целочисленности), а остальные процессоры находятся в ожидании.
Для загрузки любого процессора в его оперативную память пересылаются следующие данные: целочисленные массивы 1 В, 1и — соответственно списки базисных, супербазисных и небазисных переменных- часть массива Н (часть списка оценочных задач) — х^ - массив значений супербазисных переменных. Значения элементов массивов 1 В, 1я, 1и, х^ и какая часть массива Н пересылается, определяется уровнем к, с которого начинается решение первой после загрузки оценочной задачи. На основе списка 1 В строится матрица, обратная к базисной, и начинается исполнение алгоритма ветвей и границ с уровня к.
Для упорядочения загрузки процессорных элементов на нулевом процессоре хранится массив пар (г, кг), где г — номер процессора, а кг — уровень непомеченной ветви на г-м процессоре. В списке пар для каждого процессора г сохраняется кг с наименьшим значением (что соответствует наиболее высокой непомеченной ветви на данном процессоре). Если некоторый г-й процессор простаивает, то это помечается присвоением кг = -1.
Рассмотрим случай, когда г-й процессор завершил выполнение алгоритма (достигнут нулевой уровень) или процессор изначально еще не загружен. Если номер процессорного элемента г & gt- 0, то нулевому процессору г-й процессор посылает запрос на загрузку. После получения запроса на нулевом процессоре выбирается пара (], kj) с минимальным значением kj & gt- 0 и на ]-й процессор отправляется сообщение о необходимости загрузки г-го процессора. В результате на ]-м процессоре ветвь уровня kj помечается и отправляются данные для загрузки на г-й процессор. Если не существует kj & gt- 0, то г-й процесс находится в ожидании. При завершении исполнения алгоритма на нулевом процессоре аналогично определяется номер ]-го процессора для загрузки.
Если на некотором процессоре получено новое значение рекорда г, то это значение пересылается всем остальным процессорам. При получении значения рекорда на каждом из процессоров проверяется выполнение неравенства Н (4, к) — Н (5, к) & gt- г, где к — уровень, с которого начиналось исполнение алгоритма на данном процессоре. Если неравенство выполняется на каком-либо из процессоров, то на этом процессоре завершается исполнение алгоритма ветвей и границ и формируется запрос на загрузку.
Полученное значение г используется для обновления информации о помеченных и непомеченных ветвях в массиве Н. Может оказаться, что на некотором из процессоров ] помечается не помеченная ветвь с минимальным значением kj. В этом случае на нулевой процессор посылается информация об изменении kj.
После решения очередной оценочной задачи на каждом из процессоров производится обращение к функции МР1_1РИ, ОВЕ [5] для проверки, имеются ли сообщения для данного процесса. Кроме того, обращение к этой функции производится после выполнения определенного числа итераций решения оценочной задачи. Частота обращений к функции за-
дается в таблице управляющих параметров, где заданы и другие значения или начальные значения управляющих параметров, такие как частота перепостроения матриц, обратных к базисным, допустимый уровень ведущих элементов и др.
Задача решена, если на всех процессорных элементах достигнут нулевой уровень.
Передача исходных данных о задаче производится с помощью блокированной функции MPI & quot-один — всем& quot-. В дальнейшем все процессы выполняются в асинхронном режиме и обмены осуществляются с помощью неблокированных функций.
4. Решение задач
Тестовые задачи взяты из [6]. В табл. 1 указаны входные параметры задач.
В табл. 2 приведены результаты расчетов в однопроцессорном режиме на системе МВС 1000 М с процессорами альфа-21 264, с тактовой частотой 830 МГц.
Таблица 1
Входные параметры задач
ь Название
задачи задачи m me n ni nb nz
1 ibell3a 106 — 122 60 31 304
2 ibienstl 576 128 505 28 28 2184
3 ield76 75 75 1898 1898 1898 19 111
4 imas284 68 — 151 150 150 9631
5 imisc07 212 35 260 259 259 8619
6 imodOll 4480 4367 10 957 97 96 22 253
7 iqiu 1192 132 840 48 48 3432
8 iran13 195 26 338 169 169 676
9 iran8 296 40 512 256 256 1024
Примечание. т — общее число ограничений, среди которых те ограничений являются равенствами- п — число переменных, среди которых Пг являются целочисленными и щ переменных — булевыми- пг — число ненулей в матрице А.
Таблица 2
Результаты расчетов в однопроцессорном режиме на системе МВС 1000 М
ь задачи it cp t Pi P2 Po
1 354 065 35 765 66. 76 2680 9648 18
2 32 317 905 71 883 45 315. 26 2500 30 294
3 42 003 240 91 423 51 143. 150 4385 10 577
4 1 054 772 53 052 380.8 409 16 672 127
5 3 848 924 54 826 3057. 202 6797 2086
6 5 500 000 5411 83 258. 15 193 —
7 32 000 000 150 246 85 163. 368 16 556 4166
8 18 496 522 1 324 478 4480. 8847 537 387 6604
9 16 239 520 985 215 5751. 4385 400 440 15 115
Примечание. И — суммарное число итераций в оценочных задачах- ср — количество оценочных задач- Ь — время решения, Р1 указывает, сколько раз в процессе решения задачи выполнилось условие fг + Рт1п & gt- гг- Р2 — число случаев выполнения условия fг + Р, & gt- гг, где Р, = Р+ или
Pj
Р- - Р0 — указывает, сколько раз в процессе решения задачи было получено Ртах = 0.
В задачах 6 и 7 в однопроцессорном режиме не удалось получить точное решение из-за слишком больших затрат машинного времени. Задачи снимались со счета после выполнения заданного количества итераций. Обозначим через f полученное приближение по функ-
f — f *
ции, через f * - оптимальное значение и определим относительную ошибку f
If *l
Значения величины 8f для задач 6 и 7 соответственно равны 0. 119 и 0. 002.
Характеристики процесса решения задач параллельным алгоритмом содержит табл. 3. В параллельном варианте расчетов использовались восемь процессоров. Во всех вариантах расчетов значения управляющих парметров выбирались одни и те же, кроме параметра, определяющего частоту обращения к функции MPI_IPROBE. Этот параметр для разных задач принимался в пределах от 10 до 30 в зависимости от числа строк в матрице ограничений.
Для задач 6 и 7, в которых в однопроцессорном режиме не удалось получить точное решение, были проведены дополнительно расчеты на 12 процессорах. В результате суммарное число итераций сократилось по сравнению с расчетами на восьми процессорах, для задачи 6 — на 4.1 млн итераций и для задачи 7 — на 9.6 млн итераций. Отношение te/t12, где te — время исполнения на восьми процессорах и t12 — на двенадцати, для задач 6 и 7 соответственно составило 1.9 и 1.7.
Анализ результатов численного эксперимента показывает, что параллельный алгоритм обеспечивает достаточно равномерную загрузку процессоров. С увеличением размерности, например, при числе строк, равном нескольким десяткам тысяч, можно ожидать возрастания простоев на начальном этапе у процессорных элементов с номерами i & gt- 0. Здесь был бы полезен вариант MPI с динамическим определением числа используемых процессоров [7].
Параллельный вариант программы включен в пакет программ ЦКП. В него кроме программ решения задач входят программы обработки данных, которые осуществляют логический контроль данных, их перевод из входного MPS-формата [2] во внутренний разреженный формат и другие преобразования. Пакет организован аналогично пакету ЦЛП [4].
Таблица 3
Характеристики процесса решения задач параллельным алгоритмом
Номер задачи
1 2 3 4 5 6 7 8 9
it1 48 243 3 650 064 7 246 301 123 501 339 065 2 319 167 15 065 957 2 018 111 1 466 794
it2 47 246 3 770 439 4 379 823 116 647 324 842 2 500 000 14 285 510 2 005 200 1 490 414
i GU 48 071 3 504 662 4 669 789 121 271 333 464 2 348 354 14 172 405 2 040 903 1 486 702
it4 48 566 3 947 817 4 523 608 115 781 315 766 2 361 807 13 030 163 2 039 631 1 423 969
it5 47 621 3 745 081 3 431 145 113 520 348 888 2 320 268 14 535 137 1 993 103 1 460 954
ite 46 654 3 707 780 3 473 066 115 535 309 939 2 326 606 14 784 738 2 019 802 1 449 547
it7 47 678 3 413 254 3 294 001 115 900 307 952 2 412 781 13 322 902 2 010 772 1 467 747
ite 45 102 3 532 041 3 281 568 119 760 302 719 2 317 997 14 299 376 1 963 402 1 529 396
Sit 379 181 29 271 138 34 299 301 941 915 2 582 635 18 906 980 113 496 188 16 090 924 11 775 523
tp 9. 18 5187. 8851 45. 49 293.4 40 439. 38 292 502.7 528. 6
t/tp 7. 27 8. 74 5. 78 8. 37 10. 42 — - 8. 91 10. 88
Примечание: г^ - суммарное число итераций в оценочных задачах, выполненных на г-м процессоре- Бц — суммарное число итераций, выполненное на всех процессорах- Ьр — время решения задач параллельным алгоритмом, Ь/Ьр — отношение времени решения задачи последовательным алгоритмом ко времени решения параллельным алгоритмом.
Заключение
Реализован параллельный алгоритм целочисленного квадратичного программирования,
который позволяет снизить суммарную трудоемкость решения задач по сравнению с последовательным алгоритмом, и за счет этого обеспечивается ускорение выше линейного.
Список литературы
[1] Ковалев М. М. Дискретная оптимизация (целочисленное программирование). Минск: Изд-во Белорус. ун-та, 1977.
[2] Муртаф Б. Современное линейное программирование. Теория и практика. М.: Мир, 1984.
[3] Уилкинсон Дж.Х., РАйнш К. Справочник алгоритмов на языке Алгол. Линейная алгебра. М.: Машиностроение, 1976.
[4] Zabinyako G.I., Kotel'-nikov E.A. Linear optimization programs // NCC Bulletin, Series Num. Anal. Novosibirsk: NCC Publisher, 2002. Issue 11. P. 103−112.
[5] Корнеев В. Д. Параллельное программирование в MPI. Новосибирск: Изд-во СО РАН, 2000.
[6] FTP: //PLATO. LA. ASU. EDU
[7] Воеводин В. В., Воеводин Вл.В. Параллельные вычисления. СПб.: БХВ Санкт-Петербург, 2002.
Поступила в редакцию 1 октября 2003 г.

ПоказатьСвернуть
Заполнить форму текущей работой