Contents

function bt_mor_DAE1_tol(k, istest)
% Computes a reduced order model via Balanced Truncation for the proper
% index-1 System BIPS98_606 from https://sites.google.com/site/rommes/software
% following the method suggested in [1]
%
% Input:
% k         select model (allowed values 1,..5, 7, .., 13, default 7)
%           note that BIPS model 6 is not stable.
% 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] F. Freitas, J. Rommes, N. Martins, Gramian-based reduction method
%    applied to large sparse power system descriptor models, IEEE Trans.
%    Power Syst. 23 (3) (2008) 1258–1270. doi:10.1109/TPWRS.2008.926693

%
% This file is part of the M-M.E.S.S. project
% (http://www.mpi-magdeburg.mpg.de/projects/mess).
% Copyright © 2009-2022 Jens Saak, Martin Koehler, Peter Benner and others.
% All rights reserved.
% License: BSD 2-Clause License (see COPYING)
%
if nargin < 1, k = 7; end
if nargin < 2, istest = 0; end

set operation manager for the Gramian computations

oper = operatormanager('dae_1');

Read problem data

if k == 6
    warning('MESS:illegal_input', ...
        ['The Juba5723 model is not stable and thus not supported. ',...
         'Using BIPS bips98_606 instead.']);
    k = 7;
end
eqn = mess_get_BIPS(k);

Turn off close to singular warnings

(this model is really badly conditioned)

orig_warnstate = warning('OFF','MATLAB:nearlySingularMatrix');

ADI tolerances and maximum iteration number

opts.adi.maxiter = 300;
opts.adi.res_tol = 1e-10;
opts.adi.rel_diff_tol = 1e-16;
opts.adi.info = 1;
opts.norm = 'fro';
opts.shifts.method = 'projection';
opts.shifts.num_desired= 20;
eqn.type='N';
t_mess_lradi = tic;
outB = mess_lradi(eqn, opts, oper);
t_elapsed1 = toc(t_mess_lradi);
fprintf(1,'mess_lradi took %6.2f seconds \n', t_elapsed1);
if istest
    if min(outB.res)>=opts.adi.res_tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result');
    end
else
    figure;
    semilogy(outB.res,'LineWidth',3);
    title('0 = BB^T + A X E^T + E X A^T');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
end

disp('size outB.Z:');
disp(size(outB.Z));
ADI step:    1 normalized residual: 1.249224e-02 relative change in Z: 1.000000e+00
ADI step:    2 normalized residual: 1.861655e-02 relative change in Z: 2.822544e-01
ADI step:    3 normalized residual: 3.466669e-03 relative change in Z: 3.542643e-01
ADI step:    4 normalized residual: 4.613706e-03 relative change in Z: 1.895522e-01
ADI step:    5 normalized residual: 4.724451e-03 relative change in Z: 2.028230e-01
ADI step:    6 normalized residual: 5.569876e-03 relative change in Z: 2.115781e-01
ADI step:    8 normalized residual: 6.176912e-03 relative change in Z: 1.348785e-01
ADI step:    9 normalized residual: 6.588077e-03 relative change in Z: 1.139927e-01
ADI step:   10 normalized residual: 4.620909e-03 relative change in Z: 1.228504e-01
ADI step:   11 normalized residual: 1.266073e-03 relative change in Z: 1.340039e-01
ADI step:   12 normalized residual: 6.020717e-04 relative change in Z: 7.264241e-02
ADI step:   13 normalized residual: 1.132884e-04 relative change in Z: 7.928436e-02
ADI step:   15 normalized residual: 6.147827e-05 relative change in Z: 5.632046e-02
ADI step:   17 normalized residual: 4.614739e-05 relative change in Z: 3.718772e-02
ADI step:   18 normalized residual: 4.156553e-05 relative change in Z: 2.158940e-02
ADI step:   19 normalized residual: 4.028177e-05 relative change in Z: 1.325031e-02
ADI step:   20 normalized residual: 3.968866e-05 relative change in Z: 1.009396e-02
ADI step:   21 normalized residual: 4.829208e-05 relative change in Z: 1.303421e-01
ADI step:   22 normalized residual: 4.670224e-05 relative change in Z: 9.636264e-02
ADI step:   24 normalized residual: 4.113879e-05 relative change in Z: 7.980326e-02
ADI step:   26 normalized residual: 3.922578e-05 relative change in Z: 3.877317e-02
ADI step:   28 normalized residual: 3.960925e-05 relative change in Z: 1.983986e-02
ADI step:   30 normalized residual: 3.792689e-05 relative change in Z: 1.289513e-02
ADI step:   32 normalized residual: 3.631034e-05 relative change in Z: 6.538058e-02
ADI step:   34 normalized residual: 3.371415e-05 relative change in Z: 2.411855e-02
ADI step:   36 normalized residual: 3.001546e-05 relative change in Z: 1.681402e-01
ADI step:   38 normalized residual: 2.894825e-05 relative change in Z: 3.509838e-02
ADI step:   40 normalized residual: 2.310367e-05 relative change in Z: 1.367781e-01
ADI step:   41 normalized residual: 2.641129e-05 relative change in Z: 7.114943e-02
ADI step:   43 normalized residual: 2.522149e-05 relative change in Z: 2.000890e-02
ADI step:   45 normalized residual: 2.189231e-05 relative change in Z: 5.117453e-02
ADI step:   46 normalized residual: 7.579906e-06 relative change in Z: 2.780314e-01
ADI step:   48 normalized residual: 7.532878e-06 relative change in Z: 3.504571e-02
ADI step:   50 normalized residual: 6.878246e-06 relative change in Z: 1.107888e-02
ADI step:   52 normalized residual: 6.552806e-06 relative change in Z: 3.908740e-02
ADI step:   54 normalized residual: 3.750783e-06 relative change in Z: 4.953786e-02
ADI step:   56 normalized residual: 3.098839e-06 relative change in Z: 2.964697e-02
ADI step:   58 normalized residual: 2.638372e-06 relative change in Z: 6.926623e-03
ADI step:   60 normalized residual: 2.290819e-06 relative change in Z: 1.051677e-02
ADI step:   61 normalized residual: 2.639404e-06 relative change in Z: 1.029437e-02
ADI step:   62 normalized residual: 2.253806e-06 relative change in Z: 2.608155e-02
ADI step:   64 normalized residual: 2.484190e-06 relative change in Z: 8.275264e-03
ADI step:   66 normalized residual: 1.567765e-06 relative change in Z: 1.965742e-02
ADI step:   68 normalized residual: 1.426942e-06 relative change in Z: 4.424302e-03
ADI step:   70 normalized residual: 1.357729e-06 relative change in Z: 5.059327e-03
ADI step:   72 normalized residual: 1.220516e-06 relative change in Z: 6.242079e-03
ADI step:   74 normalized residual: 1.359000e-06 relative change in Z: 6.644560e-03
ADI step:   76 normalized residual: 1.035372e-07 relative change in Z: 3.050838e-02
ADI step:   78 normalized residual: 8.951783e-08 relative change in Z: 2.735113e-03
ADI step:   80 normalized residual: 7.523501e-08 relative change in Z: 2.760003e-03
ADI step:   81 normalized residual: 6.096580e-08 relative change in Z: 2.606757e-03
ADI step:   82 normalized residual: 4.348605e-08 relative change in Z: 6.827184e-03
ADI step:   84 normalized residual: 4.098094e-08 relative change in Z: 7.487300e-04
ADI step:   86 normalized residual: 1.827304e-08 relative change in Z: 4.848195e-03
ADI step:   88 normalized residual: 1.799273e-08 relative change in Z: 5.127323e-04
ADI step:   90 normalized residual: 1.695302e-08 relative change in Z: 1.522590e-03
ADI step:   92 normalized residual: 1.608422e-08 relative change in Z: 1.102288e-03
ADI step:   94 normalized residual: 1.572221e-08 relative change in Z: 7.361704e-04
ADI step:   96 normalized residual: 1.653842e-08 relative change in Z: 4.994851e-04
ADI step:   98 normalized residual: 1.649109e-08 relative change in Z: 9.588159e-04
ADI step:  100 normalized residual: 1.334391e-08 relative change in Z: 1.362960e-03
ADI step:  101 normalized residual: 1.317343e-08 relative change in Z: 1.568726e-03
ADI step:  103 normalized residual: 1.258336e-08 relative change in Z: 5.042037e-04
ADI step:  105 normalized residual: 1.238122e-08 relative change in Z: 2.133441e-04
ADI step:  107 normalized residual: 1.187104e-08 relative change in Z: 4.004058e-04
ADI step:  109 normalized residual: 2.064275e-09 relative change in Z: 2.421535e-03
ADI step:  111 normalized residual: 1.755973e-09 relative change in Z: 8.847296e-04
ADI step:  113 normalized residual: 1.707741e-09 relative change in Z: 1.263139e-04
ADI step:  115 normalized residual: 1.664406e-09 relative change in Z: 3.079688e-04
ADI step:  116 normalized residual: 2.718718e-10 relative change in Z: 1.869124e-03
ADI step:  118 normalized residual: 2.519715e-10 relative change in Z: 1.154816e-04
ADI step:  120 normalized residual: 2.341509e-10 relative change in Z: 1.979802e-04
ADI step:  121 normalized residual: 9.006480e-11 relative change in Z: 6.201689e-05
mess_lradi took   3.53 seconds 
size outB.Z:
   606   484

eqn.type='T';
t_mess_lradi = tic;
outC = mess_lradi(eqn, opts, oper);
t_elapsed2 = toc(t_mess_lradi);
fprintf(1,'mess_lradi took %6.2f seconds \n', t_elapsed2);

if istest
    if min(outC.res)>=opts.adi.res_tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result');
    end
else
    figure;
    semilogy(outC.res,'LineWidth',3);
    title('0 = C^T C + A^T X E + E^T X A');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
disp('size outC.Z:');
disp(size(outC.Z));
ADI step:    1 normalized residual: 8.341012e-01 relative change in Z: 1.000000e+00
ADI step:    2 normalized residual: 8.329581e+01 relative change in Z: 9.997990e-01
ADI step:    3 normalized residual: 3.145341e+02 relative change in Z: 9.640121e-01
ADI step:    4 normalized residual: 4.728576e+01 relative change in Z: 6.318905e-01
ADI step:    5 normalized residual: 4.332043e+01 relative change in Z: 2.389648e-01
ADI step:    6 normalized residual: 4.335982e+01 relative change in Z: 2.920793e-01
ADI step:    8 normalized residual: 5.617203e+01 relative change in Z: 4.165169e-01
ADI step:   10 normalized residual: 2.202578e+01 relative change in Z: 4.429715e-01
ADI step:   12 normalized residual: 1.489192e+01 relative change in Z: 3.616397e-01
ADI step:   13 normalized residual: 1.126917e+01 relative change in Z: 1.290182e-01
ADI step:   15 normalized residual: 1.359594e+01 relative change in Z: 1.667615e-01
ADI step:   17 normalized residual: 1.450070e+01 relative change in Z: 1.621300e-01
ADI step:   18 normalized residual: 8.266242e+00 relative change in Z: 1.081806e-01
ADI step:   19 normalized residual: 4.262993e+00 relative change in Z: 6.317378e-02
ADI step:   20 normalized residual: 4.169228e+00 relative change in Z: 9.120033e-03
ADI step:   21 normalized residual: 3.511480e+00 relative change in Z: 5.324245e-02
ADI step:   23 normalized residual: 3.385841e+00 relative change in Z: 6.771030e-02
ADI step:   24 normalized residual: 3.393270e+00 relative change in Z: 3.185809e-03
ADI step:   26 normalized residual: 3.818003e+00 relative change in Z: 1.867800e-01
ADI step:   28 normalized residual: 3.970619e+00 relative change in Z: 1.354925e-02
ADI step:   29 normalized residual: 3.749534e+00 relative change in Z: 3.378714e-02
ADI step:   31 normalized residual: 4.149615e+00 relative change in Z: 4.677622e-02
ADI step:   33 normalized residual: 2.238886e+00 relative change in Z: 1.918269e-01
ADI step:   35 normalized residual: 3.073296e+00 relative change in Z: 7.906918e-02
ADI step:   37 normalized residual: 3.763925e+00 relative change in Z: 3.402160e-02
ADI step:   39 normalized residual: 1.577345e+00 relative change in Z: 2.664323e-01
ADI step:   41 normalized residual: 1.585400e+00 relative change in Z: 2.971082e-02
ADI step:   42 normalized residual: 3.859236e+00 relative change in Z: 4.807399e-02
ADI step:   43 normalized residual: 3.860839e+00 relative change in Z: 2.856641e-03
ADI step:   45 normalized residual: 3.827653e+00 relative change in Z: 1.768378e-02
ADI step:   47 normalized residual: 2.072817e+00 relative change in Z: 1.065736e-01
ADI step:   49 normalized residual: 1.981590e+00 relative change in Z: 1.876702e-02
ADI step:   51 normalized residual: 2.062949e+00 relative change in Z: 2.454722e-02
ADI step:   53 normalized residual: 3.614905e-01 relative change in Z: 8.638273e-02
ADI step:   55 normalized residual: 3.211747e-01 relative change in Z: 2.636745e-02
ADI step:   57 normalized residual: 2.232300e-01 relative change in Z: 2.194743e-02
ADI step:   58 normalized residual: 2.160605e-01 relative change in Z: 4.054070e-03
ADI step:   59 normalized residual: 2.179565e-01 relative change in Z: 5.599015e-03
ADI step:   61 normalized residual: 1.996228e-01 relative change in Z: 9.207628e-03
ADI step:   62 normalized residual: 2.131091e-01 relative change in Z: 7.186096e-03
ADI step:   64 normalized residual: 1.962125e-01 relative change in Z: 6.852630e-03
ADI step:   66 normalized residual: 4.975793e-02 relative change in Z: 3.535560e-02
ADI step:   68 normalized residual: 5.615909e-02 relative change in Z: 4.170434e-03
ADI step:   70 normalized residual: 1.738090e-02 relative change in Z: 1.156919e-02
ADI step:   72 normalized residual: 1.608667e-02 relative change in Z: 3.952712e-03
ADI step:   74 normalized residual: 1.231934e-02 relative change in Z: 3.630013e-03
ADI step:   76 normalized residual: 1.086041e-02 relative change in Z: 5.459894e-03
ADI step:   78 normalized residual: 7.922740e-03 relative change in Z: 1.711902e-02
ADI step:   79 normalized residual: 7.564029e-03 relative change in Z: 2.300073e-03
ADI step:   80 normalized residual: 7.763314e-03 relative change in Z: 6.381721e-04
ADI step:   82 normalized residual: 5.834638e-03 relative change in Z: 6.179926e-03
ADI step:   83 normalized residual: 2.608336e-02 relative change in Z: 3.601881e-03
ADI step:   85 normalized residual: 2.643208e-02 relative change in Z: 6.472311e-03
ADI step:   86 normalized residual: 2.639864e-02 relative change in Z: 3.019126e-04
ADI step:   88 normalized residual: 2.543823e-02 relative change in Z: 1.657503e-03
ADI step:   90 normalized residual: 2.488711e-02 relative change in Z: 2.708028e-03
ADI step:   92 normalized residual: 2.989310e-03 relative change in Z: 1.208542e-02
ADI step:   94 normalized residual: 2.783350e-03 relative change in Z: 4.979101e-04
ADI step:   96 normalized residual: 2.033616e-03 relative change in Z: 1.774386e-03
ADI step:   98 normalized residual: 1.889935e-03 relative change in Z: 2.327530e-03
ADI step:   99 normalized residual: 2.322765e-03 relative change in Z: 6.304912e-04
ADI step:  101 normalized residual: 3.004050e-03 relative change in Z: 1.058910e-03
ADI step:  103 normalized residual: 1.780888e-03 relative change in Z: 1.199370e-03
ADI step:  104 normalized residual: 3.742034e-03 relative change in Z: 1.406561e-03
ADI step:  106 normalized residual: 3.795710e-03 relative change in Z: 5.532008e-04
ADI step:  107 normalized residual: 3.737064e-03 relative change in Z: 2.025597e-04
ADI step:  109 normalized residual: 3.832716e-03 relative change in Z: 6.362630e-04
ADI step:  111 normalized residual: 3.041474e-04 relative change in Z: 6.036291e-03
ADI step:  113 normalized residual: 3.045180e-04 relative change in Z: 6.337521e-04
ADI step:  115 normalized residual: 2.177559e-04 relative change in Z: 1.758688e-03
ADI step:  117 normalized residual: 2.549174e-04 relative change in Z: 3.502542e-04
ADI step:  119 normalized residual: 1.738552e-04 relative change in Z: 2.741697e-04
ADI step:  121 normalized residual: 2.539824e-04 relative change in Z: 6.609090e-04
ADI step:  123 normalized residual: 3.651742e-04 relative change in Z: 3.565011e-04
ADI step:  124 normalized residual: 1.079403e-04 relative change in Z: 3.533319e-04
ADI step:  126 normalized residual: 1.175526e-04 relative change in Z: 2.872537e-04
ADI step:  128 normalized residual: 1.547044e-04 relative change in Z: 2.069268e-04
ADI step:  130 normalized residual: 5.076777e-05 relative change in Z: 6.064107e-04
ADI step:  132 normalized residual: 6.310836e-05 relative change in Z: 1.010909e-04
ADI step:  133 normalized residual: 5.999603e-05 relative change in Z: 6.654815e-05
ADI step:  135 normalized residual: 5.070026e-05 relative change in Z: 1.301886e-04
ADI step:  137 normalized residual: 6.649690e-05 relative change in Z: 2.727154e-04
ADI step:  139 normalized residual: 8.769217e-06 relative change in Z: 1.112640e-03
ADI step:  141 normalized residual: 6.319251e-06 relative change in Z: 1.186939e-04
ADI step:  143 normalized residual: 3.193556e-06 relative change in Z: 5.811154e-05
ADI step:  144 normalized residual: 1.313842e-06 relative change in Z: 3.654803e-05
ADI step:  145 normalized residual: 1.296647e-06 relative change in Z: 8.065006e-06
ADI step:  147 normalized residual: 1.258418e-06 relative change in Z: 1.391970e-05
ADI step:  149 normalized residual: 1.121113e-06 relative change in Z: 4.110717e-05
ADI step:  151 normalized residual: 1.065404e-06 relative change in Z: 3.147982e-05
ADI step:  152 normalized residual: 1.098363e-06 relative change in Z: 7.989842e-06
ADI step:  154 normalized residual: 8.858622e-08 relative change in Z: 7.968126e-05
ADI step:  156 normalized residual: 8.749477e-08 relative change in Z: 1.078510e-05
ADI step:  158 normalized residual: 1.774616e-07 relative change in Z: 1.775152e-05
ADI step:  160 normalized residual: 5.637692e-08 relative change in Z: 2.683975e-05
ADI step:  162 normalized residual: 1.784761e-08 relative change in Z: 1.892640e-05
ADI step:  164 normalized residual: 2.090382e-08 relative change in Z: 1.507129e-06
ADI step:  165 normalized residual: 2.739810e-08 relative change in Z: 7.017073e-06
ADI step:  167 normalized residual: 2.689303e-08 relative change in Z: 9.712160e-06
ADI step:  168 normalized residual: 2.385078e-08 relative change in Z: 1.047469e-06
ADI step:  169 normalized residual: 2.433743e-08 relative change in Z: 5.897054e-07
ADI step:  171 normalized residual: 7.880353e-09 relative change in Z: 4.839816e-06
ADI step:  173 normalized residual: 8.932618e-09 relative change in Z: 4.235514e-06
ADI step:  175 normalized residual: 2.948845e-09 relative change in Z: 7.629617e-06
ADI step:  177 normalized residual: 2.558604e-09 relative change in Z: 1.106287e-06
ADI step:  179 normalized residual: 2.775637e-09 relative change in Z: 5.514441e-07
ADI step:  181 normalized residual: 1.927496e-09 relative change in Z: 3.106171e-06
ADI step:  183 normalized residual: 3.694518e-09 relative change in Z: 1.737800e-06
ADI step:  185 normalized residual: 3.439907e-10 relative change in Z: 3.700856e-06
ADI step:  186 normalized residual: 2.450340e-10 relative change in Z: 3.674148e-07
ADI step:  188 normalized residual: 2.084278e-10 relative change in Z: 3.792256e-07
ADI step:  190 normalized residual: 2.012881e-10 relative change in Z: 1.104120e-07
ADI step:  192 normalized residual: 1.825901e-10 relative change in Z: 2.768216e-07
ADI step:  194 normalized residual: 2.801960e-10 relative change in Z: 7.520497e-07
ADI step:  196 normalized residual: 1.778026e-10 relative change in Z: 2.609004e-07
ADI step:  197 normalized residual: 1.925840e-10 relative change in Z: 7.852756e-08
ADI step:  199 normalized residual: 2.255852e-11 relative change in Z: 1.013337e-06
mess_lradi took   5.27 seconds 
size outC.Z:
   606   796

Compute reduced system matrices

Perform Square Root Method (SRM)

% BT tolerance and maximum order for the ROM
opts.srm.tol=1e-3;
opts.srm.max_ord=250;

% SRM verbosity
if istest
    opts.srm.info=1;
else
    opts.srm.info=2;
end

% The actual SRM
[TL,TR,hsv] = mess_square_root_method(eqn,opts,oper,outB.Z,outC.Z);

% compute ROM matrices
B1 = TL'*(eqn.A_(1:eqn.st,1:eqn.st))*TR;
B2 = TL'*(eqn.A_(1:eqn.st,eqn.st+1:end));
A1 = eqn.A_(eqn.st+1:end,1:eqn.st)*TR;

ROM.A = B1 - B2*(eqn.A_(eqn.st+1:end,eqn.st+1:end)\A1);
ROM.B = TL'*eqn.B(1:eqn.st,:) - ...
     B2*(eqn.A_(eqn.st+1:end,eqn.st+1:end)\eqn.B(eqn.st+1:end,:));
ROM.C = eqn.C(:,1:eqn.st)*TR - ...
     eqn.C(:,eqn.st+1:end)*(eqn.A_(eqn.st+1:end,eqn.st+1:end)\A1);
ROM.D = -eqn.C(:,eqn.st+1:end)*(eqn.A_(eqn.st+1:end,eqn.st+1:end)\...
    eqn.B(eqn.st+1:end,:));
ROM.E = eye(size(ROM.A));
reduced system order: 96  (max possible/allowed: 484/250)

Evaluate the ROM quality

while the Gramians are computed on the hidden manifold, we need to do the frequency domain computations without (implicitly) using the Schur complement (due to the construction of the function handles)

oper = operatormanager('default');

if istest
    opts.sigma.info=0;
else
    opts.sigma.info=2;
end

opts.sigma.fmin=-3;
opts.sigma.fmax=4;

out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err;

if istest
    if max(err)>5e-3
        error('MESS:TEST:accuracy','unexpectedly inaccurate result');
    end
else
    figure;
    semilogy(hsv,'LineWidth',3);
    title('Computed Hankel singular values');
    xlabel('index');
    ylabel('magnitude');
end
Computing TFMs of original and reduced order systems and MOR errors
 Step  10 / 100 Step  20 / 100 Step  30 / 100 Step  40 / 100 Step  50 / 100 Step  60 / 100 Step  70 / 100 Step  80 / 100 Step  90 / 100 Step 100 / 100

reset warning state

warning(orig_warnstate);