Contents

function LQR_DAE3_SO(model, istest)
% Computes a Riccati feedback control for the constrained vibrating 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. https://doi.org/10.1007/3-540-27909-1_3

%
% This file is part of the M-M.E.S.S. project
% (http://www.mpi-magdeburg.mpg.de/projects/mess).
% Copyright (c) 2009-2025 Jens Saak, Martin Koehler, Peter Benner and others.
% All rights reserved.
% License: BSD 2-Clause License (see COPYING)
%
narginchk(0, 2);

if nargin < 1
    model = 'Stykel_small';
end

if nargin < 2
    istest = false;
end

set operation

opts = struct();
[oper, opts] = operatormanager(opts, '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 = true;
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 = false;
opts.nm.linesearch = true;
opts.nm.projection.freq = 0;
opts.nm.projection.ortho = true;
opts.nm.res = struct('maxiter', 10, 'tol', 1e-6, 'orth', 0);

The actual Newton call

eqn.type = 'T';

t_mess_lrnm = tic;
outnm = mess_lrnm(eqn, opts, oper);
t_elapsed1 = toc(t_mess_lrnm);
mess_fprintf(opts, 'mess_lrnm took %6.2f seconds \n', t_elapsed1);

if istest
    if min(outnm.res) >= opts.nm.res_tol
        mess_err(opts, 'TEST:accuracy', 'unexpectedly inaccurate result');
    end
else
    figure(1);
    semilogy(outnm.res, 'LineWidth', 3);
    title('0= C^T C + A^T X E + E^T X A -E^T X BB^T X E');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
mess_fprintf(opts, '2-Norm of the resulting feedback matrix: %g\n', ...
             norm(outnm.Z, 2));

mess_fprintf(opts, 'size outnm.Z: %d x %d\n\n', ...
             size(outnm.Z, 1), size(outnm.Z, 2));
Warning: Non-stable Ritz values were detected.
These will be removed from the set for the subsequent computations. 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/DAE3_SO/LQR_DAE3_SO.m',101)">LQR_DAE3_SO (line 101)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_lrnm.m',804)">mess_lrnm (line 804)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/shifts/mess_para.m',279)">mess_para (line 279)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/usfs/dae_3_so/get_ritz_vals_dae_3_so.m',74)">get_ritz_vals_dae_3_so (line 74)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/shifts/mess_get_ritz_vals.m',113)">mess_get_ritz_vals (line 113)</a>
  
 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 for the subsequent computations. 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/DAE3_SO/LQR_DAE3_SO.m',101)">LQR_DAE3_SO (line 101)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_lrnm.m',804)">mess_lrnm (line 804)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/shifts/mess_para.m',279)">mess_para (line 279)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/usfs/dae_3_so/get_ritz_vals_dae_3_so.m',74)">get_ritz_vals_dae_3_so (line 74)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/shifts/mess_get_ritz_vals.m',113)">mess_get_ritz_vals (line 113)</a>
  
 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   
  
mess_lrnm took   0.55 seconds   
2-Norm of the resulting feedback matrix: 2.27051  
size outnm.Z: 1200 x 50  
  

Lets try the RADI method 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 = true;
opts.radi.get_ZZt = true;
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;

t_mess_lrradi = tic;
outradi = mess_lrradi(eqn, opts, oper);
t_elapsed2 = toc(t_mess_lrradi);
mess_fprintf(opts, 'mess_lrradi took %6.2f seconds \n', t_elapsed2);

if istest
    if min(outradi.res) >= opts.radi.res_tol
        mess_err(opts, 'TEST:accuracy', 'unexpectedly inaccurate result');
    end
else
    figure(2);
    semilogy(outradi.res, 'LineWidth', 3);
    title('0= C^T C + A^T X E + E^T X A -E^T X BB^T X E');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
end

mess_fprintf(opts, 'size outradi.Z: %d x %d\n', ...
             size(outradi.Z, 1), size(outradi.Z, 2));
Warning: Non-stable Ritz values were detected.
These will be removed from the set for the subsequent computations. 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/DAE3_SO/LQR_DAE3_SO.m',148)">LQR_DAE3_SO (line 148)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_lrradi.m',874)">mess_lrradi (line 874)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/shifts/mess_lrradi_get_shifts.m',60)">mess_lrradi_get_shifts (line 60)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/shifts/mess_para.m',279)">mess_para (line 279)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/usfs/dae_3_so/get_ritz_vals_dae_3_so.m',74)">get_ritz_vals_dae_3_so (line 74)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/shifts/mess_get_ritz_vals.m',113)">mess_get_ritz_vals (line 113)</a>
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  
mess_lrradi took   0.28 seconds   
size outradi.Z: 1200 x 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-', ...
             'LineWidth', 3);
    title('0= C^TC + A^T X E + E^T X A - E^T X BB^T X E');
    xlabel('number of solves with A+p*M');
    ylabel('normalized residual norm');
    legend('LR-NM', 'RADI');
end