MATLAB Examples

Contents

HINFBT for descriptor systems

This demo script contains the application of the generalized H-infinity balanced truncation method (ml_hinfbt_dss) on a test descriptor system with finite stable, finite unstable and infinite eigenvalues of the form

$$
\setlength\arraycolsep{2pt}
\begin{array}{rl}
E\dot{x}(t) & = Ax(t) + Bu(t),\\
y(t) & = Cx(t) + Du(t).
\end{array}
$$

After loading the demo data, the optional parameters are assigned here explicitly and the ml_hinfbt_dss function is called with its different input interfaces. The dss object version is only called if the System Control Toolbox (Matlab) or the Control Package (Octave) is installed/loaded.

To show the performance of the model reduction method, the sigma error of the full-order and reduced-order model is plotted. In case the System Control Toolbox is installed, also a bode magnitude plot of the error system is shown. The plots of the scaled coprime factorizations like in the standard case cannot be shown, since there is no Matlab implementation of the corpime factorization for descriptor systems.

%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU Affero General Public License as published
% by the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Affero General Public License for more details.
%
% You should have received a copy of the GNU Affero General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) 2006-2017 Peter Benner, Steffen W. R. Werner
%

Initialization

For demonstration, a random descriptor system example was generated by the script morlab_data_desc_infunstab.m and saved in morlab_data_desc_infunstab.mat. The number of finite stable, finite unstable and infinite eigenvalues as well as the complete size of the system is saved in the datainfo structure.

if exist('OCTAVE_VERSION', 'builtin')
    warning('off', 'Octave:data-file-in-path');
    load morlab_data_desc_infunstab.mat;
    warning('on', 'Octave:data-file-in-path');
else
    load morlab_data_desc_infunstab.mat;
end

% Get information about installed/loaded toolboxes.
hasControlPkg = size(license('inuse', 'control'), 2);
hasControlTbx = license('test', 'control_toolbox');

Construction of the descriptor system structure

To test the different input-output formats, the struct and state-space object shapes of the descriptor system are formulated here.

sys_struct = struct('A', A, 'B', B, 'C', C, 'D', D, 'E', E);

if hasControlPkg || hasControlTbx
    sys_ss = dss(A, B, C, D, E);
end

Construction of the normalized LQG system and GAM computation

The optimal gamma from the normalized LQG system couldn't be computed in Matlab due to the miss of the necessary subroutines for descriptor systems. Therefore, the optimal GAM number is set to infinity. If you have an appropriate subroutine for the computation, the normalized LQG system is commented below.

% if hasControlTbx || hasControlPkg
%     n = size(A, 1);
%     m = size(B, 2);
%     p = size(C, 1);
%
%     sys_normal = dss(A, ...
%         [B, zeros(n, p), B], ...
%         [C; zeros(m, n); C], ...
%         [zeros(p, m), zeros(p, p), zeros(p, m);
%         zeros(m, m), zeros(m, p), eye(m, m);
%         zeros(p, m), eye(p, p), D], ...
%         E);
% end

GAM = Inf;

Set of optional parameters

The default values are mainly taken here, which can be modified. Alternative values depending on the system are commented out. Usually for using default values the corresponding parameters are not set or empty. Also, the function call "opts = ml_morlabopts('ml_hinfbt_dss')" generates an empty option struct of the following form.

% Option struct for generalized discrete-time Lyapunov equation solver.
gdlyapopts = struct(...
    'AbsTol' , 0, ...
    'Index'  , Inf, ... % Index = 3
    'Info'   , 0, ... % Info = 1
    'MaxIter', 100, ...
    'RelTol' , 1.0e+01 * datainfo.ni * eps);

% Option struct for Lyapunov solver inside the Riccati solver.
lyapopts = struct(...
    'AbsTol' , 0, ...
    'CompTol', log(datainfo.n)*eps, ...
    'Info'   , 0, ... % Info = 1
    'MaxIter', 100, ...
    'RelTol' , 1.0e+01 * datainfo.n * eps);

