ИЗВЕСТИЯ РАН. ТЕОРИЯ И СИСТЕМЫ УПРАВЛЕНИЯ, 2015, № 2, с. 3-21
РЕДУКЦИЯ БОЛЬШИХ ДИНАМИЧЕСКИХ СИСТЕМ МЕТОДОМ ПОДПРОСТРАНСТВ А.Н. КРЫЛОВА. АНАЛИЗ ПОДХОДОВ* © 2015 г. Н. Е. Зубов, Е. А. Микрин, М. Ш. Мисриханов,
г. Королев МО, ОАО РКК "Энергия", Москва, МГТУим. Н.Э. Баумана, МФТИ Поступила в редакцию 24.07.14 г., после доработки 27.10.14 г.
Описываются подходы и численные алгоритмы понижения порядка (редукции) математических моделей многомерных динамических систем в пространстве состояний на основе метода подпространств А.Н. Крылова. Для вычисления матриц в представлении редуцированной модели в пространстве состояний применены: метод Ланцоша и метод Арнольди. Приведены решения практических примеров редукции больших систем.
DOI: 10.7868/S0002338815020146
Введение. В 1931 г. акад. А.Н. Крылов опубликовал работу [1] на тему, известную теперь как подпространства Крылова или метод подпространств Крылова. Статья была посвящена вычислению коэффициентов характеристического полинома заданной числовой матрицы. В этом труде акад. Крылов коснулся эффективности расчетов и подсчитал вычислительные затраты как количество "отдельных операций перемножения" (явление не типичное для математической публикации 1931 г.). Он представил метод, который был лучшим из известных к тому времени вычислительных методов и с тех пор широко используется во всем мире, особенно при итеративном решении систем линейных матричных уравнений больших порядков.
Метод подпространств Крылова нашел также широкое применение при анализе устойчивости, управляемости и наблюдаемости линейных динамических систем в пространстве состояний, а также при синтезе законов управления с обратной связью (см., например, [2—14] и др.).
Другим направлением использования метода подпространств Крылова является редукция математических моделей динамических систем, заданных, например, в виде марковских процессов. Редукция математической модели динамической системы, основанная на методе подпространств Крылова, служит альтернативой известному методу редукции, основанному на 8УО-разложении. Она применяется для больших динамических систем с размерностью пространства состояний п > 1, когда матрицы, входящие в описание системы, бывают достаточно разреженными.
Ключевым компонентом редукции на основе подпространств Крылова является так называемая подгонка моментов (параметров Маркова). Ее идея заключается в установлении соответствия между моментами передаточной функции (матрицы) исходной модели системы £ высокого порядка при разложении ее в ряд Лорана и моментами модели системы низкого порядка 2 г. При этом используются управляемое и наблюдаемое подпространства системы.
Данная работа посвящена анализу известных подходов к редукции математических моделей динамических систем и численных алгоритмов, построенных на основе метода подпространств Крылова и получению с использованием приведенных алгоритмов решений для трудно разрешимых задач.
1. Описание динамических систем. Рассмотрим динамическую систему 2 (линейную стационарную систему или сокращенно ЬТТ-систему) следующего вида:
А. В. Пролетарский, В. Н. Рябченко
(1.1)
*Исследование выполнено за счет гранта Российского научного фонда (проект № 14-11-00046).
где x = (( ¡•••jxn)T е Rn — вектор состояния; u = (щ ¡•••jum)T е Rm — вектор входа; y = (y1 yp)T е е № — вектор выхода; x0 — начальное условие.
Если m = p = 1, тогда (1.1) именуется SISO-системой (Single Input Single Output). Если m,p > 1, тогда (1.1) называется MIMO-системой (Multi Input Multi Output). При m = 1, p > 1 (1.1) называется SIMO-системой (Single Input Multi Output), а в случае m > 1, p = 1 — соответственно MISO-системой (Multi Input Single Output).
Определим для системы 2 (1.1) при нулевых начальных условиях оператор свертки У от u(t) к y(t) [12]:
да
У : u(t) ^ y(t) = G * u = jG(t - x)u(x)dт (1.2)
о
с ядром в виде весовой (импульсной) функции (матрицы) G(t). Преобразование Лапласа этой функции G(t) дает передаточную функцию (матрицу)
F(s) = Ж (G)(5) = jG(t) e~s'dt = C (sIn - A)-1 B + D. (1.3)
о
Здесь In — единичная матрица размера n x n.
В дальнейшем будем считать, что система £ (1.1) является полностью управляемой [15], т.е. матрица управляемости
% = (b|ab|a 2B |-| A n-1B) (1.4)
имеет полный ранг
rank% = n, (1.5)
и полностью наблюдаемой [15], т.е. матрица наблюдаемости
О = (CT |ATCt I (AT)2CT |-| (AT)n-1CT) (1.6)
также имеет полный ранг
rankO = n. (1.7)
Пусть система £ (1.1) — устойчивая по Ляпунову, тогда существуют матрицы-грамианы 9 и & , определенные как [5]
да да
9 = \e ATBBTe ATVx, & = \e A^CTCe AVx (1.8)
о
и имеющие свойства 9 = 9т > 0, & = &т > 0. Более того, матрицы (1.8) являются решениями соответствующих алгебраических уравнений Ляпунова
А9 + 9 Ат + ВВт = 0, А т& + &Л + С тС = 0. (1.9)
Вернемся к оператору свертки У (1.2) и, ограничивая его области определения и значений, введем понятие оператора Ганкеля [12]:
0
Ж : и_(0 ^ у+(0 = Ж (и_) = | в(? -т)и_(т) йт, г > 0. (1.10)
—да
Здесь, нестрого говоря, и _($) — "прошлые" входы, у+(г) — "будущие" выходы. Другими словами, оператор Ганкеля отображает "прошлые" входы в "будущие" выходы. В отличие от оператора свертки оператор Ганкеля имеет конечный ранг [12], не превосходящий п, и, следовательно, конечное множество сингулярных значений ст,(£), определенных следующим образом.
Пусть система £ (1.1) — полностью управляемая, полностью наблюдаемая и асимптотически устойчивая. Тогда сингулярные значения CTi(£) есть положительные квадратные корни собственных значений Xi произведения грамианов ФИ из (1.8) или (1.9) [16], т.е.
ст;(2) = i (ФИ), i = й (1.11)
Максимальное сингулярное число оператора Ганкеля a max(2) определяет ганкелеву норму системы £ (1.1) [12], т.е.
И Ж = ° max(2). (1.12)
При этом норма Ж „ системы £ (1.1) равна
||£|| Ж m = sup ^ max
(С(/ш)}, (1.13)
™ шеЖ
а норма Ж 2 —
/ » Л°.5
|1ж 2 _
J tr(G*(y ®)G(/ rn))d ю
(1.14)
Здесь символ * означает сопряжение.
2. Подпространства Крылова. Метод подпространств Крылова нашел широкое применение при решении разнообразных вычислительных задач и прежде всего в процессе итеративного решения больших матричных уравнений вида [17]
Ax = b, (2.1)
где A, b — заданные матрица n х n и вектор n х 1 соответственно.
Пусть x0 е Жп — начальная аппроксимация (приближение) решения уравнения (2.1), где
Ж = К или Ж = С,
Го = b - Ax0 (2.2)
— начальная невязка и пусть
Ж m (A, Го) = span {го, Aro, A 2Го, ..., A и-1Го| (2.3)
— подпространство размерности m, определенное с помощью матрицы A и вектора го. Для подпространств (2.3) справедливо включение
Ж m ^ Ж m+1.
Методы, основанные на подпространствах Ж, являются итеративными [18]. С их помощью аппроксимация решения уравнения (2.1) на m -м шаге находится в виде
Xm = Х0 + Qm-1 (A) Го, (2.4)
где qm-1 — матричный полином от A степени не больше чем m - 1. Если матричное уравнение (2.1) задано над К, тогда qm-1 также имеет действительные коэффициенты.
Согласно (2.4), невязка решения уравнения (2.1) на m -м шаге может быть записана в виде
rm = b - Axm = Го - Aqm-1 (A) Го = Pm (A) Го, (2.5)
где pm — остаточный полином. Аналогично если x* — решение (2.1), то
Хm - Х* = pm (A) (X0 - X*) .
В процессе итеративного построения базиса Ж m каждый метод, основанный на подпространствах (2.3), на каждой итерации использует одно или два умножения матрицы на вектор в форме
z = Av (или y = A w). Поэтому алгоритмы данных методов могут применяться как для решения задач, в которых матрица системы задана в явной форме, так и для решения задач, в которых мат-
рица доступна только через операцию умножения на вектор [16]. Следует заметить, что сходимость итерационного процесса обеспечивается применением того или иного алгоритма.
Во всех методах на основе подпространств Ж m инициализация итеративных процедур осуществляется путем задания начального приближения x0 и соответствующей начальной невязки (2.2). Аппроксимация решения на m -м шаге находится в виде
x о + Ж m (A, Го ) .
Без ограничения общности можно полагать, что хо = о и го = b. Исходя из сказанного, введем определение.
Определение 1. Для системы £ (1.1), заданной парой (A, b), где A е [Rnxn и b е Rn, к-я последовательность Крылова образуется следующими к вектор-столбцами:
Жk (A,b) = {b, Ab, A2b,..., Ak-1b}. (2.6)
Линейное пространство, образованное оболочкой, натянутой на столбцы векторов (2.6), т.е.
Ж k (A,b) = span {b, Ab, A2b,..., Ak-1b}, (2.7)
называется к-м подпространством Крылова.
Согласно теореме Гамильтона—Кэли, для подпространств Крылова справедливо тождество
V k > n : Жk (A,b) = Ж n (A,b). (2.8)
Свойство (2.8) определяет непосредственную связь подпространств Крылова с матрицами управляемости % (1.4) и наблюдаемости О (1.6) системы £ (1.1).
Обобщениями данного подхода являются так называемые частичная матрица достижимости
Ж k (A, B, а) = ((ctIn - A)-1 BI (ctIn - A)-2 B | • • • | (aIn - A)-k B) (2.9)
и частичная матрица наблюдаемости
Оk (C, A,a) = ((aIn - AT)-1CT (aIn - AT)-kCT). (2.10)
3. Проекторы. Проектор P e С "x" в подпространство Ж с Сn определяется как линейное идем-потентное отображение в Сn [19], т.е.
P2 = P, Px е Ж, Vx е Cn. (3.1) В силу идемпотентности проектора P выполняется следующее равенство:
ker (P) = range (I - P), (3.2) т.е. ядро проектора P совпадает с областью значений отображения (I - P). Каждый вектор x е Сn может быть представлен в виде разложения
x = Px + (I - P) x, (3.3)
и, следовательно, пространство Сn может быть декомпозировано в прямую сумму подпространств
Сп = ker (P) © range (P). (3.4)
Из равенства (3.4) вытекает, что для каждой пары подпространств Ж с С", Ж с С", которые в
виде прямой суммы типа (3.4) определяют Сn, и для каждого вектора x е Сn вектор Px соответствует условиям
Px е Ж, (3.5)
x - Px е Ж. (3.6)
При этом говорят, что линейное отображение P е Сnxn проектирует x е Сn в Ж с Сn в том же направлении или параллельно подпространству Я с Сn.
Ортогональный проектор P в подпространство Я с Сп также определяется как линейное отображение Сn с условием (3.1) и дополнительным требов
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.