Contents

function bilinear_BT_rail(refinements)
% Demonstrates the basics of truncated Gramian based balanced truncation
% for the bilinear formulation of the rail model.
%
% Call: bilinear_BT_rail(refinements)
%
% Input:
%   refinement   The grid refinement level determining the problem size.
%                See `mess_get_bilinear_rail` for details.  The default
%                uses the size n = 5177 model.
%                (optional, default: 3)

%
% 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)
%

%
% 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)
%

Octave version check

We need some specific Octave as it seems to ignore mass matrices in old versions and many ode functions.

if exist('OCTAVE_VERSION', 'builtin') && ...
        (not(exist('verLessThan', 'file')) || verLessThan('octave', '5.2'))

    disp('This demo needs at least octave 5.2 with ode15s support.');

    return
end

Set the input argument, if it is not given.

if nargin < 1
    refinements=3;
end

Fetch model data

eqn = mess_get_bilinear_rail(refinements);

Set process parameters

verbosity of the Lyapunov-plus-positive solver

opts.blyap.info = 2;

% Set tolerances and iteration limit for the Lyapunov-plus-positive solver
opts.blyap.res_tol  =  1e-6;
opts.blyap.rel_diff_tol  =  1e-6 ;
opts.blyap.maxiter = 10;

% Set options for the ADI-based inner Lyapunov solver
% and related shift computation
opts.shifts.num_Ritz = 50;
opts.shifts.num_hRitz = 25;
opts.shifts.method = 'projection';
opts.shifts.num_desired = 6;

opts.adi.maxiter = 100;
opts.adi.res_tol = 1e-12;
opts.adi.rel_diff_tol = 0;
opts.adi.info = 2;
opts.adi.norm = 'fro';
opts.norm = 'fro';

% Set mess_res2_norm options
opts.resopts.res.maxiter = 10;
opts.resopts.res.tol = 1e-6;
opts.resopts.res.orth = 0;

opts.fopts.LDL_T = 0;
opts.fopts.norm = 'fro';

opts.srm.tol = 1e-10;
opts.srm.info = 1;

choose USFS set

oper = operatormanager('default');

Perform system reduction

