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 © 2009-2022 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 = 0; 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  = 0;

% set operation
oper = operatormanager('default');

global options

opts.norm = 'fro';

RADI parameters

opts.shifts.naive_update_mode = 0;
opts.radi.compute_sol_fac     = 1;
opts.radi.compute_res         = 0;
opts.radi.get_ZZt             = 1;
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   = 0;
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
        error('MESS: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.494748e+03 + 0.000000e+00i normalized residual: 1.089110e-02 relative change in Z: 9.575758e-03
RADI step:    8 pc: -1.116444e+00 + 0.000000e+00i normalized residual: 1.083860e-02 relative change in Z: 1.157451e-02
RADI step:    9 pc: -1.077506e+03 + 0.000000e+00i normalized residual: 1.040756e-03 relative change in Z: 4.392728e-03
RADI step:   10 pc: -4.471525e+00 + 0.000000e+00i normalized residual: 1.027884e-03 relative change in Z: 1.394280e-03
RADI step:   11 pc: -2.794687e+01 + 0.000000e+00i normalized residual: 9.855012e-04 relative change in Z: 1.034566e-03
RADI step:   12 pc: -6.526616e+03 + 0.000000e+00i normalized residual: 1.362372e-04 relative change in Z: 1.074209e-03
RADI step:   13 pc: -1.117897e+00 + 0.000000e+00i normalized residual: 1.359461e-04 relative change in Z: 1.115155e-04
RADI step:   14 pc: -1.041276e+03 + 0.000000e+00i normalized residual: 2.361170e-05 relative change in Z: 4.130731e-04
RADI step:   15 pc: -4.471514e+00 + 0.000000e+00i normalized residual: 2.338554e-05 relative change in Z: 1.571960e-04
RADI step:   16 pc: -2.794687e+01 + 0.000000e+00i normalized residual: 2.243961e-05 relative change in Z: 1.170509e-04
RADI step:   18 pc: -6.227725e+03 + -3.597591e+02i normalized residual: 5.256882e-07 relative change in Z: 1.623687e-04
RADI step:   19 pc: -1.117898e+00 + 0.000000e+00i normalized residual: 5.242264e-07 relative change in Z: 1.259725e-05
RADI step:   20 pc: -1.091629e+03 + 0.000000e+00i normalized residual: 6.797841e-08 relative change in Z: 2.649700e-05
RADI step:   21 pc: -1.216101e+03 + 0.000000e+00i normalized residual: 1.115192e-08 relative change in Z: 1.029177e-05
RADI step:   22 pc: -1.006091e+01 + 0.000000e+00i normalized residual: 9.091590e-09 relative change in Z: 2.052900e-05
RADI step:   24 pc: -6.274322e+03 + -1.064010e+03i normalized residual: 4.449329e-10 relative change in Z: 3.234767e-06
RADI step:   25 pc: -1.788603e+01 + 0.000000e+00i normalized residual: 2.292637e-10 relative change in Z: 5.641699e-06
RADI step:   26 pc: -1.222273e+03 + 0.000000e+00i normalized residual: 2.796066e-11 relative change in Z: 5.144872e-07

Size outReg.Z: 905 x 32
solving the regulator Riccati equation took   0.18 seconds 
Residual computation -- RADI: 2.796066e-11 | mess_res2_norms: 8.155300e-11 | eigs: 8.155296e-11 

Solve filter Riccati equation.

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

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

% Some additional scaling for non-trivial LDL' term.
eqn.S        = diag(1:n_unstable);
eqn.B        = eqn.B * diag(1 ./ sqrt(1:n_unstable));

% Set corresponding residual.
opts.radi.W0 = WC0 * diag(1 ./ sqrt(1:n_unstable));
opts.radi.S0 = 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.S) * eqn.B');
eqn.S_diag = diag(eqn.S);
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.S * (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
        error('MESS: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.014781e+03 + 0.000000e+00i normalized residual: 1.165367e-02 relative change in Z: 1.054087e-02
RADI step:    8 pc: -1.116444e+00 + 0.000000e+00i normalized residual: 1.162216e-02 relative change in Z: 8.254271e-03
RADI step:   10 pc: -6.220259e+02 + -5.816253e+02i normalized residual: 1.580948e-02 relative change in Z: 9.500773e-03
RADI step:   11 pc: -2.794687e+01 + 0.000000e+00i normalized residual: 1.452220e-02 relative change in Z: 3.124170e-03
RADI step:   12 pc: -5.945117e+03 + 0.000000e+00i normalized residual: 1.231813e-02 relative change in Z: 5.578680e-03
RADI step:   13 pc: -1.117897e+00 + 0.000000e+00i normalized residual: 1.228722e-02 relative change in Z: 3.865830e-04
RADI step:   15 pc: -6.293714e+02 + -1.465612e+03i normalized residual: 5.689496e-04 relative change in Z: 5.881489e-03
RADI step:   16 pc: -2.794688e+01 + 0.000000e+00i normalized residual: 5.596639e-04 relative change in Z: 7.459946e-04
RADI step:   18 pc: -9.118950e+02 + -3.038626e+03i normalized residual: 4.300191e-04 relative change in Z: 1.540990e-03
RADI step:   19 pc: -1.117897e+00 + 0.000000e+00i normalized residual: 4.279174e-04 relative change in Z: 1.146041e-04
RADI step:   21 pc: -1.054695e+03 + -1.036002e+02i normalized residual: 6.087122e-05 relative change in Z: 1.277466e-03
RADI step:   22 pc: -1.107375e+03 + 0.000000e+00i normalized residual: 6.078395e-05 relative change in Z: 2.997363e-04
RADI step:   23 pc: -4.472160e+00 + 0.000000e+00i normalized residual: 6.058693e-05 relative change in Z: 7.767462e-05
RADI step:   25 pc: -9.195671e+02 + -3.884581e+03i normalized residual: 2.210906e-05 relative change in Z: 2.110296e-04
RADI step:   26 pc: -1.788599e+01 + 0.000000e+00i normalized residual: 2.192515e-05 relative change in Z: 5.462723e-05
RADI step:   27 pc: -9.430082e+02 + 0.000000e+00i normalized residual: 1.096999e-05 relative change in Z: 1.693018e-04
RADI step:   28 pc: -4.471521e+00 + 0.000000e+00i normalized residual: 1.091638e-05 relative change in Z: 1.433180e-05
RADI step:   30 pc: -8.604730e+02 + -4.586325e+03i normalized residual: 5.182548e-06 relative change in Z: 1.179524e-04
RADI step:   31 pc: -1.788603e+01 + 0.000000e+00i normalized residual: 5.054292e-06 relative change in Z: 1.764528e-05
RADI step:   32 pc: -1.669957e+03 + 0.000000e+00i normalized residual: 1.571658e-06 relative change in Z: 1.212049e-04
RADI step:   33 pc: -9.514423e+00 + 0.000000e+00i normalized residual: 1.573388e-06 relative change in Z: 7.461725e-06
RADI step:   35 pc: -1.020407e+03 + -5.345356e+03i normalized residual: 6.440315e-07 relative change in Z: 4.423136e-05
RADI step:   37 pc: -9.965652e+02 + -3.708808e+03i normalized residual: 4.251783e-07 relative change in Z: 3.217101e-05
RADI step:   38 pc: -2.058298e+03 + 0.000000e+00i normalized residual: 2.126008e-07 relative change in Z: 3.622314e-05
RADI step:   39 pc: -1.006090e+01 + 0.000000e+00i normalized residual: 2.116580e-07 relative change in Z: 1.349319e-06
RADI step:   41 pc: -1.319288e+03 + -6.671412e+03i normalized residual: 1.285451e-07 relative change in Z: 8.883680e-06
RADI step:   43 pc: -1.224419e+03 + -3.118088e+03i normalized residual: 6.423242e-09 relative change in Z: 1.133498e-05
RADI step:   44 pc: -2.626777e+03 + 0.000000e+00i normalized residual: 4.906290e-09 relative change in Z: 3.409488e-06
RADI step:   45 pc: -1.455835e+01 + 0.000000e+00i normalized residual: 4.810056e-09 relative change in Z: 4.333556e-07
RADI step:   47 pc: -1.154675e+03 + -7.261498e+03i normalized residual: 3.292168e-09 relative change in Z: 1.656312e-06
RADI step:   49 pc: -1.068699e+03 + 2.897627e+03i normalized residual: 4.663999e-10 relative change in Z: 2.159784e-06
RADI step:   51 pc: -3.389061e+03 + -4.327642e+02i normalized residual: 1.207785e-10 relative change in Z: 8.067382e-07
RADI step:   53 pc: -1.282176e+03 + -8.496731e+03i normalized residual: 1.921659e-11 relative change in Z: 3.357523e-07

Size outFil.Z: 905 x 73
solving the filter Riccati equation took   0.34 seconds 
Residual computations -- RADI: 1.921659e-11 | mess_res2_norms: 1.107772e-10 | eigs: 1.107844e-10