CC BY
17
ТЕОР1Я I МЕТОДИ АВТОМАТИЧНОГО УПРАВЛШНЯ
ТЕОРИЯ И МЕТОДЫ АВТОМАТИЧЕСКОГО УПРАВЛЕНИЯ
THEORY AND METHODS OF AUTOMATIC CONTROL
УДК 62 - 506
А. К. Новиков
СИНТЕЗ ОПТИМАЛЬНЫХ УПРАВЛЕНИЙ НА ОСНОВЕ АНАЛИТИЧЕСКОГО РЕШЕНИЯ ЛИНЕЙНОЙ СТАЦИОНАРНОЙ ДВУХТОЧЕЧНОЙ КРАЕВОЙ ЗАДАЧИ
Часть 1. Постановка и решение задачи. Базовый алгоритм
Рассматривается построение аналитического решения краевой задачи для линейных стационарных систем без ограничений на управляющие функции. В основу предлагаемого подхода положено, как и в классической теории управления, задание распределения корней управляемой системы. Решение задачи представлено в виде универсального (базового) алгоритма.
Не решай сложную задачу, не решив простую (принцип простоты).
Первая заповедь математического моделирования [1]
ВВЕДЕНИЕ
Задачи, рассматриваемые в теории управления, разбиваются на два класса. К первому из них относятся задачи стабилизации движения относительно программных траекторий или заданных стационарных состоя© Новиков А. К., 2006
ний, ко второму - задачи отыскания этих программных траекторий на основе решения краевых задач. Все успехи в развитии теории управления связаны, в основном, с решением задач первого класса. Что касается задач второго класса, то здесь достижения являются намного менее впечатляющими. Основные методы построения оптимальных программных траекторий (принцип максимума и метод динамического программирования) были сформулированы около пятидесяти лет назад. За это время с использованием указанных методов были решены многие важные частные задачи, однако какого-либо общего конструктивного подхода не найдено до сих пор [2, 3]. Как и полвека назад, непреодолимым препятствием на пути широкого применения метода динамического программирования лежит пресловутое «проклятие размерности». Высокая чувствительность
краевой задачи к погрешностям задания начальных значений сопряженных переменных в принципе максимума сильно затрудняет построение оптимальных траекторий методами «проб и ошибок» даже без учета ограничений на управляющие функции. Так как в любой реальной задаче указанные ограничения всегда присутствуют, подход с использованием принципа максимума оказывается здесь и вовсе бесперспективным. Все это дает основания говорить о наличии определенного кризиса в развитии науки об управлении [2].
По мнению автора, одна из гипотез о возможных причинах данного кризиса может быть сформулирована следующим образом. Множества методов решения задач каждого класса (и внутри классов) являются весьма специфичными и не имеют общих элементов. При этом в каждом множестве существуют задачи, для которых внутри самого этого множества методы решения отсутствуют. Положение, как говорится, патовое. Между тем, история науки знает немало примеров, когда новые идеи и даже новые области знаний возникали как раз на стыке наук. Опираясь на этот факт, в данной работе сделана попытка построения нового подхода к проблеме решения краевых задач, органически сочетающего в себе основные принципы классической теории управления, принципа максимума и метода динамического программирования. Основной особенностью данного подхода является аналитическое решение краевой задачи для специального класса линейных стационарных систем без ограничений на управляющие функции.
Предлагаемый в работе подход существенно опирается на возможность проведения символьных преобразований в средах современных систем программирования MATLAB, Maple V и др. Для понимания работы требуются лишь минимальные знания в области матричного анализа, поэтому она вполне доступна широкому кругу инженеров, занимающихся разработкой систем управления.
1 ПОСТАНОВКА ЗАДАЧИ
Движение системы рассматривается на интервале времени Ь е [ Ьо, . Начальный £0 и конечный ^ моменты времени считаются заданными. Граничные условия для системы (1) задаются соотношениями
x(t0) co,
Mx( tf)
cf.
(2)
(3)
Постоянные векторы с0 е К , е^ е К и постоянная матрица М е К х п ранга 0 < £ < п считаются заданными. В качестве критерия оптимальности перевода системы (1) из состояния (2) в состояние (3) принимается классический функционал вида
J = 2 J(xTQx + u Ru)dt,
(4)
где О е Кп х п, К е Кт х т; О = diag(ц2,.., дп)> 0, К > 0. В соответствии с формализмом принципа максиму-
1 т т
ма [4] составляется гамильтониан Н = 2- (х Ох + и Ки) +
+ рТ(Гх + Ои), где р е Кп - сопряженный вектор, и записываются необходимые условия оптимальности
X = ^^ = Fx + Gu; op
• OH y?t „ P = --т- = -F p - Qx; ox
0 = O-H- = Ru + GTp. ou
(5)
(6)
(7)
Уравнения (5)-(7) являются уравнениями краевой задачи. Граничные условия для нее задаются соотношениями (2), (3) и выражениями
Чем проще модель, тем реже она обманет (принцип надежности) .
Четвертая заповедь математического моделирования
p(to) = ^o;
p( tf) = M'^f,
(8)
(9)
Рассматривается задача терминального управления в следующей максимально простой постановке. Пусть дана полностью управляемая линейная стационарная система
x = Fx + Gu; x e R , u e R ,
F e Rn
G e Rn
(1)
вытекающими из условий трансверсальности. В этих выражениях ^ е К и ^ е К - неопределенные множители Лагранжа, которые должны выбираться так, чтобы выполнялись условия (2), (3).
Задача формулируется следующим образом: построить аналитическое решение полученной линейной стационарной двухточечной краевой задачи.
2 РЕШЕНИЕ ЗАДАЧИ
или в развернутом виде
Из уравнения (7) найдем выражение для оптимального управления
-1 т и = -К О р,
(10)
подставим его в уравнение (5) и запишем полученное соотношение и уравнение (6) в виде следующей расширенной системы:
где
А =
г = Аг,
-1 т
Г -ОК О1
т
-О -г1
(11)
Г ,
; г = X . (12)
Р
Пусть ¿1, ¿2 е [¿0, ^] - два любых момента времени и Ь = ¿2 - ¿1 есть соответствующий им интервал. Тогда решение однородной линейной стационарной системы (11) может быть представлено в виде [5]
2 п 2 2 2 п - 1
ае1 (А) = (/)- ^(Ц.-.д;)(^) + + 54(4 ...,к2п)(З2) -2... ±Б2;(4 ...,к2п), (16)
22
где .--Дп), к = 1, ..., п - элементарные симмет-
рические функции, являющиеся суммами всех Сп к-членных произведений различных собственных значений, выбираемых среди 4 ..., 1?п [7]. С другой стороны, в матрице А для каждого к = 1, ..., 2п существует
различных главных миноров порядка к, суммы которых обозначаются через Ек(А) [7]. Напомним, что главным минором порядка к называется определитель подматрицы, расположенной на пересечении к строк и столбцов с одинаковыми множествами номеров. Вследствие гамильтоновости матрицы А все суммы Ек (А) с нечетными индексами оказываются равными нулю, так что ее определитель, выраженный через указанные суммы, принимает вид
г(¿2) = Ф(¿1)г(¿1),
(13)
Аг
где Ф(¿2, ¿1) = Ф(^2 - ¿1) = Ф(¿) = е - переходная матрица системы. Построение переходной матрицы является первым и наиболее трудным этапом решения краевой задачи. Второй этап, собственно решение этой задачи с учетом краевых условий (2), (3) и (8), (9), особых трудностей не представляет. Рассмотрим оба эти этапа последовательно.
Построение переходной матрицы. Точное выражение для переходной матрицы стационарной системы (11) дается формулой [5]
е = Ь 1 [(з1 - А) 1 ],
(14)
где символом Ь [.] обозначено обратное преобразование Лапласа. При практическом применении этого преобразования основной трудностью является разло-
1
жение знаменателя резольвенты (э1-А) матрицы А на элементарные множители, то есть решение полной проблемы собственных значений для матриц произ-
т^ л х-,2п х 2п
вольного вида. Гамильтонова матрица А е К имеет 2п собственных значений (с учетом их кратности), из которых всегда ровно половина имеют отрицательные вещественные части, а вторая половина - положительные, причем чисто мнимых собственных значений нет [6]. Следовательно, ее определитель может быть записан следующим образом:
ае1;(А) = (52- 4( з2-
■Л^ТЖ-Г), (15)
аеЦА) =
2 п 2 п 1 2 п 2
= (з1) -Е2(А)(з2) + Е4(А)(З2) ...±Е2п(А). (17)
Сравнивая в уравнениях (16) и (17) выражения при одинаковых степенях з , можем записать следующие соотношения:
Б2к (4-Дп) = Е2к(А), к = 1,..., п, (18)
существование которых свидетельствует о тесной взаимосвязи между элементарными симметрическими функциями собственных значений матрицы А и суммами ее главных миноров. Считая соотношения (18) алгебраическими уравнениями относительно неизвестных ..., Хп с правыми частями, зависящими от элементов матрицы А, мы приходим к одной из труднейших задач численного анализа - классической задаче на собственные значения в ее традиционной постановке. Имеется много подходов к решению этой задачи, однако ни один из них не является абсолютно надежным и универсальным. Вычисление собственных значений непосредственно из уравнения ае1( э1 -А) = 0 «...является одной из хорошо известных глупостей в численном анализе» [8]. Приведение матриц к диагональному виду не всегда возможно, так как далеко не всякая матрица является диагонализуемой. Процесс приведения матрицы к форме Жордана является и технически сложным, и исключительно неустойчивым. Малейшие изменения в матрице могут приводить к большим изменениям элементов жордановой матрицы и даже ее
структуры. Обеспечение численной устойчивости этого процесса - дело безнадежное [7]. Даже применение одного из самых эффективных алгоритмов вычислительной математики - ОК-алгоритма приведения матриц к форме Шура - далеко не всегда обеспечивает требуемую точность вычисления собственных значений и собственных векторов матриц общего вида, особенно при наличии у них кратных или близких по модулю собственных значений [6, 7].
Таким образом, в настоящее время не существует надежных способов преодоления самой первой и, по-видимому, самой главной трудности в решении краевой задачи - трудности решения полной проблемы собственных значений. Единственный корректный выход из сложившегося положения - вообще не решать данную проблему, а считать, как и в классической теории управления, все собственные значения заданными. Такое простое, на первый взгляд, решение меняет ситуацию самым радикальным образом. Все то, что для традиционного подхода было неразрешимой проблемой, с принятием нового подхода осуществляется легко и просто. Появляется возможность, например, осуществлять точное вычисление матриц преобразования к диагональной или жордановой форме, а также приведение матрицы А к форме Шура. Причем, во многих случаях необходимые для этого вычисления могут быть доведены до конца даже в символьном (!) виде. В этом и состоит основная идея рассматриваемого в работе подхода. Соотношения (18) при этом, разумеется, по-прежнему должны выполняться. Чтобы при новом подходе это стало возможным, в правых частях указанных соотношений должны существовать некоторые свободные параметры, выбором которых можно было бы обеспечить выполнение этих соотношений. В краевой задаче такие параметры существуют, и как раз в нужном количестве, - это суть п диагональных элементов ..., Яп матрицы О. Принимая это во внимание, преобразуем уравнения (18) к виду
Е2к(Чь .••' Яп) = .••' ^X (19)
где неизвестными являются теперь не параметры ^1,...Дп, а параметры Я1>.--> Яп. Эти параметры входят в левые части уравнений (19) в виде линейных комбинаций, так называемых, полилинейных функций, т. е. р-членных произведений различных величин, выбираемых среди Я1>.--> Яп. Здесь р = 1,..., т, где т - размерность вектора управления. В общем случае, уравнения (19), будучи записанными в форме
/к(Я)= Е2к(Я1 .• Яп) -32к(х1.--х2п) = 0
к = 1, ..., п, (20)
легко решаются относительно компонент вектора я = т
= (Я1,..., Яп) любым численным методом (например, итерационным методом Ньютона). Вычисление требующегося в методе Ньютона якобиана 5 f/Ъ я облегчается тем, что выражения для /(я) могут быть получены в символьном виде. Последующее дифференцирование этих функций также легко осуществляется в символьной среде МАТЬАБ. В частном случае, когда управление в задаче является скалярным (т = 1), уравнения (20) оказываются линейными, и их решение легко может быть получено в символьном виде.
Процесс построения переходной матрицы может осуществляться несколькими способами. Можно, например, представить детерминант матрицы А в форме (15) и выражения для элементов переходной матрицы получать непосредственно из уравнения (14). Символь-
-1
ное представление резольвенты (з1 - А) матрицы А в этом случае легко вычисляется с использованием алгоритма Леверье [5]. Кстати, получаемые в ходе реализации этого рекуррентного алгоритма выражения для коэффициентов а2к разложения детерминанта мат-2
рицы А по степеням 5 , совпадают с выражениями для коэффициентов Таким образом, применение алгоритма Леверье можно рассматривать в качестве еще одного способа формирования уравнений (20). Подобное непосредственное вычисление элементов переходной матрицы является, пожалуй, единственным способом, позволяющим, как говорится, «воочию» увидеть выражения для этих элементов в символьном виде. Однако на этом достоинства этого способа и заканчиваются. Указанные выражения в большинстве случаев оказываются такими невообразимо громоздкими, что никакого практического значения этот способ не имеет и применять его, вообще говоря, не следует (за исключением, может быть, самых простых случаев).
Гораздо более конструктивным является другой подход к построению переходной матрицы, основанный на методах теории подобных преобразований. Суть этого подхода состоит в представлении матрицы А в следующем факторизованном виде [5]:
А = ТУТ(21)
где т - матрица подобного преобразования (трансформирующая матрица), У - одна из матричных канонических форм. Выражение (14) для переходной матрицы принимает при этом следующий вид:
еА = ТЬ [(51 - У) ']Т . (22)
В зависимости от выбранного распределения корней и способа приведения матрицы А к виду (21), матрица У может оказаться диагональной, жордановой матрицей
или треугольной формой Шура. Смысл введения представления (22) заключается в том, что, по сравнению с матрицей А, матрица У имеет более простую структуру, благодаря чему обратное преобразование Лапласа в формуле (22) легко вычисляется. При этом выражение для переходной матрицы принимает вид
Уt -1
Ф(¿) = ТеУ Т .
(23)
Если в выбранном распределении Х^, ...Хп нет крат-
1
ных корней, символьные представления матриц т и т формируются в соответствии с алгоритмом [5]. Столбцами матрицы т являются 2п правых собственных векторов матрицы А, строки матрицы т формируются из соответствующих левых собственных векторов матрицы А. Поскольку собственные векторы всегда определяются неоднозначно, для полученных таким способом
1
матриц вместо условия Т Т = I выполняется обычно 1
лишь условие Т Т = В, где В - диагональная матрица. Поэтому на этапе численных расчетов вместо мат-1
риц Т и Т , вычисленных по их символьным представ-
1
" 2
лениям, должны использоваться матрицы Т = ТВ 1
"-1 2 -1
и Т = ВТ . Если среди корней ...,Хп имеются кратные, то формируется только символьное представ-
1
ление матрицы Т. Матрица Т вычисляется затем численным методом уже на этапе численных расчетов. При заданном распределении корней и неизменных па-
1
раметрах управляемой системы матрицы Т и Т являются постоянными. Таким образом, при вычислении переходной матрицы в соответствии с выражением (23) для каждого нового значения t заново приходится
вычислять только матрицу простой структуры е .
Как всякая переходная матрица, матрица Ф, определяемая уравнением (23), является симплектической, то есть удовлетворяет условию [6]
способе реализуются следующим образом. Вначале для заданных значений Х^, . ..,Хп и элементов матрицы К вычисляются значения элементов матрицы О, обеспечивающие сохранение инвариантов подобия, и формируется (в численном виде) матрица А. Затем эта матрица приводится к верхней форме Хессенберга Н. Наконец, с использованием ОК-алгоритма (со сдвигами) выполняется приведение матрицы Н к верхней форме Шура Б. В данном случае все сдвиги известны точно -это суть заданные нами собственные значения Х,1, ..., Хп, -Х1, ..., -Хп. Вследствие этого, теоретически бесконечный итерационный процесс приведения, характерный для стандартного ОК-алгоритма, в рассматриваемом случае превращается в конечный - приведение к форме Шура осуществляется всего за (2п - 1) шагов. Данная модификация ОК-алгоритма обладает многими привлекательными свойствами. Прежде всего, этот алгоритм одинаково хорошо работает при любом распределении корней и любых значениях элементов матрицы К. Кроме того, сам алгоритм является идеально обусловленным, т. к. матрицы преобразования к форме Хессенберга и форме Шура - всегда ортогональны. Наконец, так как сдвиги могут задаваться сразу в требуемом порядке, исключается операция переформирования трансформирующих матриц и самой формы Шура, обычно присутствующая в стандартном ОК-алгоритме. Недостаток у этого алгоритма всего один, но довольно существенный. При любом изменении параметров модели, значений корней распределения или элементов матрицы К вычисления с использованием этого довольно громоздкого алгоритма должны выполняться заново, начиная с вычисления новых значений элементов матрицы О, формирования матрицы А и предварительного приведения ее к форме Хессенберга. В связи с этим, применение формы Шура в ряде случаев может оказаться затруднительным.
Решение краевой задачи. Разобьем переходную матрицу на блоки
Ф(¿2, ¿1) =
Ф'/Ф=/,
(24)
Ф11(^ ¿1) Ф12(^ ¿1) Ф21 (^ ¿1) Ф22(^ ¿1)
где ] =
0 I -I 0
матричный аналог мнимой единицы.
Указанное свойство значительно облегчает обращение матрицы Ф, а равенство (24) может служить хорошей проверкой правильности ее вычисления.
Хотя в данной статье рассматриваются, в основном, способы аналитического конструирования переходной матрицы, следует отметить, что существует также способ численного решения указанной задачи, основанный на использовании в качестве канонической формы У треугольной формы Шура Б [6, 7]. Вычисления в этом
и запишем соотношение (13) в виде двух уравнений
х(¿2) = Фц(¿2, ¿1)х(11) + Ф12(¿2, ¿1)р(¿1); Р(¿2) = Ф21 (¿2, ¿1)X(11) + Ф22(¿2, ¿1 )р(¿1) .
Эти уравнения справедливы при любых значениях ¿1, ¿2 е [¿0, ¿у]. В частности, полагая в первом из них ¿2 = ¿0, ¿1 = , можем записать
X(¿0) = Фц(¿0, ух(У + Ф12(¿0, уР(у). (25)
Полагая во втором уравнении = Ь/, ^ = получим
р(/) = Ф21 (/ х(^) + Ф22(/ Ур(У. (26)
Для сокращения записи введем обозначения
Ф(/ ^) = и
Ф(/) = V
и11 и12 и21 и22
V11 V12 V21 ^^22
(27)
(28)
Симплектичность матриц и и V позволяет обойтись вычислением по уравнению (23) только одной (любой) из них. Если, например, вычислена матрица V, то обратная ей матрица и может быть получена из соотношений
и = V =
ТТ
ТТ -^1 VT1
(29)
Выразим х(Ь/) из уравнения (25) и, с учетом граничных условий (2), (9), запишем
1 Т
х(/) = (С0- Vl2Mí 4/).
Умножим это уравнение на матрицу М слева
1 1 Т
Мх( /) = MV11c0-MV11V12Mi^f
(30)
и введем обозначение
1 Т
Н = MV11V12Mi .
Тогда из уравнения (30), с учетом граничного условия (3), получим
1 1
4/ = Н 1 (MVl1Co- С/).
(31)
Перепишем уравнение (26), с учетом (8),(9), в виде
Mт 4/ = и21С0 + и2240.
Отсюда с использованием соотношений (28), (29) получим окончательное выражение для множителя Ла-1 т
гранжа 40 = и22(M 4/- и2С) или, с учетом (27)-(29), (31),
-Т Т -1 -1 Т Т -1
40 = V1\[(MTH 1MV11 + VT1) С0-M1H 1/ (32)
Выражения для вычисления в момент Ь е [ /] программных значений переменных краевой задачи принимают при этом следующий вид:
х (Ь) = Фц(Ь, ^)С0 + Ф12(Ь,
*
р (= Ф21 (Ь0)С0 + Ф22(и Ь0)40.
(33)
Как правило, исходные физические уравнения управляемой системы являются нелинейными, их вид определяется соответствующими законами природы. Опыт успешных применений методов оптимизации к решению практических задач [10, 11, 12] показывает, что для решения задач оптимизации использовать эти уравнения в таком виде никогда не следует. Для постановки оптимизационной задачи на основе исходных уравнений всегда должна разрабатываться более простая линейная стационарная модель управляемой системы. Если управление и(Ь) в этой линейной модели соответствует физическому управлению в исходной модели, то уравнение
* 1 Т *
и (Ь) = -К 1 о'р (Ь)
(34)
представляет собой искомое аналитическое решение краевой задачи в форме программного управления. Если положить в приведенных выше формулах ^ равным текущему моменту времени Ь, то решение краевой задачи, с учетом (32), может быть получено в форме закона управления с обратной связью
и (Ь) = -К( Ь) х (Ь) + К/( Ь) С/,
1 Т Т Т 1 1 Т где К (Ь) = К 1 G1V1¡(M1H 1MV11 + VT1), К/( Ь) =
-1 Т -Т Т -1 = К 1 G1V11iM1H .
Размерность линейной модели не обязательно должна совпадать с размерностью исходной системы. Если используется расширенная модель системы, управление и(Ь) в которой не связано с физическим управлением в исходной задаче, решение задачи программного управления определяется соотношением
и (Ь) = хк(,
(35)
где к - номер координаты в линейной системе, связанной с физическим управлением в исходной задаче.
Для компенсации внешних возмущений и погрешностей модели суммарное управление в линейной системе должно формироваться в виде
и (Ь) = и*( Ь) + 5 и (Ь),
(36)
где 5и(Ь) - некоторая поправка, обеспечивающая стабилизацию движения объекта относительно программной траектории х*(Ь). Вычисление величины этой
поправки может осуществляться как методами классической теории регулирования, так и с использованием базового алгоритма.
Анализ обусловленности задачи. Решение краевой задачи происходит в два этапа и на каждом из них существует опасность встретиться с плохо обусловленными вычислениями. Из сказанного выше, уже ясно, что для реализации вычислений на первом этапе - этапе построения переходной матрицы - не существует алгоритма, наилучшего во всех отношениях. Вычисления с использованием формы Шура всегда идеально обусловлены, и в этом смысле их применение является весьма желательным. Однако эти вычисления очень громоздки и во многих случаях попытки их реализации в адаптивных системах в реальном масштабе времени могут натолкнуться на серьезные трудности. Напротив, вычисления по аналитическим формулам очень просты и, в этом смысле, идеально вписываются в концепцию построения адаптивных систем. Однако, если не предпринимать никаких специальных мер, эти же вычисления оказываются, как правило, и чрезвычайно плохо обусловленными. Связано это с тем, что трансформирующие матрицы преобразования (21) формируются из собственных векторов, а плохая обусловленность вычислений с собственными векторами - дело обычное [6].
Как показывают результаты моделирования, характер процессов в задаче терминального управления слабо зависит от выбранного распределения корней и элементов матрицы К, однако сам процесс вычисления переходной матрицы при неудачном выборе этих параметров может оказаться чрезвычайно плохо обусловленным. Поэтому начинать нужно, прежде всего, с выбора таких значений указанных параметров, при которых как сама матрица А краевой задачи, так и трансформирующие матрицы преобразования (21) имели бы как можно меньшее число обусловленности (напомним, что число обусловленности - это тот коэффициент, на который умножаются ошибки в начальных данных при формировании конечного результата). Использование этого подхода позволяет в некоторых случаях добиться даже идеальной обусловленности. Например, если матрица А принадлежит к классу так называемых нормальных матриц, то есть, если для нее выполняется условие ТТ
(АА = А А), то, независимо от наличия у нее кратных собственных значений, она всегда диагонализуема, причем диагонализация всегда осуществляется унитарным преобразованием [7]. Таким образом, нормальная матрица А - это тот идеал, к которому следует всячески стремиться в ходе решения краевой задачи. В некоторых случаях этого идеала действительно можно достичь путем рационального выбора распределения корней характеристического уравнения и элементов матрицы К. Например, если система (1) представляет собой «-мерный интегратор, наилучшим оказывается
выбор К = 1 и распределение корней по типу распределения Баттерворта [5]
Хк = гехр[(2к=-1 + З) Хп + к =
к = 1, ..., п (37)
с радиусом круга распределения г = 1. При таком выборе все диагональные элементы матрицы О, кроме становятся равными нулю, ^ оказывается равным 1, матрица А переходит в класс нормальных, а матрица преобразования Т превращается в унитарную с числом обусловленности, равным 1.
Если бы мы выбрали другое распределение корней, например, биномиальное, число обусловленности матрицы Т в рассматриваемом примере составило бы величину порядка 5 108. Это свидетельствует о том, что обусловленность вычислений в краевой задаче сильно зависит от распределения корней и элементов матрицы К. Поэтому вопросам выбора этих параметров необходимо уделять самое серьезное внимание и не жалеть для этого никаких усилий. Даже если всегда не удается найти такое идеальное их сочетание, как в рассмотренном выше примере, всегда можно выбрать указанные параметры таким образом, чтобы существенно улучшить обусловленность вычислений по сравнению с их обусловленностью при каком-либо случайном выборе данных параметров.
Необходимо стремиться также к хорошей обусловленности еще одной матрицы, входящей в разложение
(21), а именно, матрицы еУ1:. Ее сингулярное число
^ 2at
обусловленности можно оценить величиной е , где а = тах[аЬ8(Ие(Хк))], к =1, ..., п, (для диагональных матриц эта оценка является точной). Величина а должна выбираться в зависимости от величины интервала времени ¿ таким образом, чтобы обеспечивалось согласование чисел обусловленности матрицы е^ и трансформирующей матрицы Т.
Таким образом, требование наилучшей обусловленности вычислений накладывает довольно жесткие ограничения на свободу выбора распределения корней и элементов матрицы К. Можно сказать даже, что указанным требованием эти параметры определяются почти однозначно. Отсюда следует, что в рассматриваемом подходе матрицы К и О более не являются свободными параметрами и им невозможно придать тот физический смысл, который обычно приписывается этим параметрам в классической теории оптимальных систем. Теперь К - это просто некоторая матрица коэффициентов, совместно с корнями характеристического уравнения обеспечивающая наилучшую численную обусловленность процесса построения переходной матрицы, а О - зависящая от них матрица, отвечающая за
сохранение инвариантов подобия. Некоторая свобода выбора элементов матрицы К появляется лишь в тех случаях, когда имеющейся разрядности оказывается достаточно для проведения вычислений и можно особо не заботиться о проблеме обусловленности.
Подводя итог, можно сказать, что при использовании на первом этапе вычислений по аналитическим формулам, за счет выбора распределения корней и элементов матрицы К всегда можно (если не радикально, то, по крайней мере, существенно) улучшить обусловленность вычислений. В критических случаях можно попытаться использовать описанную выше модификацию ОК-алгоритма. Другими словами, обусловленность вычислений на первом этапе является в той или иной мере управляемой.
Совсем иная ситуация складывается на втором этапе решения краевой задачи - этапе формирования значений адаптивных управлений. В классической и современной теории управления для обеспечения устойчивости движения в качестве полюсов замкнутой системы всегда выбираются только те п характеристических чисел, которые расположены в левой полуплоскости. Для рассматриваемого в данной статье подхода характерно, что в решении краевой задачи участвуют все 2п характеристических чисел матрицы А расширенной системы (11), половина из которых расположена в левой полуплоскости, а вторая половина является их зеркальным отражением относительно мнимой оси, т. е. расположена в правой полуплоскости. Это означает что как матрица А, так и связанная с ней краевая задача являются неустойчивыми. К тому же, в силу специфики самой задачи, в ходе ее решения приходится оперировать не с самой переходной матрицей (как единым объектом), а с ее блоками. Даже в том случае, когда выбором распределения корней и элементов матрицы К удается добиться идеальной обусловленности самой переходной матрицы, ее внедиагональные блоки могут оказаться крайне плохо обусловленными. Это связано с тем, что характеристические числа матрицы А в выражения для вычисления элементов этих блоков всегда входят парами - наряду с числом Хк используется и число с противоположным знаком -Хк. Это приводит к тому, что выражения для вычисления элементов указанных матриц оказываются знакопеременными, вследствие чего итоговые значения этих элементов вычисляются всегда как малые разности близких (но отнюдь не малых) чисел.
Вот типичный пример подобной ситуации, когда система (1) является «-мерным интегратором (п = 8). Для этой системы при Ь = 0, 025 с и распределении корней в соответствии с (37) число обусловленности переходной матрицы и в уравнении (27) равно 1,05. Если бы вычисления проводились с самой этой матрицей, столь близкое к единице число обусловленности, без сомнения, свидетельствовало бы об очень хо-
рошей обусловленности вычислений. Однако, как мы уже убедились, в решении краевой задачи участвует не сама матрица и, а ее блоки. И если диагональные блоки и 11 и и22 обусловлены так же хорошо, как и сама матрица и, то обусловленность внедиагональных блоков в данном случае оказывается просто катастрофически плохой. Число обусловленности блока и^, на-
30
пример, характеризуется величиной 1, 25 • 10 . Вдумайтесь только - это число более чем в миллиард миллиардов миллиардов (!) раз больше числа обусловленности матрицы и. Причину столь плохой обусловленности можно понять, если детально рассмотреть процесс вычисления какого-либо элемента этой матрицы, например, элемента и^( 1, 1). В выражениях для его вычисления положительные и отрицательные слагаемые чередуются самым регулярным образом. Вычисленная средствами МАТЬАБ сумма положительных слагаемых равна
0,3203292768698897028728241250754498983474819975, сумма отрицательных
-0,3203292768698897028728241250754798976352843760, точное решение
0,0000000000000000000000000000000000007121976215.
Легко видеть, что, если проводить вычисления с разрядностью, меньшей 36, в ответе не окажется вообще ни одного верного знака.
Во многих случаях ситуация оказывается не столь катастрофической, как в приведенном примере, однако общая тенденция всегда сохраняется. Поэтому нет ничего удивительного в том, что попытки решения краевых задач традиционными численными методами (да еще в рамках стандартной разрядности) очень часто оканчиваются неудачей. В таких случаях обычно отмечается, что задача обладает слишком высокой чувствительностью конечных значений фазовых переменных к погрешностям вычисления начальных значений сопряженных переменных [9]. Как видим, причина этой высокой чувствительности крайне проста и заключается она в недостаточной разрядности проводимых вычислений.
Для линейных краевых задач самым радикальным и единственным способом преодоления указанной трудности является проведение вычислений с повышенной разрядностью. Для нелинейных задач, по мнению автора, вряд ли помогут вычисления с любой, даже самой высокой, разрядностью, так как к этим задачам неприменим краеугольный принцип теории линейных систем - принцип суперпозиции.
Из приведенного примера видно, что требуемое повышение разрядности может иногда оказаться довольно значительным, особенно для систем высокого порядка. В этом смысле, образно можно сказать, что над методами оптимизации, основанными на принципе максимума, тяготеет некое «проклятие», которое, по аналогии со знаменитым беллмановским «проклятием размерности» в динамическом программировании, можно было бы назвать «проклятием разрядности». К счастью, кажется, преодолеть это второе «проклятие» намного легче, чем первое. Для этого достаточно разработать новые версии компиляторов, допускающих вычисления с повышенной разрядностью. Задача эта не представляется такой уж неразрешимой. По крайней мере, один успешный опыт такой разработки известен - созданный много лет назад компилятор с языка PL/1 для ЕС ЭВМ позволял проводить вычисления с разрядностью мантиссы, равной 33. Стандартной разрядностью повсеместно распространенных в настоящее время IBM PC является разрядность 16, и возможностей ее увеличения в рамках существующих систем программирования компилятивного типа, к сожалению, не существует.
Базовый алгоритм. Данный алгоритм не является вычислительным алгоритмом в обычном понимании этого слова. Наряду с чисто вычислительными этапами он содержит этапы принятия решений, которые не могут быть выполнены без участия человека. Кроме того, алгоритм самым существенным образом опирается на возможность проведения символьных преобразований, которые для получения достаточно простых конечных выражений необходимо проводить в интерактивном режиме, т. е. также с участием человека. Таким образом, можно сказать, что базовый алгоритм - это последовательность действий человека, проектирующего адаптивную систему для своего конкретного приложения и находящегося с этой целью в некоторой среде символьного программирования.
Алгоритм формулируется в виде следующей последовательности шагов.
1. На основе исходных (как правило, нелинейных нестационарных) уравнений, вытекающих из законов природы, строится линейная стационарная модель управляемого процесса. Для этого могут использоваться различные методы, среди которых особого внимания заслуживает подход, основанный на идее [10] о возможности точного представлении уравнений исходной задачи в виде нескольких независимых 2-мерных интеграторов. Эта идея допускает естественное обобщение на случай интеграторов произвольного числа измерений. При этом вектор состояния может расширяться как «вверх» (введением интегралов - например, для придания системе свойств астатизма), так и «вниз» (введением высших производных). Подобное расширение во многих случаях оказывается весьма эффек-
тивным средством придания управляемому процессу новых свойств, которые невозможно получить при равенстве размеров исходной и оптимизационной моделей.
Часто в качестве оптимизационной модели может быть принято даже довольно грубое приближение к уравнениям исходной задачи. При этом параметры этой упрощенной модели необходимо выбирать таким образом, чтобы в ходе управления с обратной связью влияние сделанных допущений уменьшалось по мере приближения управляемого объекта к заданному терминальному состоянию Cf. Удачные примеры такого подхода содержатся в [11, 12].
Важным является также требование, чтобы оптимизационная задача с использованием построенной модели не содержала ограничений на управляющие функции и, тем более, на фазовые координаты. Только в этом случае краевая задача, оставаясь линейной, может быть решена предлагаемым в данной работе способом. Если в исходной задаче ограничения на управления все-таки присутствуют, они могут быть учтены a posteriori в ходе управления с обратной связью.
2. Для построенной линейной модели в соответствии с выражением (12) формируется символьное представление матрицы A. При этом для параметров модели, относительно которых предполагается, что они могут изменяться, вводятся краткие символьные обозначения.
3. Принимается решение о характере распределения собственных значений ..., Xn матрицы A: являются ли все они различными или же среди них имеется определенное количество кратных собственных значений.
4. Для принятого характера распределения собственных значений определяется структура канонической формы Y, к которой преобразованием подобия (21) приводится матрица A. Самый простой случай, когда все собственные значения - различны. Каноническая форма в этом случае является диагональной матрицей с элементами, равными этим собственным значениям, а матрицы преобразования T и T 1 составляются из правых и левых собственных векторов, соответствующих этим собственным значениям.
Если в выбранном распределении присутствуют кратные собственные значения, в число столбцов и строк матриц преобразования входят соответствующие этим значениям обобщенные собственные векторы. Канонической формой в этом случае является жорданова матрица, структура которой определяется числом и размерами жордановых блоков. Для того чтобы определить эти параметры, необходимо для каждого различного Хр
определить степени (XpI-A)k, где k = 0, 1, ..., n, проанализировать последовательность рангов этих степеней и по результатам анализа определить размеры и число жордановых клеток, отвечающих данному собственному значению [7].
5. В соответствии с алгоритмом [5] формируется символьное представление матрицы преобразования Т. Ее столбцами являются 2п правых собственных векторов матрицы А (включая и обобщенные). Если в выбранном распределении нет кратных корней, символьное представление строк матрицы Т 1 формируется из соответствующих левых собственных векторов матрицы А. Для этого указанный алгоритм применяется к матрице АТ.
6. С использованием обратного преобразования Лапласа вычисляется символьное представление матрицы
УЬ -1 -1
е = Ь 1 [(^1 -У) 1 ].
7. Формируются уравнения для определения элементов д^ • ••, Яп. Получить эти уравнения можно одним из следующих способов:
- приравниванием коэффициентов при одинаковых
2
степенях 5 в выражениях для детерминантов матриц (51-А) и (51-У);
- с использованием инвариантов подобия Е2к в уравнении (20);
- с использованием коэффициентов а2к, получаемых в ходе реализации рекуррентного алгоритма Ле-верье [5];
- путем вычисления (2пх2п) - матрицы М = ТУ -- УТ с последующим приравниванием нулю тех п из 2п ее ненулевых элементов, которые образуют систему независимых уравнений.
8. В случае, если полученные уравнения оказываются линейными относительно они решаются в символьном виде. В результате находятся символьные выражения для д^ •.., дп. В противном случае в символьном виде формируются система нелинейных уравнений
/(д) = 0, где д = [д^ •.., дп] , и якобиан этой системы
_5£
Символьные выражения для элементов якобиана
находятся с использованием соответствующих символьных средств системы МАТЬАБ.
9. В соответствии с выбранным распределением корней задаются их численные значения
итераций значения дк должны быть неотрицательными. Добиться этого можно выбором распределения корней и/или значений элементов матрицы К. В необходимых случаях пункты 9-11 приходится повторить несколько раз.
12. По символьным представлениям трансформи-
1
рующих матриц Т и Т определяются их численные
1
значения, находится диагональная матрица В = Т Т,
1 1 " 2 " -1 2 -1 и вычисляются матрицы Т = ТВ и Т = В Т ,
" -1 -
удовлетворяющие условию Т Т = I. Если среди кор-
1
ней распределения имеются кратные, матрица Т вычисляется как обратная матрице Т. В этом случае указанная выше балансировка строк и столбцов не производится.
13. Задается величина интервала Ь = ^ - Г и по сим-
УЬ
вольному представлению матрицы е вычисляется ее численное значение.
14. По формуле (23) вычисляется переходная матрица V = Ф(и по формуле (24) осуществляется проверка правильности ее вычисления.
15. Задаются граничные условия (2), (3) краевой задачи и по формуле (32) вычисляется начальное значение ^0 сопряженного вектора р.
16. Для требуемого момента времени Ь по формулам (33) вычисляются значения параметров х*(Ь) и р*(Ь) программной траектории и по формулам (34) или (35) -значение программного управления м*(Ь).
17. С использованием соответствующих алгоритмов формируется величина поправки 5м(Ь) и по формуле (36) вычисляется управление м(Ь) в линейной задаче.
18. При необходимости, по управлению м(Ь) формируется физическое управление т( Ь), подаваемое на исполнительные органы исходной системы.
Возможности данного алгоритма ограничиваются только возможностями существующих пакетов символьного программирования.
Хк = ак + V; Хп+к = -Хк;к =1,;п.
10. Задаются численные значения элементов матрицы К.
11. Вычисляются значения параметров д^-., дп по формулам п. 8 или с использованием итерационного алгоритма Ньютона
др+ 1 = д
/(др), р = 0,1,;.
д = др
В ходе вычислений необходимо контролировать знаки параметров дк; все полученные после окончания
ВЫВОДЫ
В работе предложен конструктивный подход к решению проблемы синтеза адаптивных линейных стационарных систем без ограничений на фазовые переменные и управляющие функции. В его основу положена идея о замене традиционного подхода с априорным заданием матриц О и К критерия качества новым подходом, в котором априорно задаются полюса замкнутой системы. Сохранение инвариантов подобия матрицы А краевой задачи при этом обеспечивается соответствующим выбором диагональных элементов матрицы О критерия качества.
В классической теории управления, в основе которой лежит тот же принцип априорного выбора полюсов, основным критерием для подобного выбора является устойчивость движения управляемой системы. Методам, основанным на принципе максимума, внутренне присущи высокие требования к разрядности вычислений («проклятие разрядности»). Поэтому в предлагаемом подходе аналогичный критерий должен основываться на понятии обусловленности вычислений.
ПЕРЕЧЕНЬ ССЫЛОК
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. Дуайер Т. А. У, 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
с использованием других подходов. Эти примеры охватывают круг задач терминального управления, оптимального регулирования и наблюдения, управления с эталонной моделью, управления с ограничениями в виде неравенств, наложенными на управляющие функции.
Остальные примеры, взятые из близкой автору области управления в космосе, иллюстрируют возможности применения базового алгоритма к решению задач регулирования и стабилизации, наведения, а также обеспечения плавной стыковки участков траектории в случае, когда общая задача управления разбивается
