CC BY
31
Математическое моделирование. Оптимальное управление Вестник Нижегородского университета им. Н.И. Лобачевского, 2011, № 3 (2), с. 13-19
УДК 519.853.4
ОЦЕНКИ ЭФФЕКТИВНОСТИ ПАРАЛЛЕЛЬНОГО ИНДЕКСНОГО МЕТОДА ГЛОБАЛЬНОЙ ОПТИМИЗАЦИИ
© 2011 г. К.А. Баркалов
Нижегородский госуниверситет им. Н.И. Лобачевского
Поступила в редакцию 27.01.2011
Рассмотрен параллельный алгоритм решения многоэкстремальных задач с невыпуклыми ограничениями, основанный на сведении исходной многомерной задачи к набору связанных одномерных задач. Схема редукции размерности основана на построении множества отображений типа кривой Пеано, при этом на одномерной шкале сохраняется часть информации о близости точек в многомерном пространстве. Приведены результаты экспериментов, позволяющие оценить эффективность алгоритма при использовании схемы построения множественных отображений.
Ключевые слова: многоэкстремальная оптимизация, невыпуклые ограничения, кривые Пеано, параллельные алгоритмы, эффективность.
Введение
Оригинальный подход к минимизации многоэкстремальных функций при невыпуклых ограничениях, описанный в работах [1-4] и получивший название индексного метода глобальной оптимизации, основан на раздельном учете каждого ограничения задачи и не связан с использованием штрафных функций. При этом решение многомерных задач сводится к решению эквивалентных им одномерных. Соответствующая редукция основана на использовании отображений единичного отрезка вещественной оси на гиперкуб. Роль таких отображений играют аппроксимации кривой Пеано, называемые также развертками (см. [4, 5]). В работе [6] предложено использовать семейство кривых Пеано (вращаемые развертки) с целью построения параллельного алгоритма для многопроцессорных вычислительных систем с распределенной память. Целью данной работы являлось определение показателей эффективности и избыточности для параллельного индексного метода глобальной оптимизации с вращаемыми развертками.
1. Постановка задачи
Рассмотрим многомерную задачу глобальной оптимизации
ф(у*) = тт{ф(у): уеБ, $(у) < 0, 1 <у < ш}, (1.1)
где область поиска Б задана как единичный гиперкуб Ж-мерного евклидова пространства
Б = {уеВ*: -2-1 <у < 2-1, 1 < I < Ж}.
Важный в прикладном отношении подкласс задач вида (1.1) характеризуется тем, что все функции, входящие в определение задачи, заданы некоторыми (программно реализуемыми) алгоритмами вычисления значений ф(у), gj(y), 1 < у < ш, в точках области поиска Б. При этом решение задачи (1.1) сводится к построению оценки у* е Б, отвечающей некоторому понятию близости к точке у* (например, чтобы у* - у* || <е или |ф(/) — ф(у* )| <е, где 8 > 0
есть заданная точность) на основе некоторого числа к значений функционалов задачи, вычисленных в точках области Б.
В задачах многоэкстремальной оптимизации возможность достоверной оценки глобального оптимума принципиально основана на наличии априорной информации о функции, позволяющей связать возможные значения минимизируемой функции с известными значениями в точках осуществленных поисковых итераций. Весьма часто такая априорная информация о задаче (1.1) представляется в виде предположения, что целевая функция ф (в дальнейшем обозначаемая также gm+l) и левые части ограничений gj(y), 1 < у < ш, удовлетворяют условию Липшица с соответствующими константами Ьр, 1 < у < ш+1, а именно
\gj(Уl) - gj(y2)1< Ь\\У1 -^
1 < у < ш+1, у1, у2еБ. (1.2)
В общем случае все эти функции могут быть многоэкстремальными.
Использование развертки Пеано у(х), однозначно отображающей единичный отрезок вещественной оси на единичный гиперкуб, позволяет свести многомерную задачу условной минимизации в области Б к одномерной задаче условной минимизации на отрезке [0,1]
ф(у(х*)) = тт{ф(у(х)):хе[0,1],
(1.3)
gj(y(x)) < 0, 1 <у < ш}.
Рассматриваемая схема редукции размерности сопоставляет многомерной задаче с липши-цевой минимизируемой функцией и липшице-выми ограничениями одномерную задачу, в которой соответствующие функции удовлетворяют равномерному условию Гельдера (см. [4]), т.е.
Ш - gj(y(x"))\ < Щ\х '- х''\ш, х', х''е [0,1], 1 <у < ш+1,
где N есть размерность исходной многомерной задачи, а коэффициенты К связаны с константами Липшица Ь исходной задачи соотношениями К) < 4Ь, у[Й. Вопросы численного построения отображений типа кривой Пеано и соответствующая теория подробно рассмотрены в работах [4, 5]. Здесь же отметим, что численно построенная кривая является приближением к теоретической кривой Пеано с точностью не хуже 2- т по каждой координате (параметр ш называется плотностью развертки).
Для целей дальнейшего изложения введем классификацию точек у из области поиска Б с помощью индекса V = у(у(х)). Указанный индекс
V определяется условиями
gj(y(x)) < 0, 1 <у < V - 1, gv(y(x)) > 0, (1.4)
где последнее неравенство несущественно, если
V = ш+1, и удовлетворяет неравенствам
1 < V = v(y(x)) < ш+1.
Данная классификация порождает функцию
ХКх)) = gv(y(x)), v = Чу(х)Х (1.5)
определенную и вычислимую всюду в Б. Ее значение в точке у есть либо значение левой
части ограничения, нарушенного в этой точке (случай, когда V < ш), либо значение минимизируемой функции (случай, когда V = ш+1). Поэтому определение значения /(у(х)) сводится к последовательному вычислению величин gj(y(x)), 1 <у < V = v(x), т.е. последующее значение gj+1(y(x)) вычисляется лишь в том случае, когда gj(y(x)) < 0. Процесс вычислений завершается либо в результате установления неравенства gj(y(x)) > 0, либо в результате достижения значения v(x) = ш+1.
Описанная процедура, называемая испытанием в точке у, автоматически приводит к определению индекса V этой точки. Пара значений
^ = gv(y), V = V(y), (1.6)
порожденная испытанием в точке уеБ, называется результатом испытания.
Заметим, что максимальный индекс М и значения констант Липшица 1 < V < М, являются неизвестными. Однако эту проблему можно преодолеть, используя вместо этих величин их адаптивные оценки, получаемые в процессе решения задачи на основании результатов испытаний.
2. Использование вращаемых разверток отображений
Редукция многомерных задач к одномерным с помощью разверток имеет такие важные свойства, как непрерывность и сохранение равномерной ограниченности разностей функций при ограниченности вариации аргумента. Однако при этом происходит потеря части информации о близости точек в многомерном пространстве, так как точка хе[0,1] имеет лишь левых и правых соседей, а соответствующая ей точка у(х)еР!м имеет соседей по 2Ж направлениям. Как результат, при использовании отображений типа кривой Пеано близким в Ж-мерном пространстве образам у', у" могут соответствовать достаточно далекие прообразы х , х" на отрезке [0,1].
Сохранить часть информации о близости точек позволяет использование множества отображений
^(х) = •••, /(х)} (2.1)
вместо применения единственной кривой Пеано у(х) (см. [4]). Каждая кривая Пеано У(х) из Уь(х) может быть получена в результате поворота развертки вокруг начала координат. При этом найдется отображение У(х), которое точкам
многомерного пространства у , у , которым при исходном отображении соответствовали достаточно далекие прообразы на отрезке [0,1], будет сопоставлять более близкие прообразы х', х".
Максимальное число различных поворотов развертки, отображающей Ж-мерный гиперкуб на одномерный отрезок, составляет 2Ж. Использование всех из них является избыточным, требуется выбрать лишь часть из всех возможных вариантов.
Предлагается осуществлять преобразование развертки в виде поворота на угол ±п/2 в каждой из координатных плоскостей. Число подобных пар поворотов определяется числом координатных плоскостей пространства, которое равно с2 = N (N -1)/2, а общее число преобразований - N(N-1). Учитывая исходное отображение, приходим к заключению, что данный способ позволяет строить до N(N - 1)+1 разверток для отображения ^мерной области на соответствующие одномерные отрезки.
Использование множества отображений Уь(х) = {у (х), ...,у (х)} приводит к формированию соответствующего множества одномерных многоэкстремальных задач
тіп{ф(у;(х)):хє[0,1], &(у(х)) < 0, 1 < ] < т}, 1 < 1 < Ь.
(2.2)
Каждая задача из данного набора может решаться независимо, при этом любое вычисленное значение 7 = £\,( у'), у' = у' (х') функции gv(y) в г-й задаче может интерпретироваться как вычисление значения 7 = gv( у'), у' = У( х") для любой другой ^-й задачи без повторных трудоемких вычислений функции gv(y). Подобное информационное единство позволяет решать исходную задачу (1.3) путем параллельного решения индексным методом Ь задач вида (2.2) на наборе отрезков [0,1]. Каждая одномерная задача решается на отдельном процессоре. Для организации взаимодействия на каждом процессоре создается Ь очередей, в которые процессоры помещают информацию о выполненных итерациях. Используемая схема не содержит какого-либо единого управляющего процессора, что увеличивает надежность выполняемых вычислений.
3. Организация параллельных вычислений
Использование множественных отображений позволяет решать исходную задачу (1.1)
путем параллельного решения индексным методом Ь задач вида (2.2) на наборе отрезков [0,1]. Каждая одномерная задача решается на
отдельном процессоре. Результаты испытания в
к
точке х , полученные конкретным процессором для решаемой им задачи, интерпретируются как результаты испытаний во всех остальных задачах (в соответствующих точках хк1, ..., хкь). При таком подходе испытание в точке хке[0,1], осуществляемое в я-й задаче, состоит в последовательности действий:
1. Определить образ ук = У(хк) при соответствии У(х).
2. Проинформировать остальные процессоры о начале проведения испытания в точке ук (блокирование точки ук).
3. Вычислить величины g1(yk), ..., gv(yk), где значения индекса V < ш определяются условиями gj(yk) < 0, 1 <у < V, gv(yk) > 0, V < ш.
Выявление первого нарушенного ограниче-
к т»
ния прерывает испытание в точке у . В случае,
к
когда точка у допустима, испытание включает вычисление значений всех функционалов задачи, при этом значение индекса принимается равным величине V = ш+1. Тройка
у\хк), V = v(xk), 7к = gv(ys(xk)) (3.1)
является результатом испытания в точке хк.
4. Определить прообразы хиє[0,1], 1 < 1 < Ь, точки ук и интерпретировать испытание, проведенное в точке укеБ, как проведение испытаний в Ь точках
к1
кЬ
(3.2)
с одинаковыми результатами
v(xk1) = ... = v(xkЬ) = v(xk),
1 к1 Ь кЬ к
gv(y (х )) = ... = gv(y (х )) = 7 .
Проинформировать остальные процессо-
к
ры о результатах испытания в точке у .
Каждый процессор имеет свою копию программных средств, реализующих вычисление функций задачи, и решающее правило алгоритма. Для организации взаимодействия на каждом процессоре создается Ь очередей, в которые процессоры помещают информацию о выполненных итерациях в виде троек: точка очередной итерации, индекс и значение из (3.1), причем индекс заблокированной точки полагается равным -1, а значение функции в ней не определено.
4. Описание алгоритма
Алгоритм выбора точек итераций для всех процессоров одинаков. Начальная итерация осуществляется в произвольной точке х1е(0,1) (начальные точки для всех процессоров задаются различными). Выбор точки х^1, д > 1, любого последующего испытания определяется следующими правилами.
Правило 1. Изъять из всех очередей, закрепленных за данным процессором, записанные для него результаты, включающие множество 7?={у?г:1 < г < точек итераций в области Б и вычисленные в них значения индекса и величин из (3.1). Определить множествоХ?={х?г:1 < г < я?} прообразов точек множества Ук при соответствии у(х).
Правило 2. Точки множества итераций
для относительных разностей функций gv, 1 < V < ш+1. Если множество IV, содержит менее двух элементов, или если ц из (4.2) оказывается равным нулю, то принять ц = 1. Из (4.2) следует, что оценки являются неубывающими начиная с момента, когда (4.2) порождает первое положительное значение ц„.
Правило 5. Для всех непустых множеств Д,,
1 < V < ш + 1, вычислить оценки
I -
1 min{gv (xi): * є Iv }
где вектор
{Xl}uXlU_UXq
(4.З)
перенумеровать нижними индексами в порядке увеличения значений координаты
0 = х0 < х1 < ... < хг-< ... < хк< хк+1 = 1, (4.1)
где к = 1 + я1 + . + яд, и сопоставить им значения 21 = gv(xг•), V = v(xг•), 1 < г < к, вычисленные в этих точках. При этом индекс блокированной точки хг- (т.е. точки, в которой начато проведение испытания другим процессором) полагается равным -1, т.е. v(xг•) = -1, значение 2г- является неопределенным. Точки х0, хк+1 введены дополнительно для удобства последующего изложения, индекс данных точек полагается равным -2, т.е. v(x0) = v(xk+1) = -2, а значения 20, 2к+1 являются неопределенными.
Правило 3. Провести классификацию номеров г, 1 < г < к, точек из ряда (4.1) по числу ограничений задачи, выполняющихся в этих точках, путем построения множеств
I-2 = {0, k+1},
I-1 = {i: 1 < i < k, v(x) = -1},
Iv = {i: 1 < i < k, v(X;) = v}, 1 < v < m + 1,
содержащих номера всех точек хг, 1 < і < к, имеющих индексы, равные одному и тому же значению V. Определить максимальное значение индекса
М = max{v = v(xг•), 1 < г < к}.
Правило 4. Вычислить текущие нижние границы
имеющий положительные координаты, называется вектором резервов.
Правило 6. Для каждого интервала (хг-1, хг),
1 < і < к+1, вычислить характеристику Я(і), где
Л(І) = д, + (2- ^)2 - 2(- 2^ >,
Г М V Д і Г М V
V = V(X- ) ^ (X),
(2- — 2*)
т = 2А; — 4^^-^, уСх-О <у(х,-) = у,
Г ^
R(i) = 2A, - 4 (z‘-1 Zv), V = v(x, j) > v(x,),
rv Mv Ai = (Xi - Xi-i)1/N.
Величины rv> 1, 1 < v < m + 1, являются параметрами алгоритма.
Правило 7. Определить интервал (xt-1, xt), которому соответствует максимальная характеристика
R(t) = max{R(i): 1 < i < k+1}. (4.4)
Правило S. Провести очередное испытание в серединной точке интервала (xt -1, xt), если индексы его концевых точек не совпадают, т.е.
xq+1 = (Xt + Xt-i)2, v(xt -i) Ф v(xt).
В противном случае провести испытание в точке
xq+1 = (xt + Xt-i) / 2 - sign(zt-zt-i) —
2rv
v = v(Xt -i) = v(Xt).
N
*
zv =
zt- zt-1
и
v
Результаты испытания занести во все очереди, закрепленные за данным процессором. Увеличив д на единицу, перейти к новой итерации.
Описанные правила можно дополнить условием остановки, прекращающим испытания, если At < е, где ? из (4.4) и е > 0 имеет порядок желаемой покоординатной точности в задаче.
5. Операционные характеристики и сравнение алгоритмов
Один из известных подходов к оценке эффективности методов безусловной глобальной оптимизации основан на численном решении этими методами всех задач из некоторой случайно генерируемой выборки большого объема. Генераторы таких выборок можно рассматривать как некоторые классы функций Ау), уеБ, с определенной на них вероятностной мерой.
Среди подобных генераторов следует отметить используемый в данной работе ОКЬ8-генератор, описанный в [7] и доступный для свободного скачивания. Отличительной особенностью ОКЬ8-генератора является прозрачный способ задания сложности генерируемых задач: параметрами генератора являются количество локальных минимумов задачи, размеры их областей притяжения, удаленность глобального минимума от локальных и многое другое. Всё это определяет востребованность ОКЬ8-генератора, который используется более чем в 40 странах для тестирования методов глобальной оптимизации. Для сравнения методов использовались описанные в [8] несколько классов тестовых задач по 100 функций каждый. Область поиска Б во всех задачах составляет гиперкуб [-1.0, 1.0]ж.
Примененная процедура оценки основана на операционных характеристиках из [9] и состоит в следующем. Пусть некоторая задача из рассматриваемой выборки решается с помощью алгоритма 5. При этом задаче сопоставляется порождаемая алгоритмом последовательность точек испытаний { ук} . Указанная последова-
Надежность Р(к)
Итерации к
тельность усекается (т.е. процесс решения прекращается) либо в связи с первым попаданием точки очередного испытания в заданную 8-окрестность решения у , либо в связи с тем, что в ходе выполнения заданного числа Ктах испытаний такое попадание не имело места. В проведенных численных экспериментах использовалось значение Ктах = 10000.
Результат решения всех задач выборки с помощью алгоритма 5 представляется функцией Р (к), характеризующей долю задач, для которых в ходе к шагов поиска имели место попадания точек испытаний в заданную 8-окрестность решения. Такую функцию будем называть операционной характеристикой алгоритма 5. Отметим, что при использовании ОКЬ8-генератора 8-окрестность решения задачи заранее известна.
Эксперименты проводились для двух реализаций алгоритма с вращаемыми развертками: последовательной и параллельной при плотности ш = 10, числе разверток Ь = 6. При этом использовались параметры надежности г = 3.8 (данные параметры являются минимальными для соответствующих методов, при которых достигается решение 100% задач из выборки); точность попадания в окрестность решения составляла 8 = 0.01.
Операционные характеристики для обоих методов, полученные на некоторых классах, порожденных ОКЬ8-генератором, представлены соответственно на рис. 1 и 2. При этом нижняя разрывная кривая характеризует последовательный метод, а верхняя непрерывная - параллельный метод с вращаемой разверткой. Указанное расположение кривых показывает, что параллельный алгоритм с вращаемыми развертками обеспечивает в среднем значительно более быстрое получение оценок, лежащих в заданной окрестности решения, чем его последовательный прототип.
В проведенных экспериментах доля правильно решенных задач определяет надеж-
Надежность Р(к)
'
с
/
/ *
/
/ /
/ у
^ /
(/
!/
0 1000 2000 3000 4000 3300 600
Итераци
Рис. 1
Рис. 2
ность поиска решения, а число выполненных испытаний целевой функции - затратность при поиске решения. Для последовательного метода число итераций k совпадает с количеством испытаний K, а при применении параллельной версии K я k-L.
Отметим, что оценка по итерациям косвенно характеризует выигрыш по времени от применения параллельных реализаций. Таким образом, выигрыш по времени на первом классе функций (при условии правильного решения всех задач) примерно равен 1.9 раза, а на втором он составил 2.4 раза, т.е. коэффициент ускорения возрастает с ростом сложности задач. Аналогичный эффект наблюдается при увеличении размерности (отчасти - за счет возможности увеличения коэффициента распараллеливания, совпадающего с числом используемых разверток L).
Дадим оценки избыточности. Анализируя результаты, приведенные на рис. 1, можно оценить общее количество вычислений функций, выполненных параллельной версией, исходя из общего числа итераций (оно составило 475), необходимого для решения 100% задач. Общее число измерений функции (в худшем случае) при этом оценивается как 475L = 28500 на задачу. При этом для полного решения всех задач в худшем случае последовательной версии потребовалось 900 итераци. Таким образом, коэффициент избыточности по количеству испытаний при L = 6 за счет параллелизма может быть оценен как 28500/900 = 3.17. Аналогичный способ оценивания для второго эксперимента (рис. 2) дает оценку избыточности, равную 2500L/6000 = 2.5.
Укажем также абсолютное время, затраченное программной реализацией алгоритма на поддержку собственного функционирования в условиях суперкомпьютерных вычислений (без учета времени вычисления функций). Оно составило в проведенных экспериментах (в среднем на одну задачу) 0.6 с на каждый процессор. Общее время вычислений зависит от затрат на одно вычисление функции, т.е. определяется их вычислительной сложностью. Указанные временные характеристики соответствуют следующим показателям использованного оборудования. При проведении экспериментов использовался высокопроизводительный кластер, приобретенный Нижегородским университетом в 2007 г., с двухпроцессорными двухъядерными серверами Intel Xeon 3.2 GHz, число разверток
соответствовало числу вычислительных ядер. Размер оперативной памяти, выделяемой каждому ядру, составлял 1 Гб.
Заключение
Можно отметить, что схема построения множества кривых Пеано, основанная на вращении разверток, ориентирована на параллельные вычислительные системы с большим числом узлов. Указанная схема реализована в параллельной системе глобальной оптимизации Global Expert. С помощью Global Expert проведены вычислительные эксперименты, результаты которых позволяют оценить эффективность параллельного индексного метода.
Работа выполнена при поддержке Совета по грантам Президента Российской Федерации (грант № МК-1536.2009.9 и № НШ-64729.2010.9).
Список литературы
1. Стронгин Р.Г., Маркин Д.Л. Минимизация многоэкстремальных функций при невыпуклых ограничениях // Кибернетика. 1986. № 4. С. 63-69.
2. Стронгин Р.Г. Поиск глобального оптимума. М.: Знание, 1990.
3. Стронгин Р.Г., Баркалов К.А. О сходимости индексного алгоритма в задачах условной оптимизации с s-резервированными решениями // Математические вопросы кибернетики. М.: Наука, 1999. С. 273-288.
4. Strongin R.G., Sergeyev Ya.D. Global optimization with non-convex constraints. Sequential and parallel algorithms. Kluwer Academic Publishers, Dordrecht, 2000.
5. Стронгин Р.Г. Численные методы в многоэкстремальных задачах. М.: Наука, 1978.
6. Стронгин Р.Г., Гергель В.П., Баркалов К. А. Параллельные методы решения задач глобальной оптимизации // Известия вузов. Приборостроение. 2009. Т. 52. № 10. С. 25-32.
7. Gaviano M., Kvasov D.E., Lera D., Sergeyev Ya.D. Software for generation of classes of test functions with known local and global minima for global optimization: [http://si.deis.unical.it/~yaro/GKLS.html] // ACM TOMS. 2003. Vol. 29. № 4. P. 469-480.
8. Lera D., Sergeyev Ya.D. Lipschitz and Holder global optimization using space-filling curves // Applied Numerical Mathematics. January 2010. Vol. 60. № 1-2. P. 115-129.
9. Гришагин В.А. Операционные характеристики некоторых алгоритмов глобального поиска // Проблемы случайного поиска. Рига: Зинатне, 1978. № 7. С. 198-206.
EFFICIENCY ESTIMATIONS OF PARALLEL INDEX METHOD OF GLOBAL OPTIMIZATION
K.A. Barkalov
A parallel algorithm for solving multiextremal problems with nonconvex constraints is considered. It is based on the reduction of a multidimensional problem to a set of related one-dimensional ones. The dimension reduction scheme is based on Peano-type space-filling curves mapping, the information about proximity of the points in the multidimensional space being preserved in a one-dimensional address space. The algorithm efficiency in multidimensional mapping constructions has been illustrated by the results of numerical experiments.
Keywords: multiextremal optimization, nonconvex constraints, Peano space-filling curves, parallel algorithms, efficiency.
