Contents

function HINFR_rail(k, istest)
% Computes the a robust suboptimal Hinf feedback via the low-rank
% Riccati iteration [1] using the RADI method [2] for the selective cooling of
% Steel profiles application described in [3,4,5].
%
% Usage: HINFR_rail(k, istest)
%
% Inputs:
%
% k           refinement level of the model to use
%             (1-4, i.e. 1357-79841Dofs)
%             (optinal, defaults to 1)
%
% istest      flag to determine whether this demo runs as a CI test or
%             interactive demo
%             (optional, defaults to 0, i.e. interactive demo)
%
% References:
% [1] A. Lanzon, Y. Feng, B. D. O. Anderson, An iterative algorithm to
%     solve algebraic Riccati equations with an indefinite quadratic term,
%     2007 European Control Conference (ECC), pp. 3033--3039, 2007.
%     https://doi.org/10.23919/ecc.2007.7068239
%
% [2] P. Benner, Z. Bujanović, P. Kürschner, J. Saak, RADI: A low-rank
%     ADI-type algorithm for large scale algebraic Riccati equations,
%     Numer. Math. 138 (2) (2018) 301–330.
%     https://doi.org/10.1007/s00211-017-0907-5.
%
% [3] J. Saak, Efficient numerical solution of large scale algebraic matrix
%     equations in PDE control and model order reduction, Dissertation,
%     Technische Universität Chemnitz, Chemnitz, Germany (Jul. 2009).
%     URL http://nbn-resolving.de/urn:nbn:de:bsz:ch1-200901642
%
% [4] P. Benner, J. Saak, A semi-discretized heat transfer model for
%     optimal cooling of steel profiles, in: P. Benner, V. Mehrmann, D.
%     Sorensen (Eds.), Dimension Reduction of Large-Scale Systems, Vol. 45
%     of Lect. Notes Comput. Sci. Eng., Springer-Verlag, Berlin/Heidelberg,
%     Germany, 2005, pp. 353–356. https://doi.org/10.1007/3-540-27909-1_19.
%
% [5] J. Saak, Efficient numerical solution of large scale algebraic matrix
%     equations in PDE control and model order reduction, Dissertation,
%     Technische Universität Chemnitz, Chemnitz, Germany (Jul. 2009).
%     URL http://nbn-resolving.de/urn:nbn:de:bsz:ch1-200901642
%

%
% 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
%
narginchk(0,2);
if nargin<1, k=1; end
if nargin<2, istest=0; end

set operation

oper = operatormanager('default');

Problem data

eqn = getrail(k);
% Reformulate as normalized Hinf control problem.
eqn.B1 = eqn.B;
eqn.B2 = eqn.B;
eqn.C1 = eqn.C;
eqn = rmfield(eqn,'B');
eqn = rmfield(eqn,'C');

Optional parameters.

opts.norm                = 2;
% Shift selection settings.
opts.shifts.history = 42;
opts.shifts.num_desired      = 5;
% choose either of the three shift methods, here
opts.shifts.method = 'gen-ham-opti';
%opts.shifts.method = 'heur';
%opts.shifts.method = 'projection';

% RADI settings
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.maxiter             = 200;
opts.radi.res_tol             = 1.0e-10;
opts.radi.rel_diff_tol        = 1.0e-16;
opts.radi.info                = 1;

% Riccati iteration settings.
opts.ri.riccati_solver = 'radi';
opts.ri.maxiter        = 10;
opts.ri.res_tol        = 1.0e-10;
opts.ri.rel_diff_tol   = 1.0e-16;
opts.ri.compres_tol    = 1.0e-16;
opts.ri.info           = 1;

Solve the equation.

