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