CC BY
111
УДК 519.6
ПРИМЕНЕНИЕ МЕТОДОВ ПЧЕЛИНОЙ КОЛОНИИ И ПОИСКА ГАРМОНИИ К ЗАДАЧЕ НАХОЖДЕНИЯ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ ДИСКРЕТНЫМИ ДЕТЕРМИНИРОВАННЫМИ СИСТЕМАМИ
А.В. ПАНТЕЛЕЕВ, Е.А. АЛЁШИНА
Рассмотрена проблема применения методов поиска гармонии, роя пчел и искусственной пчелиной колонии к задачам оптимизации нелинейных дискретных детерминированных систем. Разработано программное обеспечение, с помощью которого произведено сравнение указанных методов на тестовом примере, содержащем систему управления сложной структуры.
Ключевые слова: метод пчелиной колонии, метод поиска гармонии, оптимальное управление.
Введение
Метаэвристические методы являются актуальным направлением теории оптимизации [1-6]. Использование таких методов целесообразно для решения прикладных задач оптимального управления летательными аппаратами (самолетами, ракетами, космическими аппаратами) на различных участках полета со сложной структурой математических моделей, ввиду того, что традиционные методы оптимизации не способны справиться с подобными проблемами. В данной работе на примере оптимизации нелинейной дискретной детерминированной системы управления [5] сравниваются метод поиска гармонии, метод искусственной пчелиной колонии и метод роя пчел [1-4].
1. Постановка задачи
Поведение модели объекта управления описывается разностным уравнением
х(г +1) = / (г, х(| X и (г )), (1)
где I - дискретное время, Iе Т = [0,1,...,N — 1]; х - вектор состояния системы, хе Я"; и - вектор управления, и еи(I) с Я ; и(I) - множество допустимых значений управления, для каждого значения I представляющее собой прямое произведение отрезков [а{(I),Ь(I)], г = 1, 2,...,q; число шагов дискретного времени N задано; /(I,х,и) = (/1(1,х,и),...,/п(I,х,и))Т - непрерывная вектор-функция.
Начальное состояние задано
х(0) = х0. (2)
Правый конец траектории х^) свободен. Предполагается, что при управлении используется информация только о дискретном времени I, т.е. применяется программное управление. Таким образом, рассматриваемая система управления является разомкнутой по состоянию.
Множество допустимых процессов .0(0, х0) - это множество пар ё = (х(), и()), включающих траекторию х() = {х0,х(1),...,х^)} и допустимое управление и() = {и(0),и(1),...,и^ — 1)} , где и(1) е и(I), удовлетворяющих уравнению состояния (1) и начальному условию (2).
На множестве допустимых процессов 0(0, х0) определен функционал качества управления
N—1
1 (ё) = I /0 (I, х(г), и (I))+^ (х( N)),
г=0
где f0 (t, x(t), u(t)), F (x) - заданные непрерывные функции. Требуется найти такую пару d* = (x * (•), u * (•)) е 0(0, x0), что
I(d*) = шт IШ) . (4)
dеD (0,!0)
Искомые элементы пары d*: траектория x*(•) = x*(1),...,x*(#)} и управление u *(•) = {а *(0),u *(1),...,u *(N-1)} называются соответственно оптимальной траекторией и оптимальным управлением.
2. Идея решения задачи
При использовании метаэвристических методов для решения задачи (1)-(4) будем оптимизировать управление u(•). Позиция частицы (пчелы) с номером ] на k -й итерации будет представлять собой вектор-строку u1’к =(uJ’k(0),u1’k(1),...,uJ’k^-1)). На каждой итерации будет формироваться новая позиция частицы ^,1с+1 =(^д+1(0),^д+1(1),...,^,k+1(N-1)), ее необходимо
будет оценить с помощью критерия (3), который будет целевой функцией. Для получения значения целевой функции потребуется: а) найти траекторию системы
x1 Л+1 =(x0,x1 д+1(1),...,x1 М1^)), соответствующую управлению uJ,к+1, из уравнения состояния
(1) с учетом условия (2); б) вычислить значение критерия (3), соответствующее ^,k+1 и x1 ’к+х. Выпишем формулу в явном виде
N-1
I (dJ,k+1) = 2 f0 (t, xJ,k+1 (t), u м+1 (t)) + F (xJ,k+1 (N)),
I=0
где d1 Д+1 =(x1 ^(О,^^(О) - процесс из множества допустимых процессов Ф(0,x0),
x]’к+1(t),t = 0,...,N определяется уравнением состояния (1) и начальным условием (2)
7 (t -1, xJ,k+1 (t -1), u J,k+1 (t -1)), t = 1,..., N,
XJ,k +1(t ) =
^0, t = 0.
3. Метод поиска гармонии
Стратегия поиска решения. Метод имитирует процесс импровизации музыканта-исполнителя [1]. В процессе поиска каждое решение из множества допустимых порождает значение функции с целью достижения глобального экстремума. При этом используются идеи метода имитации отжига, метода частиц в стае, стохастического градиента.
Метод поиска гармонии сконструирован таким образом, что процедура генерирования нового решения позволяет выходить из области притяжения локального экстремума.
Алгоритм
Шаг 1. Задать параметры метода:
hms - размер памяти гармонии;
hmcr - частоту выбора значений из памяти гармонии;
par - частоту выбора соседнего значения;
^ - вектор максимального изменения приращения;
K - максимальное число итераций.
Положить число итераций k = 0 .
Шаг 2. На множестве допустимых решений генерировать hms решений и1,...,ыНт!1, и1 =(и1 (0),и1 (1),...,и1 (М-1)). Используя уравнение состояния (1), получить допустимые процессы ё1 =(х1,и1), 1 = 1,...,hms и найти соответствующие значения целевой функции
I (ё1),..., I (ёНт). Сохранить их в памяти гармонии ( НМ )
' и\ • К f (и1) Л
„ hms и, V 1 hms • ип f (иhms) ,
HM =
Найти наихудшее решение в памяти uworst.
Шаг 3. Генерировать новый вектор unew. Для всех i = 1,...,n :
Шаг 3.1. Получить значения и' следующим образом:
- с вероятностью hmcr выбрать элемент из памяти HM с номером int[r (0,1) • hms] +1,
и' = и int[r (0,1) hms]+1
где r (0,1) - равномерно распределенное число от 0 до 1; int[] - операция нахождения целого числа;
- с вероятностью 1 - hmcr выбрать значение внутри промежутка [a{, bt ]: и' е [a{, Ь ].
Шаг 3.2. Если значение и' выбрано из памяти, то
- с вероятностью par найти
и- = и' + fWi • r(-1,1), где r (-1,1) - равномерно распределенное число на [-1,1];
- с вероятностью 1 - par положить
new /
и1 = и1.
ТТТ J Т"1 new worst
Шаг 4. Если новое решение и лучше наихудшего решения и в памяти, то заменить nworst на ипew. Иначе - не заменять.
Шаг 5. Проверка условий окончаний процесса поиска:
- если k < K , положить к = к +1 и перейти к шагу 3;
- если к = K, процесс завершить, найти наилучшее решение иЬеst в памяти и положить
и* = иЬы.
Для сравнения с другими методами использовалась модификация метода поиска гармонии, в которой величины fwi экспоненциально убывают.
4. Метод роя пчел
Стратегия поиска решения. Идея метода состоит в моделировании поведения пчёл при поиске нектара [2]. На основе информации, полученной от других пчёл, пчела может двигаться к предполагаемой точке экстремума; выбраться из области притяжения локального экстремума; исследовать ранее неисследованные области.
Алгоритм
Шаг 1. Задать максимальное число итераций K, число пчёл-разведчиков s, пороговое значение расстояния между пчёлами в начальный момент threshold, параметр области локального поиска А, числа пчёл B и P, посылаемых в «наилучшие» и «перспективные» области; числа b
и p = s - b отбираемых наилучших и перспективных значений целевой функции. Положить число итераций к = 0 .
Шаг 2. Сформировать начальное положение s пчёл-разведчиков при помощи равномерного распределения на множестве U. Если расстояние (евклидово) между двумя пчёлами-разведчиками меньше порогового значения threshold , остаётся пчела с большим значением функции. Недостающие пчёлы-разведчики генерируются заново.
Для всех пчёл-разведчиков подсчитываются значения функционала качества управления.
Шаг 3. Все значения критерия (для всех пчёл, находящихся на множестве U ) упорядочить по убыванию: I (d1 )> /(d2) > ... > I (ds). Из них выбрать b наилучших значений и p перспективных значений, наиболее близких к наилучшим.
Шаг 4. Выбранным на шаге 3 значениям ставятся в соответствие области локального поиска, каждая из которых представляет собой гиперкуб; центр области определяется координатами пчелы-разведчика, длины сторон равны 2А. Результатом шага 4 является нахождение b «наилучших областей» и p перспективных областей.
Шаг 5. В каждой из «наилучших областей» случайным образом (используя равномерное распределение) генерируется B новых решений, а в каждой из «перспективных областей» генерируется P новых решений. Для всех новых решений подсчитываются значения целевой функции.
Шаг 6. Проверяется выполнение условия окончания:
- если к < K -1, то положить к = к +1 и перейти к шагу 3;
- если к = K -1, процесс завершить и в качестве и * выбрать наилучшее решение из достигнутых.
При решении задачи поиска оптимального управления использовалась модификация метода, в которой области сжимаются с ростом числа итераций А = А • R, где 0 < R < 1. В начальный момент времени значение А = А0 задано.
Важная особенность алгоритма состоит в том, что новые решения, генерируемые в процессе поиска, не могут выйти за границы множества допустимых решений. Для этого достаточно при задании «наилучших» и «перспективных» областей корректировать их границы с учетом границ множества допустимых решений.
5. Метод искусственной пчелиной колонии
Стратегия поиска решения. В методе искусственной пчелиной колонии (ABC - artificial bee colony algorithm) [3; 4] в качестве основы взяты особенности поведения пчёл при поиске источников пищи (нектара). Основное отличие от метода роя пчел - другой способ генерации нового решения.
Алгоритм
Шаг 1. Инициализация алгоритма и генерация начальных позиций источников питания.
Шаг 1.1. Задать параметры метода: максимальное число итераций K; количество источников питания, количество рабочих пчёл, количество пчёл-наблюдателей SN; предельное количество итераций, по прошествии которых решение считается забытым, если оно не было улучшено, limit .
Шаг 1.2. Для каждой рабочей пчелы находятся координаты её положения
и/ = at + rand• (b -at), i = 1,...,n, j = 1,...,SN, где rand - равномерное распределённое случайное число на отрезке [0; 1].
Шаг 1.3. Для всех полученных решений (для всех пчёл-разведчиков) находятся соответствующие допустимые процессы ё1 =(х1,и1) и значения критерия I (ё1), у = 1,...,8Ы.
Шаг 1.4. Положить к = 0 .
Шаг 2. Получение новых решений для рабочих пчёл Шаг 2.1. Сгенерировать новые решения V1
V1. = и/ +ф1 • (и/ - и'), . = 1,...,п, 1 = 1,..., ,
где ф1 - равномерное распределённое случайное число на отрезке [-1; 1], номер г выбирается случайным образом из {1,...,SN} , но при этом г отличается от у.
Шаг 2.2. Рассчитать значения критерия I (ё}), где ё} =( х1, V1), у = 1,..., .
Шаг 3. Произвести селекцию среди новых решений V1 и существующих решений и1:
- если I (ё} )< I (ё}), заменить и1 на новое решение V1
,>j •
- иначе и j не заменять.
Шаг 4. Рассчитать вероятности выбора р1 для решений и1
I (d1)
р'= ■ j = 1,..,SN.
Я (<!■)
s=1
Шаг 5. Получение новых решений для пчёл-наблюдателей аналогично решениям рабочих пчёл на шаге 2.
Шаг 5.1. Сгенерировать новые решения v1
vj = и/+фг; • (и/- и'), i = 1,..., n, j = 1,..., SN, где ф1 - равномерное распределённое случайное число на отрезке [-1; 1], номер r выбирается случайным образом на основе вероятностей р1, но при этом r отличается от j .
Шаг 5.2. Рассчитать значения критерия I (d1), d} =( x1, v1), 1 = 1,..., SN.
Шаг 6. Произвести селекцию среди новых решений v1 и существующих решений и1, таким же образом, как на шаге 3.
Шаг 7. Генерация новых решений вместо «забытых». Решения, для которых количество итераций без изменения значения превысило limit , считаются «забытыми». Среди всех забытых решений выбирается одно с наихудшим значением целевой функции. Это решение исключается из рассмотрения, а вместо него генерируется новое случайным образом на множестве U :
и! = ai + rand• (bi -ai), i = 1,...,n, 1 = 1,...,SN, где rand - равномерное распределённое случайное число на отрезке [0; 1].
Шаг 8. Определить и сохранить в памяти лучшее к данному моменту решение иbest.
Шаг 9. Проверить условие окончания процесса поиска:
- если к < K -1 , то перейти к шагу 2.
- иначе процесс завершить. В качестве и * выбрать наилучшее решение иь<!st.
6. Программное обеспечение
На основе разработанного алгоритма сформировано программное обеспечение на языке C# в среде Microsoft Visual Studio 2005, включающее графический пользовательский интерфейс. С помощью программного обеспечения пользователь может выполнять следующие действия:
вводить параметры постановки задачи; задавать параметры метода; просматривать результат решения задачи в виде графиков; анализировать полученный результат; сохранять результат в памяти компьютера для последующего анализа.
7. Решение тестовых примеров
Пример 1. Рассмотрим задачу оптимального управления, описанную в [5].
Модель дискретной системы управления описывается конечно-разностными уравнениями:
хі0 )
x1(t +1):
1 + 0,01 u1(t) (3 + u2(t))
x2(t +1) = -
= x2 (t) + u1 (t)x1 (t + 1) ; 1 + u1(t) (l + u2(t) )
x3(t)
x3(t +1) = 3W
1 + 0,01 u2 (t) (l + u3 (t))
г = 0,1,..., N -1.
Начальные условия х(0) = [2 5 7]Г .
Ограничения на управление: 0 < и1(г) < 4, 0 < и2(г) < 4, 0 < и3(г) < 0,5. Критерий оценки качества управления
I = У x2(tN)+
У x2(t-1) + x22(t-1) + 2u32(t-1) її У x32(t-1) + 2u12(t-1) + 2u22(t-1)
Необходимо минимизировать критерий.
В [5] рассмотрены случаи N = 20, N = 50 и N = 100.
Разработанное программное обеспечение позволяет провести серию опытов и рассчитать статистические характеристики получаемой выборки значений. При проведении опытов запускалась серия из 100 решений (для N = 20 ) или 30 решений (для N = 50 и N = 100) одной и той
же задачи с одними и теми же значениями параметров. Для полученной выборки {/', 12,..., 1100}
_ 1 Ю0 _ ,=--
вычислялись выборочное среднее / =-------У / и среднеквадратическое отклонение О/ =у1 £ 100 ,
100
i =1
_ 1 100 _ 2
S100 = — У (Г -1) .
99 >
В табл. 1 сведены результаты серий опытов при разных значениях времени функционирования системы N, а также для сравнения приведены результаты из [5], полученные с использованием метода итерационного динамического программирования. Использовались следующие обозначения: HS (Harmony Search) - метод поиска гармонии; BA (Bees Algorithm) - метод роя пчел; ABC (Artificial Bee Colony) - метод искусственной пчелиной колонии; nf - требуемое количество вычислений значения целевой функции. Параметры методов подбирались таким образом, чтобы реальное количество расчетов значений целевой функции было приблизительно равно nf .
По результатам, представленным в табл. 1, делаем вывод, что в целом наиболее близкие к оптимальным решения получаются при использовании метода роя пчел, однако в случае ограниченных вычислительных ресурсов иногда оказывается лучше использовать метод поиска гармонии (случаи N = 50 и N = 100, когда количество вычислений значений функции не должно превышать 25000).
2
i =1
t=1
t=1
Таблица 1
N nf Метод Г Г наим Sr Значение из [5]
20 5000 HS 209,563 210,6242 0,524269 209,26937
BA 209,3048 209,9043 0,311183
ABC 212,8299 215,3463 0,887809
20 25000 HS 208,672 209,2445 0,216007
BA 208,4255 208,6678 0,124801
ABC 209,3878 209,9216 0,315238
20 100000 HS 208,4796 208,802 0,157133
BA 208,4369 208,6077 0,097741
ABC 208,471 208,6318 0,053726
50 25000 HS 243,0316 243,781 0,484859 240,917
BA 246,4545 247,9878 0,617528
ABC 253,5078 258,0076 1,948513
50 100000 HS 241,2325 241,8219 0,27996
BA 240,2685 240,4086 0,071366
ABC 248,8643 250,3288 0,687646
50 250000 HS 240,5671 240,9298 0,189076
BA 239,9008 239,975 0,034119
ABC 242,3185 242,8575 0,325716
100 25000 HS 275,5493 278,9272 1,515548 258,33922
BA 297,3412 301,2077 1,720014
ABC 313,3766 320,3505 3,291568
100 100000 HS 267,8703 269,732 0,936816
BA 259,8398 260,1668 0,233011
ABC 303,2857 307,3917 1,860283
HS 264,0494 265,1362 0,583137
100 250000 BA 258,1714 258,3441 0,09873
ABC 281,2116 284,0097 1,72876
Заключение
Рассмотренные метаэвристические методы позволяют решать задачи поиска оптимального управления нелинейными дискретными детерминированными системами сложной структуры. Наилучшие результаты дает применение метода роя пчел, однако следует отметить, что в [4] приведены данные, полученные одним из вариантов метода искусственной пчелиной колонии (описанным недостаточно подробно), который несколько отличается от алгоритма, используемого авторами в данной статье.
При решении тестовой задачи подбор параметров осуществлялся эмпирическим образом. Такой подход не гарантирует оптимальность подобранных значений. В связи с этим интерес представляет реализация адаптивных версий метаэвристических методов, которая позволит провести более корректное сравнение.
ЛИТЕРАТУРА
1. Lee K.S., Geem Z.W. A new meta-heuristic algorithm for continuous engineering optimization: harmony search theory and practice // Computer Methods in Applied Mechanics and Engineering, 194, pp.3902-3933, 2005.
2. Pham D.T., Ghanbarzadeh A., Koc E., Otri S., Rahim S., Zaidi M. The Bees Algorithm - A Novel Tool for Complex Optimization Problems // Proceedings of IPROMS 2006 Conference, pp. 454-461, 2006.
3. Basturk B. An artificial bee colony (abc) algorithm for numeric function optimization / D. Karaboga, In IEEE
Swarm Intelligence Symposium 2006, Indianapolis, Indiana, USA, May 2006.
4. Karaboga D., Akay B. A comparative study of artificial bee colony algorithm. Applied Mathematics and Computation, In Press, 2009.
5. Luus R. Iterative Dynamic Programming. - Monographs and Surveys in Pure and Applied Mathematics, London, UK: Chapman & Hall/CRC, 2000.
6. Пантелеев А.В. Метаэвристические алгоритмы поиска глобального экстремума. - М.: Изд-во МАИ-ПРИНТ, 2009.
BEES OPTIMIZATION AND HARMONY SEARCH APPLIED TO OPTIMAL CONTROL OF DISCRETE TIME DETERMINISTIC SYSTEM
Panteleyev A.V., Alyoshina E.A.
We consider the application of bees, artificial bee colony and harmony search algorithms to process of searching the optimal control of discrete deterministic system. The special software for solving such problems was developed. We performed the comparison of the methods considered on test problem containing control system of complex structure.
Key words: artificial bee colony algorithm, harmony search, optimal control.
Сведения об авторах
Пантелеев Андрей Владимирович, 1955 г.р., окончил МГТУ им. Н.Э. Баумана (1978), доктор физико-математических наук, профессор, заведующий кафедрой математической кибернетики МАИ, автор более 150 научных работ, область научных интересов - методы синтеза оптимальных нелинейных систем управления, методы оптимизации.
Алёшина Екатерина Александровна, окончила МАИ (2010), аспирантка МАИ, автор 14 научных работ, область научных интересов - методы оптимизации, численные методы.
