Особенности численного решения матричного алгебраического уравнения Риккати методом установления Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ Ц А Г И

Т о м X 19 7 9 № 1

УДК 62-501.12

ОСОБЕННОСТИ ЧИСЛЕННОГО РЕШЕНИЯ МАТРИЧНОГО АЛГЕБРАИЧЕСКОГО УРАВНЕНИЯ РИККАТИ МЕТОДОМ УСТАНОВЛЕНИЯ

В. М. Кувшинов

Рассматривается условие сходимости и скорость сходимости итерационного процесса решения алгебраического уравнения Рик-кати методом установления. Исследуются особенности численного решения, возникающие при нарушении условия сходимости. Проводится сравнение эффективности применения в рассматриваемом случае методов Эйлера и Рунге—Кутта численного интегрирования дифференциальных уравнений.

При решении ряда задач синтеза оптимальных систем автоматического управления летательных аппаратов возникает необходимость нахождения положительно определенного решения матричного алгебраического уравнения Риккати

~ А^Р-РА + РВР~^В^Р-^ = 0. (1)

Уравнение (1) появляется в задачах синтеза оптимальных регуляторов (матрица Р определяет коэффициенты обратной связи) и в задаче оптимальной фильтрации в квазистационарном случае (Р — ковариационная матрица вектора фазового состояния).

Методам решения этого уравнения посвящено большое количество работ, например, [1-—4]. В ряде работ [3 — 5] для его решения предлагается использовать метод установления, основанный на том, что решение матричного дифференциального уравнения Риккати

-А'РМ — РЮА+РЮВК-'В'Р^-ЦГ^ — РУ) (2)

с начальным условием Р(0)=э0 сходится к положительно определенной матрице Р, удовлетворяющей матричному алгебраическому

уравнению Риккати, т. е. при />(0) = 0, НтЯ(0 = Р. Этот факт был

*-►00

установлен Калманом [6].

Однако вопрос о рациональном выборе метода и шага численного интегрирования уравнения (2) при практическом использовании данного алгоритма на ЭЦВМ остается открытым. Для того, чтобы

иметь возможность произвести такой выбор, необходимо исследовать особенности поведения численного решения уравнения (2), возникающие из-за конечности шага интегрирования.

Основные особенности метода установления решения уравнения Риккати могут быть выявлены при анализе численного решения уравнения Риккати первого порядка

л: = — Ь2 х2 + %ах + Я2> х(0) = х0, (3)

которое можно также представить в виде

х = — Ь‘2(х — г^){х — г2),

в"Ь 10910 13 ГЛ ТУ

где 21,2 = —----р-----—корни уравнения Ь2 х* -}- 2ах + ц1 = 0. Как

показывают численные расчеты, эти особенности сохраняются и для уравнений более высокого порядка.

Рассмотрим численное решение уравнения (3) методом Эйлера:

*/+1 = х1,х, И или Х1+1~Х1 +{?2 + 2ах1 — №х]) к, (4)

где Л — шаг интегрирования.

Процесс численного интегрирования (4) можно представить на фазовой плоскости (х, х) так, как показано на фиг. 1, на которой приведен пример решения, сходящегося к искомому установившемуся значению гх. Можно построить на фазовой плоскости пример решения, расходящегося после второго шага интегрирования при

и \ + ^г\ ~ 4г1 г2 т-

----------^2-------• Если пользоваться данной графической интерпретацией, а также примерами расчетов, то можно установить, что возможны и другие типы решений при различных значениях шага интегрирования. Так, на фиг. 2 приведен пример решения, выходящего на режим установившихся колебаний. Рассмотрим этот тип решения.

Как видно из построения на фиг. 2, установившиеся колебания возможны при выполнении следующих условий:

АВ = АР, ВС — АВ-к. (5)

Пусть Аха= хА—гх, Ах0=х0—ги где хА=хв и х0 —хЕ— значения, которые поочередно принимает величина х при установившихся колебаниях. Очевидно,

АВ = — й2 &хА (Аха + к), АР = — Ь2 Ахв (Ахв + к), ВС = Ах0 — Ахл,

где А = — г2,

и условия (5) примут вид

Ах])-^- Ах\ + + Ал;л)& = 0; |

Ахв — Аха — — Ь2 Аха (Длл + к) А. | ^

