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!