Contents

function LQG_FDM_unstable_radi(n0, n_unstable, istest)
% Computes stabilizing and detecting solutions of the regulator and filter
% Riccati equations via the RADI, respectively. The solution of the
% regulator equation is approximated via a ZZ' factorization, and for the
% filter equation an LDL' approximation is used.
% Note: The LDL' case is only used for demonstration, but with no
% underlying necessity here.
%
% Inputs:
%
% n0          n0^2 gives the dimension of the original model, i.e. n0 is
%             the number of degrees of freedom, i.e. grid points, per
%             spatial direction
%             (optional; defaults to 30)
%
% n_unstable  number of unstable eigenvalues by construction
%             (optional; defaults to 5)
%
% istest      flag to determine whether this demo runs as a CI test or
%             interactive demo
%             (optional, defaults to 0, i.e. interactive demo)
%

%
% This file is part of the M-M.E.S.S. project
% (http://www.mpi-magdeburg.mpg.de/projects/mess).
% Copyright (c) 2009-2025 Jens Saak, Martin Koehler, Peter Benner and others.
% All rights reserved.
% License: BSD 2-Clause License (see COPYING)
%

narginchk(0, 3);
if nargin < 1
    n0 = 30;
end
if nargin < 2
    n_unstable = 5;
end
if nargin < 3
    istest = false;
end

construct unstable matrix

Generate the stable part.

A0         = fdm_2d_matrix(n0, '10*x', '1000*y', '0');
n_stable   = size(A0, 1);
n          = n_stable + n_unstable;
B          = [eye(n_stable, n_unstable); diag(1:n_unstable)];

C                                 = [transpose(B); zeros(2, n)];
C(n_unstable + 1, n_unstable + 1) = 1;
C(n_unstable + 2, n_unstable + 2) = 1;

% The instability is added by fixing the solution Y of the partial
% stabilization problem via Bernoulli and solving the resulting Lyapunov
% equation for the unstable matrix Au:
%
%   Au'*Y + Y*Au - Y*Bt*Bt'*Y = 0.
%
% The resulting stabilizing feedback is defined such that the unstable
% eigenvalues of the system are mirrored into the left open half-plane.

Y   = eye(n_unstable);
Bt  = B(n_stable + 1:n, :); % Part of B acting on unstable part.
Vl  = [zeros(n_stable, n_unstable); eye(n_unstable)];
if exist('lyap', 'file') % Solve Lyapunov equation.
    Au = lyap(Y, -Bt * Bt');
else
    Au = lyap2solve(Y, (-Bt * Bt')');
end

ZB0 = Vl;
YB0 = Y \ eye(n_unstable);
WB0 = C';

ZC0 = ZB0;
YC0 = YB0;
WC0 = B;

% The full A is constructed via additive decomposition (block-diagonal).
eqn.A_     = blkdiag(A0, Au);
eqn.B      = B;
eqn.C      = C;
eqn.haveE  = false;

% set operation
opts = struct();
[oper, opts] = operatormanager(opts, 'default');

global options

opts.norm = 'fro';

RADI parameters

opts.shifts.naive_update_mode = false;
opts.radi.compute_sol_fac     = true;
opts.radi.compute_res         = false;
opts.radi.get_ZZt             = true;
opts.radi.maxiter             = 200;
opts.radi.res_tol             = 1.0e-10;
opts.radi.rel_diff_tol        = 1.0e-16;
opts.radi.info                = 1;

% Only for mess_res2_norms (not needed in RADI):
opts.radi.res = struct('maxiter', 10, ...
                       'tol', 1.0e-06, ...
                       'orth', 0);

shift parameters via projection

opts.shifts.num_desired = 5;
opts.shifts.method      = 'projection';

for shift banned eigenvalues

opts.shifts.banned     = -eig(Au);
opts.shifts.banned_tol = 1e-6;

Solve regulator Riccati equation.

fprintf('Solve regulator Riccati equation\n');
fprintf('--------------------------------\n');
eqn.type     = 'T';
opts.LDL_T = false;
opts.radi.Z0 = ZB0;
opts.radi.Y0 = YB0;
opts.radi.W0 = WB0;

time_radi_regulator = tic;
[outReg, eqn, opts, oper] = mess_lrradi(eqn, opts, oper);
t_elapsed1 = toc(time_radi_regulator);

% Size of solution factor.
fprintf('\nSize outReg.Z: %d x %d\n', ...
        size(outReg.Z, 1), size(outReg.Z, 2));

% Check residuals.
res0 = norm(eqn.C * eqn.C');
res1 = mess_res2_norms(outReg.Z, 'riccati', ...
                       eqn, opts, oper, opts.radi, []) / res0;
res2 = abs(eigs(@(x) eqn.A_' * (outReg.Z * (outReg.Z' * x)) + ...
                (outReg.Z * (outReg.Z' * (eqn.A_ * x)))  + ...
                eqn.C' * (eqn.C * x) - outReg.K' * (outReg.K * x), ...
                n, 1, 'LM')) / res0;

fprintf('solving the regulator Riccati equation took %6.2f seconds \n', ...
        t_elapsed1);
fprintf(['Residual computation -- RADI: %e | ' ...
         'mess_res2_norms: %e | eigs: %e \n'], ...
        outReg.res(end), res1, res2);

% Print convergence behavior.
if istest
    if min(outReg.res) >= opts.radi.res_tol
        mess_err(opts, 'TEST:accuracy', 'unexpectedly inaccurate result');
    end
else
    figure(1);
    semilogy(outReg.res, 'LineWidth', 3);
    title('0 = C^T C + A^T X + X A - X B B^T X');
    xlabel('number of iteration steps');
    ylabel('normalized residual norm');
    pause(1);
end

fprintf('\n');
Solve regulator Riccati equation
--------------------------------
RADI step:    1 pc: -1.250003e+01 + 0.000000e+00i normalized residual: 2.504853e-01 relative change in Z: 6.334954e-01  
RADI step:    2 pc: -8.000025e+00 + 0.000000e+00i normalized residual: 1.248903e-01 relative change in Z: 2.307932e-01  
RADI step:    3 pc: -4.500025e+00 + 0.000000e+00i normalized residual: 9.837321e-02 relative change in Z: 1.554579e-01  
RADI step:    4 pc: -2.000025e+00 + 0.000000e+00i normalized residual: 9.051185e-02 relative change in Z: 1.122526e-01  
RADI step:    5 pc: -5.000250e-01 + 0.000000e+00i normalized residual: 8.921764e-02 relative change in Z: 3.235217e-02  
RADI step:    6 pc: -2.759799e+01 + 0.000000e+00i normalized residual: 7.456561e-02 relative change in Z: 4.044819e-02  
RADI step:    7 pc: -6.520281e+03 + 0.000000e+00i normalized residual: 1.095846e-02 relative change in Z: 9.569335e-03  
RADI step:    8 pc: -1.116444e+00 + 0.000000e+00i normalized residual: 1.090575e-02 relative change in Z: 1.157454e-02  
RADI step:    9 pc: -1.077506e+03 + 0.000000e+00i normalized residual: 1.053820e-03 relative change in Z: 4.404639e-03  
RADI step:   10 pc: -4.471525e+00 + 0.000000e+00i normalized residual: 1.040873e-03 relative change in Z: 1.394344e-03  
RADI step:   11 pc: -2.794687e+01 + 0.000000e+00i normalized residual: 9.980822e-04 relative change in Z: 1.034868e-03  
RADI step:   12 pc: -6.526749e+03 + 0.000000e+00i normalized residual: 1.373945e-04 relative change in Z: 1.081194e-03  
RADI step:   13 pc: -1.117897e+00 + 0.000000e+00i normalized residual: 1.371014e-04 relative change in Z: 1.115304e-04  
RADI step:   14 pc: -1.041276e+03 + 0.000000e+00i normalized residual: 2.386307e-05 relative change in Z: 4.146519e-04  
RADI step:   15 pc: -4.471514e+00 + 0.000000e+00i normalized residual: 2.363546e-05 relative change in Z: 1.572055e-04  
RADI step:   16 pc: -2.794687e+01 + 0.000000e+00i normalized residual: 2.268074e-05 relative change in Z: 1.171042e-04  
RADI step:   18 pc: -6.258103e+03 + -3.947523e+02i normalized residual: 5.405305e-07 relative change in Z: 1.630778e-04  
RADI step:   19 pc: -1.117898e+00 + 0.000000e+00i normalized residual: 5.390402e-07 relative change in Z: 1.259984e-05  
RADI step:   20 pc: -1.091178e+03 + 0.000000e+00i normalized residual: 7.032812e-08 relative change in Z: 2.683536e-05  
RADI step:   21 pc: -1.216261e+03 + 0.000000e+00i normalized residual: 1.149965e-08 relative change in Z: 1.042322e-05  
RADI step:   22 pc: -1.006091e+01 + 0.000000e+00i normalized residual: 9.455159e-09 relative change in Z: 2.053063e-05  
RADI step:   24 pc: -6.290886e+03 + -1.037961e+03i normalized residual: 4.509840e-10 relative change in Z: 3.297730e-06  
RADI step:   25 pc: -1.788603e+01 + 0.000000e+00i normalized residual: 2.371104e-10 relative change in Z: 5.642218e-06  
RADI step:   26 pc: -1.222450e+03 + 0.000000e+00i normalized residual: 2.903984e-11 relative change in Z: 5.232219e-07  

Size outReg.Z: 905 x 32
Warning: resopts.res takes precedence over resopts; please use resopts as
documented 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/FDM/LQG_FDM_unstable_radi.m',135)">LQG_FDM_unstable_radi (line 135)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/norms/mess_res2_norms.m',113)">mess_res2_norms (line 113)</a>
solving the regulator Riccati equation took   0.19 seconds 
Residual computation -- RADI: 2.903984e-11 | mess_res2_norms: 8.155423e-11 | eigs: 8.155418e-11 

Solve filter Riccati equation.

fprintf('Solve filter Riccati equation\n');
fprintf('-----------------------------\n');

eqn.type     = 'N';
opts.LDL_T   = true;
opts.radi.Z0 = ZC0;
opts.radi.Y0 = YC0;

% Some additional scaling for non-trivial LDL' term.
eqn.B        = eqn.B * diag(1 ./ sqrt(1:n_unstable));
eqn.R = diag(1:n_unstable);
eqn.Q = eye(size(eqn.C, 1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set corresponding residual.
opts.radi.W0 = WC0 * diag(1 ./ sqrt(1:n_unstable));
opts.radi.T0 = diag(1:n_unstable);

time_radi_filter = tic;
[outFil, eqn, opts, oper] = mess_lrradi(eqn, opts, oper);
t_elapsed2 = toc(time_radi_filter);

% Size of solution factor.
fprintf('\nSize outFil.Z: %d x %d\n', ...
        size(outFil.Z, 1), size(outFil.Z, 2));

% Check residuals.
res0       = norm((eqn.B * eqn.R) * eqn.B');

res3       = mess_res2_norms(outFil.Z, 'riccati', ...
                             eqn, opts, oper, opts.radi, outFil.D) / res0;
ZD         = outFil.Z * outFil.D;
res4       = abs(eigs(@(x) eqn.A_ * (ZD * (outFil.Z' * x)) + ...
                      (ZD * (outFil.Z' * (eqn.A_' * x))) + ...
                      eqn.B * (eqn.R * (eqn.B' * x)) - (outFil.K)' * ...
                      ((outFil.K) * x), n, 1, 'LM')) / res0;

fprintf('solving the filter Riccati equation took %6.2f seconds \n', ...
        t_elapsed2);
fprintf(['Residual computations -- RADI: %e | ' ...
         'mess_res2_norms: %e | eigs: %e \n'], ...
        outFil.res(end), res3, res4);

% Print convergence behavior.
if istest
    if min(outFil.res) >= opts.radi.res_tol
        mess_err(opts, 'TEST:accuracy', 'unexpectedly inaccurate result');
    end
else
    figure(2);
    semilogy(outFil.res, 'LineWidth', 3);
    title('0 = B S B^T + A Y + Y A^T - Y C^T C Y');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
Solve filter Riccati equation
-----------------------------
RADI step:    1 pc: -1.250003e+01 + 0.000000e+00i normalized residual: 2.471852e-01 relative change in Z: 8.425445e-01  
RADI step:    2 pc: -8.000025e+00 + 0.000000e+00i normalized residual: 1.179628e-01 relative change in Z: 2.442724e-01  
RADI step:    3 pc: -4.500025e+00 + 0.000000e+00i normalized residual: 8.950024e-02 relative change in Z: 1.303531e-01  
RADI step:    4 pc: -2.000025e+00 + 0.000000e+00i normalized residual: 8.083576e-02 relative change in Z: 8.548480e-02  
RADI step:    5 pc: -5.000250e-01 + 0.000000e+00i normalized residual: 7.940022e-02 relative change in Z: 2.739268e-02  
RADI step:    6 pc: -2.759799e+01 + 0.000000e+00i normalized residual: 6.388342e-02 relative change in Z: 6.269812e-02  
RADI step:    7 pc: -5.014896e+03 + 0.000000e+00i normalized residual: 1.165382e-02 relative change in Z: 1.054083e-02  
RADI step:    8 pc: -1.116444e+00 + 0.000000e+00i normalized residual: 1.162231e-02 relative change in Z: 8.254271e-03  
RADI step:   10 pc: -6.219560e+02 + -5.812256e+02i normalized residual: 1.580419e-02 relative change in Z: 9.501003e-03  
RADI step:   11 pc: -2.794687e+01 + 0.000000e+00i normalized residual: 1.451894e-02 relative change in Z: 3.122320e-03  
RADI step:   12 pc: -5.945165e+03 + 0.000000e+00i normalized residual: 1.231246e-02 relative change in Z: 5.577753e-03  
RADI step:   13 pc: -1.117897e+00 + 0.000000e+00i normalized residual: 1.228162e-02 relative change in Z: 3.862549e-04  
RADI step:   15 pc: -6.293837e+02 + -1.465574e+03i normalized residual: 5.688371e-04 relative change in Z: 5.880379e-03  
RADI step:   16 pc: -2.794688e+01 + 0.000000e+00i normalized residual: 5.596881e-04 relative change in Z: 7.449497e-04  
RADI step:   18 pc: -9.118909e+02 + -3.038603e+03i normalized residual: 4.287436e-04 relative change in Z: 1.540944e-03  
RADI step:   19 pc: -1.117897e+00 + 0.000000e+00i normalized residual: 4.266501e-04 relative change in Z: 1.144259e-04  
RADI step:   21 pc: -1.054685e+03 + -1.036088e+02i normalized residual: 6.088963e-05 relative change in Z: 1.274954e-03  
RADI step:   22 pc: -1.107356e+03 + 0.000000e+00i normalized residual: 6.080957e-05 relative change in Z: 2.998477e-04  
RADI step:   23 pc: -4.472160e+00 + 0.000000e+00i normalized residual: 6.061245e-05 relative change in Z: 7.767297e-05  
RADI step:   25 pc: -9.195846e+02 + -3.884575e+03i normalized residual: 2.211912e-05 relative change in Z: 2.110609e-04  
RADI step:   26 pc: -1.788599e+01 + 0.000000e+00i normalized residual: 2.193512e-05 relative change in Z: 5.462805e-05  
RADI step:   27 pc: -9.428568e+02 + 0.000000e+00i normalized residual: 1.097522e-05 relative change in Z: 1.693259e-04  
RADI step:   28 pc: -4.471521e+00 + 0.000000e+00i normalized residual: 1.092157e-05 relative change in Z: 1.433438e-05  
RADI step:   30 pc: -8.604715e+02 + -4.586307e+03i normalized residual: 5.185893e-06 relative change in Z: 1.179754e-04  
RADI step:   31 pc: -1.788603e+01 + 0.000000e+00i normalized residual: 5.057562e-06 relative change in Z: 1.764861e-05  
RADI step:   32 pc: -1.670155e+03 + 0.000000e+00i normalized residual: 1.572012e-06 relative change in Z: 1.212355e-04  
RADI step:   33 pc: -9.512802e+00 + 0.000000e+00i normalized residual: 1.573741e-06 relative change in Z: 7.462918e-06  
RADI step:   35 pc: -1.020421e+03 + -5.345332e+03i normalized residual: 6.442747e-07 relative change in Z: 4.423560e-05  
RADI step:   37 pc: -9.965875e+02 + -3.708754e+03i normalized residual: 4.253481e-07 relative change in Z: 3.217756e-05  
RADI step:   38 pc: -2.058356e+03 + 0.000000e+00i normalized residual: 2.126599e-07 relative change in Z: 3.623041e-05  
RADI step:   39 pc: -1.006090e+01 + 0.000000e+00i normalized residual: 2.117168e-07 relative change in Z: 1.349507e-06  
RADI step:   41 pc: -1.319279e+03 + -6.671383e+03i normalized residual: 1.285788e-07 relative change in Z: 8.884593e-06  
RADI step:   43 pc: -1.224416e+03 + -3.118018e+03i normalized residual: 6.424611e-09 relative change in Z: 1.133633e-05  
RADI step:   44 pc: -2.626692e+03 + 0.000000e+00i normalized residual: 4.907051e-09 relative change in Z: 3.409770e-06  
RADI step:   45 pc: -1.453144e+01 + 0.000000e+00i normalized residual: 4.810975e-09 relative change in Z: 4.329954e-07  
RADI step:   47 pc: -1.154669e+03 + -7.261482e+03i normalized residual: 3.292811e-09 relative change in Z: 1.656437e-06  
RADI step:   49 pc: -1.068737e+03 + -2.897572e+03i normalized residual: 4.664451e-10 relative change in Z: 2.159986e-06  
RADI step:   51 pc: -3.388976e+03 + -4.327969e+02i normalized residual: 1.207816e-10 relative change in Z: 8.067781e-07  
RADI step:   53 pc: -1.282176e+03 + -8.496703e+03i normalized residual: 1.921627e-11 relative change in Z: 3.357603e-07  

Size outFil.Z: 905 x 73
Warning: resopts.res takes precedence over resopts; please use resopts as
documented 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/FDM/LQG_FDM_unstable_radi.m',193)">LQG_FDM_unstable_radi (line 193)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/norms/mess_res2_norms.m',113)">mess_res2_norms (line 113)</a>
solving the filter Riccati equation took   0.32 seconds 
Residual computations -- RADI: 1.921627e-11 | mess_res2_norms: 1.107844e-10 | eigs: 1.107815e-10