ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2013 Управление, вычислительная техника и информатика № 1(22)
УДК 517.977
А.И. Рубан
МЕТОД ГЛОБАЛЬНОЙ ОПТИМИЗАЦИИ,
ОСНОВАННЫЙ НА СЕЛЕКТИВНОМ УСРЕДНЕНИИ КООРДИНАТ, ПРИ НАЛИЧИИ ОГРАНИЧЕНИЙ
Изложены идеи конструирования алгоритмов недифференцируемой глобальной оптимизации, в основе которых лежит: разнесение во времени пробных и рабочих шагов, селективное усреднение координат по результатам пробных движений, адаптивная перестройка размеров прямоугольной области пробных движений и учёт ограничений типа неравенств и типа равенств.
Ключевые слова: глобальная оптимизация, селективное усреднение координат, ограничения типа неравенств и равенств.
Среди существующих методов поиска глобального минимума перспективны рандомизированные подходы. Они часто построены на сочетании случайного просмотра области поиска X и локальных детерминированных алгоритмов поиска минимума.
В отличие от алгоритмов метода стохастической аппроксимации существенным продвижением при конструировании алгоритмов глобальной оптимизации является осознание того факта, что для более гарантированного движения к глобальному экстремуму необходимо проводить усреднение оптимизируемой функции во всей заданной области, где нужно отыскивать экстремум. На этом базируются подходы А. А. Красовского, А. И. Каплинского, А. И. Пропоя. Авторы используют потенциальные функции.
Другой способ отыскания глобального экстремума основан на селективном усреднении (также во всей заданной области на начальном этапе поиска) не самой оптимизируемой функции, а искомых координат. Такой способ предложен в работах [1-3]; его развивает и автор [4-11] данной статьи. Производится параллельный расчёт всех искомых переменных.
Решаем задачу поиска глобального минимума функции при изменении непрерывных переменных внутри допустимой области:
f (x) = min. (1)
xeX
1. Метод синтеза алгоритмов глобальной оптимизации
Вводим последовательность {ps(y), s = 1,2,...} непрерывных, положительных в R+ функций, таких, что для любых y < z (где y, z с R1) последовательность {Ps (У)/Ps (z), s = 1, 2,...} монотонно нарастает с ростом s :
lim ps (y)/ps (z) = ж; y < z . Считаем, что известны точная верхняя
fmax = sup f (x) и точная нижняя fmin = inf f (x) границы минимизируемой
xeX xeX
функции f (x).
Справедлив следующий результат:
Ps (gmin( -)) _ .ч f (-) - fm
Xmln =
limf xps,min(—)d-, Ps,mn(x) = ^s^ V gm,n(x) = . (2)
X J Ps (gmin (t)) dt fmax - fmin
X
Безразмерная неотрицательная переменная gmin( -) лежат в интервале [0; 1].
Заметим, что в (2) можно и не переходить к безразмерной переменной gmin (-). Достаточно положить gmin (-) = f (-). Предельные свойства сохраняются. Использование безразмерных переменных придает алгоритмам универсальность.
В качестве ядер ps (g) могут использоваться, например: 1) экспоненциальное
в степени s : exp(-sg), 2) гиперболическое в степени s : g~s, 3) линейное
(r =1), параболическое (r =2), кубическое (r =3) и т. д. в степени s: Ps(g) = (1 -gr)s, r = 1, 2, 3, •••. Эти виды ядер конструируются по единому образцу. За основу берется убывающая в интервале [0; 1] функция p(g) (с максимумом в точке g = 0), а затем эта функция возводится в степень s:
Ps (g) = (p(g))s. Для таких p(g) всегда 1 < (p(y)/ p(z)) при y < z . Следовательно, ядра ps (g) удовлетворяют основному условию для ядер.
Перечисленный набор ядер ps (g) пока достаточен для решения задач глобальной минимизации. При увеличении s растет «селективность» (способность к локализации положения глобального экстремума) нормированного ядра ps min (-). В пределе (при s ) ядро стремятся к многомерной дельта-функции с особой точкой - = -min . Усреднение в правой части формулы (2) позволяют предсказывать положение глобального экстремума при фиксированной степени s , которую будем называть «степенью селективности ядра».
Приведенный результат позволяет построить целый спектр практически реализуемых быстросходящихся алгоритмов, обеспечивающих поиск глобального экстремума с вероятностью, близкой к единице. Укажем возможный вариант класса алгоритмов (при различных параметрах 0 < у , q е {1, 2, • • •} , 0 < s и различных формах ядра) поиска глобального минимума:
( V/q
xl+І =
Xv - XV 1 q P. ,mln (X) dx
v = І, m , (3)
X1 = П1 пX , П1 =[х!-Ат!, X +Ат!]х---х[х1т-Ах1п, х1т + Ах1т],
I = 0, 1, 2, •••; [0 < уд, д е {1, 2, •••}, 0 < 5].
Здесь I - номер итерации; уд, д, 5 - подбираемые фиксированные параметры; П1 - прямоугольная область в окрестности точки х1. При 0 < у < 1 всегда Ах^1 < удАх(,, V = 1, т . Отсюда следует конечно-шаговое (и достаточно быстрое) сжатие области поиска до заданных малых размеров.
Алгоритмам (3) можно придать более привычный вид и осуществить в них приближенное вычисление интеграла. В области X, если она задана ограничениями неравенствами и не является узким «жгутом», размещается п пробных точек х(г), . = 1, п . В них вычисляется минимизируемая функция /(г) = /(х(г)), . = 1, п , и операция интегрирования заменяется суммированием.
При получении пробных точек х(г), . = 1, п , последовательно генерируются равномерно распределённые точки в прямоугольной области П1 с центром в точке х1 :
х^ = х(, + Ах(, • ««, иу е [-1; 1], V = 1, т, . = 1, 2, — , (4)
и из них оставляется п точек, попадающих в допустимую область X. Пробные точки лежат теперь в области X1 =Пг П X .
Исходная точка х0 и размеры Ах0 прямоугольной области П0 выбираются так, чтобы П0 охватывала допустимую область X или ту её часть, где расположен искомый глобальный экстремум.
В результате указанного перехода от интегрирования к соответствующему суммированию алгоритмы (3) приобретают следующий вид:
V/ ч
XV +1 1 ■ 1 — ■ 1 + 1 ■ 1 I ^ . (.') ч — (.')
I п у/4 _
^ ^ • ^,шт , ^+1 = Уч • ^ • I X I «V ) I , V = 1 т , (5)
11П =Х^г)Й
-(.) = Рз („(.) = /(0 - /т
5,Ш1П
. =1
X Рз(^ П)
1=1
/пах /т
I = 0, 1, 2, •••, [0 < уч, ч е {1,2,•••}, 0 < з].
Здесь /тах = тах{/(.),.=1,п}, /т1п = тт{/(0,. =1,п}; в переменных е[-1;1], Р«п для упрощения записи опущен номер итерации I. Весовые коэффициенты
п
(ядра) р«п нормированы на системе п пробных точек: XР3т1п = 1.
. =1
Прекращение работы алгоритмов можно осуществлять, как обычно, либо по заданной величине размера области пробных движений: тах{| АxV |, V = 1, т} < е1, либо по величине наибольшего уклонения минимизируемой функции на множестве пробных точек: /тах - /тт < е2.
В полученных алгоритмах пробные и рабочие движения разнесены. На каждом рабочем шаге осуществляется переход в новую точку, в среднем более близкую к глобальному минимуму, и производится изменение размеров области поиска.
Тестирование большого числа многоэкстремальных функций показало, что алгоритм (5) имеет устойчивые достаточно высокие окончательные результаты (по скорости сходимости, точности и близкой к 1 оценке вероятности отыскания глобального экстремума) при следующих условиях:
1) 0,8<уд < 1,2; д = 2;
2) степенных типах ядер: экспоненциальном, гиперболическом, линейном, параболическом, с целочисленными параметрами селективности 5, лежащими в широких пределах;
3) 50 < п ;
4) равномерном распределении пробных точек (случайном с независимыми координатами либо с использованием ЛП т последовательности) внутри области пробных движений.
Скорость сходимости достаточно высокая: 5-12 итераций. Обычно коэффициент уд берется единичным.
2. Учёт ограничений на искомые переменные
Допустимая область X формируется за счет существующих ограничений типа неравенств
/(х) = шш, ¥} (х) < 0, у = 1, р, (6а)
и типа равенств
/(х) = шш, Фк (х) = 0, к = 1, д . (6б)
Они могут присутствовать по отдельности и вместе.
Ограничения неравенствами относятся к менее жестким ограничениям по сравнению с ограничениями равенствами. При наличии только ограничений неравенствами в достаточно широкой допустимой области X удается сравнительно просто реализовать вышеуказанную процедуру размещения (на каждом рабочем шаге) пробных точек, не выходящих за пределы области X.
В остальных случаях необходимо использовать штрафы на этапе рабочих движений. Пробные точки при этом равномерно размещаются в прямоугольной области П1 с центром в точке X . Большая часть этих точек (или все) находится вне допустимой области. Для этих точек формируются штрафы. Они двух видов.
1. Нормированное ядро сроится в виде произведения ядер для минимизируемой функции, для функций нарушаемых ограничений неравенств (6а) и модулей всех функций ограничений равенств (6б):
- (о _ ^.щщ т _ р („ (і ) )
-г*,тіп п > -г^.шіп ,тт-'
У -(г).
/ ^ г *,тіп
I _1
т(і) _ /00 - /тш „(і) _Фі (Х0)) -Фі
( ^1 ( т, ^
І 2
П -*(о П -*()
V іє (і) У
К.
V к _1 У
(7а)
где я/)тт _ “-----^, яфі _ —-----------------^-, і _ 1 ті ; при і є 3() є 1, т все-
/тах /тіп Фі,тах Ф ],щіп
гда 0 <ф;(х(і)); £« _ Ік(^ и , * _ЇТт, ; 1 <^,1 <р,.
1 к к |тах - 1 к к |тіп
2. Минимизируется не исходная функция /(х), а функция I (х) со штрафными добавками; она в пробных точках имеет вид
г(0 _ 7
.(і)) _ ^тіп
+ тах <! р1
(■
, к _ 1, ш2, (7б)
/ = 1, и, 1 < Р1; 1 < Р2 .
При каждом фиксированном значении I для пробной точки х() в расчетах участвуют только те функции фу (х(г)), которые больше нуля (т. е. те функции, для ко-
торых ограничения-неравенства нарушены): 0 <ф. (х(г)) при і є 3(і) є 1, ш1 .
Для третьей группы алгоритмов применяются обе идеи: формирование пробных точек внутри допустимой области (когда это принципиально возможно) и создание штрафов (на основе модулей функций для ограничений равенств) при совершении рабочих движений.
Для глобальной минимизации функций синтезировано достаточное число алгоритмов: 1) 9 - при наличии только ограничений типа неравенств, 2) 8 - при наличии ограничений типа равенств, 3) 16 - при наличии обоих типов ограничений.
По указанным схемам учёта ограничений типа неравенств и типа равенств создаются группы алгоритмов для: 1) поиска главных минимумов многоэкстремальных функций [7, 8]; 2) многокритериальной оптимизации [9]; 3) решения систем нелинейных уравнений. Вторая проблема кратко затронута в заключении.
15(Х_с5) _ 6| х _ 2|1Л +4| х2 _ 2|17 +7, І6( х _с6) _ 5| х1 _ 4|1Л +5| х2|18 +9,
17(Х_с7)_6|х1 _4|06 +7|х2_4|06 +4 , 18(Х_с8)_6|х1 + 4|06 +6|х2_4|16 +3 ,
19(Х _ с9) _ 3|х1 + 4|12 +3|х2 + 4|05 +7.5, 110(х _ с10) _ 2| х1 _ 3|09 +4| х2 + 5|03 +8.5,
которая имеет 10 минимумов (см. рис. 1) в точках: (0; 0), (-2; 0), (0; -2), (0; 4), (2;
2), (4; 0), (4; 4), (-4; 4), (-4; -4), (3; -5). Глобальный минимум находится в начале координат: I(0) _ 0 .
На рис. 1 показана минимизируемая функция, а на фоне линий равных уровней функции (8) нанесены первые шаги (один или два) алгоритма (5) при четырех различных начальных условиях. Область X прямоугольная и достаточно широкая. Скорость сходимости высокая.
На рис. 2 представлены изменения от шага к шагу значения координат, размеров области поиска и минимизируемой функции. В алгоритме ядро параболическое Р6- (8) _ (1 _ 82 )5 при 5 _ 10 , число пробных точек п _ 50, а также д _ 2 и у _ 1.
3. Пример
Рассмотрим многоэкстремальную функцию
Іг( х _ с,) _ 6| хх|2 +7| х212,
І3(х_с3)_5|х Г3 +5|х2 + 2|13 +5,
12(х_с2)_5|х + 2|05 +5|х2|05 +6 , 14( х _с4) _ 4| х1|08 +3| х2 _ 4|12 +8,
I (х) _ тіп(Іг- (х _ сі), і _ 1,10},
(8)
хг
б
х1
Рис. 1. Функция (8), линии равных уровней её и первые шаги алгоритма
_Дх # _ Дх2
1 (х1, х2)
Рис. 2. Изменение по итерациям координат, размеров области поиска и значений функции
х
Минимизируем ту же 10-экстремальную функцию (8) при наличии двух ограничений неравенств:
ф1 (х) = -х1 - 6 + х2 < 0 (или х2 < х1 + 6 ), ф2 (х) = х1 - 6 - х2 < 0 (или х1 - 6 < х2),
выделяющих полосу в окрестности диагонали х2 _ х1, и одного ограничения равенства
у1(х) = х1 + 4.25Бтх1 _ х2 _ 0 ,
которое устанавливает гармоническую зависимость (х2 _ х1 + 4.258ІПх1) между координатами.
На рис. 3 на фоне линии равных уровней показана допустимая область и линия ограничения равенства.
Поведение минимизируемой функции (8) вдоль линии ограничения равенства представлено на рис. 4. Глобальный минимум её внутри допустимой области приходится на начало координат.
Рис. 3. Линии равных уровней функции (8) и ограничения
/(х1, х1 + 4.25 sin х1)
Рис. 4. Функция (8) при х2 = х1 + 4.25Бт х1
Изменение основных переменных при поиске глобального минимума на основе применения многомерного ядра (7 а) представлено на рис. 5. Использованы следующие ядра и параметры: ядро линейное р, (g) = (1 - g)ж с коэффициентом селективности , =50, а также п =100, д = 2, у = 1, Р1=1, Р2=1.
Качество работы алгоритмов оптимизации достаточно высокое.
/(х1) -О- ф1 О- ф2 ^1
Рис. 5. Изменение по итерациям координат, минимизируемой функции и левых частей ограничений неравенств и равенств; реализован алгоритм с многомерным нормированным ядром (7а)
/среди (-^Ь -^2 ) ^К( - ^| - -'^2 )
Рис. 6. Изменение по итерациям координат, усреднённого значения (по пробным точкам) минимизируемой функции и левой части ограничения равенства; алгоритм с нормированным ядром (7 а); 100 % помеха
Алгоритмы сохраняют работоспособность при высоком уровне аддитивных помех измерения оптимизируемой функции. На рис. 6. приведены результаты поиска минимума функции (8) при 100 %-й равномерно распределённой помехе. В параметрах алгоритма увеличено лишь число пробных точек: п = 200.
Многокритериальная глобальная оптимизация имеет следующую форму записи:
Для попадания в область Парето, образуемую положением глобальных минимумов всех минимизируемых функций, при совершении каждого рабочего шага используем либо произведение ядер
либо переход к одной минимизируемой функции во всех пробных точках:
Здесь а(, / = 1, г , положительные коэффициенты, отражающие приоритеты (лица, принимающего решение) минизируемых функций, и наборы этих параметров дают различные предельные точки внутри области Парето.
Указанные выше способы учета ограничений порождают другие алгоритмы многокритериальной глобальной оптимизации. В настоящее время их 38.
При наличии сравнительно простых ограничений неравенств, обеспечивающих малые вычислительные затраты при получении пробных точек внутри допустимой области, построены также алгоритмы поиска глобальной седловой точки [6, 10]: /(х, у) = тт тах. Учёт более сложных видов ограничений ждёт своего
Для исследования синтезированных алгоритмов и решения на основе них реальных задач глобальной оптимизации разработан комплекс программ [11]. Этот комплекс тоже адаптивно перестраивается по мере расширения спектра алгоритмов, более полно учитывающих наличие ограничений.
Самостоятельная проблема при синтезе алгоритмов возникает, если осуществлять оптимизацию не только по непрерывным переменным (как было рассмотрено выше), а и по смешанным переменным, включающим в себя и дискретные переменные. Они могут быть трёх типов: 1) дискретные неупорядоченные, 2) дискретные упорядоченные и 3) дискретные числовые переменные. Для первого типа переменных не остается ничего другого, как перебирать их возможные значения и распараллеливать процесс вычислений глобальных минимумов всех возникающих функций непрерывных переменных. При втором и третьем типах дискретных переменных возможный путь состоит в переходе от соответствующей дискретной переменной к непрерывной. Учёт однозначной функции соответствия между ними необходим просто для синхронизации процесса вычислений (при этом мы не
Заключение
І} (х) = ШІП, ] = 1, Р, X = хт )Т .
І ,ШІП
п
і ,шах
:,ШІП
І(1) = І (х(і)) = шах
хєХ уєі
решения.
отходим от дискретных значений соответствующих переменных) и для обратного перехода от рассчитанных непрерывных значений к истинным дискретным значениям.
ЛИТЕРАТУРА
1. Медведев А.В, Цыкунова И.М. Об алгоритмах случайного поиска // Применение вычислительных машин в системах управления непрерывными производствами: сб. статей. Фрунзе: Изд-во «Илим», 1975. С. 81-92.
2. Рубан А.И. Метод непараметрической оптимизации стохастических объектов // Системы управления: сб. научных работ. Вып. 1. Томск: Изд-во Том. ун-та, 1975. С. 101-107.
3. Экстремальная радионавигация / В. И. Алексеев и др. М.: Наука, 1978. 280 с.
4. Рубан А.И. Метод непараметрической поисковой глобальной оптимизации // Кибернетика и вуз: Сб. научных работ. Вып. 28. Томск: ТПУ, 1994. С. 107-114.
5. Рубан А. И. Метод непараметрической поисковой оптимизации // Изв. вузов. Физика. 1995. Т. 38. № 9. С. 65-73.
6. Рубан А.И. Глобальная оптимизация методом усреднения координат: монография. Красноярск: ИПЦ КГТУ. 2004. 303 с.
7. Кузнецов А.В., Рубан А.И. Поиск главных минимумов многоэкстремальных функций при активном учёте ограничений неравенств // Техника и технологии: журн. СФУ. 2010. Т. 3. № 3. С. 335-346.
8. Кузнецов А.В., Рубан А.И. Алгоритмы метода усреднения координат при поиске главных минимумов многоэкстремальных функций // Вестник СибГАУ. 2010. Вып. 5(31). С. 36-41.
9. Рубан А.И. Глобальная многокритериальная минимизация функций методом усреднения координат // Информационные технологии в территориальном управлении, промышленности, образовании: сб. статей. Томск: ТУСУР, 2005. С. 142-156.
10. Rouban A.I. Method for global minimax optimization in continuous space // Advances in Modelling & Analysis: Series A. Mathematical, general mathematical modelling. France: A.M.S.E. 2007. V. 44. No. 2. P. 46-65.
11. Кузнецов А.В., Рубан А.И. Комплекс программ «Global Optimizer v2.0» // Свидетельство о государственной регистрации программ для ЭВМ № 2011617970 «Global Optimizer v2.0». Зарегистрировано в Реестре программ для ЭВМ 12 октября 2011 г.
Рубан Анатолий Иванович Сибирский федеральный университет
E-mail: rouban@mail.ru Поступила в редакцию 30 мая 2012 г.
Rouban Anatoliy I. (Siberian Federal University). Global optimization method based on the selective averaging coordinate with restrictions.
Keywords: global optimization, selective averaging of coordinates, restrictions of the form of inequalities and equalities.
There are described ideas of design of non-differentiable global optimization algorithms, which are based on: separation in time of exploratory and pattern steps, selective averaging of coordinates on the results of test movements, adaptive reconstruction the size of rectangular region of test motions and taking into account the restrictions in the form of inequalities and equalities.
Inequality restrictions are less restrictive than equality constraints. If there is only inequality restrictions and a fairly wide feasible region one can (before every working step) relatively simple implement the procedure of placing the sampling points in the admissible region.
In other cases, penalties are used. Sampling points with are uniformly placed in a rectangular area centered at the point from which the algorithm performs the pattern step. Most of the sampling points (or all) are out of the admissible area. For these points are formed penalties. They are of two types: 1) the calculation of the normalized core pattern steps built in the form of the product cores for function to be minimized, for functions with violated inequalities and for modules of
all functions with equality restrictions, and 2) minimizing the penalty function. In test points the penalty function has several forms built on combinations of operations of maximization and summation.
In all global optimization algorithms the transformations of optimized functions and functions of restrictions are performing for dimensionless variables. This increases accuracy and reduces the number of adjustable parameters in the algorithms.
Convergence rate of the algorithm is rather high: 5-12 pattern steps in the absence and in the presence of additive noise of high intensity for optimized functions.