Contents

function LQG_FDM_unstable_nwt(n0, n_unstable, istest)
% Computes stabilizing and detecting solutions of the regulator and filter
% Riccati equations via the Newton-ADI method, 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_unstable, n_stable), eye(n_unstable)];
KB0 = Bt' * Y * Vl;
if exist('lyap','file') % Solve Lyapunov equation.
    Au = lyap(Y, -Bt * Bt');
else
    Au = lyap2solve(Y, -Bt * Bt');
end
KC0 = C(:, n_stable+1:n) * Y * Vl;

% 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';

ADI tolerances and maximum iteration number

opts.adi.maxiter          = 200;
opts.adi.res_tol          = 1e-10;
opts.adi.rel_diff_tol     = 1e-16;
opts.adi.info             = 0;
opts.adi.accumulateDeltaK = 1;
opts.adi.accumulateK      = 0;
opts.adi.compute_sol_fac  = 1;

shift parameters via projection

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

Newton tolerances and maximum iteration number

opts.nm.maxiter       = 15;
opts.nm.res_tol       = 1e-8;
opts.nm.rel_diff_tol  = 1e-16;
opts.nm.info          = 1;
opts.nm.linesearch    = 1;
opts.nm.accumulateRes = 1;
opts.nm.res           = struct('maxiter',10,'tol',1e-6,'orth',0);

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.nm.K0 = KB0;

time_nwt_reguator = tic;
[outReg, eqn, opts, oper]=mess_lrnm(eqn, opts, oper);
t_elapsed1 = toc(time_nwt_reguator);

% 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.nm, []) / 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 -- Newton: %e | ' ...
    'mess_res2_norms: %e | eigs: %e \n'], ...
    outReg.res(end), res1, res2);

% Print convergence behavior.
if istest
    if min(outReg.res) >= opts.nm.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  BB^T X');
    xlabel('number of iteration steps');
    ylabel('normalized residual norm');
    pause(1);
end

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

 NM step:    1  normalized residual: 9.444737e-01
               relative change in K: 4.999150e-01
               number of ADI steps: 20 


 NM step:    2  normalized residual: 1.049167e-01
               relative change in K: 1.999310e-01
               number of ADI steps: 21 


 NM step:    3  normalized residual: 2.139893e-03
               relative change in K: 2.939241e-02
               number of ADI steps: 20 


 NM step:    4  normalized residual: 9.675242e-07
               relative change in K: 6.253776e-04
               number of ADI steps: 25 


 NM step:    5  normalized residual: 1.057951e-09
               relative change in K: 2.830132e-07
               number of ADI steps: 19 


Size outReg.Z: 905 x 32
solving the regulator Riccati equation took   0.69 seconds 
Residual computation -- Newton: 1.057951e-09 | mess_res2_norms: 8.313153e-10 | eigs: 8.313153e-10 

Solve filter Riccati equation.

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

eqn.type   = 'N';
opts.LDL_T = 1;
opts.nm.K0 = KC0;

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

time_nwt_filter = tic;
[outFil, eqn, opts, oper]=mess_lrnm(eqn, opts, oper);
t_elapsed2 = toc(time_nwt_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.nm, 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 -- Newton: %e | ' ...
    'mess_res2_norms: %e | eigs: %e \n'], ...
    outFil.res(end), res3, res4);

% Print convergence behavior.
if istest
    if min(outFil.res) >= opts.nm.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
-----------------------------

 NM step:    1  normalized residual: 9.453716e-01
               relative change in K: 4.999209e-01
               number of ADI steps: 59 


 NM step:    2  normalized residual: 1.050126e-01
               relative change in K: 1.999301e-01
               number of ADI steps: 50 


 NM step:    3  normalized residual: 2.141849e-03
               relative change in K: 2.939225e-02
               number of ADI steps: 51 


 NM step:    4  normalized residual: 9.684086e-07
               relative change in K: 6.253743e-04
               number of ADI steps: 62 


 NM step:    5  normalized residual: 9.470465e-09
               relative change in K: 2.830117e-07
               number of ADI steps: 41 


Size outFil.Z: 905 x 492
solving the filter Riccati equation took   1.66 seconds 
Residual computations -- Newton: 9.470465e-09 | mess_res2_norms: 1.186182e-08 | eigs: 1.186182e-08