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