Contents

function LQR_FDM_unstable(n0, n_unst, istest)
% Computes a stabilizing feedbacks for primal and dual systems by
% solving the algebraic Riccati equations with both the ZZ' and LDL'
% approximations of the solutions.
%
% Inputs:
%
% n0          n0^2 gives the dimension of the original model, i.e. n0 is
%             the number of degrees of freedom, i.e. grid points, per
%             spatial direction
%             (optional; defaults to 20)
%
% n_unst      number of unstable shifts by construction
%             (optional; defaults to 5)
%
% istest      flag to determine whether this demo runs as a CI test or
%             interactive demo
%             (optional, defaults to 0, i.e. interactive demo)
%

%
% 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,3);
if nargin<1, n0 = 20; end
if nargin<2, n_unst = 5; end
if nargin<3, istest = 0; end

construct unstable matrix

first generate the stable part

A0 = fdm_2d_matrix(n0, '10*x', '1000*y', '0');
n_stable = size(A0, 1);
n_unstable = n_unst; % no. of unstable eigenvalues
n = n_stable+n_unstable;
B = ones(n, 5);

% now add the instability by constructing a known antistable part and
% corresponding initial stabilizing feedback
Y = eye(n_unstable); % fix "projected feedback" (the Bernoulli/Lyap. sol)

