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 Transactions % on Power Systems 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 (c) 2009-2025 Jens Saak, Martin Koehler, Peter Benner and others. % All rights reserved. % License: BSD 2-Clause License (see COPYING) %
if nargin < 1 istest = false; end
set operation
opts = struct();
[oper, opts] = operatormanager(opts, '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 = true; opts.nm.accumulateRes = true;
t_mess_lrnm = tic; outnm = mess_lrnm(eqn, opts, oper); t_elapsed1 = toc(t_mess_lrnm); mess_fprintf(opts, 'mess_lrnm took %6.2f seconds \n', t_elapsed1); if istest if min(outnm.res) >= opts.nm.res_tol mess_err(opts, '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 [mZ, nZ] = size(outnm.Z); mess_fprintf(opts, 'size outnm.Z: %d x %d\n', mZ, nZ);
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: 188 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: 236 NM step: 3 normalized residual: 8.550157e-03 relative change in K: 9.423954e-01 number of ADI steps: 238 NM step: 4 normalized residual: 2.080652e-03 relative change in K: 7.739811e-01 number of ADI steps: 219 NM step: 5 normalized residual: 5.951448e-04 relative change in K: 4.726820e-01 number of ADI steps: 208 NM step: 6 normalized residual: 5.145986e-04 relative change in K: 3.286749e-01 number of ADI steps: 198 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: 224 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: 214 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: 198 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: 167 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: 199 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: 175 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: 152 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: 187 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: 166 NM step: 16 normalized residual: 1.410808e-04 relative change in K: 8.932130e-02 number of ADI steps: 180 NM step: 17 normalized residual: 5.011963e-06 relative change in K: 1.902161e-02 number of ADI steps: 148 NM step: 18 normalized residual: 4.131444e-07 relative change in K: 5.234902e-03 number of ADI steps: 133 NM step: 19 normalized residual: 7.625077e-08 relative change in K: 2.143992e-03 number of ADI steps: 127 NM step: 20 normalized residual: 5.861314e-09 relative change in K: 6.077819e-04 number of ADI steps: 122 NM step: 21 normalized residual: 6.364230e-11 relative change in K: 6.659441e-05 number of ADI steps: 85 mess_lrnm took 137.79 seconds size outnm.Z: 606 x 119

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 = true; opts.radi.get_ZZt = true; opts.radi.compute_res = false; 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); mess_fprintf(opts, 'mess_lrradi took %6.2f seconds \n', t_elapsed2); if istest if min(outradi.res) >= opts.radi.res_tol mess_err(opts, '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 [mZ, nZ] = size(outradi.Z); mess_fprintf(opts, 'size outradi.Z: %d x %d\n', mZ, nZ);
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.000264e-05 RADI step: 25 normalized residual: 2.046697e-05 RADI step: 26 normalized residual: 2.078785e-05 RADI step: 27 normalized residual: 2.266193e-05 RADI step: 28 normalized residual: 2.410330e-05 RADI step: 30 normalized residual: 2.385302e-05 RADI step: 32 normalized residual: 1.333946e-07 RADI step: 33 normalized residual: 1.329083e-07 RADI step: 35 normalized residual: 1.573559e-07 RADI step: 37 normalized residual: 1.590883e-07 RADI step: 38 normalized residual: 1.601287e-07 RADI step: 40 normalized residual: 1.626166e-07 RADI step: 41 normalized residual: 1.270577e-07 RADI step: 43 normalized residual: 1.273634e-07 RADI step: 44 normalized residual: 1.143290e-07 RADI step: 46 normalized residual: 1.153258e-07 RADI step: 48 normalized residual: 1.028288e-07 RADI step: 50 normalized residual: 9.266728e-08 RADI step: 51 normalized residual: 9.228526e-08 RADI step: 53 normalized residual: 9.576228e-08 RADI step: 54 normalized residual: 8.870237e-08 RADI step: 56 normalized residual: 8.791572e-08 RADI step: 58 normalized residual: 8.459842e-08 RADI step: 59 normalized residual: 7.546439e-09 RADI step: 60 normalized residual: 1.875140e-08 RADI step: 61 normalized residual: 1.892236e-08 RADI step: 63 normalized residual: 1.932520e-08 RADI step: 65 normalized residual: 1.947682e-08 RADI step: 67 normalized residual: 1.980736e-08 RADI step: 68 normalized residual: 2.037698e-08 RADI step: 70 normalized residual: 1.961033e-08 RADI step: 71 normalized residual: 2.015205e-08 RADI step: 73 normalized residual: 2.062464e-08 RADI step: 75 normalized residual: 2.058992e-08 RADI step: 77 normalized residual: 2.040077e-08 RADI step: 79 normalized residual: 2.066021e-08 RADI step: 81 normalized residual: 2.063680e-08 RADI step: 83 normalized residual: 2.063209e-08 RADI step: 84 normalized residual: 2.050078e-08 RADI step: 86 normalized residual: 5.172236e-09 RADI step: 87 normalized residual: 5.146273e-09 RADI step: 89 normalized residual: 4.552584e-09 RADI step: 91 normalized residual: 4.407804e-09 RADI step: 93 normalized residual: 4.175370e-09 RADI step: 94 normalized residual: 4.021581e-09 RADI step: 96 normalized residual: 4.078232e-09 RADI step: 98 normalized residual: 4.075851e-09 RADI step: 100 normalized residual: 4.028576e-09 RADI step: 102 normalized residual: 3.618141e-09 RADI step: 104 normalized residual: 3.512252e-09 RADI step: 106 normalized residual: 3.500713e-09 RADI step: 107 normalized residual: 3.473549e-09 RADI step: 109 normalized residual: 3.413547e-09 RADI step: 111 normalized residual: 3.340438e-09 RADI step: 112 normalized residual: 3.315342e-09 RADI step: 114 normalized residual: 3.081033e-09 RADI step: 116 normalized residual: 3.037409e-09 RADI step: 118 normalized residual: 3.001774e-09 RADI step: 120 normalized residual: 2.990168e-09 RADI step: 121 normalized residual: 1.142086e-10 RADI step: 123 normalized residual: 1.137821e-10 RADI step: 124 normalized residual: 1.118964e-10 RADI step: 125 normalized residual: 1.153001e-10 RADI step: 127 normalized residual: 1.176298e-10 RADI step: 129 normalized residual: 1.182581e-10 RADI step: 130 normalized residual: 1.678357e-10 RADI step: 132 normalized residual: 1.967571e-10 RADI step: 134 normalized residual: 2.053741e-10 RADI step: 136 normalized residual: 2.017562e-10 RADI step: 137 normalized residual: 2.019593e-10 RADI step: 138 normalized residual: 2.023465e-10 RADI step: 140 normalized residual: 2.033768e-10 RADI step: 142 normalized residual: 2.034074e-10 RADI step: 143 normalized residual: 2.034336e-10 RADI step: 145 normalized residual: 2.034683e-10 RADI step: 146 normalized residual: 2.046968e-10 RADI step: 147 normalized residual: 2.059176e-10 RADI step: 148 normalized residual: 2.063508e-10 RADI step: 150 normalized residual: 2.073416e-10 RADI step: 152 normalized residual: 2.079803e-10 RADI step: 154 normalized residual: 2.082592e-10 RADI step: 155 normalized residual: 2.082119e-10 RADI step: 157 normalized residual: 2.077113e-10 RADI step: 158 normalized residual: 2.039846e-10 RADI step: 160 normalized residual: 2.033916e-10 RADI step: 162 normalized residual: 1.987409e-10 RADI step: 164 normalized residual: 2.001868e-10 RADI step: 165 normalized residual: 1.938879e-10 RADI step: 166 normalized residual: 1.902885e-10 RADI step: 168 normalized residual: 1.854150e-10 RADI step: 169 normalized residual: 1.830819e-10 RADI step: 171 normalized residual: 1.797678e-10 RADI step: 172 normalized residual: 1.748152e-10 RADI step: 174 normalized residual: 1.734806e-10 RADI step: 176 normalized residual: 1.680378e-10 RADI step: 178 normalized residual: 1.618832e-10 RADI step: 180 normalized residual: 1.617886e-10 RADI step: 182 normalized residual: 1.612449e-10 RADI step: 184 normalized residual: 1.606989e-10 RADI step: 185 normalized residual: 1.574595e-10 RADI step: 186 normalized residual: 1.574352e-10 RADI step: 188 normalized residual: 1.567034e-10 RADI step: 190 normalized residual: 1.558677e-10 RADI step: 191 normalized residual: 1.538670e-10 RADI step: 193 normalized residual: 1.463787e-10 RADI step: 195 normalized residual: 1.299298e-10 RADI step: 197 normalized residual: 1.260510e-10 RADI step: 199 normalized residual: 1.243959e-10 RADI step: 200 normalized residual: 1.157281e-10 RADI step: 202 normalized residual: 1.082525e-10 RADI step: 203 normalized residual: 9.735267e-11 mess_lrradi took 10.90 seconds size outradi.Z: 606 x 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');