Алгоритм RADI для решения крупномасштабных уравнений Рикатти

Если m > n, R имеет только n строк. Если m <= n, это то же самое, что R = qr(A).[A, S, B ] = svd( inner_piece ) – сингулярное разложение матрицыif any(~isreal(d)) Z= Z*diag(d); elseds = sign(d);d = sqrt(abs(d));Z = Z*diag(d); Y = (Y*diag(d.*ds));end%[QZ, RZ] = qr( Z, 0); [QY, RY] = qr( Y, 0);% % XX * YY^T = (QX*RX) * (QY*RY)^T = QX*(RX*RY^T)*QY^Tinner_piece = RZ*RY'; [A, S, B ] = svd( inner_piece ); if rel ==1 tol = tol*S(1,1); endidx = find(diag(S)>tol, 1, 'last');A = A(:,1:idx); B = B(:,1:idx); S = S(1:idx,1:idx); A = QZ * A; B = QY * B; endЧисленный алгоритмы нахождения сингулярного разложения реализован в модуле lowrank_svdСингулярное разложение матрицы M позволяет вычислять сингулярные числа данной матрицы, а также левые и правые сингулярные векторы матрицы M:левые сингулярные векторы матрицы M — это собственные векторы матрицы MM*;правые сингулярные векторы матрицы M — это собственные векторы матрицы M*M.Сингулярное разложение является удобным при вычислении ранга матрицы, ядра матрицы и псевдообратной матрицы.Сингулярное разложение также используется для приближения матриц матрицами заданного ранга.Неотрицательное вещественное число σ называется сингулярным числом матрицы M тогда и только тогда, когда существуют два вектора единичной длины u∊Km и v∊Kn такие, что:Mv = σu, M*u = σv.Такие векторы u и v называются, соответственно, левым сингулярным вектором и правым сингулярным вектором, соответствующим сингулярному числу σ.Для матрицы M порядка n х m при необходимости приближения матрицей ранга r меньшего, чем n часто используют компактное представление разложения:M = UrΣrVr*.Вычисляются только r столбцов U и r строк V*. Остальные столбцы U и строки V* не вычисляются. Это экономит большое количество памяти при r« n. При работе с матрицей меньшей размерности нужно вычислить сингулярное разложение: M = UΣV* при этом матрица Σ является диагональной. После чего, если нужно сохранить только 90%информации, то нужно взять r, таким образом, чтобы сумма квадратов первых элементов Σ была 90% от общей суммы всех квадратов диагональных элементов Σ.Таким образом получаем U размерностью [mхr] (взяв r столбцов), Σ с размерностью [rх r] и V* с [rх n]. После, вместо матрицы M можно манипулировать матрицей M'=UΣ с меньшей размерностью [mх r], которую часто интерпретируют, как матрицу оценок.ВМатЛабфункция [U,S,V] = svd(A) выполняет сингулярное разложение матрицы A, так что A = U*S*V'.function [U, S, V] = lowrank_svd(C,R, varargin)[Q1, R1] = qr(C,0); [Q2, R2] = qr(R,0); % Q1*R1*R2'*Q2' = C*R': [U, S, V] = svd(R1*R2'); if ~isempty(varargin)if numel(varargin)==3 p = varargin{1}; S = S(1:p, 1:p); U = U(:,1:p); V = V(:,1:p); elseif nargin(varargin)==4 p = varargin{2}; idx = find(diag(S)>p*S(1,1), 1, 'last'); S = S(1:idx, 1:idx); U = U(:,1:idx); V = V(:,1:idx); endendU = Q1*U;V = Q2*V; end3.2. Уравнение Рикатти возможно решить напрямую через встроенные функции МатЛабШаг 1. Инициализация переменных и массивов.На первом шаге инициализируем параметры формируемых массивов: параметр n определяет размер матрицы А nхn, nи m – размер матрицы В nхm, Cnх1и Еmх1 –матрицы размера nи m соответственно.Необходимо создать матрицу размером (n×n) из случайных чисел, равномерно распределенных в диапазоне от 0 до 1: A заполняется случайными числами с помощью функции rand(n, n), аналогично создаем матрицу размером (n×m) из случайных чисел В - rand(n, m) ; C и Е – единичные матрицы. % A'X + XA - XB(R-1)B'X+C'C=0clear alln = 100; m = 5; p =1;A = rand(n, n); B = rand(n, m); C = eye(n); E = eye(m); R = 3;Шаг 2. Подготовка матриц. Матрица Qформируется ортогонально-треугольным разложением с помощью функции qr(A,0). Затем формируются A и В как нормальные и имеющие желаемые собственные значения:%A = rand(n, n); [Q, ~] = qr(A,0); t = -pi + rand(n,1)*2*pi; rad = rand(n,1);eigsA = rad*(r1).*exp(1i*t)+o1;A = Q*diag(eigsA)*Q'; %A является нормальным и имеет желаемые собственные значенияШаг 3. Решение алгебраического уравнения Риккати с непрерывным временем. Функция care (A, B, Q)находит единственное симметричное стабилизирующее решение в виде матрицы X обобщенного алгебраического уравнения Риккати с непрерывным временем.Решение алгебраического уравнения РиккатиX замкнутой системы A, собственные значениякоторой расположены в левой полуплоскости комплексной переменной:[x,l,g] = care(A,B,C'*C);Это дает решение и вычисляет единственное решение X алгебраического уравнения Риккати с непрерывным временем. [eig(A) eig(A-B*g)];Переменная l содержит собственные значения с обратной связьюl=eig(A-B*G,E)Функция также возвращает матрицу усиления,G =R−1(BTXE+ST)Процедура также вычисляет относительную невязку решения, т. е. погрешность подстановки решения в исходное уравнение (погрешность вычислений).Шаг 4. Можно провести диагностику [x,l,g,report] = care(A,B,C'*C);которая возвращает -1, когда связанный гамильтонов пучок имеет собственные значения на мнимой оси или очень близко к ней (отказ); -2, когда нет конечного стабилизирующего решения X; Фробениусову норму относительной невязки, если X существует и конечно.Этот синтаксис не выдает никаких сообщений об ошибках, когда X не существует.В приведенном примере report = 1.6329e-14Шаг 5. Можно убедиться, что полученное решение действительно стабилизируется:plot(1:100, x(1:4,:))Анализируя этот график, можно видеть, что решение стабилизируется уже на 4-ой итерации. Можно убедиться, что это решение действительно стабилизируется, сравнив собственные значения А и А - В*G.Список использованной литературыBenner, Peter; Köhler, Martin; Saak, Jens: Matrix equations, sparse solvers: \textttM-M.E.S.S.-2.0.1 -- philosophy, features, and application for (parametric) model order reduction (2021)Benner, Peter; Palitta, Davide: On the solution of the nonsymmetric T-Riccati equation (2021)Chen, Hao; Wang, Ying: Modified Douglas splitting method for differential matrix equations (2021)Dorschky, Ines; Reis, Timo; Voigt, Matthias: Balanced truncation model reduction for symmetric second order systems -- a passivity-based approach (2021)Benner, Peter; Bujanović, Zvonimir; Kürschner, Patrick; Saak, Jens: A numerical comparison of different solvers for large-scale, continuous-time algebraic Riccati equations and LQR problems (2020)Зеликин М. И. Однородные пространства и уравнение Риккати в вариационном исчислении, — Факториал, Москва, 1998.Егоров А. И. Уравнения Риккати, — Физматлит, Москва, 2001.Лауфер М. Я. О решении уравнений Риккати // Лауфер М. Я. Избранные задачи математической физики. Сб. статей.— Северодвинск: НТО кораблестроителей им. акад. А. Н. Крылова, Севмашвтуз, Северодв. отд-ние Ломоносов. фонда, 2005.— стр. 137—140.— ISBN 5-7723-0605-9.