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 program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
%               2009-2020
%

if nargin<1, istest=0; end

set operation

oper = operatormanager('dae_1');

Problem data

fname  = sprintf('%s/../models/BIPS/bips98_606.mat',...
    fileparts(mfilename('fullpath')));
Bips = load(fname);
% from https://sites.google.com/site/rommes/software
p = find(diag(Bips.E));
np = find(diag(Bips.E) == 0);
pp = [p;np];
eqn.A_ = Bips.A(pp, pp);
eqn.E_ = Bips.E(pp, pp);
eqn.B = Bips.b(pp, :);
eqn.C = 0.01*Bips.c( : , pp);
eqn.st = length(p);
eqn.haveE = 1;
clear Bips;

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;

tic;
outnm = mess_lrnm(eqn, opts, oper);
toc;
if istest
    if min(outnm.res)>=opts.nm.res_tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result in LRNM');
    end
else
    figure(1);
    semilogy(outnm.res);
    title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
disp('size outnm.Z:');
disp(size(outnm.Z));

		Using line search (res: 6.145277e+09)
		lambda: 9.490913e-06

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


NM step:    2  normalized residual: 2.947417e-04
               relative change in K: 2.468724e-02
               number of ADI steps: 268 


NM step:    3  normalized residual: 7.562885e-05
               relative change in K: 1.208700e-02
               number of ADI steps: 248 


NM step:    4  normalized residual: 2.632597e-06
               relative change in K: 2.230511e-03
               number of ADI steps: 226 


NM step:    5  normalized residual: 1.265741e-09
               relative change in K: 5.006001e-05
               number of ADI steps: 300 


NM step:    6  normalized residual: 1.506303e-11
               relative change in K: 5.134368e-08
               number of ADI steps: 239 

Elapsed time is 44.423862 seconds.
size outnm.Z:
   606    72

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;

tic;
outradi = mess_lrradi( eqn, opts, oper );
toc;

if istest
    if min(outradi.res)>=opts.radi.res_tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result in RADI');
    end
else
    figure(2);
    semilogy(outradi.res);
    title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
end
disp('size outradi.Z:');
disp(size(outradi.Z));

