Contents

function out = LQR_LTV_smallscale_splitting(method, istest)
% Computes the optimal feedback via low-rank splitting schemes [1, 2] for
% a small-scale LTV system, where the matrices are time-varying in the
% sense that A(t) = a(t) \hat{A}.
%
% Inputs:
%
% 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
%
% 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] 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
%     differential Riccati equations, Numer. Algorithms,  (2017).
%     https://doi.org/10.1007/s11075-017-0416-8.

%
% 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
    method.order = 2;
    method.additive = false;
    method.symmetric = false;
end
if nargin < 2
    istest = 0;
end

Nx = 10;
Nx2 = Nx^2;

% Finite-difference approximation of Laplacian on [0,1]^2 with Dirichlet
% boundary conditions
xx = linspace(0,1, Nx+2);
x = xx(2:end-1);
dx = 1/(Nx+1);
e = ones(Nx, 1);
A_1D = 1/dx^2*spdiags([e, -2*e, e], -1:1, Nx, Nx);

[X, Y] = meshgrid(x,x);
I_1D = eye(Nx, Nx);
A_2D = kron(A_1D, I_1D) + kron(I_1D, A_1D);
M = speye(Nx2, Nx2); % mass matrix = I in this case

% 3 inputs, control on the squares [j/4, j/4 + 1/8] x [j/4, j/4 + 1/8]
Nb = 3;
B = zeros(Nx2, Nb);
for j = 1:Nb
    B(:, j) = reshape((X > j/4) & (X < j/4 + 1/8) & (Y > j/4) & (Y < j/4 + 1/8), Nx2, 1);
end

% 1 output, average of all states
C = 1/Nx2 * ones(1, Nx2);

% Time dependency through factors
alpha = @(t) 1 + 10*sin(2*pi*t); % A
mu = @(t) 2 + 7.1*sin(2*pi*t);   % M
dmu = @(t) 7.1*2*pi*cos(2*pi*t); % dM/dt
beta = @(t) 3 + cos(t);          % B
gamma = @(t) 1 - min(t, 1);      % C

% Alternative, autonomous case:
% alpha = @(t) 1;        % A
% mu = @(t) 1;           % M
% dmu = @(t) 0;          % dM/dt
% beta = @(t) 1;         % B
% gamma = @(t) 1;        % C

At = @(t) alpha(t) * A_2D;
Mt = @(t) mu(t) * M;
dMt = @(t) dmu(t) * M;
Bt = @(t) beta(t) * B;
Ct = @(t) gamma(t) * C;

eqn.A_time = At;
eqn.E_time = Mt;
eqn.dt_E_time = dMt;
eqn.B_time = Bt;
eqn.C_time = Ct;

eqn.haveE = 1;

eqn.type = 'T';

% LDLT-factorization of initial condition P(0) = 0
L0 = zeros(Nx2, 1);
eqn.L0 = L0;
eqn.D0 = eye(size(L0, 2));

% R^{-1}, inverse of weighting factor for input in cost functional
eqn.Rinv = 1;

oper = operatormanager('default');

eqn.LTV = 1; % Specify that this is a time-varying problem

% Time interval [0, 0.1] and 100 time steps
t0 = 0;
tend = 0.1;
Nt = 100;

General splitting parameters

opts.splitting.time_steps = linspace(t0, tend, Nt+1);
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 = eps;

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

Matrix exponential actions

opts.exp_action.method = 'LTV';
opts.exp_action.tol = 1e-8;

Compute the approximation

