Contents

function BT_DAE3_SO(model, tol, max_ord, maxiter, istest)

BT_DAE3_SO demonstrates the use of balanced truncation for index-3 second order models. It reports the results presented in [3].

Usage:
 BT_DAE3_SO(model,tol,max_ord,maxiter,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
         'Truhar_Veselic' a version of the triple chain mass spring damper
                          test example from [2] with holonomic constraints
                          coupling certain masses. Dimension 451 in second
                          order form.
         'TV'             another version of the example from [2] with
                          other constraints
         'TV2'            another version of the example from [2] with
                          yet another set of constraints
tol       Balanced Truncation tolerance (optional, default: 1e-4)
max_ord   maximal reduced order allowed (optional, default: 500)
maxiter   maximal number of ADI iterations (optional, default: 300)
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

[2] N. Truhar and K. Veselic, An efficient method for estimating the optimal dampers’ viscosity for linear vibrating systems using Lyapunov equation, SIAM J. Matrix Anal. Appl., 31 (2009), pp. 18–39. https://doi.org/10.1137/070683052

[3] J. Saak, M. Voigt, Model reduction of constrained mechanical systems in M-M.E.S.S., IFAC-PapersOnLine 9th Vienna International Conference on Mathematical Modelling MATHMOD 2018, Vienna, Austria, 21–23 February 2018 51 (2) (2018) 661–666. https://doi.org/10.1016/j.ifacol.2018.03.112

%
% 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)
%

Input checks

narginchk(0, 5);

if nargin < 1
    model = 'Stykel_small';
end

if nargin < 2
    tol = 1e-4;
end

if nargin < 3
    max_ord = 500;
end

if nargin < 4
    maxiter = 300;
end

if nargin < 5
    istest = false;
end

set operation manager for the structured computations of Gramians

opts = struct();
[oper, opts] = operatormanager(opts, 'dae_3_so');

load problem data

switch lower(model)
    case {'stykel_small', 'stykel_large'}
        if strcmp(model, 'stykel_small')
            sys = load(sprintf('%s/../models/ms_ind3_by_t_stykel/g600.mat', ...
                               fileparts(mfilename('fullpath'))));
        else
            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 = full(sys.B(1:2 * nv, :));
        eqn.C = full(sys.C(:, 1:2 * nv));
        clear E A B C M D K G;
    case 'tv2'
        n1 = 151; % make sure it is odd!
        alpha = 0.01;
        v = 5e0;

        [eqn.M_, eqn.E_, eqn.K_] = triplechain_MSD(n1, alpha, alpha, v);
        eqn.E_ = -eqn.E_;
        eqn.K_ = -eqn.K_;
        eqn.G_ = sparse(3, 3 * n1 + 1);
        n12 = ceil(n1 / 2); % n1 is odd so this is the index of
        % the center of the string
        eqn.G_(1, 1) = -1;
        eqn.G_(1, n12) = 2;
        eqn.G_(1, n1) = -1;
        eqn.G_(2, n1 + 1) = -1;
        eqn.G_(2, n1 + n12) = 2;
        eqn.G_(2, 2 * n1) = -1;
        eqn.G_(3, 2 * n1 + 1) = -1;
        eqn.G_(3, 2 * n1 + n12) = 2;
        eqn.G_(3, 3 * n1) = -1;

        nv = size(eqn.M_, 1);
        np = size(eqn.G_, 1);
        eqn.B = zeros(6 * n1 + 2, 3);
        eqn.B(3 * n1 + 6, 1) = 1;
        eqn.B(5 * n1 + n12 + 6, 2) = 1;
        eqn.B(end - 6, 3) = 1;
        eqn.C = zeros(3, 6 * n1 + 2);
        eqn.C(1, 3 * n1 + 1) = 1;
        eqn.C(2, 2 * n1) = 1;
        eqn.C(3, 2 * n1 + floor(n12 / 2)) = 1;
        eqn.haveE = true;
        eqn.alpha = -0.02;
    case 'tv'
        n1 = 151; % make sure it is odd!
        alpha = 0.01;
        v = 5e0;

        [eqn.M_, eqn.E_, eqn.K_] = triplechain_MSD(n1, alpha, alpha, v);
        eqn.E_ = -eqn.E_;
        eqn.K_ = -eqn.K_;
        eqn.G_ = sparse(3, 3 * n1 + 1);
        n12 = ceil(n1 / 2); % n1 is odd so this is the index of
        % the center of the string
        eqn.G_(1, 1) = 1 / 3;
        eqn.G_(1, n12) = 1 / 3;
        eqn.G_(1, n1) = 1 / 3;
        eqn.G_(2, n1 + 1) = 1 / 3;
        eqn.G_(2, n1 + n12) = 1 / 3;
        eqn.G_(2, 2 * n1) = 1 / 3;
        eqn.G_(3, 2 * n1 + 1) = 1 / 3;
        eqn.G_(3, 2 * n1 + n12) = 1 / 3;
        eqn.G_(3, 3 * n1) = 1 / 3;

        nv = size(eqn.M_, 1);
        np = size(eqn.G_, 1);
        eqn.B = [zeros(3 * n1 + 1, 1); ones(3 * n1 + 1, 1)];
        eqn.C = [ones(3 * n1 + 1, 1); zeros(3 * n1 + 1, 1)]';
        eqn.haveE = true;
        eqn.alpha = -0.02;
    case 'truhar_veselic'
        n1 = 1500; % make sure it is even!
        alpha = 0.01;
        v = 5e0;

        [eqn.M_, eqn.E_, eqn.K_] = triplechain_MSD(n1, alpha, alpha, v);
        eqn.E_ = -eqn.E_;
        eqn.K_ = -eqn.K_;
        eqn.G_ = sparse(6, 3 * n1 + 1);
        eqn.G_(1, 1) = 1;
        eqn.G_(1, n1 / 2) = -1;
        eqn.G_(2, n1 / 2 + 1) = 1;
        eqn.G_(2, n1) = -1;
        eqn.G_(3, n1 + 1) = 1;
        eqn.G_(3, n1 + n1 / 2) = -1;
        eqn.G_(4, n1 + n1 / 2 + 1) = 1;
        eqn.G_(4, 2 * n1) = -1;
        eqn.G_(5, 2 * n1 + 1) = 1;
        eqn.G_(5, 2 * n1 + n1 / 2) = -1;
        eqn.G_(6, 2 * n1 + n1 / 2 + 1) = 1;
        eqn.G_(6, 3 * n1) = -1;
        nv = size(eqn.M_, 1);
        np = size(eqn.G_, 1);
        eqn.B = [zeros(3 * n1 + 1, 1); ones(3 * n1 + 1, 1)];
        eqn.C = [ones(3 * n1 + 1, 1); zeros(3 * n1 + 1, 1)]';
        eqn.haveE = true;
        eqn.alpha = -0.02;
    otherwise
        mess_err(opts, 'unknown model requested!\n');
        return
