Contents

function LQR_DAE1(istest)
% Computes a Riccati feedback control for the proper
% index-1 System BIPS98_606 from https://sites.google.com/site/rommes/software
% following the ideas introduced in [1] for Lyapunov equations using he
% Newton-ADI iteration.
%
% Input:
% istest  decides whether the function runs as an interactive demo or a
%         continuous integration test. (optional; defaults to 0, i.e.
%         interactive demo)
%
% References:
%[1] F. Freitas, J. Rommes, N. Martins, Gramian-based reduction method
%    applied to large sparse power system descriptor models, IEEE Trans.
%    Power Syst. 23 (3) (2008) 1258–1270. doi:10.1109/TPWRS.2008.926693

%
% 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)
%
if nargin<1, istest=0; end

set operation

oper = operatormanager('dae_1');

Problem data

eqn = mess_get_BIPS(7);

Turn off close to singular warnings

(this model is really badly conditioned)
orig_warnstate = warning('OFF','MATLAB:nearlySingularMatrix');
opts.norm = 'fro';

% 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 = 0;
opts.adi.projection.freq=0;

eqn.type='N';
opts.shifts.num_desired=25;
opts.shifts.num_Ritz=50;
opts.shifts.num_hRitz=25;
opts.shifts.method = 'projection';
opts.shifts.num_desired= 9;

Newton tolerances and maximum iteration number

opts.nm.maxiter = 30;
opts.nm.res_tol = 1e-10;
opts.nm.rel_diff_tol = 1e-16;
opts.nm.info = 1;
opts.nm.linesearch = 1;
opts.nm.accumulateRes = 1;
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 istest
    if min(outnm.res)>=opts.nm.res_tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result in LRNM');
    end
else
    figure();
    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 iterations');
    ylabel('normalized residual norm');
    pause(1);
end
disp('size outnm.Z:');
disp(size(outnm.Z));
		 Using line search (res: 2.838281e-01)
		 lambda: 8.784350e-01

 NM step:    1  normalized residual: 2.420598e-01
               relative change in K: 1.000000e+00
               number of ADI steps: 187 


		 Using line search (res: 5.458145e-02)
		 lambda: 1.201723e+00

 NM step:    2  normalized residual: 4.363637e-02
               relative change in K: 3.145467e+00
               number of ADI steps: 207 


 NM step:    3  normalized residual: 8.550157e-03
               relative change in K: 9.423954e-01
               number of ADI steps: 222 


 NM step:    4  normalized residual: 2.080652e-03
               relative change in K: 7.739811e-01
               number of ADI steps: 236 


 NM step:    5  normalized residual: 5.951448e-04
               relative change in K: 4.726820e-01
               number of ADI steps: 192 


 NM step:    6  normalized residual: 5.145986e-04
               relative change in K: 3.286749e-01
               number of ADI steps: 188 


		 Using line search (res: 1.037793e-03)
		 lambda: 2.717878e-01

 NM step:    7  normalized residual: 4.430466e-04
               relative change in K: 2.805264e-01
               number of ADI steps: 194 


		 Using line search (res: 9.712399e-04)
		 lambda: 2.456325e-01

 NM step:    8  normalized residual: 3.877868e-04
               relative change in K: 2.440394e-01
               number of ADI steps: 217 


		 Using line search (res: 8.876159e-04)
		 lambda: 2.317425e-01

 NM step:    9  normalized residual: 3.423524e-04
               relative change in K: 2.149219e-01
               number of ADI steps: 207 


		 Using line search (res: 7.917199e-04)
		 lambda: 2.265002e-01

 NM step:   10  normalized residual: 3.032567e-04
               relative change in K: 1.907198e-01
               number of ADI steps: 222 


		 Using line search (res: 6.873039e-04)
		 lambda: 2.286794e-01

 NM step:   11  normalized residual: 2.683577e-04
               relative change in K: 1.698967e-01
               number of ADI steps: 221 


		 Using line search (res: 5.774723e-04)
		 lambda: 2.387039e-01

 NM step:   12  normalized residual: 2.361637e-04
               relative change in K: 1.513977e-01
               number of ADI steps: 170 


		 Using line search (res: 4.649922e-04)
		 lambda: 2.589076e-01

 NM step:   13  normalized residual: 2.054636e-04
               relative change in K: 1.344069e-01
               number of ADI steps: 186 


		 Using line search (res: 3.525780e-04)
		 lambda: 2.951567e-01

 NM step:   14  normalized residual: 1.750367e-04
               relative change in K: 1.181660e-01
               number of ADI steps: 175 


		 Using line search (res: 2.432579e-04)
		 lambda: 3.624454e-01

 NM step:   15  normalized residual: 1.432222e-04
               relative change in K: 1.017240e-01
               number of ADI steps: 158 


 NM step:   16  normalized residual: 1.410808e-04
               relative change in K: 8.932130e-02
               number of ADI steps: 147 


 NM step:   17  normalized residual: 5.011963e-06
               relative change in K: 1.902161e-02
               number of ADI steps: 128 


 NM step:   18  normalized residual: 4.131441e-07
               relative change in K: 5.234900e-03
               number of ADI steps: 142 


 NM step:   19  normalized residual: 7.625096e-08
               relative change in K: 2.143991e-03
               number of ADI steps: 131 


 NM step:   20  normalized residual: 5.860942e-09
               relative change in K: 6.077827e-04
               number of ADI steps: 142 


 NM step:   21  normalized residual: 7.850643e-11
               relative change in K: 6.660950e-05
               number of ADI steps: 84 