tic;
[out, ~,opts, ~] = mess_splitting_dre(eqn,opts,oper);
toc;
Step: 1 of 100, time = 0 
Step: 2 of 100, time = 1.000000e-03 
Step: 3 of 100, time = 2.000000e-03 
Step: 4 of 100, time = 3.000000e-03 
Step: 5 of 100, time = 4.000000e-03 
Step: 6 of 100, time = 5.000000e-03 
Step: 7 of 100, time = 6.000000e-03 
Step: 8 of 100, time = 7.000000e-03 
Step: 9 of 100, time = 8.000000e-03 
Step: 10 of 100, time = 9.000000e-03 
Step: 11 of 100, time = 1.000000e-02 
Step: 12 of 100, time = 1.100000e-02 
Step: 13 of 100, time = 1.200000e-02 
Step: 14 of 100, time = 1.300000e-02 
Step: 15 of 100, time = 1.400000e-02 
Step: 16 of 100, time = 1.500000e-02 
Step: 17 of 100, time = 1.600000e-02 
Step: 18 of 100, time = 1.700000e-02 
Step: 19 of 100, time = 1.800000e-02 
Step: 20 of 100, time = 1.900000e-02 
Step: 21 of 100, time = 2.000000e-02 
Step: 22 of 100, time = 2.100000e-02 
Step: 23 of 100, time = 2.200000e-02 
Step: 24 of 100, time = 2.300000e-02 
Step: 25 of 100, time = 2.400000e-02 
Step: 26 of 100, time = 2.500000e-02 
Step: 27 of 100, time = 2.600000e-02 
Step: 28 of 100, time = 2.700000e-02 
Step: 29 of 100, time = 2.800000e-02 
Step: 30 of 100, time = 2.900000e-02 
Step: 31 of 100, time = 3.000000e-02 
Step: 32 of 100, time = 3.100000e-02 
Step: 33 of 100, time = 3.200000e-02 
Step: 34 of 100, time = 3.300000e-02 
Step: 35 of 100, time = 3.400000e-02 
Step: 36 of 100, time = 3.500000e-02 
Step: 37 of 100, time = 3.600000e-02 
Step: 38 of 100, time = 3.700000e-02 
Step: 39 of 100, time = 3.800000e-02 
Step: 40 of 100, time = 3.900000e-02 
Step: 41 of 100, time = 4.000000e-02 
Step: 42 of 100, time = 4.100000e-02 
Step: 43 of 100, time = 4.200000e-02 
Step: 44 of 100, time = 4.300000e-02 
Step: 45 of 100, time = 4.400000e-02 
Step: 46 of 100, time = 4.500000e-02 
Step: 47 of 100, time = 4.600000e-02 
Step: 48 of 100, time = 4.700000e-02 
Step: 49 of 100, time = 4.800000e-02 
Step: 50 of 100, time = 4.900000e-02 
Step: 51 of 100, time = 5.000000e-02 
Step: 52 of 100, time = 5.100000e-02 
Step: 53 of 100, time = 5.200000e-02 
Step: 54 of 100, time = 5.300000e-02 
Step: 55 of 100, time = 5.400000e-02 
Step: 56 of 100, time = 5.500000e-02 
Step: 57 of 100, time = 5.600000e-02 
Step: 58 of 100, time = 5.700000e-02 
Step: 59 of 100, time = 5.800000e-02 
Step: 60 of 100, time = 5.900000e-02 
Step: 61 of 100, time = 6.000000e-02 
Step: 62 of 100, time = 6.100000e-02 
Step: 63 of 100, time = 6.200000e-02 
Step: 64 of 100, time = 6.300000e-02 
Step: 65 of 100, time = 6.400000e-02 
Step: 66 of 100, time = 6.500000e-02 
Step: 67 of 100, time = 6.600000e-02 
Step: 68 of 100, time = 6.700000e-02 
Step: 69 of 100, time = 6.800000e-02 
Step: 70 of 100, time = 6.900000e-02 
Step: 71 of 100, time = 7.000000e-02 
Step: 72 of 100, time = 7.100000e-02 
Step: 73 of 100, time = 7.200000e-02 
Step: 74 of 100, time = 7.300000e-02 
Step: 75 of 100, time = 7.400000e-02 
Step: 76 of 100, time = 7.500000e-02 
Step: 77 of 100, time = 7.600000e-02 
Step: 78 of 100, time = 7.700000e-02 
Step: 79 of 100, time = 7.800000e-02 
Step: 80 of 100, time = 7.900000e-02 
Step: 81 of 100, time = 8.000000e-02 
Step: 82 of 100, time = 8.100000e-02 
Step: 83 of 100, time = 8.200000e-02 
Step: 84 of 100, time = 8.300000e-02 
Step: 85 of 100, time = 8.400000e-02 
Step: 86 of 100, time = 8.500000e-02 
Step: 87 of 100, time = 8.600000e-02 
Step: 88 of 100, time = 8.700000e-02 
Step: 89 of 100, time = 8.800000e-02 
Step: 90 of 100, time = 8.900000e-02 
Step: 91 of 100, time = 9.000000e-02 
Step: 92 of 100, time = 9.100000e-02 
Step: 93 of 100, time = 9.200000e-02 
Step: 94 of 100, time = 9.300000e-02 
Step: 95 of 100, time = 9.400000e-02 
Step: 96 of 100, time = 9.500000e-02 
Step: 97 of 100, time = 9.600000e-02 
Step: 98 of 100, time = 9.700000e-02 
Step: 99 of 100, time = 9.800000e-02 
Step: 100 of 100, time = 9.900000e-02 
Elapsed time is 6.927538 seconds.
if not(istest)
    t = opts.splitting.time_steps;
    figure;
    plot(t, out.ms);
    title('Ranks of approximations over time');

    y = zeros(1,length(out.Ks));
    for i=1:length(out.Ks)
        y(i) = out.Ks{i}(1,1);
    end

    figure;
    plot(t, y);
    title('evolution of component (1,1) of the optimal feedback');
else

    if abs(norm(out.Ds{1}) / 6.078945091766749e-05 - 1) >= 1e-10
       error('MESS:TEST:accuracy','unexpectedly inaccurate result');
   end
end
ans = 

  struct with fields:

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