Contents

MORLAB Demo: Model Reduction for Discrete-Time Descriptor Systems

This demo script contains the application of the model reduction methods for discrete-time descriptor systems and some hints for proper usage.

%
% 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-2019 Peter Benner, Steffen W. R. Werner
%

General Handling

The MORLAB toolbox implements model reduction methods for discrete-time descriptor systems, i.e., for systems of the form

E*x(t+1) = A*x(t) + B*u(t),
    y(t) = C*x(t) + D*u(t).

In practice, those systems are described by a large number of difference(-algebraic) equations, which makes the usage in optimization and controller design difficult. The goal is to construct an approximation given by

Er*xr(t+1) = Ar*xr(t) + Br*u(t),
     yr(t) = Cr*xr(t) + Dr*u(t),

which is described by less difference(-algebraic) equations than the original system. The new system has to approximate the input-output behavior of the original system, i.e., the same inputs u(t) the outputs are close in a certain sense |y(t) - yr(t)| < tol.

To get started with the model reduction in MORLAB, we load a prepared data file containing a stable descriptor system with no infinite eigenvalues.

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

For the application of the model reduction methods, we just need to give the system to the model reduction function we want to use. This can be done in different ways. First, MORLAB supports the single matrices given to the function. As example we want to apply the classical modal truncation to the loaded matrices.

[Ar, Br, Cr, Dr, Er, info] = ml_dt_dss_mt(A, B, C, D, E);

The output matrices Ar, Br, Cr, Dr, Er are now the system matrices of the approximation. The info struct contains general information about the underlying used methods.

disp(info);
     infoADTF: [1x1 struct]
    infoSIGNM: [1x1 struct]
            N: 0
           Ni: 0

A more advanced way of using the model reduction function is to construct first the system as a struct. This is recommended, since most of the evaluation and analysis tools in MORLAB assume the system to be either in a struct or an ss object. Note that MORLAB doesn't need a sampling time Ts for the model reduction but for the evaluation we will use it.

sys = struct( ...
    'A', A, ...
    'B', B, ...
    'C', C, ...
    'D', D, ...
    'E', E, ...
    'Ts', Ts);

Now we want to apply the balanced truncation method to this system.

rom = ml_dt_dss_bt(sys);

The reduced-order model is now given in the rom struct with the same naming of the matrices as in sys.

disp(rom);
     A: [37x37 double]
     B: [37x3 double]
     C: [2x37 double]
     D: [2x3 double]
     E: [37x37 double]
    Ts: 1

To take a closer look on the result, we can use the analysis tools of MORLAB, e.g., let's compute the transfer functions of the original system and the reduction in the frequency domain by the sigma plot function.

ml_sigmaplot(sys, rom);

Or we can compare the two systems in a time simulation.

ml_dt_dss_simulate(sys, rom);

For creating other plots we are closing the subplot above again.

close;

All model reduction methods allow customization by optional parameters. option contructor. Here we will just call the function for simplicity and adjust the parameters as we want to.

opts = ml_morlabopts('ml_dt_dss_bt');

Now opts is an empty option struct with all the fields that the balanced truncation routine accepts. Note that MORLAB functions doesn't need the full option struct in general, non-existent or empty fields are replaced by default values.

disp(opts);
           DecompEig: []
           DecompTol: []
         dlyapdlopts: [1x1 struct]
        gdlyapdlopts: [1x1 struct]
       ImproperTrunc: []
               Index: []
          infdecopts: [1x1 struct]
              Method: []
               Order: []
    OrderComputation: []
         stabdecopts: [1x1 struct]
     StoreProjection: []
           Tolerance: []

As a simple start, we want to compute a reduced-order model of order 7. Therefor, we are arranging the opts struct in the following way.

opts.Order            = 7;
opts.OrderComputation = 'Order';

Now we can start the model reduction method again but with opts as additional input.

rom = ml_dt_dss_bt(sys, opts);

disp(rom);
     A: [7x7 double]
     B: [7x3 double]
     C: [2x7 double]
     D: [2x3 double]
     E: [7x7 double]
    Ts: 1

And we see that rom now contains the matrices of an order 7 system.

Special Features: Vectorization and Projection Matrices

An important feature of many model reduction methods in MORLAB is the vectorization, i.e., one can compute a bunch of reduced-order models efficiently by one call of the model reduction routine.

As example, we want to compute 3 different reduced-order models using the balanced truncation method. Again we have to adjust the option struct to tell the method how the different reduced-order models should be constructed. For simplicity, we just set 3 different orders.

opts = ml_morlabopts('ml_dt_dss_bt');

opts.OrderComputation = 'Order';
opts.Order            = {10, 15, 20};

rom = ml_dt_dss_bt(sys, opts);

disp(rom);
    [1x1 struct]    [1x1 struct]    [1x1 struct]

rom is now a cell array with the 3 entries corresponding to our chosen order.

