Contents

function LQR_DAE2(problem, level, re, istest)
% Computes a stabilizing feedback by implicitly solving the generalized
% Riccati equation for the equivalent projected system on the hidden manifold.
%
% Inputs:
% problem       either 'Stokes' or 'NSE' to choose the Stokes demo or the
%               linearized Navier-Stokes-Equation.
%               (optional, defaults to 'Stokes')
%
% level         discretization level 1 through 5
%               (optional, only used in 'NSE' case, default: 1)
%
% re            Reynolds number 300, 400, or 500
%               (optional, only used in 'NSE' case, default: 500)
%
% istest        flag to determine whether this demo runs as a CI test or
%               interactive demo
%               (optional, defaults to 0, i.e. interactive demo)
%
% Note that the 'NSE' option requires additional data available in a
% separate 270MB archive and at least the 5th discretization level needs a
% considerable amount of main memory installed in your machine.

%
% 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)
%

Set operations

oper = operatormanager('dae_2');

Problem data

if nargin<1, problem='stokes'; end
if nargin<2, level=1; end
if nargin<3, re=500; end
if nargin<4, istest=0; end

switch lower(problem)
    case 'stokes'
        nin = 5;
        nout = 5;
        nx = 10;
        ny = 10;
        [eqn.E_,eqn.A_,eqn.B,eqn.C]=stokes_ind2(nin,nout,nx,ny);
        eqn.haveE=1;
        st=trace(eqn.E_); % Stokes is FDM discretized, so so this is
        % the dimension of the velocity space
        eqn.st=st;
        eqn.B=eqn.B(1:st,:);
        eqn.C=eqn.C(:,1:st);
    case 'nse'
        [eqn, K0, ~] = mess_get_NSE(re, level);
        opts.nm.K0 = K0;
        opts.radi.K0 = K0;
    otherwise
        error('input ''problem'' must be either ''NSE'' or ''Stokes''');
end
degrees of freedom: 
------------------------------
 total         :   280
 velocity      :   180
 pressure      :   100
 n_finite      :    81
------------------------------
Generating FVM matrices...
 -> Laplacians...
 -> gradient and divergence operator...
 -> B1, C1 and v1 velocity nodes...
 -> B2, C2 and v2 velocity nodes...
Setting up system matrices ...

First we run the Newton-ADI Method

opts.norm = 2;

% ADI tolerances and maximum iteration number
opts.adi.maxiter = 300;
opts.adi.res_tol = 1e-12;
opts.adi.rel_diff_tol = 1e-16;
opts.adi.info = 1;
opts.adi.LDL_T=0;
eqn.type='T';
n=size(eqn.A_, 1);
opts.shifts.num_desired=5;%*nout;
opts.shifts.num_Ritz=50;
opts.shifts.num_hRitz=25;
opts.shifts.method = 'projection';
opts.shifts.b0=ones(n,1);

Newton tolerances and maximum iteration number

opts.nm.maxiter = 20;
opts.nm.res_tol = 1e-10;
opts.nm.rel_diff_tol = 1e-16;
opts.nm.info = 1;
opts.nm.projection.freq=0;
opts.nm.projection.ortho=1;
%opts.nm.projection.meth='care_nwt_fac';
opts.nm.res=struct('maxiter',10,'tol',1e-6,'orth',0);
opts.nm.linesearch = 1;
opts.nm.inexact = 'superlinear';
opts.nm.tau = 0.1;
opts.nm.accumulateRes = 1;

use low-rank Newton-Kleinman-ADI

t_mess_lrnm = tic;
outnm = mess_lrnm(eqn, opts, oper);
t_elapsed1 =toc(t_mess_lrnm);
fprintf(1,'mess_lrnm took %6.2f seconds \n' , t_elapsed1);
if not(istest)
    figure(1);
    disp(outnm.res);
    semilogy(outnm.res,'LineWidth',3);
    title('0 = C^T C + A^T X E + E^T X A - E^T X BB^T X E');
    xlabel('number of newton iterations');
    ylabel('normalized residual norm');
    pause(1);
