CC BY
149
УДК 519.8
ПРИМЕНЕНИЕ МЕТОДА ДИФФЕРЕНЦИАЛЬНОЙ ЭВОЛЮЦИИ И ЕГО МОДИФИКАЦИЙ В ЗАДАЧЕ ПОИСКА ОПТИМАЛЬНОГО УПРАВЛЕНИЯ ДИСКРЕТНЫМИ ДЕТЕРМИНИРОВАННЫМИ СИСТЕМАМИ
А.В. ПАНТЕЛЕЕВ, И.Ф. ДМИТРАКОВ
Сформирован алгоритм применения метода дифференциальной эволюции и его модификаций в задаче поиска оптимального управления дискретными детерминированными системами, создана программная среда визуализации его работы, позволяющая оценить его сравнительную эффективность.
Ключевые слова: метод дифференциальной эволюции, оптимальное управление, дискретные системы, мета-эвристические методы.
Введение
Отличительной особенностью решения актуальных задач анализа и синтеза сложных детерминированных систем управления является высокая трудоемкость вычислений при использовании прямых методов оптимизации. Поэтому в настоящее время большое внимание уделяется метаэвристическим методам глобальной оптимизации [1-5], позволяющим получить решение «высокого качества» за приемлемое (с практической точки зрения) время. К метаэвристическим методам относятся и алгоритмы дифференциальной эволюции [1-5], основанные на анализе эволюционных процессов. Особенностью алгоритмов дифференциальной эволюции является использование различий между индивидами (значениями аргумента целевой функции), реализованное линейным оператором, называемым «дифференциация».
Целью данной работы является создание алгоритмического и программного обеспечения применения метода дифференциальной эволюции и его модификаций в задаче оптимального управления дискретными детерминированными системами, разработка комплекса программных средств, позволяющего получить характеристики тестируемого метода оптимизации и делающего процесс сравнения его с другими методами удобным и наглядным.
Методы решения поставленной задачи основываются на алгоритме, изложенном в [1], и модификациях, изложенных в [5].
1. Постановка задачи
Поведение модели объекта управления описывается разностным уравнением
х(к +1) = / (к, х(к), и(к)), к = 0,1,..., N -1, (1)
где х - вектор состояния системы, х е Кп; и - вектор управления, и(к) еи(к) с К9; и (к) = [а^к, Ьук ]х [а2к, Ь к ]х...х [адк, Ь>чк ] - множество допустимых значений управления;
к - дискретное время, к е Т = [0,1,..., N -1], число шагов N задано; / (к, х, и) - непрерывная
вектор-функция, /(к,х,и) = (У1 (к,х,и),...,/п (к,х,и))Т .
Начальное состояние системы (1) задано х(0) = хо. Правый конец траектории свободен. Предполагается, что при управлении используется информация только о дискретном времени к , т.е. применяется так называемое программное управление.
Множество допустимых процессов Ф (0, хо) - это множество пар ё = ( х(-), и(-)), включаю-
щих траекторию х(-) = {х0,х(1),...,х^)] и управление и(-) = {и(0),и(1),...,и^-1)} , удовлетворяющих уравнению (1) с начальным условием х(0) = х0 .
На множестве Ф(0, х0) определен функционал качества управления
N-1
I(ё) = ^ /’ (к, х(к), и(к))+F (х( N)), (2)
к=0
где /0 (к, х, и), F (х) - заданные непрерывные функции.
* / * * \
Требуется найти такую пару ё = (х (•), и (•)) е Ф(0, х0), что
I(ё ) = шт I(ё) . (3)
йеФ (0,х0)
2. Стратегия поиска решения
При решении задачи (3) будем искать оптимальное управление и *(•) и, как следствие, оптимальную траекторию х *(•). Пусть вектор V - блочный столбец, состоящий из значений
ТТ управления на каждом шаге V = (V0,...,Vk ,...,VN-1) = (и(0),..., и (к),..., и(М -1)) , где блок
Т Т
^к = <7иУг£, .. ,^9,к) = и(к) = (и1(к)^^Х..-ид(к)) - управление на к -м шаге,
к = 0,1,..., N -1.
Будем обозначать через I(V) значение критерия (2), соответствующее V. Для его вычисления сначала ищется траектория системы х() = { х(0), х(1),..., х^)} из уравнения состояния (1)
с известными значениями управления - компонентами вектора V и начальным условием. Далее вычисляется значение критерия (2).
Поясним основную стратегию метода. Сначала на множестве допустимых управлений ио = и(0)хи(1)х...хи(N -1) случайным образом генерируется конечный набор
С = {Ут,т = 1,2,...,векторов, называемый начальной популяцией, где ^ - размер популяции. Далее идет циклический процесс замены текущей популяции новой. Он заканчивается, когда количество сформированных популяций оказывается равным заданному максимальному числу популяций.
Для формирования новой популяции последовательно выбирается каждый элемент текущей популяции (он называется вектором-мишенью V1), и принимается решение, остается ли он в новой популяции или его заменяет сгенерированный специальным образом вектор-образец.
3. Алгоритм
Шаг 1. Задать размер популяции ^; весовой коэффициент F; параметр операции скрещивания СЯ; максимальное количество популяций М .
Шаг 2. Сгенерировать с помощью датчика равномерного распределения на множестве
ио = и(0)хи(1)х...хи^-1) начальную популяцию - ^ векторов V2,...,VNP . Положить т = 0.
Шаг 3. Положить к = 1 (вектором-мишенью выбран вектор V1 = V1).
Шаг 4. Из числа оставшихся членов популяции случайным образом выбрать три различных члена Vе1 V, Vе, отличающихся от вектора-мишени V1.
Шаг 5. Сформировать вектор Vе1: Vе1 = Vе + F(Vе - VЬ). Если в полученном векторе Vе1 координата VClk £ [а5 к, Ъ8к ], то сгенерировать с помощью равномерного распределения на отрезке [а^к, Ъ8к ] случайную величину и VClk приравнять ей.
Шаг 6. Сформировать вектор-образец VР.
Для этого положить г = 1, ] = 0 и выполнить следующие действия:
а) сгенерировать при помощи равномерного распределения на [0; 1] случайное число 3г ]-;
б) если г = 9 и ] = N -1, то ViРj = Vе; если 3г ] < СЯ, то ViР]■ = Vе; если 3г ] > СЯ, то VР- = V* ■■
в) проверить выполнение неравенства г > 9 :
- если оно выполнено, то положить г = 1 и перейти к шагу 6 г;
- если нет, положить г = г +1 и перейти к шагу 6 а;
г) проверить выполнение неравенства ] > N -1:
- если оно выполнено, то процедуру формирования вектора VР завершить;
- если нет, положить ] = ] +1 и перейти к шагу 6 а.
Шаг 7. Формирование новой популяции. Подсчитать значения функционала для вектора-образца VР и вектора-мишени V*. Сравнить полученные величины I(VР) и I(V*) :
а) если I(VР) < I(V*), то поместить вектор-образец VР в новую популяцию;
б) если I(VР) > I(V*), то поместить вектор-мишень V* в новую популяцию.
Шаг 8. Проверить выполнение неравенства к > ^ :
а) если оно выполнено, то перейти к шагу 9;
б) если нет, то положить к = к +1, V* = Vк и перейти к шагу 4.
Шаг 9. Проверить условие окончания:
а) если число сформированных популяций т < М , то положить т = т + 1 и перейти к шагу 3;
б) если т = М , то процесс закончить. Выбрать в полученной популяции вектор, которому соответствует наименьшее значение целевой функции, и считать его приближенным решением задачи.
Замечания
1. Обычно полагают ^ = (5 15) • N • 9; 0,2 < F < 1; 0 < СЯ < 1. При N • 9 > 10 для быстрого получения приблизительного решения рекомендуется брать 0,5 < F < 0,8, для улучшения точности - увеличить М совместно с понижением F : 0,2 < F < 0,4 . Также рекомендуется начинать с СЯ = 0,9 или СЯ = 1 перед тем как пробовать СЯ = 0,1 или СЯ = 0 .
2. Если при использовании рекомендованных параметров полученные результаты не обеспечивают достаточной точности, рекомендуется применить метод еще раз с теми же параметрами или задать большее значение ^ или (и) М .
3. В пунктах 3-5 осуществляется операция мутации, а в пункте 6 - скрещивания.
4. На шаге 5 алгоритма реализуется оператор «дифференциация».
5. Изложенное в [5] обобщение схем метода дифференциальной эволюции также можно распространить на решение задачи (3). На шаге 5 предлагается использовать формулу
Vе1 = р + F 3,
где р - базисный вектор; 3 - вектор различий (направление поиска); F - коэффициент влияния.
Вектор получается двумя способами:
а) в популяции случайным образом выбирается элемент без учета информации о значении критерия (2);
б) в популяции находится элемент Vbest с наилучшим значением критерия (2).
Вектор 3 определяется также двумя способами:
а) различиями между случайно выбранными элементами (двумя, четырьмя и т.д.);
б) направлением (оно принимается за «хорошее» направление поиска, имитирующее градиентное, поскольку характеризуется различиями между наилучшим элементом Vbest и случайно выбранным из популяции).
В результате получены четыре группы методов.
Первая группа. Базисный вектор и направление поиска определяются по случайно выбранным элементам текущей популяции без учета значений критерия (2):
а) в популяции выделяются два класса Си С , содержащих пи п элементов соответственно;
1 п 1 П
б) находятся барицентры каждого класса: БЯ(' = — '*ТУ1 , БЯ(" = — '*ТУ1, где V1 - эле-
п г=1 п г=1
менты популяции, входящие в данные классы;
в) базисный вектор р находится одним из двух способов: Ре С' и С" или р£ С' и С";
г) сгенерированный вектор находится по формуле
Vе1 =р + F (ВЯс» - ВЯс'), где F выбирается из множества (-1, 0) и (0, 1) .
Вторая группа. По сравнению с формулой п. «г» первой группы методов в качестве базисного вектора выбирается наилучший элемент Vbestc' из класса С'
Vе1 = Пе*с + F (ВЯс- - ВЯ,у ).
Третья группа. Вектор различий формируется с учетом информации о величине критерия (2):
а) в популяции выбираются два класса С+ и С- с числом элементов п+ и п- соответственно так, что значение критерия для каждого элемента из класса С+ меньше значения критерия
для любого элемента из класса С- : I(У1)<I(У]), г = 1,...,п+;] = 1,...,п-;
тах
б) в каждом классе находятся «минимальный» и «максимальный» элементы: УC+m, Vc V(-m , V(ШaX , для которых выполняются условия
I^0“) < IУ) < I(VC+ax) "Vг е С+; I^0™) < IУ) < I(УCшax) "Vг е С-;
в) в каждом классе находятся приращения ВЯ+ = УC+1n - VC+ax, ВЯ- = V(C—n - V|C—X и сред-
0 п ВЯ+ + ВЯ-
нее приращение ВЯ8 = -
2
1 п+ 1 п-
г) находятся барицентры классов С+ и С- : ВЯс =— '^Vг , ВЯС- =— '^Vг ;
П + . , п . ,
+ 1=1 - 1=1
д) сгенерированный вектор находится по формуле
Vе1 = ВЯс+ + F (ВЯс+ - ВЯС- + ВЯ8 ) .
Заметим, что по сравнению с классической формулой параметр F может принимать как положительные, так и отрицательные значения.
Четвертая группа. В формуле п. «д» третьей группы в качестве базисного вектора используется наилучший элемент Vbestc+ в классе С+
Vе1 = ^«0+ + F ( ВЯс+ - ВЯс_ + ВЯ5).
Выделение описанных четырех групп позволяет в рамках каждой группы выделять еще несколько различных стратегий поиска.
4. Программное обеспечение
Разработана программа поиска оптимального управления, реализующая вышеописанный метод и основные его модификации. Среда разработки - Microsoft Visual Studio, язык программирования - C#.
Разработан интерфейс пользователя, позволяющий: вводить параметры постановки задачи; легко управлять параметрами тестируемого метода; анализировать эффективность его работы на модельных примерах; «видеть» процесс его работы; анализировать сходимость метода; просмотреть промежуточные результаты счета; удобно проводить сравнение выбранного метода с другими методами глобальной оптимизации.
Уравнение состояния оптимизируемых систем в большинстве случаев полагалось линейным
x(k +1) = Ax(k) + Bu(k), k = 0,1,..., N -1, (4)
где A - матрица размером n х n; B - матрица размером n х q .
Функционал качества управления считался квадратичным или линейно-квадратичным (для случая n = 1):
N-1/
I1(d) = 2 (xT (k)Sx(k) + uT (k)Qu(k)) + Lx(N) ( n = 1), (5)
k=0 N-1
I2 (d) = 2 (xT (k)Sx(k) + uT (k)Qu(k)) + xT (N)Lx(N), (6)
k=0
где S - матрица размером n хn; Q - матрица размером qх q; L - матрица размером nх n или число при n =1 .
5. Тестовые примеры
Пример 1. Даны модель объекта управления
x1 (k + 1) = x1 (k) + u(k) , x2 (k +1) = 2x1 (k) + x2 (k) , k = 0,1,
2 5
где x e R , | u(k) |< 10 , начальные условия x^0) = 2, x2(0) = 1 и функционал
1
I = 2 x12(k) + x22(k) + u2(k)
к=0
* *
Требуется найти минимальное значение функционала и оптимальный процесс (х (•), и (•)) на котором это значение достигается.
**
л _ 1 1 -Л х1 (°) = 2, х1 (1) = х2 (1) = i, х2
Аналитическое решение задачи: u (0) = -1, u (1) = 0 x1*(0) = 2, x1*(1) = x2*(1) = 1, x2*(0) = 1,
% * x2 (1) = 5 , x2 (2) = 7 , min I = 32 .
Результат работы программы представлен на рис. 1: и (0) = -1, и (1) = 0, min I = 32. В данном случае найденное с помощью программного обеспечения оптимальное управление совпадает с точным значением.
Рис. 1. Ввод параметров постановки задачи
Пример 2. Даны модель объекта управления x(k +1) = x(k) + u(k), начальное состояние
x(0) = 3, k = 0,1,..., N -1, где x е R, |u(k )| < 102, и функционал
N-1
I = ^ x2(k) + u2(k) ®min.
k=0
* *
Требуется найти минимальное значение функционала и оптимальный процесс (x (•), и (•)), на котором это значение достигается.
Будем искать численное решение, совпадающее с точным до 3-го знака после запятой, с помощью классического метода дифференциальной эволюции и четырех его модификаций, описанных в п. 5 замечаний к алгоритму. Приведем табл. 1, в которой покажем, сколько раз в процессе решения избранным методом понадобилось вычислять функционал качества управления для достижения требуемой точности.
Пусть N = 10. Аналитическое решение задачи min I = 60885/4181 @ 14,562.
Всеми методами (табл. 1) получено точное решение с требуемой точностью. Высокое, по сравнению с обычным методом, количество вычислений критерия при использовании групп
модификаций 3 и 4 обусловлено тем, что каждый раз при формировании вектора Vc1 требуется упорядочивать элементы популяции по возрастанию значения критерия, что вызывает многократное вычисление критерия. При использовании модификаций группы 2 для формирования
вектора Vc1 требуется находить элемент класса C', соответствующий наименьшему значению критерия, что также увеличивает количество вычислений критерия.
Таблица 1
Метод Параметры метода Количество вычислений функционала
NP F CR M дополнительные
Классический 45 0,4 0,9 190 1,7-104
1-я группа 70 -0,9 0,5 300 n = 10, п= 10 4,0-104
2-я группа 70 -0,9 0,5 200 п = 10, п = 10 17,0-104
3-я группа 70 -0,9 0,5 300 п+ = 10, п- = 10 46,0-104
4-я группа 50 0,5 0,5 500 п+ = 20, п- = 20 100,0-104
Можно заключить, что использование групп модификаций 1-4 связано с большими вычислительными затратами. Поэтому, если условия задачи не требуют высокой точности численного решения, можно рекомендовать использовать классический метод.
Пример 3. Даны модель объекта управления
xi(k )
x2(k)+
Ui (k ) Xi (k )
Xi (k +1) =
1 + 0,01- Ui (k )(3 + U2 (k)) x3(k +1) =
x2(k +1) =-x3(k )
1 + 0,01- Ui (k )(3 + U2 (k))
1 + Ui (k )(1 + U2 (k))
1 + 0,01-U2(k)(1 + U3(k)) ’
начальное состояние x(0) = [ 2; 5; 7] , к = 0,1,..., N -1, где x e R ; ограничения на управление: 0 < ^i(k) < 4, 0 < U2(k) < 4, 0 < щ(к) < 0,5; и функционал
I = x2 (N) + x22 (N) + x32 (N) +
+
f N-1
V N
2 x2 (k) + x2 (k) + 2u2 (k) 2 X3 (k) + 2u2 (k) + 2u| (k)
V k=0 Jv k =1
Требуется найти минимальное значение функционала и оптимальный процесс (x (•), и (•)), на котором это значение достигается. Для N = 100 аналитическое решение задачи, найденное в [8] minI @ 258,339; численные решения, полученные с помощью разработанного программного обеспечения, представлены в табл. 2.
Таблица 2
1
2
№ п/п Параметры метода min I
NP F CR M
1 2000 0,7 0,90 2000 358,787
2 2000 0,3 0,30 4000 258,341
3 2000 0,3 0,20 4000 258,340
4 2000 0,3 0,15 4000 258,340
Для выбранного значения N с помощью разработанного программного обеспечения получены решения, близкие к точному. В случае 1 задано большое значение параметра Е , что привело к увеличению погрешности.
Заключение
В данной работе создан алгоритм и программное обеспечение применения метода дифференциальной эволюции и его модификаций в задаче поиска оптимального управления дискретными детерминированными системами. Эффективность разработанного программного обеспечения продемонстрирована на модельных примерах.
ЛИТЕРАТУРА
1. Пантелеев А. В., Дмитраков И. Ф. Применение метода дифференциальной эволюции для оптимизации параметров аэрокосмических систем // Электронный журнал “Труды МАИ” - 2010. - № 37.
2. Пантелеев А.В. Метаэвристические алгоритмы поиска глобального экстремума. - М.: МАИ-Принт, 2009.
3. Дмитраков И.Ф. Анализ эффективности метаэвристических методов поиска глобального экстремума функций многих переменных // Актуальные проблемы авиационных и аэрокосмических систем: процессы, модели, эксперимент. - 2009. - Т. 14. - Вып. 2 (29).
4. Price K., Storn R., Lampinen J. Differential évolution - a practical approach to global optimization, Springer, 2005.
5. Feoktistov V. Differential Evolution In Search of Solutions, Springer, 2006.
6. Пантелеев А.В., Летова Т. А. Методы оптимизации в примерах и задачах. - М.: Высшая школа, 2008.
7. Пантелеев А.В., Бортаковский А.С. Теория управления в примерах и задачах. - М. : Высшая школа, 2003.
8. Luus R. Iterative Dynamic Programming, London, UK: Chapman & Hall/CRC, 2000.
APPLICATION OF DIFFERENTIAL EVOLUTION METHOD TO THE DISCRETE TIME DETERMINISTIC SYSTEMS OPTIMAL CONTROL PROBLEM
Panteleyev A.V., Dmitrakov I.F.
Detailed algorithms of application of the method of differential evolution and its modifications to the problem of optimal control of the discrete time deterministic system are suggested; the software environment for visualization of algorithms, allowing to estimate its comparative efficiency, is developed.
Key words: the differential evolution method, optimal control, the discrete time systems, the methaeuristic methods.
Сведения об авторах
Пантелеев Андрей Владимирович, 1955 г.р., окончил МГТУ им. Н.Э. Баумана (1978), доктор физико-математических наук, профессор, заведующий кафедрой математической кибернетики МАИ, автор более 130 научных работ, область научных интересов - методы синтеза оптимальных нелинейных систем управления, методы оптимизации.
Дмитраков Игнатий Филиппович, 1988 г.р., студент МАИ, автор 5 научных работ, область научных интересов - анализ и оценка эффективности сложных систем, методы оптимизации.
