CC BY
48
Математическое моделирование. Оптимальное управление Вестник Нижегородского университета им. Н.И. Лобачевского, 2012, № 1 (1), с. 144-155
УДК 519.853.4
ТРИАНГУЛЯЦИОННЫЕ МЕТОДЫ ПАРАБОЛОИДОВ В ЗАДАЧАХ МНОГОЭКСТРЕМАЛЬНОЙ ОПТИМИЗАЦИИ С ОГРАНИЧЕНИЯМИ ДЛЯ КЛАССА ФУНКЦИЙ С ЛИПШИЦЕВЫМИ ПРОИЗВОДНЫМИ ПО НАПРАВЛЕНИЯМ
© 2012 г. С.Ю. Городецкий
Нижегородский госуниверситет им. Н.И. Лобачевского
gorosyu@gmail. com
Поступила в редакцию 17.07.2011
Предложен новый метод решения многоэкстремальных задач с невыпуклыми ограничениями для класса функций с липшицевыми производными по направлениям. Измеряются только значения функций. Применена адаптивная триангуляция области поиска нерегулярными симплексами. Миноранты функций в симплексах строятся в виде параболоидов. Метод редуцирует задачу с ограничениями к задаче на гиперинтервале с перестраиваемой целевой функцией. Получены условия сходимости процесса редуцирования. Приведены результаты вычислительных экспериментов.
Ключевые слова: многоэкстремальная оптимизация, невыпуклые ограничения, компонентный подход, триангуляция, симплексы, липшицевы производные по направлениям.
Введение
Компонентный подход является одним из способов построения вычислительно реализуемых алгоритмов многоэкстремальной оптимизации, использующих оценки поведения функций по конечному числу проведенных испытаний. В методах такого типа начальная область поиска динамически разбивается на подобласти-компоненты определенной структуры (как правило — гиперинтервалы) таким образом, чтобы в каждой компоненте было размещено определенное количество точек испытаний функций задачи (как правило — одна или две). Структура компонент и расположение точек измерений в них выбираются так, чтобы существовали простые (желательно, конечные) методы вычисления по этим измерениям оценок поведения функций задачи в пределах компоненты. На основе полученных оценок выбирается компонента, подвергаемая делению с проведением дополнительных измерений. Методам такого типа посвящен ряд работ разных авторов [1-7].
В данной статье применяется компонентный подход нетрадиционной структуры, предложенный в [3, 4], основанный на общей идее адаптивной триангуляции области поиска D, задаваемой в виде N-мерного гиперинтервала. Множество D последовательно разбивается на не пересекающиеся по внутренним областям нерегулярные полноразмерные симплексы S, которые в совокупности покрывают D. Первый из методов такого типа, названных SM-методами (от Simplex Metods), был предложен и
исследован в работе [3], опубликованной в 1999 году, применительно к задаче без функциональных ограничений с липшицевой целевой функцией. В [4] впервые было приведено описание аналогичного метода для класса дифференцируемых функций с липшицевыми производными по направлениям для случая испытаний, включающих измерения градиента. Методы многоэкстремальной оптимизации, использующие симплексы, рассматривались позднее в работах других авторов для липшицевых задач без функциональных ограничений. Ряд таких методов предложен и исследован в работах А.Г. Жилинскаса (см. [8, 9] и библиографию в [10]). Метод, основанный на триангуляции Делоне, предложен в [11]. В [8] вместо триангуляционного разбиения использовалось покрытие области D правильными симплексами, имеющими существенные пересечения по внутренним областям, однако в [9,10] применяется триангуляционное разбиение, аналогичное использованному в [3, 4] .
В рассматриваемых в данной работе SM-методах вершины симплексов размещаются в точках испытаний (т.е. измерений функций, входящих в постановку решаемой задачи). Начальная триангуляция строится по точкам измерений, проводимых в вершинах исходного гиперинтервала D, а также в его геометрическом центре и центрах тяжести всех граней D размерности от 1 до N - 1 (всего 3^ точек). Новое испытание проводится в симплексе £, являющемся «лучшим» по отношению к введенным характеристикам-приоритетам симплексов, и
размещается в центре его наибольшего ребра. Выбранный «лучший» симплекс, а также другие симплексы, включающие делимое ребро симплекса S, разделяются (каждый из них) на пару новых симплексов меньшей меры с использованием точки нового испытания. Таким образом, в отличие от [9-11], каждое испытание может приводить к делению не одного, а группы симплексов (такой эффект можно назвать неотложенным делением). Приоритеты симплекса S определяются на основе оценок поведения функций задачи в пределах симплекса. Оценки на множестве S строятся по имеющимся изме-
Л, 0 1 N
рениям функций в его вершинах V , V ,...^ . Таким образом, одновременно учитываются результаты N + 1 испытания, что должно приводить к более качественному прогнозу поведения функций, чем при использовании других компонентных схем.
Одним из центральных вопросов при построении SM-методов является разработка вычислительно эффективных способов определения характеристик-приоритетов симплексов для задач с различными свойствами. Выбор подходящего класса функций всецело определяет возможность и эффективность такого вычисления и является чрезвычайно важным. В данной работе приводится описание нового SM-метода для задач с ограничениями-неравенствами, в котором эта проблема решается. Метод построен в предположении, что функции задачи относятся к некоторому специальному классу, который далее будем обозначать С^1р (О). Предложенный класс функций оказывается особенно удобным при его использовании в SM-методах. Это связано с простой квадратичной (параболо-идной) формой минорант функций в пределах каждого симплекса S, построенных по измерениям в его вершинах. Методы, использующие такие миноранты, далее будем называть пара-болоидными и обозначать SMP-методы. Новыми в работе являются: использование указанного класса функций (при измерениях только значений функций задачи) в сочетании с триангуляционным подходом, а также построенный для данного класса алгоритмов метод, учитывающий функциональные ограничения. Предварительные результаты были представлены в [12].
1. Постановка задачи и построение минорант
Рассмотрим задачу (1)-(3) с ограничениями-неравенствами
I (У) ^ т1п, у е Y, (1)
У = {у е D : Vi = 1,...,т gl(у) < 0}, (2)
D = {у е ^ : а] < у] < Ь} (1 = 1,..., N)}. (3)
Для удобства обозначений объединим функции, входящие в постановку задачи, в вектор-функцию Q(y) с компонентами ч1(у) (7 = 0,...,т), где qo = I, а qi = gi ^ = 1,..., т).
Определим класс С^р(Е) как специальное
подмножество класса непрерывно-дифференцируемых функций с липшицевыми производными по направлениям. А именно, для скалярной
функции q из класса С^^) должна существовать константа L > 0 , чтобы V у1, у2 е D,
у1 Ф у2 при и = (у2 - у 1)/\\у 2 - у1 || выполнялось неравенство
\ дд(у2)/ди-дд(у')/ди \< L \\ у2 - у1 \\. (4) Будем предполагать, что функции qj (7 =0, ...,т) задачи (1)-(3) удовлетворяют (4) со своими значениями констант Ь, величины которых считаются неизвестными и подлежат оцениванию. При построении вычислительных методов принимается, что в любой точке области поиска D измерению доступны значения всех функций задачи, а их производные не измеряются.
Вопрос о виде нижней оценки функции q этого класса рассмотрен в [6; 12] в следующей постановке. Пусть выполнено N + 1 измерение q в наборе точек V0, v1,...,vN области D, образующих полноразмерный симплекс в . Выпуклую линейную оболочку этого набора - обозначим ее S (невырожденный выпуклый полиэдр размерности Ы) - также будем называть симплексом. Пусть Юдч - совокупность измерений
0 N.
V ,...^ функции q в вершинах симплекса:
“N+1 = “N+l(q,^ = {(^,ч): 1 = 0,...,N} (5)
Лемма 1. Пусть фунация ч принадлежит
алассу с известным значением аон-
станты Ь и для симплеаса S с D получены результаты измерений (5), тогда точная нижняя оценаа ч ~(у) значений ч(у) для у е S имеет вид параболоида:
ч (у) = ф(у, С,w) = С + 0.5 • Ь \\ у - w \\2,( 6) где параметры С = С(£>) = С(юлч) и w = w(S) = =w(юN+l) однозначно определяются из системы уравнений:
ф^1,С,w) = ч], (V1,ч])е“N+1 О' = 0,....^.(7) Доказательство приведено в [6; 12]. Нетрудно видеть, что полная нелинейная система (7) из N + 1 уравнения относительно неизвестных С, w1,...,wN сводится к решению линейной. Действительно, вычитая из всех уравнений (7) с 1 > 0 уравнение с 7 = 0, получим невырожденную линейную систему относительно неизвестных w:
.(В)
(9)
\Ь(w, V1 -Vй) = ч -ч° + 0.5Ь^((IV0 \\2 -\\ V1 \\2) [/• =1,...Д
Далее из уравнения с 7 = 0 в (7) находим С = ч0 -0.5Ь \\ V0 -w \\2.
2. Вычисление оценок констант класса функций
В (8), (9) необходимо иметь значение константы Ь класса функций (4). Поскольку Ь априори неизвестно, должно выполняться его оценивание. Для каждой функции чг-, i = 0,., т, должна быть построена отдельная оценка соответствующей константы Ь Метод оценивания значения Ь опишем на примере некоторой функции ч, опуская индекс.
Алгоритм оптимизации после выполнения очередного (к-го) шага будет вычислять и использовать несколько типов оценок для Ь: общую по всем симплексам разбиения текущую глобальную оценку I к, ее завышенное ступенчато изменяемое значение Ь , локальные оценки 1(5) для каждого симплекса 5 текущего разбиения, а также завышенные локализованные
(смешанные) оценки Ь(5) для каждого из симплексов. Для вычисления нижней оценки ч (у) в симплексе 5 из (6) в дальнейшем будет использована завышенная локализованная оценка, т.е. в (8), (9) будет принято Ь = Ь(5).
Все типы оценок строятся на основе получения локальных оценок 1(5). Вначале опишем способ их вычисления. При определении 1(5) в симплексах учитывается особенность их возникновения за счет деления симплексов-пред-ков. Каждый симплекс 5 поддерживает очередь из фиксированного количества локальных оценок 1Ь12,...,1Г, где выбрано г = N + 1. При делении симплекса 5 на две части каждая из них наследует от 5 копию его очереди оценок и осуществляет в этой копии сдвиг в сторону больших номеров: 4 := 4^ (у = г,...,2). Наиболее «старый» элемент 1г при этом удаляется из очередей и забывается. В освободившиеся элементы 11 очередей оценок новых симплексов помещается новая локальная оценка. Она вычисляется по измерениям функции, выполненным в трех точках, лежащих на делимом ребре: на его концах и в центре - точке нового испытания. В качестве 11 принимается модуль второй производной полинома Лежандра, построенного по этим трем измерениям. Окончательно, в качестве локальной оценки 1(5) для симплекса 5 принимается максимум элементов очереди этого симплекса и малого 0 < 5 << 1, гарантированно отделяющего локальную оценку от нуля:
l(S) = max{l1v.., lr, 5 }. (10)
Для начального заполнения очередей симплексов на шаге k = 0 используются начальные измерения функций (в вершинах D, его геометрическом центре и центрах граней размерностей от 1 до N - 1). Для каждой из триад точек измерений, размещённых на одной прямой, вычисляется (аналогично предыдущему) значение оценки l и добавляется в начало очередей тех симплексов S, которым принадлежат любые две точки использованной триады. Добавление происходит в элемент /1 с предварительным сдвигом существующих элементов очередей в направлении больших номеров. Незаполненные «старшие» элементы очередей обнуляются.
Общая текущая глобальная оценка l k вычисляется по правилу:
l *k = max {l(S,): i = 1,..., nk }, (11)
где nk — количество симплексов текущего разбиения на шаге k.
Поскольку значение lk, в общем случае, меньше неизвестной истинной величины L, вы-
*k 1—г
числяется завышенное значение L . Причем применяется специальный прием для уменьшения числа коррекций этой величины, исполь-
— * k
зующий вспомогательное L .
1. При k = 0 строится вспомогательная за-
1*k
вышенная оценка для l :
L * k=0 =Y1 . l * k=0 (1 <Y1). (12)
2. На последующих шагах организуется пересчет:
fLk, при l *k+1 < L *k
lY1 • l *k+1, при l *k+1 > L *k . (13) Заметим, что при yj = 1 правило (13) упро-
- * k + 1 7 * k + 1
щается до обычного равенства L = l
(величины изменяются одновременно), а при Y1>
Т* k fk
> 1 значение L изменяется реже, чем l .
L
* k+l
Построение завышенной глобальной оценки
*к
Ь зависит от двух параметров: номера шагаМ>> >> 1 и коэффициента надежности у2 > у1 - и выполняется следующим образом:
Lk = J(Y2/Ті ^L *k, при k — M,
\L*k + A ,
(14)
при k > М, где Д вычисляется на шаге k = М по правилу Д = (у2/у1 -1)- L М = LM - L М . Аналогично
скорректируем локальные оценки 1(5,) из (10), получив их завышенные аналоги Ь (5,):
г(5)=1(ї2/у1)-1(5^ при k<M,
( ') =[1(5,.)+Д-(?(5)/Ь*k), при k>м.(15)
Учитывая вариабельность константы L для разных подобластей области поиска D, метод использует в качестве ее оценки для симплекса
S смешанные (локализованные) значения L(S) , определяемые по двум завышенным значениям
(14), (15), путем их смешивания с весами, зависящими от диаметра симплекса d(S) =
( II > п\\ > >> О 'ї •
= max { у - у ■ у, у є S }■
Щ)=
L
L
*(Si)\l-dSi)ld*)+L*k •(d(Si)/dj, d(Si)—d*,
(16)
где й — пороговое значение диаметра, при котором начинается смешивание. Обычно й определяется как заданная доля 0 < ц < 1 от диаметра исходной области: й = ц • diam(D).
Описанный выше метод получения оценок констант зависит от следующего набора параметров: 0 < 5 << 1, 1 <у1 < у2, М >> 1, 0 <ц <1.
Заметим, что иные формы смешанных оценок констант классов функций широко использовались в [7] применительно к другим типам компонентных методов, основанным на разбиении множества D на гиперинтервалы.
3. Общая структура SMP-метода
Приведенное ниже описание является лишь принципиальной схемой алгоритма. При его программной реализации для уменьшения вычислительных затрат на шаге 2 не следует заново вычислять приоритеты всех симплексов, а только новых. Полный перерасчет необходимо выполнять только в моменты изменения общей
т* к
завышенной оценки Ь- одной из констант Ь- ,
а также при изменении достигнутой оценки fk глобально-оптимального значения целевой функции.
Принципиальная схема алгоритма
Шаг 1. Задаются параметры метода и параметр останова е, а также максимальное число итераций Мах_к. Область поиска D преобразуется к [0,1^. Начальные испытания проводятся в вершинах D, его геометрическом центре и в геометрических центрах граней меньших размерностей. Строится начальная триангуляция для D, где точки испытаний отвечают вершинам симплексов 5г-. Полагается к = 0. Заполняются очереди локальных оценок констант всех функций задачи qj для симплексов 5,, согласно (10) определяются локальные оценки 1-(5,), по (11), (12), (14) вычисляются для к = 0 глобальные
оценки констант
l
симплексов разбиения, рк=0 - число выполненных испытаний функций задачи.
Шаг 2. Определяется новое значение оценки
Р *
/к . Для всех функций задачи (ч-,- = 0,...,т) в каждом из симплексов 5- вычисляются локализованные оценки Ь- (5-) констант Ь- согласно описанному в разделе 2 алгоритму по формулам
(15), (16).
Шаг 3. Для компонент 5- триангуляции с использованием локализованных оценок констант класса Ь- (5-) вычисляются параметры минорант-параболоидов для функций задачи и на их основе строятся характеристики-приоритеты Л(5г). Определяется наиболее приоритетный симплекс 5Ь где
Я(5,) = тт {Я(5-): - пк }. (17)
Шаг 4. Если диаметр й(5) < е или число итераций к превысило заданное значение Мах_к, метод останавливается, иначе - переходит на шаг 5.
Шаг 5. Точка очередного испытания ук1 размещается в центре наибольшего ребра г/: симплекса St. Проводится испытание в выбранной точке с вычислением набора значений @ функций задачи, пополняется множество точек испытаний и их результатов.
Шаг 6. Компонента 5/ разделяется на две новые с помощью точки ук1. Отыскиваются все остальные симплексы 5 триангуляции, содержащие ребро г (эта операция выполняется без полного перебора компонент за счет включения в симплексы ссылок на соседей через грани размерности N - 1). Производится разделение
найденных 5 на две части с использованием в ^ к+1 /
качестве их новой вершины точки у (неотложенное деление). Возникает новая триангуляция области D.
Шаг 7. В новых симплексах вычисляются новые значения локальных оценок констант /,(5г), исходя из правил (11), (13), (14) определяется необходимость коррекции глобальных
k+l
оценок констант
l.
и их завышенных значе-
* k
ний Lj и Lj . Далее полагается k = k+1,
и их завышенные значе-
T~*k T*k r-,
ния Lj и Lj . Запоминается nk=0 - количество
г *k
корректируется число компонент щ и испытаний ріс , выполняется переход на шаг 2.
Заметим, что эффективная программная реализация алгоритма основана на использовании достаточно сложно организованных динамических структур данных [13]. Текущая экспериментальная реализация выполнена на основе левосторонних куч - приоритетных очередей, представленных в виде левостороннего бинарного дерева. В качестве ключей элементов кучи
Рис. 1. Построение начальной триангуляции: слева - базовая триангуляция для N = 2; справа - один из этапов рекурсивного построения триангуляции гиперкуба размерности п + 1 по триангуляции одной из граней размерности п (на примере п = 2)
Принцип деления симплексов: по точке нового испытания У+1 в центре большего ребра г, наиболее приоритетного (первично делимого) симплекса 5,; вторично делимые симплексы (отмечены цифрами 1, 2 и 3) включают общее с 5, ребро г,
использованы приоритеты Я(5) симплексов, связанных с этими узлами [14]. Для построения на шаге 1 начальной триангуляции множества D с ^ применяется рекурсивная по размерности п = 2,...Л процедура, основанная на использовании для граней размерности два базовой триангуляции по девяти точкам (см. рис. 1). Общее число порождаемых начальных испыта-
~^N
ний р0 = 3 , а число симплексов начального разбиения п0 = N^2N. Включение в схему триангуляции центральных точек для граней любых размерностей необходимо для применения метода оценивания констант Ь-, описанного в разделе 2. Возможно использование более экономичных алгоритмов триангуляций [15], но это потребует иного способа оценивания констант.
В качестве дополнительной информации к описанию алгоритма (шаги 4-6) укажем на известную оценку (см., например, [16]) скорости уменьшения диаметров последовательности { 5(п) }п=0 взаимно вложенных симплексов, получаемых в результате их деления «пополам» точкой в середине наибольшего ребра делимого симплекса:
й(5(п)) < (Тз/2) ^^ • й(5(0)), где п - количество выполненных последовательных делений, а \rnlN] - наибольшее целое,
не превосходящее пШ.
Поддержка операции неотложенного деления «соседних» с 5/ симплексов на шаге 6 может приводить к уплощению части симплексов и образованию значительного количества новых. Однако этот механизм способствует скорейшему учету информации о проведенном испытании не только в наиболее приоритетном (делимом первично) симплексе 5/, но и во всей
группе вторично делимых симплексов, содержащих разделяемое ребро г/ (см. рис. 2).
Приведенная принципиальная схема алгоритма требует конкретизации способа вычисления приоритетов-характеристик Я(5) симплексов 5. Этот вопрос непосредственно связан с выбором метода учета функциональных ограничений.
4. Редуцирование задачи условной оптимизации
Перейдём к вопросу о вычислении приоритетов симплексов при наличии ограничений-неравенств (2). В [6; 12] для задачи (1)-(3) на рассматриваемом классе функций были предложены два подхода, приводящие к двум разным методам.
Один из них, исследуемый в данной статье, основан на редуцировании задачи с функциональными ограничениями к задаче без функциональных ограничений с перестраиваемой целевой функцией и является универсальным. Он применим, так же как и известные индексные методы [17; 18], даже в той ситуации, когда допустимая область исходной задачи (1)-(3) пуста. В последнем случае решение понимается в расширенном смысле, как глобальный минимум невязки в нормированных ограничениях исходной задачи.
Пусть / - глобально-оптимальное значение целевой функции в задаче (1)-(3), если допустимая область не пуста, и / = +<», в случае ее пустоты.
Заменим неизвестное / его оценкой. Пусть
*
после к испытаний fk — наименьшее вычисленное значение целевой функции в допустимой области У, а если таких значений нет, то
*
]к = +<», т.е.
^ _jmmf(ys):s _ 1,... k,/ еГ},если{у1,...,yk}f|Y* 0,^ k l+K, в противном случае.
Введем для функций задачи нормирующие множители Lt (i _ 0,...,m), выбираемые на первом этапе поиска решения задачи (при k < M < K) равными оценкам L tk из (13). Эти множители
в процессе дальнейшего поиска решения (при k>
> K) не будут изменяться.
Построим редуцированную задачу без функциональных ограничений с перестраиваемой целевой функцией:
^) (У) — min У e D, n(k) _ fk*, (19)
р1(k)(y) _ max{(f (y) - n(k»/~0;G(y)}, (20)
G(y) _ max {gt(у)/1, : i _ 1,...,m }. (21)
Легко видеть, что при замене n(k) значением f редуцированная задача (19)-(21) эквивалентна исходной (1)—(3), решение которой понимается в расширенном смысле.
В разделе 5 будет конкретизирован метод решения задачи (19)—(21), существенно учитывающий информацию об особенности ее структуры и опирающийся на поведение минорант всех функций исходной задачи.
Далее в текущем разделе будут рассмотрены условия, обеспечивающие сходимость процесса редуцирования к глобальному минимуму исходной задачи. Для упрощения записи в (20),
(21) будем опускать множители Lt (i _ 0,...,m),
а также знак тильды над функцией G , поскольку это несущественно для проводимого анализа.
Временно абстрагируемся от конкретной
k+1
процедуры порождения точек у поисковой последовательности и наложим формальные ограничения на их выбор, определяемые положениями A, B.
Положение А. Пусть выбор очередной точки испытания yk+1 удовлетворяет следующим условиям. Для функции P^(k )(у) _ max {f (у) -
--n(k);G(у)} где G(у)_max{g,-(у):j_1,...,m},
выполнено:
P,(k)(yk+1) <P,(k)(yk), (22)
P„(k) (yk+1) ^ min{pn(k) (у): у e d}+ sk, sk > 0 ,(23)
lim s k _ 0 при k —— да, ,(24)
где r[(k +1) _ f (у^1), если уk+1 e Y, и n(k+1) = = да в противном случае.
Обратим внимание, что в силу (22) последовательность n(k) невозрастающая.
Положение В. Если допустимая область
Y ^ 0, то cov(intY) _ Y и Vу eintY: G^) < 0.
Заметим, что из положения В непосредственно следует, что в сколь угодно малой окрестности любой допустимой точки найдется допустимая точка у е У , в которой G(у) < 0 .
Теорема 1. Пусть для исходной задачи выполнены положения В, функции в (1)-(3) по крайней мере непрерывны, а выбор точек ук+1 последовательности испытаний удовлетворяет положениям А, тогда при У ^ 0 предельными точками этой последовательности будут являться только точки глобальных минимумов у * е У *, если же допустимая область У в (2) не пуста, предельными будут точки глобальных минимумов невязки в ограничениях G(y).
Доказательство. Рассмотрим случай У = 0, тогда />Л(к) (у) = G(y). В силу компактности D и непрерывности G(y) существует и достигается ее конечное глобально-оптимальное значение G на D. Из (22) и (23) следует монотонная сходимость G(у к) к G , что в рассматриваемых условиях означает, что все предельные точки последовательности ук глобально оптимальны для функции невязки.
Рассмотрим случай, когда допустимая область в (2) не пуста. Докажем вначале несколько лемм.
Лемма 2. Если у допустима и ук+1, удовлетворяющая (22), существует, то у к+1 е т; У и
РЛ(к)(ук+1) < 0.
Действительно, из допустимости ук следует, что G(ук) < 0, а "л(к) = f (ук), поэтому РЛ(к)(ук) = тах{ 0; 0(ук) }< 0 . Таким образом, из (22) получаем Р^к)(ук+1) < 0, что возможно только при G(у к+1) < 0, следовательно,
у к+1 е Ы У . Лемма доказана.
Лемма 3. Пусть для допустимой области задачи (1)-(3) выполнены положения В. Если у допустимо и значение Г[(к) = f (ук) отделено от глобально-оптимального f (у*) так, что "Л(к) > f (у*) + й, й > 0, то УДе (0, й) Зу е D и те (0,Д/4), что е[/(у*) + Д,"л(к)] выполнится неравенство Рл (у) < .
Доказательство. Выберем произвольное Д из условия 0 < Д< й. Найдется у е У, что
f (У) < f (у *) + Д/ 2. Поскольку у - допустима, то G(y) < 0 . Из положений В следует существование (в сколь угодно малой окрестности уу ) внутренней для У точки у, удовлетворяющей условиям:
G(у) < -т < 0, /(у) < /(у*) + Д/2 + 8(т), (25) где 5(х) > 0 и 5(х) ^ 0 при т ^ 0, а положительная величина т выбрана настолько малой, что т < Д/4 и 5(т) < Д/4.
Пусть п выбрано из промежутка, указанного в условиях леммы, тогда из (25):
/(у) - Л < /(у*) + Д/2 + 8(т) - (/(у*) + Д) =
= 8(т)-Д/2 <-Д/ 4. .
Учитывая ранее наложенные ограничения на величину т, получаем, что значение
Рл (у) < тах(-Д/4, -т) < -т, что и требовалось.
Лемма 4. Если у допустима, то либо ук еУ *, т.е. является глобально-оптимальной, либо существует ук+1, удовлетворяющая (22), т.е. данное правило корректно определено.
Доказательство. Пусть ук €У *, ук еУ, тогда, применяя лемму 3 с й = /(ук) - /(у ) > 0 и
учитывая, что Рп(к)(ук) = 0, непосредственно по-
к+1
лучаем существование точки у , удовлетворяющей условию (22). Доказательство завершено.
Леммы доказаны. Продолжим доказательство теоремы для случая непустой допустимой области. Заметим, что до обнаружения первой допустимой точки Рп(к)(у) = G(y) и на этой стадии поиска, в силу положений А , будет выполняться оценка G(yk) < G* + ек при ек ^ 0 с ростом к. Но поскольку, в силу положения В, будет выполнено неравенство G < 0, то на некотором шаге будет найдена допустимая точка ук. После этого, в силу леммы 2, все последующие точки этой последовательности останутся допустимыми, и, если еще л(к) Ф f , в силу леммы 4, выбор каждой следующей точки по правилу (22) будет корректно определен. Заметим, что
(22) означает одновременное выполнение условия (23) при некотором ек > 0. Пусть теперь, в силу условий теоремы, используемый метод выбора ук1 обеспечивает выполнение (24).
Заметим, что из леммы 2 следует
рл(к)(ук+1) < 0 , поэтому f (ук+1) - Л(к) < 0. Если ещё учесть допустимость рассматриваемых к+1 ^
точек у , получим двойное неравенство:
/(у *) <л(к +1) <л(к).
Таким образом, последовательность п(к) сходится к некоторому значению п .
Предположим, что л > /(у ). Применяя лемму 3 при й = л* - /(у*) > 0 и Д = й/2, получим, что существует точка у е Б и значение 0 < т < Д/4 , что для всякого к (поскольку из
Л* < л(к) следует л(к) е [/(у*) + Д, л(к)]), будет выполняться
Рл(к)(у) <-т. (26)
Фиксируя значение т, выберем Е: 0 < £ < т/2. Для него найдется К(£), что У к > К (Е): 0 < п(к) -- п* < £ и, следовательно, 0 < п(к) - п(к+1) < £. Учитывая правило выбора л (к +1), получаем,
что /(ук+1) -л(к) >-Е , следовательно, при всех достаточно больших к выполнится неравенство:
Рл( к)(ук+1) >4>-т/ 2.
Тогда из (23) следует, что при всех достаточно больших к
уеБ1>-т/2-ек >-т. (27)
С учетом (26), (27) приходим к противоречивому утверждению: З у е Б, что при больших к
Рл(к) (у) < ^Кк)(у): у е °).
Это доказывает, что л = /(у ), следовательно, / (ук) ^ / (у ) при к ^ да и допустимости точек ук, что определяет глобальную оптимальность предельных точек последовательности испытаний. Теорема доказана.
5. Определение приоритетов симплексов редуцированной задачи
Для построения вычислительного метода решения редуцированной задачи (19)-(21) используем текущую триангуляцию области поиска Б. Пусть 5 - один из симплексов. В пределах 5 заменим все функции qj (у) (и = 0,...,т) исходной задачи, входящие в (20), (21), на их нижние оценки qИ (у, 5), имеющие вид параболоидов (6). При этом для каждой функции на текущем шаге к будем использовать свои смешанные оценки Ьк(5) констант Липшица (16) для производных по направлениям. Получим вспомогательную оптимизационную задачу на симплексе 5. Оптимальное значение в этой задаче является оценкой минимума функции (20) на 5. Эту оценку (с учетом нормирующих делителей Lj, полагаемых в реализации метода рав-*к
ными Ь' ) примем в качестве значения характеристики-приоритета ^(5).
Получим следующее описание способа вычисления R(5):
Я(5) = шт {тах{с°(5) + 0.5 • М 0(5) х
х11 у- ™0(5)( ) ^-(у)}: уе 5} ( )
где
G- (у) = шах{С' (5) + 0.5 • М} (5) х
II и 2 (29)
х у - wи (5) : И = 1,..., т },
С0(5) = (С0(5) - /к*)/Ь, С (5) = Си (5)/Ь
(30)
(И = 1,...5 т)
М] (5) = Ьки (5)/~ (И = 0,..., т). (31)
Условие у е 5 эквивалентно системе из N + 1 линейного неравенства
(у -у ^ а) < 0 0' = и.N) (32)
(у-у7",а0) < 0 , ( )
0 N о 0 N
где V ,^,у — вершины 5, а а ,...,а — внешние нормали к граням размерности N - 1 симплекса
5, противолежащим вершинам V0,... ,^. Каждый
' N
из векторов а, например а , вычисляется с помощью процедуры, подобной ортогонализации
с- 10 N 0
базиса V - V ,..., V - V .
Для вычисления R(5) задача (28)-(32) эквивалентно приводится к задаче сепарабельного программирования за счёт введения искусственной переменной z е R1 :
Я(5) = шт 2, (33)
°.5 • Л~0 ^|у - W0of + С00 < 2;
0.5 • Ми -||у - + СИ < 2 (И = 1,..., т),(34)
(у-у0,а) < 0 (' = 1,...,N),
(у-уN,а0) < 0
(35)
где в целях сокращения записи не указана зависимость величин С', и', М' от 5.
Все функции в (33)—(35) сепарабельны и выпуклы по у, г. Если бы целевая функция была строго выпуклой, то в силу известной теоремы (см., например, [19] ) можно было бы приблизить эту задачу к аппроксимационной кусочнолинейной, эквивалентной некоторой задаче линейного программирования по отношению к специально вводимым новым переменным. В нашем случае данная теорема неприменима, но в [6; 12] доказан ее аналог, в силу которого указанная технология определения R(5) все же может быть использована из-за специфики структуры задачи (33)-(35). Аппроксимационная задача строится по значениям функций в (34), (35) от переменных у, вычисляемых в узлах сетки из
N
п^а точек, покрывающей наименьший гиперинтервал [а(5), й(5)] (в исходной системе координат), содержащий симплекс 5.
Таким образом, вычисление R(5) приближенно сводится к решению специальной задачи линейного программирования, вид которой приведен в [6; 12].
Заметим, что в описании SMP-метода условной оптимизации в разделах 2, 3, 5 нумерация точек испытаний ук не совпадает с нумерацией, принятой в разделе 4. А именно, очередной ите-
рации, как она понимается в разделе 4, соответствует момент улучшения текущего рекордного значения целевой функции, достигаемого в задаче (19)—(21) при выполнении вычислений по SMP-методу.
Отметим также, что анализ выполнения положения А при использовании SMP-метода выходит за пределы данной статьи и в последующем будет проведен дополнительно.
6. Результаты вычислительных экспериментов
Приведены результаты предварительных вычислительных экспериментов для размерности N = 2 на трех тестовых задачах, описания которых представлены в табл. 1, 2. В пятой колонке табл. 2 указано количество локальных
минимумов в задаче с учетом ограничений, а
0
также величина / , характеризующая наименьшее локально-оптимальное значение целевой функции, отличное от глобального. Параметры глобально-оптимальных решений, приведенные в шестой колонке этой таблицы, получены с использованием локальных методов пакета MatCad.
Для каждой из задач были проведены вычислительные эксперименты при нескольких значениях параметров SMP-метода с ограничениями. Использованные наборы параметров и результаты вычислений представлены в табл. 3.
Параметры расчетов в табл. 3 соответствуют обозначениям, введенным в разделах 2, 3. В блок данных о процессе решения включены следующие показатели: р^ - общее число испытаний к моменту выполнения критерия оста-
У
нова и рк(е) - число испытаний к этому моменту в допустимой области; ик(е) - число симплексов в построенной триангуляции к моменту останова; к - число испытаний до момента нахождения первой допустимой точки, со значением целевой функции, меньшим /0 (с этого момента задача правильного оценивания глобального минимума может считаться решенной); Р - число полных перевычислений приоритетов симплексов текущих триангуляций.
В блоке оценок решения табл. 3 подчеркнуты разряды, значения в которых определены верно. Реальная точность определения глобально-оптимального решения в большинстве проведенных расчетов оказывается ниже заданного значения параметра останова е = 10-4. В то же время во всех случаях задача определения глобального минимума оказалась корректно решена за число испытаний к0 от 49 до 75 для задач 1 и 2, а для задачи 3 (с малым значением разности /0 -/ ~ 0.00097) за число испытаний от 134 до 153.
Таблица 1
Набор функций, использованных в вычислительных экспериментах
Обозначение функции Формула Источник
/і - І.5 у2 ехр(1 - у2 - 20.25(у - У2)2) - - (0.5(уі - 1)(у2 - І))4ехр(2 - (0.5(уі - І))4 - (у2 - І)4) [11]
/2 (4 - 2.І УІ + Уі4/3) УІ + У1У2 + (4 У2 - 4) У2 [1]
/з 0.0І (уіУ2 + (Уі - п)2 + 3(У2 - п)2) - ^іп(уі>іп(2у2))2 [1], с модификацией
/4 0.00і((у1 - 2.2)2 + (у2 - І.2)2 - 2.25) [11]
/5 100(1 - ((Уі - 2)/1.2)2 - (О.5у2)2) [11]
/6 10 (у2 - 1.5 - І.5sin(2п(yІ - 1.15))) [11]
Її - (1.5уІ -у2 - 0.2)2- (2sin(2y2) + 0.2)2 + 1 —
/8 - 14 + уІ + О.І|3 + 2|у2 - 0.2|3 —
/9 - (уІ - п + О.І)2 - ^іп(у2) + О.2)2 + 6 —
Таблица 2
Описания тестовых задач
Задача Постановка Функции Область поиска Б Количество локальных минимумов; / 0 Глобальный минимум / *
Задача 1 тіп /(у), gi(y) < О (і = 1, 2, 3), у є D /=/1; gl = и ^ = А ^ = /6 [ О, 4 ] [-1, 3 ] 5 -І.ООООО / * = - 1.48968 (О.94248; О.94526)
Задача 2 тіп /(у), gi(y) < О (і = 1, 2), у є D /=/2; Й = /1, g2 = /8 [-2.5, 2.5 ] [-1.5, 1.5 ] 1 -О.2І546 / * = - О.8О461 (-О.3252О; О.18І91)
Задача 3 тіп /(у), gi(y) < О, у є D /=/3; gi = /9 [ О, 2п] [ О, 2п] 8 -О.8І8І4 / * = - О.8І9ІІ (І.3О499; 2.21249)
Таблица 3
Результаты вычислительных экспериментов
№ Задача Параметры Оценка решения Данные о процессе решения
Ті Ї2 м є / Л к (є) У*(є) р«.) ( р«.) ) пкє) кО Р
1-І 1 1 2.5 ІОО О.25 ІО'4 -1.4896 О.942О; О.9442 161(52) 319 58 23
1-2 1 1 3.О ІОО О.25 ІО'4 -1.4896 О.942О; О.9442 114(53) 332 64 25
1-3 1 1.2 3.О ІОО О.25 ІО'4 -1.4896 О.942О; О.9442 112(51) 328 68 21
1-4 1 1.5 3.О ІОО О.25 ІО'4 -1.4891 О.9424; О.945І 116(53) 336 68 19
1-5 1 1.5 3.О 50 О.25 ІО'4 -1.4894 О.94І2; О.94І2 131(30) 258 61 18
1-6 1 1.5 3.О 50 О.І25 ІО‘4 -1.4894 О.94І2; О.94І2 195(50) 314 84 18
1-1 1 1.5 3.О І5О О.25 ІО‘4 -1.4891 О.9426; О.9458 212(11) 4О8 69 20
1-8 1 1.5 3.О 5ОО О.25 ІО‘4 -1.4891 О.9426; О.9458 212(11) 4О8 69 21
2-І 2 1 2.5 ІОО О.25 ІО‘4 -О.8О46 -О.3253; О.18І9 225(128) 42О 16 34
2-2 2 1 2.О ІОО О.25 ІО‘4 -О.8О46 -О.3253; О.18І9 198(111) 311 11 34
2-3 2 1.5 3.О ІОО О.25 ІО‘4 -О.8О46 -О.3253; О.18І9 218(115) 4О5 69 19
2-4 2 1.5 3.О 50 О.25 ІО‘4 -О.8О46 -О.3253; О.18І9 202(115) 331 69 11
2-5 2 1.5 3.О 5ОО О.25 ІО‘4 -О.8О46 -О.3253; О.18І9 216(115) 4ОІ 69 22
3-О 3 1 3.О ІОО О.25 ІО‘4 -О.8І9О І.3О48; 2.2126 115(92) 319 159 14
3-І 3 1.5 3.О ІОО О.25 ІО‘4 -О.8І9О І.3О5І; 2.2122 168(89) 305 151 11
3-2 3 1.5 3.О ІОО О.5О ІО‘4 -О.8І9О І.3О5І; 2.2122 146(18) 262 143 11
3-3 3 1.5 3.О 50 О.25 ІО‘4 -О.8І9О І.3О5І; 2.2122 161(89) 292 151 11
3-4 3 1.5 3.О 5ОО О.25 ІО‘4 -О.8І9О І.3О5О; 2.2124 181(96) 331 162 14
Рис. 3. Размещение испытаний и вид триангуляции в эксперименте 1-4 для задачи 1
Рис. 4. Размещение испытаний и вид триангуляции в эксперименте 2-3 для задачи 2
Рис. 5. Размещение испытаний и вид триангуляции в эксперименте 3-2 для задачи 3
На рис. 3-5 показаны характерные виды размещения точек испытаний при решении каждой из трех тестовых задач, а также построенные триангуляции. Более темным цветом выделена допустимая область. Приведенные рисунки демонстрируют адаптивный характер проведения испытаний и получаемых триангуляций в зависимости от задачи. Следует обратить внимание на присутствие эффекта уплощения симплексов за счет механизма вторичных делений.
Общее число испытаний, проведенных методом до останова, в выполненных экспериментах колебалось от 131 до 225, что характеризует экономичность метода (табл. 3). Повышение коэффициента надежности у2 (завышения оце-
нок констант Lj), а также увеличение параметра М приводит к некоторому росту затрат по числу испытаний и несколько повышает точность определения решения к моменту останова. Следует отметить, что существенное повышение реальной точности определения решения самим SMP-методом за счет увеличения значений параметров у2, М представляется нецелесообразным. Более эффективным является использование процедур локальной оптимизации на стадии уточнения найденной оценки глобальнооптимального решения.
Результаты в колонке Р табл. 3 показывают возможность (за счет использования значений Уі > 1) снижения числа наиболее вычислительно
затратных операций - полных пересчетов приоритетов симплексов без заметного снижения качества работы метода.
Заключение
В статье рассмотрен новый компонентный метод многоэкстремальной оптимизации для задач с ограничениями-неравенствами, основанный на адаптивной триангуляции области поиска по точкам измерений. Предложена схема сведения исходной задачи к задаче без функциональных ограничений с перестраиваемой целевой функцией, имеющей вид верхней огибающей набора ограничений и целевой функции исходной задачи с регулируемым смещением. Получены условия сходимости процесса редуцирования. В случае пустоты допустимой области метод определяет глобальный минимум невязки в нормированных ограничениях. Показано удобство использования класса функций с липшицевыми производными по направлениям в качестве модели поведения многоэкстремальных функций в симплексе (в предположении измеримости только значений функций): для этого класса миноранты функций имеют в каждом симплексе вид параболоидов. Планирование размещения очередных измерений основано на вычислении приоритетов симплексов. С учетом вида минорант это вычисление сводится к численному решению специальных задач выпуклого сепарабельного программирования, приближенно приводимых к задачам линейного программирования. Представлены результаты предварительных вычислительных экспериментов, показывающие экономичность построенного метода по числу проводимых испытаний.
Автор выражает благодарность выпускнику ВМК HH^ А.О. Парфенову за помощь в подготовке экспериментальной версии программной реализации метода.
Работа выполнена в рамках ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы, госконтракт М 02. 740.11.5018.
Список литературы
1. Евтушенко Ю.Г., Ратькин В.А. Метод ноло-винных делений для глобальной оптимизации функций многих переменных // Техническая кибернетика. 1987. № 1. С. 119-127.
2. Pinter J. Global optimization in Action. Dordrecht: Kluwer Academic Publishers, 199б. 478 p.
3. Городецкий С.Ю. Многоэкстремальная оптимизация на основе триангуляции области // Вестник Нижегородского государственного университета. Математическое моделирование и оптимальное управление. Н.Новгород: Изд-во ННГУ, 1999. № 2 (21). С. 249-269.
4. Городецкий С.Ю. Методы многоэкстремальной оптимизации на основе триангуляции области поиска // Первая Всероссийская научно-практическая конференция по вопросам решения научно-практических задач в промышленности ОПТИМ-2001. Сборник докладов. СПб.: ЦНИИТС, 2001. С. 191-196.
5. Сергеев Я.Д., Квасов Д.Е. Адаптивные диагональные кривые и их программная реализация // Вестник Нижегородского университета. Математическое моделирование и оптимальное управление. Н.Новгород: Изд-во ННГУ, 2001. № 2 (24). С. 300317.
6. Городецкий С.Ю., Гришагин В.А. Нелинейное программирование и многоэкстремальная оптимизация. Н.Новгород: Изд-во ННГУ, 2007. 489 с.
7. Сергеев Я.Д., Квасов Д.Е. Диагональные методы глобальной оптимизации. М.: Физматлит, 2008. 352 с.
8. Clausen J., Zilinskas A. Subdivision, Sampling, and Initialization Strategies for Simplical Branch and Bound in Global Optimization // Computers and Mathematics with Applications. 2002. № 44. P. 957-967.
9. Zilinskas J. Branch and bound with simplicial partitions for global optimization // Math. Model. Anal. 2008. V. 13, № 1. P. 145-159.
10. Paulavicius R., Zilinskas J., Grothey A. Investigation of selection strategies in branch and bound algorithm with simplicial partitions and combination of Lip-schitz bounds // Optim. Lett. 2010. № 4. P. 173-183.
11. Елсаков С. М., Ширяев В.И. Однородные алгоритмы многоэкстремальной оптимизации // Журн. ВМ и МФ. 2010. Т. 50. № 10. С. 1-14.
12. Городецкий С. Ю. SM-методы в задачах многоэкстремальной оптимизации с ограничениями для класса функций с липшицевыми производными по направлениям. Н.Новгород: ННГУ, 2007. 22 с. Деп. в ВИНИТИ 11.01.07, № 18-В2007.
13. Алексеев В.Е., Таланов В.А. Графы. Модели вычислений. Структуры данных. Н.Новгород: Изд-во ННГУ, 2005. 307 c.
14. Городецкий С.Ю., Осминин М.К. Использование динамических структур типа куч в многомерных триангуляционных методах многоэкстремальной оптимизации // Матер. Всерос. конф. «Технологии Microsoft в теории и практике программирования», Н.Новгород, 13-14 мая 2010 г. Н.Новгород: ННГУ, 2010. С. 89-92 [Электронный ресурс]. - Режим доступа: http://www.agora.guru.ru/ display.php? conf=msconf-2010.
15. Шевченко В.Н., Груздев Д.В. Модификация алгоритма Фурье-Моцкина для построения триангуляции //Дискретный анализ и исследование операций. 2003. Сер. 2. Т. 10. № 1. С. 53-64.
16. Horst R. On generalized bisection of N-simplices //Mathimatics of Computation. 1997. V. 66. № 218. P. 691-698.
17. Стронгин Р.Г., Баркалов К.А. О сходимости индексного алгоритма в задачах условной оптимизации с е-резервированными решениями // Математические вопросы кибернетики. М.: Наука, 1999. С. 273-278.
18. Баркалов К.А., Стронгин Р.Г. Метод глобальной оптимизации с адаптивным порядком проверки ограничений // Журн. ВМ и МФ. 2о02. Т. 42. № 9. С. 1338-1350.
19. Базара М., Шетти К. Нелинейное программирование. Теория и алгоритмы. М.: Мир, 1982. 583 с.
PARABOLOID TRIANGULATION METHODS IN SOLVING MULTIEXTREMAL OPTIMIZATION PROBLEMS WITH CONSTRAINTS FOR A CLASS OF FUNCTIONS WITH LIPSCHITZ DIRECTIONAL DERIVATIVES
S.Yu. Gorodetsky
A new method is proposed to solve multiextremal problems with nonconvex constraints for a class of functions with Lipschitz directional derivatives. Only function values are measured. Adaptive triangulation of the search region with irregular simplexes is applied. Minorants of functions in simplexes are constructed as paraboloids. The method reduces the problem with constraints to a problem on a hyper interval with a configurable objective function. Convergence conditions for the reduction process have been obtained. Computational experiment results are presented.
Keywords: multiextremal optimization, nonconvex constraints, partition algorithms, triangulation, simplexes, lipschitzian directional derivatives.