[ROM, ~, eqn, opts, oper]= mess_balanced_truncation_bilinear(eqn, opts ,oper);
ROM.haveE = 0;
% only needed for function bilinear systems to run for both FOM and ROM:
ROM.A_=ROM.A;
ROM.N_=ROM.N;
ADI step:    1 normalized residual: 5.850156e-01 
ADI step:    2 normalized residual: 4.753125e-01 
ADI step:    3 normalized residual: 3.994897e-01 
ADI step:    4 normalized residual: 3.195283e-01 
ADI step:    5 normalized residual: 2.667712e-01 
ADI step:    6 normalized residual: 2.291903e-01 
ADI step:    7 normalized residual: 1.975979e-01 
ADI step:    8 normalized residual: 1.970844e-01 
ADI step:    9 normalized residual: 7.849211e-02 
ADI step:   10 normalized residual: 7.726303e-02 
ADI step:   11 normalized residual: 6.349575e-02 
ADI step:   12 normalized residual: 5.299628e-02 
ADI step:   13 normalized residual: 5.120739e-02 
ADI step:   14 normalized residual: 5.120563e-02 
ADI step:   15 normalized residual: 7.224509e-03 
ADI step:   16 normalized residual: 7.204177e-03 
ADI step:   17 normalized residual: 5.259276e-03 
ADI step:   18 normalized residual: 5.258073e-03 
ADI step:   19 normalized residual: 5.032815e-03 
ADI step:   20 normalized residual: 1.805588e-03 
ADI step:   21 normalized residual: 1.805535e-03 
ADI step:   22 normalized residual: 5.774697e-04 
ADI step:   23 normalized residual: 5.744826e-04 
ADI step:   24 normalized residual: 5.742695e-04 
ADI step:   25 normalized residual: 5.097725e-04 
ADI step:   26 normalized residual: 5.097670e-04 
ADI step:   27 normalized residual: 3.441204e-04 
ADI step:   28 normalized residual: 1.704472e-04 
ADI step:   29 normalized residual: 1.681480e-04 
ADI step:   30 normalized residual: 1.681177e-04 
ADI step:   31 normalized residual: 1.629152e-04 
ADI step:   32 normalized residual: 1.629124e-04 
ADI step:   33 normalized residual: 1.556947e-05 
ADI step:   34 normalized residual: 2.688670e-06 
ADI step:   35 normalized residual: 2.683910e-06 
ADI step:   36 normalized residual: 4.663566e-07 
ADI step:   37 normalized residual: 3.758180e-07 
ADI step:   38 normalized residual: 3.757462e-07 
ADI step:   39 normalized residual: 3.210477e-07 
ADI step:   40 normalized residual: 6.184098e-08 
ADI step:   41 normalized residual: 1.764867e-08 
ADI step:   42 normalized residual: 1.651598e-08 
ADI step:   43 normalized residual: 1.594091e-08 
ADI step:   44 normalized residual: 1.642639e-09 
ADI step:   45 normalized residual: 1.642037e-09 
ADI step:   46 normalized residual: 1.632330e-09 
ADI step:   47 normalized residual: 2.312456e-10 
ADI step:   48 normalized residual: 1.635641e-10 
ADI step:   49 normalized residual: 1.516371e-10 
ADI step:   50 normalized residual: 6.503436e-11 
ADI step:   51 normalized residual: 6.500396e-11 
ADI step:   52 normalized residual: 1.102681e-11 
ADI step:   53 normalized residual: 1.094060e-11 
ADI step:   54 normalized residual: 7.624897e-12 
ADI step:   55 normalized residual: 7.202863e-12 
ADI step:   56 normalized residual: 1.992392e-12 
ADI step:   57 normalized residual: 1.991335e-12 
ADI step:   58 normalized residual: 8.499744e-13 
step:    1 residual: 0.000000e+00 rank(Z): 144 
ADI step:    1 normalized residual: 9.255309e-01 
ADI step:    2 normalized residual: 7.907249e-01 
ADI step:    3 normalized residual: 5.871261e-01 
ADI step:    4 normalized residual: 4.381557e-01 
ADI step:    5 normalized residual: 3.143962e-01 
ADI step:    6 normalized residual: 2.417279e-01 
ADI step:    7 normalized residual: 2.225004e-01 
ADI step:    8 normalized residual: 1.029552e-01 
ADI step:    9 normalized residual: 1.028412e-01 
ADI step:   10 normalized residual: 7.941353e-02 
ADI step:   11 normalized residual: 7.928711e-02 
ADI step:   12 normalized residual: 7.060856e-02 
ADI step:   13 normalized residual: 6.924555e-02 
ADI step:   14 normalized residual: 6.924681e-02 
ADI step:   15 normalized residual: 3.334129e-02 
ADI step:   16 normalized residual: 3.329523e-02 
ADI step:   17 normalized residual: 2.570410e-02 
ADI step:   18 normalized residual: 2.525764e-02 
ADI step:   19 normalized residual: 2.466967e-02 
ADI step:   20 normalized residual: 2.466948e-02 
ADI step:   21 normalized residual: 8.711176e-03 
ADI step:   22 normalized residual: 8.706068e-03 
ADI step:   23 normalized residual: 4.300220e-03 
ADI step:   24 normalized residual: 4.282946e-03 
ADI step:   25 normalized residual: 4.154593e-03 
ADI step:   26 normalized residual: 4.154566e-03 
ADI step:   27 normalized residual: 2.696639e-03 
ADI step:   28 normalized residual: 2.694473e-03 
ADI step:   29 normalized residual: 1.323225e-03 
ADI step:   30 normalized residual: 6.537163e-05 
ADI step:   31 normalized residual: 3.477391e-05 
ADI step:   32 normalized residual: 3.476172e-05 
ADI step:   33 normalized residual: 3.405674e-05 
ADI step:   34 normalized residual: 3.351649e-05 
ADI step:   35 normalized residual: 1.658817e-06 
ADI step:   36 normalized residual: 1.441807e-06 
ADI step:   37 normalized residual: 1.344603e-06 
ADI step:   38 normalized residual: 1.392039e-07 
ADI step:   39 normalized residual: 1.389558e-07 
ADI step:   40 normalized residual: 6.056773e-08 
ADI step:   41 normalized residual: 5.970674e-08 
ADI step:   42 normalized residual: 5.951934e-08 
ADI step:   43 normalized residual: 5.834090e-08 
ADI step:   44 normalized residual: 2.946857e-09 
ADI step:   45 normalized residual: 2.941690e-09 
ADI step:   46 normalized residual: 1.955680e-09 
ADI step:   47 normalized residual: 1.952057e-09 
ADI step:   48 normalized residual: 3.110592e-10 
ADI step:   49 normalized residual: 3.030109e-10 
ADI step:   50 normalized residual: 3.029456e-10 
ADI step:   51 normalized residual: 2.305157e-11 
ADI step:   52 normalized residual: 2.300702e-11 
ADI step:   53 normalized residual: 1.605965e-11 
ADI step:   54 normalized residual: 2.381142e-12 
ADI step:   55 normalized residual: 2.327500e-12 
ADI step:   56 normalized residual: 2.323640e-12 
ADI step:   57 normalized residual: 8.416387e-14 
step:    2 residual: 1.362079e-26 rank(Z): 199 
Converged at step:    2 with residual: 1.362079e-26 rank(Z):199 
Elapsed time is 203.571209 seconds.
ADI step:    1 normalized residual: 3.042920e-01 
ADI step:    2 normalized residual: 1.461599e-01 
ADI step:    3 normalized residual: 8.921908e-02 
ADI step:    4 normalized residual: 6.057646e-02 
ADI step:    5 normalized residual: 4.482928e-02 
ADI step:    6 normalized residual: 3.315240e-02 
ADI step:    7 normalized residual: 2.611895e-02 
ADI step:    8 normalized residual: 2.422521e-02 
ADI step:    9 normalized residual: 3.321623e-03 
ADI step:   10 normalized residual: 3.111151e-03 
ADI step:   11 normalized residual: 1.992631e-03 
ADI step:   12 normalized residual: 1.952416e-03 
ADI step:   13 normalized residual: 1.741349e-03 
ADI step:   14 normalized residual: 1.732727e-03 
ADI step:   15 normalized residual: 2.868363e-04 
ADI step:   16 normalized residual: 2.817366e-04 
ADI step:   17 normalized residual: 1.777327e-04 
ADI step:   18 normalized residual: 1.774337e-04 
ADI step:   19 normalized residual: 1.655852e-04 
ADI step:   20 normalized residual: 6.511069e-05 
ADI step:   21 normalized residual: 6.508239e-05 
ADI step:   22 normalized residual: 5.180204e-05 
ADI step:   23 normalized residual: 5.177502e-05 
ADI step:   24 normalized residual: 5.120464e-05 
ADI step:   25 normalized residual: 5.093405e-05 
ADI step:   26 normalized residual: 5.093376e-05 
ADI step:   27 normalized residual: 6.138270e-06 
ADI step:   28 normalized residual: 3.888903e-06 
ADI step:   29 normalized residual: 3.888128e-06 
ADI step:   30 normalized residual: 3.880681e-06 
ADI step:   31 normalized residual: 3.816710e-06 
ADI step:   32 normalized residual: 5.646882e-07 
ADI step:   33 normalized residual: 5.646461e-07 
ADI step:   34 normalized residual: 9.994580e-08 
ADI step:   35 normalized residual: 9.928148e-08 
ADI step:   36 normalized residual: 9.185417e-08 
ADI step:   37 normalized residual: 2.052937e-08 
ADI step:   38 normalized residual: 2.051883e-08 
ADI step:   39 normalized residual: 1.850808e-08 
ADI step:   40 normalized residual: 1.808808e-08 
ADI step:   41 normalized residual: 1.158856e-08 
ADI step:   42 normalized residual: 9.522061e-09 
ADI step:   43 normalized residual: 3.545576e-09 
ADI step:   44 normalized residual: 2.670415e-09 
ADI step:   45 normalized residual: 2.669277e-09 
ADI step:   46 normalized residual: 2.632094e-09 
ADI step:   47 normalized residual: 1.730874e-09 
ADI step:   48 normalized residual: 5.906592e-10 
ADI step:   49 normalized residual: 2.855743e-10 
ADI step:   50 normalized residual: 2.712909e-10 
ADI step:   51 normalized residual: 2.708522e-10 
ADI step:   52 normalized residual: 1.215551e-10 
ADI step:   53 normalized residual: 1.082785e-10 
ADI step:   54 normalized residual: 5.812570e-13 
step:    1 residual: 0.000000e+00 rank(Z): 185 
ADI step:    1 normalized residual: 7.970838e-01 
ADI step:    2 normalized residual: 6.634989e-01 
ADI step:    3 normalized residual: 5.507982e-01 
ADI step:    4 normalized residual: 4.306935e-01 
ADI step:    5 normalized residual: 2.833548e-01 
ADI step:    6 normalized residual: 1.588064e-01 
ADI step:    7 normalized residual: 1.385566e-01 
ADI step:    8 normalized residual: 5.476261e-02 
ADI step:    9 normalized residual: 5.463959e-02 
ADI step:   10 normalized residual: 2.851662e-02 
ADI step:   11 normalized residual: 2.840280e-02 
ADI step:   12 normalized residual: 2.111007e-02 
ADI step:   13 normalized residual: 2.034876e-02 
ADI step:   14 normalized residual: 2.034683e-02 
ADI step:   15 normalized residual: 1.002348e-02 
ADI step:   16 normalized residual: 9.999533e-03 
ADI step:   17 normalized residual: 5.596786e-03 
ADI step:   18 normalized residual: 3.904233e-03 
ADI step:   19 normalized residual: 3.901650e-03 
ADI step:   20 normalized residual: 3.901591e-03 
ADI step:   21 normalized residual: 1.196946e-03 
ADI step:   22 normalized residual: 1.195790e-03 
ADI step:   23 normalized residual: 6.627942e-04 
ADI step:   24 normalized residual: 5.719791e-04 
ADI step:   25 normalized residual: 5.658318e-04 
ADI step:   26 normalized residual: 5.658138e-04 
ADI step:   27 normalized residual: 3.619835e-04 
ADI step:   28 normalized residual: 3.614669e-04 
ADI step:   29 normalized residual: 1.789279e-04 
ADI step:   30 normalized residual: 1.789146e-04 
ADI step:   31 normalized residual: 1.784227e-04 
ADI step:   32 normalized residual: 1.784207e-04 
ADI step:   33 normalized residual: 3.031264e-07 
ADI step:   34 normalized residual: 2.751264e-07 
ADI step:   35 normalized residual: 5.690918e-08 
ADI step:   36 normalized residual: 3.364051e-08 
ADI step:   37 normalized residual: 1.643682e-08 
ADI step:   38 normalized residual: 1.608495e-08 
ADI step:   39 normalized residual: 1.571400e-08 
ADI step:   40 normalized residual: 1.349732e-08 
ADI step:   41 normalized residual: 7.041645e-09 
ADI step:   42 normalized residual: 5.684835e-10 
ADI step:   43 normalized residual: 4.559154e-10 
ADI step:   44 normalized residual: 4.555300e-10 
ADI step:   45 normalized residual: 1.697340e-10 
ADI step:   46 normalized residual: 1.668208e-10 
ADI step:   47 normalized residual: 3.749107e-11 
ADI step:   48 normalized residual: 2.307260e-11 
ADI step:   49 normalized residual: 2.048191e-11 
ADI step:   50 normalized residual: 9.304975e-12 
ADI step:   51 normalized residual: 9.297416e-12 
ADI step:   52 normalized residual: 2.373141e-12 
ADI step:   53 normalized residual: 2.328949e-12 
ADI step:   54 normalized residual: 9.202005e-14 
step:    2 residual: 1.017419e-11 rank(Z): 679 
Converged at step:    2 with residual: 1.017419e-11 rank(Z):679 
Elapsed time is 342.511548 seconds.
reduced system order: 107  (max possible/allowed: 199/5177)

