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.

[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.

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

Input checks

narginchk(0,5);
if nargin==0
    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 = 0;
end

set operation manager for the structured computations of Gramians

oper = operatormanager('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=1;
        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;
        beta=alpha;
        v=5e0;

        [eqn.M_,eqn.E_,eqn.K_]=triplechain_MSD(n1,alpha,beta,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)=-1;   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=1;
        eqn.alpha = -0.02;
    case 'tv'
        n1=151; % make sure it is odd!
        alpha=0.01;
        beta=alpha;
        v=5e0;

        [eqn.M_,eqn.E_,eqn.K_]=triplechain_MSD(n1,alpha,beta,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=1;
        eqn.alpha = -0.02;
   case 'truhar_veselic'
        n1=1500; % make sure it is even!
        alpha=0.01;
        beta=alpha;
        v=5e0;

        [eqn.M_,eqn.E_,eqn.K_]=triplechain_MSD(n1,alpha,beta,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=1;
        eqn.alpha = -0.02;
    otherwise
        fprintf('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;
tic;
opts.shifts.p = 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= opts.shifts.p-0.5;
end

outC = mess_lradi(eqn, opts, oper);
toc;

% observability Gramian
eqn.type = 'N';
tic;
outB = mess_lradi(eqn, opts, oper);
toc;
Warning: Could not compute initial projection shifts. Going to retry with random
right hand side. 
ADI step:    1 normalized residual: 1.053678e+01 relative change in Z: 1.000000e+00
ADI step:    2 normalized residual: 3.012160e+00 relative change in Z: 5.506374e-01
ADI step:    3 normalized residual: 1.697071e+00 relative change in Z: 1.913104e-01
ADI step:    4 normalized residual: 1.852880e+00 relative change in Z: 1.309126e-01
ADI step:    5 normalized residual: 2.112116e+00 relative change in Z: 1.550425e-01
ADI step:    7 normalized residual: 1.902741e-01 relative change in Z: 4.235858e-01
ADI step:    9 normalized residual: 7.157019e-02 relative change in Z: 1.763331e-01
ADI step:   11 normalized residual: 1.084505e-02 relative change in Z: 6.989401e-02
ADI step:   12 normalized residual: 1.142861e-02 relative change in Z: 1.831825e-02
ADI step:   13 normalized residual: 1.836154e-03 relative change in Z: 7.700851e-03
ADI step:   15 normalized residual: 1.445044e-04 relative change in Z: 1.823306e-02
ADI step:   17 normalized residual: 1.060878e-04 relative change in Z: 3.933767e-03
ADI step:   19 normalized residual: 1.178264e-04 relative change in Z: 3.259449e-03
ADI step:   21 normalized residual: 8.750229e-05 relative change in Z: 2.257048e-03
ADI step:   23 normalized residual: 3.022722e-05 relative change in Z: 1.388048e-03
ADI step:   25 normalized residual: 4.943924e-06 relative change in Z: 6.985186e-04
ADI step:   27 normalized residual: 3.878335e-07 relative change in Z: 5.062565e-04
ADI step:   29 normalized residual: 1.495482e-07 relative change in Z: 2.855021e-04
ADI step:   31 normalized residual: 8.338307e-09 relative change in Z: 8.570238e-05
ADI step:   33 normalized residual: 1.242339e-11 relative change in Z: 1.508670e-05
Elapsed time is 0.319768 seconds.
ADI step:    1 normalized residual: 1.558938e-01 relative change in Z: 1.000000e+00
ADI step:    2 normalized residual: 4.401792e-02 relative change in Z: 5.072038e-01
ADI step:    3 normalized residual: 7.352801e-02 relative change in Z: 5.006779e-01
ADI step:    4 normalized residual: 5.158659e-02 relative change in Z: 3.600499e-01
ADI step:    6 normalized residual: 5.117526e-03 relative change in Z: 4.068991e-01
ADI step:    8 normalized residual: 1.161967e-03 relative change in Z: 1.250397e-01
ADI step:   10 normalized residual: 1.973618e-04 relative change in Z: 5.280161e-02
ADI step:   12 normalized residual: 1.472261e-05 relative change in Z: 3.107334e-02
ADI step:   14 normalized residual: 1.073122e-05 relative change in Z: 7.185637e-03
ADI step:   16 normalized residual: 6.855427e-06 relative change in Z: 2.382091e-03
ADI step:   18 normalized residual: 1.873297e-06 relative change in Z: 1.566160e-03
ADI step:   20 normalized residual: 1.311802e-07 relative change in Z: 1.011573e-03
ADI step:   22 normalized residual: 1.071308e-08 relative change in Z: 6.292104e-04
ADI step:   24 normalized residual: 7.695432e-11 relative change in Z: 1.391128e-04
Elapsed time is 0.144138 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;

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));
toc;
Computing reduced order model via square root method

reduced system order: 9  (max possible/allowed: 24/500)

Elapsed time is 0.010480 seconds.

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

ROM and comparison to the original model.

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

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

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 = 1;

operu = operatormanager('so_1');

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

toc;
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

Elapsed time is 2.708875 seconds.

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
    % wether they are "good enough"...
    if (max(err) > 50*tol)
        error('MESS:TEST:accuracy', ['unexpectedly inaccurate result ' ...
                            'for %s %g %d %d (%g)'], model, tol, ...
              max_ord, maxiter,max(err));
    end
end