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