Evaluation of reduction results

Set ode45 parameters

t0 = 0;
tf = 4500;
tvals = t0:1:tf;

% Start ode45 for original and reduced order systems
u = ones(size(eqn.B,2),1);
x0 = zeros(size(eqn.A_,1),1);
options=odeset(...
    'RelTol', 1e-6, ...
    'Mass', eqn.E_, ...
    'MStateDependence', 'none',...
    'MassSingular', 'no');

[~, x] = ...
    ode15s(@(t,x)bilinear_system(t,x, u, eqn, opts, oper), tvals, x0, options);

options=odeset('RelTol', 1e-6);
u = ones(size(ROM.B,2),1);
x0 = zeros(1,size(ROM.A,1));
[~,x_r] = ...
    ode15s(@(t,x)bilinear_system(t, x, u, ROM, opts, oper),  tvals, x0, options);

% Calculate and plot system outputs
y   = eqn.C*x';
y_r = ROM.C*x_r';

y = y';
y_r = y_r';

leg_comp = {'FOM out 1', 'FOM out 2', 'FOM out 3', 'FOM out 4', 'FOM out 5', ...
    'FOM out 6', 'ROM out 1', 'ROM out 2', 'ROM out 3', 'ROM out 4',...
    'ROM out 5', 'ROM out 6'};
