function BT_sym_TripleChain(version,istest)
% % Computes a reduced order model (ROM) for the triple chain example of % Truhar and Veselic [1] via Balanced truncation, e.g. [2], exploiting % symmetry of the control system, i.e. equality of the system Gramians. % In comparison to BT_TripleChain this demo exploits, that the first order % form of the system has symmetric E and A and B=C', such that the two % Lyapunov equations coincide. % % Usage: BT_sym_TripleChain(version) % % Input: % % version Decides the Balanced Truncation version to use. % Possible values: % 'FO' for reduction of the first order form to first order form % 'VV' velocity-velocity balancing of the second order form to % second order form. % 'PP' position-position balancing of the second order form to % second order form. % 'PV' position-velocity balancing of the second order form to % second order form. % 'VP' velocity-position balancing of the second order form to % second order form.% % % istest flag to determine whether this demo runs as a CI test or % interactive demo % (optional, defaults to 0, i.e. interactive demo) % % References: % % [1] 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. % % [2] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems, Vol. % 6 of Adv. Des. Control, SIAM Publications, Philadelphia, PA, 2005. % https://doi.org/10.1137/1.9780898718713. % % 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 %
narginchk(0,2); if nargin==0 version = 'FO'; end if nargin<2 istest=0; end format long e; % set operation oper = operatormanager('so_2'); % Problem data n1=500; alpha=.002; beta=alpha; v=5; [eqn.M_,eqn.E_,eqn.K_]=triplechain_MSD(n1,alpha,beta,v); s = size(eqn.K_,1); B = ones(s,1); O = zeros(s,1); Cp = ones(1,s); Cv = O'; eqn.B = [B; O]; eqn.C = [Cp, Cv]; eqn.haveE=1;
ADI tolerances and maximum iteration number
opts.adi.maxiter = 300; opts.adi.res_tol = 1e-10; opts.norm = 'fro'; %opts.adi.rel_diff_tol = 1e-16; opts.adi.rel_diff_tol = 0; opts.adi.info = 1; %opts.adi.accumulateK = 0; %opts.adi.accumulateDeltaK = 0; opts.adi.compute_sol_fac = 1; eqn.type = 'N';
%Heuristic shift parameters via projection opts.shifts.num_desired=5; opts.shifts.info=0; opts.shifts.method = 'projection';
Compute Gramian Factor (one is enough, since the Gramians are equal)
tic; outB = mess_lradi(eqn, opts, oper); toc; if istest if min(outB.res)>=1e-1 error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(1); semilogy(outB.res); title('0= AXM^T + MXA^T -BB^T'); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end disp('size outB.Z:'); disp(size(outB.Z));
ADI step: 1 normalized residual: 1.270108e+00 ADI step: 2 normalized residual: 3.427418e+00 ADI step: 4 normalized residual: 4.477922e+02 ADI step: 6 normalized residual: 3.943500e+02 ADI step: 8 normalized residual: 6.085237e+02 ADI step: 9 normalized residual: 7.293151e+03 ADI step: 11 normalized residual: 7.332737e+03 ADI step: 13 normalized residual: 7.483298e+03 ADI step: 14 normalized residual: 5.877023e+02 ADI step: 16 normalized residual: 7.272517e+02 ADI step: 18 normalized residual: 8.594421e+02 ADI step: 20 normalized residual: 2.043734e+02 ADI step: 21 normalized residual: 8.008062e+02 ADI step: 23 normalized residual: 7.946889e+02 ADI step: 25 normalized residual: 6.376582e+01 ADI step: 27 normalized residual: 6.430936e+01 ADI step: 28 normalized residual: 9.394251e+01 ADI step: 29 normalized residual: 7.753989e+01 ADI step: 31 normalized residual: 8.002104e+00 ADI step: 33 normalized residual: 7.641858e+00 ADI step: 35 normalized residual: 6.705939e+00 ADI step: 36 normalized residual: 5.442093e+00 ADI step: 38 normalized residual: 5.539283e+00 ADI step: 40 normalized residual: 4.767961e+00 ADI step: 42 normalized residual: 4.926510e+00 ADI step: 43 normalized residual: 5.954219e+00 ADI step: 44 normalized residual: 4.564030e+00 ADI step: 46 normalized residual: 5.300509e+00 ADI step: 48 normalized residual: 2.859722e+00 ADI step: 50 normalized residual: 1.514706e+00 ADI step: 51 normalized residual: 3.250111e+00 ADI step: 53 normalized residual: 2.174333e+00 ADI step: 55 normalized residual: 1.813607e+00 ADI step: 57 normalized residual: 1.797216e+00 ADI step: 58 normalized residual: 1.571941e+00 ADI step: 60 normalized residual: 1.090202e+00 ADI step: 62 normalized residual: 1.067706e+00 ADI step: 63 normalized residual: 1.210877e+00 ADI step: 65 normalized residual: 1.156980e+00 ADI step: 66 normalized residual: 7.188757e-01 ADI step: 68 normalized residual: 1.172664e+00 ADI step: 70 normalized residual: 1.129371e+00 ADI step: 71 normalized residual: 8.961108e-01 ADI step: 73 normalized residual: 8.713580e-01 ADI step: 75 normalized residual: 7.378370e-01 ADI step: 77 normalized residual: 5.279300e-01 ADI step: 78 normalized residual: 8.063868e-01 ADI step: 80 normalized residual: 6.666620e-01 ADI step: 81 normalized residual: 7.099150e-01 ADI step: 83 normalized residual: 6.526702e-01 ADI step: 85 normalized residual: 4.615182e-01 ADI step: 87 normalized residual: 4.450105e-01 ADI step: 88 normalized residual: 5.715083e-01 ADI step: 90 normalized residual: 5.039076e-01 ADI step: 92 normalized residual: 3.200358e-01 ADI step: 93 normalized residual: 4.142990e-01 ADI step: 95 normalized residual: 4.217200e-01 ADI step: 97 normalized residual: 3.333007e-01 ADI step: 98 normalized residual: 3.777807e-01 ADI step: 100 normalized residual: 3.617455e-01 ADI step: 101 normalized residual: 3.912320e-01 ADI step: 103 normalized residual: 3.370354e-01 ADI step: 105 normalized residual: 3.317509e-01 ADI step: 107 normalized residual: 3.355219e-01 ADI step: 108 normalized residual: 3.173621e-01 ADI step: 110 normalized residual: 2.591779e-01 ADI step: 111 normalized residual: 2.713245e-01 ADI step: 113 normalized residual: 2.680856e-01 ADI step: 115 normalized residual: 2.582790e-01 ADI step: 117 normalized residual: 2.575180e-01 ADI step: 118 normalized residual: 1.824067e-01 ADI step: 120 normalized residual: 2.048517e-01 ADI step: 121 normalized residual: 2.396445e-01 ADI step: 123 normalized residual: 2.420309e-01 ADI step: 124 normalized residual: 2.250282e-01 ADI step: 126 normalized residual: 2.267304e-01 ADI step: 128 normalized residual: 2.233017e-01 ADI step: 130 normalized residual: 2.275887e-01 ADI step: 132 normalized residual: 2.323182e-01 ADI step: 133 normalized residual: 2.233032e-01 ADI step: 135 normalized residual: 1.474431e-01 ADI step: 137 normalized residual: 1.422854e-01 ADI step: 138 normalized residual: 1.711563e-01 ADI step: 139 normalized residual: 1.590640e-01 ADI step: 141 normalized residual: 1.009243e-01 ADI step: 143 normalized residual: 9.804670e-02 ADI step: 144 normalized residual: 1.874438e-01 ADI step: 146 normalized residual: 1.877888e-01 ADI step: 148 normalized residual: 1.842147e-01 ADI step: 149 normalized residual: 2.021664e-01 ADI step: 151 normalized residual: 2.020127e-01 ADI step: 153 normalized residual: 1.791933e-01 ADI step: 155 normalized residual: 1.766039e-01 ADI step: 157 normalized residual: 1.592584e-01 ADI step: 158 normalized residual: 1.749497e-01 ADI step: 160 normalized residual: 1.720773e-01 ADI step: 162 normalized residual: 1.675864e-01 ADI step: 163 normalized residual: 1.559620e-01 ADI step: 165 normalized residual: 1.566123e-01 ADI step: 167 normalized residual: 1.591720e-01 ADI step: 168 normalized residual: 1.588365e-01 ADI step: 169 normalized residual: 1.641616e-01 ADI step: 171 normalized residual: 1.639970e-01 ADI step: 173 normalized residual: 1.451533e-01 ADI step: 175 normalized residual: 1.383603e-01 ADI step: 177 normalized residual: 1.234213e-01 ADI step: 178 normalized residual: 1.053449e-01 ADI step: 179 normalized residual: 1.262741e-01 ADI step: 181 normalized residual: 1.191211e-01 ADI step: 183 normalized residual: 1.181099e-01 ADI step: 184 normalized residual: 6.863730e-02 ADI step: 186 normalized residual: 6.824738e-02 ADI step: 188 normalized residual: 6.207974e-02 ADI step: 189 normalized residual: 9.478628e-02 ADI step: 191 normalized residual: 9.501622e-02 ADI step: 193 normalized residual: 9.264676e-02 ADI step: 194 normalized residual: 9.925469e-02 ADI step: 196 normalized residual: 9.861238e-02 ADI step: 198 normalized residual: 9.929500e-02 ADI step: 199 normalized residual: 1.043346e-01 ADI step: 201 normalized residual: 1.040538e-01 ADI step: 203 normalized residual: 1.035106e-01 ADI step: 204 normalized residual: 7.429947e-02 ADI step: 206 normalized residual: 7.597314e-02 ADI step: 208 normalized residual: 7.518157e-02 ADI step: 209 normalized residual: 9.231182e-02 ADI step: 211 normalized residual: 8.773386e-02 ADI step: 213 normalized residual: 8.725790e-02 ADI step: 214 normalized residual: 6.214829e-02 ADI step: 216 normalized residual: 6.135600e-02 ADI step: 218 normalized residual: 5.987776e-02 ADI step: 219 normalized residual: 1.004795e-01 ADI step: 221 normalized residual: 1.003918e-01 ADI step: 223 normalized residual: 9.785594e-02 ADI step: 224 normalized residual: 9.320641e-02 ADI step: 226 normalized residual: 9.283261e-02 ADI step: 228 normalized residual: 8.572871e-02 ADI step: 229 normalized residual: 9.553589e-02 ADI step: 231 normalized residual: 9.657063e-02 ADI step: 233 normalized residual: 9.249383e-02 ADI step: 234 normalized residual: 7.123043e-02 ADI step: 236 normalized residual: 6.958544e-02 ADI step: 238 normalized residual: 7.414030e-02 ADI step: 239 normalized residual: 5.894823e-02 ADI step: 241 normalized residual: 6.038939e-02 ADI step: 243 normalized residual: 5.916657e-02 ADI step: 244 normalized residual: 8.411748e-02 ADI step: 246 normalized residual: 8.425901e-02 ADI step: 248 normalized residual: 8.243929e-02 ADI step: 249 normalized residual: 7.637549e-02 ADI step: 251 normalized residual: 7.632137e-02 ADI step: 253 normalized residual: 7.705582e-02 ADI step: 255 normalized residual: 7.636337e-02 ADI step: 256 normalized residual: 7.639694e-02 ADI step: 258 normalized residual: 7.923463e-02 ADI step: 259 normalized residual: 4.060082e-02 ADI step: 261 normalized residual: 4.059216e-02 ADI step: 263 normalized residual: 3.871716e-02 ADI step: 264 normalized residual: 4.788324e-02 ADI step: 266 normalized residual: 4.702063e-02 ADI step: 268 normalized residual: 4.721126e-02 ADI step: 269 normalized residual: 5.602991e-02 ADI step: 271 normalized residual: 5.538506e-02 ADI step: 273 normalized residual: 5.433885e-02 ADI step: 274 normalized residual: 3.250868e-02 ADI step: 276 normalized residual: 3.236078e-02 ADI step: 278 normalized residual: 3.151680e-02 ADI step: 279 normalized residual: 4.724736e-02 ADI step: 281 normalized residual: 4.753528e-02 ADI step: 283 normalized residual: 4.491330e-02 ADI step: 284 normalized residual: 4.464692e-02 ADI step: 286 normalized residual: 4.490678e-02 ADI step: 287 normalized residual: 5.450309e-02 ADI step: 288 normalized residual: 5.716043e-02 ADI step: 289 normalized residual: 7.086736e-02 ADI step: 291 normalized residual: 7.097008e-02 ADI step: 293 normalized residual: 6.753440e-02 ADI step: 294 normalized residual: 4.793382e-02 ADI step: 296 normalized residual: 4.775572e-02 ADI step: 298 normalized residual: 4.622585e-02 ADI step: 299 normalized residual: 5.323988e-02 ADI step: 301 normalized residual: 5.316499e-02 Elapsed time is 0.364077 seconds. size outB.Z: 3002 301

