function BT_sym_TripleChain(variant,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: % % variant 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. % https://doi.org/10.1137/070683052 % % [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 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) %
narginchk(0,2); if nargin==0 variant = '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)
t_mess_lradi = tic; outB = mess_lradi(eqn, opts, oper); t_elapsed = toc(t_mess_lradi); fprintf(1,'mess_lradi took %6.2f seconds \n',t_elapsed); if istest if min(outB.res)>=1e-1 error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(1); semilogy(outB.res,'LineWidth',3); title('0 = A X ^T + E X A^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.002105e+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.813608e+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.188758e-01 ADI step: 68 normalized residual: 1.172664e+00 ADI step: 70 normalized residual: 1.129371e+00 ADI step: 71 normalized residual: 8.961118e-01 ADI step: 73 normalized residual: 8.713570e-01 ADI step: 75 normalized residual: 7.378360e-01 ADI step: 77 normalized residual: 5.279292e-01 ADI step: 78 normalized residual: 8.063871e-01 ADI step: 80 normalized residual: 6.666642e-01 ADI step: 81 normalized residual: 7.099135e-01 ADI step: 83 normalized residual: 6.526697e-01 ADI step: 85 normalized residual: 4.615173e-01 ADI step: 87 normalized residual: 4.450099e-01 ADI step: 88 normalized residual: 5.715111e-01 ADI step: 90 normalized residual: 5.039108e-01 ADI step: 92 normalized residual: 3.200317e-01 ADI step: 93 normalized residual: 4.143070e-01 ADI step: 95 normalized residual: 4.217265e-01 ADI step: 97 normalized residual: 3.333128e-01 ADI step: 98 normalized residual: 3.777849e-01 ADI step: 100 normalized residual: 3.617540e-01 ADI step: 101 normalized residual: 3.912230e-01 ADI step: 103 normalized residual: 3.370220e-01 ADI step: 105 normalized residual: 3.317460e-01 ADI step: 107 normalized residual: 3.355299e-01 ADI step: 108 normalized residual: 3.173677e-01 ADI step: 110 normalized residual: 2.591600e-01 ADI step: 111 normalized residual: 2.713324e-01 ADI step: 113 normalized residual: 2.680951e-01 ADI step: 115 normalized residual: 2.582824e-01 ADI step: 117 normalized residual: 2.575210e-01 ADI step: 118 normalized residual: 1.824294e-01 ADI step: 120 normalized residual: 2.049555e-01 ADI step: 121 normalized residual: 2.397734e-01 ADI step: 123 normalized residual: 2.421403e-01 ADI step: 124 normalized residual: 2.247502e-01 ADI step: 126 normalized residual: 2.264447e-01 ADI step: 128 normalized residual: 2.230960e-01 ADI step: 130 normalized residual: 2.271362e-01 ADI step: 132 normalized residual: 2.325429e-01 ADI step: 133 normalized residual: 2.222306e-01 ADI step: 135 normalized residual: 1.468320e-01 ADI step: 137 normalized residual: 1.418196e-01 ADI step: 138 normalized residual: 1.719007e-01 ADI step: 139 normalized residual: 1.616633e-01 ADI step: 141 normalized residual: 1.078697e-01 ADI step: 143 normalized residual: 1.049704e-01 ADI step: 144 normalized residual: 1.892963e-01 ADI step: 146 normalized residual: 1.895386e-01 ADI step: 148 normalized residual: 1.778348e-01 ADI step: 149 normalized residual: 2.060114e-01 ADI step: 151 normalized residual: 2.060300e-01 ADI step: 153 normalized residual: 1.763440e-01 ADI step: 155 normalized residual: 1.804357e-01 ADI step: 157 normalized residual: 1.739263e-01 ADI step: 158 normalized residual: 1.921374e-01 ADI step: 160 normalized residual: 1.461589e-01 ADI step: 162 normalized residual: 1.433190e-01 ADI step: 163 normalized residual: 1.485836e-01 ADI step: 165 normalized residual: 1.227408e-01 ADI step: 167 normalized residual: 1.224209e-01 ADI step: 168 normalized residual: 1.295107e-01 ADI step: 170 normalized residual: 1.159963e-01 ADI step: 172 normalized residual: 1.100151e-01 ADI step: 173 normalized residual: 1.154309e-01 ADI step: 174 normalized residual: 9.360058e-02 ADI step: 176 normalized residual: 9.177159e-02 ADI step: 178 normalized residual: 8.333038e-02 ADI step: 179 normalized residual: 7.422130e-02 ADI step: 181 normalized residual: 7.365669e-02 ADI step: 183 normalized residual: 7.292316e-02 ADI step: 185 normalized residual: 7.307011e-02 ADI step: 187 normalized residual: 7.096115e-02 ADI step: 188 normalized residual: 8.845479e-02 ADI step: 190 normalized residual: 9.420441e-02 ADI step: 192 normalized residual: 9.470731e-02 ADI step: 193 normalized residual: 9.967799e-02 ADI step: 194 normalized residual: 9.501825e-02 ADI step: 196 normalized residual: 9.489766e-02 ADI step: 198 normalized residual: 9.454401e-02 ADI step: 200 normalized residual: 9.479503e-02 ADI step: 202 normalized residual: 9.439792e-02 ADI step: 203 normalized residual: 1.034206e-01 ADI step: 204 normalized residual: 1.074042e-01 ADI step: 206 normalized residual: 1.052745e-01 ADI step: 208 normalized residual: 1.052426e-01 ADI step: 210 normalized residual: 1.051079e-01 ADI step: 212 normalized residual: 1.042939e-01 ADI step: 213 normalized residual: 1.027145e-01 ADI step: 215 normalized residual: 1.020585e-01 ADI step: 217 normalized residual: 1.017780e-01 ADI step: 218 normalized residual: 1.010942e-01 ADI step: 220 normalized residual: 1.002537e-01 ADI step: 222 normalized residual: 1.020959e-01 ADI step: 223 normalized residual: 9.459831e-02 ADI step: 225 normalized residual: 9.295336e-02 ADI step: 226 normalized residual: 7.709465e-02 ADI step: 228 normalized residual: 7.670378e-02 ADI step: 229 normalized residual: 6.854768e-02 ADI step: 231 normalized residual: 6.711718e-02 ADI step: 233 normalized residual: 6.611763e-02 ADI step: 234 normalized residual: 5.228905e-02 ADI step: 236 normalized residual: 5.226955e-02 ADI step: 238 normalized residual: 5.362616e-02 ADI step: 239 normalized residual: 9.468063e-02 ADI step: 241 normalized residual: 9.402312e-02 ADI step: 243 normalized residual: 9.365705e-02 ADI step: 244 normalized residual: 6.628590e-02 ADI step: 246 normalized residual: 6.501904e-02 ADI step: 248 normalized residual: 6.482750e-02 ADI step: 249 normalized residual: 7.485161e-02 ADI step: 251 normalized residual: 7.483962e-02 ADI step: 253 normalized residual: 7.504306e-02 ADI step: 254 normalized residual: 6.256013e-02 ADI step: 256 normalized residual: 6.236421e-02 ADI step: 258 normalized residual: 6.001722e-02 ADI step: 259 normalized residual: 5.742161e-02 ADI step: 261 normalized residual: 5.713425e-02 ADI step: 263 normalized residual: 5.971008e-02 ADI step: 265 normalized residual: 6.031225e-02 ADI step: 267 normalized residual: 6.027308e-02 ADI step: 268 normalized residual: 6.064504e-02 ADI step: 270 normalized residual: 6.397064e-02 ADI step: 272 normalized residual: 6.372679e-02 ADI step: 273 normalized residual: 6.674236e-02 ADI step: 274 normalized residual: 3.606668e-02 ADI step: 276 normalized residual: 3.602915e-02 ADI step: 278 normalized residual: 3.561930e-02 ADI step: 279 normalized residual: 4.706202e-02 ADI step: 281 normalized residual: 4.623834e-02 ADI step: 283 normalized residual: 4.618432e-02 ADI step: 284 normalized residual: 3.994970e-02 ADI step: 286 normalized residual: 3.984032e-02 ADI step: 288 normalized residual: 3.966268e-02 ADI step: 289 normalized residual: 3.880402e-02 ADI step: 291 normalized residual: 3.878791e-02 ADI step: 293 normalized residual: 3.884510e-02 ADI step: 295 normalized residual: 3.844234e-02 ADI step: 297 normalized residual: 3.868438e-02 ADI step: 298 normalized residual: 3.985688e-02 ADI step: 299 normalized residual: 5.061408e-02 ADI step: 301 normalized residual: 5.055444e-02 mess_lradi took 0.55 seconds size outB.Z: 3002 301

switch upper(variant) 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(variant,'FO')) max_ord = 75; tol = eps; inform = 1; [TL,TR] = square_root_method_SO(eqn.M_, max_ord, tol, inform, 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

