Contents
- Input checks
- set operation manager for the structured computations of Gramians
- load problem data
- options
- LRADI for the two Gramian factors
- Reduced Order Model computation via square root method (SRM)
- Frequency-domain evaluation of the (transfer function of the)
- final accuracy test used in the continuous integration system or
function BT_DAE3_SO(model, tol, max_ord, maxiter, istest)
BT_DAE3_SO demonstrates the use of balanced truncation for index-3 second order models. It reports the results presented in [3].
Usage: BT_DAE3_SO(model,tol,max_ord,maxiter,test)
Input:
model choice of the example system. Possible values 'Stykel_small' Stykel's mass spring damper system from [1] of dimension 600 in second order form (default) 'Stykel_large' Stykel's mass spring damper system from [1] of dimension 6000 in second order form 'Truhar_Veselic' a version of the triple chain mass spring damper test example from [2] with holonomic constraints coupling certain masses. Dimension 451 in second order form. 'TV' another version of the example from [2] with other constraints 'TV2' another version of the example from [2] with yet another set of constraints
tol Balanced Truncation tolerance (optional, default: 1e-4)
max_ord maximal reduced order allowed (optional, default: 500)
maxiter maximal number of ADI iterations (optional, default: 300)
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] V. Mehrmann and T. Stykel. Balanced truncation model reduction for large-scale systems in descriptor form. In P. Benner, V. Mehrmann, and D. Sorensen, editors, Dimension Reduction of Large-Scale Systems, volume 45 of Lecture Notes in Computational Science and Engineering, pages 83–115. Springer-Verlag, Berlin/Heidelberg, 2005. https://doi.org/10.1007/3-540-27909-1_3
[2] 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
[3] J. Saak, M. Voigt, Model reduction of constrained mechanical systems in M-M.E.S.S., IFAC-PapersOnLine 9th Vienna International Conference on Mathematical Modelling MATHMOD 2018, Vienna, Austria, 21–23 February 2018 51 (2) (2018) 661–666. https://doi.org/10.1016/j.ifacol.2018.03.112
% % 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) %
Input checks
narginchk(0, 5); if nargin < 1 model = 'Stykel_small'; end if nargin < 2 tol = 1e-4; end if nargin < 3 max_ord = 500; end if nargin < 4 maxiter = 300; end if nargin < 5 istest = false; end
set operation manager for the structured computations of Gramians
opts = struct();
[oper, opts] = operatormanager(opts, 'dae_3_so');
load problem data
switch lower(model) case {'stykel_small', 'stykel_large'} if strcmp(model, 'stykel_small') sys = load(sprintf('%s/../models/ms_ind3_by_t_stykel/g600.mat', ... fileparts(mfilename('fullpath')))); else sys = load(sprintf('%s/../models/ms_ind3_by_t_stykel/g6000.mat', ... fileparts(mfilename('fullpath')))); end eqn.M_ = sys.M; eqn.E_ = sys.D; eqn.K_ = sys.K; eqn.G_ = sys.G; eqn.haveE = true; eqn.alpha = -0.02; nv = size(eqn.M_, 1); np = size(eqn.G_, 1); eqn.B = full(sys.B(1:2 * nv, :)); eqn.C = full(sys.C(:, 1:2 * nv)); clear E A B C M D K G; case 'tv2' n1 = 151; % make sure it is odd! alpha = 0.01; v = 5e0; [eqn.M_, eqn.E_, eqn.K_] = triplechain_MSD(n1, alpha, alpha, v); eqn.E_ = -eqn.E_; eqn.K_ = -eqn.K_; eqn.G_ = sparse(3, 3 * n1 + 1); n12 = ceil(n1 / 2); % n1 is odd so this is the index of % the center of the string eqn.G_(1, 1) = -1; eqn.G_(1, n12) = 2; eqn.G_(1, n1) = -1; eqn.G_(2, n1 + 1) = -1; eqn.G_(2, n1 + n12) = 2; eqn.G_(2, 2 * n1) = -1; eqn.G_(3, 2 * n1 + 1) = -1; eqn.G_(3, 2 * n1 + n12) = 2; eqn.G_(3, 3 * n1) = -1; nv = size(eqn.M_, 1); np = size(eqn.G_, 1); eqn.B = zeros(6 * n1 + 2, 3); eqn.B(3 * n1 + 6, 1) = 1; eqn.B(5 * n1 + n12 + 6, 2) = 1; eqn.B(end - 6, 3) = 1; eqn.C = zeros(3, 6 * n1 + 2); eqn.C(1, 3 * n1 + 1) = 1; eqn.C(2, 2 * n1) = 1; eqn.C(3, 2 * n1 + floor(n12 / 2)) = 1; eqn.haveE = true; eqn.alpha = -0.02; case 'tv' n1 = 151; % make sure it is odd! alpha = 0.01; v = 5e0; [eqn.M_, eqn.E_, eqn.K_] = triplechain_MSD(n1, alpha, alpha, v); eqn.E_ = -eqn.E_; eqn.K_ = -eqn.K_; eqn.G_ = sparse(3, 3 * n1 + 1); n12 = ceil(n1 / 2); % n1 is odd so this is the index of % the center of the string eqn.G_(1, 1) = 1 / 3; eqn.G_(1, n12) = 1 / 3; eqn.G_(1, n1) = 1 / 3; eqn.G_(2, n1 + 1) = 1 / 3; eqn.G_(2, n1 + n12) = 1 / 3; eqn.G_(2, 2 * n1) = 1 / 3; eqn.G_(3, 2 * n1 + 1) = 1 / 3; eqn.G_(3, 2 * n1 + n12) = 1 / 3; eqn.G_(3, 3 * n1) = 1 / 3; nv = size(eqn.M_, 1); np = size(eqn.G_, 1); eqn.B = [zeros(3 * n1 + 1, 1); ones(3 * n1 + 1, 1)]; eqn.C = [ones(3 * n1 + 1, 1); zeros(3 * n1 + 1, 1)]'; eqn.haveE = true; eqn.alpha = -0.02; case 'truhar_veselic' n1 = 1500; % make sure it is even! alpha = 0.01; v = 5e0; [eqn.M_, eqn.E_, eqn.K_] = triplechain_MSD(n1, alpha, alpha, v); eqn.E_ = -eqn.E_; eqn.K_ = -eqn.K_; eqn.G_ = sparse(6, 3 * n1 + 1); eqn.G_(1, 1) = 1; eqn.G_(1, n1 / 2) = -1; eqn.G_(2, n1 / 2 + 1) = 1; eqn.G_(2, n1) = -1; eqn.G_(3, n1 + 1) = 1; eqn.G_(3, n1 + n1 / 2) = -1; eqn.G_(4, n1 + n1 / 2 + 1) = 1; eqn.G_(4, 2 * n1) = -1; eqn.G_(5, 2 * n1 + 1) = 1; eqn.G_(5, 2 * n1 + n1 / 2) = -1; eqn.G_(6, 2 * n1 + n1 / 2 + 1) = 1; eqn.G_(6, 3 * n1) = -1; nv = size(eqn.M_, 1); np = size(eqn.G_, 1); eqn.B = [zeros(3 * n1 + 1, 1); ones(3 * n1 + 1, 1)]; eqn.C = [ones(3 * n1 + 1, 1); zeros(3 * n1 + 1, 1)]'; eqn.haveE = true; eqn.alpha = -0.02; otherwise mess_err(opts, 'unknown model requested!\n'); return end
options
Adi options
opts.adi.maxiter = maxiter; opts.adi.res_tol = 1e-10; opts.adi.rel_diff_tol = 1e-11; opts.norm = 'fro'; opts.shifts.method = 'projection'; opts.shifts.num_desired = 25;
LRADI for the two Gramian factors
controllability Gramian
eqn.type = 'T'; opts.adi.info = 1; t_mess_lradi1 = tic; [p, ~, eqn, opts, oper] = mess_para(eqn, opts, oper); % use an additional alpha-shift to improve convergence and ROM quality for % the triple chain model if strcmp(model, 'Truhar_Veselic') || strcmp(model, 'TV') || strcmp(model, 'TV2') opts.shifts.p = p - 0.5; else opts.shifts.p = p; end outC = mess_lradi(eqn, opts, oper); t_elapsed1 = toc(t_mess_lradi1); mess_fprintf(opts, 'mess_lradi took %6.2f seconds \n\n', t_elapsed1); % observability Gramian eqn.type = 'N'; t_mess_lradi2 = tic; outB = mess_lradi(eqn, opts, oper); t_elapsed2 = toc(t_mess_lradi2); mess_fprintf(opts, 'mess_lradi took %6.2f seconds \n', t_elapsed2);
Warning: Could not compute initial projection shifts. Going to retry with random right hand side. ↳ 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/DAE3_SO/BT_DAE3_SO.m',214)">BT_DAE3_SO (line 214)</a> ↳ In <a href="matlab:opentoline('/builds/mess/mmess/shifts/mess_para.m',379)">mess_para (line 379)</a> ADI step: 1 normalized residual: 1.215480e+01 relative change in Z: 1.000000e+00 ADI step: 2 normalized residual: 3.525543e+00 relative change in Z: 5.520973e-01 ADI step: 3 normalized residual: 1.953222e+00 relative change in Z: 1.898487e-01 ADI step: 4 normalized residual: 2.131585e+00 relative change in Z: 1.324398e-01 ADI step: 5 normalized residual: 2.435279e+00 relative change in Z: 1.566566e-01 ADI step: 7 normalized residual: 2.191473e-01 relative change in Z: 4.227025e-01 ADI step: 9 normalized residual: 8.229267e-02 relative change in Z: 1.761189e-01 ADI step: 11 normalized residual: 1.247249e-02 relative change in Z: 6.986264e-02 ADI step: 12 normalized residual: 1.309187e-02 relative change in Z: 1.825205e-02 ADI step: 13 normalized residual: 2.110282e-03 relative change in Z: 7.668970e-03 ADI step: 15 normalized residual: 1.666456e-04 relative change in Z: 1.817214e-02 ADI step: 17 normalized residual: 1.223639e-04 relative change in Z: 3.929405e-03 ADI step: 19 normalized residual: 1.359440e-04 relative change in Z: 3.259148e-03 ADI step: 21 normalized residual: 1.010439e-04 relative change in Z: 2.257748e-03 ADI step: 23 normalized residual: 3.494455e-05 relative change in Z: 1.389120e-03 ADI step: 25 normalized residual: 5.723896e-06 relative change in Z: 6.988380e-04 ADI step: 27 normalized residual: 4.484029e-07 relative change in Z: 5.062538e-04 ADI step: 29 normalized residual: 1.727537e-07 relative change in Z: 2.856449e-04 ADI step: 31 normalized residual: 9.661411e-09 relative change in Z: 8.582632e-05 ADI step: 33 normalized residual: 1.434879e-11 relative change in Z: 1.510282e-05 mess_lradi took 0.47 seconds ADI step: 1 normalized residual: 3.133718e-01 relative change in Z: 1.000000e+00 ADI step: 2 normalized residual: 8.672847e-02 relative change in Z: 5.044804e-01 ADI step: 3 normalized residual: 1.460907e-01 relative change in Z: 5.012199e-01 ADI step: 4 normalized residual: 1.035263e-01 relative change in Z: 3.622584e-01 ADI step: 6 normalized residual: 1.030719e-02 relative change in Z: 4.069064e-01 ADI step: 8 normalized residual: 2.348522e-03 relative change in Z: 1.250323e-01 ADI step: 10 normalized residual: 3.960126e-04 relative change in Z: 5.286041e-02 ADI step: 12 normalized residual: 2.940750e-05 relative change in Z: 3.112750e-02 ADI step: 14 normalized residual: 2.145534e-05 relative change in Z: 7.172585e-03 ADI step: 16 normalized residual: 1.375092e-05 relative change in Z: 2.387882e-03 ADI step: 18 normalized residual: 3.776225e-06 relative change in Z: 1.570942e-03 ADI step: 20 normalized residual: 2.663151e-07 relative change in Z: 1.012300e-03 ADI step: 22 normalized residual: 2.130823e-08 relative change in Z: 6.299396e-04 ADI step: 24 normalized residual: 1.533744e-10 relative change in Z: 1.394598e-04 ADI step: 26 normalized residual: 3.735751e-11 relative change in Z: 2.349302e-05 mess_lradi took 0.22 seconds
Reduced Order Model computation via square root method (SRM)
fprintf('\nComputing reduced order model via square root method\n\n'); opts.srm.tol = tol; opts.srm.max_ord = max_ord; opts.srm.info = 1; t_SRM_ROM = tic; [TL, TR] = mess_square_root_method(eqn, opts, oper, outB.Z, outC.Z); % compute ROM matrices ROM.A = TL' * oper.mul_A(eqn, opts, 'N', TR, 'N'); ROM.B = TL' * eqn.B; ROM.C = eqn.C * TR; ROM.E = eye(size(ROM.A)); t_elapsed3 = toc(t_SRM_ROM); mess_fprintf(opts, 'ROM matrices computation took %6.2f seconds \n\n', ... t_elapsed3);
Computing reduced order model via square root method reduced system order: 9 (max possible/allowed: 26/500) ROM matrices computation took 0.02 seconds
Frequency-domain evaluation of the (transfer function of the)
ROM and comparison to the original model.
We feed the mess_tf_plot with usfs that do not exploit the DAE structure:
t_FD_eval = tic; opts.tf_plot.nsample = 200; if istest opts.tf_plot.info = 0; else opts.tf_plot.info = 2; end if strcmp(model, 'TV2') opts.tf_plot.fmin = -2; opts.tf_plot.fmax = 1; else opts.tf_plot.fmin = -4; opts.tf_plot.fmax = 4; end opts.tf_plot.type = 'sigma'; NG = sparse(np, nv); NS = sparse(np, np); eqnu.M_ = [eqn.M_ NG'; NG NS]; eqnu.E_ = [-eqn.E_ NG'; NG NS]; eqnu.K_ = [-eqn.K_ -eqn.G_'; -eqn.G_ NS]; eqnu.C = [zeros(size(eqn.C, 1), np + nv), eqn.C(:, 1:nv), zeros(size(eqn.C, 1), np)]; eqnu.B = [eqn.B(nv + 1:2 * nv, :); zeros(2 * np + nv, size(eqn.B, 2))]; eqnu.haveE = true; [operu, opts] = operatormanager(opts, 'so_1'); out = mess_tf_plot(eqnu, opts, operu, ROM); err = out.err; t_elapsed4 = toc(t_FD_eval); mess_fprintf(opts, 'frequency-domain evaluation took %6.2f \n\n', t_elapsed4);
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 frequency-domain evaluation took 3.13


final accuracy test used in the continuous integration system or
plot of the computed
if istest % the errors are not always perfect in this example, but let's see % whether they are "good enough"... if max(err) > (50 * tol) mess_err(opts, 'TEST:accuracy', ... 'unexpectedly inaccurate result for %s %g %d %d (%g)', ... model, tol, max_ord, maxiter, max(err)); end end