Contents
function bt_mor_DAE1_tol(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: % 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 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 %
if nargin<1, istest=0; end
set operation manager for the Gramian computations
oper = operatormanager('dae_1');
Read problem data
fname = sprintf('%s/../models/BIPS/bips98_606.mat',... fileparts(mfilename('fullpath'))); Bips = load(fname); % from https://sites.google.com/site/rommes/software % Note that we here use the alpha shift suggested by Rommes and coauthors % [1], but do not perform the shift back. p = find(diag(Bips.E)); np = find(diag(Bips.E) == 0); pp = [p;np]; eqn.A_ = Bips.A(pp, pp)-0.05*Bips.E(pp,pp); eqn.E_ = Bips.E(pp, pp); eqn.B = Bips.b(pp, :); eqn.C = Bips.c( : , pp); eqn.st = length(p); eqn.haveE = 1;
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;
opts.shifts.p=mess_para(eqn,opts,oper);
disp(opts.shifts.p);
-332.4833 -198.8877 -66.7167 -1.3436
eqn.type='N'; tic; outB = mess_lradi(eqn, opts, oper); toc; if istest if min(outB.res)>=opts.adi.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure; semilogy(outB.res); title('0= BB^T + AXM^T + MXA^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.829220e-05 relative change in Z: 1.303414e-01 ADI step: 22 normalized residual: 4.671708e-05 relative change in Z: 9.591799e-02 ADI step: 24 normalized residual: 4.115405e-05 relative change in Z: 7.980469e-02 ADI step: 26 normalized residual: 3.924151e-05 relative change in Z: 3.876808e-02 ADI step: 28 normalized residual: 3.763448e-05 relative change in Z: 1.287818e-02 ADI step: 30 normalized residual: 3.794276e-05 relative change in Z: 1.986166e-02 ADI step: 32 normalized residual: 3.632613e-05 relative change in Z: 6.538409e-02 ADI step: 34 normalized residual: 3.372976e-05 relative change in Z: 2.412149e-02 ADI step: 36 normalized residual: 3.002802e-05 relative change in Z: 1.682046e-01 ADI step: 38 normalized residual: 2.896024e-05 relative change in Z: 3.510276e-02 ADI step: 40 normalized residual: 2.311519e-05 relative change in Z: 1.367862e-01 ADI step: 41 normalized residual: 2.642194e-05 relative change in Z: 7.117977e-02 ADI step: 43 normalized residual: 2.523100e-05 relative change in Z: 2.002000e-02 ADI step: 45 normalized residual: 2.190225e-05 relative change in Z: 5.117811e-02 ADI step: 46 normalized residual: 7.579046e-06 relative change in Z: 2.781234e-01 ADI step: 48 normalized residual: 7.532084e-06 relative change in Z: 3.504476e-02 ADI step: 50 normalized residual: 6.877403e-06 relative change in Z: 1.107903e-02 ADI step: 52 normalized residual: 6.551939e-06 relative change in Z: 3.908824e-02 ADI step: 54 normalized residual: 3.750068e-06 relative change in Z: 4.953681e-02 ADI step: 56 normalized residual: 3.098153e-06 relative change in Z: 2.964644e-02 ADI step: 58 normalized residual: 2.637938e-06 relative change in Z: 6.923949e-03 ADI step: 60 normalized residual: 2.290169e-06 relative change in Z: 1.051753e-02 ADI step: 61 normalized residual: 2.640161e-06 relative change in Z: 1.029665e-02 ADI step: 62 normalized residual: 2.254367e-06 relative change in Z: 2.609252e-02 ADI step: 64 normalized residual: 2.484287e-06 relative change in Z: 8.268925e-03 ADI step: 66 normalized residual: 1.567875e-06 relative change in Z: 1.965620e-02 ADI step: 68 normalized residual: 1.427024e-06 relative change in Z: 4.423741e-03 ADI step: 70 normalized residual: 1.357797e-06 relative change in Z: 5.060376e-03 ADI step: 72 normalized residual: 1.220496e-06 relative change in Z: 6.241354e-03 ADI step: 74 normalized residual: 1.359185e-06 relative change in Z: 6.646571e-03 ADI step: 76 normalized residual: 1.035576e-07 relative change in Z: 3.050473e-02 ADI step: 78 normalized residual: 8.954692e-08 relative change in Z: 2.735657e-03 ADI step: 80 normalized residual: 7.524479e-08 relative change in Z: 2.760664e-03 ADI step: 81 normalized residual: 6.097447e-08 relative change in Z: 2.607093e-03 ADI step: 82 normalized residual: 4.347984e-08 relative change in Z: 6.829676e-03 ADI step: 84 normalized residual: 4.097442e-08 relative change in Z: 7.486790e-04 ADI step: 86 normalized residual: 1.827831e-08 relative change in Z: 4.847401e-03 ADI step: 88 normalized residual: 1.799834e-08 relative change in Z: 5.127881e-04 ADI step: 90 normalized residual: 1.695841e-08 relative change in Z: 1.522774e-03 ADI step: 92 normalized residual: 1.608588e-08 relative change in Z: 1.103021e-03 ADI step: 94 normalized residual: 1.572546e-08 relative change in Z: 7.364181e-04 ADI step: 96 normalized residual: 1.653972e-08 relative change in Z: 4.994718e-04 ADI step: 98 normalized residual: 1.649536e-08 relative change in Z: 9.589022e-04 ADI step: 100 normalized residual: 1.334355e-08 relative change in Z: 1.363124e-03 ADI step: 101 normalized residual: 1.317651e-08 relative change in Z: 1.567924e-03 ADI step: 103 normalized residual: 1.264063e-08 relative change in Z: 4.818100e-04 ADI step: 105 normalized residual: 1.243398e-08 relative change in Z: 2.174048e-04 ADI step: 107 normalized residual: 2.035920e-09 relative change in Z: 2.454738e-03 ADI step: 109 normalized residual: 2.064934e-09 relative change in Z: 1.394831e-04 ADI step: 111 normalized residual: 2.016980e-09 relative change in Z: 1.353135e-04 ADI step: 113 normalized residual: 1.597638e-09 relative change in Z: 1.029640e-03 ADI step: 115 normalized residual: 1.555359e-09 relative change in Z: 3.025770e-04 ADI step: 116 normalized residual: 3.268467e-10 relative change in Z: 1.752376e-03 ADI step: 118 normalized residual: 3.074456e-10 relative change in Z: 1.186513e-04 ADI step: 120 normalized residual: 2.790767e-10 relative change in Z: 2.570926e-04 ADI step: 121 normalized residual: 1.383010e-10 relative change in Z: 7.258464e-05 ADI step: 123 normalized residual: 1.120660e-11 relative change in Z: 5.737296e-04 Elapsed time is 2.687318 seconds. size outB.Z: 606 492

eqn.type='T'; tic; outC = mess_lradi(eqn, opts, oper); toc; if istest if min(outC.res)>=opts.adi.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure; semilogy(outC.res); title('0= C^TC + A^TXE + E^TXA'); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end disp('size outC.Z:'); disp(size(outC.Z));
ADI step: 1 normalized residual: 1.202824e+00 relative change in Z: 1.000000e+00 ADI step: 2 normalized residual: 4.050102e+00 relative change in Z: 8.985337e-01 ADI step: 3 normalized residual: 2.305832e+01 relative change in Z: 9.605920e-01 ADI step: 4 normalized residual: 5.601395e+01 relative change in Z: 9.938077e-01 ADI step: 5 normalized residual: 1.877752e+02 relative change in Z: 5.296185e-01 ADI step: 7 normalized residual: 2.232009e+02 relative change in Z: 3.184176e-01 ADI step: 9 normalized residual: 1.957431e+02 relative change in Z: 4.332124e-01 ADI step: 11 normalized residual: 4.320781e+01 relative change in Z: 6.347510e-01 ADI step: 13 normalized residual: 4.222705e+01 relative change in Z: 7.523328e-02 ADI step: 15 normalized residual: 1.594059e+01 relative change in Z: 2.391749e-01 ADI step: 16 normalized residual: 1.799811e+01 relative change in Z: 1.509212e-01 ADI step: 17 normalized residual: 1.686195e+01 relative change in Z: 9.987618e-02 ADI step: 18 normalized residual: 1.466933e+01 relative change in Z: 7.872198e-02 ADI step: 19 normalized residual: 1.355524e+01 relative change in Z: 5.282422e-02 ADI step: 20 normalized residual: 1.346740e+01 relative change in Z: 7.777360e-03 ADI step: 21 normalized residual: 9.582863e+00 relative change in Z: 1.172079e-01 ADI step: 22 normalized residual: 9.582422e+00 relative change in Z: 6.048739e-03 ADI step: 24 normalized residual: 8.738475e+00 relative change in Z: 9.509954e-02 ADI step: 25 normalized residual: 9.009804e+00 relative change in Z: 3.304692e-02 ADI step: 27 normalized residual: 5.350630e+00 relative change in Z: 3.667354e-01 ADI step: 29 normalized residual: 4.928253e+00 relative change in Z: 7.713453e-02 ADI step: 31 normalized residual: 5.083533e+00 relative change in Z: 1.772831e-01 ADI step: 33 normalized residual: 5.223994e+00 relative change in Z: 1.916352e-02 ADI step: 35 normalized residual: 4.795910e+00 relative change in Z: 2.148451e-01 ADI step: 37 normalized residual: 3.612973e+00 relative change in Z: 8.987243e-02 ADI step: 39 normalized residual: 2.880429e+00 relative change in Z: 3.577668e-02 ADI step: 41 normalized residual: 4.065443e+00 relative change in Z: 3.772076e-02 ADI step: 42 normalized residual: 3.759057e+00 relative change in Z: 2.395453e-02 ADI step: 44 normalized residual: 1.745338e+00 relative change in Z: 1.033159e-01 ADI step: 46 normalized residual: 1.816830e+00 relative change in Z: 1.790056e-02 ADI step: 48 normalized residual: 1.595539e+00 relative change in Z: 3.174118e-02 ADI step: 50 normalized residual: 1.764331e+00 relative change in Z: 1.729769e-02 ADI step: 51 normalized residual: 1.760178e+00 relative change in Z: 2.012452e-03 ADI step: 53 normalized residual: 2.048000e+00 relative change in Z: 4.492451e-02 ADI step: 55 normalized residual: 1.720457e+00 relative change in Z: 3.451409e-02 ADI step: 57 normalized residual: 1.933750e+00 relative change in Z: 2.474784e-02 ADI step: 59 normalized residual: 8.012666e-01 relative change in Z: 7.041308e-02 ADI step: 61 normalized residual: 5.471749e-01 relative change in Z: 1.411295e-02 ADI step: 62 normalized residual: 1.069122e+00 relative change in Z: 1.651164e-02 ADI step: 64 normalized residual: 8.730749e-01 relative change in Z: 2.104850e-02 ADI step: 66 normalized residual: 1.041290e+00 relative change in Z: 2.881269e-02 ADI step: 68 normalized residual: 1.101326e+00 relative change in Z: 8.375999e-03 ADI step: 70 normalized residual: 5.379291e-02 relative change in Z: 7.054036e-02 ADI step: 71 normalized residual: 5.479420e-02 relative change in Z: 1.304036e-03 ADI step: 73 normalized residual: 3.150090e-02 relative change in Z: 8.460189e-03 ADI step: 75 normalized residual: 3.119396e-02 relative change in Z: 3.336979e-03 ADI step: 77 normalized residual: 3.705368e-02 relative change in Z: 5.943746e-03 ADI step: 79 normalized residual: 1.659676e-02 relative change in Z: 1.145817e-02 ADI step: 81 normalized residual: 3.436526e-02 relative change in Z: 6.480687e-03 ADI step: 82 normalized residual: 5.716515e-02 relative change in Z: 4.231095e-03 ADI step: 84 normalized residual: 1.299322e-02 relative change in Z: 1.901395e-02 ADI step: 86 normalized residual: 7.385314e-03 relative change in Z: 2.882285e-03 ADI step: 88 normalized residual: 1.061139e-02 relative change in Z: 2.874863e-03 ADI step: 90 normalized residual: 5.297361e-03 relative change in Z: 5.596457e-03 ADI step: 91 normalized residual: 5.356198e-03 relative change in Z: 4.258366e-04 ADI step: 93 normalized residual: 6.844082e-03 relative change in Z: 4.372806e-03 ADI step: 95 normalized residual: 7.357109e-03 relative change in Z: 6.812892e-04 ADI step: 97 normalized residual: 4.883417e-03 relative change in Z: 3.660378e-03 ADI step: 99 normalized residual: 9.244692e-04 relative change in Z: 2.520909e-03 ADI step: 101 normalized residual: 2.123864e-03 relative change in Z: 1.403914e-03 ADI step: 102 normalized residual: 4.907279e-03 relative change in Z: 1.304758e-03 ADI step: 104 normalized residual: 3.164948e-03 relative change in Z: 1.854023e-03 ADI step: 106 normalized residual: 3.946063e-03 relative change in Z: 2.351810e-03 ADI step: 108 normalized residual: 3.009306e-03 relative change in Z: 1.192170e-03 ADI step: 109 normalized residual: 3.067043e-03 relative change in Z: 1.601425e-04 ADI step: 111 normalized residual: 3.189752e-03 relative change in Z: 2.317468e-04 ADI step: 113 normalized residual: 3.324455e-04 relative change in Z: 3.747250e-03 ADI step: 115 normalized residual: 3.216530e-04 relative change in Z: 5.056479e-04 ADI step: 117 normalized residual: 1.471250e-04 relative change in Z: 5.127442e-04 ADI step: 119 normalized residual: 3.339876e-04 relative change in Z: 7.928982e-04 ADI step: 121 normalized residual: 4.484060e-04 relative change in Z: 1.374526e-03 ADI step: 122 normalized residual: 3.337753e-04 relative change in Z: 1.828815e-04 ADI step: 124 normalized residual: 3.758983e-04 relative change in Z: 2.605702e-04 ADI step: 126 normalized residual: 2.283871e-04 relative change in Z: 4.771613e-04 ADI step: 128 normalized residual: 2.849426e-04 relative change in Z: 5.813636e-04 ADI step: 130 normalized residual: 2.746930e-04 relative change in Z: 5.280990e-05 ADI step: 132 normalized residual: 7.194683e-06 relative change in Z: 1.147138e-03 ADI step: 134 normalized residual: 4.515688e-06 relative change in Z: 6.038787e-05 ADI step: 135 normalized residual: 4.567368e-06 relative change in Z: 6.436784e-06 ADI step: 137 normalized residual: 1.187164e-06 relative change in Z: 7.891308e-05 ADI step: 139 normalized residual: 6.274526e-07 relative change in Z: 5.774421e-05 ADI step: 141 normalized residual: 5.229292e-07 relative change in Z: 1.491434e-05 ADI step: 142 normalized residual: 4.260937e-07 relative change in Z: 8.344926e-06 ADI step: 144 normalized residual: 4.256090e-07 relative change in Z: 6.206398e-07 ADI step: 146 normalized residual: 3.262073e-07 relative change in Z: 1.697903e-05 ADI step: 148 normalized residual: 4.159036e-07 relative change in Z: 1.765883e-05 ADI step: 150 normalized residual: 3.656251e-07 relative change in Z: 8.175008e-06 ADI step: 152 normalized residual: 8.496967e-08 relative change in Z: 2.028409e-05 ADI step: 153 normalized residual: 8.367522e-08 relative change in Z: 1.347172e-06 ADI step: 155 normalized residual: 5.117775e-08 relative change in Z: 1.202413e-05 ADI step: 157 normalized residual: 1.901574e-07 relative change in Z: 8.645037e-06 ADI step: 159 normalized residual: 1.123623e-07 relative change in Z: 8.918116e-06 ADI step: 161 normalized residual: 2.179426e-08 relative change in Z: 6.309749e-06 ADI step: 162 normalized residual: 1.032776e-07 relative change in Z: 5.836519e-06 ADI step: 164 normalized residual: 9.696707e-08 relative change in Z: 2.576380e-06 ADI step: 165 normalized residual: 9.619087e-08 relative change in Z: 6.244908e-07 ADI step: 167 normalized residual: 8.844280e-08 relative change in Z: 4.475350e-06 ADI step: 169 normalized residual: 6.090795e-08 relative change in Z: 4.698606e-06 ADI step: 171 normalized residual: 5.319742e-08 relative change in Z: 1.811500e-06 ADI step: 173 normalized residual: 4.491929e-08 relative change in Z: 4.997787e-06 ADI step: 175 normalized residual: 4.379374e-08 relative change in Z: 1.693454e-06 ADI step: 177 normalized residual: 3.720800e-08 relative change in Z: 2.667125e-06 ADI step: 179 normalized residual: 3.580057e-08 relative change in Z: 4.361012e-07 ADI step: 181 normalized residual: 3.337012e-09 relative change in Z: 5.901256e-06 ADI step: 182 normalized residual: 2.840442e-09 relative change in Z: 4.957047e-07 ADI step: 184 normalized residual: 2.719953e-09 relative change in Z: 6.021206e-07 ADI step: 186 normalized residual: 2.955935e-09 relative change in Z: 5.144543e-07 ADI step: 188 normalized residual: 2.353533e-09 relative change in Z: 1.293506e-06 ADI step: 190 normalized residual: 2.102016e-09 relative change in Z: 1.139851e-06 ADI step: 192 normalized residual: 2.080104e-09 relative change in Z: 1.002696e-07 ADI step: 194 normalized residual: 7.179784e-10 relative change in Z: 1.160817e-06 ADI step: 196 normalized residual: 8.789659e-10 relative change in Z: 6.508999e-07 ADI step: 198 normalized residual: 6.898603e-10 relative change in Z: 2.174102e-07 ADI step: 200 normalized residual: 5.120167e-11 relative change in Z: 9.074723e-07 Elapsed time is 4.310821 seconds. size outC.Z: 606 800

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: 492/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); 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);