% Option struct for Riccati equation solver.
careopts = struct(...
    'AbsTol'  , 0, ...
    'Info'    , 0, ... % Info = 1
    'lyapopts', lyapopts, ...
    'MaxIter' , 100, ...
    'RelTol'  , 1.0e+02 * datainfo.n * eps);

% Option struct for decomposition of infinite part.
infdecopts = struct(...
    'AbsTol'   , 0, ...
    'Dimension', -1, ... % Dimension = datainfo.ni
    'Info'     , 0, ... % Info = 1
    'MaxIter'  , 100, ...
    'RelTol'   , 1.0e+02 * datainfo.n * eps, ...
    'RankTol'  , log(datainfo.n) * eps);

% Option struct for decomposition of unstable part.
stabdecopts = struct(...
    'AbsTol'   , 0, ...
    'Dimension', -1, ... % Dimension = 10
    'Info'     , 0, ... % Info = 1
    'MaxIter'  , 100, ...
    'RelTol'   , 1.0e+02 * datainfo.n * eps, ...
    'RankTol'  , log(datainfo.n) * eps);

% Option struct for the partial stabilization method.
stabmethodopts = struct(...
    'AbsTol' , 0, ...
    'Info'   , 0, ... % Info = 1
    'MaxIter', 100, ...
    'RelTol' , 1.0e+01 * datainfo.n * eps);

% Option struct for the complete function.
opts = struct(...
    'Beta'            , 0.1, ...
    'careopts'        , careopts, ...
    'DecompTol'       , log(datainfo.n) * eps, ...
    'Gamma'           , GAM, ...
    'gdlyapopts'      , gdlyapopts, ...
    'ImproperTrunc'   , log(datainfo.n) * eps, ...
    'Index'           , Inf, ... % Index = 3
    'infdecopts'      , infdecopts, ...
    'Method'          , 'sr', ... % Method = 'bfsr'
    'Order'           , 10, ...
    'OrderComputation', 'tolerance', ... % OrderComputation = 'order'
    'stabdecopts'     , stabdecopts, ...
    'stabmethodopts'  , stabmethodopts, ...
    'StabMethod'      , 'abe', ... % StabMethod = 'lyap'
    'Tolerance'       , 1.0e-02);

Application of the function

Here the application of the ml_hinfbt_dss function is shown for different interfaces and input-data. The default calls are commented out.

% Application with single matrices.
% [Ar, Br, Cr, Dr, Er, info] = ml_hinfbt_dss(A, B, C, D, E);
[Ar, Br, Cr, Dr, Er, info] = ml_hinfbt_dss(A, B, C, D, E, opts);

% Application with structure.
% [rom_struct, info_struct] = ml_hinfbt_dss(sys_struct);
[rom_struct, info_struct] = ml_hinfbt_dss(sys_struct, opts);

% Application with state-space object.
if hasControlPkg || hasControlTbx
    % [rom_ss, info_ss] = ml_hinfbt_dss(sys_ss);
    [rom_ss, info_ss] = ml_hinfbt_dss(sys_ss, opts);
end

Report

As visualization, a sigmaplot of the error system of the (right) coprime factors is made for the standard system structures and a bode magnitude plot for the state-space objects. Note: Since there is no MATLAB implementation of the generalized coprime factorization for descriptor systems, the computed error bound is ignored.

% Sigmaplot of the error system.
figure;
ml_sigmaplot(sys_struct, rom_struct, -4, 4, 100, [], 'b.');
legend('Sigma error');
title({'HINFBT\_DSS (sigmaplot, error system)'; ...
    ['Full order = ' int2str(size(A, 1)) '; ' ...
    'Reduced-order = ' int2str(size(Ar, 1))]});

if hasControlTbx
    % Bode magnitude plot of the error system.
    bodeopts          = bodeoptions('cstprefs');
    bodeopts.MagUnits = 'abs';

    figure;
    bodemag(sys_ss - rom_ss, bodeopts);
    title({'HINFBT\_DSS (Bode magnitude plot, error system)'; ...
        ['Full order = ' int2str(size(A, 1)) '; ' ...
        'Reduced-order = ' int2str(size(Ar, 1))]});
end