Исключив из этой системы Ахв и отбрасывая решения Аха = 0 и Дхл = — А, не представляющие интереса для анализа установившихся колебаний, получим квадратное уравнение для определения

Дх,

Ах\+(к-^.)АхА-

&

= 0.

Действительные корни этого уравнения при ^>-^-дают решения системы (6)

2

Ад:,

4

Л2 62

- Л

Л62

4

лг г>2

(7)

Таким образом, при интегрировании уравнения Риккати (3) мето-

2

дом Эйлера с шагом Л>ЛА, где /гА = — критический размер

шага, могут иметь место установившиеся колебания решения с амплитудой, определяемой соотношениями (7), которые после введения обозначения га = Л/й* можно переписать следующим образом:

.(8)

Однако реализация таких колебаний зависит от устойчивости этого-режима численного решения.

Если при установившихся колебаниях величина х поочередно принимает значение хв и хв (см. фиг. 2), то должны выполняться следующие условия:

ХЕ “ ХВ + (?* + 2аХВ — 1)2 Х%) к>

хв> =ХЕ + (я2 + 2 ахЕ — Ь2 х|) /г;

(9)

Для исследования устойчивости колебательного численного решения отклоним точку В от режима колебаний на некоторую малую величину Ьхв. Обозначим новое значение абсциссы точки В хв + 4- Ьхв, а точек Ей В' — хЕ + ЪхЕ и хв,-}-Ьхв' соответственно. Очевидно, что для устойчивости колебательного решения необходимо

<1. Найдем Ьхв, из системы уравнений

выполнение условия

Ьхв,

Ьхв

ХЕ + ЬХЕ=ХВ + Ьхв + [Я2 + 2« (хв + ЬхВ) — Ь% (ХВ + 8*в)2]

ХВ' + ЬхВ’ = ХЕ + ЬхЕ + № + 2а (ХЕ + 3*я) ~ 1,2 (ХЕ + Н-

Учитывая условия (9) и линеаризуя эти уравнения по 8л;в и ЬхЕ, получим следующее выражение для Ьхв,:

Ьхв, = Ъхв [1 -|- (2а — 2Ьъхв)к\ [1 4- (2а — 2Ь2 хЕ)Ь].

Тогда условие устойчивости колебательного решения в малом примет вид

|[1 +{2а-2Ь2хв)Н\[\ (2а - 26* дся) А] | < 1. (10)

Преобразовав неравенство (10) к переменной /г. = Л/Л* с учетом соотношений (7), получим условие устойчивости колебательного решения в малом

|5 — 4д!|<1 или 1<«<>/1,5.

Следовательно, при «>У^ 1,5 установившиеся колебания реализоваться не могут. Можно показать, что при «> 1,5 решение является расходящимся. В интервале же 1/1,5<«< 1,5, как видно из примеров расчетов, имеют место неустановившиеся колебания.

Заметим, что искомое установившееся решение является частным случаем колебательного решения при хв = г1 и хЕ = ги Подставляя эти значения хв, хЕ в неравенство (10), получим условие устойчивости установившегося решения или, что то же самое, условие сходимости процесса численного решения

(1 — 2л)2 < 1 или 0<«<1-Области различного поведения численного решения уравнения Риккати первого порядка (3), полученного методом Эйлера, при различных значениях шага интегрирования приведены на фиг. 3.

Результаты аналитического исследования подтверждаются примерами численного интегрирования уравнений Риккати первого порядка. На фиг. 4 приведены результаты интегрирования уравнения

* = - 0.28ЛГ2 — 0,56л; + 6,72, (11)

для которого Ай = 0,7143, при л = 0,1; «=0,2; п—0,3; « = 0,5; «=0,9; л = 1,1, «=1,3. При рассмотрении интегрирования уравнения Риккати (3) методом Рунге — Кутта начнем с примеров численного

интегрирования. Расчеты проводились при вариации параметра п = Щкк. Некоторые характерные результаты интегрирования уравнения (11) приведены на фиг. 5. Для значений п= 1,4ч- 1,7 проявляется новая особенность поведения численного решения, а именно, численное решение стремится к установившемуся значению которое не удовлетворяет алгебраическому уравнению Риккати.

Для метода Рунге — Кутта Для метода Эйлера

Р1 — область сходимости; Э1 — область сходимости;

Р2 — область фиктивной сходимости; Э2 — область установившихся колебаний;

РЗ— область установившихся колебаний; ЭЗ — область неустановившихся колебаний;

Р4 — область неустановившихся колебаний; Э4 — область расходимости

Р5 — область расходимости

Фиг. 3

Численное решение уравнения (1.7) методом Эйлера

Фиг. 4

Численное решение уравнения (1.7) методом Рунге — Кутта

Фиг. 5

При одинаковых значениях п решение имеет один и тот же характер для уравнений Риккати с различными значениями коэффициентов. Покажем, что параметр п является в определенном смысле параметром подобия численных решений уравнений Риккати. Сделаем в уравнении (3)замену переменных

АХ'

k ’

Тогда уравнение (3) примет вид

dy

dx

*0

к

(12)

Численное интегрирование уравнения (3) в реальном времени t с шагом h эквивалентно интегрированию уравнения (12) в безразмерном времени х с шагрм =/г/ЛА = га. Таким образом, роль шага интегрирования h в безразмерном времени играет параметр га.

Уравнение (12) не содержит никаких параметров исходного уравнения (3) и зависимость точного или численного его решения от параметров уравнения (3) входит только через начальные условия. Условия же сходимости, колебательности, устойчивости колебаний, фиктивного установления не включает в себя начальных условий, а содержат только параметры рассматриваемого уравнения и размер шага интегрирования. Значит, для уравнения (12) аналитически (как для метода Эйлера) или с помощью численного эксперимента (для метода Рунге — Кутта или других методов) можно определить функции yf(n), yD{n), уА(п), удовлетворяющие таким условиям, как условие фиктивного установления, условие колебательности решения, т. е. условиям в виде равенств, а также получить области параметров га, удовлетворяющих таким условиям,^ как условие сходимости, условие реализации режима фиктивного установления, условие устойчивости режима установившихся колебаний, т. е. условиям в виде неравенств. Для метода Эйлера роль функций AxDlk=y3D(n), AxA/k=y3A(n) играют зависимости (8) (см. фиг. 3), а границы реализации различных режимов nk = 1, гаср = =1/1,5, «р=1,5. Для метода Рунге — Кутта по результатам проведенных расчетов можно построить „экспериментальные" зависимости Axf/k=yj(n), ДдГд/А («\ AxAlk=y^(n) и определить границы и области различных режимов. Эти зависимости и границы (см. фиг. 3) будут иметь универсальный характер для уравнений Риккати с различными значениями коэффициентов.

Как показывают примеры расчетов, для уравнений Риккати, имеющих порядок TV > 1, характер поведения численного решения, получаемого методами Эйлера и Рунге — Кутта при различных значениях шага интегрирования, тот же, что и для уравнения первого порядка. При интегрировании уравнения Риккати (2)третьего порядка с матрицами

'-2,66 -1,57 —24,3 ' '-52,6 —16,3 5,55 "

А = —0,09 -0,66 —14,4 ; В = 1,79 - 7,52 3,82

0,042 1 - 0,318 0 — 0,056 -0,026

'10 0 0 ■ '10 0 0‘

W= 0 1 0 ; я = 0 4 0

0 0 100 0 0 10

методом Эйлера в диапазоне 0< А <0,016 имеет место сходимость, в диапазоне 0,017 < А <0,020 — установившиеся, а в диапазоне 0,021 <А<0,025— неустановившиеся колебания, и при А>0,026-расходимость решения. При интегрировании того же уравнения методом Рунге — Кутта при 0<А<0,023 решение сходится к искомому значению, при 0,024 < Л < 0,029 наблюдается фиктивное установление, при 0,030<Л< 0,031 — установившиеся, а при 0,032< <А<0,035 —неустановившиеся колебания, и при А>-0,036 решение расходится. При интегрировании уравнения Риккати третьего порядка с матрицами

'-1 -0,001 -0,001* '—52,6 -16,3 5,55 "

А - -0,008 — 1,167 -0,28 ; в = 1,79 - 7,52 3,82

0,031 0,72 —0,481 0 — 0,056 -0,026

'0,131 3,06 2,22 ' ' 2767 856 —29 Г

3,06 72,0 51,8 ; Я = 856 271 - 93

2,22 51,8 39,7 — 291 - 93 32

методом Эйлера при 0<Л<0,11 решение сходится, при 0,12<Л< <0,14 реализуется режим установившихся колебаний, при 0,15 < < А < 0,17 — режим неустановившихся колебаний, и при Л>-0,18 решение расходится. Отметим, что в этом случае все собственные числа X* матрицы оптимальной замкнутой системы А = А ■—Р действительные, а критическое значение шага оказывается

I

Наличие режимов фиктивного установления и установившихся колебаний представляет не только теоретический интерес, но имеет и практическое значение, точнее говоря, таит в себе определенную опасность при практическом использовании метода установления. Так, если шаг интегрирования лежит в области режимов установившихся колебаний, при выводе на печать или на контроль значений решения, например, только после четного числа шагов (2, 4, 6... или 10, 20, 30...) может создаться иллюзия сходимости решения к искомому решению алгебраического уравнения Риккати. Если же шаг интегрирования лежит в области режимов фиктивного установления, такая иллюзия может возникнуть при любой печати значений решения. Следовательно, контроль сходимости нужно проводить не по разности соседних значений решения, а по матрице производных (правых частей) дифференциального уравнения Риккати.

Выбор метода и шага численного интегрирования для метода установления, очевидно, должен производиться так, чтобы получить решение алгебраического матричного уравнения Риккати (1) с требуемой точностью и за минимальное машинное время. Как показывает анализ процессов получения численных решений уравнений Риккати первого и более высоких порядков, наименьшее число итераций требуется для сходимости при достаточно большом размере шага интегрирования, близком к критическому. При таком размере шага наблюдается более сильное отклонение численного решения дифференциального уравнения Риккати от точного его решения, чем при малых шагах интегрирования, однако это численное решение быстрее выходит на искомое установившееся значение

(фиг. 4, 5). Таким образом, Для получения искомого установившегося значения решения дифференциального уравнения Риккати (2) не требуется большая точность интегрирования. Заметим, что метод Рунге — Кутта, за счет увеличения числа машинных операций {при том же размере шага интегрирования) как раз дает дополнительную, в данном случае излишнюю точность. Сильно же увеличить шаг интегрирования методом Рунге —Кутта по сравнению с методом Эйлера оказывается невозможным, поскольку критические значения шагов интегрирования для этих методов близки. Примеры численного интегрирования уравнений Риккати также показывают, что при выборе для обоих методов оптимальных с точки зрения экономии машинного времени шагов интегрирования число необходимых машинных операций для метода Эйлера намного меньше, чем для метода Рунге —Кутта. Так, в первом из приведенных примеров для получения четырех верных знаков решения при использовании метода Рунге — Кутта потребовалось 120 раз вычислять правую часть дифференциального уравнения Риккати, а при использовании метода Эйлера — всего 40 раз. Применение в данном случае других одношаговых методов интегрирования, более сложных, чем метод Эйлера, по-видимому, также нецелесообразно по тем же причинам, что и указанные выше для метода Рунге — Кутта.

Таким образом, наиболее рациональным методом численного интегрирования дифференциального уравнения Риккати для получения решения алгебраического уравнения Риккати является метод Эйлера, причем шаг интегрирования должен выбираться достаточно большим (а не Л<СЛК). Как показывают примеры решения

различных уравнений Риккати, при удачном выборе шага А этот метод является достаточно быстрым и может конкурировать с другими методами [1, 2] по быстродействию. По простоте же и надежности данный алгоритм превосходит другие алгоритмы.

Автор выражает благодарность И. В. Колину, сделавшему ряд критических замечаний.

ЛИТЕРАТУРА

1. Potter J. Е. Matrix quadratic solutions. .SIAM J. Applied Math.“, vol. 14, N 3, May 1966.

2. Kleinman D. L. On an iterative technique for Riccati equation computations. „IEEE Transactions on Automatic Control", vol. AC—13, N 4, 1968.

3. Kalman R. E., Englar T. A user’s manual for the automatic synthesis program. NASA CR-475, 1966.

4. Репин Ю. М., Третьяков В. Е. Решение задачи об аналитическом конструировании регуляторов на электронных моделирующих устройствах. „Автоматика и телемеханика*, т, 24, № 6, 1963.

5. Брайсон А., Хо-Ю-Ши. Прикладная теория оптимального управления. М., „Мир”, 1972.

6. К а л м а н Р. Об общей теории систем управления. Труды I Конгресса ИФАК, т. 2, Изд-во АН СССР, 1961.

Рукопись поступила 24jVIII 1977 г.