% grab the B part corresponding to the antistable subsystem
Bt = B(n_stable+1 : n, : );
% initial feedback mirrors unstable eigenvals to negative halfplane, these
% may be computed as shifts and make (A+pI) singular
Vl = [zeros(n_stable,n_unstable); eye(n_unstable)];
K0 = Vl*Y*Bt;
if exist('lyap','file')% extra antistable part of A
    Au = lyap(Y, -Bt * Bt');
else
    Au = lyap2solve(Y, (-Bt * Bt')');
end


%set the full A and corresponding C together with the initial feedback
eqn.A_ = blkdiag(A0,Au);
opts.nm.K0 = K0';
C = ones(10, n);

eqn.haveE = 0;

% set operation
oper = operatormanager('default');

global options

opts.norm = 'fro';

ADI tolerances and maximum iteration number

opts.adi.maxiter = 200;
opts.adi.res_tol = 1e-10;
opts.adi.rel_diff_tol = 1e-16;
opts.adi.info = 0;
opts.adi.accumulateDeltaK = 1;
opts.adi.accumulateK = 0;
opts.adi.compute_sol_fac = 1;

shift parameters via projection

opts.shifts.num_desired = 5;
opts.shifts.method = 'projection';

Newton tolerances and maximum iteration number

opts.nm.maxiter = 15;
opts.nm.res_tol = 1e-8;
opts.nm.rel_diff_tol = 1e-16;
opts.nm.info = 1;
opts.nm.linesearch = 1;
opts.nm.accumulateRes = 1;
opts.nm.res=struct('maxiter',10,'tol',1e-6,'orth',0);

for shift banned eigenvalues

opts.shifts.banned = -eig(Au);
opts.shifts.banned_tol = 1e-6;

solve Riccati equation

for type=['T','N']
    eqn.type=type;
    eqn.B = B;
    eqn.C = C;
    tic;
    [outB, eqn1, opts1, oper1]=mess_lrnm(eqn, opts, oper);
    toc;

    disp('size outB.Z:');
    disp(size(outB.Z));
		Using line search (res: 9.067072e+01)
		lambda: 6.886611e-02

NM step:    1  normalized residual: 8.571260e-01
               relative change in K: 9.706812e-01
               number of ADI steps: 141 


		Using line search (res: 2.592935e+00)
		lambda: 3.842189e-01

NM step:    2  normalized residual: 5.018087e-01
               relative change in K: 7.920006e-01
               number of ADI steps: 106 


NM step:    3  normalized residual: 9.867003e-02
               relative change in K: 2.678522e-01
               number of ADI steps: 106 


NM step:    4  normalized residual: 2.025621e-03
               relative change in K: 3.911352e-02
               number of ADI steps: 96 


NM step:    5  normalized residual: 6.303703e-07
               relative change in K: 6.901770e-04
               number of ADI steps: 125 


NM step:    6  normalized residual: 8.081437e-09
               relative change in K: 1.567180e-07
               number of ADI steps: 86 

Elapsed time is 1.450450 seconds.
size outB.Z:
   405    58

		Using line search (res: 4.718175e+01)
		lambda: 1.218982e-01

NM step:    1  normalized residual: 5.671448e-01
               relative change in K: 9.673735e-01
               number of ADI steps: 99 


NM step:    2  normalized residual: 2.506721e-01
               relative change in K: 4.148393e-01
               number of ADI steps: 89 


NM step:    3  normalized residual: 1.602879e-02
               relative change in K: 1.131273e-01
               number of ADI steps: 90 


NM step:    4  normalized residual: 1.068051e-04
               relative change in K: 9.290895e-03
               number of ADI steps: 82 


NM step:    5  normalized residual: 1.078169e-08
               relative change in K: 9.335850e-05
               number of ADI steps: 98 


NM step:    6  normalized residual: 1.371685e-09
               relative change in K: 1.881478e-08
               number of ADI steps: 78 

Elapsed time is 1.340320 seconds.
size outB.Z:
   405    48

check residual norm

    if eqn.type == 'T'
        res1 = mess_res2_norms(outB.Z,'riccati',eqn1,opts1,oper1,opts1.nm,[]) ...
            / norm(eqn.C*eqn.C');
        res2 = abs(eigs(@(x) eqn.A_'*(outB.Z*((outB.Z'*x)))+(outB.Z*((outB.Z'*(eqn.A_*x))))...
            +eqn.C'*(eqn.C*x)-(outB.K)'*((outB.K)*x),n,1,'LM'))...
            /norm(eqn.C*eqn.C');
    else
        res1 = mess_res2_norms(outB.Z,'riccati',eqn1,opts1,oper1,opts1.nm,[]) ...
            / norm(eqn.B'*eqn.B);
        res2 = abs(eigs(@(x) eqn.A_*(outB.Z*((outB.Z'*x)))+(outB.Z*((outB.Z'*(eqn.A_'*x))))...
            +eqn.B*(eqn.B'*x)-(outB.K)'*((outB.K)*x),n,1,'LM'))...
            /norm(eqn.B'*eqn.B);
    end
    fprintf('Newton: %e \t mess_res2_norms: %e \t eigs: %e \n', ...
        [outB.res(end), res1, res2]);
Newton: 8.081437e-09 	 mess_res2_norms: 7.995688e-09 	 eigs: 7.995688e-09 
Newton: 1.371685e-09 	 mess_res2_norms: 1.359802e-09 	 eigs: 1.359802e-09 

print output

    if istest
        if min(outB.res)>=opts.nm.res_tol
            error('MESS:TEST:accuracy','unexpectedly inaccurate result');
        end
    else
        figure(1);
        semilogy(outB.res);
        title('0= C^TC + A^TXE + E^TXA -E^TXBB^TXE');
        xlabel('number of iterations');
        ylabel('normalized residual norm');
        pause(1);
    end

Test with LDL_T formulation

    fprintf('\n\n');
    disp('Repeat computations with LDL^T representation of data and solution');
    opts.LDL_T = 1;
    if eqn.type == 'T'
        eqn.S = diag([4,4,9,9,16,16,64,64,81,81]);
        eqn.C = diag([4,4,9,9,16,16,64,64,81,81].^(-0.5)) * eqn.C;
    else
        eqn.S = diag([4,4,9,9,16]);
        eqn.B = eqn.B *diag([4,4,9,9,16].^(-0.5));
    end

Repeat computations with LDL^T representation of data and solution

Repeat computations with LDL^T representation of data and solution

solve Riccati equation

    tic;
    [outB2, eqn2, opts2, oper2]=mess_lrnm(eqn, opts, oper);
    toc;
		Using line search (res: 9.067072e+01)
		lambda: 6.886611e-02

NM step:    1  normalized residual: 8.571260e-01
               relative change in K: 9.706812e-01
               number of ADI steps: 186 


		Using line search (res: 2.592935e+00)
		lambda: 3.842189e-01

NM step:    2  normalized residual: 5.018087e-01
               relative change in K: 7.920006e-01
               number of ADI steps: 108 


NM step:    3  normalized residual: 9.867003e-02
               relative change in K: 2.678522e-01
               number of ADI steps: 106 


NM step:    4  normalized residual: 2.025622e-03
               relative change in K: 3.911352e-02
               number of ADI steps: 96 


NM step:    5  normalized residual: 6.303703e-07
               relative change in K: 6.901770e-04
               number of ADI steps: 125 


NM step:    6  normalized residual: 8.081521e-09
               relative change in K: 1.567180e-07
               number of ADI steps: 86 

Elapsed time is 1.549733 seconds.
		Using line search (res: 4.718175e+01)
		lambda: 1.218982e-01

NM step:    1  normalized residual: 5.671448e-01
               relative change in K: 9.673735e-01
               number of ADI steps: 118 


NM step:    2  normalized residual: 2.506721e-01
               relative change in K: 4.148393e-01
               number of ADI steps: 93 


NM step:    3  normalized residual: 1.602879e-02
               relative change in K: 1.131273e-01
               number of ADI steps: 90 


NM step:    4  normalized residual: 1.068051e-04
               relative change in K: 9.290895e-03
               number of ADI steps: 88 


NM step:    5  normalized residual: 1.078169e-08
               relative change in K: 9.335850e-05
               number of ADI steps: 98 


NM step:    6  normalized residual: 1.372881e-09
               relative change in K: 1.881476e-08
               number of ADI steps: 78 

Elapsed time is 1.396881 seconds.

print output

    if istest
        if min(outB2.res)>=opts.nm.res_tol
            error('MESS:TEST:accuracy','unexpectedly inaccurate result');
        end
    else
        figure(2);
        semilogy(outB2.res);
        title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
        xlabel('number of iterations');
        ylabel('normalized residual norm');
        pause(1);
    end
    disp('size outB2.Z:');
    disp(size(outB2.Z));
size outB2.Z:
         405        1290

size outB2.Z:
         405        1170

check residual norm

    if eqn.type == 'T'
        res0 = norm(eqn.C'* eqn.S * eqn.C);
        eqn2.S_diag = diag(outB2.S);
        res1 = mess_res2_norms(outB2.Z,'riccati',eqn2,opts2,oper2,opts1.nm,outB2.D)/res0;
        outB2.ZD = outB2.Z*outB2.D;
        res2 = abs(eigs(@(x) eqn.A_'*(outB2.ZD*((outB2.Z'*x)))+(outB2.ZD*((outB2.Z'*(eqn.A_*x))))...
               +eqn.C'*eqn.S*(eqn.C*x)-(outB2.K)'*((outB2.K)*x),n,1,'LM'))...
               /res0;
    else
        res0 = norm((eqn.B * eqn.S) * eqn.B');
        eqn2.S_diag = diag(outB2.S);
        res1 = mess_res2_norms(outB2.Z,'riccati',eqn2,opts2,oper2,opts1.nm,outB2.D)/res0;
        outB2.ZD = outB2.Z*outB2.D;
        res2 = abs(eigs(@(x) eqn.A_*(outB2.ZD*((outB2.Z'*x)))+(outB2.ZD*((outB2.Z'*(eqn.A_'*x))))...
               +eqn.B*eqn.S*(eqn.B'*x)-(outB2.K)'*((outB2.K)*x),n,1,'LM'))...
               /res0;
    end
    fprintf('inner: %e \t mess_res2_norms: %e \t eigs: %e \n', ...
        [outB2.res(end), res1, res2]);
    % only needed to make the next iteration work with ZZ' again
    if opts.LDL_T, opts=rmfield(opts,'LDL_T'); end
inner: 8.081521e-09 	 mess_res2_norms: 7.995767e-09 	 eigs: 7.995767e-09 
inner: 1.372881e-09 	 mess_res2_norms: 1.360896e-09 	 eigs: 1.360896e-09 
end