function LQR_rail_BDF(k)
% Computes the optimal feedback via the low-rank Rosenbrock[1,2] methods for
% the selective cooling of Steel profiles application described in [3,4,5].
%
% Usage: LQR_Rail(k,shifts,inexact,Galerkin,istest)
%
% Inputs:
%
% k           k-step BDF method
%             possible values: 1, ..., 6
%             (optinal, defaults to 2)
%
% References:
%
% [1] N. Lang, H. Mena, J. Saak, On the benefits of the LDLT factorization
%     for largescale differential matrix equation solvers, Linear Algebra
%     Appl. 480 (2015) 4471.  https://doi.org/10.1016/j.laa.2015.04.006.
%
% [2] N. Lang, Numerical methods for large-scale linear time-varying
%     control systems and related differential matrix equations,
%     Dissertation, Technische Universitt Chemnitz, Chemnitz, Germany,
%     logos-Verlag, Berlin, ISBN 978-3-8325-4700-4 (Jun. 2017).
%     URL https://www.logos-verlag.de/cgi-bin/buch/isbn/4700
%
% [3] J. Saak, Effiziente numerische Lsung eines
%     Optimalsteuerungsproblems fr die Abkhlung von Stahlprofilen,
%     Diplomarbeit, Fachbereich 3/Mathematik und Informatik, Universitt
%     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. 353356. 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 Universitt 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,1);
if nargin<1, k =2; end

set operation

oper = operatormanager('default');
% Problem data

eqn=getrail(0);
opts.norm = 'fro';

% 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;
opts.adi.compute_sol_fac = 1;
opts.cc_info = 0;

eqn.type = 'T';
%Heuristic shift parameters via basic Arnoldi
n=oper.size(eqn, opts);
opts.shifts.num_desired=7;
opts.shifts.num_Ritz=50;
opts.shifts.num_hRitz=25;
opts.shifts.method = 'heur';

opts.shifts.b0=ones(n,1);

Newton tolerances and maximum iteration number

opts.nm.maxiter = 8;
opts.nm.res_tol = 1e-10;
opts.nm.rel_diff_tol = 1e-16;
opts.nm.info = 0;
opts.norm = 'fro';
opts.nm.accumulateRes = 1;
opts.nm.linesearch = 1;

BDF parameters

opts.bdf.time_steps = 0 : 50: 4500;
opts.bdf.step = k;
opts.bdf.info = 1;
opts.bdf.save_solution = 0;
opts.bdf.startup_iter = 7;
tic;
[out_bdf]=mess_bdf_dre(eqn,opts,oper);
toc;
Warning: Initial condition factor L0 is not defined or corrupted. Setting it to
the zero vector. 
Warning: Initial condition factor D0 is not defined or corrupted. Setting it to
the identity matrix. 
BDF step 4450 s
	 Newton steps:  2 
	 Rank 1
BDF step 4400 s
	 Newton steps:  2 
	 Rank 79
BDF step 4350 s
	 Newton steps:  2 
	 Rank 79
BDF step 4300 s
	 Newton steps:  2 
	 Rank 82
BDF step 4250 s
	 Newton steps:  2 
	 Rank 84
BDF step 4200 s
	 Newton steps:  2 
	 Rank 85
BDF step 4150 s
	 Newton steps:  2 
	 Rank 86
BDF step 4100 s
	 Newton steps:  2 
	 Rank 88
BDF step 4050 s
	 Newton steps:  2 
	 Rank 89
BDF step 4000 s
	 Newton steps:  2 
	 Rank 90
BDF step 3950 s
	 Newton steps:  2 
	 Rank 92
BDF step 3900 s
	 Newton steps:  2 
	 Rank 92
BDF step 3850 s
	 Newton steps:  2 
	 Rank 92
BDF step 3800 s
	 Newton steps:  2 
	 Rank 94
BDF step 3750 s
	 Newton steps:  2 
	 Rank 94
BDF step 3700 s
	 Newton steps:  2 
	 Rank 95
BDF step 3650 s
	 Newton steps:  2 
	 Rank 97
BDF step 3600 s
	 Newton steps:  2 
	 Rank 97
BDF step 3550 s
	 Newton steps:  2 
	 Rank 97
BDF step 3500 s
	 Newton steps:  2 
	 Rank 98
BDF step 3450 s
	 Newton steps:  2 
	 Rank 98
BDF step 3400 s
	 Newton steps:  2 
	 Rank 98
BDF step 3350 s
	 Newton steps:  2 
	 Rank 99
BDF step 3300 s
	 Newton steps:  2 
	 Rank 100
BDF step 3250 s
	 Newton steps:  2 
	 Rank 100
BDF step 3200 s
	 Newton steps:  2 
	 Rank 101
BDF step 3150 s
	 Newton steps:  2 
	 Rank 101
BDF step 3100 s
	 Newton steps:  2 
	 Rank 102
BDF step 3050 s
	 Newton steps:  2 
	 Rank 102
