CC BY
25
УДК 519.1
ГРАВИТАЦИОННЫЕ АНАЛОГИИ В ЗАДАЧЕ ОПТИМИЗАЦИИ
Н.Ф. Жихалкина, Р.Т. Файзуллин
The most spreaded optimization methods are oriented on the looking for a local minimum and results depend from an initial approximation. Therefore it is frequently nessesary to use additional methods but usually the payment for it is the complicating a structure of the algorithm and reducing stability of its work.
The new method of looking for an extremum is suggested in this paper. This algorithm we have named gravitational. It is based on mechanics namely gravitational laws and genetics analogies.
Besides it is represented some theoretical and numerical results.
1. Введение
Нелинейная задача безусловной оптимизации с действительными переменными математически формулируется следующим образом:
Дано: F : Rk —>■ R
Найти: г * £ Rk : F(r*) > F(r), Vr £ Rk, где Rk - Ачмерное евклидово пространство.
Существуют хорошо известные алгоритмы решения поставленной проблемы в предположении гладкости целевой функции: метод Ньютона, методы секущих, в частности формула Пауэлла, и т.п [1]. Данные методы не гарантируют глобальность найденного максимума. Здесь важную роль играет начальное приближение, с которого стартует оптимизационный алгоритм. Для многомерных пространств близость начального приближения к оптимуму представляет собой существенное ограничение.
Наряду с традиционными градиентными методами при решении оптимизационных задач применяются эволюционные и генетические алгоритмы. Данный подход наиболее характерен для дискретной оптимизации.
Следующий класс алгоритмов поиска экстремумов функции включает в себя методы, связанные с квантово-механическими, в частности, гравитационными аналогиями.
0 1998 Н.Ф. Жихалкина, Р.Т. Файзуллин
E-mail: zhihal@univer.omsk.su, faizulin@univer.omsk.su Омский государственный университет
2. Классический градиентный алгоритм
Традиционно для поиска экстремума функции используется гибридный кеази-нъютоновский алгоритм. При этом помимо целевой функции F задается начальное приближение х0 [1]. На каждой итерации строится локальная модель, описывающая поведение целевой функции F вблизи текущего приближения ж&. Если найденное решение (или точка близкая к решению) модельной задачи неудовлетворительно, то для поиска Xk+\ применяется глобальная стратегия.
Применительно к рассматриваемым в данном разделе алгоритмам термин «глобальный» используется для обозначения метода, который сходится к локальному экстремуму, но почти из любой начальной точки. Тогда как для локально сходящихся методов необходима достаточная близость начального приближения к точке оптимума. Задача поиска глобального экстремума является более сложной и при ее решении наряду с классическими используются эвристические, эволюционные и т.п. методы.
Рассмотрим одну из реализаций предложенного алгоритма. Пусть задана целевая функция F : R -У R. На каждом шаге итерационного процесса используется метод Ньютона: в некоторой окрестности максимума целевой функции
F\х*) = 0 —> F(х*) = maxF(x) .
Пусть F'(x) = G(x).
G'(x+) ss
Ay
Ax
G(x) - G(x+)
При G(x+) = 0 получаем: x+ = x — ^т^+у- Таким образом, шаг метода Ньютона осуществляется по следующей схеме:
Щ+1
= хг
F'(xn) F"{xn) '
В качестве глобальной стратегии можно взять метод деления пополам:
while \F(xn+1)\ < \F(xn)\ do
r, __xn+l-\~xn
^n+l •— 2
Если Ньютоновская точка xn+i не приводит к увеличению |Е(ж)|, то разумная стратегия заключается в дроблении шага в обратном направлении от xn+i к хП} пока не встретится точка х+ : |Е(ж+)| > \F(xn)\.
В методах глобальной оптимизации возникает две проблемы:
1) направление шага;
2) длина шага.
При дополнительных ограничениях на целевую функцию поставленные задачи успешно решаются [1], но здесь важную роль играют размерность пространства и гладкость функции F.
3. Генетические алгоритмы
3.1. Основные операторы ГА
При работе с генетическими алгоритмами используется следующая терминология: хi (в простейшем случае битовая строка длины п) - хромосома; каждая позиция в строке - локус; значение переменной в локусе - аллель; сама переменная - ген. Множество хромосом называется генотипом, который определяет фенотип или конкретного индивида.
Эволюционные алгоритмы, моделирующие естественные процессы, применялись для решения оптимизационных задач еще в 1960-е годы. Основная работа была направлена на изучение воздействия мутации, скрещивания и селекции - основных операторов этих алгоритмов, на эволюцию генотипов в случае нелинейной целевой функции или функции пригодности. Наряду с теоретическим анализом активно использовались и продолжают использоваться результаты, полученные в ходе численных экспериментов.
Еще в 1966 г. Bremermann ввел термин «скрещивание» (mating) при рекомбинации двух (и более) родительских строк для создания строки-потомка.
В 70-е годы независимо друг от друга появились генетический алгоритм [Holland, 1975] и эволюционная стратегия [Rechenberg, 1973 и Schwefel, 1981].
Теоретический анализ этих алгоритмов базируется на популяционной генетике и статистике и при некоторых предположениях позволяет определить:
- приемлемую схему селекции;
- прогресс популяции (улучшение среднего значения целевой функции) при данной схеме селекции и рекомбинации;
- прогресс популяции при данной схеме селекции и мутации;
- соотношение селекции, мутации и рекомбинации в конкретном алгоритме.
Как упоминалось выше, селекция, мутация и рекомбинация являются основными операторами генетических алгоритмов. Два последних (мутацию и рекомбинацию или скрещивание) следует оценивать в соответствии с вероятностью, с которой они находят лучшее решение. Мутация базируется на случайности, шаг мутации почти непредсказуем: на данном шаге мутации с вероятностью m изменяется каждый бит строки. Рекомбинация (или crossover) - это некоторый глобальный поиск, который в принципе ограничен заданной популяцией. Обычно crossover оперирует двумя строками. Существует множество вариантов работы этого оператора, в частности:
1) однородная рекомбинация:
ж = (жь ..., ж„), у = (yi,..., Уп) -У Z = (щ, . . . , Zn),
где Zi = жi либо Zi = уг-. Вероятность выбора жг- или уг- равна 0.5;
2) в результате воздействия оператора на две строки вновь получаем два потомка по следующей схеме:
(Ж1, • • • 1 Жга)? (У1, • • • ? Уп) t (Ж1, . . • , Ж^, Ук+1, • • • 1 Уп) (У17 • • • 7 Vki %k+17 • • • 7 %п))
где к = 1,. . . , п.
Остается вопрос: как совмещать эти стратегии (селекцию, мутацию и рекомбинацию)? Большинство результатов, полученных в этом направлении, являются эмпирическими и основываются на численных экспериментах и компьютерном моделировании.
3.2. Обзор эволюционных алгоритмов
В общем случае генетический алгоритм является эвристическим оптимизационным алгоритмом, который на каждой итерации использует популяцию точек вместо одного начального приближения, выбираемого в классических методах оптимизации. В процессе работы вычисляются лишь значения целевой функции и не требуются данные о производной. Пространство оптимизации дискретно и состоит из конечного числа точек. Случайный характер поиска ведет к тому, что вероятность оптимальности найденного решения на итерации t не равна единице. Но при ?, стремящемся к бесконечности, вероятность того, что оптимальное решение не попадет ни в одно из t поколений, равная (1 — (m/2)n)tM, стремится к нулю. Здесь п - длина хромосомы, (га/2)п - оценка вероятности перехода фиксированного представителя популяции Pit) в точку оптимума (вероятность одновременной мутации во всех позициях генотипа с удачными исходами), М - размер популяции.
1. (р + X) - эволюционная стратегия
Шаг 1: создание начальной популяции размера Л.
Шаг 2: вычисление целевой функции или пригодности С(жг), г = 1,. . . , А.
Шаг 3: выбор р < А лучших индивидов.
Шаг 4: создание - потомков из каждого из р индивидов с малой вариацией координат. Возврат на шаг 2.
По сути, данный алгоритм представляет собой случайный поиск, использующий селекцию и мутацию.
2. Генетический алгоритм
Шаг 0: представление исходной задачи в терминах ГА.
Шаг 1: создание начальной популяции Р(0) = ж°,. . . , ж^, t = 0.
Шаг 2: вычисление средней пригодности текущей популяции:
N
и «нормализованной» пригодности каждого индивида:
Шаг 3: вычисление для каждого индивида вероятности, с которой он может быть выбран в качестве родителя для создания потомка, пропорционально его
пригодности (пропорциональная селекция):
f(xi)
N
Используя это распределение, выбирается N векторов из популяции P(t), тем самым создается набор родителей.
Шаг 4: случайным образом формируется N/2 пар. Применяя к каждой паре, с определенной вероятностью,
1) оператор скрещивания (crossover),
2) оператор мутации и т.п.,
формируется новая популяция P(t + 1). Возврат на шаг 2.
3. Параллельный генетический алгоритм
Шаг 0: представление исходной задачи в терминах ГА.
Шаг 1: создание начальной популяции.
Шаг 2: каждый индивидуум выбирает партнера для скрещивания среди
Шаг 3: используя операторы рекомбинации и мутации, генерируется потомок.
Шаг 4: в следующей популяции потомок заменяет родителя, если он лучше с точки зрения критерия оптимизации. Предварительно может производиться
некоторый локальный поиск в окрестности потомка с целью улучшить значение целевой функции. Возврат на шаг 2.
В этом алгоритме информационный обмен напоминает диффузию.
4- Селекционный генетический алгоритм (Breeder GA)
Шаг 0: представление исходной задачи в терминах ГА.
Шаг 1: создание начальной популяции -Р(О) размера N, t = 0.
Шаг 2: каждый индивидуум может выполнить локальный поиск.
Шаг 3: выбирается Т% популяции для скрещивания, тем самым создается набор родителей.
Шаг 4: случайным образом формируется N/2 пар. Применяя операторы рекомбинации и мутации, переходим к новой популяции P(t-\-1), t —У t+1. Возврат на шаг 2.
Главное отличие этого алгоритма от классического ГА - в методах селекции. Данный алгоритм может использовать различные селекционные стратегии в зависимости от решаемой задачи.
3.3. Теоретические и эмпирические результаты
В теории, связанной с эволюционными алгоритмами, наиболее значимой является теорема о схемах или фундаментальная теорема генетических алгоритмов [Schema theorem, Goldberg, 1989], которая позволяет предсказать число строк
своих ближайших соседей.
со схемой Н в поколении t + 1 по значению этой характеристики для поколения t:
m(H, t + 1) = m(H, t) ,
где f(H,t) - средняя пригодность схемы Н в поколении t. К сожалению, этот результат справедлив только для пропорциональной селекции. В общем случае данная теорема неприменима.
Неудачи в теоретических расчетах привели к активному развитию экспериментальных исследований. Так, эмпирически было найдено оптимальное соотношение основных параметров генетического алгоритма [Grefenstette, 1986; Schaffer, 1989]:
ln(N) + 0.93/n(m) + 0.45/n(n) = 0.56,
где N - размер популяции, m - темп мутации, п - длина хромосомы. Для расчетов удобна приближенная формула:
= 1.7 .
Рассмотрим более подробно искусственную селекцию. На каждой итерации ГА в качестве родителей выбирается 10 — 15% лучших, с точки зрения пригодности, особей. Их число определяет так называемый порог, равный trunc%. Пусть R(t) = fit + 1) — f(t) - результат селекции на шаге t и пусть S(t) = /s(t + l) — f(t) - селекционная разность, где fs(t-\-1) - средняя полезность родителей, отобранных для селекции.
Задача: зная S(t), найти R{t).
Будем считать [Falconer, 1981], что R(t) = b{t)S(t)} где b{t) - реализованная наследственность. Эта величина либо переносится из предыдущей популяции, либо оценивается на текущей итерации различными методами. Но обычно b{t) берется за константу для некоторого числа поколений. Интенсивность селекции / = ^у, где cr{t) - отклонение пригодности индивидов, может быть вычислена аналитически, если пригодность имеет нормальное распределение [Bulmer, 1980]. Численные эксперименты показывают, что данное уравнение применимо для многих практических приложений и в том случае, когда распределение не является нормальным.
Для произвольного распределения получена следующая оценка [Nagaraja, 1982]:
м
<r{t)
Таким образом, R(t) = blcr{t). Очевидно, что для увеличения R(t) требуется увеличивать b или cr{t).
Заканчивая обзор генетических алгоритмов, необходимо отметить, что имеются различные оценки R(t)} например, для бинарной одноэкстремальной функции с однородной рекомбинацией, для некоторых дискретных функций. Кроме того, установлено, что оптимальный темп мутации m для дискретных
<
(100 — trunc)
trunc
функций равен 1 /га, причем как для больших, так и для малочисленных популяций [2]. Если вернуться к вопросу об ожидаемом числе строк с заданной схемой Н в поколении t, то представляет интерес оценка числа различных схем в популяции размера п на строках длины /, а также скорость обработки этих схем. Некоторые результаты для бинарных строк постоянной длины (в предположении пропорциональной селекции) могут быть найдены в статье [3]. Нетрудно показать, что ожидаемое число схем в популяции
I
S(n3) = 3'-^c'2-(l-(l/2)T .
%=О
Известно, что независимо от п, при t -У оо, S(t) -У 21 [3]. Знание этих величин дает возможность вычислить темп обработки схем М, в частности, если все схемы имеют равный приоритет, обрабатываются бинарные структуры постоянной длины, то
М(га, /, (3) =
S(n,l)-2l
ПСП1 Р ’
где пс - число поколений; (3 - степень параллелизма процессора. Если целевая функция задана на алфавите из к элементов, начальная популяция выбирается случайно, выбор г-го элемента алфавита производится с вероятностью 1/к, с - требуемое число дубликатов одной схемы, тогда ожидаемое число схем в популяции
I п
5'(n, I, с, к) =(k+ 1)' - £ С‘к- Y, С?(1 - {1/к)У(1/кр~» .
г=0 j—п — с-\-1
Последняя формула позволяет расширить полученные результаты и на более общий класс функций.
4. Переход к гравитационному алгоритму
4.1. Постановка задачи и решение проблемы N тел
Рассмотрим несколько иной подход к поставленной проблеме, который использует физические аналогии и вероятностную схему преобразования координат взаимодействующих частиц.
Пусть в начальный момент времени t = 0 в пространстве оптимизации, представляющем собой связную замкнутую область, случайным образом генерируется N точек, взаимодействующих с силой, обратно пропорциональной квадрату расстояния (в данном случае это гравитационные силы). В качестве области распределения берется Ачмерный куб, который затем замыкается в тор, чтобы, в частности, обеспечить постоянное число точек внутри области. Масса г-й частицы зависит от ее координат и равна значению целевой функции в этой точке. Таким образом, функция F аппроксимируется некоторой равномерной сеткой, в узлах которой расположены материальные точки с массами тг = F{fi).
На каждом шаге алгоритма решается задача N тел [4], иными словами, происходит взаимодействие материальных точек под действием гравитационных сил на основе закона тяготения.
Рассматривается модель гравитационного взаимодействия материальных точек без столкновений. Согласно законам движения Ньютона и при некоторых дополнительных предположениях [5], формулы расчета координат взаимодействующих частиц принимают вид:
3? = £f=
,(г7-гТ)
j=l,j^i If*™—г™ rT+i = yn + C(J(t)^n ^
где C > 0 - константа; a(t) - некоторая функция, зависящая от шага по времени; к - размерность пространства (в первом приближении к = 2).
Одна из трудностей изучения системы с дальнодействующими силами состоит в том, что каждая материальная точка может эффективно взаимодействовать с любой другой. Поэтому в системе из N частиц существует N^N2 ^ взаимодействий. Следовательно, число N ограничено сверху величиной порядка 103, так как при этом приходится следить за миллионом взаимодействий. Однако 103 частиц достаточно для проявления статистических свойств ансамбля, и это дает возможность «экспериментального» изучения поведения системы.
Кроме того, всю систему целесообразно разделить по отношению к текущей точке fi на два объема VI и V2, где VI - шар радиуса R с центром в точке у, a V2 - оставшаяся часть системы. Если радиус R достаточно велик, то частицы внешнего объема по отношению к выделенной точке ведут себя как непрерывная среда, а частицы, находящиеся внутри шара VI, проявляются как дискретная структура, поэтому их нужно описывать отдельно. Величина R называется дебаевским радиусом. Гравитационное поле в точке у можно представить в виде суммы двух слагаемых: первое слагаемое Е1 - локальное поле, создаваемое частицами объема VI, второе слагаемое Е2 - среднее поле, создаваемое кажущимся континуумом зарядов внешнего объема V2, таким образом, Е = El + Е2 [4]. Зависимость суммарной гравитационной силы, действующей на частицу со стороны соседей вне шара VI, от расстояния между взаимодействующими частицами близка к экспоненциальной (в непрерывном случае порядка 1 /Rk~1). Следовательно, выбирая дебаевский радиус определенным образом, слагаемым Е2 можно пренебречь.
В рассматриваемой системе для каждой частицы вводится дебаевский радиус, точнее, более технический элемент - радиус взаимодействия R. Если полный объем, занятый системой, равен V, то с каждой частицей связан свободный объем:
1 = — /?3 = — п 3 N ’
здесь п - плотность частиц [к = 3). Характерное расстояние взаимодействия:
R = ез 1пш^.
В случае произвольной размерности пространства радиус взаимодействия R вычисляется по новой формуле. Во так как важна лишь зависимость R от к,
то, полагая V = 1 и оставляя объем Агмерного шара равным |7Г, получаем:
R = ДгпчД
И наконец, шаг по времени t следует выбирать достаточно малым, чтобы для всех i выполнялось условие:
R
Итак, под действием гравитационных сил частицы стягиваются к центру масс системы, т.е. к тем подобластям, где значение функции максимально. В зависимости от целевой функции и от расположения оптимума внутри рассматриваемой области, центр масс может не совпадать с областью концентрации наиболее тяжелых частиц. Для решения этой проблемы на каждом шаге алгоритма «лучшие», с точки зрения критерия оптимальности, частицы фиксируются и наделяются фиктивной массой, в несколько раз превышающей значение целевой функции.
4.2. Необходимость малых возмущений системы
Использование только гравитационных аналогий делает данный алгоритм незащищенным от зацикливания, в частности, от скопления частиц в окрестности седловой точки целевой функции. Поэтому на каждой итерации осуществляется так называемое «отражение» в одной из координатных плоскостей ф, где г-я координата частицы вычисляется через j-e координаты «соседей». Выбор плоскости «отражения» носит вероятностный характер.
По аналогии с гравитацией формулы для вычисления координат г-й частицы имеют вид:
71+1
71+1
- xi + Ci(Ti{t)
= - Cpr^t) Ef=ij#
| jfn _jfn |2 | jfn _jfn |2
5
остальные координаты на данном шаге не меняются.
Очевидно, что расстояние ггд между частицами г и j в плоскости «отражения» увеличится:
[(гдП2 = (Д - +Д + (у? - У])2 [+»r+1]2 = ДД1 - Д+1)2 + (уД1 - уД1)2 =
д
= [(гдТ\2 + + •
Как показывают численные эксперименты, данный процесс, который интерпретируется как малые возмущения системы, позволяет не только успешно решать проблемы, связанные с зацикливанием, но и компенсировать ошибки выбора шага по времени.
4.3. Сходимость метода
При работе гравитационного алгоритма генерируется последовательность {г*}, где 1 < * < и t = 1,2,.... Все элементы последовательности принадлежат некоторому тору Т. Пусть t = Д, среди N точек ту1 выбирается г*1, такая, что F(r)1) — ^(д1)- Зафиксировав эту частицу, осуществляется переход к t = Д. Вновь отбирается лучшая, с точки зрения критерия оптимизации, точка и фиксация переносится на нее. В силу неподвижности г*1, справедливо неравенство F(r^) > F(r£). Рассмотрим последовательность {г))}
Определение 1. Метод сходится, если HniF^r*3) = maxF(r), при j стремящемся к бесконечности и г принадлежащем Т.
Замечание 1. Будем считать, что в некоторой окрестности предельной точки последовательности {г Д } функция F имеет один максимум.
Лемма 1. Flycmb функция F непрерывна. Тогда данный метод является сходящимся.
Доказательство. Так как все точки последовательности принадлежат некоторой связной замкнутой области (тору), то из последовательности {гг( можно выбрать сходящуюся подпоследовательность, такую, что ту* —>■ г* £ Т.
Из непрерывности функции F следует : limF(r))) = F(r*). Последовательность {F(ryp}j=iv.. монотонно возрастает (по построению) и ограничена сверху, в силу существования максимума (функция, непрерывная на компакте, достигает своей верхней границы). Таким образом, эта последовательность также сходится, причем limF(r^) = F(r*).
Осталось показать, что F(r*) = maxF(r), г* £ Т. Допустим, что это не так, тогда F(r*) < maxF(r), г* £ Т.
В силу определения г* и непрерывности F существует г, такая, что F(r) > F(r*) (*) и F возрастает в направлении г*г. Пусть F(r) — F(r*) = г.
На шаге I фиксируется частица г)), такая, что F(r)() > F(r)(), j = 1,. . . , N. Если выделенная точка удовлетворяет условию (*), то все доказано. Пусть C(r‘;)<F(r>).
Рассмотрим сферу с центром в точке г*, радиус которой равен г. Под действием оператора «отражения» расстояние между частицами увеличивается, и в плоскости действия алгоритма они движутся по траекториям, близким к раскручивающейся спирали. Таким образом, на следующем шаге гравитационного алгоритма движение к центру масс системы будет происходить под другим углом. Так как точка г* является предельной точкой последовательности {г*3} и в силу независимости процесса «отражения» от гравитационного взаимодействия, вероятность того, что найдется элемент последовательности, расположенный сколь угодно близко к отрезку г*г, стремится к единице при j -У оо. Следовательно, учитывая непрерывность функции F, можно утверждать, что найдется к > /, т.ч. F(r*) < F(r\kk). Данное противоречие и доказывает лемму.■
Замечание 2. Лемма 1 не дает оценки скорости сходимости данного метода. В режиме реального времени мы можем не достичь точки точного экстремума.
4.4. Общие характеристики алгоритма
Подобно генетическим алгоритмам, гравитационный алгоритм является методом нулевого порядка: поиск оптимума построен только на вычислении значений критерия оптимальности и не требует дополнительных знаний о производных и константе Липшица. На функцию F не накладываются ограничения дифференцируемости, одноэкстремальности.
4.5. Связь гравитационного и генетического алгоритмов
Итак, пусть в пространстве оптимизации задана целевая функция Г(ж), х £ Rk и решается задача на максимум. Проведем аналогии между гравитационным и генетическими алгоритмами.
1) В предположении непрерывности функции F, в начальный момент времени t = 0 в пространстве оптимизации случайным образом (равномерно) генерируется N материальных точек. Тем самым получаем начальную популяцию Р(0) = {ж°,. . ., х%\х° £ Rk}. Каждой частице i ставится в соответствие масса, равная F{ж°).
2) Вычисление Fit), F{ж*), р(ж(Л) происходит по стандартным формулам ГА.
3) На этом этапе возможно реализовать два подхода. Первый: если придерживаться классического генетического алгоритма, то
а) в соответствии с вероятностью р(ж(Л) выбирается N точек из текущей популяции, либо, следуя селекционному генетическому алгоритму, рассматривается Т% самых «тяжелых» частиц; далее из этих векторов формируется N/m наборов по m элементов;
б) теперь к каждому набору необходимо применить операторы рекомбинации и мутации; возможно выделение так называемой элиты, т.е. тех материальных точек, на которых значение целевой функции максимально в данной популяции Pit), и перенос этих частиц в следующее поколение без изменений.
Второй подход заключается в следующем: предлагается перейти от стандартных генетических операторов рекомбинации и мутации к одному шагу гравитационного алгоритма. В этом случае целесообразно рассмотреть другой вариант ГА: параллельный генетический алгоритм, в котором каждый индивидуум использует для создания своего потомка (ж) —у x\+l) ближайших соседей. По-прежнему выделяется элита, которая не меняется при переходе в следующую популяцию, но участвует в создании потомков соседей.
Идея использования более двух родителей для рекомбинации неоднократно применяется в классических ГА. Muhlenbein [1989] использовал 8 родителей, мультирекомбинацию применяли Rane и Ruttkay [1994] и т.д. Суть такой ре-
комбинации заключается в том, что два родительских аллеля случайно выбираются из целого множества родителей.
Итак, шаг гравитационного алгоритма формально состоит из двух частей:
1) решение задачи N тел (задачи гравитационного взаимодействия частиц) для всей системы в целом или для каждой точки и подмножества ее ближайших «соседей»;
2) локальное возмущение системы («отражение» в одной из координатных плоскостей).
В сравнении с операторами ГА можно провести следующие параллели:
- рекомбинация УД «отражение»;
- мутация УД гравитационное взаимодействие.
При «отражении» в качестве родительских аллелей выбираются координаты ближайших соседей, и все они участвуют в расчете новой координаты текущей частицы. Кроме того, гравитационное взаимодействие не является случайным поиском, как в случае оператора мутации, и связано с критерием оптимизации - с производной целевой функции.
Главное отличие заключается в том, что если в генетических алгоритмах мутация играет роль малых возмущений, тогда как рекомбинация выступает основным оператором получения новых особей в популяции, то в гравитационном алгоритме - наоборот: решение задачи N тел является инструментом получения новых точек, близких к экстремуму, а «отражение» вносит в систему элемент случайности.
Представленная связь алгоритмов позволяет обобщить некоторые результаты из теории ГА на гравитационный метод оптимизации. Так, экспериментально полученное оптимальное соотношение между основными параметрами ГА (обратная зависимость темпа мутации от длины хромосомы при постоянном размере популяции) подтверждается в ходе численных экспериментов и для гравитационного метода, если в качестве темпа мутации рассматривать величину шага по времени при решении задачи N тел.
4.6. Градиент целевой функции в гравитационных алгоритмах
Как упоминалось выше, одной из основных проблем в классических методах глобальной оптимизации является поиск направления наискорейшего роста целевой функции, т.е. Vf(r). Перейдем теперь к гравитационному алгоритму и попытаемся проследить его связь с вышеназванными методами.
Отличительной особенностью гравитационных алгоритмов является то, что на каждом шаге итерационного процесса:
а) вместо одного текущего приближения гп рассматривается набор точек тф,. . ., Гдг, в которых известны значения целевой функции;
б) вычисление градиента заменяется решением задачи N тел.
Таким образом, делая шаг из каждой точки на данной итерации алгоритма, используется информация о поведении функции в некоторой окрестности этой
точки. Сравним насколько хорошо при гравитационном взаимодействии аппроксимируется градиент целевой функции.
Пусть F : Rk —> R и число взаимодействующих материальных точек равно N Т 1. Новые координаты текущей точки г^+1 вычисляются по схеме:
N
У1 = е0 + Ca(t)
рднд - ?,
=5)
%=1
I iytfI ГрУI I k
где t - шаг по времени, С = |G, cr(t) = t2 механики. Оценим
для стандартных формул из
N
Е
%=1
F(P;)(P;-y_VF{rZn
I z?n
Г: — Г,
I к
Пусть F(rY) = Mi7 rf = гд \rf — r^ \ = Ri, H = Мгr°^, тогда проекция H
на ось Ож вычисляется по формуле:
=
N
Е
г = 1
Mi(xi - х0)
R§
Задача сводится к оценке \НХ — F'x\ в точке г0. Для произвольной функции G(x):
d2G
dx2
dG
dx
Xq
G(xi) — G(xq)
x0 Xl ~ X°
G(x1) — 2G(x0) + G(x2) h2
где h = x2 — x0 = x0 — X\. Очевидно, что порядок аппроксимации таких схем равен 1, кроме того, обе схемы являются устойчивыми, следовательно, сходятся к точным значениям производных.
Пусть
1) N = 2
2) R\ = R2
3) и м0.
Тогда
_Mi{xi-xq) М2(х2 — ж0) _ х~ Rk + “
М\Х\ — (Mi Т- \12 !./'ц Т- М2х2
Rk
Mi ж Щ-2
Mqxq _i_ М?х?
Rk-2
X]__9М0Ж I
-2 ^ Т2к—2 I
Rk-2
R2
д2
:М
дх2 V Rk~2
= *
Из этой формулы, в частности, следует, что, какова бы ни была размерность пространства к, степень модуля расстояния между взаимодействующими частицами целесообразно брать равной двум. Тот же результат получен в ходе численных экспериментов.
(*) = Т(М + *юи, (2АД + .тЛОЦ. « •
Вообще говоря, мы получаем производные по R} но, предполагая линейную зависимость х от R} нетрудно показать, что в этом случае формулы будут отличаться лишь на некоторую константу. Таким образом,
Нх ~ 2 М'х .
Отсюда
К(го) ~ •
Приведенные расчеты справедливы, если слагаемым, содержащим вторую производную, можно пренебречь. В этом случае гравитационные формулы представляют собой аппроксимацию градиента, индуцированного ближайшими «соседями». Но в окрестности экстремума производная целевой функции М'х близка к нулю, следовательно, именно слагаемое хМхх делает основной вклад в результирующую сумму Нх. Этим объясняется относительно низкая эффективность алгоритма при локальном поиске (в сравнении с традиционным методом Ньютона) и высокая скорость сходимости при использовании алгоритма в качестве глобальной стратегии.
Обращаясь к физике данного процесса, заметим, что суммарная сила, действующая на частицу г*0 со стороны «соседей», не зависит от ее массы. Перенося начало локальной системы координат в данную точку, получаем, что предположение (3) оказывается излишним, так как
— Tq) Т М2^Х2 — Ж0) = А41\Х\ AI2X2 — М1\Х\ — 2A1qXq Т AI2X2 •
Интересно, что в этом случае при R = (ж0 — Xj) формула для расчета Нх представляет собой центральную разность:
\ /1 ./• | — М2х2 К1
Ml - м>
Аж
2М'х + АжМ”х .
Перейдем к общему случаю, когда N - произвольно. Пусть в качестве области оптимизации рассматривается тор. Для любой точки г0 найдутся «соседи», проекции которых на координатную ось будут как меньше, так и больше ж0. Тогда
Г\
Мг
Е
ier
- Vi
cardl ’ Еiei~ ^
cardl
где 1~ - множество индексов г, такое, что жг- < х0. Аналогично вычисляются г2 и М2.
Компактификация пространства (замена Ачмерного куба на тор) может привести к нарушению условия непрерывности целевой функции. Но если значения функции F на границе совпадают (близки к нулю), а этого можно добиться, используя так называемые «сглаживающие» функции, то условие непрерывности сохраняется. Примером такой вспомогательной функции для случая, когда F : R -У R задана на отрезке [ад, х2\, является функция
{fl(x) = у(1 + sin^(X — Xi — !)), Xi < X < Xi + £
(3, .r I • • • X ■
Д(ж) = f (1 + sin^(—x + x2 — !)), x2 — e < x < x2 .
Другой вариант обобщения полученных результатов - это представление системы гравитирующих частиц в виде связного графа. Для каждой материальной точки i выделяется подмножество ближайших соседей Д,. . .,Д, что в графовом представлении задачи соответствует существованию ребер ijm} где m = 1,. . . , I. Задав направление и учитывая, что в качестве области оптимизации рассматривается тор, получаем замкнутую сеть.
Оценка значения целевой функции в узле г, аналогичная допущению (3),
Мп + • • • + Мп
I
~ АД (^ иД
накладывает на функцию F довольно жесткие ограничения и в общем случае не всегда корректна.
Иной подход к решению данной проблемы возможен в том случае, когда целевая функция F гармоническая, т.е. описывает некоторый физический процесс, заданный однородным уравнением Лапласа с условием периодичности на границе:
A F F
О
F
г2
Пусть даны функции АД i = l,...,s, такие, что Д равна единице в узле с номером i и нулю - во всех остальных узлах, на ребре ij функция линейна. Таким образом, если все ребра имеют длину R, то V- = const = —А. Кроме того, имеется некоторая интерполяционная формула Fh = ДЖ'-
Домножим AF на функцию V) и проинтегрируем по ребру щ, если поток направлен из узла i в узел j (иначе - по ребру Д). Интегрируя по частям, получаем:
AFVj = VF
6 узле j
. + —(FI
1 Ry 1
6 узле j
— FI
16 узле г>
= * * *
подставляя вместо F интерполянт и учитывая, что в итоговой сумме первые слагаемые дают нуль (так как задан поток и на границе выполнено условие периодичности, т.е. возможен переход к замкнутой системе), имеем:
О * *) ~ - Fi) ■
Повторяя описанную процедуру для всех Vj j = l,...,s и по всем ребрам, получаем матричное уравнение:
где / Т = (Fi,. . . , Fs), элемент агд матрицы А вычисляется по формуле:
Следовательно, при данных предположениях относительно целевой функции F, допустимость оценки (**) определяется погрешностью интерполяционной формулы.
Остановимся на вопросе аппроксимации функции и ее производной применительно к гравитационному алгоритму. В частности, если имеется некоторая интерполяционная формула Fh(r) = ^СгК(г), то представляют интерес погрешность вычисления производной:
Итак, в рассмотренном примере гравитационные формулы представляют собой интерполяцию производной при помощи фундаментальных решений оператора Лапласа [6].
4.7. Заключение
Вместе с традиционными оптимизационными методами гравитационный алго-
риментах физики высоких энергий при поиске глобального экстремума функции, зависящей от одной переменной. Предложенный метод прошел тестирование на специальных тестовых функциях, в ходе численных экспериментов дано обоснование полученных результатов [5].
Af = 0 ,
dF
дх
Рассмотрим
N
N
Fh(r) = ^Р{гг)Щг) = ^МгЩг) ,
где г £ R3. Пусть Щ = | у — г \,
Rt у7(хг - ж)2 + (уг - уУ + (щ - z)2
dVi(r) х i — х
дх Rf
ритм был применен к задаче геометрической реконструкции событий в экспе-
Представляется возможным использование аналога гравитационного алгоритма при решении задачи гидравлического расчета оптимального режима работы разветвленного нефтепровода.
Литература
1. Дэннис Дж., Шнабель Р. Численные методы безусловной оптимизации и решения нелинейных уравнений. - М.: Мир, 1988.
2. Miihlenbein Heinz Genetic algorithms // local Search in Optimization / Edited by E. AArts and J. K. Lenstra. 1997.
3. Goldberg David E. Sizing Populations for Serial and Parallel Genetic Algorithms II The Proceedings of the Third International Conference on Genetic Algorithms. -California, 1989.
4. Поттер П. Вычислительные методы в физике. - М.: Мир, 1975.
5. Жихалкина Н. Ф., Кисель И. В., Назаренко М. А., Файзуллин Р. Т. Гравитационный метод безусловной глобальной оптимизации // Сообщения ОИЯИ Р5-97-255. Дубна, 1997.
6. Тихонов А. Н., Самарский А. А. Уравнения математической физики. - М.: Наука, 1977.