RADI step:    1 normalized residual: 9.599173e-03
RADI step:    2 normalized residual: 1.908370e-02
RADI step:    4 normalized residual: 8.986170e-04
RADI step:    6 normalized residual: 8.511832e-05
RADI step:    7 normalized residual: 8.569602e-05
RADI step:    9 normalized residual: 2.228144e-04
RADI step:   11 normalized residual: 1.448097e-04
RADI step:   13 normalized residual: 1.920097e-04
RADI step:   14 normalized residual: 2.195916e-04
RADI step:   15 normalized residual: 1.554772e-03
RADI step:   16 normalized residual: 4.855531e-03
RADI step:   18 normalized residual: 7.028810e-03
RADI step:   19 normalized residual: 2.203283e-02
RADI step:   20 normalized residual: 8.762692e-03
RADI step:   21 normalized residual: 1.908620e-03
RADI step:   22 normalized residual: 4.711082e-04
RADI step:   23 normalized residual: 4.198868e-05
RADI step:   25 normalized residual: 5.179948e-04
RADI step:   26 normalized residual: 3.313485e-04
RADI step:   27 normalized residual: 1.135597e-04
RADI step:   29 normalized residual: 1.459155e-04
RADI step:   31 normalized residual: 1.098862e-03
RADI step:   33 normalized residual: 3.098942e-03
RADI step:   35 normalized residual: 2.206764e-03
RADI step:   36 normalized residual: 1.962263e-03
RADI step:   37 normalized residual: 1.115376e-03
RADI step:   38 normalized residual: 3.425401e-04
RADI step:   40 normalized residual: 2.202839e-04
RADI step:   42 normalized residual: 2.198304e-04
RADI step:   44 normalized residual: 3.992017e-05
RADI step:   45 normalized residual: 3.454886e-05
RADI step:   46 normalized residual: 3.612336e-05
RADI step:   47 normalized residual: 9.800702e-06
RADI step:   48 normalized residual: 1.082880e-05
RADI step:   49 normalized residual: 1.396232e-05
RADI step:   51 normalized residual: 1.403041e-05
RADI step:   52 normalized residual: 1.487255e-05
RADI step:   53 normalized residual: 1.178772e-05
RADI step:   54 normalized residual: 6.166076e-05
RADI step:   56 normalized residual: 6.625390e-05
RADI step:   57 normalized residual: 6.704507e-05
RADI step:   59 normalized residual: 6.477090e-05
RADI step:   60 normalized residual: 8.393733e-05
RADI step:   61 normalized residual: 4.385472e-05
RADI step:   63 normalized residual: 6.373227e-05
RADI step:   64 normalized residual: 1.319460e-04
RADI step:   66 normalized residual: 8.949816e-05
RADI step:   67 normalized residual: 9.625756e-05
RADI step:   69 normalized residual: 1.967213e-05
RADI step:   71 normalized residual: 7.661990e-06
RADI step:   72 normalized residual: 7.481344e-06
RADI step:   73 normalized residual: 9.008916e-06
RADI step:   74 normalized residual: 1.013113e-05
RADI step:   76 normalized residual: 3.858499e-06
RADI step:   77 normalized residual: 2.870654e-06
RADI step:   78 normalized residual: 2.723241e-06
RADI step:   79 normalized residual: 2.687573e-06
RADI step:   80 normalized residual: 9.181773e-06
RADI step:   82 normalized residual: 1.054376e-05
RADI step:   84 normalized residual: 6.527040e-06
RADI step:   86 normalized residual: 2.379435e-06
RADI step:   88 normalized residual: 1.562171e-06
RADI step:   89 normalized residual: 1.680017e-06
RADI step:   91 normalized residual: 1.910162e-06
RADI step:   93 normalized residual: 2.612364e-06
RADI step:   95 normalized residual: 3.892949e-06
RADI step:   97 normalized residual: 2.534717e-06
RADI step:   99 normalized residual: 6.643923e-07
RADI step:  101 normalized residual: 5.481648e-07
RADI step:  102 normalized residual: 5.242353e-07
RADI step:  103 normalized residual: 5.477457e-07
RADI step:  104 normalized residual: 5.612204e-07
RADI step:  106 normalized residual: 5.153755e-07
RADI step:  107 normalized residual: 5.519777e-07
RADI step:  109 normalized residual: 3.818941e-07
RADI step:  110 normalized residual: 3.818050e-07
RADI step:  112 normalized residual: 3.814077e-07
RADI step:  113 normalized residual: 3.814069e-07
RADI step:  114 normalized residual: 3.786058e-07
RADI step:  116 normalized residual: 3.677029e-07
RADI step:  117 normalized residual: 3.542040e-07
RADI step:  119 normalized residual: 3.084748e-07
RADI step:  121 normalized residual: 3.005680e-07
RADI step:  122 normalized residual: 3.001007e-07
RADI step:  124 normalized residual: 2.845097e-07
RADI step:  126 normalized residual: 2.841175e-07
RADI step:  127 normalized residual: 2.837284e-07
RADI step:  129 normalized residual: 2.794714e-07
RADI step:  131 normalized residual: 2.770217e-07
RADI step:  132 normalized residual: 2.459212e-07
RADI step:  134 normalized residual: 2.366464e-07
RADI step:  136 normalized residual: 2.314890e-07
RADI step:  138 normalized residual: 2.293491e-07
RADI step:  139 normalized residual: 2.264931e-07
RADI step:  141 normalized residual: 2.257476e-07
RADI step:  143 normalized residual: 2.245017e-07
RADI step:  145 normalized residual: 2.242471e-07
RADI step:  146 normalized residual: 2.155094e-07
RADI step:  148 normalized residual: 2.135228e-07
RADI step:  149 normalized residual: 2.121524e-07
RADI step:  150 normalized residual: 2.120637e-07
RADI step:  152 normalized residual: 2.089365e-07
RADI step:  153 normalized residual: 1.983920e-07
RADI step:  155 normalized residual: 1.856248e-07
RADI step:  156 normalized residual: 1.855993e-07
RADI step:  158 normalized residual: 1.701099e-07
RADI step:  159 normalized residual: 1.680616e-07
RADI step:  161 normalized residual: 1.653339e-07
RADI step:  162 normalized residual: 1.567953e-07
RADI step:  163 normalized residual: 1.501578e-07
RADI step:  164 normalized residual: 1.497228e-07
RADI step:  166 normalized residual: 1.452866e-07
RADI step:  167 normalized residual: 1.112695e-09
RADI step:  168 normalized residual: 1.111093e-09
RADI step:  170 normalized residual: 1.116877e-09
RADI step:  171 normalized residual: 9.473453e-10
RADI step:  172 normalized residual: 9.228011e-10
RADI step:  174 normalized residual: 7.391312e-10
RADI step:  176 normalized residual: 7.057959e-10
RADI step:  177 normalized residual: 6.803305e-10
RADI step:  179 normalized residual: 6.646255e-10
RADI step:  181 normalized residual: 6.537771e-10
RADI step:  183 normalized residual: 6.300925e-10
RADI step:  185 normalized residual: 6.270189e-10
RADI step:  187 normalized residual: 6.057439e-10
RADI step:  189 normalized residual: 5.826613e-10
RADI step:  191 normalized residual: 5.595750e-10
RADI step:  192 normalized residual: 5.479129e-10
RADI step:  194 normalized residual: 5.435029e-10
RADI step:  195 normalized residual: 5.334675e-10
RADI step:  197 normalized residual: 5.313044e-10
RADI step:  199 normalized residual: 4.924531e-10
RADI step:  201 normalized residual: 4.729925e-10
RADI step:  203 normalized residual: 4.625162e-10
RADI step:  205 normalized residual: 4.610462e-10
RADI step:  207 normalized residual: 4.510313e-10
RADI step:  209 normalized residual: 4.500268e-10
RADI step:  211 normalized residual: 4.498292e-10
RADI step:  212 normalized residual: 4.494223e-10
RADI step:  214 normalized residual: 4.472030e-10
RADI step:  216 normalized residual: 4.471251e-10
RADI step:  218 normalized residual: 4.340936e-10
RADI step:  220 normalized residual: 4.289454e-10
RADI step:  222 normalized residual: 4.083883e-10
RADI step:  223 normalized residual: 4.076139e-10
RADI step:  225 normalized residual: 3.923050e-10
RADI step:  227 normalized residual: 3.869945e-10
RADI step:  228 normalized residual: 3.715376e-10
RADI step:  229 normalized residual: 3.713599e-10
RADI step:  231 normalized residual: 3.710322e-10
RADI step:  232 normalized residual: 3.709626e-10
RADI step:  234 normalized residual: 3.615990e-10
RADI step:  236 normalized residual: 3.606716e-10
RADI step:  238 normalized residual: 3.573253e-10
RADI step:  239 normalized residual: 3.572622e-10
RADI step:  241 normalized residual: 3.555258e-10
RADI step:  242 normalized residual: 3.526018e-10
RADI step:  244 normalized residual: 3.485781e-10
RADI step:  246 normalized residual: 3.475126e-10
RADI step:  247 normalized residual: 3.474811e-10
RADI step:  249 normalized residual: 3.420404e-10
RADI step:  250 normalized residual: 3.420300e-10
RADI step:  251 normalized residual: 3.361750e-10
RADI step:  252 normalized residual: 3.237146e-10
RADI step:  254 normalized residual: 3.134949e-10
RADI step:  255 normalized residual: 3.134490e-10
RADI step:  256 normalized residual: 2.871997e-10
RADI step:  258 normalized residual: 2.736017e-10
RADI step:  259 normalized residual: 2.730073e-10
RADI step:  260 normalized residual: 2.729066e-10
RADI step:  262 normalized residual: 2.727443e-10
RADI step:  263 normalized residual: 2.389986e-10
RADI step:  264 normalized residual: 2.389982e-10
RADI step:  266 normalized residual: 2.368738e-10
RADI step:  267 normalized residual: 2.368652e-10
RADI step:  269 normalized residual: 2.202510e-10
RADI step:  270 normalized residual: 2.172891e-10
RADI step:  272 normalized residual: 2.086569e-10
RADI step:  274 normalized residual: 1.936065e-10
RADI step:  275 normalized residual: 1.935796e-10
RADI step:  276 normalized residual: 1.763168e-10
RADI step:  277 normalized residual: 9.347948e-13
Elapsed time is 13.731614 seconds.
size outradi.Z:
   606    72

compare

if not(istest)
    figure(3);
    ls_nm=[outnm.adi.niter];
    ls_radi=1:outradi.niter;

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

reset warning state

warning(orig_warnstate');

Published with MATLAB® R2018a