disp(rom{1});
disp(rom{2});
disp(rom{3});
     A: [10x10 double]
     B: [10x3 double]
     C: [2x10 double]
     D: [2x3 double]
     E: [10x10 double]
    Ts: 1

     A: [15x15 double]
     B: [15x3 double]
     C: [2x15 double]
     D: [2x3 double]
     E: [15x15 double]
    Ts: 1

     A: [20x20 double]
     B: [20x3 double]
     C: [2x20 double]
     D: [2x3 double]
     E: [20x20 double]
    Ts: 1

Note that the evaluation tools also allow cell arrays of structs and ss objects to be evaluated. So we could just call

ml_sigmaplot(rom)

To see the sigma plots of all 3 computed roms.

The second feature we want to emphasize here is the possibility to get the projection matrices. This feature can be activated for projection-based methods by the optional parameter StoreProjection.

As example, we want to get the projection matrices for using the balanced truncation method. Since we only need this one field, we construct the option struct this time by hand. Remember that the non-existent fields in the struct are replaced by default values.

opts = struct('StoreProjection', 1);

[rom, info] = ml_dt_dss_bt(sys, opts);

disp(rom);
disp(info);
     A: [37x37 double]
     B: [37x3 double]
     C: [2x37 double]
     D: [2x3 double]
     E: [37x37 double]
    Ts: 1

    AbsErrBound: 0.0084
           Hsvi: [0x1 double]
           Hsvp: [80x1 double]
       infoADTF: [1x1 struct]
    infoDLYAPDL: [1x1 struct]
             Ni: 0
             Np: 37
             Nu: 0
              T: [100x37 double]
              W: [37x100 double]

Info now contains the matrices T and W, that can be used to contruct the reduced-order model from the original data:

rom.A = info.W * sys.A * info.T
rom.B = info.W * sys.B
rom.C = sys.C * info.T
rom.D = sys.D
rom.E = info.W * sys.E * info.T

Note that for example the Hankel-norm approximation has not this option, since this is no projection method.

Unstable Systems

In principle, most of the model reduction methods can also handle unstable systems. There are two approaches to this.

For the demonstration we load an unstable standard system, which has 10 anti-stable eigenvalues, 10 infinite and 80 stable ones.

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

Now we perform the classical balanced truncation method on the system.

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

[rom, info] = ml_dt_dss_bt(sys);

If we take a look into info, then we see that the method found the anti-stable and infinite parts and we will see that rom.A is a block diagonal matrix with the anti-stable part in the middle and the infinite part in the lower right corner.

disp(rom);
disp(info);
    A: [58x58 double]
    B: [58x3 double]
    C: [2x58 double]
    D: [2x3 double]
    E: [58x58 double]

     AbsErrBound: 0.0094
            Hsvi: [6x1 double]
            Hsvp: [80x1 double]
        infoADTF: [1x1 struct]
     infoDLYAPDL: [1x1 struct]
    infoGDLYAPDL: [1x1 struct]
              Ni: 3
              Np: 47
              Nu: 8
               T: []
               W: []

But MORLAB also implements model reduction methods designed for unstable systems. The LQG balanced truncation is a methods, which solves Riccati equations to reduce unstable systems. As example, we apply the LQGBT method on the unstable system.

[rom, info] = ml_dt_dss_lqgbt(sys);

disp(rom);
disp(info);
    A: [53x53 double]
    B: [53x3 double]
    C: [2x53 double]
    D: [2x3 double]
    E: [53x53 double]

     AbsErrBound: 0.0077
            Hsvi: [6x1 double]
            Hsvp: [88x1 double]
        infoADTF: [1x1 struct]
      infoDAREDL: [1x1 struct]
    infoGDLYAPDL: [1x1 struct]
              Nf: 50
              Ni: 3
               T: []
               W: []

Note that those two methods can also be used for systems with eigenvalues on the unit disk boundary, which is not possible by classical ones.

If it is in advance known how many unstable eigenvalues exist, one can turn of the additive decomposition for the anti-stable part by opts.stabdecopts.Dimension = 0.

Interface to the System Control Toolbox / Control Package

MORLAB also gives a basic interface for the state-space (ss) object of the System Control Toolbox in MATLAB or the Control Package in Octave. All first-order model reduction methods as well as the system evaluation routines allow the use of ss objects. The following lines of code will only be executed if the System Control Toolbox or the Control Package is installed.

We load again the unstable system as example.

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

The following lines are just testing if the toolboxes are installed and loaded.

hasControlPkg = size(license('inuse', 'control'), 2);
hasControlTbx = license('test', 'control_toolbox');

if hasControlPkg || hasControlTbx

Now we create our system as state-space object.

    sys = dss(A, B, C, D, E, Ts);

And we perform the classical balanced truncation for this system.

    [rom, info] = ml_dt_dss_bt(sys);

The resulting reduced-order model is also saved in an ss object and we can use rom with the MORLAB routines as well as the system control routines.

For the demonstration, we compute the sigma relative error plot for our approximation with a MORLAB function.

    ml_sigmaplot(sys, rom, struct('DiffMode', 'rel'));

And a bode magnitude plot with the function from the System Control Tooblox.

    if hasControlTbx
        bodemag(sys, rom);
    end
end

Remarks

See Also

ml_dt_dss_bt | ml_dt_dss_lqgbt | ml_dt_dss_mt