leg_err = {'out 1', 'out 2', 'out 3', 'out 4', 'out 5', 'out 6'};

figure()
plot(tvals, y, '-', 'LineWidth', 3)
hold on
plot(tvals, y_r, '--', 'LineWidth', 3)
xlabel('time [s]')
ylabel('magnitude')
legend(leg_comp, 'Location','northeastoutside')
title('FOM (solid) versus ROM (dashed) outputs')
hold off
% Evaluate absolute error
absErr = abs(y-y_r);
% and corresponding relative error
relErr = abs(absErr ./ y);

% plot relative error
figure()
semilogy(tvals, relErr, 'LineWidth', 3);
xlabel('time [s]')
ylabel('magnitude')
legend(leg_err, 'Location','northeastoutside')
title('pointwise relative output errors')
% Check the relative error
for i = 60:67
    test_condition = relErr(i) < 1e-2;
    if test_condition == 0
        error('limit for relative Error exceeded')
    end
end

fprintf('Everything looks good!\n')
end

helper function for the ode integrator

function f = bilinear_system(~, x ,u ,eqn, opts, oper)

    f = eqn.A_*x + eqn.B*u;

    [eqn, opts, oper] = oper.mul_N_pre(eqn, opts, oper);
    numberOf_N_matrices = length(eqn.N_);

    for currentN_k = 1:numberOf_N_matrices
        f = f + oper.mul_N(eqn, opts, 'N', x, 'N', currentN_k)*u(currentN_k);
    end

    [~, ~, ~] = oper.mul_N_post(eqn, opts, oper);

end
Everything looks good!