CC BY
20
В классической теории управления, в основе которой лежит тот же принцип априорного выбора полюсов, основным критерием для подобного выбора является устойчивость движения управляемой системы. Методам, основанным на принципе максимума, внутренне присущи высокие требования к разрядности вычислений («проклятие разрядности»). Поэтому в предлагаемом подходе аналогичный критерий должен основываться на понятии обусловленности вычислений.
ПЕРЕЧЕНЬ ССЫЛОК
1. Андриевский Б. Р., Фрадков А. Л. Элементы математического моделирования в программных средах MAT-LAB и Scilab. - СПб.: Наука, 2001.
2. Красовский А. А. Науковедение и состояние теории процессов управления // АиТ. - 2000. - № 4. - С. 3-19.
3. Габасов Р., Дмитрук Н. М., Кириллова Ф. М. Оптимизация многомерных систем управления с параллелепипед-ными ограничениями. // АиТ. - 2002. - № 3. - С. 3-26.
4. СейджЭ.П, Уайт Ч. С., III. Оптимальное управление системами: Пер. с англ. - М.: Радио и связь, 1982. -392 с.
5. Квакернаак X., Сиван Р. Линейные оптимальные системы управления. - М.: Наука, 1977. - 652 с.
6. Икрамов X. Д. Численное решение матричных уравнений. - М.: Наука, 1984. - 192 с.
7. Хорн Р., Джонсон Ч. Матричный анализ. - М.: Мир, 1989. - 656 с.
8. Хайрер Э, Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. - М.: Мир, 1990. -512 с.
9. Токарев В. В. О знаках импульсов в задачах оптимального управления с закрепленными концами траектории. // АиТ. - 2001. - № 8. - С. 46-55.
10. Дуайер Т. А. V, III. Точное нелинейное управление быстрыми вращениями КЛА посредством внутренней передачи количества движения // Аэрокосмическая техника. - 1987. - № 3. - С. 151-159.
11. Синха С. К., Шривастава С. К. Оптимальный закон наведения многоступенчатой ракеты-носителя на трехмерной траектории, использующий выражения в явном виде // Аэрокосмическая техника. - 1991. - № 3. -С. 74-87.
12. Чендлер Д. С., Смит И. Е. Схема итеративного управления и ее использование для различных аппаратов и космических операций // Астронавтика и ракетоди-намика. - 1967. - № 8. - С. 1-24.
Надшшла 15.07.06
Розглядаетъся побудова анал1тичного розв'язку крайо-eoï задач1 для л1ншних стацюнарних систем без обме-женъ на керуюч1 функци. В основу запропонованого nid-ходу покладено, як i в класичнш теори керування, за-дання рoзnoдiлу коретв керoванoï системи. Розв'язок за-дачi подано у виглядi утверсалъного (базового) алгоритму.
It is consider the analytical decision of edge problem for linear stationary systems without restriction of control functions. For basis of the proposed approach it was taken, as in the classical control theory, the forming of control system roots distribution. The decision it is produced in the form of universal (base) algorithm.
УДК 62 - 506
А. К. Новиков
СИНТЕЗ ОПТИМАЛЬНЫХ УПРАВЛЕНИЙ НА ОСНОВЕ АНАЛИТИЧЕСКОГО РЕШЕНИЯ ЛИНЕЙНОЙ СТАЦИОНАРНОЙ ДВУХТОЧЕЧНОЙ КРАЕВОЙ ЗАДАЧИ
Часть 2. Примеры применения базового алгоритма
На конкретных примерах рассматриваются особенности применения базового алгоритма, построенного в первой части статьи. Для систем с ограничениями на управляющие функции используется принцип поэтапного планирования метода динамического программирования.
ВВЕДЕНИЕ
Издавна в научных кругах принято считать хорошей традицией проверку любого нового метода или алгоритма на известных примерах. Не отступая от этой традиции, автор включил в данный раздел несколько известных из литературы примеров, решенных ранее
© Новиков А. К., 2006
с использованием других подходов. Эти примеры охватывают круг задач терминального управления, оптимального регулирования и наблюдения, управления с эталонной моделью, управления с ограничениями в виде неравенств, наложенными на управляющие функции.
Остальные примеры, взятые из близкой автору области управления в космосе, иллюстрируют возможности применения базового алгоритма к решению задач регулирования и стабилизации, наведения, а также обеспечения плавной стыковки участков траектории в случае, когда общая задача управления разбивается
на несколько элементарных подзадач. Внимания заслуживает также пример использования базового алгоритма для решения такой нетрадиционной для теории оптимального управления задачи, как высокоточное численное интегрирование уравнений небесной механики.
Автор надеется, что приведенные примеры помогут читателям найти свои подходы к решению интересующих их задач в других областях применения.
ПРИМЕРЫ
Пример 1. Рассматривается задача синтеза алгоритмов оптимального регулятора с обратной связью по выходной переменной у = Сх. Как известно, такой регулятор является комбинацией наблюдателя, в котором восстанавливается состояние системы, и закона управления, представляющего собой линейную функцию этого восстановленного состояния [2].
Решение задачи оптимального регулирования для системы X = Гх + Ои и критерия качества
tf
= 1 хТ(Ь{)Ргх(Ьр + 11(xTQx + и Ки)ёЬ
Ьо
имеет вид закона управления с обратной связью и(Ь) = = -К(Ь)х(Ь), где К( Ь) = ВГХОТР (Ь), Р(Ь) - решение уравнения Риккати [2]
• Т -1 Т
-Р(Ь) = ¥1 Р(Ь) + Р(Ь)Г + Q - Р( Ь)ВК Б'Р( Ь),
Р (Ь/) = Рf, Ь < Ь^
В явном виде это уравнение решается в исключительно редких случаях, применение численных методов также является не вполне корректным, т. к. процесс численного интегрирования данного уравнения в прямом направлении является неустойчивым. Даже нахождение его установившегося решения наталкивается на серьезные трудности, поскольку эта задача тесно связана с решением уже упоминавшейся проблемы собственных значений.
Применение базового алгоритма позволяет обойти все эти трудности. Вычислив с использованием этого алгоритма переходную матрицу Ф(Ь, Ь^, для регулятора с переменной настройкой можем записать [2]
Р(Ь) = [Ф21(ь, Ь^ + Ф22(ь, Ь^][Ф11 (Ь, Ьf) +
-1
+ Ф12( Ь, ьf) Pf г, Ь < Ьг_
Матрица коэффициентов закона управления для этого случая имеет вид
Установившееся решение уравнения Риккати находится из выражения Р = Т22Т12, где Т12 и Т22 - блоки матрицы преобразования к канонической форме
Т =
Т11 Т12 Т21 Т22
. При этом матрица коэффициентов уси-
ления оптимального регулятора с постоянной настройкой определяется соотношением К к = К 1ОТР. В тех случаях, когда возможностей вычислительных средств оказывается достаточно для символьного обращения матрицы Т12, для элементов матрицы Кк могут быть получены даже символьные выражения, часто помогающие глубже понять суть рассматриваемой задачи. Например, для управляемой системы [2]
х (Ь) =
О 1 О -а
х (Ь) +
ц( Ь),
(1)
элементы матрицы Q, отвечающие за сохранение инвариантов подобия, определяются выражениями
= (^2^2Г1 /Я2 = (^1 + ^2 - а2)Г1 /£2 (2)
где Г1 - элемент скалярной матрицы К, а выражения для элементов матрицы коэффициентов усиления Кк = [&1 £2] оптимального регулятора с постоянной настройкой имеют вид
£1 = ^1^2/£; £2 = (^1 + ^2 - а)/£.
(3)
Анализ выражений (3) показывает, что коэффициенты регулятора не зависят от элементов матриц Q и К, а определяются только параметрами управляемой системы и выбранным распределением корней (ср. с [2]). На основании формул (2),(3) могут быть сделаны также обоснованные выводы о параметрах этого распределения. Ясно, например, что для значений должны выполняться условия
222 + > а; + ^2 > а .
(4)
С учетом того, что ^ = а + Ы, ^ = а - Ы условия
(4) преобразуются к виду а > л/я^^+Ь2.
Решение задачи синтеза алгоритмов оптимального наблюдателя может быть получено на основе свойства дуальности задач оптимального регулирования и оптимального наблюдения. Как известно [2], задачу оптимального наблюдения можно представить в виде задачи оптимального регулирования для системы
1 ^ Т,
КК(Ь) = К О Р(Ь)
• тТ Т х = - Г х - С и
с критерием качества
/0 = £0) БоХ (¿о) + 1( хТОх + иТ Яи) М.
Ч
Комбинируя полученные для регулятора и наблюдателя результаты, можем для системы (1) записать уравнения оптимального регулятора с обратной связью по измеренным значениям выходной переменной у = = СХ + и [2]
Матрица расширенной системы (12) для данной задачи имеет вид
А =
Т Т -1
-гТ с'к С
О
г
В соответствии с базовым алгоритмом для нее может быть вычислена переходная матрица Ф(Ь, ¿о). С использованием этой матрицы решение матричного уравнения Риккати в задаче наблюдения
5 (Ь) = ГБ( Ь) + Б (Ь) ГТ + ОТ -1
- Б(Ь)СТЕ 1 СБ(Ь), Б(Ь0) = Бо
Х= (Г - ОКк-Кд С)х + КдУ, Х ( Ьд) = Хд,
и = -Ккх.
Пример 2. Решается задача оптимального управления движением мостового крана. Уравнения движения, все обозначения и величины параметров приняты в соответствии с работой [6]. Дана линейная стационарная система вида Х = Гх + Ои, где
о 1 о о о
Г = -а о о о , О = е
о о о 1 о
-С о о о е
,е = Ь/ип
записывается в виде
Б(Ь) = [Ф21(Ь, Ьо) + Ф22(Ь, Ьо)Бо]х х[Фц(Ь, Ьо) + Ф12(Ь, Ьо)Бо]Л Ь > Ьо.
Для нахождения установившегося решения этого уравнения используется соотношение
- -1 Б = Т22Т12.
Оптимальные коэффициенты наблюдателя опреде-
Т -1
ляются при этом выражением Ко(Ь) = Б(Ь)С Я , Ь > Ьо
или Ко = БСТЯ 1
Для рассматриваемого примера (1) элементы матрицы О в задаче наблюдения определяются выражениями
ql = (^2 + }?2 -д2)г\> Я2 = ((^1 + )а2 + )*1.
Выражения для коэффициентов Ко = ^ к^ оптимального наблюдателя с постоянной настройкой имеют вид
2
&1 = (^1 + Х^)-а; ^2 = а -(^ + ^)а + ^1^2.
Отметим, что эти коэффициенты также не зависят от матриц О и Я (ср. с [2]). На основании анализа полученных выражений могут быть также сформулированы условия вида (4), которым должны удовлетворять параметры распределения корней наблюдателя.
Граничные условия (2), (3) в рассматриваемом примере имеют вид
со = [о, о, о, оТ, СЧ = [о, о, sf/I, оТ, М = I.
Ограничение, наложенное на управление и, учитывается в данном случае путем выбора такой величины Ьч , при которой условие \и\ < 1 заведомо выполняется
на всем интервале управления. Другие способы учета ограничений рассматриваются в следующих примерах.
Продемонстрируем на этом простом примере все этапы решения краевой задачи. Прежде всего, выберем в качестве корней характеристического полинома замкнутой системы корни распределения Баттерворта (37) для п = 4 и представим каноническую форму матрицы А в виде диагональной символьной матрицы П = = 4,-^1,^,-^4). Уравнения для опреде-
ления параметров ..., Я4 в данном случае имеют вид
2 2 2 2 2 (2Г1 а-е (у2 + Я4)) = -(^1 + + ^3 + ^4)г1>
(Я1 + Яз-2°Я4)е + а Г1 = (^2^2 + ^3^4 + + ^2^3 + Х2 + ^4)^*1 ;
2 222 222 222 222 (2^3 - ОЯ4)ье = -(^1 ^3^4 + ^2^4 + ^1^2^3 + ^2^3^4)1;
®2е2^з = .
Ь
Для сокращения записи в этих формулах введено обозначение V = а — с. В результате решения этой системы линейных уравнений находим
= Х^^^Х^ / V / в2;
д4 = ((Х^Х3Х^ + Х^Х2;Х^ + Х^Х^ + Х^Х^Х^)г^/V/в2-
+2д3)/V;
д^ = ((Х2Х4 + Х;Х4 + Х1 Х; + + Х; Х3 + Х.^4 — а )Х
х Г! — дз + 2 vд4)/в
д2 = (Х; + Х4 + Х2 + + 2а) Г1 / в2 — д4
Обозначим столбцы матрицы преобразования Т к канонической форме через к = 1, ..., 2п. В соответствии с алгоритмом [2], символьные выражения для этих столбцов находятся путем решения 2п систем линейных уравнений вида
(А — Хк1 )Ьк = 0, к = 1, ..., 2п.
(5)
Проще всего решение этих уравнений может быть получено следующим образом. Вначале для каждого к путем вычеркивания в символьной матрице (А — Х^) строки и столбца с номером к формируется квадратная матрица Вк порядка 2п - 1. Вычеркиванием элемента с номером к в векторе ^ вводится новый вектор Vk. Полученная в результате система уравнений с невырожденной матрицей Вк легко решается относительно Vk средствами символьного пакета МАТЬАВ. По результатам решения формируется искомый вектор причем вначале элементу в позиции к этого вектора присваивается значение, равное 1.
Полученные при этом выражения зачастую оказываются очень громоздкими. Чтобы получить более про-
стые выражения, можно воспользоваться тем, что векторы Ьк из уравнений (5) определяются неоднозначно, и умножить каждый из них на любое подходящее выражение. Например, таким способом в выражениях для элементов ^ можно избавиться от знаменателей, которые как раз и придают полученным первоначально выражениям упомянутый необозримо громоздкий вид. В рассматриваемом примере окончательные выражения для Ьк могут иметь, например, вид, определяемый уравнениями
2 2 2 в (Хк + v)(— дз + д4Хк)
2 2 2 Хкв (Хк + — дз+д4Хк)
2 2 2 2 2 Хк((Хк +а) Г1 +в (д1—Хкд2))
з 2 2 2 2
Хк((Хк +а) Г1 +в (д1—Хкд2))
2 2 2 —Хк(с(Хк +а)Г1 +в (д1 + vд2))(— дз + д{Хк)
2 2 2 2 (с(Хк + а)Г1 + в (д1 — Хкд2))(— дз + д4Хк)
2 2 2 2 —Хкдз((Хк +а) Г1 +в (д1—Хкд2))
2 2 2 2 2 — ((Хк + а) Г1 +в (д1— Хкд2))(— дз + д4Хь)
Ч =
к = 1, ..., 2п.
(6)
Вычисления на последующих этапах осуществляются в соответствии с базовым алгоритмом.
Результаты моделирования движения крана при ^ = = 4с приведены на рис. 1, 2.
Пример 3. Рассматривается задача регулирования из статьи [1]. Дана линейная система х = ¥х + Оы, матрицы которой имеют вид
^ =
0 0 1 0 0 0 0
0 0 0 1 , О = 0 0 0
—1 1 0 0 1 0 —1
0,1 —1, 02 0 0 0 —1 1_
1.5
0.5
-0.5
;
1 г 1 1 х4 : ■""Мь- - - »----.....
/Л
I, с
Рисунок 1
Рисунок 2
Движение системы рассматривается на интервале Ь е
Т
е [о, 15]. Граничные условия имеют вид: Со = [2, 1, 7, 5] ;
Т
Сч = [о, о, о, о] ; М = I. В качестве канонической формы в данном примере принята диагональная матрица.
На управляющие функции и^(к = 1, 2, 3) наложены ограничения вида итШ(Ь)< ик(Ь) < итах(Ь). Учесть эти ограничения в постановке оптимизационной задачи невозможно, т. к. краевая задача при этом станет нелинейной. Поэтому используется способ апостериорного учета этих ограничений, заключающийся в том, что для каждого Ь величины ик (Ь) вычисляются вначале в соответствии с базовым алгоритмом (т. е. так, как будто никаких ограничений нет), а затем они подвергаются нелинейному преобразованию
иН(Ь) = т1п{итах(Ь), тах[итт(Ь),и£(Ь)] } (7)
в результате чего формируются величины управлений ик (Ь), фактически подаваемые на исполнительные органы. Результаты моделирования для случая итШ = = -1, итах = 1 приведены на рис. 3-6.
Пример 4. Рассматривается задача построения алгоритмов автопилота с эталонной моделью, обеспечивающего управление боковым движением летательного аппарата с полной развязкой каналов крена и рыс-
кания [7]. Угловое движение аппарата описывается уравнениями
ю = а11( Ь)ю + ац( Ь)Р + а1з( Ь )Р + Ь1 (Ь) и^
|3 = аз1( Ь )ю + аз2( Ь )Р + азз( Ь )(3 + Ь2( Ь )и2, (8)
где ю - угловая скорость по крену; р - угол скольжения; и1, и2 - углы отклонения рулевых органов; коэффициенты а- (г = 1, 3; / = 1, 2, 3), Ьк (А = 1, 2) изменяются по времени с изменением высоты и скорости полета.
В качестве эталонной модели принимается система уравнений вида
г = Сг + д,
(9)
• Т Т
где г = [юш, Pш, Рш] ; д = [gl, ^ д2] ; gl, д2 - управляющие воздействия по крену и рысканию, поступающие от ручки и педали управления соответственно. Элементы матрицы
с11 о о
о о 1
о С32 С33
являются заданными постоянными величинами.
Рисунок 3
Рисунок 4
С=
Решение уравнения (9) дается формулой Коши
£
г(£) = £, £0)г(£0) + £, т)д(т)Л. (10)
Переходная матрица £, £о) легко находится с использованием средств символьного пакета МАТЬАВ
£, £о) =
0 ¿2е3- ¿3е2 е2
0 с32( е2 — е3) ¿2е2- ¿3е3
где
ек = ехр[£-£0)], к = 1, 2, 3;
¿1 = Сц; ¿2 = (С33 + с)/2; ¿3 = (С33 -с)/2;
с = л/С33 + 4с32.
Система (8) являются нестационарной, поэтому непосредственное применение к ней базового алгоритма невозможно. Мы поступим следующим образом. Разобьем время управления на циклы постоянной длины к. Величина к определяет время реакции системы на действия пилота, поэтому она должна выбираться минимально возможной (например, к < 0,1 с). На малом интервале времени от момента £0 = £ до момента £к = £ + к исходную модель управляемой системы (8) заменим более простой моделью, представляющей собой два независимых интегратора вида ю = »1 и в = »2, где 01, »2 - новые управляющие переменные. В таком подходе нет ничего необычного - замена на малом интервале времени сложных моделей более простыми математическими объектами (полиномами, сплайнами и т. п.) широко используется в вычислительной математике, например, при построении алгоритмов численного интегрирования.
Представим новую расширенную модель системы в виде
где
х = Гх + ОV,
х = [ю,ю,в, в, в] ; V = [»1, »2]
Г =
0 1 0 0 0 0 0
0 0 0 0 0 1 0
0 0 0 1 0 ; о = 0 0
0 0 0 0 1 0 0
0 0 0 0 0 0 1
В начальный момент £0 = £ координаты ю, в, в вектора с0 в уравнении (2) могут быть найдены по результатам измерений фактических значений параметров < и в путем их численного дифференцирования либо получены в качестве выходных параметров некоторого наблюдателя. Координаты вектора с^ в конечный момент £! = £к определим следующим образом. Располагая информацией о значениях параметров ^, д2 эталонной модели в к моментов времени £0 = £, £1 = £-к, ..., £к = £-(к- 1)к, ..., мы можем вычислить экстраполированные значения этих параметров в момент £к = £ + к. Тогда значения параметров в к, вк вычисляются по формуле (10), а параметров ю^,вк - по уравнениям эталонной модели (9). Сформированные подобным образом граничные условия С0 и с^ используются в базовом алгоритме для вычисления по формулам (33) программной траектории х*(£) перехода от фактического углового положения летательного аппарата к задаваемому эталонной моделью требуемому положению. Вычислять вектор р*(£) при этом нет необходимости - полученное с использованием этого вектора управление (34) для нас совершенно бесполезно, т. к. оно не связано ни с каким физически реализуемым управлением в исходной задаче. Поэтому в данном случае используется другой способ вычисления требуемых управляющих функций «1 и «2. Заметим, что в числе программ х*(£) вычисляются функции х2(£) = ю*(£) и х*(£) = в*(£). Полагая в уравнениях
(8) фактической модели ю( £) = ю * (£), в( £) = в*( £), значения параметров управления найдем из уравнений
«1 = (со *-(вцю + а^в + а^в))/&1;
П2 = (в* -(^Ю + ^32в + в33в))/
При таком подходе к выбору параметров управления движение летательного аппарата осуществляется в соответствии с уравнениями эталонной модели (9), т. е. оказывается не зависящим от любых конкретных (но известных) значений параметров его модели. При этом в качестве модели летательного аппарата (8) могут быть использованы уравнения любой сложности, а не только линеаризованные уравнения.
Моделирование процесса управления проводилось для случая, когда коэффициенты модели изменяются во времени произвольным образом. В данном случае рассматривался вариант задания этих функций в виде отрезков синусоид со случайно выбранными частотами, амплитудами и фазами.
Результаты моделирования, приведенные на рис. 7, 8, свидетельствуют о том, что при описанном подходе обеспечивается практически идеальное отслеживание
0
се
е
3
эталонной траектории, задающей независимое движение летательного аппарата в каналах крена и рыскания. На рис. 9 приведены отклонения фактических значений угла скольжения от его эталонных значений в процессе маневра. Всплески погрешностей отслеживания эталонных значений на графике соответствуют точкам, в которых производная функции д2 претерпевает разрывы, вследствие чего возрастают погрешности экстраполяции значений этой функции на моменты th. На рис. 10 представлен характер изменения в ходе маневра величин углов отклонения элеронов («1) и руля направления («2).
Пример 5. Рассматривается задача перевода КА с ма-ховичной системой управления, находящегося на круговой орбите, из произвольного неориентированного положения в положение, соответствующее полету в орбитальной системе координат (ОСК). Цель данного примера - показать, каким образом базовый алгоритм (алгоритм терминального управления с фиксированным временем) может использоваться для решения задач регулирования и стабилизации, где время заранее не задается или даже может быть вообще неограниченным.
Угловое движение КА относительно ОСК описывается уравнениями
/ю = Н х (ю + Юо) + /(ю х Юо - юо) - т;
У = Гю, (11)
где Юо - вектор орбитальной угловой скорости, остальные обозначения совпадают с [4]. Для решения задачи воспользуемся известным преобразованием [4], суть которого заключается в переходе от множества переменных {ю, у, т} исходной задачи к новому множеству {у, у, и}. Система (11) при этом точно преобразуется в эквивалентную линейную систему вида х = ¥х + Ои,
где х =
Таким образом, система (11) представляется в виде двумерного интегратора. Граничные условия Со и е^ для этой новой системы задаются следующим образом. В текущий момент времени t, принимаемый за начальный tо, координаты вектора Со определяются по данным системы определения ориентации и показаниям измерителей угловой скорости. Вектор е^ задается в виде е! = о, матрица М = I.
х; У ; ^ = о 1 ; о = о
х2 x о о _1
■4
; ;
\ > /
\ \® г /
«Д / \/
\ ■ [г ■
* ;
6 I, с £
Рисунок 7:
д1 - управляющее воздействие от ручки управления; ю - эталонная и фактическая угловые скорости крена
Рисунок 8:
д2 - управляющее воздействие от педали управления; в - эталонный и фактический узлы скольжения
И;, Щ 0.6
0.4 02 О -Й 2 -0.4
:
(\1 ■
И) * * # 1 н / 1 г > "А Л....... -О
\ Л? \ / Ч X л у/ „„V,/! —
Я I, с
Рисунок 9
Рисунок 10
Для использования базового алгоритма требуется также задание времени достижения заданного терминального состояния е^. В задачах регулирования, где это время обычно неизвестно, значение может выступать в качестве одного из параметров, регулирующих скорость протекания переходного процесса.
Возврат от управления и в линейной задаче к управлению т в исходной задаче осуществляется по формуле, аналогичной [4]
-1 -1 т
т = /Г и + к х (ю + ю0) - 1 /2у0 (ю ю)/у +
0.6
+ /(ю х ю0 - юо),
0.4
0.2
-0.2
;
11 1 У =17 5с \
Ъ Л л | Ч !
,А/=5'0с
200
400
600
^00, 1000 1, с
Рисунок
где
то = V1 - хТх1; У = Х1;
Ю 2 (Уо + То х 1 х 1 х2 — х 1 х х2).
Учет ограничений на величины управляющих моментов маховиков осуществляется путем использования нелинейного преобразования (7).
Результаты моделирования процесса переориентации для двух вариантов задания приведены на рисунках 11, 12. Видно, что процесс регулирования, длящийся примерно 400-800 с плавно переходит в процесс стабилизации, который теоретически может продолжаться неограниченно долго.
Пример 6. Для многих управляемых процессов характерной является ситуация, когда общая задача управления разбивается на несколько элементарных подзадач. Траектория движения в этом случае оказывается состоящей из нескольких участков. При этом важной является проблема обеспечения плавной стыковки этих участков. Рассмотрим один из способов решения данной проблемы на примере задачи кадровой фотосъемки участков земной поверхности с борта КА, находящегося на околоземной орбите.
На рис. 13 приведены результаты моделирования процесса управления съемкой участка поверхности Земли размером 15 х 15 км (3 строки по 3 кадра, размер кадра 5 км) в предположении, что, как и в предыдущем примере, для перехода в пространство линейных параметров используется модель движения в виде совокупности независимых интеграторов. Как видно из графиков, на границах смежных участков производные угловых скоростей в каналах крена (Ю1) и тангажа (Ю2) терпят разрывы. При этом создаются неблагоприятные условия для системы стабилизации углового движения КА относительно программных траекторий. Плавная стыковка отдельных участков может быть обеспечена следующим образом. Добавим к вектору состояния эквивалентной линейной системы (6) координату хз = у и изменим соответствующим
-2
-4
1 ^¿■=175/:
\ 1 1 \ 1
\
\£у=50с\
О 200 400 600 800, 1000
% С
Рисунок 12
аз, градус/с 0.6
0.4 0.2 О -0.2 -0.4
10 20 Рисунок 13
30 4 с
образом матрицы Г и О. В исходной задаче граничные условия на эту координату не наложены, поэтому мы можем распоряжаться ими по своему усмотрению.
В частности, мы можем потребовать, чтобы в моменты стыковки смежных участков совпадали не только значения координат х1, х2, но и значения координаты х3. Тем самым, будет обеспечена требуемая плавность изменения параметров ю1 и ю2. Результаты моделирования для этого случая приведены рис. 14, 15 (пунктирные линии на графиках соответствуют линейной модели в виде 2-мерного интегратора, сплошные - расширенной модели в виде 3-мерного интегратора). Графики управляющих функций показаны
Рисунок 14
10
Рисунок 15
10'
10 20
Рисунок 16
2 1 и -1
-2 -3
¡0
: /Л [| ||
..........1/____
т—1 -ч 1-----
я
! ь 3 -Ц-. шМ
1
20
Рисунок 17
30 I, с
на рис. 16, 17. Из рисунков видно, что предложенный подход обеспечивает более благоприятные условия для работы системы стабилизации. При дальнейшем расширении вектора состояния может быть обеспечена любая желаемая степень гладкости.
Пример 7. Для традиционных подходов к задачам оптимального управления характерно почти полное отсутствие возможностей учета априорной информации об управляемом процессе. Между тем, зачастую, еще только приступая к решению задачи, мы уже достаточно много знаем о ней. Например, при решении задач наведения в космосе мы можем принять во внимание хорошо известный результат, полученный ранее многими исследователями, - для маневров в гравитационном поле оптимальными являются режимы полета с максимальной тягой. С учетом сказанного, рассмотрим задачу терминального управления на примере выведения КА в заданную точку эллиптической орбиты в заданное время. Для простоты будем считать, что движение происходит в центральном поле (как мы увидим далее, вид гравитационного поля особого значения не имеет). Уравнения движения КА относительно инерциального базиса имеют вид
единичный вектор направления
СО80СО8^ где е = 8Ш0СО8^
силы тяги; а = ис/(т - Ь) - модуль ускорения силы тяги; 0 - угол тангажа; у - угол рыскания; ис -скорость истечения; т = /с - время полного «выгорания» топлива; - начальная масса; с - секундный расход; д = -ог - вектор гравитационного ускорения; V = Г3
Чтобы привести исходное нелинейное уравнение (12) к эквивалентной линейной системе введем новые переменные х1 = г, х2 = Г, х3 = Г, х4 = г и новое управление и. Запишем новые уравнения движения в виде трех (по числу каналов управления) независимых четырехмерных интеграторов. Уравнения движения в каждом канале управления при этом приводятся к уравнениям вида х = Гх + Ои с матрицами
0 10 0 0
Г = 0 0 10 , О = 0
0 0 0 1 0
0 0 0 0 _1
Г = ае + д,
(12)
Теперь для того чтобы применить базовый алгоритм, осталось сформировать граничные условия с0
и С/. Будем считать начальные условия движения г0, г0 известными, конечные условия т/, т/ и время наведения - to) - заданными. Заданными будем считать также начальные и конечные значения углов тангажа 00, 0/ и рыскания V/. Располагая этими данными, мы можем непосредственно сформировать граничные условия для переменных х\ и %2. Для переменной х3 граничные условия определяются из уравнения (12). Чтобы задать указанные условия для переменной х4, найдем по уравнению (12) производную
т = ае + ае - Vт - Vт,
где
-8Ш 0ТО8у ТО808Ш у а =
е = 0 сО80ТО8у - Т ; V
_ 0 _ сО8у р =
т• ,, ,2
Будем предполагать, что на значения параметров 0, у в моменты времени tf никаких условий не наложено, и мы можем распоряжаться ими по своему усмотрению. Выберем их наиболее естественным способом:
0о = 0/ = (0/-0о)/(/- to),
Р = Т/ = (Т/ - То)/(/- tо).
Наконец, зададим матрицу М в уравнении (3) равной единичной матрице и сформируем векторы Со и С/. В заключение, выберем корни характеристического уравнения в соответствии с распределением Баттер-ворта и зададим радиус круга этого распределения, равным 1.
В процессе работы базового алгоритма в числе программ х*(^ вычисляются функции х1 (^ = т*(t) и х3(t) = т*(t). На основании этих данных программные значения величин ускорения силы тяги а*(0 и углов ориентации вектора тяги 0*(t), у*(t) в момент t определяются следующим образом.
На основании уравнения (12) можем записать х3( t) = а (t) е (t) + д [ х1( t)].
Отсюда получаем
а*( t) = »¡[^Л-д^у^^^)-^^]; е*(t) = [х*(t) - д*(t)]/а*(t).
Окончательно имеем
у*( t) = -агсзт[ е3( t)]; 0*( t) = агйап[ е2( t)/e* (t)].
На рис. 18, 19 приведены результаты моделирования процесса наведения для следующих значений параметров КА и участка наведения [5]:
^ = 0; / = 387 с; т0 = 14000 кг; е = 25 кг / с; ис = 4, 5 км/с; 00 = 47, 35°; 0/ = 96, 23°; у0 = 55, 4°; у/ = 0; т0 = [6372, 9;-1570, 0;-299, 4] км; т0 = [ 1, 502;5, 477;1, 885] км/с; т! = [6510, 663510927766;1238, 442814248355;0] км; т/ = [-1,105899452214530;10,13460122429019;0] км/с.
При вычислениях с разрядностью мантиссы 16 десятичных знаков погрешности выведения КА составляют: по координатам ~3,110-3 м, по скоростям ~2,910-8 м/с. Вычисления с более высокой разрядностью обеспечивают еще большую точность выведения. В этом, конечно же, нет ничего удивительного, так как мы решили краевую задачу аналитически, т. е. совершенно точно. Погрешности при этом определяются только погрешностями численной реализации процесса вычислений программного движения по этим формулам и дей-
Рисунок 18: 1 - заданная орбита; 2 - траектория выведения
100 80 60 40
20
О
— *
100
200
100 I, с
Рисунок 19:
1 - угол тангажа; 2 - угол рыскания
ствующими возмущениями. Для компенсации возмущений управление в линейной задаче должно формироваться в виде u (t) = u*( t) + 5 u (t). Поэтому KA должен быть оснащен двигателем с регулируемой в небольших пределах тягой.
Интересно сравнить рассмотренный в данном примере метод наведения с лучшим на сегодняшний день явным методом [5, 8], который впервые был применен в 60-70-х годах прошлого века для наведения ракеты-носителя Saturn V в ходе реализации программы полетов KK Apollo на Луну. В настоящее время этот, несколько модифицированный, метод применяется для наведения большинства современных ракет-носителей (Space Shuttle, Arian и др.). Главное отличие заключается в том, что с использованием явного метода обеспечивается выполнение заданных конечных условий лишь для пяти параметров движения (отсутствует управление продольной координатой). При этом время выхода на заданную орбиту и положение точки на этой орбите, в которой оказывается KA в момент конца участка выведения, зависят от возмущений и заранее неизвестны. Метод наведения с использованием базового алгоритма обеспечивает выполнение заданных конечных условий по всем шести параметрам движения KA точно в заданное время. Данное свойство позволяет, дополнительно к тем возможностям, которыми обладает явный метод наведения, решать многие качественно новые задачи, возникающие в ходе выполнения маневров в космосе. Например, при выполнении операции по спасению космонавтов с терпящего бедствие KA решающим может оказаться фактор времени. С использованием предлагаемого метода оказывается возможным непосредственное решение задачи встречи на орбите, минуя занимающие значительное время операции фазирования орбит. При этом обеспечивается столь высокая точность наведения, что практически осуществимым может оказаться выведение KA непосредственно в ту область пространства, из которой далее могут начинаться этапы причаливания и стыковки.
^оме того, в рассматриваемом подходе краевая задача решена точно, без каких бы то ни было упрощений исходной модели движения, которые надо было бы компенсировать дополнительным расходом топлива. Поэтому данный подход обеспечивает более высокую степень оптимальности управления по сравнению с явным методом наведения. Наконец, что тоже немаловажно, алгоритмы рассматриваемого метода значительно проще, чем алгоритмы явного метода наведения.
Пример 8. Рассматривается применение базового алгоритма для высокоточного численного интегрирования уравнений небесной механики. Традиционно методы численного интегрирования являются предметом вычислительной математики. Подход к задачам интегрирования как к задачам оптимального управления позволяет по-новому взглянуть на многие проблемы, свя-
занные с численным интегрированием. Рассмотрим особенности предлагаемого подхода на примере интегрирования периодических орбит в круговой ограниченной задаче трех тел. Уравнения движения имеют вид [3]
.. п . , 50 .. п . , 50 л: = 2у + —; у = 2х +—, 5х ёу
где
0.'[ (i - ^ + 2) + „(Р2 + -2)}
2 , , -2 , 2 2 ,, ,-,2,2
Р1 = (х + + У ; Р2 = [х - (1 - 1)] + У I = 0, 012277471.
Начальные условия интегрирования:
х0 = 0, 994; х0 = 0;
у0=2,0317326295573368357302057924;
период Р = 11,124340337266085134999734047. Орбита, соответствующая этим начальным условиям, приведена на рис. 20.
Как видно из рисунка, в окрестности точки 2 наблюдается тесное сближение движущегося тела с Луной. Именно на этом участке при численном интегрировании происходит основная потеря точности. Поэтому интегрирование подобных орбит является хорошим испытанием для любого метода интегрирования. Уравнения (13) являются нелинейными. Ясно, что непосредственно к ним базовый алгоритм применить невозможно. Поэтому мы вновь воспользуемся возможностью
1
□.в □.6 0,4 0.2 0 0.2 ■0,4 -0.6 -0.В -1
1 1
■
¡
Л ■ ---
1 __ 2
-0.5
0.6
X
Рисунок 20: 1 - Земля; 2 - Луна
точного приведения этих уравнений к эквивалентной линейной системе. Введем новый «-мерный вектор состояния
т т т т г = , г2,..., г«] ,
Г , г- -п г («-!) («-
где г1 = [ху]; г2 = [х>у]; ...> г« = [х ,у ].
В этих новых переменных уравнения (13) представляются в виде «-мерного интегратора г. = г2, г2 = = гз, ..., г«-1 = г«, г« = и, где и - некоторое фиктивное управление. Ввиду недостатка места символьные выражения для координат г4, ...,г« здесь не приводятся, однако они легко могут быть получены путем дифференцирования уравнений (13) с использованием средств символьного пакета МАТЬАВ.
Граничные условия с0 = г (£0) в начале каждого шага интегрирования вычисляются с использованием уравнений (13) и полученных символьных выражений для высших производных. Граничные условия в конце шага интегрирования задаются следующим образом.
Поскольку в момент Ь/ = £0 + к требуемые значения параметров г^Ь/), г2(Ь/) неизвестны, в связи с чем граничные условия для них задать невозможно, выберем матрицу М в уравнении (3) в виде единичной матрицы порядка « с вычеркнутыми из нее двумя первыми строками. Интегрирование будем вести по обычной для неявных методов интегрирования схеме «прогноз плюс коррекция». На самом первом шаге интегрирования прогнозированные значения г1 (Ь/), г2(Ь/) вычисляются с использованием разложения в ряд Тейлора
2 з
г1 (Ь/) = г^Ьд) + г2(£0)к + гз(Ь0)к /2 + г4(Ь0)к /6 +...;
г2 (Ь/) = г2( Ьд) + гз( Ь0)к + г4 (Ь0) к2 / 2 + ....
На последующих шагах интегрирования прогноз осуществляется с результатами использования базового алгоритма, причем в качестве ^ берется его значение, полученное на предыдущем шаге.
По вычисленным таким образом значениям г1 (Ь/), г2( Ь/) с использованием полученных ранее символьных выражений находятся начальные приближения для оставшихся координат вектора г(Ь/) и формируется вектор с/ = [гз(Ь/), ..., г«(Ь/)]. На этом подготовка к ис пользованию базового алгоритма завершается.
В результате работы базового алгоритма на этапе коррекции вычисляются некоторые новые значения г1(Ь/) и г2(Ь/), по ним вычисляется новые значения координат вектора с/ и т. д. Как показывают численные эксперименты, процесс этот сходится за 2-з итерации.
Моделирование процесса интегрирования проводилось на интервале времени, равном одному периоду Р, для двух значений порядков расширенной системы: « =6 и « = 7. При « = 6 точность интегрирования характеризуется следующими величинами погрешностей:
-1з
- по координатам Ът = 4 • 10 ,
-11
- по составляющим скорости Ът = 6 • 10 .
Такая точность сопоставима с точностью, обеспечиваемой лучшими на сегодняшний день методами интегрирования [9]. При той же последовательности шагов интегрирования для случая « = 7 погрешности составляют соответственно Ът = з • 10 18, Ът = 5 • 10 16. Как видим, увеличение порядка метода на единицу приводит к уменьшению величин погрешностей на пять порядков. Это свойство предлагаемого подхода к задачам численного интегрирования, казалось бы, указывает путь к возможности практически неограниченного повышения точности интегрирования. Однако, в настоящее время этому препятствует уже упоминавшееся выше «проклятие разрядности».
Пример 9. При рассмотрении предыдущих примеров мы убедились в том, что в случае, когда уравнения исходной задачи являются нелинейными, расширение вектора состояния путем включения в него нескольких высших производных позволяет построить эффективные законы управления, обладающие рядом ценных свойств. Если исходная система уравнений уже является линейной, к ней может быть применен точно такой же подход с расширением вектора состояния.
Пусть исходная система задана в виде системы второго порядка
х 1 = /11х1 + /12х2; х 2 = /21 х1 + /22 х2 + 9и.
Для моделирования в качестве такой исходной системы возьмем систему (1), уже рассматривавшуюся в примере 1. Зададим интервал управления, равным 10 с, и граничные условия в виде С0 = [ 5 ;-5 ]; с/ = = [-1 ;1 ]. Результаты моделирования процесса управления для этой системы второго порядка, полученные с использованием базового алгоритма, приведены на рис. 21, 22.
Как видно из рис. 21, для процесса управления в исходной задаче характерным является приложение максимальных значений и в начале и в конце интервала управления. Между тем ясно, что такой характер изменения управляющей функции по ряду причин на практике может оказаться неприемлемым. Например, на управляющую функцию могут быть наложены ограничения вида ит;„ < и < ит„При наличии подобных
ннп шах.
ограничений решение нелинейной краевой задачи получить вообще невозможно.
Рисунок 21:
1 - управление в исходний задаче; 2 - то же в расширенной задаче
1, 2
2 3
Рисунок 22:
координата и скорость в исходной задаче; з, 4 - то же в расширенной задаче
Расширим вектор состояния системы (14), для чего введем обозначение /2з = 9 и добавим к вектору состояния новую переменную хз = и. В исходной задаче никаких ограничений на эту переменную не наложено. Следовательно, мы можем распорядиться ее граничными значениями, а также граничными значениями нескольких ее высших производных, по своему усмотрению. Ограничиваясь для простоты первой и второй производными, запишем вместо системы (14) расширенную систему
х = Гх + О»,
где » - новое управление, а матрицы Г и О имеют вид
0 1 0 0 0
Г = 0 -а £ 0 ; о = 0
0 0 0 1 0
0 0 0 0 _1
Пользуясь свободой выбора граничных условий для переменных хз и х4, сформируем граничные условия для этой новой модели управляемой системы в виде
0
= [5; -5; 2е - з; 0]; с/ = [-1; 1; 6е - з; - 5е - 4}
Результаты моделирования процесса управления с данными граничными условиями, представленные на рис. 22, показывают, что управляющая функция и = = хз имеет в данном случае гораздо более приемлемый для практической реализации вид. Это свидетельствует о том, что подход с расширением вектора состояния линейных систем в сочетании с использованием базового алгоритма позволяет оказывать активное влияние на характеристики управляемых процессов. С использованием этого подхода в большинстве случаев можно добиться существенного улучшения указанных характеристик.
выводы
Сочетание основных принципов классической теории управления (свободный выбор распределения корней системы и расширение вектора состояния для улучшения свойств замкнутой системы), принципа максимума (базовый алгоритм) и метода динамического программирования (поэтапное планирование операций) в рамках единого подхода, разработанного в первой части, во многих случаях позволяет получить решение задачи синтеза оптимальных управлений как для линейных, так и для нелинейных нестационарных систем, в том числе с ограничениями на управляющие функции.
БЛАГОДАРНОСТИ
Данная работа, состоящая из двух частей, своими корнями уходит в далекие 60-е годы прошлого века. Тогда автору, только что окончившему ХАИ, выпала редкая удача - в самом начале своего профессионального пути поработать в уникальном научном коллективе на предприятии, которое теперь называется «Харт-рон». В те годы на этом вновь организованном предприятии все были молоды, веселы и очень талантливы. Уже тогда они искали свои подходы к решению проблем управления, закладывая фундамент научной школы, получившей со временем всеобщее признание. Именно этим коллективом впоследствии, в числе многих уникальных проектов, была разработана система управления стратегической ракеты, известной теперь под названием «Сатана».
Общение с такими выдающимися специалистами как Я. Е. Айзенберг, |А. С. Гончар, С. С. Корума, |В. Н. Романенко, В. Д. Стадник и др. наложило глу-
бокий отпечаток на мировоззрение автора, круг его научных интересов и, в конечном счете, на его судьбу. Всем этим замечательным людям - живым, а также тем, кого, к величайшему сожалению, уже нет с нами, -автор выражает свою глубокую искреннюю благодарность.
ПЕРЕЧЕНЬ ССЫЛОК
1. Габасов Р., Дмитрук Н. М., Кириллова Ф. М. Оптимизация многомерных систем управления с параллелепи-педными ограничениями // АиТ. - 2002. - № 3. - С. 3-26.
2. Квакернаак X., Сиван Р. Линейные оптимальные системы управления. - М.: Наука, 1977. - 652 с.
3. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. - М.: Мир, 1990. -512 с.
4. Дуайер Т. А. V., III. Точное нелинейное управление быстрыми вращениями КЛА посредством внутренней передачи количества движения // Аэрокосмическая техника. - 1987. - № 3. - С. 151-159.
5. Синха С. К., Шривастава С. К. Оптимальный закон наведения многоступенчатой ракеты-носителя на трехмерной траектории, использующий выражения в явном
виде // Аэрокосмическая техника. - 1991. - № 3. -С. 74-87.
6. Кабанов С. А. Управление системами на самоорганизующихся моделях // АиТ. - 2001. - № 7. - С. 122-128.
7. Земляков С. Д., Рутковский В. Ю. О некоторых результатах развития теории и практического применения беспоисковых адаптивных систем. // АиТ. - 2001. -№ 7. - С. 103-121.
8. Чендлер Д. С., Смит И. Е. Схема итеративного управления и ее использование для различных аппаратов и космических операций // Астронавтика и ракетоди-намика. - 1967. - № 8. - С. 1-24.
9. Everhart E. Implicit Single Sequence Methods for Integrating Orbits // Celest. Mech. - 1974. - V. 10. - P. 35-55.
Надшшла 15.07.06
На конкретных прикладах розглядаються особливост1 застосування базового алгоритму, побудованого у першш частит статт1. Для систем з обмеженнями на керуюч1 функцИ використовуеться принцип поетапного плануван-ня методу динам1чного програмування.
It is considered on concrete exemples the peculiarity of use of the base algorithm, wich was synthesized in the furst part of the articl. For system with restrictions of the control function it is utilized the principle of stage by stage planning of dynamic programming method.
