Contents

function LQR_DAE3_SO(model,istest)
% Computes a Riccati feedback control for the constrained vibrational model
% from [1]
%
%
%   Usage:
%    LQR_DAE3_SO(model,test)
%
% Input:
%
%  model    choice of the example system. Possible values
%           'Stykel_small'   Stykel's mass spring damper system from [1]
%                            of dimension 600 in second order form (default)
%           'Stykel_large'   Stykel's mass spring damper system from [1]
%                            of dimension 6000 in second order form
%
%  istest    decides whether the function runs as an interactive demo or a
%            continuous integration test.
%            (optional; defaults to 0, i.e. interactive demo)
%
% References:
% [1] V. Mehrmann and T. Stykel. Balanced truncation model reduction for
%     large-scale systems in descriptor form.
%     In P. Benner, V. Mehrmann, and D. Sorensen, editors, Dimension
%     Reduction of Large-Scale Systems, volume 45 of Lecture Notes in
%     Computational Science and Engineering, pages 83–115. Springer-Verlag,
%     Berlin/Heidelberg, 2005.

% 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,2);

if (nargin < 1)
    model = 'Stykel_small';
end

if ( nargin < 2 )
    istest = 0;
end

set operation

oper = operatormanager('dae_3_so');

load problem data

switch lower(model)
    case 'stykel_small'
             sys = load(sprintf('%s/../models/ms_ind3_by_t_stykel/g600.mat',...
                fileparts(mfilename('fullpath'))));
    case 'stykel_large'
            sys = load(sprintf('%s/../models/ms_ind3_by_t_stykel/g6000.mat',...
                fileparts(mfilename('fullpath'))));
end
eqn.M_=sys.M;
eqn.E_=sys.D;
eqn.K_=sys.K;
eqn.G_=sys.G;
eqn.haveE=1;
eqn.alpha = -0.02;
nv = size(eqn.M_,1);
np = size(eqn.G_,1);
eqn.B = sys.B(1:2*nv,:);
eqn.C = sys.C(:,1:2*nv);
eqn.type = 'T';

First we run the Newton-ADI Method

opts.norm = 'fro';

% ADI options
opts.adi.maxiter = 1000;
opts.adi.res_tol = 1e-15;
opts.adi.rel_diff_tol = 1e-16;
opts.adi.info = 0;
opts.shifts.num_desired=25;
opts.shifts.num_Ritz=50;
opts.shifts.num_hRitz=25;
opts.shifts.b0=ones(2*nv+np,1);

% Newton options and maximum iteration number
opts.nm.maxiter = 20;
opts.nm.res_tol = 1e-10;
opts.nm.rel_diff_tol = 1e-16;
opts.nm.info = 1;
opts.nm.accumulateRes = 0;
opts.nm.linesearch = 1;
opts.nm.projection.freq=0;
opts.nm.projection.ortho=1;
opts.nm.res=struct('maxiter',10,'tol',1e-6,'orth',0);

The actual Newton call

eqn.type = 'T';

tic;
outnm = mess_lrnm(eqn, opts, oper);
toc;

if istest
    if min(outnm.res)>=opts.nm.res_tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result');
    end
else
    figure(1);
    semilogy(outnm.res);
    title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
fprintf('2-Norm of the resulting feedback matrix: %g\n', norm(outnm.Z,2));

disp('size outnm.Z:');
disp(size(outnm.Z));
Warning: Non-stable Ritz values were detected.
 These will be removed from the set in further computations. 

NM step:    1  normalized residual: 1.415415e-01
               relative change in K: 1.000000e+00
               number of ADI steps: 46 

Warning: Non-stable Ritz values were detected.
 These will be removed from the set in further computations. 

NM step:    2  normalized residual: 2.631219e-05
               relative change in K: 1.382088e-02
               number of ADI steps: 46 


NM step:    3  normalized residual: 4.572588e-13
               relative change in K: 1.822007e-06
               number of ADI steps: 46 

Elapsed time is 0.453428 seconds.
2-Norm of the resulting feedback matrix: 2.27051
size outnm.Z:
        1200          50

Lets try the RADI metod and compare

RADI-MESS settings

opts.shifts.history = opts.shifts.num_desired*size(eqn.C,1);
opts.shifts.num_desired=25;
opts.shifts.num_Ritz=50;
opts.shifts.num_hRitz=25;
opts.shifts.b0=ones(2*nv+np,1);

% choose either of the three shift methods, here
%opts.shifts.method = 'gen-ham-opti';
opts.shifts.method = 'heur';
%opts.shifts.method = 'projection';

opts.shifts.naive_update_mode = false; % .. Suggest false (smart update is faster; convergence is the same).
opts.radi.compute_sol_fac = 1;
opts.radi.get_ZZt = 1;
opts.radi.maxiter = opts.adi.maxiter;
opts.norm = 2;
opts.radi.res_tol = opts.nm.res_tol;
opts.radi.rel_diff_tol = 0;
opts.radi.info = 1;

tic;
outradi = mess_lrradi(eqn, opts, oper);
toc;

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

disp('size outradi.Z:');
disp(size(outradi.Z));
Warning: Non-stable Ritz values were detected.
 These will be removed from the set in further computations. 
RADI step:    1 normalized residual: 1.001499e+00
RADI step:    2 normalized residual: 4.884287e+00
RADI step:    4 normalized residual: 2.281476e+00
RADI step:    6 normalized residual: 3.593905e+00
RADI step:    8 normalized residual: 1.014958e+00
RADI step:   10 normalized residual: 1.425091e-01
RADI step:   12 normalized residual: 1.952333e-01
RADI step:   14 normalized residual: 7.656443e-02
RADI step:   16 normalized residual: 4.215812e-03
RADI step:   18 normalized residual: 2.983008e-03
RADI step:   20 normalized residual: 1.777674e-05
RADI step:   22 normalized residual: 7.529772e-07
RADI step:   24 normalized residual: 9.044695e-09
RADI step:   26 normalized residual: 5.459748e-10
RADI step:   27 normalized residual: 5.372400e-10
RADI step:   28 normalized residual: 9.524957e-10
RADI step:   30 normalized residual: 1.102023e-10
RADI step:   32 normalized residual: 1.669552e-10
RADI step:   34 normalized residual: 1.471539e-10
RADI step:   36 normalized residual: 1.558063e-11
Elapsed time is 0.191757 seconds.
size outradi.Z:
        1200          50

compare

if not(istest)
    figure(3);
    ls_nm=cumsum([outnm.adi.niter]);
    ls_radi=1:outradi.niter;

    semilogy(ls_nm,outnm.res,'k--',ls_radi,outradi.res,'b-');
    title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
    xlabel('number of solves with A+p*M');
    ylabel('normalized residual norm');
    legend('LR-NM','RADI');
end