Параллельные вычисления в задачах многоэкстремальной оптимизации Текст научной статьи по специальности «Математика»

Математическое моделирование. Оптимальное управление Вестник Нижегородского университета им. Н.И. Лобачевского, 2009, № 6 (1 ), с. 1 71-1 77

УДК 519.853.4

ПАРАЛЛЕЛЬНЫЕ ВЫЧИСЛЕНИЯ В ЗАДАЧАХ МНОГОЭКСТРЕМАЛЬНОЙ ОПТИМИЗАЦИИ

© 2009 г. КА. Баркалов, С.В. Сидоров, В.В. Рябов

Нижегородский госуниверситет им. Н.И. Лобачевского [email protected]

Поступилн врвднкцию 15.05.2009

Рассмотрен параллельный алгоритм решения многоэкстремальных задач с невыпуклыми ограничениями, основанный на сведении исходной многомерной задачи к набору связанных одномерных задач. Предложена новая схема построения множества отображений типа кривой Пеано, сохраняющая часть информации о близости точек в многомерном пространстве. Приведены результаты экспериментов, подтверждающие эффект ускорения сходимости алгоритма при использовании новой схемы построения множественных отображений.

Клюиввыв словн: многоэкстремальная оптимизация, невыпуклые ограничения, кривые Пеано, параллельные алгоритмы.

Введение

Данная работа продолжает развитие известного подхода к минимизации многоэкстремальных функций при невыпуклых ограничениях, описанного в [1-6] и получившего название индексного метода глобальной оптимизации. Подход основан на раздельном учете каждого ограничения задачи и не связан с использованием штрафных функций. При этом решение многомерных задач сводится к решению эквивалентных им одномерных. Соответствующая редукция основана на использовании разверток единичного отрезка вещественной оси на гиперкуб. Роль таких разверток играют непрерывные однозначные отображения типа кривой Пеано, называемые также кривыми, заполняющими пространство. Предложена новая схема построения множества кривых Пеано («вращаемые развертки»), достоинством которой (в отличие от схемы «сдвиговых разверток») является простота построения и использования. Описан параллельный индексный алгоритм, основанный на одновременном использовании множества вращаемых разверток. Приведены результаты экспериментов, убедительно подтверждающие достоинство новой схемы построения множественных отображений.

1. Постановка задачи

Рассмотрим многомерную задачу глобальной оптимизации

ф(у*) = тт{ф(у): yeD, gj(y) < 0, 1 <] <т}, (1.1)

где область поиска D задана как единичный гиперкуб ^-мерного евклидова пространства

D={yеR: -2-1<уг<2-1, 1</<^}.

Важный в прикладном отношении подкласс задач вида (1.1) характеризуется тем, что все функции, входящие в определение задачи, заданы некоторыми (программно реализуемыми) алгоритмами вычисления значений ф(у), gj(y), 1 <] <т, в точках области поиска D. При этом решение задачи (1.1) сводится к построению оценки у* е D , отвечающей некоторому понятию близости к точке у (например, чтобы ||у* -у*|| <8 или |ф(у*) — ф(у*)| <8 , где

е > 0 есть заданная точность), на основе некоторого числа k значений функционалов задачи, вычисленных в точках области D.

В задачах многоэкстремальной оптимизации возможность достоверной оценки глобального оптимума принципиально основана на наличии априорной информации о функции, позволяющей связать возможные значения минимизируемой функции с известными значениями в точках осуществленных поисковых итераций. Весьма часто такая априорная информация о задаче (1.1) представляется в виде предположения, что целевая функция ф (в дальнейшем обозначаемая также gm+1) и левые части ограничений gj(y),1<j<m, удовлетворяют условию Липшица с соответствующими константами Lj, 1<<т+1, а именно

I gj(Уl)—gj(У2) I <^-|У1—У2||, 1</'<т+1, уь y2еD. (1.2)

В общем случае все эти функции могут быть

многоэкстремалъными.

