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 (c) 2009-2025 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'))

    mess_fprintf(struct, ...
                 '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 = 2;
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 = false;

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

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

choose USFS set

[oper, opts] = operatormanager(opts, 'default');

Perform system reduction

[ROM, ~, eqn, opts, oper] = mess_balanced_truncation_bilinear(eqn, opts, oper);
ROM.haveE = false;
% 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.852602e-01   
ADI step:    2 normalized residual: 4.725143e-01   
ADI step:    3 normalized residual: 3.957490e-01   
ADI step:    4 normalized residual: 3.154090e-01   
ADI step:    5 normalized residual: 2.629962e-01   
ADI step:    6 normalized residual: 2.253168e-01   
ADI step:    7 normalized residual: 1.733749e-01   
ADI step:    8 normalized residual: 1.730036e-01   
ADI step:    9 normalized residual: 8.852700e-02   
ADI step:   10 normalized residual: 8.592194e-02   
ADI step:   11 normalized residual: 4.866952e-02   
ADI step:   12 normalized residual: 4.863480e-02   
ADI step:   13 normalized residual: 4.796868e-02   
ADI step:   14 normalized residual: 4.796515e-02   
ADI step:   15 normalized residual: 2.409286e-02   
ADI step:   16 normalized residual: 2.399171e-02   
ADI step:   17 normalized residual: 1.107968e-02   
ADI step:   18 normalized residual: 1.374140e-04   
ADI step:   19 normalized residual: 3.396044e-05   
ADI step:   20 normalized residual: 3.385299e-05   
ADI step:   21 normalized residual: 3.074579e-05   
ADI step:   22 normalized residual: 2.279768e-05   
ADI step:   23 normalized residual: 1.329990e-05   
ADI step:   24 normalized residual: 1.305970e-05   
ADI step:   25 normalized residual: 3.911814e-06   
ADI step:   26 normalized residual: 3.844116e-06   
ADI step:   27 normalized residual: 3.840302e-06   
ADI step:   28 normalized residual: 2.154556e-06   
ADI step:   29 normalized residual: 1.809099e-06   
ADI step:   30 normalized residual: 1.790384e-06   
ADI step:   31 normalized residual: 1.191187e-06   
ADI step:   32 normalized residual: 1.190726e-06   
ADI step:   33 normalized residual: 1.098747e-06   
ADI step:   34 normalized residual: 1.358373e-07   
ADI step:   35 normalized residual: 1.346303e-07   
ADI step:   36 normalized residual: 1.743221e-08   
ADI step:   37 normalized residual: 1.065640e-08   
ADI step:   38 normalized residual: 1.064839e-08   
ADI step:   39 normalized residual: 9.510056e-09   
ADI step:   40 normalized residual: 8.907265e-10   
ADI step:   41 normalized residual: 8.582273e-10   
ADI step:   42 normalized residual: 8.550998e-10   
ADI step:   43 normalized residual: 7.610448e-10   
ADI step:   44 normalized residual: 2.045672e-10   
ADI step:   45 normalized residual: 2.043110e-10   
ADI step:   46 normalized residual: 6.444107e-11   
ADI step:   47 normalized residual: 6.342724e-11   
ADI step:   48 normalized residual: 6.321517e-11   
ADI step:   49 normalized residual: 5.551946e-11   
ADI step:   50 normalized residual: 5.545529e-11   
ADI step:   51 normalized residual: 1.471981e-11   
ADI step:   52 normalized residual: 1.449097e-11   
ADI step:   53 normalized residual: 2.853497e-12   
ADI step:   54 normalized residual: 2.838058e-12   
ADI step:   55 normalized residual: 2.519011e-12   
ADI step:   56 normalized residual: 2.517386e-12   
ADI step:   57 normalized residual: 1.670909e-13   
step:    1 residual: 0.000000e+00 rank(Z): 124   
ADI step:    1 normalized residual: 9.262721e-01   
ADI step:    2 normalized residual: 7.939284e-01   
ADI step:    3 normalized residual: 5.952297e-01   
ADI step:    4 normalized residual: 4.475445e-01   
ADI step:    5 normalized residual: 3.241449e-01   
ADI step:    6 normalized residual: 2.518143e-01   
ADI step:    7 normalized residual: 2.288776e-01   
ADI step:    8 normalized residual: 2.286211e-01   
ADI step:    9 normalized residual: 1.410828e-01   
ADI step:   10 normalized residual: 1.399839e-01   
ADI step:   11 normalized residual: 1.109378e-01   
ADI step:   12 normalized residual: 1.109215e-01   
ADI step:   13 normalized residual: 1.065639e-01   
ADI step:   14 normalized residual: 1.065540e-01   
ADI step:   15 normalized residual: 3.494491e-02   
ADI step:   16 normalized residual: 3.486746e-02   
ADI step:   17 normalized residual: 2.101068e-02   
ADI step:   18 normalized residual: 2.100745e-02   
ADI step:   19 normalized residual: 1.996968e-02   
ADI step:   20 normalized residual: 1.996918e-02   
ADI step:   21 normalized residual: 7.280126e-03   
ADI step:   22 normalized residual: 7.227859e-03   
ADI step:   23 normalized residual: 7.345134e-05   
ADI step:   24 normalized residual: 7.313596e-05   
ADI step:   25 normalized residual: 6.521968e-05   
ADI step:   26 normalized residual: 6.518320e-05   
ADI step:   27 normalized residual: 4.239147e-05   
ADI step:   28 normalized residual: 4.166790e-05   
ADI step:   29 normalized residual: 2.231052e-05   
ADI step:   30 normalized residual: 1.568582e-06   
ADI step:   31 normalized residual: 6.123582e-07   
ADI step:   32 normalized residual: 5.944533e-07   
ADI step:   33 normalized residual: 5.444474e-07   
ADI step:   34 normalized residual: 4.473177e-07   
ADI step:   35 normalized residual: 3.962567e-08   
ADI step:   36 normalized residual: 3.881816e-08   
ADI step:   37 normalized residual: 4.531329e-09   
ADI step:   38 normalized residual: 4.323515e-09   
ADI step:   39 normalized residual: 2.158421e-09   
ADI step:   40 normalized residual: 1.440968e-09   
ADI step:   41 normalized residual: 3.093612e-10   
ADI step:   42 normalized residual: 2.590839e-10   
ADI step:   43 normalized residual: 1.733240e-10   
ADI step:   44 normalized residual: 9.625777e-11   
ADI step:   45 normalized residual: 9.571664e-11   
ADI step:   46 normalized residual: 1.727017e-12   
ADI step:   47 normalized residual: 1.630507e-12   
ADI step:   48 normalized residual: 1.404758e-12   
ADI step:   49 normalized residual: 1.110139e-12   
ADI step:   50 normalized residual: 1.103525e-12   
ADI step:   51 normalized residual: 4.596788e-13   
step:    2 residual: 9.080593e-27 rank(Z): 168   
Converged at step:    2 with residual: 9.080593e-27 rank(Z): 168   
Elapsed time is 8.916579 seconds.
ADI step:    1 normalized residual: 3.040808e-01   
ADI step:    2 normalized residual: 1.460416e-01   
ADI step:    3 normalized residual: 8.956017e-02   
ADI step:    4 normalized residual: 5.835787e-02   
ADI step:    5 normalized residual: 3.986876e-02   
ADI step:    6 normalized residual: 2.591982e-02   
ADI step:    7 normalized residual: 1.906199e-02   
ADI step:    8 normalized residual: 6.402335e-03   
ADI step:    9 normalized residual: 2.676688e-03   
ADI step:   10 normalized residual: 1.442800e-03   
ADI step:   11 normalized residual: 1.371090e-03   
ADI step:   12 normalized residual: 8.279338e-04   
ADI step:   13 normalized residual: 7.228789e-04   
ADI step:   14 normalized residual: 7.041602e-04   
ADI step:   15 normalized residual: 2.579749e-04   
ADI step:   16 normalized residual: 2.545852e-04   
ADI step:   17 normalized residual: 1.809466e-04   
ADI step:   18 normalized residual: 1.790073e-04   
ADI step:   19 normalized residual: 1.748071e-04   
ADI step:   20 normalized residual: 1.747866e-04   
ADI step:   21 normalized residual: 7.890061e-05   
ADI step:   22 normalized residual: 3.676306e-05   
ADI step:   23 normalized residual: 3.667969e-05   
ADI step:   24 normalized residual: 8.127336e-07   
ADI step:   25 normalized residual: 3.793238e-07   
ADI step:   26 normalized residual: 2.575559e-07   
ADI step:   27 normalized residual: 2.353880e-07   
ADI step:   28 normalized residual: 1.175716e-07   
ADI step:   29 normalized residual: 1.109699e-07   
ADI step:   30 normalized residual: 9.511325e-08   
ADI step:   31 normalized residual: 3.672500e-08   
ADI step:   32 normalized residual: 3.611592e-08   
ADI step:   33 normalized residual: 3.340876e-08   
ADI step:   34 normalized residual: 1.542078e-08   
ADI step:   35 normalized residual: 1.492578e-08   
ADI step:   36 normalized residual: 1.323623e-08   
ADI step:   37 normalized residual: 7.390861e-09   
ADI step:   38 normalized residual: 7.089013e-09   
ADI step:   39 normalized residual: 6.210005e-09   
ADI step:   40 normalized residual: 6.028678e-09   
ADI step:   41 normalized residual: 6.957560e-10   
ADI step:   42 normalized residual: 6.134987e-10   
ADI step:   43 normalized residual: 5.129509e-10   
ADI step:   44 normalized residual: 5.036310e-10   
ADI step:   45 normalized residual: 7.570670e-11   
ADI step:   46 normalized residual: 7.280566e-11   
ADI step:   47 normalized residual: 1.306112e-11   
ADI step:   48 normalized residual: 1.276137e-11   
ADI step:   49 normalized residual: 1.189979e-11   
ADI step:   50 normalized residual: 7.664440e-13   
step:    1 residual: 0.000000e+00 rank(Z): 157   
ADI step:    1 normalized residual: 8.323664e-01   
ADI step:    2 normalized residual: 6.527352e-01   
ADI step:    3 normalized residual: 4.772238e-01   
ADI step:    4 normalized residual: 3.266597e-01   
ADI step:    5 normalized residual: 2.133796e-01   
ADI step:    6 normalized residual: 1.597953e-01   
ADI step:    7 normalized residual: 1.324924e-01   
ADI step:    8 normalized residual: 1.323701e-01   
ADI step:    9 normalized residual: 6.277101e-02   
ADI step:   10 normalized residual: 6.169991e-02   
ADI step:   11 normalized residual: 2.749827e-02   
ADI step:   12 normalized residual: 2.547262e-02   
ADI step:   13 normalized residual: 2.355286e-02   
ADI step:   14 normalized residual: 1.204765e-02   
ADI step:   15 normalized residual: 1.204749e-02   
ADI step:   16 normalized residual: 5.118936e-03   
ADI step:   17 normalized residual: 5.098945e-03   
ADI step:   18 normalized residual: 4.388951e-04   
ADI step:   19 normalized residual: 1.278458e-04   
ADI step:   20 normalized residual: 1.275838e-04   
ADI step:   21 normalized residual: 1.259613e-04   
ADI step:   22 normalized residual: 1.225461e-04   
ADI step:   23 normalized residual: 6.125543e-05   
ADI step:   24 normalized residual: 5.489724e-05   
ADI step:   25 normalized residual: 2.072566e-05   
ADI step:   26 normalized residual: 2.043932e-05   
ADI step:   27 normalized residual: 2.037695e-05   
ADI step:   28 normalized residual: 1.946105e-05   
ADI step:   29 normalized residual: 9.373888e-06   
ADI step:   30 normalized residual: 8.350686e-06   
ADI step:   31 normalized residual: 7.166759e-06   
ADI step:   32 normalized residual: 1.183588e-06   
ADI step:   33 normalized residual: 1.174956e-06   
ADI step:   34 normalized residual: 3.052019e-07   
ADI step:   35 normalized residual: 3.016719e-07   
ADI step:   36 normalized residual: 2.129306e-07   
ADI step:   37 normalized residual: 1.906235e-07   
ADI step:   38 normalized residual: 3.212059e-08   
ADI step:   39 normalized residual: 3.186044e-08   
ADI step:   40 normalized residual: 9.036800e-09   
ADI step:   41 normalized residual: 8.935192e-09   
ADI step:   42 normalized residual: 6.277655e-09   
ADI step:   43 normalized residual: 5.827587e-09   
ADI step:   44 normalized residual: 5.814391e-09   
ADI step:   45 normalized residual: 4.403359e-10   
ADI step:   46 normalized residual: 4.319751e-10   
ADI step:   47 normalized residual: 2.097258e-10   
ADI step:   48 normalized residual: 2.007870e-11   
ADI step:   49 normalized residual: 1.915797e-11   
ADI step:   50 normalized residual: 2.080172e-12   
ADI step:   51 normalized residual: 1.786382e-12   
ADI step:   52 normalized residual: 1.172227e-12   
ADI step:   53 normalized residual: 7.276950e-13   
step:    2 residual: 1.078911e-11 rank(Z): 501   
Converged at step:    2 with residual: 1.078911e-11 rank(Z): 501   
Elapsed time is 13.195502 seconds.
reduced system order: 92  (max possible/allowed: 168/1357)  
  

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 not(test_condition)
        mess_err(opts, 'relative_error', 'limit for relative Error exceeded');
    end
end

mess_fprintf(opts, '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!