end
disp('size outnm.Z:');
disp(size(outnm.Z));
ADI step:    1 normalized residual: 4.749604e-01 relative change in Z: 1.000000e+00
		 normalized outer residual: 4.748502e-01
ADI step:    2 normalized residual: 2.313480e-01 relative change in Z: 6.115557e-01
		 normalized outer residual: 2.309720e-01
ADI step:    3 normalized residual: 1.160895e-01 relative change in Z: 4.086903e-01
		 normalized outer residual: 1.154641e-01
ADI step:    4 normalized residual: 3.274932e-02 relative change in Z: 3.384842e-01
		 normalized outer residual: 3.186747e-02
ADI step:    5 normalized residual: 4.408328e-03 relative change in Z: 2.031583e-01
		 normalized outer residual: 3.477590e-03

 outer tolerance: 1.188528e-02

 NM step:    1  normalized residual: 3.477590e-03
               relative change in K: 1.000000e+00
               number of ADI steps: 5 

ADI step:    1 normalized residual: 4.751369e-01 relative change in Z: 1.000000e+00
		 normalized outer residual: 4.746418e-01
ADI step:    2 normalized residual: 2.310230e-01 relative change in Z: 6.109172e-01
		 normalized outer residual: 2.308408e-01
ADI step:    3 normalized residual: 1.156490e-01 relative change in Z: 4.073025e-01
		 normalized outer residual: 1.155924e-01
ADI step:    4 normalized residual: 3.274385e-02 relative change in Z: 3.354002e-01
		 normalized outer residual: 3.273882e-02
ADI step:    5 normalized residual: 4.365684e-03 relative change in Z: 2.003457e-01
		 normalized outer residual: 4.365677e-03
ADI step:    6 normalized residual: 1.459723e-03 relative change in Z: 6.921308e-02
		 normalized outer residual: 1.459682e-03
ADI step:    7 normalized residual: 9.263155e-05 relative change in Z: 5.492726e-02
		 normalized outer residual: 9.253389e-05

 outer tolerance: 3.863989e-04

 NM step:    2  normalized residual: 9.253389e-05
               relative change in K: 1.322216e-02
               number of ADI steps: 7 

ADI step:    1 normalized residual: 5.908896e-01 relative change in Z: 1.000000e+00
		 normalized outer residual: 5.902469e-01
ADI step:    2 normalized residual: 2.801793e-01 relative change in Z: 6.987959e-01
		 normalized outer residual: 2.799246e-01
ADI step:    3 normalized residual: 7.784030e-02 relative change in Z: 5.237500e-01
		 normalized outer residual: 7.780849e-02
ADI step:    4 normalized residual: 6.236314e-04 relative change in Z: 3.257720e-01
		 normalized outer residual: 6.234901e-04
ADI step:    5 normalized residual: 5.863780e-05 relative change in Z: 3.252883e-02
		 normalized outer residual: 5.863776e-05
ADI step:    6 normalized residual: 2.832583e-05 relative change in Z: 5.135531e-03
		 normalized outer residual: 2.832582e-05
ADI step:    7 normalized residual: 2.435343e-05 relative change in Z: 4.159017e-03
		 normalized outer residual: 2.435343e-05
ADI step:    8 normalized residual: 2.563349e-07 relative change in Z: 1.651384e-03
		 normalized outer residual: 2.563341e-07

 outer tolerance: 3.304782e-06

 NM step:    3  normalized residual: 2.563341e-07
               relative change in K: 1.145651e-04
               number of ADI steps: 8 

ADI step:    1 normalized residual: 5.952306e-01 relative change in Z: 1.000000e+00
		 normalized outer residual: 5.945829e-01
ADI step:    2 normalized residual: 3.403191e-01 relative change in Z: 6.609782e-01
		 normalized outer residual: 3.399864e-01
ADI step:    3 normalized residual: 1.078645e-01 relative change in Z: 5.650374e-01
		 normalized outer residual: 1.078086e-01
ADI step:    4 normalized residual: 2.082073e-02 relative change in Z: 3.428347e-01
		 normalized outer residual: 2.081777e-02
ADI step:    5 normalized residual: 1.759708e-04 relative change in Z: 1.767522e-01
		 normalized outer residual: 1.759642e-04
