Contents

function out =  LQR_rail_splitting(k, exp_action, method,istest)
% Computes the optimal feedback via low-rank splitting schemes [1, 2] for
% the selective cooling of Steel profiles application described in [3,4,5].
%
% Inputs:
%
% k           refinement level of the model to use
%             (1-4, i.e. 1357-79841 DOFs)
%             (optional, defaults to 1)
%
% exp_action  parameters for computing actions of matrix exponentials, see
%             mat-eqn-solvers/private/mess_exp_action.m for details
%             (optional, defaults to exp_action.method = 'Krylov',
%             exp_action.tol = 1e-8)
%
% method      choice of splitting scheme; structure with fields 'order',
%             'additive' and 'symmetric'
%             (optional, defaults to order = 2, additive = false, symmetric
%             = false)
%             NOTE: the additive schemes typically need a running parallel
%             pool in order to be competitive
%
%
% References:
%
% [1] T. Stillfjord, Low-rank second-order splitting of large-scale
%     differential Riccati equations, IEEE Trans. Autom. Control, 60
%     (2015), pp. 2791-2796. https://doi.org/10.1109/TAC.2015.2398889.
%
% [2] T. Stillfjord, Adaptive high-order splitting schemes for large-scale
%     diffederential Riccati equations, Numer. Algorithms,  (2017).
%     https://doi.org/10.1007/s11075-017-0416-8.
%
% [3] J. Saak, Effiziente numerische L\"{o}sung eines
%     Optimalsteuerungsproblems f\"{u}r die Abk\"{u}hlung von
%     Stahlprofilen, Diplomarbeit, Fachbereich 3/Mathematik und Informatik,
%     Universit\"{a}t Bremen, D-28334 Bremen (Sep. 2003).
%
% [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\"{a}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
%

if nargin < 1
    k = 1;
end
if nargin < 2
    exp_action.method = 'Krylov';
    exp_action.tol = 1e-8;
end
if nargin < 3
    method.order = 2;
    method.additive = false;
    method.symmetric = false;
end
if nargin < 4
    istest = 0;
end

Equation parameters

Default (E, A, B, C) system

oper = operatormanager('default');

eqn = getrail(k);
eqn.Rinv = 1;

eqn.type = 'T';

eqn.L0 = rand(size(eqn.A_, 1), 1);
eqn.D0 = 1;

General splitting parameters

opts.splitting.time_steps = 0 : 50 : 4500;
opts.splitting.order = method.order;
opts.splitting.additive = method.additive;
opts.splitting.symmetric = method.symmetric;
opts.splitting.info = 2;
opts.splitting.intermediates = 1;
opts.splitting.trunc_tol = 1e-10;

% Quadrature (for integral term) parameters
opts.splitting.quadrature.type = 'adaptive';
opts.splitting.quadrature.tol = 1e-4 ;

Matrix exponential action parameters

opts.exp_action = exp_action;

Compute the approximation

tic;
[out, ~, opts, ~] = mess_splitting_dre(eqn, opts, oper);
toc;
Step: 1 of 90, time = 0 
Step: 2 of 90, time = 50 
Step: 3 of 90, time = 100 
Step: 4 of 90, time = 150 
Step: 5 of 90, time = 200 
Step: 6 of 90, time = 250 
Step: 7 of 90, time = 300 
Step: 8 of 90, time = 350 
Step: 9 of 90, time = 400 
Step: 10 of 90, time = 450 
Step: 11 of 90, time = 500 
Step: 12 of 90, time = 550 
Step: 13 of 90, time = 600 
Step: 14 of 90, time = 650 
Step: 15 of 90, time = 700 
Step: 16 of 90, time = 750 
Step: 17 of 90, time = 800 
Step: 18 of 90, time = 850 
Step: 19 of 90, time = 900 
Step: 20 of 90, time = 950 
Step: 21 of 90, time = 1000 
Step: 22 of 90, time = 1050 
Step: 23 of 90, time = 1100 
Step: 24 of 90, time = 1150 
Step: 25 of 90, time = 1200 
Step: 26 of 90, time = 1250 
Step: 27 of 90, time = 1300 
Step: 28 of 90, time = 1350 
Step: 29 of 90, time = 1400 
Step: 30 of 90, time = 1450 
Step: 31 of 90, time = 1500 
Step: 32 of 90, time = 1550 
Step: 33 of 90, time = 1600 
Step: 34 of 90, time = 1650 
Step: 35 of 90, time = 1700 
Step: 36 of 90, time = 1750 
Step: 37 of 90, time = 1800 
Step: 38 of 90, time = 1850 
Step: 39 of 90, time = 1900 
Step: 40 of 90, time = 1950 
Step: 41 of 90, time = 2000 
Step: 42 of 90, time = 2050 
Step: 43 of 90, time = 2100 
Step: 44 of 90, time = 2150 
Step: 45 of 90, time = 2200 
Step: 46 of 90, time = 2250 
Step: 47 of 90, time = 2300 
Step: 48 of 90, time = 2350 
Step: 49 of 90, time = 2400 
Step: 50 of 90, time = 2450 
Step: 51 of 90, time = 2500 
Step: 52 of 90, time = 2550 
Step: 53 of 90, time = 2600 
Step: 54 of 90, time = 2650 
Step: 55 of 90, time = 2700 
Step: 56 of 90, time = 2750 
Step: 57 of 90, time = 2800 
Step: 58 of 90, time = 2850 
Step: 59 of 90, time = 2900 
Step: 60 of 90, time = 2950 
Step: 61 of 90, time = 3000 
Step: 62 of 90, time = 3050 
Step: 63 of 90, time = 3100 
Step: 64 of 90, time = 3150 
Step: 65 of 90, time = 3200 
Step: 66 of 90, time = 3250 
Step: 67 of 90, time = 3300 
Step: 68 of 90, time = 3350 
Step: 69 of 90, time = 3400 
Step: 70 of 90, time = 3450 
Step: 71 of 90, time = 3500 
Step: 72 of 90, time = 3550 
Step: 73 of 90, time = 3600 
Step: 74 of 90, time = 3650 
Step: 75 of 90, time = 3700 
Step: 76 of 90, time = 3750 
Step: 77 of 90, time = 3800 
Step: 78 of 90, time = 3850 
Step: 79 of 90, time = 3900 
Step: 80 of 90, time = 3950 
Step: 81 of 90, time = 4000 
Step: 82 of 90, time = 4050 
Step: 83 of 90, time = 4100 
Step: 84 of 90, time = 4150 
Step: 85 of 90, time = 4200 
Step: 86 of 90, time = 4250 
Step: 87 of 90, time = 4300 
Step: 88 of 90, time = 4350 
Step: 89 of 90, time = 4400 
Step: 90 of 90, time = 4450 
Elapsed time is 74.856233 seconds.
if not(istest)
    t = opts.splitting.time_steps;
    figure;
    plot(t, out.ms);
    title('Ranks of approximations over time');
end
ans = 

  struct with fields:

    Ls: {1×91 cell}
    Ds: {1×91 cell}
    ms: [91×1 double]
    Ks: {1×91 cell}