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