ADI step:    6 normalized residual: 8.406516e-05 relative change in Z: 1.521109e-02
		 normalized outer residual: 8.406342e-05
ADI step:    7 normalized residual: 9.211191e-06 relative change in Z: 1.267083e-02
		 normalized outer residual: 9.211191e-06
ADI step:    8 normalized residual: 3.101855e-06 relative change in Z: 1.302219e-03
		 normalized outer residual: 3.101855e-06
ADI step:    9 normalized residual: 2.053260e-08 relative change in Z: 1.875062e-03
		 normalized outer residual: 2.053260e-08
ADI step:   10 normalized residual: 6.170189e-09 relative change in Z: 8.933796e-05
		 normalized outer residual: 6.170189e-09
ADI step:   11 normalized residual: 1.864742e-09 relative change in Z: 6.988243e-05
		 normalized outer residual: 1.864742e-09
ADI step:   12 normalized residual: 1.354435e-09 relative change in Z: 2.042627e-05
		 normalized outer residual: 1.354435e-09
ADI step:   13 normalized residual: 9.506928e-11 relative change in Z: 4.068726e-05
		 normalized outer residual: 9.506928e-11

 outer tolerance: 3.943602e-09

 NM step:    4  normalized residual: 9.506928e-11
               relative change in K: 1.797500e-07
               number of ADI steps: 13 

mess_lrnm took   0.50 seconds 
    0.0035
    0.0001
    0.0000
    0.0000

size outnm.Z:
   180    54

Lets try the RADI method and compare

opts.norm = 2;
% RADI-MESS settings
opts.shifts.history = opts.shifts.num_desired*size(eqn.C,1);
opts.shifts.num_desired = opts.shifts.num_desired;

% choose either of the three shift methods, here
opts.shifts.method = 'gen-ham-opti';
%     opts.shifts.method = 'heur';
%     opts.shifts.method = 'projection';

opts.shifts.naive_update_mode = false; % .. Suggest false (smart update is faster; convergence is the same).
opts.shifts.info = 0;
opts.radi.compute_sol_fac = 1; % Turned on for numerical stability reasons.
opts.radi.get_ZZt = 0;
opts.radi.maxiter = opts.adi.maxiter;
opts.radi.res_tol = opts.nm.res_tol;
opts.radi.rel_diff_tol = 0;
opts.radi.info = 1;


t_mess_lrradi = tic;
outradi = mess_lrradi(eqn, opts, oper);
t_elapsed2 = toc(t_mess_lrradi);
fprintf(1,'mess_lrradi took %6.2f seconds \n' ,t_elapsed2);
if not(istest)
    figure();
    semilogy(outradi.res,'LineWidth',3);
    title('0 = C^TC + A^T X E + E^T X A - E^T X BB^T X E');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
end
RADI step:    1 normalized residual: 7.035353e-02
RADI step:    2 normalized residual: 1.197302e-02
RADI step:    3 normalized residual: 3.148166e-03
RADI step:    4 normalized residual: 2.370135e-03
RADI step:    5 normalized residual: 6.534372e-05
RADI step:    6 normalized residual: 5.589811e-05
RADI step:    7 normalized residual: 2.820919e-06
RADI step:    8 normalized residual: 2.015254e-07
RADI step:    9 normalized residual: 1.213884e-07
RADI step:   10 normalized residual: 7.433566e-09
RADI step:   11 normalized residual: 4.192365e-09
RADI step:   12 normalized residual: 7.152029e-10
RADI step:   13 normalized residual: 3.266705e-11
mess_lrradi took   0.10 seconds 

compare

if istest
    if min(outnm.res)>=opts.nm.res_tol, error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end
    if min(outradi.res)>=opts.radi.res_tol, error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end
else
    figure();
    ls_nm=[outnm.adi.niter];
    ls_radi=1:outradi.niter;

    semilogy(cumsum(ls_nm),outnm.res,'k--',ls_radi,outradi.res,'b-','LineWidth',3);
    title('0 = C^T C + A^T X E + E^T X A - E^T X BB^T X E');
    xlabel('number of solves with A+p*M');
    ylabel('normalized residual norm');
    legend('LR-NM','RADI');
end