УДК 681.518.001.33
ИСПОЛЬЗОВАНИЕ ГЕНЕТИЧЕСКОГО АЛГОРИТМА ДЛЯ ПОИСКА ОПТИМАЛЬНОЙ ТРАЕКТОРИИ НАБЛЮДАТЕЛЯ
Д.В. Степанов, А.А. Шалыто
Предложена модификация генетического алгоритма, предназначенная для поиска оптимальной траектории наблюдателя в задаче оценивания параметров движущейся цели по угловым измерениям (Bearing-Only Target Motion Analysis, BOTMA). Возникающие проблемы сходимости и попадания в локальные экстремумы, характерные для генетических алгоритмов, решены с помощью идей, заложенных в оптимизационные методы случайного восхождения и имитации отжига. Для ускорения работы алгоритма отбор особей производился на уровне хорошо приспособленных схем. Ключевые слова: генетические алгоритмы, метод случайного восхождения, метод имитации отжига, теорема схем Холланда, гипотеза строительных блоков Гольдберга, расстояние Хэмминга, задача BOTMA.
Введение
В задаче определения параметров движения цели по угловым измерениям BOTMA одной из наиболее важных проблем является выбор движения наблюдателя, приводящего к получению качественных оценок [1, 2]. Параметры движения цели оцениваются с помощью линейных регрессионных моделей. Качество получаемых оценок характеризуется их ковариационной матрицей, выражающейся через матрицу плана регрессионной модели. Элементы матрицы плана, в свою очередь, зависят от параметров движения наблюдателя. Задача планирования эксперимента в случае регрессионной модели может быть сформулирована как задача выбора матрицы плана, доставляющей оценки с минимальной дисперсией. Таким образом, задачу планирования эксперимента в рассматриваемом случае можно представить в виде оптимизационной задачи поиска варианта движения наблюдателя, при котором минимизируется ковариационная матрица регрессионных оценок [3].
Будем рассматривать вариант кусочно-линейного движения наблюдателя с постоянной скоростью. При этом траектория наблюдателя может быть представлена в виде последовательности углов его поворота. Подобную оптимизационную задачу имеет смысл решать с помощью генетического алгоритма (ГА) [3], однако в случае большого числа «изломов» (галсов) траектории проявляется одна из частых проблем ГА, состоящая в попадании решения в локальные экстремумы.
Как показал анализ траекторий [3], в оптимальных многогалсовых траекториях часть галсов «сливается» вместе, и фактическое их число оказывается значительно меньше первоначально заданного. Объяснением данному факту является то, что, помимо осуществления маневров, для получения хороших оценок наблюдатель должен также производить измерения из точек, как можно более далеко разнесенных друг от друга. Траектория, удовлетворяющая данному требованию, не может быть сильно изломана. Оптимальные траектории, построенные из большого числа частично «сливающихся» галсов, имеют одно важное свойство - они могут содержать галсы различной длительности, что в ряде случаев приводит к улучшению оценок.
В настоящей работе предлагается модификация ГА, позволяющая находить оптимальные решения в задаче поиска траекторий наблюдателя с различными по длительности галсами.
Планирование эксперимента в задаче BOTMA
В активном эксперименте, каким является задача BOTMA, матрица плана регрессии выбирается так, чтобы повысить качество оценок. Ковариационная матрица управляется параметрами движения наблюдателя. Элементы ковариационной матрицы тем меньше, чем «больше» информационная матрица. В качестве характеристики матрицы будем рассматривать определитель матрицы.
Обозначим через 5 функцию от определителя информационной матрицы. В рассматриваемой постановке задачи в качестве управляющих параметров ограничимся только выбором углов поворота наблюдателя Cj eQC = [aMIN, aMAX), 0°<aMIN < aMAX < 360°, на заданном числе NLEG галсов. При отсутствии ограничений на угол поворота на галсе естественно выбрать для рассмотрения aMIN = 0°, aMAX = 360°. Формально оптимизационная задача состоит в поиске
argmax (5). (1)
(Ci,...,CNLBa ))leg
Пространство QCiLBG представляет собой пространство траекторий. Выбор возможных Cj ограничим дискретным множеством, состоящим из точек вида Cj = j • ДС, где ДС - шаг дискретизации. Число раз-
вариант°в углов поворота наблюдателя определяется как NALPH Обозначим ЧИсло
галсов Nleg . Тогда общее число траекторий, состоящих из NLEG галсов с NALPH вариантами возможных углов поворота, равно NALpH . Каждой траектории сопоставим угол, рассчитываемый по формуле
. . М1ЕО Д
ф^,...,^ _1)=С0 + х с,. • ы:1РН +-м-^н.
1=1 2
Интерпретируя 5 как функцию от ф , изобразим график функции 5 (рис. 1).
-20
Рис. 1. Способ изображения функции 5 при шаге дискретизации 30°
Сформулируем (1) как задачу поиска угла ф , в направлении которого функция 5 наиболее сильно удаляется от нуля. Как видно из рис. 1, график функции 5 имеет МА1РН корней. Они соответствуют прямолинейным траекториям наблюдателя, на которых определитель информационной матрицы обнуляется. Корни делят весь график функции 5 на МА1РН сильно изрезанных «лепестков». Каждый из них отвечает
выбору первого угла поворота траектории наблюдателя. Лепестки обладают множеством локальных максимумов, многие из которых имеют близкие значения. Рис. 1 дает представление о том, насколько сложной оказывается задача оптимизации (1).
Генетические алгоритмы и гипотеза строительных блоков
ГА представляет собой оптимизационный алгоритм, имитирующий процесс биологической эволюции [4]. Основной теоремой теории ГА, обосновывающей эффективность алгоритма, является теорема схем (шаблонов) Дж. Холланда, доказанная в 1975 г. [5, 6].
Схемой (шаблоном) Н называется подмножество множества генотипов, допустимое в данной популяции, заданное в виде хромосомы с зафиксированными значениями некоторых генов. Число «зафиксированных» генов схемы называется порядком, а расстояние между крайними зафиксированными генами - определяющей длиной. Функцией приспособленности ^ (Н) схемы называется среднее значение
функций приспособленности всех ее генотипов [5, 6]. Схемы с функцией приспособленности выше среднего по популяции, малым порядком и малой определяющей длиной принято называть строительными блоками.
Из теоремы схем следует, что шанс увеличить свое представительство в популяции следующего поколения имеется у схем малого порядка, малой определяющей длины и высокой приспособленности -у строительных блоков. В 1989 г. Д. Гольдбергом [6] была высказана гипотеза строительных блоков, состоящая в том, что ГА ведет поиск решения на множестве строительных блоков [7].
С практической точки зрения гипотеза Д. Гольдберга состоит в том, что наибольшая эффективность от применения ГА ожидается в тех задачах, в которых возможно выделение строительных блоков.
В задаче поиска оптимальной траектории (1) хромосомой длины Ышк = Ы1ЕО является последовательность углов поворота наблюдателя в моменты смены галса, например, при МУАК = 6 хромосома может принять вид С = (С0,С1,С2, С3,С4,С5) = (270°,0о,180о,0о,0°,210о). Одна из возможных схем Н , которой принадлежит приведенная выше хромосома, будет задаваться выражением Н = (*,0, *,0,0, *). Порядок приведенной схемы равен трем, а определяющая длина схемы также равна трем: 5 - 2 = 3.
В следующем разделе приводится модификация ГА, обеспечивающая решение оптимизационной задачи (1), использующая идею гипотезы строительных блоков.
Генетический алгоритм с эмиссией интронов
Сформулируем основные отличия предлагаемого метода от классической реализации ГА.
- Интенсивность применения генетических операторов привязывается к мере разнообразия популяции, меняющейся от поколения к поколению (аналогия с методом имитации отжига).
- Вероятность сохранения изменений, возникающих при применении генетических операторов, ставится в зависимость от наблюдаемого изменения функции приспособленности (аналогия с методом случайного восхождения).
- Выделяется ген (интрон), высокое содержание которого характеризует приспособленные хромосомы. На начальном этапе работы алгоритма осуществляется масштабное искусственное добавление данного гена в хромосомы.
- Производится поиск хорошо приспособленных схем (в согласии с гипотезой Д. Гольдберга) путем анализа расположения интронов в хромосомах - выделения интронных островов и их дополнений -экзонных островов.
- Вводится ряд новых операторов на экзонных островах, реализующих перебор траекторий-кандидатов.
Для измерения степени разнообразия популяции будем использовать расстояние Хэмминга dн (А, B), определяемое как число позиций, в которых соответствующие символы двух строк A и B
одинаковой длины различны. Обозначим через dн (О) среднее расстояние Хэмминга по популяции О . Используя идею метода имитации отжига [8], интенсивность применения генетических операторов поставим в зависимость от показателя dн (О).
Предполагая, что выпрямленные траектории будут иметь большую приспособленность, чем сильно изломанные, выскажем гипотезу о существовании строительных блоков (или хорошо приспособленных схем), состоящих из последовательностей нулей. При таком рассмотрении нас будет интересовать разделение генов на нулевые (интроны) и все остальные (экзоны).
Сделаем необходимые замечания относительно терминов «интрон» и «экзон». Термин «интрон» заимствован из биологии, где он обозначает некодирующий участок хромосомы, задающий пространственную структуру хромосомы. Термин «экзон» в биологии обозначает кодирующий участок хромосомы. Термин «интрон» используется также в области, родственной ГА - в генетическом программировании, где он обозначает неэффективный фрагмент кода, изменения в котором не влекут изменения в работе программы [9].
Оператор эмиссии интронов Ei создадим по аналогии с оператором мутации. Каждый ген Ck хромосомы C может мутировать в интрон (соответствующую новую хромосому обозначим C 0'к) с вероятностью pi, зависящей от среднего расстояния Хэмминга по популяции и изменения функции приспособленности. Оператор эмиссии интронов может быть записан в виде
Е,( +1,Ck е C) = I 0 р ^,0р, + = 1,
1 ' К, Ч, (к,0) "
где = 1 - р, - вероятность того, что данный ген не будет изменен оператором Е,: Г dн (Ор
(C ,о) =
1—
N
1 v VAR У
(C )> F (c0
Обозначим через MXQ популяционную матрицу, строками которой будут хромосомы особей, прошедших отбор. Если через NX обозначить число особей, проходящих отбор, то размерность матрицы Mx Q записывается как NX х NVAR . Далее, матрице Mx Q сопоставим матрицу M^Q , элементами которой являются индикаторы принадлежности к группе интронов «0» или экзонов «1». Отбор строительных блоков будем производить с помощью матрицы M^Q. Введем еще одно определение: назовем локусом
loci (j, Q) позицию гена в хромосоме. Задача отбора блоков состоит в определении локусов, в которых должны быть расположены интроны. Такие локусы будем называть далее интронными. Рассмотрим j -й столбец матрицы M^Q. Он содержит все варианты заполнения локуса loci (j, Q) в текущей популяции Q . Введем в рассмотрение среднее расстояние Хэмминга в локусе j популяции Q , характеризующее степень определенности относительно того, интрон или экзон должен располагаться в данном локусе:
dH (j, Q) = N (N 1) Z dH (AВ) .
NX v X _ ^ AiBsloci(j, Q)
Максимальное значение dH (у, О) равно NX/2 и соответствует ситуации, при которой в данном локусе находится поровну интронов и экзонов («0» и «1»). Минимальное значение dH (у, О) равно нулю и отвечает ситуации, в которой в данном локусе находятся либо только интроны, либо только экзоны. Введем в рассмотрение расстояние dH ,+1 (у, О), учитывающее информацию о среднем расстоянии Хэм-минга в локусе, содержащуюся в предшествующих поколениях:
t• Nx N -1)-dн,((, О) + X dн (А,В)
dн м (У, О) =-----, ^■О)-.
н' ^ ' (+l)•Nx N -1)
Используем аналогию, подсказанную моделями статистической физики: будем интерпретировать
dн ,+1 (У, О) как энергию системы частиц, а dн (О) - как температуру. Обратим внимание на то, что в
рассматриваемой системе с ростом времени (числа поколений) температура dн (О) (при правильно работающем алгоритме) будет снижаться. Наиболее устойчивое состояние системы характеризуется минимумом энергии, что в данном случае соответствует ситуации определенности в распределении интронов и экзонов в локусах. Дискретизированная энергия системы Ег =|^н ,+1 (у, О^ (где ] - целая часть
числа) в каждом локусе может принимать значения из отрезка [0; Nx/2]. Обозначим через Nz число частиц, находящихся на энергетическом уровне г , тогда значения Nz/Nx описываются распределением Максвелла-Больцмана:
( Ч_ N^ _ gz exp (-EjdH (Q)) Pmb (z) _ ,, _
Nx Zg^ exp(-Ew/dH (Q))
w
Будем называть энергетически неустойчивым значение элемента матрицы M^Q (к, j), соответствующее локусу loci (j, Q) в хромосоме с номером к е [1, NX ], если данное значение относится к менее многочисленному множеству при разбиении на интроны-экзоны. Например, если в локусе loci (j, Q) большинство составляют интроны, а MX Q (к, j) является экзоном или наоборот, то MX Q) (к, j) является энергетически неустойчивым значением. Рассматривая последовательно локусы loci (j, Q) в хромосомах с номерами к е [1, NX ], будем изменять энергетически неустойчивые значения MX Q (к, j) (находящиеся на энергетическом уровне z ) на противоположные с вероятностью 1 - pMB (z )• pCF (к), где pCF (к) -дискретное распределение вида
Pcf (к)_8(фк8(ф„)] . (2)
Это распределение отдает приоритет в сохранении энергетически неустойчивых состояний более приспособленным хромосомам. Вышеописанная операция, проведенная для всех энергетически неустойчивых элементов матрицы MX (к, j), преобразует исходную матрицу в некоторую новую матрицу
MX Q) (к, j), с помощью которой далее будет построена маска mQ (вектор размерности NVAR) для оценки интронных локусов:
m ( ■ V к е[1; Nx ] Mi^(к, j) _ 0; [1, 3 к е [1; Nx ]: MMХ£(к, j) * 0.
Далее, как только определены интронные острова, поиск оптимального решения будет происходить путем изменения генов в экзонных локусах. Определим четверку операторов, производящих следующие изменения:
1. оператор очистки - заменяет одиночный экзон на экзонном острове на интрон;
2. оператор сдвига - сдвигает весь экзонный остров на некоторое случайное число позиций вправо или влево;
3. оператор склейки - суммирует управляющие воздействия экзонного острова и помещает суммарное воздействие в один из локусов данного экзонного острова, согласно распределению (2);
4. оператор перемешивания - суммирует управляющие воздействия экзонного острова и распределяет суммарное воздействие по случайному числу локусов данного экзонного острова, согласно распределению (2).
После каждого цикла применения операторов 1-4 запускается оператор уборки мусора, который будет удалять противоположные по направлению управляющие воздействия, например, поворот на 30° и
следующий за ним поворот на минус 30°. В отличие от операторов 1-4, действие оператора уборки мусора не ограничивается границами экзонного острова. Изменения, производимые данными операторами, будем сохранять только в том случае, когда это приводит к увеличению приспособленности особи. Для ускорения работы алгоритма будем применять операторы 1-4 с вероятностью, пропорциональной частоте принятых изменений, произведенных соответствующим оператором.
Укажем изменения, которые необходимо внести в алгоритмы работы операторов рекомбинации и мутации, чтобы их действия были согласованы с алгоритмом отбора строительных блоков. В классическом ГА сохранению блоков большего порядка препятствует оператор мутации. Соответственно, чтобы блоки смогли выжить в ранних поколениях, мутация должна быть низкой в начале и может увеличиваться по мере снижения среднего расстояния Хэмминга по популяции. Действительно, при появлении в ней выживающих из поколения в поколение блоков dн (О) будет уменьшаться (тем больше, чем выше порядок блоков). В более поздних поколениях (в популяциях с малым dн (О)) увеличение числа мутаций
может помочь ускорить отбор более приспособленных блоков, правда, при следующем условии: будут разрешены только мутации, приводящие к увеличению значения функции приспособленности особи. Данное условие аналогично принципу, лежащему в основе метода случайного восхождения [8]. Зададим формулу вероятности мутации особи C (теперь зависящую от времени):
dн (О)
Pm ( + 1, C ) =
1—
N
VAR
'7{F,+1 (C)>Ft (C)} •
Рассмотрим механизм, который позволяет производить эмиссию интронов непосредственно в ходе применения оператора рекомбинации. Этот оператор, применяемый совместно с оператором рекомбинации, во избежание путаницы будем называть оператором выпрямления R . Он должен обнулять часть генов в хромосомах. Наиболее приспособленная особь в популяции (после применения генетических операторов) с большей вероятностью будет содержать искомые строительные блоки. По этой причине имеется заинтересованность в распространении ее генотипа. Из всех возможных вариантов оператора выбора пары более всего подходит вариант вожака стада (herd leader breeding [5]), при котором наиболее приспособленная особь образует пары со всеми остальными. Поскольку при скрещивании желательно не разрушать уже полученные блоки, выберем одноточечный вариант скрещивания, при котором родительские хромосомы рвутся в случайной точке и затем комбинируются для образования хромосом потомков.
Суммируя изложенное выше относительно операторов одноточечного скрещивания X и выпрямления R , приведем формальное описание совместного действия этих операторов:
. dn И.
(XR ) =
X,
R, a<Katt
nvar
du (A)
nvar '
где a - равномерно распределенная на [0,1] величина, а км - коэффициент затухания, принимающий
неотрицательные значения.
Описанный выше алгоритм будет быстро сходиться к экстремуму, отвечающему найденным блокам с высокой приспособленностью. Для того чтобы увеличить число перебираемых вариантов, реализуем алгоритм как островной [4]. Когда на острове с не самой высокой приспособленностью происходит вырождение, будем удалять популяцию данного острова и переселять наилучшую особь с лучшего острова, а остальные особи будем генерировать случайным образом, оставляя нетронутыми те гены, которые были нулевыми у наилучшей особи. Тем самым будут перенесены на обновленный остров выращенные строительные блоки с лучшего острова, обеспечивая запуск генетического отбора на данном острове с более высокого уровня приспособленности.
Предложенную модификацию ГА будем называть ГА с эмиссией интронов (IEGA - Intron Emission GA). При этом исходный алгоритм ГА будем обозначать GA.
Результаты моделирования
Приведем результаты работы GA и IEGA. На рис. 2 изображены графики функций приспособленности лучших особей (по поколениям) для GA (без мутации и с мутацией) и IEGA, из которых видно, что модифицированный алгоритм попадает на более высокий уровень приспособленности особей. Результаты работы алгоритмов приведены в таблице. Ниже представлено их сравнение:
- IEGA генерирует в качестве оптимальных траектории со значительно меньшим числом галсов, чем GA (4 против 18-19);
- функция приспособленности особей IEGA выше;
- разброс значений функции приспособленности особей IEGA ниже (среднеквадратическое отклонение (СКО) результатов IEGA в три-пять раз ниже, чем у результатов GA);
- 1БвЛ требует для сходимости в среднем больше поколений, чем вЛ без мутации. вЛ с мутацией не сходится за 100 поколений;
- разброс числа различных решений 1БвЛ ниже, чем у вЛ (16 против 30).
21,50 21,00 20,50 20,00 19,50 19,00 18,50
0
■е S
т °
1 °
s ^
,——*
1 3 5 7 9 11 13 15 17 19 21 23 25 27 йА Г,-!-; г.-,уг.;-д.1-.1; ОЛ,-.1у1 и',Ъ; пил
номер поколения
Рис. 2. Графики функций приспособленности лучших особей для ЗА и !ЕЗА
Алгоритм Число галсов 5 СКО(5) Среднее число поколений Общее число решений Число различных решений
GA (без мутации) 19 19,95 0,26 7 30 30
GA (частота мутации 0,05) 18 20,41 0,15 100 30
IEGA 4 21,06 0,05 20 16
Таблица. Сравнение результатов работы GA и IEGA Заключение
Предложенный генетический алгоритм с эмиссией интронов позволяет находить близкие к оптимальным траектории маневрирования наблюдателя, решающего задачу оценивания параметров движущейся цели. Изменения, внесенные в логику работы генетических операторов рекомбинации и мутации, и осуществление отбора на уровне высоко приспособленных схем, увеличили вероятность нахождения глобального оптимума функции приспособленности по сравнению с классическим вариантом генетического алгоритма.
Литература
1. Ince L., Sezen B., Saridogan E., Ince H. An evolutionary computing approach for the target motion analysis (TMA) problem for underwater tracks // Expert Systems with Applications. - 2009. - V. 36. - № 1. -Р. 3866-3879.
2. Cadre J.P., Jauffret C. Discrete-Time Observability and Estimability for Bearings-Only Target Motion Analysis // IEEE Transactions on Aerospace and Electronic Systems. - 1997. - V. 33. - № 1. - Р. 178-201.
3. Степанов Д.В. Использование генетического алгоритма для нахождения оптимального маневра в задаче N-пеленгов // III Всероссийский конкурс молодых ученых. - Миасс, 2011. - 11 с.
4. Haupt R.L., Haupt S.E. Practical genetic algorithms. - John Wiley & Sons, 2004. - 253 р.
5. Reeves C.R., Rowe J.E. Genetic Algorithms: Principles and Perspectives. - Kluwer Academic Publishers, 2003. - 332 р.
6. Goldberg D.E. Genetic Algorithms in Search, Optimization and Machine Learning. - Addison-Wesley, 1989. - 412 с.
7. Spears W.M. Evolutionary Algorithms: The Role of Mutation and Recombination. - Springer, 2000. - 222 р.
8. Russell S.J., Norvig P. Artificial Intelligence: A Modern Approach. - Prentice Hall, 2003. - 1081 р.
9. Brameier M., Banzhaf W. Linear Genetic Programming. - Springer, 2007. - 315 р.
Степанов Денис Вячеславович - ОАО «Концерн «НПО «Аврора», мл. научный сотрудник,
denis_v_stepanov@hotmail.com Шалыто Анатолий Абрамович - Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, доктор технических наук, профессор, зав. кафедрой, shalyto@mail.ifmo.ru