mess_lrnm took 152.71 seconds 
size outnm.Z:
   606   118

Lets try RADI

opts.norm = 2;

% RADI-MESS settings
opts.shifts.history = opts.shifts.num_desired*size(eqn.C,1);
opts.shifts.method  = 'gen-ham-opti';

opts.shifts.naive_update_mode = false;
% .. Suggest false (smart update is faster; convergence is the same).

opts.radi.compute_sol_fac = 1;
opts.radi.get_ZZt = 1;
opts.radi.compute_res = 0;
opts.radi.maxiter = 500;
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 istest
    if min(outradi.res)>=opts.radi.res_tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result in RADI');
    end
else
    figure();
    semilogy(outradi.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 iterations');
    ylabel('normalized residual norm');
end
disp('size outradi.Z:');
disp(size(outradi.Z));
RADI step:    1 normalized residual: 9.571872e-03
RADI step:    2 normalized residual: 1.049748e-02
RADI step:    4 normalized residual: 1.037405e-03
RADI step:    5 normalized residual: 9.851085e-04
RADI step:    6 normalized residual: 9.321639e-04
RADI step:    7 normalized residual: 8.933725e-04
RADI step:    8 normalized residual: 2.603056e-06
RADI step:    9 normalized residual: 6.172319e-06
RADI step:   10 normalized residual: 5.845786e-06
RADI step:   12 normalized residual: 2.321815e-06
RADI step:   14 normalized residual: 2.331021e-06
RADI step:   15 normalized residual: 2.457001e-06
RADI step:   17 normalized residual: 3.033602e-06
RADI step:   19 normalized residual: 3.013795e-06
RADI step:   21 normalized residual: 3.089809e-06
RADI step:   23 normalized residual: 2.000302e-05
RADI step:   25 normalized residual: 2.046737e-05
RADI step:   26 normalized residual: 2.078825e-05
RADI step:   27 normalized residual: 2.266239e-05
RADI step:   28 normalized residual: 2.410377e-05
RADI step:   30 normalized residual: 2.385350e-05
RADI step:   32 normalized residual: 1.333943e-07
RADI step:   33 normalized residual: 1.329079e-07
RADI step:   35 normalized residual: 1.573517e-07
RADI step:   37 normalized residual: 1.590843e-07
RADI step:   38 normalized residual: 1.601255e-07
RADI step:   40 normalized residual: 1.626163e-07
RADI step:   41 normalized residual: 1.270573e-07
RADI step:   43 normalized residual: 1.273628e-07
RADI step:   44 normalized residual: 1.143198e-07
RADI step:   46 normalized residual: 1.153163e-07
RADI step:   48 normalized residual: 1.028228e-07
RADI step:   50 normalized residual: 9.267039e-08
RADI step:   51 normalized residual: 9.228848e-08
RADI step:   53 normalized residual: 9.575965e-08
RADI step:   54 normalized residual: 8.869360e-08
RADI step:   56 normalized residual: 8.791654e-08
RADI step:   58 normalized residual: 8.459751e-08
RADI step:   59 normalized residual: 7.556533e-09
RADI step:   60 normalized residual: 1.876330e-08
RADI step:   61 normalized residual: 1.893301e-08
RADI step:   63 normalized residual: 1.933470e-08
RADI step:   65 normalized residual: 1.948601e-08
RADI step:   67 normalized residual: 1.981708e-08
RADI step:   68 normalized residual: 2.038561e-08
RADI step:   70 normalized residual: 1.963424e-08
RADI step:   71 normalized residual: 2.018131e-08
RADI step:   73 normalized residual: 2.065537e-08
RADI step:   75 normalized residual: 2.061558e-08
RADI step:   77 normalized residual: 2.042812e-08
RADI step:   79 normalized residual: 2.067530e-08
RADI step:   81 normalized residual: 2.066508e-08
RADI step:   83 normalized residual: 2.066456e-08
RADI step:   84 normalized residual: 2.058455e-08
RADI step:   86 normalized residual: 5.872939e-09
RADI step:   87 normalized residual: 5.819833e-09
RADI step:   88 normalized residual: 5.812557e-09
RADI step:   90 normalized residual: 5.240826e-09
RADI step:   91 normalized residual: 5.158949e-09
RADI step:   92 normalized residual: 5.017203e-09
RADI step:   94 normalized residual: 4.759440e-09
RADI step:   96 normalized residual: 4.671644e-09
RADI step:   98 normalized residual: 4.751987e-09
RADI step:   99 normalized residual: 4.598024e-09
RADI step:  100 normalized residual: 4.494784e-09
RADI step:  101 normalized residual: 4.313083e-09
RADI step:  103 normalized residual: 3.744081e-09
RADI step:  104 normalized residual: 1.320733e-10
RADI step:  105 normalized residual: 1.320634e-10
RADI step:  107 normalized residual: 1.349657e-10
RADI step:  109 normalized residual: 1.358398e-10
RADI step:  110 normalized residual: 1.417912e-10
RADI step:  112 normalized residual: 1.328264e-10
RADI step:  113 normalized residual: 1.249755e-10
RADI step:  114 normalized residual: 1.114197e-10
RADI step:  116 normalized residual: 1.393720e-10
RADI step:  117 normalized residual: 1.396652e-10
RADI step:  118 normalized residual: 1.481514e-10
RADI step:  119 normalized residual: 1.492920e-10
RADI step:  121 normalized residual: 1.703151e-10
RADI step:  122 normalized residual: 1.878151e-10
RADI step:  123 normalized residual: 1.911557e-10
RADI step:  125 normalized residual: 1.906492e-10
RADI step:  127 normalized residual: 1.864185e-10
RADI step:  128 normalized residual: 1.990756e-10
RADI step:  129 normalized residual: 2.219637e-10
RADI step:  131 normalized residual: 2.225979e-10
RADI step:  133 normalized residual: 2.237630e-10
RADI step:  134 normalized residual: 2.281696e-10
RADI step:  136 normalized residual: 2.279055e-10
RADI step:  138 normalized residual: 2.280573e-10
RADI step:  140 normalized residual: 1.452224e-10
RADI step:  142 normalized residual: 1.476485e-10
RADI step:  143 normalized residual: 1.469076e-10
RADI step:  145 normalized residual: 1.370782e-10
RADI step:  147 normalized residual: 1.301370e-10
RADI step:  148 normalized residual: 1.308642e-10
RADI step:  149 normalized residual: 1.309574e-10
RADI step:  151 normalized residual: 1.301693e-10
RADI step:  152 normalized residual: 1.256852e-10
RADI step:  154 normalized residual: 1.191989e-10
RADI step:  156 normalized residual: 1.191493e-10
RADI step:  158 normalized residual: 1.198222e-10
RADI step:  159 normalized residual: 1.186296e-10
RADI step:  160 normalized residual: 1.178697e-10
RADI step:  161 normalized residual: 1.175573e-10
RADI step:  162 normalized residual: 1.173704e-10
RADI step:  164 normalized residual: 1.169534e-10
RADI step:  166 normalized residual: 1.169540e-10
RADI step:  168 normalized residual: 1.105529e-10
RADI step:  170 normalized residual: 3.476434e-12
mess_lrradi took  10.58 seconds 
size outradi.Z:
   606   120

compare

if not(istest)
    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 * E');
    ylabel('normalized residual norm');
    legend('LR-NM','RADI');
end

reset warning state

warning(orig_warnstate');