switch upper(version) case 'FO'
Compute first order ROM
opts.srm.max_ord = 150; opts.srm.tol = eps; opts.srm.info = 1; [TL,TR] = mess_square_root_method(eqn,opts,oper,outB.Z,outB.Z); ROM.E = eye(size(TL,2)); ROM.A = TL'*oper.mul_A(eqn, opts, 'N', TR, 'N'); ROM.B = TL'*eqn.B; ROM.C = eqn.C*TR; ROM.D = [];
reduced system order: 150 (max possible/allowed: 301/150)
case 'VV' U = outB.Z(1:s,:); V = outB.Z(1:s,:); case 'PP' U = outB.Z(s+1:end,:); V = outB.Z(s+1:end,:); case 'PV' U = outB.Z(s+1:end,:); V = outB.Z(1:s,:); case 'VP' U = outB.Z(1:s,:); V = outB.Z(s+1:end,:); end if not(strcmp(version,'FO')) max_ord = 75; tol = eps; info = 1; [TL,TR] = square_root_method_SO(eqn.M_, max_ord, tol, info, U, V); ROM.M = eye(size(TL,2)); ROM.E = TL'*(eqn.E_*TR); ROM.K = TL'*(eqn.K_*TR); ROM.B = TL'*B; ROM.Cv = Cv*TR; ROM.Cp = Cp*TR; end
plot results
opts.sigma.fmin = 1e-4; opts.sigma.fmax = 1e0; opts.sigma.nsample = 200; if istest opts.sigma.info = 1; else opts.sigma.info = 2; end out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err; if istest if max(err) > 1000 error('MESS:TEST:accuracy','unexpectedly inaccurate result %g', max(err)); end end
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