Использование развертки Пеано у(х), однозначно отображающей единичный отрезок вещественной оси на единичный гиперкуб, позволяет свести многомерную задачу условной минимизации в области D к одномерной задаче условной минимизации на отрезке [0,1]

ф(у(х*))=тт{ф{у(х)): хе[0,1], gj(y(x'))<0, 1</<т}.(1.3)

Рассматриваемая схема редукции размерности сопоставляет многомерной задаче с липши-цевой минимизируемой функцией и липшице-выми ограничениями одномерную задачу, в которой соответствующие функции удовлетворяют равномерному условию Гельдера (см. [1]), т.е.

^(х'

X, х"е[0,1], 1</<т+1, где N есть размерность исходной многомерной задачи, а коэффициенты К связаны с константами Липшица Lj исходной задачи соотношениями Kj<4Ljл[N . Вопросы численного построения отображений типа кривой Пеано и соответствующая теория подробно рассмотрены в работах [1, 6]. Здесь же отметим, что численно построенная кривая является приближением к теоретической кривой Пеано с точностью не хуже 2 т по каждой координате (параметр т называется плотностью развертки).

Для целей дальнейшего изложения введем классификацию точек у из области поиска D с помощью индекса v=v(y(x)). Указанный индекс

V определяется условиями

gj(y(x))<0, 1<7<л—1, g1Vy(x))>0, (1.4)

где последнее неравенство несущественно, если v=m +1, и удовлетворяет неравенствам

Данная классификация порождает функцию

Яу&^М^Х ^(уС^Х (1.5)

определенную и вычислимую всюду в D. Ее значение в точке у есть либо значение левой части ограничения, нарушенного в этой точке (случай, когда v<m), либо значение минимизируемой функции (случай, когда v=m + 1). Поэтому определение значения Ау^)) сводится к последовательному вычислению величин gj(y(x)), 1<j'<v=v(x), т.е. последующее значение gj+1(y(x)) вычисляется лишь в том случае, когда gj(y(x))<0. Процесс вычислений завершается либо в результате установления неравенства gj(y(x))>0, либо в результате достижения значения v(x)=m + 1.

Описанная процедура, называемая испытанием в точке у, автоматически приводит к определению индекса V этой точки. Пара значений

^Лу), v=v(y), (1.6)

порожденная испытанием в точке yеD, называется результатом испытания.

Заметим, что максимальный индекс М и значения констант Липшица LV, 1<л<М, являются неизвестными. Однако эту проблему можно преодолеть, используя вместо этих величин их адаптивные оценки, получаемые в процессе решения задачи на основании результатов испытаний.

2. Использование множественных отображений

Редукция многомерных задач к одномерным с помощью разверток имеет такие важные свойства, как непрерывность и сохранение равномерной ограниченности разностей функций при ограниченности вариации аргумента. Однако при этом происходит потеря части информации о близости точек в многомерном пространстве, так как точка xе[0,1] имеет лишь левых и правых соседей, а соответствующая ей точка y(x)еRN имеет соседей по 2 направлениям. Как результат, при использовании отображений типа кривой Пеано близким в #-мерном пространстве образам у, у" могут соответствовать достаточно далекие прообразы x', х" на отрезке [0,1].

Сохранить часть информации о близости точек позволяет использование множества отображений

7L(x)={y1(x),..., у^)} (2.1)

вместо применения единственной кривой Пеано у^) (см. [4, 6]). Каждая кривая Пеано у^) из YL(x) может быть получена в результате некоторого сдвига вдоль главной диагонали гиперинтервала D. Таким образом сконструированное множество кривых Пеано позволяет получить для любых близких образов у, у", отличающихся только по одной координате, близкие прообразы x', x" для некоторого отображения у^).

3. Вращаемые развертки

К числу недостатков ставшей уже классической схемы построения множественных разверток (далее будем называть их сдвиговыми развертками или С-развертками) можно отнести наличие дополнительного ограничения, порождающего сложную допустимую область на одномерных отрезках.