eqn.type = 'T';
gam      = 10;
eqn.B1   = 1/gam * eqn.B1;
tic;
out = mess_lrri(eqn, opts, oper);
toc;
RADI step:    1 pc: -1.606663e-01 + 0.000000e+00i normalized residual: 2.938578e-01 relative change in Z: 1.000000e+00
RADI step:    2 pc: -5.869096e-02 + 0.000000e+00i normalized residual: 1.055053e-01 relative change in Z: 4.923292e-01
RADI step:    3 pc: -3.108071e-01 + 0.000000e+00i normalized residual: 9.156909e-02 relative change in Z: 2.239340e-01
RADI step:    4 pc: -1.851154e-02 + 0.000000e+00i normalized residual: 2.423910e-02 relative change in Z: 3.289936e-01
RADI step:    5 pc: -6.392326e-03 + 0.000000e+00i normalized residual: 1.289624e-02 relative change in Z: 2.445329e-01
RADI step:    6 pc: -1.912665e-03 + 0.000000e+00i normalized residual: 1.109388e-02 relative change in Z: 1.720236e-01
RADI step:    7 pc: -7.918107e-01 + 0.000000e+00i normalized residual: 1.318018e-03 relative change in Z: 8.532289e-02
RADI step:    8 pc: -9.356477e-04 + 0.000000e+00i normalized residual: 1.058302e-03 relative change in Z: 7.176042e-02
RADI step:    9 pc: -9.411802e-05 + 0.000000e+00i normalized residual: 1.050650e-03 relative change in Z: 1.058962e-01
RADI step:   10 pc: -3.344197e-05 + 0.000000e+00i normalized residual: 1.047080e-03 relative change in Z: 6.538866e-02
RADI step:   11 pc: -1.217272e+00 + 0.000000e+00i normalized residual: 1.300052e-04 relative change in Z: 2.417587e-02
RADI step:   12 pc: -1.807529e-05 + 0.000000e+00i normalized residual: 1.267997e-04 relative change in Z: 1.252595e-02
RADI step:   13 pc: -3.377699e-04 + 0.000000e+00i normalized residual: 1.140249e-04 relative change in Z: 1.099804e-02
RADI step:   14 pc: -3.841773e-02 + 0.000000e+00i normalized residual: 5.691451e-05 relative change in Z: 7.682463e-03
RADI step:   15 pc: -1.715145e+00 + 0.000000e+00i normalized residual: 2.032024e-05 relative change in Z: 5.169212e-03
RADI step:   16 pc: -3.916306e-03 + 0.000000e+00i normalized residual: 6.847217e-06 relative change in Z: 4.287567e-03
RADI step:   17 pc: -1.539283e-02 + 0.000000e+00i normalized residual: 4.471060e-06 relative change in Z: 1.740819e-03
RADI step:   18 pc: -2.656241e-04 + 0.000000e+00i normalized residual: 4.406977e-06 relative change in Z: 1.265239e-03
RADI step:   19 pc: -1.326839e-01 + 0.000000e+00i normalized residual: 1.100128e-06 relative change in Z: 1.239047e-03
RADI step:   20 pc: -2.033744e+00 + 0.000000e+00i normalized residual: 3.104845e-07 relative change in Z: 6.054749e-04
RADI step:   21 pc: -1.110897e-03 + 0.000000e+00i normalized residual: 2.770622e-07 relative change in Z: 4.534670e-04
RADI step:   22 pc: -4.827316e-05 + 0.000000e+00i normalized residual: 2.765238e-07 relative change in Z: 3.819660e-04
RADI step:   23 pc: -2.162945e-01 + 0.000000e+00i normalized residual: 5.297653e-08 relative change in Z: 3.058878e-04
RADI step:   24 pc: -8.579284e-03 + 0.000000e+00i normalized residual: 3.615335e-08 relative change in Z: 1.738656e-04
RADI step:   25 pc: -4.307844e-01 + 0.000000e+00i normalized residual: 4.633696e-09 relative change in Z: 1.101018e-04
RADI step:   26 pc: -3.194418e-05 + 0.000000e+00i normalized residual: 4.428805e-09 relative change in Z: 9.026628e-05
RADI step:   27 pc: -5.593559e-04 + 0.000000e+00i normalized residual: 3.991738e-09 relative change in Z: 6.418440e-05
RADI step:   28 pc: -3.742716e-02 + 0.000000e+00i normalized residual: 9.233411e-10 relative change in Z: 4.370642e-05
RADI step:   29 pc: -1.929914e-04 + 0.000000e+00i normalized residual: 8.959335e-10 relative change in Z: 1.677796e-05
RADI step:   30 pc: -2.944286e-01 + 0.000000e+00i normalized residual: 3.458257e-10 relative change in Z: 1.675726e-05
RADI step:   31 pc: -1.614106e-03 + 0.000000e+00i normalized residual: 3.011897e-10 relative change in Z: 1.443999e-05
RADI step:   32 pc: -9.074990e-01 + 0.000000e+00i normalized residual: 1.117139e-10 relative change in Z: 1.090478e-05
RADI step:   33 pc: -1.498212e-02 + 0.000000e+00i normalized residual: 5.956410e-11 relative change in Z: 7.049216e-06
RI step:    1  normalized residual: 6.978943e-07 relative change in Z: 1.000000e+00
               number of RADI steps:   33