end

options

Adi options

opts.adi.maxiter = maxiter;
opts.adi.res_tol = 1e-10;
opts.adi.rel_diff_tol = 1e-11;
opts.norm = 'fro';
opts.shifts.method = 'projection';
opts.shifts.num_desired = 25;

LRADI for the two Gramian factors

controllability Gramian
eqn.type = 'T';
opts.adi.info = 1;
t_mess_lradi1 = tic;
[p, ~, eqn, opts, oper] = mess_para(eqn, opts, oper);
% use an additional alpha-shift to improve convergence and ROM quality for
% the triple chain model
if strcmp(model, 'Truhar_Veselic') || strcmp(model, 'TV') || strcmp(model, 'TV2')
    opts.shifts.p = p - 0.5;
else
    opts.shifts.p = p;
end

outC = mess_lradi(eqn, opts, oper);

t_elapsed1 = toc(t_mess_lradi1);
mess_fprintf(opts, 'mess_lradi took %6.2f seconds \n\n', t_elapsed1);

% observability Gramian
eqn.type = 'N';
t_mess_lradi2 = tic;
outB = mess_lradi(eqn, opts, oper);
t_elapsed2 = toc(t_mess_lradi2);
mess_fprintf(opts, 'mess_lradi took %6.2f seconds \n', t_elapsed2);
Warning: Could not compute initial projection shifts. Going to retry with random
right hand side. 
↳ 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/BT_DAE3_SO.m',214)">BT_DAE3_SO (line 214)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/shifts/mess_para.m',379)">mess_para (line 379)</a>
ADI step:    1 normalized residual: 1.215480e+01 relative change in Z: 1.000000e+00  
ADI step:    2 normalized residual: 3.525543e+00 relative change in Z: 5.520973e-01  
ADI step:    3 normalized residual: 1.953222e+00 relative change in Z: 1.898487e-01  
ADI step:    4 normalized residual: 2.131585e+00 relative change in Z: 1.324398e-01  
ADI step:    5 normalized residual: 2.435279e+00 relative change in Z: 1.566566e-01  
ADI step:    7 normalized residual: 2.191473e-01 relative change in Z: 4.227025e-01  
ADI step:    9 normalized residual: 8.229267e-02 relative change in Z: 1.761189e-01  
ADI step:   11 normalized residual: 1.247249e-02 relative change in Z: 6.986264e-02  
ADI step:   12 normalized residual: 1.309187e-02 relative change in Z: 1.825205e-02  
ADI step:   13 normalized residual: 2.110282e-03 relative change in Z: 7.668970e-03  
ADI step:   15 normalized residual: 1.666456e-04 relative change in Z: 1.817214e-02  
ADI step:   17 normalized residual: 1.223639e-04 relative change in Z: 3.929405e-03  
ADI step:   19 normalized residual: 1.359440e-04 relative change in Z: 3.259148e-03  
ADI step:   21 normalized residual: 1.010439e-04 relative change in Z: 2.257748e-03  
ADI step:   23 normalized residual: 3.494455e-05 relative change in Z: 1.389120e-03  
ADI step:   25 normalized residual: 5.723896e-06 relative change in Z: 6.988380e-04  
ADI step:   27 normalized residual: 4.484029e-07 relative change in Z: 5.062538e-04  
ADI step:   29 normalized residual: 1.727537e-07 relative change in Z: 2.856449e-04  
ADI step:   31 normalized residual: 9.661411e-09 relative change in Z: 8.582632e-05  
ADI step:   33 normalized residual: 1.434879e-11 relative change in Z: 1.510282e-05  
mess_lradi took   0.47 seconds   
  
