Contents

function ARE_rail(k, shifts, inexact, Galerkin, type, istest)
% Computes the optimal feedback via the low-rank Newton-ADI [1] and RADI
% [2] methods for the selective cooling of Steel profiles application
% described in [3,4,5].
%
% Usage: ARE_Rail(k,shifts,inexact,Galerkin,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)
%
% shifts      ADI shift selection strategy. Possible values:
%              'heur'        Penzl's heuristic shifts
%              'Wachspress'  Wachspress shifts, optimally solving the dense
%                          shift selection problem.
%             (optional, defaults to 'heur')
%
% inexact     use inexact Newton method
%             (optional, defaults to 0, i.e. false)
%
% Galerkin    activate Galerkin projection acceleration in Newton method.
%             This supersedes inexact Newton selection, i.e, disables it in
%             case both are on.
%             (optional, defaults to 0, i.e. no Galerkin acceleration)
%
% type        selector for the type of equation solved:
%               'LQR'  classic LQR ARE with Q=I, R=I and opts.LDL_T=false
%               'Hinf' robust control type ARE with Q=I,
%                      R=diag([1,1,1,1,-4,-4,-4]) and opts.LDL_T=true
%               'BR'   ARE with Q=I and R=-0.2401 I and opts.LDL_T=true as
%                      appearing in bounded real balanced truncation
%             (optional, defaults to 'LQR')
%
% 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] P. Benner, J.-R. Li, T. Penzl, Numerical solution of large-scale
%     Lyapunov equations, Riccati equations, and linear-quadratic optimal
%     control problems, Numer. Lin. Alg. Appl. 15 (9) (2008) 755–777.
%     https://doi.org/10.1002/nla.622
%
% [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, Effiziente numerische Lösung eines
%     Optimalsteuerungsproblems für die Abkühlung von Stahlprofilen,
%     Diplomarbeit, Fachbereich 3/Mathematik und Informatik, Universität
%     Bremen, D-28334 Bremen (Sep. 2003).
%     https://doi.org/10.5281/zenodo.1187040
%
% [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 Lecture Notes in Computational Science and Engineering,
%     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 (c) 2009-2025 Jens Saak, Martin Koehler, Peter Benner and others.
% All rights reserved.
% License: BSD 2-Clause License (see COPYING)
%
narginchk(0, 6);
if nargin < 1
    k = 2;
end
if nargin < 2
    shifts = 'heur';
end
if nargin < 3
    inexact = false;
end
if nargin < 4
    Galerkin = false;
end
if nargin < 5
    istest = false;
end
if nargin < 6
    type = 'LQR';
end

set operation

[oper, opts] = operatormanager(struct(), 'default');

Problem data

eqn = mess_get_linear_rail(k);
switch type
    case 'LQR'
        opts.LDL_T = false;
        eqn.R = eye(size(eqn.B, 2)); % only used in here

    case 'Hinf'
        opts.LDL_T = true;
        eqn.Q = eye(size(eqn.C, 1));
        eqn.R = diag([1, 1, 1, 1, -4, -4, -4]);

    case 'BR'
        opts.LDL_T = true;
        eqn.Q = eye(size(eqn.C, 1));
        eqn.R = -0.2401 * eye(size(eqn.B, 2));

end

if opts.LDL_T
    mytitle = '0 = C^T Q C + A^T X E + E^T X A - E^T X B R^{-1} B^T X E';
else
    mytitle = '0 = C^T C + A^T X E + E^T X A - E^T X BB^T X E';
end

First we run the Newton-ADI Method

% ADI tolerances and maximum iteration number
opts.adi.maxiter = 100;
opts.adi.res_tol = 1e-14;
opts.adi.rel_diff_tol = 1e-16;
opts.adi.info = 0;

eqn.type = 'T';

Heuristic shift parameters via basic Arnoldi

opts.shifts.num_desired = 25;
opts.shifts.num_Ritz = 50;
opts.shifts.num_hRitz = 25;
n = oper.size(eqn, opts);
opts.shifts.b0 = ones(n, 1);
switch lower(shifts)

    case 'heur'
        opts.shifts.method = 'heur';

    case 'wachspress'
        opts.shifts.method = 'wachspress';
        opts.shifts.wachspress = 'T';

    case 'projection'
        opts.shifts.method = 'projection';
end

Newton tolerances and maximum iteration number

opts.nm.maxiter = 20;
opts.nm.res_tol = 1e-12;
opts.nm.rel_diff_tol = 1e-16;
if istest
    opts.nm.info = 0;
else
    opts.nm.info = 1;
end

opts.nm.accumulateRes = true;
opts.norm = 'fro';

if Galerkin
    opts.nm.linesearch = false;
    opts.nm.inexact = false;
    opts.nm.projection.freq = 2;
    opts.nm.projection.ortho = true;
elseif inexact
    opts.nm.linesearch = true;
    opts.nm.inexact = 'quadratic';
    opts.nm.projection.freq = 0;
    opts.nm.projection.ortho = false;
else
    opts.nm.linesearch = false;
    opts.nm.inexact = false;
    opts.nm.projection.freq = 0;
    opts.nm.projection.ortho = false;
end
opts.nm.res = struct('maxiter', 10, 'tol', 1e-6, 'orth', 0);
t_mess_lrnm = tic;
outnm = mess_lrnm(eqn, opts, oper);
t_elapsed1 = toc(t_mess_lrnm);
if not(istest)
    mess_fprintf(opts, 'mess_lrnm took %6.2f seconds \n', t_elapsed1);
end

if istest
    if min(outnm.res) >= opts.nm.res_tol
        mess_err(opts, 'TEST:accuracy', 'unexpectedly inaccurate result');
    end
else
    figure(1);
    semilogy(outnm.res, 'LineWidth', 3);
    title(mytitle);
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
    mess_fprintf(opts, 'size outnm.Z: %d x %d\n\n', ...
                 size(outnm.Z, 1), size(outnm.Z, 2));
end
  
 NM step:    1 normalized residual: 	3.962591e-05  
               relative change in K: 	1.000000e+00  
               number of ADI steps: 	59   
  
  
 NM step:    2 normalized residual: 	3.111624e-08  
               relative change in K: 	2.196831e-02  
               number of ADI steps: 	54   
  
  
 NM step:    3 normalized residual: 	5.824810e-13  
               relative change in K: 	6.069980e-05  
               number of ADI steps: 	42   
  
mess_lrnm took   1.03 seconds   
size outnm.Z: 1357 x 128  
  

Lets try the RADI method and compare

% RADI-MESS settings
opts.shifts.history = opts.shifts.num_desired * size(eqn.C, 1);
opts.shifts.num_desired = opts.shifts.num_desired;

% choose either of the three shift methods, here
opts.shifts.method = 'gen-ham-opti';
% opts.shifts.method = 'heur';
% opts.shifts.method = 'projection';

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.maxiter = opts.adi.maxiter;
opts.norm = 2;
opts.radi.res_tol = opts.nm.res_tol;
opts.radi.rel_diff_tol = 0;
if istest
    opts.radi.info = 0;
else
    opts.radi.info = 1;
end

t_mess_lrradi = tic;
outradi = mess_lrradi(eqn, opts, oper);
t_elapsed2 = toc(t_mess_lrradi);

if not(istest)
    mess_fprintf(opts, 'mess_lrradi took %6.2f seconds \n', t_elapsed2);
end

if istest
    if min(outnm.res) >= opts.nm.res_tol
        mess_err(opts, 'TEST:accuracy', 'unexpectedly inaccurate result');
    end
else
    figure(2);
    semilogy(outradi.res, 'LineWidth', 3);
    title(mytitle);
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    mess_fprintf(opts, 'size outradi.Z: %d x %d\n\n', ...
                 size(outradi.Z, 1), size(outradi.Z, 2));
end
RADI step:    1 normalized residual: 3.554203e-01  
RADI step:    2 normalized residual: 1.204906e-01  
RADI step:    3 normalized residual: 1.075639e-01  
RADI step:    4 normalized residual: 3.131546e-02  
RADI step:    5 normalized residual: 9.608572e-03  
RADI step:    6 normalized residual: 7.127521e-03  
RADI step:    7 normalized residual: 1.141946e-03  
RADI step:    8 normalized residual: 6.808772e-04  
RADI step:    9 normalized residual: 5.115429e-04  
RADI step:   10 normalized residual: 5.103771e-04  
RADI step:   11 normalized residual: 5.095332e-04  
RADI step:   12 normalized residual: 1.762591e-04  
RADI step:   13 normalized residual: 1.490625e-04  
RADI step:   14 normalized residual: 2.908057e-05  
RADI step:   15 normalized residual: 1.542939e-05  
RADI step:   16 normalized residual: 4.677597e-06  
RADI step:   17 normalized residual: 4.403393e-06  
RADI step:   18 normalized residual: 2.089934e-06  
RADI step:   19 normalized residual: 7.897530e-07  
RADI step:   20 normalized residual: 6.011826e-07  
RADI step:   21 normalized residual: 5.967522e-07  
RADI step:   22 normalized residual: 1.597695e-07  
RADI step:   23 normalized residual: 1.475153e-07  
RADI step:   24 normalized residual: 3.544115e-08  
RADI step:   25 normalized residual: 3.528200e-08  
RADI step:   26 normalized residual: 1.567965e-08  
RADI step:   27 normalized residual: 1.491984e-08  
RADI step:   28 normalized residual: 2.923633e-09  
RADI step:   29 normalized residual: 1.639199e-09  
RADI step:   30 normalized residual: 3.348057e-10  
RADI step:   31 normalized residual: 3.248839e-10  
RADI step:   32 normalized residual: 2.888543e-11  
RADI step:   33 normalized residual: 2.047663e-11  
RADI step:   34 normalized residual: 1.111063e-11  
RADI step:   35 normalized residual: 1.578114e-12  
RADI step:   36 normalized residual: 1.237623e-12  
RADI step:   37 normalized residual: 4.262751e-13  
mess_lrradi took   0.97 seconds   
size outradi.Z: 1357 x 128  
  

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-', ...
             'LineWidth', 3);

    title(mytitle);
    xlabel('number of solves with A+p*M');
    ylabel('normalized residual norm');
    legend('LR-NM', 'RADI');
end

if opts.LDL_T
    [ZN, DN] = mess_column_compression(outnm.Z, 'N', outnm.D, eps, not(istest));
    ZN = ZN * diag(sqrt(diag(DN)));

    [ZR, DR] = mess_column_compression(outradi.Z, 'N', outradi.D, eps, ...
                                       not(istest));
    ZR = ZR * diag(sqrt(diag(DR)));

else
    ZN = outnm.Z;
    ZR = outradi.Z;
end

D = blkdiag(eye(size(ZN, 2)), -eye(size(ZR, 2)));

Z = [ZN ZR];
f = @(x) (Z * (D * (Z' * x)));

if k < 4
    if exist('icare', 'file')
        X = icare(full(eqn.A_), eqn.B, eqn.C' * eqn.C, eqn.R, [], full(eqn.E_), []);
    elseif exist('care', 'file')
        X = care(full(eqn.A_), eqn.B, eqn.C' * eqn.C, eqn.R, [], full(eqn.E_));
    else
        X = ZN * ZN';
    end

    relerr = norm(ZN * ZN' - ZR * ZR') / norm(ZN * ZN');
    relerr(2) = max(abs(eigs(f, size(Z, 1)))) / norm(ZN' * ZN);
    relerr(3) = norm(X - ZN * ZN') / norm(X);
    relerr(4) = norm(X - ZR * ZR') / norm(X);
else

    relerr = max(abs(eigs(f, size(Z, 1)))) / norm(ZN' * ZN);

end

if istest
    if any(relerr(1:2) > 1e-6)
        % for large examples relerr computations appear to become unstable,
        % hence the comparably large tolerance
        mess_fprintf(opts, '%s %g %g %g %g\n', type, relerr);
        mess_err(opts, 'TEST:accuracy', ...
                 'Newton and RADI solution approximations deviate too much.');

    end
else
    if k < 4
        mess_fprintf(opts, ...
                     ['relative deviation of the two solution ', ...
                      'approximations:\n%g (dense 2 norm) ', ...
                      '%g (low-rank approx 2-norm)\n'], ...
                     relerr(1), relerr(2));
        if exist('icare', 'file') || exist('care', 'file')
            mess_fprintf(opts, ...
                         ['relative deviation of the two solution ', ...
                          'approximations from (i)care''s solution:\n', ...
                          '%g (NM) %g (RADI)\n'], ...
                         relerr(3), relerr(4));
        end
    else
        mess_fprintf(opts, ['relative deviation of the two solution ', ...
                            'approximations: % g\n'], ...
                     relerr);
    end
end
relative deviation of the two solution approximations:  
8.58996e-11 (dense 2 norm) 8.58996e-11 (low-rank approx 2-norm)  
relative deviation of the two solution approximations from (i)care's solution:  
2.36705e-10 (NM) 1.59839e-10 (RADI)