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 (c) 2009-2025 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 = false; end format long e; % set operation opts = struct(); [oper, opts] = operatormanager(opts, 'so_2'); % Problem data n1 = 500; alpha = 0.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 = true;
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 = false; % opts.adi.accumulateDeltaK = false; opts.adi.compute_sol_fac = true; 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); mess_fprintf(opts, 'mess_lradi took %6.2f seconds \n', t_elapsed); if istest if min(outB.res) >= 1e-1 mess_err(opts, '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 [mZ, nZ] = size(outB.Z); mess_fprintf(opts, 'size outB.Z: %d x %d\n\n', mZ, nZ);
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.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.188757e-01 ADI step: 68 normalized residual: 1.172664e+00 ADI step: 70 normalized residual: 1.129371e+00 ADI step: 71 normalized residual: 8.961106e-01 ADI step: 73 normalized residual: 8.713581e-01 ADI step: 75 normalized residual: 7.378372e-01 ADI step: 77 normalized residual: 5.279300e-01 ADI step: 78 normalized residual: 8.063868e-01 ADI step: 80 normalized residual: 6.666618e-01 ADI step: 81 normalized residual: 7.099150e-01 ADI step: 83 normalized residual: 6.526703e-01 ADI step: 85 normalized residual: 4.615184e-01 ADI step: 87 normalized residual: 4.450107e-01 ADI step: 88 normalized residual: 5.715073e-01 ADI step: 90 normalized residual: 5.039067e-01 ADI step: 92 normalized residual: 3.200372e-01 ADI step: 93 normalized residual: 4.142961e-01 ADI step: 95 normalized residual: 4.217176e-01 ADI step: 97 normalized residual: 3.332969e-01 ADI step: 98 normalized residual: 3.777805e-01 ADI step: 100 normalized residual: 3.617469e-01 ADI step: 101 normalized residual: 3.912357e-01 ADI step: 103 normalized residual: 3.370402e-01 ADI step: 105 normalized residual: 3.317560e-01 ADI step: 107 normalized residual: 3.355225e-01 ADI step: 108 normalized residual: 3.173618e-01 ADI step: 110 normalized residual: 2.591786e-01 ADI step: 111 normalized residual: 2.713283e-01 ADI step: 113 normalized residual: 2.680890e-01 ADI step: 115 normalized residual: 2.582752e-01 ADI step: 117 normalized residual: 2.575143e-01 ADI step: 118 normalized residual: 1.824167e-01 ADI step: 120 normalized residual: 2.048510e-01 ADI step: 121 normalized residual: 2.396396e-01 ADI step: 123 normalized residual: 2.420217e-01 ADI step: 124 normalized residual: 2.250442e-01 ADI step: 126 normalized residual: 2.267417e-01 ADI step: 128 normalized residual: 2.233053e-01 ADI step: 130 normalized residual: 2.276675e-01 ADI step: 132 normalized residual: 2.321377e-01 ADI step: 133 normalized residual: 2.234471e-01 ADI step: 135 normalized residual: 1.473972e-01 ADI step: 137 normalized residual: 1.421866e-01 ADI step: 138 normalized residual: 1.709875e-01 ADI step: 139 normalized residual: 1.586470e-01 ADI step: 141 normalized residual: 9.900877e-02 ADI step: 143 normalized residual: 9.614724e-02 ADI step: 144 normalized residual: 1.857306e-01 ADI step: 146 normalized residual: 1.861647e-01 ADI step: 148 normalized residual: 1.863501e-01 ADI step: 149 normalized residual: 1.995936e-01 ADI step: 151 normalized residual: 1.992495e-01 ADI step: 153 normalized residual: 1.787107e-01 ADI step: 155 normalized residual: 1.668872e-01 ADI step: 157 normalized residual: 1.534941e-01 ADI step: 158 normalized residual: 1.706631e-01 ADI step: 160 normalized residual: 1.645401e-01 ADI step: 162 normalized residual: 1.596682e-01 ADI step: 163 normalized residual: 1.581303e-01 ADI step: 165 normalized residual: 1.586407e-01 ADI step: 167 normalized residual: 1.645811e-01 ADI step: 168 normalized residual: 1.726095e-01 ADI step: 170 normalized residual: 1.726114e-01 ADI step: 172 normalized residual: 1.692694e-01 ADI step: 173 normalized residual: 1.691214e-01 ADI step: 175 normalized residual: 1.703841e-01 ADI step: 177 normalized residual: 1.661153e-01 ADI step: 178 normalized residual: 1.627408e-01 ADI step: 180 normalized residual: 1.603131e-01 ADI step: 182 normalized residual: 1.609220e-01 ADI step: 183 normalized residual: 1.654816e-01 ADI step: 185 normalized residual: 1.652809e-01 ADI step: 187 normalized residual: 1.635957e-01 ADI step: 188 normalized residual: 1.638068e-01 ADI step: 190 normalized residual: 1.645263e-01 ADI step: 192 normalized residual: 1.648147e-01 ADI step: 193 normalized residual: 1.640379e-01 ADI step: 195 normalized residual: 1.610816e-01 ADI step: 197 normalized residual: 1.279600e-01 ADI step: 198 normalized residual: 1.263336e-01 ADI step: 199 normalized residual: 1.137194e-01 ADI step: 201 normalized residual: 1.132135e-01 ADI step: 203 normalized residual: 1.045018e-01 ADI step: 204 normalized residual: 1.395207e-01 ADI step: 206 normalized residual: 1.353857e-01 ADI step: 208 normalized residual: 1.356013e-01 ADI step: 209 normalized residual: 7.226735e-02 ADI step: 211 normalized residual: 7.212180e-02 ADI step: 213 normalized residual: 7.147190e-02 ADI step: 214 normalized residual: 1.019081e-01 ADI step: 216 normalized residual: 1.019346e-01 ADI step: 218 normalized residual: 1.015503e-01 ADI step: 219 normalized residual: 9.629021e-02 ADI step: 221 normalized residual: 9.773470e-02 ADI step: 223 normalized residual: 9.797324e-02 ADI step: 224 normalized residual: 7.189089e-02 ADI step: 226 normalized residual: 7.215250e-02 ADI step: 228 normalized residual: 7.139122e-02 ADI step: 229 normalized residual: 7.682273e-02 ADI step: 231 normalized residual: 7.550024e-02 ADI step: 233 normalized residual: 7.565256e-02 ADI step: 234 normalized residual: 8.643170e-02 ADI step: 236 normalized residual: 8.579344e-02 ADI step: 238 normalized residual: 8.526905e-02 ADI step: 239 normalized residual: 7.558314e-02 ADI step: 241 normalized residual: 7.507395e-02 ADI step: 243 normalized residual: 7.458020e-02 ADI step: 244 normalized residual: 9.377813e-02 ADI step: 246 normalized residual: 9.378112e-02 ADI step: 248 normalized residual: 9.372578e-02 ADI step: 249 normalized residual: 9.792557e-02 ADI step: 251 normalized residual: 9.748247e-02 ADI step: 253 normalized residual: 9.725190e-02 ADI step: 254 normalized residual: 8.359994e-02 ADI step: 256 normalized residual: 8.365475e-02 ADI step: 258 normalized residual: 8.316054e-02 ADI step: 259 normalized residual: 7.309216e-02 ADI step: 261 normalized residual: 7.298501e-02 ADI step: 263 normalized residual: 7.186703e-02 ADI step: 265 normalized residual: 7.128906e-02 ADI step: 266 normalized residual: 6.935223e-02 ADI step: 268 normalized residual: 6.973855e-02 ADI step: 270 normalized residual: 6.975180e-02 ADI step: 272 normalized residual: 6.951310e-02 ADI step: 273 normalized residual: 7.076928e-02 ADI step: 274 normalized residual: 7.909230e-02 ADI step: 276 normalized residual: 7.930569e-02 ADI step: 278 normalized residual: 7.928774e-02 ADI step: 279 normalized residual: 4.677406e-02 ADI step: 281 normalized residual: 4.644850e-02 ADI step: 283 normalized residual: 4.551248e-02 ADI step: 284 normalized residual: 6.268944e-02 ADI step: 286 normalized residual: 6.265200e-02 ADI step: 288 normalized residual: 6.149662e-02 ADI step: 289 normalized residual: 6.020292e-02 ADI step: 291 normalized residual: 5.981246e-02 ADI step: 293 normalized residual: 5.946440e-02 ADI step: 294 normalized residual: 4.795946e-02 ADI step: 296 normalized residual: 4.789738e-02 ADI step: 298 normalized residual: 4.765180e-02 ADI step: 299 normalized residual: 5.558819e-02 ADI step: 301 normalized residual: 5.557317e-02 Warning: LR-ADI reached maximum iteration number. Results may be inaccurate! ↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a> ↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a> ↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a> ↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a> ↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a> ↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/TripleChain/BT_sym_TripleChain.m',107)">BT_sym_TripleChain (line 107)</a> ↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_lradi.m',865)">mess_lradi (line 865)</a> mess_lradi took 0.58 seconds size outB.Z: 3002 x 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.tf_plot.fmin = 1e-4; opts.tf_plot.fmax = 1e0; opts.tf_plot.nsample = 200; if istest opts.tf_plot.info = 1; else opts.tf_plot.info = 2; end opts.tf_plot.type = 'sigma'; out = mess_tf_plot(eqn, opts, oper, ROM); err = out.err; if istest if max(err) > 1000 mess_err(opts, '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