Преодолеть этот недостаток, сохранив информацию о близости точек в #-мерном пространстве, позволяет новая схема построения множественных отображений. Отличительной чертой предложенной схемы является построе-

"1 0 0 0 0" "1 0 0 0 0"

0 0 0 -1 0 0 0 0 1 0

0 0 1 0 0 0 0 1 0 0

0 1 0 0 0 0 -1 0 0 0

0 0 0 0 1 0 0 0 0 1

Рис. 2

ние множества кривых Пеано не с помощью сдвига вдоль главной диагонали гиперкуба, а поворотом развертки вокруг начала координат. При этом найдется отображение У^), которое точкам многомерного пространства у, у", которым при исходном отображении соответствовали достаточно далекие прообразы на отрезке [0,1], будет сопоставлять более близкие прообразы x', x".

Развертки, порождаемые в соответствии с новой схемой, будем называть вращаемыми развертками или В-развертками.

Для иллюстрации на рис. 1 изображены две кривые, являющиеся приближением к развертке Пеано для случая N = 2. Узлы сетки в #-мерном пространстве отмечены темными кружками, стрелкой указана одна из пар точек, прообразы которых далеки на одномерной оси при одном отображении и близки при другом.

Максимальное число различных поворотов развертки, отображающей #-мерный гиперкуб на одномерный отрезок, составляет 2 . Использование всех из них является избыточным, требуется выбрать лишь часть из всех возможных вариантов.

Предлагается осуществлять преобразование развертки в виде поворота на угол ±п/2 в каждой из координатных плоскостей. В качестве

примера на рис. 2 приведены матрицы поворота в плоскости (у2,у4) при N = 5.

Число подобных пар поворотов определяется числом координатных плоскостей простран-

^2 N(N -1) г

ства, которое равно LN =------^----, а общее

число преобразований будет равно N(N-1). Учитывая исходное отображение, приходим к заключению, что данный способ позволяет строить до N(N-1)+! развертки для отображения ^мерной области на соответствующие одномерные отрезки. При этом дополнительное ограничение, которое возникало при построении сдвиговых разверток [4], отсутствует.

В случае необходимости данный способ построения множества отображений может быть легко «отмасштабирован» для получения большего (вплоть до 2^ числа разверток.

Использование множества отображений

YL(x) = (у(х),...,у (х)} приводит к формированию соответствующего множества одномерных многоэкстремальных задач

тт{ф(у(х)): хє[0,1], gJ(yl(x))<0, 1</<т}, 1<1<Ь. (3.1) Каждая задача из данного набора может решаться независимо, при этом любое вычисленное значение z=gv(y'), у'=У(х|) функции gl^(y) в г-й задаче может интерпретироваться как вычисление значения z=gxiyt), у'=у(х") для любой другой 5-й

задачи без повторных трудоемких вычислений функции gVy). Подобное информационное единство позволяет решать исходную задачу (1.3) путем параллельного решения индексным методом L задач вида (3.1) на наборе отрезков [0,1]. Каждая одномерная задача решается на отдельном процессоре. Для организации взаимодействия на каждом процессоре создается L очередей, в которые процессоры помещают информацию о выполненных итерациях. Используемая схема не содержит какого-либо единого управляющего процессора, что увеличивает надежность выполняемых вычислений.

Проинформировать остальные процессоры о

k

результатах испытания в точке у .

Каждый процессор имеет свою копию программных средств, реализующих вычисление функций задачи, и решающее правило алгоритма. Для организации взаимодействия на каждом процессоре создается L очередей, в которые процессоры помещают информацию о выполненных итерациях в виде троек: точка очередной итерации, индекс и значение из (4.1), причем индекс заблокированной точки полагается равным -1, а значение функции в ней не определено.

4. Организация параллельных вычислений

Использование множественных отображений позволяет решать исходную задачу (1.1) путем параллельного решения индексным методом L задач вида (3.1) на наборе отрезков [0,1]. Каждая одномерная задача решается на

отдельном процессоре. Результаты испытания в

k

точке x, полученные конкретным процессором

