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
%             (0 - 5, i.e. 109 - 79841 Dofs)
%             (optional, defaults to 2, i.e. 1357 Dofs)
%
% 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 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)
%
narginchk(0,2);
if nargin<1, k=2; end
if nargin<2, istest=0; end

set operation

oper = operatormanager('default');

Problem data

eqn = mess_get_linear_rail(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;
t_mess_lrri = tic;
out = mess_lrri(eqn, opts, oper);
t_elapsed = toc(t_mess_lrri);
fprintf(1,'mess_lrri took %6.2f seconds \n' , t_elapsed );
RADI step:    1 pc: -2.049304e-01 + 0.000000e+00i normalized residual: 3.554203e-01 relative change in Z: 1.000000e+00
RADI step:    2 pc: -7.291802e-02 + 0.000000e+00i normalized residual: 1.204906e-01 relative change in Z: 4.677958e-01
RADI step:    3 pc: -3.795185e-01 + 0.000000e+00i normalized residual: 1.075639e-01 relative change in Z: 2.204062e-01
RADI step:    4 pc: -2.213539e-02 + 0.000000e+00i normalized residual: 3.131546e-02 relative change in Z: 3.164382e-01
RADI step:    5 pc: -7.308352e-03 + 0.000000e+00i normalized residual: 9.608572e-03 relative change in Z: 2.439606e-01
RADI step:    6 pc: -2.149665e-03 + 0.000000e+00i normalized residual: 7.127521e-03 relative change in Z: 1.750621e-01
RADI step:    7 pc: -8.255243e-01 + 0.000000e+00i normalized residual: 1.141946e-03 relative change in Z: 7.001013e-02
RADI step:    8 pc: -9.609574e-04 + 0.000000e+00i normalized residual: 6.808772e-04 relative change in Z: 7.544949e-02
RADI step:    9 pc: -9.490797e-05 + 0.000000e+00i normalized residual: 5.114959e-04 relative change in Z: 1.063467e-01
RADI step:   10 pc: -2.151612e-05 + 0.000000e+00i normalized residual: 5.103270e-04 relative change in Z: 9.510694e-02
RADI step:   11 pc: -1.096346e-05 + 0.000000e+00i normalized residual: 5.094826e-04 relative change in Z: 2.151030e-02
RADI step:   12 pc: -1.184978e+00 + 0.000000e+00i normalized residual: 1.751154e-04 relative change in Z: 1.654888e-02
RADI step:   13 pc: -3.256245e-04 + 0.000000e+00i normalized residual: 1.480689e-04 relative change in Z: 1.265424e-02
RADI step:   14 pc: -3.061260e-02 + 0.000000e+00i normalized residual: 2.841652e-05 relative change in Z: 7.993551e-03
RADI step:   15 pc: -1.656896e-03 + 0.000000e+00i normalized residual: 2.159098e-05 relative change in Z: 4.400037e-03
RADI step:   16 pc: -1.516658e+00 + 0.000000e+00i normalized residual: 1.312933e-05 relative change in Z: 2.772934e-03
RADI step:   17 pc: -1.199449e-02 + 0.000000e+00i normalized residual: 4.689772e-06 relative change in Z: 2.423689e-03
RADI step:   18 pc: -7.677708e-05 + 0.000000e+00i normalized residual: 4.599755e-06 relative change in Z: 1.252420e-03
RADI step:   19 pc: -7.048505e-02 + 0.000000e+00i normalized residual: 8.909293e-07 relative change in Z: 1.393627e-03
RADI step:   20 pc: -2.535844e-04 + 0.000000e+00i normalized residual: 8.477049e-07 relative change in Z: 6.752940e-04
RADI step:   21 pc: -1.396408e-02 + 0.000000e+00i normalized residual: 6.240996e-07 relative change in Z: 5.542046e-04
RADI step:   22 pc: -1.782274e-01 + 0.000000e+00i normalized residual: 1.629353e-07 relative change in Z: 4.298276e-04
RADI step:   23 pc: -2.428481e-05 + 0.000000e+00i normalized residual: 1.620478e-07 relative change in Z: 2.152853e-04
RADI step:   24 pc: -4.104582e-03 + 0.000000e+00i normalized residual: 1.296351e-07 relative change in Z: 2.985269e-04
RADI step:   25 pc: -5.829792e-01 + 0.000000e+00i normalized residual: 2.615142e-08 relative change in Z: 2.172291e-04
RADI step:   26 pc: -1.241860e+00 + 0.000000e+00i normalized residual: 3.866336e-09 relative change in Z: 1.080091e-04
RADI step:   27 pc: -5.439338e-04 + 0.000000e+00i normalized residual: 2.944237e-09 relative change in Z: 1.118151e-04
RADI step:   28 pc: -1.993953e+00 + 0.000000e+00i normalized residual: 2.087229e-09 relative change in Z: 3.047792e-05
RADI step:   29 pc: -1.844446e-04 + 0.000000e+00i normalized residual: 2.005066e-09 relative change in Z: 2.389045e-05
RADI step:   30 pc: -1.348064e-02 + 0.000000e+00i normalized residual: 8.371991e-10 relative change in Z: 2.531048e-05
RADI step:   31 pc: -5.362995e-02 + 0.000000e+00i normalized residual: 1.368347e-10 relative change in Z: 1.924337e-05
RADI step:   32 pc: -2.602382e-05 + 0.000000e+00i normalized residual: 1.356064e-10 relative change in Z: 1.161058e-05
RADI step:   33 pc: -8.525210e-04 + 0.000000e+00i normalized residual: 1.194191e-10 relative change in Z: 9.972628e-06
RADI step:   34 pc: -1.530985e-01 + 0.000000e+00i normalized residual: 2.368312e-11 relative change in Z: 7.127470e-06
RI step:    1 normalized residual: 4.992919e-07 relative change in Z: 1.000000e+00
               number of RADI steps:   34

RADI step:    1 pc: -3.678355e-02 + 0.000000e+00i normalized residual: 5.997529e-01 relative change in Z: 1.000000e+00
RADI step:    2 pc: -1.688497e-03 + 0.000000e+00i normalized residual: 4.306746e-01 relative change in Z: 8.867517e-01
RADI step:    3 pc: -5.372079e-05 + 0.000000e+00i normalized residual: 1.730035e-01 relative change in Z: 9.080842e-01
RADI step:    4 pc: -1.734928e-05 + 0.000000e+00i normalized residual: 1.241180e-01 relative change in Z: 4.903299e-01
RADI step:    5 pc: -1.096132e-05 + 0.000000e+00i normalized residual: 1.229144e-01 relative change in Z: 1.117713e-01
RADI step:    6 pc: -2.366190e-03 + 0.000000e+00i normalized residual: 4.718237e-02 relative change in Z: 6.153874e-02
RADI step:    7 pc: -3.256053e-04 + 0.000000e+00i normalized residual: 4.264110e-02 relative change in Z: 3.295457e-02
RADI step:    8 pc: -2.271895e-02 + 0.000000e+00i normalized residual: 1.856432e-02 relative change in Z: 2.537972e-02
RADI step:    9 pc: -2.890827e-01 + 0.000000e+00i normalized residual: 4.331742e-03 relative change in Z: 1.661998e-02
RADI step:   10 pc: -1.047223e-03 + 0.000000e+00i normalized residual: 2.792026e-03 relative change in Z: 1.201013e-02
RADI step:   11 pc: -2.089118e-05 + 0.000000e+00i normalized residual: 2.728619e-03 relative change in Z: 9.274648e-03
RADI step:   12 pc: -8.482697e-03 + 0.000000e+00i normalized residual: 1.139733e-03 relative change in Z: 8.409392e-03
RADI step:   13 pc: -2.032835e-01 + 0.000000e+00i normalized residual: 1.884239e-04 relative change in Z: 3.985238e-03
RADI step:   14 pc: -1.840586e-04 + 0.000000e+00i normalized residual: 1.814357e-04 relative change in Z: 2.273666e-03
RADI step:   15 pc: -7.937958e-01 + 0.000000e+00i normalized residual: 9.848090e-05 relative change in Z: 1.870338e-03
RADI step:   16 pc: -5.164413e-03 + 0.000000e+00i normalized residual: 4.750463e-05 relative change in Z: 1.139780e-03
RADI step:   17 pc: -1.150821e-01 + 0.000000e+00i normalized residual: 5.316456e-06 relative change in Z: 8.806282e-04
RADI step:   18 pc: -7.767812e-05 + 0.000000e+00i normalized residual: 5.067080e-06 relative change in Z: 5.861472e-04
RADI step:   19 pc: -5.437650e-04 + 0.000000e+00i normalized residual: 4.572434e-06 relative change in Z: 3.170414e-04
RADI step:   20 pc: -1.255460e+00 + 0.000000e+00i normalized residual: 2.923089e-06 relative change in Z: 2.626213e-04
RADI step:   21 pc: -4.341419e-02 + 0.000000e+00i normalized residual: 5.173011e-07 relative change in Z: 2.148451e-04
RADI step:   22 pc: -3.724736e-03 + 0.000000e+00i normalized residual: 2.181036e-07 relative change in Z: 8.481306e-05
RADI step:   23 pc: -4.983805e-01 + 0.000000e+00i normalized residual: 1.367721e-07 relative change in Z: 4.731425e-05
RADI step:   24 pc: -7.035243e-03 + 0.000000e+00i normalized residual: 3.606375e-08 relative change in Z: 4.266166e-05
RADI step:   25 pc: -2.577514e-04 + 0.000000e+00i normalized residual: 3.472662e-08 relative change in Z: 2.337819e-05
RADI step:   26 pc: -7.906304e-01 + 0.000000e+00i normalized residual: 2.607334e-08 relative change in Z: 1.677362e-05
RADI step:   27 pc: -1.379080e-02 + 0.000000e+00i normalized residual: 7.775028e-09 relative change in Z: 1.585718e-05
RADI step:   28 pc: -8.610394e-04 + 0.000000e+00i normalized residual: 7.328513e-09 relative change in Z: 9.808014e-06
RADI step:   29 pc: -2.147849e-05 + 0.000000e+00i normalized residual: 7.277039e-09 relative change in Z: 8.053074e-06
RADI step:   30 pc: -1.670089e-01 + 0.000000e+00i normalized residual: 1.331754e-09 relative change in Z: 9.901658e-06
RADI step:   31 pc: -1.841684e+00 + 0.000000e+00i normalized residual: 1.002571e-09 relative change in Z: 4.793498e-06
RADI step:   32 pc: -1.664347e-01 + 0.000000e+00i normalized residual: 1.943842e-10 relative change in Z: 3.339850e-06
RADI step:   33 pc: -1.522917e-01 + 0.000000e+00i normalized residual: 4.284752e-11 relative change in Z: 1.434721e-06
RI step:    2 normalized residual: 3.945027e-14 relative change in Z: 1.069088e-05
               number of RADI steps:   33

mess_lrri took   0.75 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,'LineWidth',3);
    hold on;
    for i = 1:length(out.radi), semilogy(out.radi(i).res,'LineWidth',3); end
    hold off;
    title(['0= C_1^T C_1 + A^T X E + E^T X A  + E^T X (\gamma^{-2}B_1 ' ...
           'B_1^T - B_2 B_2^T) X E']);
    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         127