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}