RADI step:    1 pc: -5.603512e-02 + 0.000000e+00i normalized residual: 5.443098e-01 relative change in Z: 1.000000e+00
RADI step:    2 pc: -3.299848e-03 + 0.000000e+00i normalized residual: 2.131345e-01 relative change in Z: 8.701085e-01
RADI step:    3 pc: -2.767165e-04 + 0.000000e+00i normalized residual: 1.398984e-01 relative change in Z: 6.953813e-01
RADI step:    4 pc: -3.110817e-05 + 0.000000e+00i normalized residual: 7.857360e-02 relative change in Z: 6.622756e-01
RADI step:    5 pc: -1.001218e-03 + 0.000000e+00i normalized residual: 5.059555e-02 relative change in Z: 1.251979e-01
RADI step:    6 pc: -1.811640e-05 + 0.000000e+00i normalized residual: 5.015026e-02 relative change in Z: 8.921797e-02
RADI step:    7 pc: -8.071368e-03 + 0.000000e+00i normalized residual: 1.377057e-02 relative change in Z: 5.459492e-02
RADI step:    8 pc: -9.163518e-02 + 0.000000e+00i normalized residual: 3.651743e-03 relative change in Z: 2.159221e-02
RADI step:    9 pc: -8.833609e-05 + 0.000000e+00i normalized residual: 3.569510e-03 relative change in Z: 1.281547e-02
RADI step:   10 pc: -3.140142e-01 + 0.000000e+00i normalized residual: 9.984065e-04 relative change in Z: 1.081364e-02
RADI step:   11 pc: -1.386174e-03 + 0.000000e+00i normalized residual: 6.363396e-04 relative change in Z: 8.847663e-03
RADI step:   12 pc: -1.958887e-02 + 0.000000e+00i normalized residual: 1.309624e-04 relative change in Z: 6.406526e-03
RADI step:   13 pc: -4.655198e-01 + 0.000000e+00i normalized residual: 3.404868e-05 relative change in Z: 2.412996e-03
RADI step:   14 pc: -5.445532e-04 + 0.000000e+00i normalized residual: 2.685013e-05 relative change in Z: 2.070759e-03
RADI step:   15 pc: -1.300257e-02 + 0.000000e+00i normalized residual: 1.329958e-05 relative change in Z: 1.164557e-03
RADI step:   16 pc: -1.518936e+00 + 0.000000e+00i normalized residual: 4.338422e-06 relative change in Z: 8.401227e-04
RADI step:   17 pc: -2.263412e-04 + 0.000000e+00i normalized residual: 3.928614e-06 relative change in Z: 5.743156e-04
RADI step:   18 pc: -6.291902e-03 + 0.000000e+00i normalized residual: 2.445629e-06 relative change in Z: 3.836634e-04
RADI step:   19 pc: -1.646897e-01 + 0.000000e+00i normalized residual: 4.038141e-07 relative change in Z: 2.850626e-04
RADI step:   20 pc: -3.078332e-05 + 0.000000e+00i normalized residual: 3.991851e-07 relative change in Z: 1.682835e-04
RADI step:   21 pc: -2.107837e-03 + 0.000000e+00i normalized residual: 3.086411e-07 relative change in Z: 1.305534e-04
RADI step:   22 pc: -1.095082e+00 + 0.000000e+00i normalized residual: 1.500644e-07 relative change in Z: 9.466912e-05
RADI step:   23 pc: -3.424694e-02 + 0.000000e+00i normalized residual: 1.296167e-08 relative change in Z: 8.206775e-05
RADI step:   24 pc: -3.835096e-04 + 0.000000e+00i normalized residual: 1.139154e-08 relative change in Z: 4.131827e-05
RADI step:   25 pc: -8.576817e-01 + 0.000000e+00i normalized residual: 4.157178e-09 relative change in Z: 2.012566e-05
RADI step:   26 pc: -1.824444e-05 + 0.000000e+00i normalized residual: 4.127396e-09 relative change in Z: 1.328583e-05
RADI step:   27 pc: -2.451899e-03 + 0.000000e+00i normalized residual: 3.305641e-09 relative change in Z: 1.310778e-05
RADI step:   28 pc: -5.721962e-01 + 0.000000e+00i normalized residual: 1.479213e-09 relative change in Z: 9.110129e-06
RADI step:   29 pc: -2.710512e-02 + 0.000000e+00i normalized residual: 5.460936e-10 relative change in Z: 7.558356e-06
RADI step:   30 pc: -3.394209e-04 + 0.000000e+00i normalized residual: 5.027410e-10 relative change in Z: 4.879596e-06
RADI step:   31 pc: -5.056266e-05 + 0.000000e+00i normalized residual: 5.005798e-10 relative change in Z: 1.363240e-06
RADI step:   32 pc: -2.041086e+00 + 0.000000e+00i normalized residual: 2.523019e-10 relative change in Z: 4.065786e-06
RADI step:   33 pc: -1.628487e-01 + 0.000000e+00i normalized residual: 4.014460e-11 relative change in Z: 2.936175e-06
RI step:    2  normalized residual: 8.421675e-15 relative change in Z: 5.820222e-06
               number of RADI steps:   33

Elapsed time is 0.633691 seconds.

Residual behavior.

if istest
    if min(out.res) >= opts.ri.res_tol
       error('MESS:TEST:accuracy','unexpectedly inaccurate result');
   end
else
    figure(1);
    semilogy(out.res);
    hold on;
    for i = 1:length(out.radi), semilogy(out.radi(i).res); end
    hold off;
    title('0= C_1^TC_1 + A^TXM + M^TXA  + M^TX(\gamma^{-2}B_1B_1^T - B_2B_2^T)XM');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    legend('Riccati Iteration', 'RADI (step 1)', 'RADI (step 2)');
end
disp('size out.Z:');
disp(size(out.Z));
size out.Z:
        1357         130