function [Ar, Br, Cr] = LQGBT_mor_FDM(tol,max_ord,n0,istest)
% LQGBT_MOR_FDM computes a reduced order model via the linear quadratic % Gaussian balanced truncation [1] for a finite difference discretized % convection diffusion model on the unit square described in [2]. % % Usage: % [Ar, Br, Cr] = LQGBT_mor_FDM(tol,max_ord,n0,test) % % Inputs % % tol truncation tolerance for the LQG characteristic values % % max_ord maximum allowed order for the reduced order model % % n0 n0^2 gives the dimension of the original model, i.e. n0 is % the number of degrees of freedom per spatial direction % % istest flag to determine whether this demo runs as a CI test or % interactive demo % (optional, defaults to 0, i.e. interactive demo) % % Outputs % % Ar, Br, Cr the reduced order system matrices. % % References % [1] D. Mustafa, K. Glover, Controller design by H∞ -balanced truncation, % IEEE Trans. Autom. Control 36 (6) (1991) 668–682. % https://doi.org/10.1109/9.86941 % % [2] T. Penzl, LyaPack Users Guide, Tech. Rep. SFB393/00-33, % Sonderforschungsbereich 393 Numerische Simulation auf massiv % parallelen Rechnern, TU Chemnitz, 09107 Chemnitz, Germany, % available from http://www.tu-chemnitz.de/sfb393/sfb00pr.html (2000). % % % 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,4);
LQGBT Tolerance and maximum order of the ROM.
if nargin<1 tol = 1e-6; end if nargin<2 max_ord = 250; end if nargin<3 n0 = 60; % n0 = number of grid points in either space direction; % n = n0^2 is the problem dimension! % (Change n0 to generate problems of different size.) end if nargin<4 istest=0; end % Problem data eqn.A_ = fdm_2d_matrix(n0,'10*x','100*y','0'); eqn.B = fdm_2d_vector(n0,'.1<x<=.3'); eqn.C = [fdm_2d_vector(n0,'.7<x<=.8'), fdm_2d_vector(n0,'.8<x<=.9')]; eqn.C = eqn.C';
Set operations.
oper = operatormanager('default');
ADI options.
opts.adi.maxiter = 200;
opts.adi.res_tol = 1.0e-10;
opts.adi.rel_diff_tol = 0;
opts.adi.info = 0;
opts.adi.accumulateK = 1;
opts.adi.accumulateDeltaK = 0;
opts.adi.compute_sol_fac = 1;
n = oper.size(eqn, opts);
opts.shifts.num_desired = 25;
opts.shifts.num_Ritz = 50;
opts.shifts.num_hRitz = 25;
opts.shifts.info = 0;
opts.shifts.method = 'heur';
opts.shifts.b0 = ones(n,1);
opts.norm = 'fro'; % Newton options. opts.nm.maxiter = 25; opts.nm.res_tol = 1e-10; opts.nm.rel_diff_tol = 1e-16; opts.nm.info = 1; opts.nm.accumulateRes = 1; opts.nm.linesearch = 1; opts.nm.inexact = 'quadratic'; opts.nm.tau = 0.1;
Solve the filter Riccati equation. A*X + X*A' - X*C'*C*X + B*B' = 0
t_mess_lrnm = tic; eqn.type = 'N'; outC = mess_lrnm(eqn, opts, oper); t_elapsed1 = toc(t_mess_lrnm); fprintf(1,'mess_lrnm took %6.2f seconds \n',t_elapsed1); if istest if min(outC.res)>=opts.nm.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(1); semilogy(outC.res,'LineWidth',3); title('A X + X A^T - X C^T C X + BB^T = 0'); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end
NM step: 1 normalized residual: 8.637261e-01 relative change in K: 1.000000e+00 number of ADI steps: 10 NM step: 2 normalized residual: 1.247428e-01 relative change in K: 7.069622e-01 number of ADI steps: 9 NM step: 3 normalized residual: 1.105634e-02 relative change in K: 2.214896e-01 number of ADI steps: 13 NM step: 4 normalized residual: 2.578772e-04 relative change in K: 3.403586e-02 number of ADI steps: 20 NM step: 5 normalized residual: 2.177720e-07 relative change in K: 9.882326e-04 number of ADI steps: 26 NM step: 6 normalized residual: 3.531935e-11 relative change in K: 6.630748e-07 number of ADI steps: 24 mess_lrnm took 3.00 seconds