для решаемой им задачи, интерпретируются как

результаты испытаний во всех остальных зада/ __ М kL\ Т—г

чах (в соответствующих точках x ,..., x ). При таком подходе испытание в точке x е [0,1], осуществляемое в 5-й задаче, состоит в последовательности действий:

1. Определить образ yk=ys(xk) при соответствии у5 ^).

2. Проинформировать остальные процессо-

k

ры о начале проведения испытания в точке у

(блокирование точки ук).

3. Вычислить величины g1(yk),..., gv(yk), где значения индекса v<m определяются условиями

gj(yk)<0, 1<7'<Л, g^Сyk)>0, Л<т. Выявление первого нарушенного ограничения прерывает испытание в точке yk. В случае k

когда точка у допустима, испытание включает вычисление значений всех функционалов задачи, при этом значение индекса принимается равным величине л=т+1. Тройка

(4.1)

5/ к\ / к\ к / 5/ К\

у (х ), у=у(х ), Z =gv(y (х )),

k

является результатом испытания в точке x .

4. Определить прообразы xие[0,1], 1</<L,

k

точки у и интерпретировать испытание, проведенное в точке у еD, как проведение испытаний в L точках

'1 " (4.2)

к1 к^

х ,..., х

с одинаковыми результатами

/ к1\ / кL\ / к\

у(х )=...=у(х )=у(х ),

/ 1/ к1\\ / L/ к^\ к

gv(y(x ))=.=gv(y(x хи?.

5. Описание алгоритма

Алгоритм выбора точек итераций для всех процессоров одинаков. Начальная итерация осуществляется в произвольной точке x1е(0,1) (начальные точки для всех процессоров задаются различными). Выбор точки xq+1, д>1, любого последующего испытания определяется следующими правилами.

Правило 1. Изъять из всех очередей, закрепленных за данным процессором, записанные для него результаты, включающие множество Уч={удг: 1</<59} точек итераций в области D и вычисленные в них значения индекса и величин из (4.1). Определить множество X(г={xqг: 1</<59} прообразов точек множества Yk при соответствии у^).

Правило 2. Точки множества итераций {x1}u^1u.uXq перенумеровать нижними индексами в порядке увеличения значений координаты

0=x0<x1<...<xг<...<xk<xk+1=1, (5.1)

где Л=1+51+.+59, и сопоставить им значения zi=gv(xi), v=v(xг), 1</<k, вычисленные в этих точках. При этом индекс блокированной точки Xi (т.е. точки, в которой начато проведение испытания другим процессором) полагается равным -1, т.е. v(xг■)=-1, значение Zi является неопределенным. Точки x0, xk+1 введены дополнительно для удобства последующего изложения, индекс данных точек полагается равным -2, т.е. v(x0)=v(xk+1)=-2, а значения z0, zk+1 являются неопределенными.

Правило 3. Провести классификацию номеров /, 1<1<к, точек из ряда (5.1) по числу ограничений задачи, выполняющихся в этих точках, путем построения множеств

1-2={0, Ж},

/_1={г: 1</<k, v(xг■)=-1},

^={1: 1<л<к, v(xг■)=v}, 1^^+^

содержащих номера всех точек хг-,1</<к, имеющих индексы, равные одному и тому же значению V. Определить максимальное значение индекса

М_тах{у=у(хі), 1<і<к}.

Правило 4. Вычислить текущие нижние гра-

ницы

(х , - х,)

(5.2)

для относительных разностей функций gv,1<v<m+1. Если множество Іу содержит менее двух элементов или если из (5.2) оказывается равным нулю, то принять ц,=1. Из

(5.2) следует, что оценки ц, являются неубывающими начиная с момента, когда (5.2) порождает первое положительное значение ц.

Правило 5. Для всех непустых множеств І,, 1<у<т+1, вычислить оценки

2 *_{-^, у< м,

2 [min{gv (хг): і є Іу}, у_М,

где вектор

єя_(єь..., Єт), (5.3)

имеющий положительные координаты, называется вектором резервов.

Правило 6. Для каждого интервала (х-х), 1<і<к+1, вычислить характеристику Я(і), где

Я(і) _&, + - 2- 2 + 2-■ - 2*)

у_Ч х-1) _Ч х);

Я(і) _ 2Д. - 4 (2 - ^), у(х. 1) <у(х.) _у;

ГуЦу

Я (і) _ 2Д,. - 4 (2-1 - ^), V _ у(хі-1) > у(х,.);

Г ц

Ді _ (хі - хі-і) .

Величины Гv>1, 1<у<т+1, являются параметрами алгоритма.

Правило 7. Определить интервал (хьЬ хг), которому соответствует максимальная характеристика

К(ґ) _ тах{Я(і): 1<і<к+1}. (5.4)

Правило 7. Провести очередное испытание в серединной точке интервала (хьЬ хґ), если индексы его концевых точек не совпадают, т.е. хч+ _ (х,+хм) / 2, у(хм) ф у(х)

В противном случае провести испытание в точке

хч+х _ х+хы) / 2^п^-2М)-1-

