Contents
MORLAB Demo: Model Reduction for Continuous-Time Standard Systems
This demo script contains the application of the model reduction methods for continuous-time standard 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 continuous-time standard systems, i.e., for systems of the form
x'(t) = 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 differential equations, which makes the usage in optimization and controller design difficult. The goal is to construct an approximation given by
xr'(t) = Ar*xr(t) + Br*u(t), yr(t) = Cr*xr(t) + Dr*u(t),
which is described by less differential 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 general stable standard system.
if exist('OCTAVE_VERSION', 'builtin') orig_warn = warning('off', 'Octave:data-file-in-path'); load morlab_data_std_br.mat; warning(orig_warn); else load morlab_data_std_br.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, info] = ml_ct_ss_mt(A, B, C, D);
The output matrices Ar, Br, Cr, Dr are now the system matrices of the approximation. The info struct contains general information about the underlying used methods.
disp(info);
infoSIGNM: [1x1 struct] infoSYLV: [1x1 struct] N: 1 T: [] W: []
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.
sys = struct( ... 'A', A, ... 'B', B, ... 'C', C, ... 'D', D);
Now we want to apply the balanced truncation method to this system.
rom = ml_ct_ss_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: -0.4978 B: [0.3681 0.3710 0.4207] C: [2x1 double] D: [2x3 double]
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_ct_ss_simulate_ss11(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_ct_ss_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);
lyapdlopts: [1x1 struct] Method: [] Order: [] OrderComputation: [] stabsignmopts: [1x1 struct] stabsylvopts: [1x1 struct] StoreProjection: [] Tolerance: [] UnstabDim: []
As a simple start, we want to compute a reduced-order model of order 5. Therefor, we are arranging the opts struct in the following way.
opts.Order = 5;
opts.OrderComputation = 'Order';
Now we can start the model reduction method again but with opts as additional input.
rom = ml_ct_ss_bt(sys, opts); disp(rom);
A: [5x5 double] B: [5x3 double] C: [2x5 double] D: [2x3 double]
And we see that rom now contains the matrices of an order 5 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 Hankel-norm approxiation 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_ct_ss_hna'); opts.OrderComputation = 'Order'; opts.Order = {1, 3, 5}; rom = ml_ct_ss_hna(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: -0.4980 B: [0.8137 0.8211 0.9310] C: [2x1 double] D: [2x3 double] A: [3x3 double] B: [3x3 double] C: [2x3 double] D: [2x3 double] A: [5x5 double] B: [5x3 double] C: [2x5 double] D: [2x3 double]
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_ct_ss_bt(sys, opts);
disp(rom);
disp(info);
A: -0.4978 B: [0.3681 0.3710 0.4207] C: [2x1 double] D: [2x3 double] AbsErrBound: 7.1059e-04 Hsv: [11x1 double] infoADTF: [1x1 struct] infoLYAPDL: [1x1 struct] Ns: 1 Nu: 0 T: [100x1 double] W: [1x100 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
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 and 90 stable ones.
if exist('OCTAVE_VERSION', 'builtin') orig_warn = warning('off', 'Octave:data-file-in-path'); load morlab_data_std_unstab.mat; warning(orig_warn); else load morlab_data_std_unstab.mat; end
Now we perform the classical balanced truncation method on the system.
sys = struct('A', A, 'B', B, 'C', C, 'D', D); [rom, info] = ml_ct_ss_bt(sys);
If we take a look into info, then we see that the method found the anti-stable part and we will see that rom.A is a block diagonal matrix with the anti-stable part in the lower right corner.
disp(rom); disp(info);
A: [27x27 double] B: [27x3 double] C: [2x27 double] D: [2x3 double] AbsErrBound: 0.0086 Hsv: [52x1 double] infoADTF: [1x1 struct] infoLYAPDL: [1x1 struct] Ns: 17 Nu: 10 T: [] W: []
But MORLAB also implements model reduction methods designed for unstable systems. The LQG and H-Infinity balanced truncation are two methods, which solve Riccati equations to reduce unstable systems. As example, we apply the LQGBT method on the unstable system.
[rom, info] = ml_ct_ss_lqgbt(sys); disp(rom); disp(info);
A: [26x26 double] B: [26x3 double] C: [2x26 double] D: [2x3 double] AbsErrBound: 0.0086 Hsv: [56x1 double] infoCAREDL: [1x1 struct] N: 26 T: [] W: []
Note that those two methods can also be used for systems with eigenvalues on the imaginary axis, 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.UnstabDim = 0.
Passive and Contractive Systems
To preserve beside stability also passivity or contractivity, the positive- and bounded-real methods are implemented in MORLAB.
As example we load a prepared data file with a passive system.
if exist('OCTAVE_VERSION', 'builtin') orig_warn = warning('off', 'Octave:data-file-in-path'); load morlab_data_std_pr.mat; warning(orig_warn); else load morlab_data_std_pr.mat; end
And we use the positive-real balanced truncation to get two reduced-order models by setting tolerances for the error bound of the method to 10^-2 and 10^-5.
sys = struct('A', A, 'B', B, 'C', C, 'D', D); opts = struct('Tolerance', {{1.0e-02, 1.0e-05}}); [rom, info] = ml_ct_ss_prbt(sys, opts); disp(rom); disp(info);
[1x1 struct] [1x1 struct] Hsv: [21x1 double] infoCAREDL: [1x1 struct] InvAbsErrBound: {[1.8384e-05] [4.3675e-06]} M: [3x3 double] N: {[4] [6]} T: [] W: []
As result, we got two passive reduced-order systems of order 4 and 6.
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_std_unstab.mat; warning(orig_warn); else load morlab_data_std_unstab.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 = ss(A, B, C, D);
And we perform the classical balanced truncation for this system.
[rom, info] = ml_ct_ss_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
- All model reduction routines are more or less equivalent to each other, which means that all usage examples here also apply to the other routines.
See Also