function [Ar, Br, Cr] = lqgbt_mor_FDM(tol,max_ord,n0,istest)
% LQGBT_MOR_FDM computes a reduced order model via the lienar 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 reuced 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 orde 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 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-2017 %
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
tic; eqn.type = 'N'; outC = mess_lrnm(eqn, opts, oper); toc; if istest if min(outC.res)>=opts.nm.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(1); semilogy(outC.res); title('AX + XA^T - XC^TCX + 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 Elapsed time is 2.584722 seconds.

Solve the regulator Riccati equation. A'*X + X*A - X*B*B'*X + C'*C = 0
tic; eqn.type = 'T'; outB = mess_lrnm(eqn, opts, oper); toc; if istest if min(outB.res)>=opts.nm.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(2); semilogy(outB.res); title('A^TX + XA - XBB^TX + C^TC = 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.263361e-11 number of ADI steps: 40 Elapsed time is 2.459279 seconds.

% Model reduction by square root method.
opts.srm.tol=tol; opts.srm.max_ord = max_ord; opts.srm.info=1; [TL, TR, hsv, 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 = Columns 1 through 3 -1.379787850218477e+02 -2.098991930101151e+02 1.945158039664558e+01 -2.098993425497630e+02 -6.657215444470442e+02 1.614371684602637e+02 -1.946667317851269e+01 -1.615077870129067e+02 -1.295849155971231e+01 -1.024027644293116e+02 -5.391996682124566e+02 -2.024614301379341e+02 -1.744205843518507e+01 -9.779198798747869e+01 -2.606563434539390e+01 -7.023755627850534e+01 -3.970330138923267e+02 -1.013966002010898e+02 -9.150460288643462e+00 -5.209958465729899e+01 -1.286698240536472e+01 -2.268014554594061e+01 -1.295551486684216e+02 -3.153584685792218e+01 8.322966600542122e+00 4.759460210521294e+01 1.152776228085584e+01 9.719557585030199e+00 5.562584502265517e+01 1.342564012423486e+01 6.939465733257562e-01 3.972197438524026e+00 9.579607739711724e-01 4.262413356774420e+00 2.439795081839154e+01 5.884498658659648e+00 6.033862060430781e-01 3.453861992598300e+00 8.329305474459288e-01 7.483603779847949e-01 4.283731030368058e+00 1.033040700245889e+00 Columns 4 through 6 -1.023861517535651e+02 -1.257265753365973e+01 -7.020968286808396e+01 -5.390823440111986e+02 -6.913639688064477e+01 -3.968588556025126e+02 2.024322765223696e+02 1.934226789051771e+01 1.013891097414623e+02 -1.336635014831780e+03 -1.896681405832564e+02 -1.588626390628468e+03 -3.875847108236453e+02 -1.128331790966220e+02 -3.076351123487116e+02 -1.588739115888475e+03 -7.769856898243755e+02 -4.294402418865620e+03 -2.244621207331687e+02 -9.202614780827804e+01 -1.214954771348800e+03 -5.724844404607552e+02 -2.808235840601703e+02 -2.378021421495675e+03 2.126703768177391e+02 1.033116404120387e+02 9.463701348687551e+02 2.505551286842463e+02 1.248710273103652e+02 1.169204825284935e+03 1.792397684706062e+01 8.924611522241307e+00 8.464215229317914e+01 1.100708029188558e+02 5.495820945795410e+01 5.190283233035052e+02 1.558633105562935e+01 7.786827973264802e+00 7.362940427766884e+01 1.933216470099979e+01 9.656576562709146e+00 9.135211730811348e+01 Columns 7 through 9 -5.762907629269254e+00 -2.255820725763411e+01 6.899733044648654e+00 -3.258465573456419e+01 -1.288411597600884e+02 3.946654083682273e+01 7.930361340571264e+00 3.147434073436109e+01 -9.302300294728083e+00 -1.239472326186882e+02 -5.702713128348192e+02 1.731975507251304e+02 -9.072347262702132e+01 -1.696823799899645e+02 1.016987789518224e+02 1.431333555576709e+02 -2.368544902722163e+03 7.363474880154822e+02 -1.073835474461003e+02 -1.711611529017456e+02 1.626087911267472e+02 -4.903154184481200e+02 -2.699668828733954e+03 5.687170287191307e+02 1.724395322287991e+02 1.654522865596895e+03 -6.328592216076318e+02 2.249870929749682e+02 2.121595681156901e+03 -1.319948145407470e+03 1.606838746635216e+01 1.654458517805755e+02 -9.903282862547407e+01 9.970848975561816e+01 9.983165364770227e+02 -6.327006238557901e+02 1.415611734941333e+01 1.430680124160949e+02 -9.166796428898252e+01 1.754733163546086e+01 1.778768773584253e+02 -1.134959682950343e+02 Columns 10 through 12 9.717602262181458e+00 9.592151167264153e-02 4.254669380276249e+00 5.561342171295966e+01 5.527709909417748e-01 2.435518294052400e+01 -1.342274925606812e+01 -9.403170850854337e-02 -5.854466430986755e+00 2.505544701611049e+02 2.134927449869502e+00 1.097349851471159e+02 8.807496306129146e+01 6.985852993597978e+00 4.189789960824827e+01 1.169165331058345e+03 8.476629068149602e+00 5.169903579454201e+02 1.351594033378031e+02 1.368531267303926e+01 6.718592439511617e+01 2.112695595640415e+03 -4.524755620151915e+00 9.844808283609590e+02 -9.998790056277833e+02 -5.978103395380602e+01 -5.400936052782504e+02 -5.459930174060410e+03 8.362096707968609e+02 -4.039196379420825e+03 -1.353328209723596e+03 -4.054835313754420e+01 3.537779672140710e+01 -4.058055607424087e+03 -5.076752522830531e+02 -5.804564925906087e+03 -6.541249424131895e+02 -7.508203931425328e+01 -1.482065562089189e+03 -8.281040083041950e+02 -8.641576252209734e+01 -1.876005251054115e+03 Columns 13 through 14 3.877705080129896e-01 6.857053679354591e-01 2.216788100672909e+00 3.926920029770919e+00 -5.608939550648020e-01 -9.288816778976313e-01 1.022678003802499e+01 1.757265460390456e+01 -5.110050536178274e-01 9.067141593922733e+00 4.885020466737890e+01 8.263724241086231e+01 -2.700274831327105e+00 1.559093172725845e+01 1.024288020593655e+02 1.546652355795167e+02 -8.738817096754818e+00 -1.116285562998081e+02 -4.078788012136729e+02 -7.494790145425161e+02 5.361441726831008e+01 -4.365927575017179e+01 -5.012711118639284e+02 -1.738258548102236e+03 -3.612341476083874e+02 7.709498948108468e+00 -8.983095857443517e+02 -1.917034827209667e+03