2гv

V _ у(хм) _ у(х)

Результаты испытания занести во все очереди, закрепленные за данным процессором. Увеличив q на единицу, перейти к новой итерации.

Описанные правила можно дополнить условием остановки, прекращающим испытания, если Дґ<є, где ґ из (5.4) и є>0 имеет порядок желаемой покоординатной точности в задаче.

6. Операционные характеристики и сравнение алгоритмов

Один из известных подходов к оценке эффективности методов безусловной глобальной оптимизации основан на численном решении этими методами всех задач из некоторой случайно генерируемой выборки большого объема (см., например, [7]).

Генераторы таких выборок можно рассматривать как некоторые классы функций Д(у), yєD, с определенной на них вероятностной мерой. Один из таких генераторов G, предложенный в работе [7], порождает многоэкстремальные функции в соответствии со схемой

Ф(у^ у2) _-'

ХХ[4а ^ у 2) + ВаЬа у2)]

1 м

ХХ^а( у 2) -,( у2)]

і_1 ,_1

где а ц(у^у2) = sin(ПУl)sin(пjУ2), ‘Муь у2) =

= cos(л/y1)cos(п/y2) определены в области 0 < у1,у2 < 1, а параметры -1<Л„, Бу, Су, Dц<1 являются независимыми случайными величинами, равномерно распределенными в указанном выше интервале. Минимизация подобных функций возникает, например, в задаче оценки максимального напряжения (определяющего прочность) в тонкой пластине при поперечной нагрузке.

Примененная процедура оценки эффективности алгоритмов использует аппарат операционных характеристик из [8] и состоит в следующем. Пусть некоторая задача вида (1.1) из рассматриваемой выборки решается с помощью алгоритма S. При этом задаче сопоставляется порождаемая алгоритмом последовательность точек испытаний {у }. Указанная последовательность усекается (т.е. процесс решения прекращается) либо в связи с первым попаданием точки очередного испытания в заданную 8-окрестность решения x*, либо в связи с тем, что

2

+

2

Р(к)

Рис. 3

в ходе выполнения заданного числа К испытаний такое попадание не имело места. В проведенных численных экспериментах использовалось значение К_1000.

Поскольку задание 8-окрестности предпола-* ^

гает координату х известной, ее значение предварительно оценивалось (для каждой функции выборки) путем перебора всех узлов равномерной (в области поиска D) сетки с шагом 10 3.

Результат решения всех задач выборки с помощью алгоритма £ представляется функцией Р£к), характеризующей долю задач, для которых в ходе к шагов поиска имели место попадания точек испытаний в заданную 8-окрестность решения. Такую функцию будем называть операционной характеристикой алгоритма &

Эксперименты проводились для описанных выше алгоритмов со сдвиговыми (С) и вращаемыми (В) развертками с плотностью т=12, число разверток L=2. При этом использовались параметры гу=2.1, є_0.01 и резервы єу_0.05, 0<у<1, из (5.3). Данные параметры являются минимальными для соответствующих методов, при которых достигается решение 100% задач из выборки.

