CC BY
49
ИНФОРМАТИКА, ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА И УПРАВЛЕНИЕ
DOI: 10.18698/0236-3933-2016-1-33-50
УДК 519.85:517.977.58
ПРИМЕНЕНИЕ ОБОБЩЕННОГО ИНВЕРСНОГО ИНТЕРВАЛЬНОГО МЕТОДА ГЛОБАЛЬНОЙ УСЛОВНОЙ ОПТИМИЗАЦИИ В ЗАДАЧЕ ПОИСКА ОПТИМАЛЬНОГО ПРОГРАММНОГО УПРАВЛЕНИЯ
А.В. Пантелеев, В.Н. Пановский
Московский авиационный институт (национальный исследовательский университет), Москва, Российская Федерация e-mail: avpanteleev@inbox.ru; panovskiy.v@yandex.ru
Разработано алгоритмическое и программное обеспечение обобщенного инверсного интервального метода глобальной условной оптимизации, а также метод его применения для решения задачи нахождения оптимального программного управления нелинейными детерминированными непрерывными динамическими системами. Разработана обобщенная модульная схема алгоритма (имеющая два заменяемых модуля проверки и сжатия), использующего операцию инвертор. Приведены доказательства теорем о сходимости метода, решения прикладных задач управления химическим процессом и преследования маневрирующей цели перехватчиком.
Ключевые слова: интервальные методы, глобальная условная оптимизация, оптимальное управление.
APPLICATION OF THE GENERALIZED INVERSE INTERVAL METHOD OF GLOBAL CONSTRAINED OPTIMIZATION FOR OPTIMAL PROGRAM CONTROL PROBLEM
A.V. Panteleev, V.N. Panovskiy
Moscow Aviation Institute (National Research University), Moscow, Russian Federation
e-mail: avpanteleev@inbox.ru; panovskiy.v@yandex.ru
The algorithmic and program software for the generalized inverse interval method of global constrained optimization as well as its application technique are developed for searching of the optimal program control of nonlinear deterministic continuous dynamical systems. The generalized module algorithm scheme (with two changeable check and compressibility modules) using the inverter operation was developed. The convergence theorems proofs, solutions of applied control problems (chemical process control and pursuit of a maneuvering target by an interceptor) are given.
Keywords: interval methods, global constrained minimization, optimal control.
Введение. В современной математике достаточно большое внимание уделяется решению задач глобальной оптимизации [1-3], возникающих в процессе проектирования конструкций самолетов, вертолетов, космических аппаратов, когда появляется необходимость оптимизации
характерных параметров (масса, дальность полета, аэродинамические характеристики). Существующие численные методы используют разнообразные подходы, однако их применение связано с большими вычислительными затратами, излишними требованиями к постановке задачи, трудностями в достижении сходимости метода.
В качестве математического аппарата, определяющего предлагаемые в настоящей работе алгоритмы, используется интервальный анализ [4-6]. Существующие интервальные методы оптимизации можно подразделить на методы безусловной (алгоритмы Мура-Скелбоу, Ичиды-Фуджи, Дюсселя, интервальный алгоритм "имитации отжига", метод случайного интервального дробления и др. [7-9]) и условной (методы Хансена, Мура и др. [10, 11]) оптимизации. Эвристические интервальные алгоритмы были применены в работе [12] для решения ряда прикладных задач. Предлагаемый интервальный метод глобальной условной оптимизации в процессе поиска использует операцию "инвертор", имеет модульную структуру и относится к группе детерминированных методов [4]. В частных случаях реализаций модулей он совпадает с методом дихотомии прямого образа, методом отсечки и стохастической отсечки виртуальных значений, которые относятся к группе методов дробления графика [9, 13]. Показано, как общие задачи поиска оптимального программного управления нелинейными непрерывными детерминированными динамическими системами могут быть сведены к задаче интервальной условной оптимизации функций многих переменных и решены приближенно с помощью разработанных интервальных алгоритмов. Приведены примеры задачи управления химическим процессом, в которой имеется как локальный, так и глобальный экстремум, а также задачи преследования маневрирующей цели перехватчиком. На указанных примерах проиллюстрирована эффективность предложенного подхода.
Используемые понятия и обозначения. Основная идея интервального анализа — окружение вещественных чисел интервалами, а вещественных векторов — интервальными векторами (параллелотопами) [4]. Для обозначения интервала используются строчные латинские буквы, заключенные в квадратные скобки ([а], [Ь], [с],...), или представление ([а;; аи], [Ь;; Ьи], [с; си],...), где указывается нижняя и верхняя границы. Для задания параллелотопа применяется то же обозначение с полужирным начертанием букв ([а], [Ь], [с],...), или прямое произведение интервалов ([1; 3] х [2; 3] х ... х [10; 11]).
Для произвольного интервала [х] определены следующие параметры: нижняя граница [х] = 1Ь ([х]) = 8ир{£ € К. и {-то, то}|У£ € [х],
£ < С}; верхняя граница [х] = иЬ ([х]) = т£{£ € К и {-то, то} € € [х],( < £}; ширина непустого интервала ш ([х]) = [х] — [х]; средняя
[x] +
точка ограниченного и непустого интервала mid ([x]) =
x
2
Перечисленные параметры определены и для параллелотопов: нижняя и верхняя границы, а также средняя точка становятся векторами, ширина рассчитывается как максимальное значение ширины всех образующих параллелотоп компонентов.
Интервальной оболочкой множества X С М называется параллелотоп с наименьшей шириной, который содержит множество X. Если множество X берется в квадратные скобки, это означает, что рассматривается интервальная оболочка [X] этого множества.
Пусть о — некоторая бинарная операция, тогда [х] о [у] = = [{£1 о £2€ [х] , £2 £ [у]}]. Обозначим через f некоторый унарный оператор, тогда f ([х]) = [^ (£) |£ € [х]}]. Описанная бинарная операция позволяет определить арифметические операции над интервалами [4].
Множество интервалов обозначается как Ж, множество интервальных векторов - Жп. Пусть имеется некоторая функция f : Мп ^ М. Функция ^] (•) называется интервальной функцией включения для f : Жп ^ Ж, если f ([х]) = ^ (£) |£ € [х]} С ] ([х]), V [х] € Жп. Функция включения позволяет получить априорную оценку множества значений функции, даже если оно не является выпуклым или связным (если вместо переменных используются интервалы и соответствующие арифметические операции, то оценка называется оценкой прямого образа функции).
Инвертор ^, [б] , [/], е) — функция, которая по заданному
интервалу [/], функции f и параметру точности е находит в области поиска [б] множество параллелотопов Р, такое, что V [х] € Р справедливо выражение [х] С [б] ([х]) > е, [f ] ([х]) С [/], или [х] С [б] ([х]) < е, ^] ([х]) П [/] = 0. Инвертор может быть реализован с помощью алгоритма БГУГА [4].
Алгоритм SIVIA. Последовательность выполнения указанного алгоритма приведена ниже.
Шаг 1. Задать функцию f, область поиска [б], интервал [/] и параметр точности е. Создать пустое множество параллелотопов Р и множество X = {[б]}.
Шаг 2. Для каждого параллелотопа [х] € X найти оценку прямого образа ^] ([х]). Проверить выполнение следующих условий (параллелотоп после проверки удаляется из множества X):
а) если ^] ([х]) П [/] = 0, то параллелотоп [х] не добавляется в множество Р;
б) если [^ ([х]) С [/], или ^] ([х]) П [1] = 0,ш ([х]) < е, то добавить параллелотоп [х] в множество Р;
в) в противном случае — задать параллелотопы [x]' = [xi] х ... х
х mid ([xk ]); [xfe ] х ... х [xn] и [x]'' = [xi] х ... х [xfc ]; mid ([xfe ]) х
х... х [xn], где k — номер компоненты параллелотопа [x] с наибольшей шириной, и добавить их в множество X. Параллелотопы [x]' и [x]'' являются левым и правым результатами бисекции, примененной к параллелотопу [x].
Шаг 3. Если X = 0, то перейти к шагу 2, в противном случае — к шагу 4.
Шаг 4. Вернуть множество параллелотопов P.
Обобщенный инверсный интервальный метод глобальной условной оптимизации. Постановка задачи интервальной е-мини-мизации и стратегия поиска решения. Пусть имеются параллелотоп [s], задающий множество допустимых решений, целевая функция f : Rn ^ R, малое число е > 0. Требуется найти параллелотоп [p]*, такой, что
[p]* С [s] ,ш ([p]*) < е, ? [x] С [s] ,ш ([x]) > е : [f] ([x]) < [f] ([p]*).
(1)
В основе метода лежит идея оценки значения функции на области поиска в виде целевого интервала. На каждой итерации алгоритма целевой интервал делится на два новых интервала [13]. Далее выполняется операция проверки того, в каком из полученных интервалов содержится значение минимума функции. Выбранный интервал вновь обозначается как целевой интервал и начинается следующая итерация. Поскольку первичная оценка значения функции на области поиска может быть неточной [4], целесообразно уточнить ее с помощью применения операции сжатия [13]. Операции проверки и сжатия выделены в алгоритме как отдельные заменяемые модули вследствие неединственности возможного описания этих процедур. Схема алгоритма представлена на рис. 1. Предлагаемый алгоритм называется обобщенным инверсным интервальным, так как при поиске наилучшего параллелотопа используются заменяемые модули сжатия и проверки, а также операция инвертор.
В качестве заменяемых модулей в алгоритме использованы следующие операторы.
1. Оператор проверки, на вход которого подается множество параллелотопов X, проверяемый интервал [/] и параметр ширины w. Результатом его работы является статус проверки "удачный". Это свидетельствует о том, что в множестве параллелотопов X есть параллелотоп [p], для которого выполнено условие ш ([p]) > w, [f] ([p]) С [/], или ш ([p]) < w, [f] ([p]) П [/] = 0). Статус проверки "неудачный" в противном случае. В алгоритме предлагается использовать:
— проверку до первого подходящего элемента First True operator (FT-operator);
Рис. 1. Схема обобщенного инверсного интервального алгоритма
— проверку до первого подходящего элемента, множество X заменяется множеством, выработанным инвертором First True with Renewal operator (FTR-operator).
2. Оператор сжатия, на вход которого подается параллелотоп [s], задающий область поиска. Результат его работы — сжатие оценки прямого образа функции на области поиска. В алгоритме предлагается применять:
— сжатие на основе разбиения области поиска, которое реализуется путем представления каждой компоненты параллелотопа, описывающего область поиска, в виде объединения непересекающихся интервалов Search Area Split operator (SAS-operator);
— сжатие на основе значений функций в случайно сгенерированных точках, которые являются интервальными векторами нулевой ширины Random Point Sample operator (RPS-operator).
Алгоритм обобщенного инверсного интервального метода. Последовательность выполнения указанного алгоритма приведена ниже.
Шаг 1. Задать область поиска [s] ; параметр точности, характеризующий размер параллелотопа в множестве P, которое будет вырабатываться инвертором б > 0; параметр точности для остановки алгоритма Z > 0. Выбрать оператор сжатия SAS с параметрами w и A или RPS с параметром A, а также оператор проверки FT или FTR с параметром w. Задать множество X = {[s]}.
Шаг 2. С помощью функции включения [f ] найти оценку прямого образа [Y] функции f на области поиска [s].
Y
, опреде-
Шаг 3. Заменить полученную оценку [Т] интервалом ленным с помощью выбранного оператора сжатия.
Шаг 4. С помощью операции бисекции разделить интервал [Т] =
= [y ; yu] на два интервала Y ] =
yi ;
Vi + Уи 2
и [Yu] =
Vi + Vu 2
-; Vu
Шаг 5. Проверить интервал [Y/] выбранным оператором проверки.
Шаг 6. Если статус проверки "удачный", то перейти к шагу 7, если нет — к шагу 8.
Шаг 7. Если ш ([Тг]) < С, то найти Р = (/, [б] , [Тг] ,е) и перейти к шагу 9. В противном случае принять интервал [Т] равным [Т;] и перейти к шагу 4.
Шаг 8. Положить интервал [Т] равным [Ти]. Если ш ([Т]) < £, то найти Р = (/, [б] , [Т], е) и перейти к шагу 9. В противном случае — перейти к шагу 4.
Шаг 9. Выбор параллелотопа.
Шаг 9.1. Принять г = 1, множество Ре = Р = {[р^}.
Шаг 9.2. Если ш ([р]^) < е, где [р]^ € Ре, то перейти к шагу 9.4, иначе — к шагу 9.3.
Шаг 9.3. Добавить в множество Ре параллелотопы [р^ х ... х
([pfeL); [pfe]» x...x [pn]»и [pi]» x...x [pfe]»;mid ([pk]»)
...
pn]i? где k — номер компоненты параллелотопа [p]^ с наибольшей шириной, удалить [p]^ из множества Pe.
Шаг 9.4. Увеличить i на единицу. Если i < |Pe|, где | • | — мощность множества, то перейти к шагу 9.2. В противном случае перейти к шагу 9.5.
Шаг 9.5. Для каждого [p]^ £ Pe найти [f] ([p]^) и выбрать [p]* =
= Arg min [f] ([p]i).
[p]»ep£-
Алгоритмы реализации операторов проверки. На вход всех операторов проверки подается множество параллелотопов X = {[х]^}, проверяемый интервал [/] и параметр ширины w.
Алгоритм FT-оператора. Последовательность выполнения указанного алгоритма приведена ниже.
Шаг 1. Принять i = 1.
Шаг 2. Выполнить следующие действия:
1) если [f] ([x]J П [/] = 0, то продолжить работу;
2) если [f]([x]i) С [/], или [f] ([x]J П [/] = 0, ш ([x]J < w, то проверка завершена "удачно";
3) в случае невыполнения действий 1 и 2 найти параллелотопы
[x]' = [xi]4 х ... х mid([xk]t); [xk ]г х ... х [xn]4 и [x]'' = [x^ х
х ... х [xk mid ([xk х ... х [xn]^, где k — номер компоненты
параллелотопа [x]^ с наибольшей шириной, и добавить их в множество X.
Шаг 3. Если i < |X|, то увеличить значение i на единицу и перейти к шагу 2, а если i = |X| , то проверка завершена "неудачно".
Алгоритм FTR-оператора. Последовательность выполнения указанного алгоритма приведена ниже.
Шаг 1. Принять i = 1. Задать множество X = X = {[x^j.
Шаг 2. Выполнить следующие действия (сам параллелотоп [x]» после проверки удаляется из множества X):
1) если [f] ([x]») П [/] = 0, то продолжить работу;
2) если [f]([x]i) С [/], или [f] ([x]J П [/] = 0, ш ([x]J < w, то проверка завершена "удачно";
3) в случае невыполнения действий 1 и 2 найти параллелотопы
[x]' = [xi]4 х ... х mid([xk]i);[xk]г х ... х [xn]» и [x]'' = [xi]» х
X
X
[xk ]»; mid ([xk ]») х ... х [xn]», где k — номер компоненты
параллелотопа [х]» с наибольшей шириной, и добавить их в множество X.
Шаг 3. Если г < |Х|, то увеличить значение г на единицу и перейти к шагу 2, а если г = |Х| — перейти к шагу 4.
Шаг 4. Заменить множество X множеством X, сообщить, что проверка завершена "неудачно".
Алгоритмы реализации операторов сжатия. На вход всех операторов сжатия подается параллелотоп [б].
Алгоритм SAS-оператора (для работы необходим параметр ширины ы и число попыток А). Последовательность выполнения указанного алгоритма приведена ниже.
Шаг 1. Найти
Y
= [f] ([s]).
Шаг 2. Для каждой компоненты вектора найти такое наимень-
ш (М)
шее целое число п, что
< w.
Шаг 3. Каждый интервал представить в виде объединения непересекающихся интервалов У . Таким образом, создается разби-
j=i
N
ение начального параллелотопа [s] = [J [s]k, где N = 1 [ n.
Шаг 4. Найти
Т
U
и
k=1
min
([s]k)
»=i
fe=jl,...,JA
T
mm
([s]k)
Y
результирующий сжатый интер-
й=л,...,.м
вал, где ^..., ^ — случайно выбранные натуральные числа из интервала [1; N ].
Алгоритм RPS-оператора (для работы необходимо число точек А). Последовательность выполнения указанного алгоритма приведена ниже.
Шаг 1. Генерировать А случайных точек € [б], г = 1,..., А.
Шаг 2. Найти сжатый интервал.
Y
[f] ([S]), Ш^А f &)
результирующий
Оценка качества решения. Теорема 1. Обобщенный инверсный интервальный метод находит решение задачи интервальной е-минимизации.
М Поскольку на шаге выбора (шаг 9) параллелотоп выбирается из множества Pe, его ширина никогда не превысит выбранного значения е (вследствие способа построения множества Pe). На каждом цикле алгоритм изначально с помощью инвертора вырабатывает множество параллелотопов P, в котором для любого параллелотопа [p] £ P справедливы соотношения [p] С [s], ш ([p]) > е и [f ] ([p]) С [Y;], или [p] С [s], ш ([p]) < е и [f ] ([p]) П [Y;] = 0. Поэтому переход к анализу [Yu] свидетельствует о том, что множество P пусто, т.е. V [x] С [s], ш ([x]) < е выполнено [f] ([p]) С [Yu]. Таким образом, в ходе каждого цикла алгоритма однозначно определяется интервал [у], такой, что $ [x] С [s], для которого выполнено условие [f] ([x]) < [у]. Это исключает возможность существования параллелотопа [x] С [s], ш ([x]) > е, для которого [f] ([x]) < [f] ([p]*), т.е. справедливо условие (1). ►
Теорема 2. Пусть x* = Arg min f (x), а [p]* — решение задачи
x€[s\
интервальной е-минимизации обобщенным инверсным интервальным методом, тогда f (x*) £ [f] ([p]*).
< Необходимо показать, что для параллелотопа [p]* одновременно выполняются неравенства f (x*) > [f] ([p]*), f (x*) < [f] ([p]*). Докажем выполнение первого неравенства. Пусть f (x*) < [f] ([p]*), тогда 3 [p] — это такое решение интервальной задачи интервальной е-минимизации, что f (x*) £ [f ] ([p]). Откуда следует [f] ([p]) < [f] ([p]*), что невозможно, так как [p]* является решением задачи интервальной
е-минимизации и при этом [p]* = Arg min [f] ([p]J. Следовательно,
[p]»ep£-
первое неравенство выполнено. Докажем выполнение второго неравенства. Пусть f (x*) > [f ] ([p]*), тогда [f ] ([p]*) = 0, что невозможно, поскольку [p]* = 0. Таким образом, второе неравенство выполнено. Два неравенства выполнены, теорема 2 доказана. ►
Следствие теоремы 2. В параллелотопе [p]* есть точка x, такая, что f (x) - f (x*) < ш ([f]([p]*)).
Оптимальное интервальное программное управление непрерывными системами. Постановка задачи. Поведение модели объекта управления описывается дифференциальным уравнением
x (t) = f (t,x (t) ,u (t)), (2)
где t — непрерывное время; t £ T = [to,ti] — заданный промежуток времени функционирования системы; x £ Rn — вектор состояния системы; u £ [u] С Rq — вектор управления; [u] — множество допустимых значений управления, представляющее собой параллело-
топ; / (¿, х, и) = (/1 (¿, х, и),..., /п (¿, х, и))т — непрерывная вектор-функция.
Начальное состояние системы (2) задано и равно
х (¿о) = хо. (3)
Конечное состояние х (¿1) должно удовлетворять условиям
Г (х (¿1)) = 0, г = 1,...,/, (4)
где 0 < / < п; Г» (х) — непрерывно дифференцируемые функции; |дг (х)/дх1,..., д (х)/<9хп|, г = 1,..., / — система векторов, линейно независимая для Ух € Кп.
При управлении используется информация только о непрерывном времени ¿, т.е. применяется программное управление. Множество допустимых процессов О (¿0, х0) определяется как множество пар d = (х (•), и (•)), включающих в себя траекторию х (•) и кусочно-непрерывное управление и (•), где У € Т и (¿) € [и], удовлетворяющих уравнению состояния (2), начальному (3) и конечному (4) условиям.
На множестве допустимых процессов определен функционал качества управления
¿1
I = у /0 (¿, х(£), и(*)) dí + ^ (х (¿1)), (5)
¿0
где /0 (¿, х, и), ^ (х) — заданные непрерывные функции.
Требуется найти пару d* = (х* (•), и* (•)) € О(£0, х0), при которой достигается минимальное значение функционала (5) на множестве допустимых процессов.
Основным элементом предлагаемого подхода к решению поставленной задачи является ее последовательное сведение к соответствующей задаче нелинейного программирования, а затем — к задаче интервальной е-минимизации и решение последней с помощью разработанного обобщенного инверсного интервального метода глобальной условной оптимизации.
Для более простого вычисления интегрального члена в критерии (5) к системе (2) добавляется уравнение хсп+1 (¿) = /0 (¿, х (¿), и (¿)) с начальным условием хг+1 (¿0) = 0, тогда значение функционала качества (5) определяется по выражению I = хг+1 (¿1) + ^ (¿1, х (¿1)).
Для сведения поставленной задачи с ограничениями (4) к задаче со свободным правым концом траектории к терминальному члену функционала добавляются либо классические штрафные функции, характеризующие степень выполнения условий (4), либо их интервальный
аналог вида
£
i=1
: Яг([Г (X (¿1)); Г (х (¿1))], [—&; &]), (6)
где Л«([а], [5]) = [ш1и(Д0оо([а], [5]),Л°,([а], [5])),шах(Д0оо([а], [5]), ([а], [5])) — метрика Хаусдорфа; ^ ([а], [5]) = г| [а] С [5] +
+ [-г; г] ,г > о| — мера близости интервалов; Яг — параметры
штрафных функций; > 0 — допустимая погрешность выполнения конечного условия, задаваемая пользователем.
Искомое оптимальное программное управление предлагается искать в двух различных классах функций: кусочно-постоянных и кусочно-линейных. Поскольку выполняется поиск наилучшего управления в каждом классе, не совпадающем с классом кусочно-непрерывных функций, найденное управление является субоптимальным, рассматриваемым в качестве кандидата в решение задачи. При решении задачи применяется интервальный метод, поэтому каждому описанному управлению сопоставляется интервальный аналог, задаваемый интервальным вектором.
Для кусочно-постоянного интервального управления (рис. 2, а) необходимо задать значения функции и (¿) в N моментах време-
¿1 — ¿0 . . П дт .
-г, г = 0,..., N — 1, т.е. управлению можно
ни T = to +
N
однозначно поставить в соответствие интервальный вектор [и] =
= [и1 (то)] х ... х [ид (то)] х ... х [и1 (т^-1)] X ... х [ид (т^-1)]. Тогда
Кто)]
/ \
/ \
/
\ / /
/
[u(tn-1 )] u(t)
о
¿0
tx t
to
tx t
to
tx t
Рис. 2. Замена кусочно-непрерывного управления интервальным кусочно-постоянным (а) и кусочно-линейным (б) управлениями
интервальное управление находится по формуле [u] (t) G [u (г»)] = = [u (т») ; u (т»)], гг < t < Ti+i, i = 0,... ,N - 1.
Для кусочно-линейного интервального управления (рис. 2, б) необходимо задать дополнительное значение управления u (t) на последнем временном интервале. Поэтому интервальный вектор, который ставится в соответствие управлению, имеет вид [u] =
= [ui (то)] х ... х [uq (то)] х ... х [ui (tn)] X ... х [uq (tn)]. Интерваль-
4-v-' 4-V-'
[u(t0)] MTN )]
ное управление определяется по формуле [u] (t) G —-[u (т)] +
Ti+1 — Ti
t — T
+-— [u (Ti+i)], Ti < t < Ti+1, i = 0,..., N — 1.
Ti+1 — Ti
Интервальной траекторией называется решение уравнения (2) с начальным условием (3), соответствующее заданному интервальному управлению. Если система (2), описывающая поведение модели объекта управления, не жесткая, то для ее численного интегрирования применяют явные методы (методы Эйлера, Эйлера-Коши, методы Рунге-Кутты различных порядков и др.). Если система (2) жесткая, то используют неявные методы численного интегрирования (неявный метод Эйлера, метод трапеций, методы Адамса-Моултона, Гира и др.). Вычисление правых частей системы (2) проводится по правилам интервальной арифметики [4].
Отметим, что описанные процедуры могут быть использованы при других способах параметризации управления, например, кусочно-полиномиальном, или при разложении по системам ортонормирован-ных базисных функций.
Решение прикладных задач. На основе приведенного алгоритма написана программа поиска субоптимального интервального программного управления. Среда разработки Microsoft Visual Studio, язык программирования C#.
Для исследования эффективности разработанного метода и программного обеспечения использованы принципы, изложенные в работе [14]: полученное приближенное решение необходимо сравнить с точным решением (если известно), с результатами применения принципа максимума, а также с результатами применения других методов.
Задача управления химическим процессом. Рассмотрим решение задачи оптимального управления химическим процессом [15, 16]. Модель объекта управления описывается системой дифференциальных уравнений
¿1 (t) = - (2 + u) (¿1 + 0,25) + (¿2 + 0,5) exp ' 25x1
¿1 + 2/
25жП (7)
¿2 (t) = 0,5 — ¿2 — (¿2 + 0,5) exp ^ ^
¿3 (t) = ¿2 + ¿2 + 0,1u2,
где ¿ь ¿2 — отклонения от безразмерной установившейся температуры и безразмерной установившейся концентрации; ¿3 — величина, определяющая значение критерия качества при t = t1, t G [0; 0,78]; начальное состояние ¿ (0) = (0,09; 0,09; 0)т, ограничение на управление [u] = [—10; 10]; критерий качества управления I = ¿3 (0,78) ^ min.
Результаты решения приведенной задачи методом динамического программирования (I = 0,13336) в виде зависимостей изменений компонент векторов состояния и управления представлены на рис. 3, а, б [15, 16]. На этом же рисунке приведены результаты решения задачи, полученные с помощью принципа максимума Понтрягина (локальный минимум I = 0,244425 и глобальный минимум I = 0,133094), алгоритма Нелдера-Мида (I = 0,1333666) и генетических алгоритмов (I = 0,13341).
Параметры обобщенного инверсного интервального метода: в качестве метода дискретизации выбран явный метод Рунге-Кутты четвертого порядка; области поиска [s] = [—10; 10] х ... х [—10; 10] для
4-v-'
N
кусочно-постоянного управления и [s] = [—10; 10] х ... х [—10; 10]
4-v-'
N+1
для кусочно-линейного управления, е = Z = 0,00001; N = 10. Результаты работы алгоритма и исследования эффективности алгоритма представлены на рис. 3, в, г ив табл. 1.
Сравнивая значения критерия оптимальности (см. табл. 1) со значениями, приведенными в работах [15, 16], можно сделать следующий вывод: разработанный алгоритм позволяет найти решение, близкое к глобальному минимуму, не застревая в окрестности локального минимума, существующего в рассматриваемой задаче. Его эффективность при заданной точности сравнима с эффективностью алгоритмов применения принципа максимума, итерационного динамического программирования, метода Нелдера-Мида и эволюционных алгоритмов. Кусочно-линейное управление предпочтительнее по функционалу, чем кусочно-постоянное. Полученные результаты иллюстрируют применимость теорем 1 и 2, поскольку находится решение задачи е-минимизации и гарантируется принадлежность глобального оптимума к получаемому интервалу значений критерия оптимальности.
Рис. 3. Зависимости изменения компонент векторов состояния Х\ (1) и х2 (2) (а) и управления и (б), найденные в работах [15, 16] и с помощью разработанного алгоритма в случае кусочно-постоянного (в) и кусочно-линейного (г) управлений
Задача преследования. Рассмотрим решение задачи преследования в трехмерном пространстве [17]. Поведение модели объекта управления описывается системой дифференциальных уравнений
X (¿) = Ух (¿), К (*) = их (¿), х* (¿) = УХ (¿), Ух* (*) = иХ (¿),
у (¿) = Уу (¿), Уу (¿) = иу (¿), у * (¿) = уу (¿), У* (¿) = иу (¿), (8)
г (¿) = У2 (¿), Уг (¿) = и (¿), г* (¿) = У* (¿), У* (¿) = и'г (¿),
где г = (х,у, г)т, г* = (х*, у*, г*)т — радиус-векторы, описывающие положение перехватчика и цели; У = (УХ, Уу, Уг)т, У* = (УХ*, У*, У*)т — ве1сторы скоростей перехватчика и цели; и = (иХ, иу, и^)т, и* = = (иХ, иу, и^) — ускорения перехватчика (управления) и цели.
Результат работы алгоритма
Операторы проверки и сжатия Функционал качества I
Кусочно-постоянное управление
FT (w = 0,001), SAS(w = 2, A = 25) [0,139861; 0,140337]
FTR (w = 0,001), SAS (w = 2, A = 25) [0,139549; 0,140376]
FT (w = 0,001), RPS (A = 100) [0,139831; 0,140791]
FTR (w = 0,001), RPS (A = 100) [0,139642; 0,140509]
Кусочно-линейное управление
FT (w = 0,001), SAS(w = 2, A = 25) [0,133087; 0,133121]
FTR (w = 0,001), SAS (w = 2, A = 25) [0,133092; 0,133101]
FT (w = 0,001), RPS (A = 100) [0,133095; 0,133118]
FTR (w = 0,001), RPS (A = 100) [0,133001; 0,133195]
Предполагается, что цель движется с постоянным ускорением П = (1; —2; 0,1)т м/с2; начальное положение цели Г = = (500; —600; 500)т м; начальная скорость цели (100; 100; 10)т м/с; начальное положение перехватчика г = (0,0,0)т м; начальная скорость перехватчика V = (150; 40; 5)т м/с. Время осуществления перехвата (координаты цели и перехватчика равны): ¿1 = 50 с. Функционал качества управления (5) задается в виде
I = (х (¿1) — X (¿1))2 + (у (¿1) — У (¿1))2 + (г (¿1) — / (¿1))2 .
Для сглаживания искомого интервального управления в функционал
м
добавим "штрафное" слагаемое Д! = 0,001 ^ (и^ + п^ + п^), где
пх,, = их (т,), = пу (т,), = п (т,), М = N — 1 для кусочно-постоянного управления, М = N для кусочно-линейного управления. В результате решается задача минимизации функционала I + Д!.
Параметры обобщенного инверсного интервального метода: в качестве метода дискретизации выбраны явная и неявная схемы Эйлера; область поиска [б] = [—2; 2] х ... х [—2; 2] для кусочно-постоянного
4-V-'
3-м
управления и [б] = [—2; 2] х ... х [—2; 2] для кусочно-линейного упра-
4-V-'
3-(м+1)
вления; е = 0,001; £ = 0,1; N = 2. Зависимость движения цели и перехватчика, полученная в работе [17], представлена на рис.4,а. Результаты работы алгоритма приведены на рис. 4, б, в ив табл. 2.
Сравнительный анализ управления и траекторий, найденных интервальным методом и методом пропорциональной навигации, свиде-
Рис. 4. Зависимость движения цели и перехватчика, полученная с помощью алгоритма пропорциональной навигации [17] (а), зависимости движения цели и перехватчика, законов изменения управления, найденные с помощью разработанного алгоритма при кусочно-постоянном (б) и кусочно-линейном (в) управлениях
тельствует о близости полученных результатов. Максимальное суммарное отклонение цели и перехватчика А = Ax + Ay + Az = 0,09 + + 0,14 + 0,07 = 0,3 для кусочно-постоянного управления и А = Ax + + Ay + Az = 0,06 + 0,08 + 0,14 = 0,28 для кусочно-линейного управления, что соответствует требованиям, предъявляемым к точности решения задачи преследования. Здесь Ay = max |y (ti) — y1 (ti)|, Az = max |z (ti) — z1 (ti )| — погрешности выполнения терминальных условий по каждой координате. Кусочно-линейное управление предпочтительнее по функционалу, чем кусочно-постоянное (см. табл.2), поскольку приводит к меньшим значениям ширины интервала неопределенностей значений критерия оптимальности. Отметим, что
Результат работы алгоритма
Операторы проверки и сжатия
Функционал качества I для схем Эйлера
Явная
Неявная
Кусочно-постоянное управление
FT (w = 0,01), SAS(w = 2, A = 25) [0; 0,0943] [0; 0,0912]
FTR (w = 0,01), SAS (w = 2, A = 25) [0; 0,0997] [0; 0,0973]
FT (w = 0,01), RPS (A = 100) [0; 0,0995] [0; 0,0996]
FTR (w = 0,01), RPS (A = 100) [0; 0,0990] [0; 0,0988]
Кусочно-линейное управление
FT (w = 0, 01), SAS(w = 2, A = 25) [0; 0,0921] [0; 0,0893]
FTR (w = 0,01), SAS (w = 2, A = 25) [0; 0,0909] [0; 0,0910]
FT (w = 0,01), RPS (A = 100) [0; 0,0887] [0; 0,0884]
FTR (w = 0,01), RPS (A = 100) [0; 0,0972] [0; 0,1087]
каждому управлению соответствуют интервалы, содержащие глобальный минимум.
Заключение. Предложен новый детерминированный обобщенный инверсный интервальный метод условной глобальной оптимизации функций многих переменных, сформулированы и доказаны теоремы о свойствах метода. Разработан алгоритм для решения задач поиска субоптимального программного управления непрерывными нелинейными детерминированными динамическими системами в классе кусочно-постоянных и кусочно-линейных управлений. Решены две прикладные задачи: управления химическим процессом и преследования цели перехватчиком. К преимуществам предложенного алгоритма относится возможность гарантии принадлежности глобального оптимума получаемому интервалу значений критерия оптимальности.
ЛИТЕРАТУРА
1. Пантелеев А.В., Летова Т.А. Методы оптимизации. Практический курс. М.: Логос, 2011. 424 с.
2. Пантелеев А.В., Метлицкая Д.В., Алешина Е.А. Методы глобальной оптимизации. Метаэвристические стратегии и алгоритмы. М.: Вузовская книга, 2013. 248 с.
3. Пантелеев А.В. Применение эволюционных методов глобальной оптимизации в задачах оптимального управления детерминированными системами. М.: Изд-во МАИ, 2013. 160 с.
4. Jaulin L., Kieffer M., Didrit O., Walter E. Applied interval analysis. London: Springer-Verlag, 2001. 384 p.
5. Moore R.E. Interval analysis. Englewood Cliffs, Prentice Hall, 1966. 145 p.
6. Шарый С.П. Конечномерный интервальный анализ. Новосибирск: XYZ, 2010. 606 с.
7. Ratschek H., Rokne J.New computer methods for global optimization. Chichester: Horwood, 1988. 229 p.
8. Shary S.P. Randomized algorithms in interval global optimization // Numerical Analysis and Applications. 2008. Vol. 1. P. 376-389.
9. Shary S.P. A surprising approach in interval global optimization // Reliable Computing. 2001. P. 497-505.
10. Hansen E. Global optimization using interval analysis. New York: Marcel Dekker, 2004. 515 p.
11. Moore R.E. Methods and applications of interval analysis. Philadelphia: SIAM, 1979. 190 p.
12. Пановский В.Н. Прикладное применение интервального метода взрывов // Труды МАИ. 2014. № 73.
URL: http://www.mai.ru/science/trudy/published.php?ID=48451
13. Пановский В.Н. Применение аппарата интервального анализа для поиска глобального экстремума функций // Труды МАИ. 2012. № 51.
URL: http://www.mai.ru/science/trudy/published.php?ID=28953
14. Федоренко Р.П. Приближенное решение задач оптимального управления. М.: Наука, 1978. 488 с.
15. Пантелеев А.В., Метлицкая Д.В. Применение генетических алгоритмов с бинарным и вещественным кодированием для приближенного синтеза субоптимального управления детерминированными системами // Автоматика и телемеханика. 2011. № 11. С. 117-129.
16. Luus R. Iterative dynamic programming. New York: Chapman & Hall, 2000. 331 p.
17. Tewari A. Advanced control of aircraft, spacecraft and rockets. New York: John Wiley & Sons, 2011.448 p.
REFERENCES
[1] Panteleev A.V., Letova T.A. Metody optimizatsii. Prakticheskiy kurs [Methods of Optimization. Practical Course]. Moscow, Logos, 2011. 424 p.
[2] Panteleev A.V., Metlitskaya D.V., Aleshina E.A. Metody global'noy optimizatsii. Metaevristicheskie strategii i algoritmy [Methods of Global Optimization. Metaheuristic Strategies and Algorithms], Moscow, Vuzovskaya kniga Publ., 2013. 248 p.
[3] Panteleev A.V. Primenenie evolyutsionnykh metodov global'noy optimizatsii v zadachakh optimal'nogo upravleniya determinirovannymi sistemami [Application of Evolutionary Methods of Global Optimization in Problems of Finding the Optimal Control of Determined Systems], Moscow, MAI Publ., 2013. 160 p.
[4] Jaulin L., Kieffer M., Didrit O., Walter E. Applied interval analysis. London, Springer-Verlag, 2001. 384 p.
[5] Moore R.E. Interval analysis. Englewood Cliffs, Prentice Hall, 1966. 145 p.
[6] Sharyy S.P. Konechnomernyy interval'nyy analiz [Finite-dimensional Interval Analysis], Novosibirsk, XYZ Publ., 2010. 606 p.
[7] Ratschek H., Rokne J. New computer methods for global optimization. Chichester, Horwood, 1988. 229 p.
[8] Shary S.P. Randomized algorithms in interval global optimization. Numerical Analysis and Applications, 2008, vol. 1, pp. 376-389.
[9] Shary S.P. A surprising approach in interval global optimization. Reliable Computing, 2001, pp. 497-505.
[10] Hansen E. Global optimization using interval analysis. N.Y., Marcel Dekker, 2004. 515 p.
[11] Moore R.E. Methods and applications of interval analysis. Philadelphia, SIAM, 1979. 190 p.
[12] Panovskiy V.N. Application of the interval explosion method. Electronnyy zhurnal "Trudy MAI", 2014, iss. 73. Available at: http://www.mai.ru/science/trudy/published.php?ID=48451
[13] Panovskiy V.N. Application of the interval analysis for the search of the global extremum of functions. Electronnyy zhurnal "Trudy MAI", 2012, iss. 51. Available at: http://www.mai.ru/science/trudy/published.php?ID=28953
[14] Fedorenko R.P. Priblizhennoe reshenie zadach optimal'nogo upravleniya [Approximate Solution of the Optimal Control Problems]. Moscow, Nauka Publ., 1978. 488 p.
[15] Panteleev A.V., Metlitskaya D.V. An application of genetic algorithms with binary and real coding for approximate synthesis of suboptimal control in deterministic systems. Avtomat. i Telemekh, 2011, no. 11, pp. 117-129 (in Russ.).
[16] Luus R. Iterative dynamic programming. N.Y.: Chapman & Hall, 2000. 331 p.
[17] Tewari A. Advanced control of aircraft, spacecraft and rockets. N.Y.: John Wiley & Sons, 2011. 448 p.
Статья поступила в редакцию 18.12.2014
Пантелеев Андрей Владимирович — д-р физ.-мат. наук, профессор, заведующий кафедрой "Математическая кибернетика" МАИ (НИУ).
Московский авиационный институт (национальный исследовательский университет), Российская Федерация, 125993, Москва, Волоколамское ш., д. 4.
Panteleev A.V. — Dr. Sci. (Phys.-Math.), Professor, Head of Mathematical Cybernetics department, Moscow Aviation Institute (National Research University). Moscow Aviation Institute (National Research University), Volokolamskoe shosse 4, Moscow, 125993 Russian Federation.
Пановский Валентин Николаевич — аспирант кафедры "Математическая кибернетика" МАИ (НИУ).
Московский авиационный институт (национальный исследовательский университет), Российская Федерация, 125993, Москва, Волоколамское ш., д. 4.
Panovskiy V.N. — post-graduate student of Mathematical Cybernetics department, Moscow Aviation Institute (National Research University).
Moscow Aviation Institute (National Research University), Volokolamskoe shosse 4, Moscow, 125993, Russian Federation.
Просьба ссылаться на эту статью следующим образом:
Пантелеев А.В., Пановский В.Н. Применение обобщенного инверсного интервального метода глобальной условной оптимизации в задаче поиска оптимального программного управления // Вестник МГТУ им. Н.Э. Баумана. Сер. Приборостроение. 2016. № 1. C. 33-50. DOI: 10.18698/0236-3933-2016-1-33-50
Please cite this article in English as:
Panteleev A.V., Panovskiy V.N. Application of the generalized inverse interval method of global constrained optimization for optimal program control problem. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Priborostr. [Herald of the Bauman Moscow State Tech. Univ., Instrum. Eng.], 2016, no. 1, pp. 33-50. DOI: 10.18698/0236-3933-2016-1-33-50
