Model order reduction example

Mario Berljafa and Stefan Güttel, May 2016Download PDF or m-file

Contents

Introduction

This script reproduces the numerical example from [2, Sec. 5.2], which relates to the INLET problem from the Oberwolfach Model Reduction Benchmark Collection [1], an active control model of a supersonic engine inlet; see also [3]. We demonstrate our near-optimal continuation strategy for approximating the transfer function of a dynamical system. In this experiment the number of parallel processors varies.

MATLAB code

Let us load the data and plot the exact transfer function.

if exist('Inlet.A') ~= 2 || exist('Inlet.B') ~= 2 || ...
   exist('Inlet.C') ~= 2 || exist('Inlet.E') ~= 2 || ...
   exist('mmread') ~= 2
   disp(['The required matrices for this problem can be '  ...
         'downloaded from https://portal.uni-freiburg.de/' ...
         'imteksimulation/downloads/benchmark/38866']);
   return
end

N = 11730;
A = mmread('oberwolfach_inlet/Inlet.A');
B = mmread('oberwolfach_inlet/Inlet.B'); B = B(:,1);
C = mmread('oberwolfach_inlet/Inlet.C');
E = mmread('oberwolfach_inlet/Inlet.E');
f = @(s) full(C*((s*E - A)\B));

f0 = 40;
s = linspace(0, f0, 12 + 11*8);
for j = 1:length(s), fs(j) = f(1i*s(j)); end

figure(1)
plot(s, abs(fs), 'k-', 'Color', [0.7,0.7,0.7], 'LineWidth', 4)
xlabel('s'), ylabel('gain |H(is)|'), hold on

Now we use the rat_krylov function simulating a varying number p of parallel processors to compute reduced order models.

P   = [1, 2, 4, 8, 12, 24];
col = {'r', 'b', 'g', 'm','c', 'k'};

ucf = @(AB, nu, mu, x, param) ...
      util_continuation_fom(AB, nu, mu, x, param);

param.continuation       = 'near-optimal';
param.continuation_m     = 5;
param.continuation_root  = inf;
param.continuation_solve = ucf;
param.orth               = 'CGS';
param.reorth             = 0;
param.waitbar            = 1;

for indp = 1:length(P)
  p = P(indp);
  if p > 1
    xi = 1i*repmat(linspace(0, f0, p), 1, 24/p);
  else
    xi = 1i*f0/2*ones(1, 24);
  end

  param.p = p;

  [V, K, H, out] = rat_krylov(A, E, full(A\B), xi, param);

  fprintf('p = %d\n', p)

  % Numerical quantities (cf. [1, Figure 5.2]).
  EV = E*V; AV = A*V; S = A\EV; S = S-V*(V\S); ss = svd(S);
  R = out.R;
  D = fminsearch(@(x) cond(R*diag(x)), ones(size(R, 2), 1), ...
                 struct('Display','off'));
  nrm = norm(V'*V - eye(size(V,2)));

  fprintf('   Cond number (scaled): %.3e\n', cond(R*diag(D)))
  fprintf('   Orthogonality check:  %.3e\n', nrm)
  fprintf('   sigma_2/sigma_1:      %.3e\n\n', ss(2)/ss(1))

  % Evaluate and plot reduced transfer function.
  Em = V'*E*V; Am = V'*A*V; Bm = V'*B; Cm = C*V;
  for j = 1:length(s)
    fsm(j) = (Cm*((1i*s(j)*Em - Am)\Bm));
  end
  plot(s, abs(fsm), '--', 'Color', col{indp})
end

title('INLET - full vs reduced models (m = 24)')
legend('full model', ...
       'p = 1', 'p = 2',  'p = 4', ...
       'p = 8', 'p = 12', 'p = 24')
p = 1
   Cond number (scaled): 1.811e+03
   Orthogonality check:  1.660e-12
   sigma_2/sigma_1:      1.027e-13

p = 2
   Cond number (scaled): 1.097e+04
   Orthogonality check:  2.198e-11
   sigma_2/sigma_1:      2.758e-11

p = 4
   Cond number (scaled): 5.496e+02
   Orthogonality check:  7.550e-13
   sigma_2/sigma_1:      2.901e-13

p = 8
   Cond number (scaled): 6.083e+03
   Orthogonality check:  1.566e-11
   sigma_2/sigma_1:      2.273e-12

p = 12
   Cond number (scaled): 1.192e+04
   Orthogonality check:  1.393e-09
   sigma_2/sigma_1:      9.424e-12

p = 24
   Cond number (scaled): 4.835e+04
   Orthogonality check:  1.092e-01
   sigma_2/sigma_1:      6.682e-11

Links to other examples

Here is a list of other numerical illustrations of parallelization strategies:

Overview of the parallelization options

Numerical illustration from [2, Sec. 3.4]

TEM example from [2, Sec. 5.1]

Waveguide example from [2, Sec. 5.3]

References

[1] Oberwolfach Model Reduction Benchmark Collection, 2003. http://www.imtek.de/simulation/benchmark.

[2] M. Berljafa and S. Güttel. Parallelization of the rational Arnoldi algorithm, SIAM J. Sci. Comput., 39(5):S197--S221, 2017.

[3] G. Lassaux and K. Willcox. Model reduction for active control design using multiple-point Arnoldi methods, AIAA Paper, 616:1--11, 2003.