ADI step:    1 normalized residual: 3.133718e-01 relative change in Z: 1.000000e+00  
ADI step:    2 normalized residual: 8.672847e-02 relative change in Z: 5.044804e-01  
ADI step:    3 normalized residual: 1.460907e-01 relative change in Z: 5.012199e-01  
ADI step:    4 normalized residual: 1.035263e-01 relative change in Z: 3.622584e-01  
ADI step:    6 normalized residual: 1.030719e-02 relative change in Z: 4.069064e-01  
ADI step:    8 normalized residual: 2.348522e-03 relative change in Z: 1.250323e-01  
ADI step:   10 normalized residual: 3.960126e-04 relative change in Z: 5.286041e-02  
ADI step:   12 normalized residual: 2.940750e-05 relative change in Z: 3.112750e-02  
ADI step:   14 normalized residual: 2.145534e-05 relative change in Z: 7.172585e-03  
ADI step:   16 normalized residual: 1.375092e-05 relative change in Z: 2.387882e-03  
ADI step:   18 normalized residual: 3.776225e-06 relative change in Z: 1.570942e-03  
ADI step:   20 normalized residual: 2.663151e-07 relative change in Z: 1.012300e-03  
ADI step:   22 normalized residual: 2.130823e-08 relative change in Z: 6.299396e-04  
ADI step:   24 normalized residual: 1.533744e-10 relative change in Z: 1.394598e-04  
ADI step:   26 normalized residual: 3.735751e-11 relative change in Z: 2.349302e-05  
mess_lradi took   0.22 seconds   

Reduced Order Model computation via square root method (SRM)

fprintf('\nComputing reduced order model via square root method\n\n');
opts.srm.tol = tol;
opts.srm.max_ord = max_ord;
opts.srm.info = 1;

t_SRM_ROM = tic;
[TL, TR] = mess_square_root_method(eqn, opts, oper, outB.Z, outC.Z);
% compute ROM matrices
ROM.A = TL' * oper.mul_A(eqn, opts, 'N', TR, 'N');
ROM.B = TL' * eqn.B;
ROM.C = eqn.C * TR;
ROM.E = eye(size(ROM.A));
t_elapsed3 = toc(t_SRM_ROM);
mess_fprintf(opts, 'ROM matrices computation took %6.2f seconds \n\n', ...
             t_elapsed3);
Computing reduced order model via square root method

reduced system order: 9  (max possible/allowed: 26/500)  
  
ROM matrices computation took   0.02 seconds   
  

Frequency-domain evaluation of the (transfer function of the)

ROM and comparison to the original model.

We feed the mess_tf_plot with usfs that do not exploit the DAE structure:

t_FD_eval = tic;
opts.tf_plot.nsample = 200;
if istest
    opts.tf_plot.info = 0;
else
    opts.tf_plot.info = 2;
end
if strcmp(model, 'TV2')
    opts.tf_plot.fmin = -2;
    opts.tf_plot.fmax = 1;
else
    opts.tf_plot.fmin = -4;
    opts.tf_plot.fmax = 4;
end

opts.tf_plot.type = 'sigma';

NG = sparse(np, nv);
NS = sparse(np, np);
eqnu.M_ = [eqn.M_ NG'; NG NS];
eqnu.E_ = [-eqn.E_ NG'; NG NS];
eqnu.K_ = [-eqn.K_ -eqn.G_'; -eqn.G_ NS];
eqnu.C =  [zeros(size(eqn.C, 1), np + nv), eqn.C(:, 1:nv), zeros(size(eqn.C, 1), np)];
eqnu.B =  [eqn.B(nv + 1:2 * nv, :); zeros(2 * np + nv, size(eqn.B, 2))];
eqnu.haveE = true;

[operu, opts] = operatormanager(opts, 'so_1');

out = mess_tf_plot(eqnu, opts, operu, ROM);
err = out.err;

t_elapsed4 = toc(t_FD_eval);
mess_fprintf(opts, 'frequency-domain evaluation took %6.2f \n\n', t_elapsed4);
Computing TFMs of original and reduced order systems and MOR errors
 Step  20 / 200 Step  40 / 200 Step  60 / 200 Step  80 / 200 Step 100 / 200 Step 120 / 200 Step 140 / 200 Step 160 / 200 Step 180 / 200 Step 200 / 200

frequency-domain evaluation took   3.13   
  

final accuracy test used in the continuous integration system or

plot of the computed
if istest
    % the errors are not always perfect in this example, but let's see
    % whether they are "good enough"...
    if max(err) > (50 * tol)
        mess_err(opts, 'TEST:accuracy', ...
                 'unexpectedly inaccurate result for %s %g %d %d (%g)', ...
                 model, tol, max_ord, maxiter, max(err));
    end
end