Solve the regulator Riccati equation. A'*X + X*A - X*B*B'*X + C'*C = 0
t_mess_lrnm = tic; eqn.type = 'T'; outB = mess_lrnm(eqn, opts, oper); t_elapsed2 = toc(t_mess_lrnm); fprintf(1,'mess_lrnm took %6.2f seconds \n' ,t_elapsed2); if istest if min(outB.res)>=opts.nm.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(2); semilogy(outB.res,'LineWidth',3); title('A^T X + X A - X BB^T X + C^T C = 0'); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end
Using line search (res: 3.850120e+00) lambda: 4.243950e-01 NM step: 1 normalized residual: 3.249656e-01 relative change in K: 1.000000e+00 number of ADI steps: 8 NM step: 2 normalized residual: 2.089833e-02 relative change in K: 1.547483e-01 number of ADI steps: 17 NM step: 3 normalized residual: 8.639340e-05 relative change in K: 1.010884e-02 number of ADI steps: 23 NM step: 4 normalized residual: 5.402815e-10 relative change in K: 2.528390e-05 number of ADI steps: 63 NM step: 5 normalized residual: 3.669164e-11 relative change in K: 9.263228e-11 number of ADI steps: 40 mess_lrnm took 3.03 seconds

% Model reduction by square root method.
opts.srm.tol=tol; opts.srm.max_ord = max_ord; opts.srm.info=1; [TL, TR, ~, eqn, opts, ~] = mess_square_root_method(eqn,opts,oper,... outB.Z,outC.Z); Ar = TR'*(eqn.A_*TL); Br = TR'*eqn.B; Cr = eqn.C*TL; opts.sigma.nsample = 200; % 200 frequency samples opts.sigma.fmin = -2; % min. frequency 1e-3 opts.sigma.fmax = 6; % max. frequency 1e4 if istest opts.sigma.info=1; % no output else opts.sigma.info = 2; % show messages and plots end ROM.A = Ar; ROM.B = Br; ROM.C = Cr; ROM.E = eye(size(ROM.A,1)); out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err;
reduced system order: 14 (max possible/allowed: 21/250) 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


Report.
if istest if max(err)>=tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end end
ans = 1.0e+03 * Columns 1 through 7 -0.1380 -0.2099 0.0195 -0.1024 -0.0126 -0.0702 -0.0058 -0.2099 -0.6657 0.1614 -0.5391 -0.0691 -0.3969 -0.0326 -0.0195 -0.1615 -0.0130 0.2024 0.0193 0.1014 0.0079 -0.1024 -0.5392 -0.2025 -1.3366 -0.1897 -1.5886 -0.1239 -0.0174 -0.0978 -0.0261 -0.3876 -0.1128 -0.3076 -0.0907 -0.0702 -0.3970 -0.1014 -1.5887 -0.7770 -4.2944 0.1431 -0.0092 -0.0521 -0.0129 -0.2245 -0.0920 -1.2150 -0.1074 -0.0227 -0.1296 -0.0315 -0.5725 -0.2808 -2.3780 -0.4903 0.0083 0.0476 0.0115 0.2127 0.1033 0.9464 0.1724 0.0097 0.0556 0.0134 0.2506 0.1249 1.1692 0.2250 0.0007 0.0040 0.0010 0.0179 0.0089 0.0846 0.0161 0.0043 0.0244 0.0059 0.1101 0.0550 0.5190 0.0997 0.0006 0.0035 0.0008 0.0156 0.0078 0.0736 0.0142 0.0007 0.0043 0.0010 0.0193 0.0097 0.0914 0.0175 Columns 8 through 14 -0.0226 0.0069 0.0097 0.0001 0.0043 0.0004 0.0007 -0.1288 0.0395 0.0556 0.0006 0.0244 0.0022 0.0039 0.0315 -0.0093 -0.0134 -0.0001 -0.0059 -0.0006 -0.0009 -0.5703 0.1732 0.2506 0.0021 0.1097 0.0102 0.0176 -0.1697 0.1017 0.0881 0.0070 0.0419 -0.0005 0.0091 -2.3685 0.7363 1.1692 0.0085 0.5170 0.0489 0.0826 -0.1712 0.1626 0.1352 0.0137 0.0672 -0.0027 0.0156 -2.6997 0.5687 2.1127 -0.0045 0.9845 0.1024 0.1547 1.6545 -0.6329 -0.9999 -0.0598 -0.5401 -0.0087 -0.1116 2.1216 -1.3199 -5.4599 0.8362 -4.0392 -0.4079 -0.7495 0.1654 -0.0990 -1.3533 -0.0405 0.0354 0.0536 -0.0437 0.9983 -0.6327 -4.0581 -0.5077 -5.8046 -0.5013 -1.7383 0.1431 -0.0917 -0.6541 -0.0751 -1.4821 -0.3612 0.0077 0.1779 -0.1135 -0.8281 -0.0864 -1.8760 -0.8983 -1.9170