Операционные характеристики обоих методов, полученные на выборке, порожденной генератором G, представлены на рис. 3 (к - число итераций, Р(к) - доля решенных задач). При

этом нижняя разрывная кривая характеризует метод со сдвиговыми развертками, средняя разрывная кривая - метод с единственной разверткой, а верхняя непрерывная линия - метод с вращаемыми развертками. Указанное расположение кривых показывает, что алгоритм с В-развертками обеспечивает в среднем более быстрое получение оценок, лежащих в заданной окрестности решения, чем алгоритм с С-развертками или же алгоритм с единственной разверткой.

С целью иллюстрации использования параллельного алгоритма с множеством В-разверток рассмотрим известную задачу минимизации функции Растригина

N . .

ф(у)=Х(у/2 -^(18у/2)), -1.5<у/ < 1.5, 1</<Ы,

/=1

N=6,

где минимальное значение ф(у ) = -N достигается в точке у =0.

Для решения данной задачи применялись последовательный и параллельный методы с вращаемой разверткой с плотностью т=10, использовался параметр метода г=2.0 и е=0.05 в критерии остановки. Число разверток L=30 соответствовало числу процессоров, задействованных при решении задачи. И последовательный, и параллельный алгоритмы нашли решение с требуемой точностью, при этом последовательный алгоритм

затратил 173116 итераций, а параллельный - 8535 (максимальное число итераций на одном процессоре). Ускорение по числу итераций составило 20.28, а ускорение по времени - 7.48.

Работа выполнена при финансовой поддержке РФФИ (грант № 07-01-00467-а) и Совета по грантам Президента Российской Федерации по государственной поддержке ведущих научных школ Российской Федерации (гранты № НШ-4694.2008.9 и № МК-1536.2009.9).

Список литературы

1. Стронгин Р.Г. Численные методы в многоэкстремальных задачах. М.: Наука, 1978.

2. Стронгин Р.Г., Маркин Д.Л. Минимизация многоэкстремальных функций при невыпуклых ограничениях // Кибернетика. 1986. № 4. С. 63-69.

3. Стронгин Р.Г. Поиск глобального оптимума. М.: Знание, 1990.

4. Стронгин Р.Г. Параллельная многоэкстремальная оптимизация с использованием множества разверток // Журн. вычисл. матем. и матем. физ. 1991. Т. 31. № 8. С. 1173-1185.

5. Стронгин Р.Г., Баркалов К.А. О сходимости индексного алгоритма в задачах условной оптимизации с е-резервированными решениями// Математические вопросы кибернетики. М.: Наука, 1999. С. 273-288.

6. Strongin R.G., Sergeyev Ya.D. Global optimization with non-convex constraints. Sequential and parallel algorithms. Kluwer Academic Publishers, Dordrecht, 2000.

7. Гергель В.П., Стронгин Р.Г. Абсолют. Программная система для исследований и изучения методов глобальной оптимизации. Н. Новгород: Изд. Нижегород. ун-та, 1998.

8. Гришагин В.А. Операционные характеристи-

PARALLEL COMPUTING IN MULTIEXTREMAL OPTIMIZATION PROBLEMS

K.A. Barkalov, S. V. Sidorov, V. V. Ryabov

A parallel algorithm for solving multiextremal optimization problems with nonconvex constraints has been considered. It is based on the reduction of an initial multidimensional problem to a set of related one-dimensional ones. A new scheme to construct a set of Peano space-filling curves has been proposed. The scheme preserves some information on the closeness of points in the multidimensional space. The results of numerical experiments have been presented to show the algorithm convergence acceleration using the new scheme of multiple curve construction.

Keywords: multiextremal optimization, nonconvex constraints, Peano curves, parallel algorithms.

ки некоторых алгоритмов глобального поиска // Проблемы случайного поиска. Рига: Зинатне, 1978. № 7. С. 198-206.