BDF step 3000 s
	 Newton steps:  2 
	 Rank 102
BDF step 2950 s
	 Newton steps:  2 
	 Rank 102
BDF step 2900 s
	 Newton steps:  2 
	 Rank 102
BDF step 2850 s
	 Newton steps:  2 
	 Rank 103
BDF step 2800 s
	 Newton steps:  2 
	 Rank 103
BDF step 2750 s
	 Newton steps:  2 
	 Rank 103
BDF step 2700 s
	 Newton steps:  2 
	 Rank 104
BDF step 2650 s
	 Newton steps:  2 
	 Rank 104
BDF step 2600 s
	 Newton steps:  2 
	 Rank 104
BDF step 2550 s
	 Newton steps:  2 
	 Rank 105
BDF step 2500 s
	 Newton steps:  2 
	 Rank 105
BDF step 2450 s
	 Newton steps:  2 
	 Rank 105
BDF step 2400 s
	 Newton steps:  2 
	 Rank 105
BDF step 2350 s
	 Newton steps:  2 
	 Rank 105
BDF step 2300 s
	 Newton steps:  2 
	 Rank 105
BDF step 2250 s
	 Newton steps:  2 
	 Rank 106
BDF step 2200 s
	 Newton steps:  2 
	 Rank 106
BDF step 2150 s
	 Newton steps:  2 
	 Rank 106
BDF step 2100 s
	 Newton steps:  2 
	 Rank 106
BDF step 2050 s
	 Newton steps:  2 
	 Rank 106
BDF step 2000 s
	 Newton steps:  2 
	 Rank 106
BDF step 1950 s
	 Newton steps:  2 
	 Rank 107
BDF step 1900 s
	 Newton steps:  2 
	 Rank 107
BDF step 1850 s
	 Newton steps:  2 
	 Rank 107
BDF step 1800 s
	 Newton steps:  2 
	 Rank 107
BDF step 1750 s
	 Newton steps:  2 
	 Rank 107
BDF step 1700 s
	 Newton steps:  2 
	 Rank 107
BDF step 1650 s
	 Newton steps:  2 
	 Rank 108
BDF step 1600 s
	 Newton steps:  2 
	 Rank 108
BDF step 1550 s
	 Newton steps:  2 
	 Rank 108
BDF step 1500 s
	 Newton steps:  2 
	 Rank 108
BDF step 1450 s
	 Newton steps:  2 
	 Rank 108
BDF step 1400 s
	 Newton steps:  2 
	 Rank 108
BDF step 1350 s
	 Newton steps:  2 
	 Rank 108
BDF step 1300 s
	 Newton steps:  2 
	 Rank 108
BDF step 1250 s
	 Newton steps:  2 
	 Rank 108
BDF step 1200 s
	 Newton steps:  2 
	 Rank 108
BDF step 1150 s
	 Newton steps:  2 
	 Rank 108
BDF step 1100 s
	 Newton steps:  2 
	 Rank 108
BDF step 1050 s
	 Newton steps:  2 
	 Rank 108
BDF step 1000 s
	 Newton steps:  2 
	 Rank 109
BDF step  950 s
	 Newton steps:  2 
	 Rank 109
BDF step  900 s
	 Newton steps:  2 
	 Rank 109
BDF step  850 s
	 Newton steps:  2 
	 Rank 109
BDF step  800 s
	 Newton steps:  2 
	 Rank 109
BDF step  750 s
	 Newton steps:  2 
	 Rank 109
BDF step  700 s
	 Newton steps:  2 
	 Rank 109
BDF step  650 s
	 Newton steps:  2 
	 Rank 109
BDF step  600 s
	 Newton steps:  2 
	 Rank 109
BDF step  550 s
	 Newton steps:  2 
	 Rank 109
BDF step  500 s
	 Newton steps:  2 
	 Rank 109
BDF step  450 s
	 Newton steps:  2 
	 Rank 109
BDF step  400 s
	 Newton steps:  2 
	 Rank 109
BDF step  350 s
	 Newton steps:  2 
	 Rank 109
BDF step  300 s
	 Newton steps:  2 
	 Rank 110
BDF step  250 s
	 Newton steps:  2 
	 Rank 110
BDF step  200 s
	 Newton steps:  2 
	 Rank 110
BDF step  150 s
	 Newton steps:  2 
	 Rank 110
BDF step  100 s
	 Newton steps:  2 
	 Rank 110
BDF step   50 s
	 Newton steps:  2 
	 Rank 110
BDF step    0 s
	 Newton steps:  2 
	 Rank 110
Elapsed time is 36.836962 seconds.
y = zeros(1,length(out_bdf.Ks));
for i=1:length(out_bdf.Ks)
    y(i) = out_bdf.Ks{i}(1,77);
end
x = opts.bdf.time_steps;
figure(1);
plot(x,y);
title('evolution of component (1,77) of the optimal feedback');