УДК 519.6
ПРИМЕНЕНИЕ МЕТОДА ЧАСТИЦ В СТАЕ К ЗАДАЧЕ ПОИСКА ОПТИМАЛЬНОГО УПРАВЛЕНИЯ НЕПРЕРЫВНЫМИ ДЕТЕРМИНИРОВАННЫМИ СИСТЕМАМИ
А.В. ПАНТЕЛЕЕВ, Е.А. АЛЁШИНА
Рассмотрена задача поиска оптимального управления непрерывными детерминированными системами с помощью метода частиц в стае. Сформирован детальный алгоритм решения поставленной задачи и на его основе создано программное обеспечение, эффективность которого продемонстрирована на примере решения прикладной задачи поиска оптимального управления химическим реактором.
Ключевые слова: непрерывная детерминированная система, оптимальное управление, метод частиц в стае.
Введение
В данной работе рассматривается применение метода частиц в стае, предложенного в [1], к задаче синтеза оптимального управления непрерывными детерминированными системами. Разработана программа, позволяющая решать поставленную задачу при различных значениях параметров метода. С помощью программы было проведено тестирование метода на задаче оптимального управления химическим реактором [2], которое показало, что при использовании модификации метода частиц в стае с локальным поиском получается лучшее решение, чем при использовании эволюционного алгоритма.
1. Постановка задачи
Поведение модели объекта управления описывается дифференциальным уравнением
х = / (7, х(г х и(г)), г ^ ], (1)
где х - вектор состояния системы, х е Кп; и - вектор управления; и = (и1,...,и9) е и с К.9;
и - некоторое замкнутое выпуклое множество допустимых значений управления; г - непрерывное время; г е Т = [г0, ]; начальный г0 и конечный моменты времени заданы,
/(г,х,и) = (/1 (г,х,и),...,/п (г,х,и))Т - непрерывно дифференцируемая вектор-функция.
Начальное состояние системы задано
х(0 = х0. (2)
Правый конец траектории свободен.
Предполагается, что при управлении используется информация только о непрерывном времени г, т.е. применяется так называемое программное управление.
Множество допустимых управлений образуют кусочно-непрерывные функции и () со значениями во множестве и. В точках разрыва значение управления определяется как предел справа.
Множество допустимых процессов О(г0, х0) - это множество пар ё = ( х(-), и(-)), включающих траекторию х() и управление и(), и (г) е и, удовлетворяющих уравнению состояния (1) с начальным условием х(г0) = х0.
На множестве допустимых процессов определен функционал качества управления
I(d) = j f0 (t, x(t), u(t)) dt + F (x(tN )), (3)
где f0 (t, x, u), F (x) - заданные непрерывно дифференцируемые функции.
Требуется найти такую пару d* = (x*( ), u*( ))е D(t0, x0), что
I(d*) = min I(d). (4)
deD(to,xo) VV
Замечание. Известно, что путем введения дополнительной переменной можно свести задачу Больца (3),(4) к задаче Майера. Поэтому в дальнейшем будем рассматривать только решение задачи Майера, в которой критерий имеет следующий вид
1 (d) = F (x(tN )). (5)
О
Задачи (1), (2), (5), (4) назовем задачей I.
2. Стратегия поиска решения
Метод частиц в стае был подробно рассмотрен авторами статьи в [3,4].
Идея метода состоит в моделировании поведения стаи животных при поиске пищи. При реализации метода используются следующие основные принципы, взятые из природы:
- для достижения цели используется множество особей (частиц);
- каждая особь получает информацию о передвижении соседей;
- движение по направлению к цели осуществляется путем подражания действиям соседа;
- учитывается информация о цели (величина целевой функции), получаемая через органы чувств (например, запах пищи);
- используется собственный прошлый опыт.
При этом моделируются тенденции живого организма повторять успешное в прошлом поведение, особенно в случае локальной неудачи; подражать успеху остальных; двигаться за лидером, выбираемым среди соседей. Сначала рассмотрим задачу поиска минимума целевой функции многих переменных на множестве допустимых решений Б .
Пусть имеется стая из ЫР частиц с позициями {х1,х2,...,хЫР} с Б с Яп . На каждом этапе
(итерации) у частицы с номером у выбирается случайное число соседей Ш1 ~ ^[Ы7тт, Ы/тах ] из членов стаи, определяемое равномерным ( Я ) распределением на отрезке [Ы/тт, Ы/тах ], где Жтт - минимальное число соседей; Жтах - максимальное число соседей (не превышающее ЫР ). Среди соседей находится частица с номером Пу и позицией хП) ’к, которой соответствует наименьшее значение целевой функции. Каждая частица с номером у хранит информацию о своей лучшей позиции хЛк, достигнутой в течение к итераций. Скорость частицы у]’к+1 с учё-
~ 1 к ~ п, ,к
том текущего положения и скорости, своей наилучшей позиции х* в прошлом, позиции х 1 лидера среди соседей, вычисляется по формуле
V1 ,к+1 = (о\],к + аг1 (х1,к - х1,к)+ сорг2 (хп1 ,к - х1 ,к), (6)
где к - номер итерации; v1,k, х1,к - скорость и позиция частицы на к -й итерации соответственно; о - весовой параметр, характеризующий инерцию (память) частицы; а и р - параметры, учитывающие степень влияния лидера среди соседей и прошлого опыта частицы (весовые
коэффициенты) на её последующую скорость; г1 и г2 - случайные параметры, заданные равно-
Новое положение частицы находится суперпозицией текущего положения и перемещения частицы в единицу времени (т.е. скорости)
Если новая позиция х1,к+1 оказалась лучше, чем х1 к, необходимо в качестве наилучшей по, к + 1 , к + 1 зиции х^ принять х1, .
Расчёт новой скорости и новой позиции по формулам (6) и (7) производится для 1 = 1,...,ЫР, т.е. для всех частиц стаи. После этого выполняется следующая итерация, затем следующая, и так до выполнения условий окончания расчета. В качестве условия окончания обычно выбирается равенство числа итераций к некоторому наперед заданному числу К (максимальному числу итераций), но такое условие не гарантирует близости к минимуму целевой функции в конце расчета.
Для решения задачи I предлагается осуществить переход к задаче поиска оптимального управления дискретной детерминированной системой, а затем построить решение исходной задачи путем интерполяции значений в узлах сетки. Предлагается искать решение задачи (управление) в виде кусочно-постоянных вектор-функций.
Множество допустимых значений управления и для простоты в дальнейшем будем полагать равным прямому произведению заданных отрезков [ а, ], I = 1,..., д.
Выберем Ы точек из отрезка [г0, гы ], в которых будем искать оптимальное дискретное управление. Пусть г0 = 0 , гк = гк-1 + к , к = 1,..., N, к = ^Ы ^— шаг. В результате получим равномерную сетку {г0, г1 = г0 + И,..., гЫ}.
Непрерывному управлению и () будет соответствовать дискретное управление и(•••) = {и(0),и(1),...,и(Ы- 1)}={иг}Ы=01, для траектории х() получим дискретный аналог:
х(---) = {х0,х(1),...,х(Ы)} = {хг}Ы=0. Дифференциальные уравнения состояния объекта будем решать приближенно с помощью численных методов [5,6]. В табл. 1 приведены формулы для некоторых численных методов применительно к поставленной задаче (учитывался тот факт, что управление предполагалось кусочно-постоянным). После численной аппроксимации уравнений объекта получим задачу поиска оптимального управления дискретной детерминированной системой [7,8].
Будем искать оптимальное управление и (•••). Позиция частицы под номером у на к -й итерации будет представлять собой вектор-строку и1,к(•••) = (и1,к(0),и1,к(1),...,и1,к(Ы-1)). На каждой итерации будет формироваться новая позиция частицы
и1 ,к+1(—) = ( и1 ,к+1(0), и1 к+1(1),..., и1 ,к+1(Ы -1)), ее необходимо будет оценить с помощью критерия
(5), который и будет целевой функцией метода частиц в стае.
Для получения значения целевой функции потребуется:
а) найти траекторию системы х1 ,к+1(---) = (х0,х1 ,к+1(1),...,х1 ,к+1(Ы)), соответствующую управлению и1 ,к+1(---) , из уравнения состояния (1) с учетом начального условия (2);
б) вычислить значение критерия (5), соответствующее и1 ,к+1(---) и х1 ,к+1(---). Выпишем соответствующую формулу в явном виде
мерным законом распределения ^[0; і].
(7)
где ё],к+1 =(х1,к+1(-.-),и],к+1(—)); х]М1(г), I = 0,...,Ы, определяется численным методом решения уравнения состояния (1) и начальным условием (2)
,к+1/ч Гх],к+1(г-1)+Дх 1, г = 1,...,Ы, х +1(г) = ^ г-1’ (9)
[ х0, г = 0,
где приращение Д хг-1 определяется в зависимости от выбранного метода (табл. 1). При дальнейшем изложении с целью сокращения записи символический аргумент (• ) будем опускать.
Таблица 1
Название метода Расчетные формулы
Метод Эйлера х,+1 = X + А X = X + h/(ї, X, и)
Метод Эйлера-Коши Х+1 = X + к/0, X, и ^ ( / (ї, X., и.) + / (ї + 1, X..,, и.+, )) X* = X + А X = X + (, ї, ї) У2( , ^ ї+1^
Метод Рунге-Кутты третьего порядка ^+1 = xí +Аxí, А xí = 4(К!+3К3), К[ = к/(/ xí, их К2 = к/(ї+3 к X+1 К^ и X 2 2 К3 = к/(ї+3 К X+3 K2, и)
Метод Рунге-Кутты четвертого порядка xí+1 = X +А X, АX =1 (К + 2К2 + 2К3 + К4) , 6 К1 = к/ ^ X, и), К 2 = к/(ї +1 к, X +1К1, и ), 2 2 2 1 * КЗ = к / (ї +1 к, X. +1К 2, и,), ^ 2 2 К4 = к / (ї + к, X* + КЗ, и)
Метод Адамса четвертого порядка к xt+l = X+А X, А X = 24(55/ - 59/-+37/-2 - ^-з^ / = /( ї, X, и); /0, /1, /2, /3 рассчитываются по одному из одношаговых методов, например, по методу Рунге-Кутты четвертого порядка
3. Алгоритм применения метода частиц в стае
Шаг 1. Задать ЫР - количество частиц в стае; Ы/тт, Ы/тах - параметры выбора числа «соседей» в стае; К - максимальное число итераций; весовой коэффициент о, характеризующий инерцию; параметры а и Ь, используемые при вычислении скорости частицы.
Шаг 2. Сформировать стаю частиц:
- положить число итераций к = 0 ;
- в области D допустимых решений задать NP частиц с позициями u1,0,u2 0,...,uNP,°. Для
этого использовать равномерный закон распределения, т.е. для каждой частицы с номером j = 1,..., NP найти позицию uj0 = (u;0(0),..., uj,0(N -1)), полагая
u!’°(0) ~ R [ a, ъ, ] v • ^ u!’°(N -1) ~ R [ а, ъ, ], * = Ъ..., q;
- положить начальные скорости равными нулю: vh° = 0, j = 1,...,NP .
Шаг 3. Для каждой частицы j = 1,..., NP найти:
- случайное число NI1 соседей при помощи равномерного распределения на отрезке [NImin , NImax ] , ^ЗДи них найти частицу с номером и позицией uHj,k, которой соответствует наименьшее среди выбранных соседей значение целевой функции;
- наилучшую позицию uj,k в течение k итераций (при к = 0 положить uj0 = uj0).
Шаг 4. Для каждой частицы в стае ( j = 1,...,NP ) найти новую скорость и новое положение
по формулам (6) и (7). Для нового положения частицы u,k+1 вычислить значение критерия I (d],k+1) по формуле (8) с использованием формулы (9) для нахождения компонент вектора
состояния xj,k+1. Обозначим допустимый процесс, в который входит uj,k, через dj,k, djk = (Xjk, u jk). Если I (djk+1) < I (djk), запоминается новая наилучшая позиция частицы:
u j,k+1 = uj,k+1, I (d j,k+1) = I (djk+1).
Шаг 5. Проверка условия окончания работы алгоритма. Если k < K, полагаем k = k +1 и переходим к Шагу 3. Если k = K , расчёт окончен.
В качестве решения задачи поиска оптимального управления дискретной системой u = (u*(0),u*(1),...,u (N -1)) выбирается позиция частицы в стае с наименьшим значением целевой функции и траектория x* = (x0, x*(1),..., x*( N)), которая рассчитывается по уравнению состояния. Приближенным решением задачи оптимального управления непрерывной детерминированной системой будет ступенчатая функция со значениями u*(0), u*(1),..., u* (N -1) на соответствующих подынтервалах отрезка T = [t0, tN ].
4. Программное обеспечение
На основе разработанного алгоритма сформировано программное обеспечение на языке C# в среде Microsoft Visual Studio 2005. Разработан интерфейс пользователя, позволяющий вводить параметры постановки задачи, задавать параметры метода, использовать модификации алгоритма, анализировать полученный результат.
С помощью разработанного программного обеспечения решим задачу оптимального управления химическим реактором, описанную в [2]. Модель непрерывной детерминированной системы описывается восемью дифференциальными уравнениями: x = u4 - qx1 -17,6x x2 - 23x1 x6u3; x2 = u1 - qx2 -17,6x1 x2 - 146x2 x3; x3 = u2 - qx3 - 73x2x3; x4 = -qx4 + 35,2x1 x2 - 51,3x4x5; x5 = -qx5 + 219x2x3 - 51,3x4x5; x6 =-qx6 +102,6x4x5 - 23x1 x6u3;
х7 = -дх7 + 46х1 х6 и3;
х8 = 5,8(qx1 -и4)-3,7и1 -4,1и2 + д(23х4 + 11х5 + 28х6 + 35х7)-5и^ -0,099, где д = и1 + и2 + и4, 0 < t < ^ = 0,2 .
Начальные условия: х(0) = [0,1883; 0,2507; 0,0467; 0,0899;0,1804; 0,1394;0,1046; 0]Г .
Ограничения на управление: 0 < ) < 20 , 0 < и2^) < 6 , 0 < u3(t) < 4 , 0 < и4 (t) < 20 .
Критерий качества управления: I = х8(^). Необходимо максимизировать критерий.
В программной реализации будем решать эквивалентную задачу с критерием I = - х8 ),
который нужно минимизировать.
Для первоначального приближения будем использовать метод частиц в стае без модификаций, так как он позволяет получить решение почти сразу же, в то время как метод частиц в стае с локальным поиском [3,4] требует для расчетов достаточно много времени. Для численного решения дифференциального уравнения выберем метод Адамса.
В программе имеется возможность решить задачу многократно, а затем для полученных значений критерия вычислить математическое ожидание и дисперсию. Используем эту возможность для оценки работы алгоритма без локального поиска. В табл. 2 приведены результаты решения задачи (выборочное среднее и выборочная дисперсия) различными численными методами.
Таблица 2
Численный метод Адамса Рунге-Кутты 4 порядка Рунге-Кутты 3 порядка
Выборочное среднее 17,8306 17,8770 17,8197
Выборочная дисперсия 0,062306 0,090277 0,092359
Из табл. 2 видно, что в среднем наилучший результат дало использование метода Рунге-Кутты 4-го порядка, однако наименьшая выборочная дисперсия получилась при использовании метода Адамса. В то же время метод частиц в стае без модификаций для всех численных методов дал решение, которое существенно хуже, чем значение критерия 21,7574, приведенное в [2]. Теперь решим задачу с использованием локального поиска. Результат работы программы представлен на рис. 1, 2.
Рис. 1. Оптимальное управление
Рис. 2. Оптимальная траектория
Наибольшее значение критерия, найденное методом частиц в стае с локальным поиском, равно 22,7, что больше, чем наилучшее значение 21,7574, полученное в [2]. Это свидетельствует о том, что метод частиц в стае с локальным поиском для данной конкретной задачи управления химическим реактором работает лучше, чем эволюционный алгоритм и метод динамического программирования.
Заключение
Предложен алгоритм решения задачи поиска оптимального управления непрерывной детерминированной системой методом частиц в стае. Сформировано программное обеспечение, с помощью которого была решена задача оптимального управления химическим реактором. Результаты решения сравнивались с результатами, полученными с использованием эволюционного алгоритма. Лучший результат 22,7 (против значения 21,7574, которое дал эволюционный алгоритм в [2]) был получен при использовании модификации метода частиц в стае с локальным поиском.
ЛИТЕРАТУРА
1. Eberhart R.C., Kennedy J. A New Optimizer using Particle Swarm Theory, Proceedings Sixth Symposium on Micro Machine and Human Science, IEEE Service Center, Piscataway, NJ, 1995
2. Lopez Cruz I. L., Van Willigenburg L. G., Van Straten G. Evolutionary Algorithms for optimal control of chemical processes. Wageningen University, Netherlands.
3. Пантелеев А.В., Алёшина Е.А. Разработка алгоритмического и программного обеспечения метода частиц в стае // Научный Вестник МГТУ ГА, серия Прикладная математика. Информатика, №145, 2009.
4. Пантелеев А.В., Алёшина Е.А. Метод частиц в стае // Гл. 4 в кн.: Пантелеев А.В. Метаэвристические алгоритмы поиска глобального экстремума. - М.: Изд-во МАИ-ПРИНТ, 2009.
5. Формалев В.Ф., Ревизников Д.Л. Численные методы. - М.: Физматлит, 2004.
6. Киреев В.И., Пантелеев А.В. Численные методы в примерах и задачах. - М.: Высшая школа, 2006.
7. Пантелеев А.В., Бортаковский А.С. Теория управления в примерах и задачах. - М.: Высшая школа, 2003.
PSO APPLICATION FOR OPTIMAL CONTROL PROBLEM OF DETERMINISTIC CONTINUOUS-TIME SYSTEM
Panteleyev A.V., Alyoshina E.A.
We consider the problem of finding optimal control for deterministic continuous-time system. The detailed algorithm that used different updates of PSO was formed. On the base of developed algorithm we made software that allows us not only to solve the problem but also to display the result in an obvious visual form. The efficiency of software was proved on a problem of finding the optimal control of high-dimensional non-linear continuous stirred tank reactor.
Key words: deterministic continuous-time system , optimal control, particle swarm method.
Сведения об авторах
Пантелеев Андрей Владимирович, 1955 г.р., окончил МГТУ им. Н.Э. Баумана (1978), доктор физико-математических наук, профессор, заведующий кафедрой математической кибернетики факультета МАИ, автор более 120 научных работ, область научных интересов - методы синтеза оптимальных нелинейных систем управления, методы оптимизации.
Алёшина Екатерина Александровна, студентка МАИ, автор 3 научных работ, область научных интересов - методы оптимизации, численные методы.