National Instruments Noise Reduction Machine 370755B 01 User Manual

TM  
MATRIXx  
TM  
Xmath Model Reduction Module  
Xmath Model Reduction Module  
April 2004 Edition  
Part Number 370755B-01  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Important Information  
Warranty  
The media on which you receive National Instruments software are warranted not to fail to execute programming instructions, due to defects  
in materials and workmanship, for a period of 90 days from date of shipment, as evidenced by receipts or other documentation. National  
Instruments will, at its option, repair or replace software media that do not execute programming instructions if National Instruments receives  
notice of such defects during the warranty period. National Instruments does not warrant that the operation of the software shall be  
uninterrupted or error free.  
A Return Material Authorization (RMA) number must be obtained from the factory and clearly marked on the outside of the package before  
any equipment will be accepted for warranty work. National Instruments will pay the shipping costs of returning to the owner parts which are  
covered by warranty.  
National Instruments believes that the information in this document is accurate. The document has been carefully reviewed for technical  
accuracy. In the event that technical or typographical errors exist, National Instruments reserves the right to make changes to subsequent  
editions of this document without prior notice to holders of this edition. The reader should consult National Instruments if errors are suspected.  
In no event shall National Instruments be liable for any damages arising out of or related to this document or the information contained in it.  
EXCEPT AS SPECIFIED HEREIN, NATIONAL INSTRUMENTS MAKES NO WARRANTIES, EXPRESS OR IMPLIED, AND SPECIFICALLY DISCLAIMS ANY WARRANTY OF  
MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. CUSTOMERS RIGHT TO RECOVER DAMAGES CAUSED BY FAULT OR NEGLIGENCE ON THE PART OF  
NATIONAL INSTRUMENTS SHALL BE LIMITED TO THE AMOUNT THERETOFORE PAID BY THE CUSTOMER. NATIONAL INSTRUMENTS WILL NOT BE LIABLE FOR  
DAMAGES RESULTING FROM LOSS OF DATA, PROFITS, USE OF PRODUCTS, OR INCIDENTAL OR CONSEQUENTIAL DAMAGES, EVEN IF ADVISED OF THE POSSIBILITY  
THEREOF. This limitation of the liability of National Instruments will apply regardless of the form of action, whether in contract or tort, including  
negligence. Any action against National Instruments must be brought within one year after the cause of action accrues. National Instruments  
shall not be liable for any delay in performance due to causes beyond its reasonable control. The warranty provided herein does not cover  
damages, defects, malfunctions, or service failures caused by owner’s failure to follow the National Instruments installation, operation, or  
maintenance instructions; owner’s modification of the product; owner’s abuse, misuse, or negligent acts; and power failure or surges, fire,  
flood, accident, actions of third parties, or other events outside reasonable control.  
Copyright  
Under the copyright laws, this publication may not be reproduced or transmitted in any form, electronic or mechanical, including photocopying,  
recording, storing in an information retrieval system, or translating, in whole or in part, without the prior written consent of National  
Instruments Corporation.  
Trademarks  
MATRIXx, National Instruments, NI, ni.com, and Xmathare trademarks of National Instruments Corporation.  
Product and company names mentioned herein are trademarks or trade names of their respective companies.  
Patents  
For patents covering National Instruments products, refer to the appropriate location: Help»Patents in your software, the patents.txtfile  
on your CD, or ni.com/patents.  
WARNING REGARDING USE OF NATIONAL INSTRUMENTS PRODUCTS  
(1) NATIONAL INSTRUMENTS PRODUCTS ARE NOT DESIGNED WITH COMPONENTS AND TESTING FOR A LEVEL OF  
RELIABILITY SUITABLE FOR USE IN OR IN CONNECTION WITH SURGICAL IMPLANTS OR AS CRITICAL COMPONENTS IN  
ANY LIFE SUPPORT SYSTEMS WHOSE FAILURE TO PERFORM CAN REASONABLY BE EXPECTED TO CAUSE SIGNIFICANT  
INJURY TO A HUMAN.  
(2) IN ANY APPLICATION, INCLUDING THE ABOVE, RELIABILITY OF OPERATION OF THE SOFTWARE PRODUCTS CAN BE  
IMPAIRED BY ADVERSE FACTORS, INCLUDING BUT NOT LIMITED TO FLUCTUATIONS IN ELECTRICAL POWER SUPPLY,  
COMPUTER HARDWARE MALFUNCTIONS, COMPUTER OPERATING SYSTEM SOFTWARE FITNESS, FITNESS OF COMPILERS  
AND DEVELOPMENT SOFTWARE USED TO DEVELOP AN APPLICATION, INSTALLATION ERRORS, SOFTWARE AND  
HARDWARE COMPATIBILITY PROBLEMS, MALFUNCTIONS OR FAILURES OF ELECTRONIC MONITORING OR CONTROL  
DEVICES, TRANSIENT FAILURES OF ELECTRONIC SYSTEMS (HARDWARE AND/OR SOFTWARE), UNANTICIPATED USES OR  
MISUSES, OR ERRORS ON THE PART OF THE USER OR APPLICATIONS DESIGNER (ADVERSE FACTORS SUCH AS THESE ARE  
HEREAFTER COLLECTIVELY TERMED “SYSTEM FAILURES”). ANY APPLICATION WHERE A SYSTEM FAILURE WOULD  
CREATE A RISK OF HARM TO PROPERTY OR PERSONS (INCLUDING THE RISK OF BODILY INJURY AND DEATH) SHOULD  
NOT BE RELIANT SOLELY UPON ONE FORM OF ELECTRONIC SYSTEM DUE TO THE RISK OF SYSTEM FAILURE. TO AVOID  
DAMAGE, INJURY, OR DEATH, THE USER OR APPLICATION DESIGNER MUST TAKE REASONABLY PRUDENT STEPS TO  
PROTECT AGAINST SYSTEM FAILURES, INCLUDING BUT NOT LIMITED TO BACK-UP OR SHUT DOWN MECHANISMS.  
BECAUSE EACH END-USER SYSTEM IS CUSTOMIZED AND DIFFERS FROM NATIONAL INSTRUMENTS' TESTING  
PLATFORMS AND BECAUSE A USER OR APPLICATION DESIGNER MAY USE NATIONAL INSTRUMENTS PRODUCTS IN  
COMBINATION WITH OTHER PRODUCTS IN A MANNER NOT EVALUATED OR CONTEMPLATED BY NATIONAL  
INSTRUMENTS, THE USER OR APPLICATION DESIGNER IS ULTIMATELY RESPONSIBLE FOR VERIFYING AND VALIDATING  
THE SUITABILITY OF NATIONAL INSTRUMENTS PRODUCTS WHENEVER NATIONAL INSTRUMENTS PRODUCTS ARE  
INCORPORATED IN A SYSTEM OR APPLICATION, INCLUDING, WITHOUT LIMITATION, THE APPROPRIATE DESIGN,  
PROCESS AND SAFETY LEVEL OF SUCH SYSTEM OR APPLICATION.  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Conventions  
The following conventions are used in this manual:  
[ ]  
Square brackets enclose optional items—for example, [response]. Square  
brackets also cite bibliographic references.  
»
The » symbol leads you through nested menu items and dialog box options  
to a final action. The sequence File»Page Setup»Options directs you to  
pull down the File menu, select the Page Setup item, and select Options  
from the last dialog box.  
This icon denotes a note, which alerts you to important information.  
bold  
Bold text denotes items that you must select or click in the software, such  
as menu items and dialog box options. Bold text also denotes parameter  
names.  
italic  
Italic text denotes variables, emphasis, a cross reference, or an introduction  
to a key concept. This font also denotes text that is a placeholder for a word  
or value that you must supply.  
monospace  
Text in this font denotes text or characters that you should enter from the  
keyboard, sections of code, programming examples, and syntax examples.  
This font is also used for the proper names of disk drives, paths, directories,  
programs, subprograms, subroutines, device names, functions, operations,  
variables, filenames, and extensions.  
monospace bold  
Bold text in this font denotes the messages and responses that the computer  
automatically prints to the screen. This font also emphasizes lines of code  
that are different from the other examples.  
monospace italic  
Italic text in this font denotes text that is a placeholder for a word or value  
that you must supply.  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 1  
Introduction  
Controllability and Observability Grammians ................................................1-7  
Hankel Singular Values...................................................................................1-8  
Singular Perturbation.......................................................................................1-11  
Chapter 2  
Algorithm ........................................................................................................2-12  
Related Functions ............................................................................................2-14  
ophank( )........................................................................................................................2-14  
Restriction........................................................................................................2-14  
Algorithm ........................................................................................................2-15  
Behaviors.........................................................................................................2-15  
© National Instruments Corporation  
v
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Contents  
Onepass Algorithm ......................................................................................... 2-18  
Multipass Algorithm ...................................................................................... 2-20  
Impulse Response Error.................................................................................. 2-22  
Chapter 3  
Hankel Singular Values of Phase Matrix of Gr............................................... 3-9  
Algorithm........................................................................................................ 3-14  
right and left.................................................................................................... 3-15  
Error Bounds................................................................................................... 3-20  
Chapter 4  
wtbalance( ) ................................................................................................................... 4-10  
Algorithm........................................................................................................ 4-12  
Related Functions............................................................................................ 4-15  
fracred( )........................................................................................................................ 4-15  
Restrictions...................................................................................................... 4-15  
Defining and Reducing a Controller ............................................................... 4-16  
Xmath Model Reduction Module  
vi  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Contents  
Algorithm ........................................................................................................4-18  
Chapter 5  
Utilities  
hankelsv( )......................................................................................................................5-1  
Related Functions ............................................................................................5-2  
Algorithm ........................................................................................................5-2  
Chapter 6  
Tutorial  
Plant and Full-Order Controller.....................................................................................6-1  
ophank( )..........................................................................................................6-9  
wtbalance.........................................................................................................6-12  
fracred..............................................................................................................6-20  
Appendix A  
Bibliography  
Appendix B  
Technical Support and Professional Services  
Index  
© National Instruments Corporation  
vii  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
1
Introduction  
This chapter starts with an outline of the manual and some useful notes. It  
also provides an overview of the Model Reduction Module, describes the  
functions in this module, and introduces nomenclature and concepts used  
throughout this manual.  
Using This Manual  
This manual describes the Model Reduction Module (MRM), which  
provides a collection of tools for reducing the order of systems.  
Readers who are not familiar with Parameter Dependent Matrices (PDMs)  
should consult the Xmath User Guide before using MRM functions and  
tools. Although several MRM functions accept both PDMs and matrices as  
input parameters, PDMs are preferable because they can include additional  
information that is useful for simulation, plotting, and signal labeling.  
Document Organization  
This manual includes the following chapters:  
Chapter 1, Introduction, starts with an outline of the manual and some  
useful notes. It also provides an overview of the Model Reduction  
Module, describes the functions in this module, and introduces  
nomenclature and concepts used throughout this manual.  
Chapter 2, Additive Error Reduction, describes additive error  
and perturbation of balanced realizations.  
Chapter 3, Multiplicative Error Reduction, describes multiplicative  
error reduction presenting considerations for using multiplicative  
rather than additive error reduction.  
Chapter 4, Frequency-Weighted Error Reduction, describes  
frequency-weighted error reduction problems, including controller  
reduction and fractional representations.  
© National Instruments Corporation  
1-1  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
             
Chapter 1  
Introduction  
Chapter 5, Utilities, describes three utility functions: hankelsv( ),  
stable( ), and compare( ).  
Chapter 6, Tutorial, illustrates a number of the MRM functions and  
their underlying ideas.  
Bibliographic References  
Throughout this document, bibliographic references are cited with  
bracketed entries. For example, a reference to [VODM1] corresponds  
to a paper published by Van Overschee and De Moor. For a table of  
bibliographic references, refer to Appendix A, Bibliography.  
Commonly Used Nomenclature  
This manual uses the following general nomenclature:  
Matrix variables are generally denoted with capital letters; vectors are  
represented in lowercase.  
G(s) is used to denote a transfer function of a system where s is the  
Laplace variable. G(q) is used when both continuous and discrete  
systems are allowed.  
H(s) is used to denote the frequency response, over some range of  
frequencies of a system where s is the Laplace variable. H(q) is used  
to indicate that the system can be continuous or discrete.  
A single apostrophe following a matrix variable, for example, x’,  
denotes the transpose of that variable. An asterisk following a matrix  
variable, for example, A*, indicates the complex conjugate, or  
Hermitian, transpose of that variable.  
Conventions  
This publication makes use of the following types of conventions: font,  
format, symbol, mouse, and note. These conventions are detailed in  
Chapter 2, MATRIXx Publications, Online Help, and Customer Support,  
of the MATRIXx Getting Started Guide.  
Xmath Model Reduction Module  
1-2  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 1  
Introduction  
Related Publications  
For a complete list of MATRIXx publications, refer to Chapter 2,  
MATRIXx Publications, Online Help, and Customer Support, of the  
MATRIXx Getting Started Guide. The following documents are particularly  
useful for topics covered in this manual:  
MATRIXx Getting Started Guide  
Xmath User Guide  
Control Design Module  
Interactive Control Design Module  
Interactive System Identification Module, Part 1  
Interactive System Identification Module, Part 2  
Model Reduction Module  
Optimization Module  
Robust Control Module  
Xµ Module  
MATRIXx Help  
Model Reduction Module function reference information is available in  
the MATRIXx Help. The MATRIXx Help includes all Model Reduction  
functions. Each topic explains a function’s inputs, outputs, and keywords  
in detail. Refer to Chapter 2, MATRIXx Publications, Online Help, and  
Customer Support, of the MATRIXx Getting Started Guide for complete  
instructions on using the help feature.  
Overview  
The Xmath Model Reduction Module (MRM) provides a collection of tools  
for reducing the order of systems. Many of the functions are based on  
state-of-the-art algorithms in conjunction with researchers at the Australian  
National University, who were responsible for the original development of  
some of the algorithms. A common theme throughout the module is the use  
of Hankel singular values and balanced realizations, although  
considerations of numerical accuracy often dictates whether these tools are  
used implicitly rather than explicitly. The tools are particularly suitable  
when, as generally here, quality of approximation is measured by closeness  
of frequency domain behavior.  
© National Instruments Corporation  
1-3  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 1  
Introduction  
As shown in Figure 1-1, functions are provided to handle four broad tasks:  
Model reduction with additive errors  
Model reduction with multiplicative errors  
Model reduction with frequency weighting of an additive error,  
including controller reduction  
Utility functions  
Functions  
Additive Error  
Model Reduction  
Multiplicative  
Model Reduction  
Frequency Weighted  
Model Reduction  
balmoore  
redschur  
ophank  
truncate  
balance  
mreduce  
bst  
mulhank  
wtbalance  
fracred  
Utility Functions  
hankelsv  
stable  
compare  
Figure 1-1. MRM Function  
Xmath Model Reduction Module  
1-4  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 1  
Introduction  
Certain restrictions regarding minimality and stability are required of the  
input data, and are summarized in Table 1-1.  
Table 1-1. MRM Restrictions  
balance( )  
A stable, minimal system  
balmoore ( )  
A state-space system must be stable and minimal,  
having at least one input, output, and state  
bst( )  
A state-space system must be linear,  
continuous-time, and stable, with full rank along  
the jω-axis, including infinity  
compare( )  
fracred( )  
Must be a state-space system  
A state-space system must be linear and continuous  
A system must be linear and stable  
hankelsv( )  
mreduce( )  
A submatrix of a matrix must be nonsingular  
for continuous systems, and variant for discrete  
systems  
mulhank( )  
A state-space system must be linear,  
continuous-time, stable and square, with full  
rank along the jω-axis, including infinity  
ophank( )  
A state-space system must be linear,  
continuous-time and stable, but can be nonminimal  
redschur( )  
A state-space system must be stable and linear,  
but can be nonminimal  
stable ( )  
truncate( )  
wtbalance( )  
No restriction  
Any full-order state-space system  
A state-space system must be linear and  
continuous. Interconnection of controller and plant  
must be stable, and/or weight must be stable.  
Documentation of the individual functions sometimes indicates how the  
restrictions can be circumvented. There are a number of model reduction  
methods not covered here. These include:  
Padé Approximation  
Methods based on interpolating, or matching at discrete frequencies  
© National Instruments Corporation  
1-5  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
                               
Chapter 1  
Introduction  
L2 approximation, in which the L2 norm of impulse response error (or,  
by Parseval’s theorem, the L2 norm of the transfer-function error along  
the imaginary axis) serves as the error measure  
Markov parameter or impulse response matching, moment matching,  
covariance matching, and combinations of these, for example,  
q-COVER approximation  
Controller reduction using canonical interactions, balanced Riccati  
equations, and certain balanced controller reduction algorithms  
Nomenclature  
This manual uses standard nomenclature. The user should be familiar with  
the following:  
sup denotes supremum, the least upper bound.  
The acute accent (´) denotes matrix transposition.  
A superscripted asterisk (*) denotes matrix transposition and complex  
conjugation.  
λmax(A) for a square matrix A denotes the maximum eigenvalue,  
presuming there are no complex eigenvalues.  
Reλi(A) and |λi(A)| for a square matrix A denote an arbitrary real part  
and an arbitrary magnitude of an eigenvalue of A.  
X( jω) for a transfer function X(·) denotes:  
sup  
max[X*(jω)X(jω)]]1/2  
∞ < ω < ∞  
An all-pass transfer-function W(s) is one where X(jω) = 1 for all ω;  
to each pole, there corresponds a zero which is the reflection through  
the jω-axis of the pole, and there are no jω-axis poles.  
An all-pass transfer-function matrix W(s) is a square matrix where  
W′(jω)W( jω) = I  
P > 0 and P 0 for a symmetric or hermitian matrix denote positive  
and nonnegative definiteness.  
P1 > P2 and P1 P2 for symmetric or hermitian P1 and P2 denote  
P1 P2 is positive definite and nonnegative definite.  
A superscripted number sign (#) for a square matrix A denotes the  
Moore-Penrose pseudo-inverse of A.  
Xmath Model Reduction Module  
1-6  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
                 
Chapter 1  
Introduction  
An inequality or bound is tight if it can be met in practice, for example  
1 + logx x 0  
is tight because the inequality becomes an equality for x = 1. Again,  
if F(jω) denotes the Fourier transform of some f (t)∈ L2 , the  
Heisenberg inequality states,  
f(t)2dt  
------------------------------------------------------------------------------------------  
1 2 4π  
1 2  
t2 f(t) 2dt  
ω2 F(jω) 2dω  
and the bound is tight since it is attained for f(t) = exp + (–kt2).  
Commonly Used Concepts  
This section outlines some frequently used standard concepts.  
Controllability and Observability Grammians  
Suppose that G(s) = D + C(sIA)–1B is a transfer-function matrix with  
Reλi(A)<0. Then there exist symmetric matrices P, Q defined by:  
PA+ AP = –BB′  
QA + AQ = –CC  
These are termed the controllability and observability grammians of the  
realization defined by {A,B,C,D}. (Sometimes in the code, WC is used for  
P and WO for Q.) They have a number of properties:  
P 0, with P > 0 if and only if [A,B] is controllable, Q 0 with Q > 0  
if and only if [A,C] is observable.  
P = eAtBBeAtdt and Q = eAtCCeAtdt  
0
0
With vec P denoting the column vector formed by stacking column 1  
of P on column 2 on column 3, and so on, and denoting Kronecker  
product  
[I A + A I]vecP = –vec(BB)  
The controllability grammian can be thought of as measuring the  
difficulty of controlling a system. More specifically, if the system is in  
the zero state initially, the minimum energy (as measured by the L2  
norm of u) required to bring it to the state x0 is x0P –1x0; so small  
eigenvalues of P correspond to systems that are difficult to control,  
while zero eigenvalues correspond to uncontrollable systems.  
© National Instruments Corporation  
1-7  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
         
Chapter 1  
Introduction  
The controllability grammian is also E[x(t)x(t)] when the system  
·
x = Ax + Bw has been excited from time by zero mean white  
noise with E[w(t)w′(s)] = Iδ(t s).  
The observability grammian can be thought of as measuring the  
information contained in the output concerning an initial state.  
·
If x = Ax, y = Cx with x(0) = x0 then:  
y′(t)y(t)dt = x0Qx0  
0
Systems that are easy to observe correspond to Q with large  
eigenvalues, and thus large output energy (when unforced).  
lyapunov(A,B*B')produces P and lyapunov(A',C'*C)  
produces Q.  
For a discrete-time G(z) = D + C(zI-A)–1B with |λi(A)|<1, P and Q are:  
P – APA= BB′  
Q – AQA = CC  
The first dot point above remains valid. Also,  
P =  
AkBBAk and Q =  
AkCCAk  
k = 0  
k = 0  
with the sums being finite in case A is nilpotent (which is the case if  
the transfer-function matrix has a finite impulse response).  
[IAA] vec P = vec (BB′)  
lyapunov( )can be used to evaluate P and Q.  
Hankel Singular Values  
If P, Q are the controllability and observability grammians of a  
transfer-function matrix (in continuous or discrete time), the Hankel  
Singular Values are the quantities λi1/2(PQ). Notice the following:  
All eigenvalues of PQ are nonnegative, and so are the Hankel singular  
values.  
The Hankel singular values are independent of the realization used to  
calculate them: when A,B,C,D are replaced by TAT–1, TB, CT–1 and D,  
then P and Q are replaced by TPTand (T–1)QT–1; then PQ is replaced  
by TPQT–1 and the eigenvalues are unaltered.  
The number of nonzero Hankel singular values is the order or  
McMillan degree of the transfer-function matrix, or the state  
dimension in a minimal realization.  
Xmath Model Reduction Module  
1-8  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 1  
Introduction  
Suppose the transfer-function matrix corresponds to a discrete-time  
system, with state variable dimension n. Then the infinite Hankel  
matrix,  
CB CAB CA2B  
CAB CA2B  
CA2B  
H =  
has for its singular values the n nonzero Hankel singular values,  
together with an infinite number of zero singular values.  
The Hankel singular values of a (stable) all pass system (or all pass matrix)  
are all 1.  
Slightly different procedures are used for calculating the Hankel singular  
values (and so-called weighted Hankel singular values) in the various  
functions. These procedures are summarized in Table 1-2.  
Table 1-2. Calculating Hankel Singular Values  
(balance( ))  
the Internally Balanced Realizations section; the  
Hankel singular values are given by  
diag(R1/2) = HSV  
balmoore( )  
For a discussion of the balancing algorithm, refer to  
the Internally Balanced Realizations section; the  
matrix SH yields the Hankel singular values through  
diag(SH)  
hankelsv( )  
ophank( )  
real(sqrt(eig(p*q)))  
Calls hankelsv( )  
redschur( )  
Computes a Schur decomposition of P*Q and then  
takes the square roots of the diagonal entries  
bst( )  
Same as redschur( )except either P or Q can be  
mulhank( )  
wtbalance( )  
fracred( )  
a weighted grammian  
© National Instruments Corporation  
1-9  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 1  
Introduction  
Internally Balanced Realizations  
Suppose that a realization of a transfer-function matrix has the  
controllability and observability grammian property that P = Q = Σ for  
some diagonal Σ. Then the realization is termed internally balanced. Notice  
that the diagonal entries σi of Σ are square roots of the eigenvalues of PQ,  
that is, they are the Hankel singular values. Often the entries of Σ are  
assumed ordered with σi ≥ σi+1  
.
As noted in the discussion of grammians, systems with small (eigenvalues  
of) P are hard to control and those with small (eigenvalues of) Q are hard  
to observe. Now a state transformation T = α I will cause P = Q to be  
replaced by α2P, α–2Q, implying that ease of control can be obtained at the  
expense of difficulty of observation, and conversely. Balanced realizations  
are those when ease of control has been balanced against ease of  
observation.  
Given an arbitrary realization, there are a number of ways of finding a  
state-variable coordinate transformation bringing it to balanced form.  
A good survey of the available algorithms for balancing is in [LHPW87].  
One of these is implemented in the Xmath function balance( ).  
The one implemented in balmoore( )as part of this module is more  
sophisticated, but more time consuming. It proceeds as follows:  
1. Singular value decompositions of P and Q are defined. Because P and  
Q are symmetric, this is equivalent to diagonalizing P and Q by  
orthogonal transformations.  
P = UcSc U′  
c
Q = UoSo Uo  
2. The matrix,  
H = S10 2UHSHV1H2  
is constructed, and from it, a singular value decomposition is obtained:  
H = UHSHVH  
3. The balancing transformation is given by:  
T = U0S01 2UHSH1 2  
The balanced realization is T–1AT, T–1B, CT.  
Xmath Model Reduction Module  
1-10  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
Chapter 1  
Introduction  
This is almost the algorithm set out in Section II of [LHPW87]. The one  
difference (and it is minor) is that in [LHPW87], lower triangular Cholesky  
factors of P and Q are used, in place of UcSc1/2 and UOSO1/2 in forming H  
in step 2. The grammians themselves need never be formed, as these  
Cholesky factors can be computed directly from A, B, and C in both  
continuous and discrete time; this, however, is not done in balmoore.  
The algorithm has the property that:  
TQT= T1P(T1)′ = SH  
Thus the diagonal entries of SH are the Hankel singular values.  
The algorithm implemented in balance( )is older, refer to [Lau80].  
A lower triangular Cholesky factor Lc of P is found, so that LcLc= P.  
Then the symmetric matrix LcQLc is diagonalized (through a singular  
value decomposition), thus LcQLc = VRU, with actually V = U. Finally,  
the coordinate basis transformation is given by T = LcVR–1/4, resulting in  
TQT = T–1P(T–1)= R1/2  
.
Singular Perturbation  
A common procedure for approximating systems is the technique of  
Singular Perturbation. The underlying heuristic reasoning is as follows.  
Suppose there is a system with the property that the unforced responses of  
the system contain some modes which decay to zero extremely fast. Then  
an approximation to the system behavior may be attained by setting state  
variable derivatives associated with these modes to zero, even in the forced  
case. The phrase “associated with the modes” is loose: exactly what occurs  
is shown below. The phrase “even in the forced case” captures a logical  
flaw in the argument: smallness in the unforced case associated with initial  
conditions is not identical with smallness in the forced case.  
Suppose the system is defined by:  
·
x1  
x1  
x2  
B1  
B2  
A11 A12  
A21 A22  
=
+
u
·
x2  
x1  
y =  
+ Du  
C1 C2  
x2  
© National Instruments Corporation  
1-11  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 1  
Introduction  
and also:  
1  
22  
Reλi(A22)<0 and Reλi(A11 A12A A21) < 0.  
Usually, we expect that,  
1  
22  
Reλi(A22) « Reλi(A11 A12A A21)  
in the sense that the intuitive argument hinges on this, but it is not necessary.  
·
Then a singular perturbation is obtained by replacing x2 by zero; this  
means that:  
1  
22  
1  
22  
A21x1 + A22x2 + B2u = 0  
Accordingly,  
x2 = –A A21x1 A B2u  
or  
1  
1  
·
x1 = (A11 = A12A22A21)x1 + (B1 A12A22B2)u  
1  
22  
1  
22  
y = (C1 C2A A21)x1 + (D C2A B2)u  
(1-2)  
Equation 1-2 may be an approximation for Equation 1-1. This means that:  
The transfer-function matrices may be similar.  
If Equation 1-2 is excited by some u(·), with initial condition x1(to), and  
if Equation 1-1 is excited by the same u(·) with initial condition given  
by,  
x1(to) and x2(to) = –A–122A21x1(to) –A22 B2u(to),  
then x1(·) and y(·) computed from Equation 1-1 and from Equation 1-2  
should be similar.  
If Equation 1-1 and Equation 1-2 are excited with the same u(·), have  
the same x1(to) and Equation 1-1 has arbitrary x2, then x1(·) and y(·)  
computed from Equation 1-1 and Equation 1-2 should be similar after  
a possible initial transient.  
As far as the transfer-function matrices are concerned, it can be verified that  
they are actually equal at DC.  
Xmath Model Reduction Module  
1-12  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 1  
Introduction  
Similar considerations govern the discrete-time problem, where,  
x1(k + 1)  
x1(k)  
B1  
B2  
A11 A12  
A21 A22  
=
+
u(k)  
x2(k + 1)  
x2(k)  
x1(k)  
y(k) =  
+ Du(k)  
C1 C2  
x2(k)  
can be approximated by:  
x1(k + 1) = [A11 + A12(I A22)1A21]x1(k) +  
[B1 + A12(I A22)1B2]u(k)  
yk = [C1 + C2(I A22)1A21]x1(k) +  
[D + C2(I A22)1B2]u(k)  
mreduce( )can carry out singular perturbation. For further discussion,  
refer to Chapter 2, Additive Error Reduction. If Equation 1-1 is balanced,  
singular perturbation is provably attractive.  
Spectral Factorization  
Let W(s) be a stable transfer-function matrix, and suppose a system S with  
transfer-function matrix W(s) is excited by zero mean unit intensity white  
noise. Then the output of S is a stationary process with a spectrum Φ(s)  
related to W(s) by:  
Φ(s) = W(s)W′(s)  
(1-3)  
Evidently,  
Φ( jω) = W( jω)W*( jω)  
so that Φ( jω) is nonnegative hermitian for all ω; when W( jω) is a scalar, so  
is Φ( jω) with Φ( jω) = |W( jω)|2.  
In the matrix case, Φ is singular for some ω only if W does not have full  
rank there, and in the scalar case only if W has a zero there.  
Spectral factorization, as shown in Example 1-1, seeks a W(jω), given  
Φ(jω). In the rational case, a W(jω) exists if and only if Φ(jω) is  
© National Instruments Corporation  
1-13  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 1  
Introduction  
nonnegative hermitian for all ω. If Φ is scalar, then Φ(jω)0 for all ω.  
Normally one restricts attention to Φ(·) with limω→∞Φ(jω)<. A key result  
is that, given a rational, nonnegative hermitian Φ(jω) with  
limω→∞Φ(jω)<, there exists a rational W(s) where,  
W()<∞.  
W(s) is stable.  
W(s) is minimum phase, that is, the rank of W(s) is constant in Re[s]>0.  
In the scalar case, all zeros of W(s) lie in Re[s]0, or in Re[s]<0 if Φ(jω)>0  
for all ω.  
In the matrix case, and if Φ(jω) is nonsingular for some ω, it means that  
W(s) is square and W–1(s) has all its poles in Re[s]0, or in Re[s]<0 if Φ(jω)  
is nonsingular for all ω.  
Moreover, the particular W(s) previously defined is unique, to within right  
multiplication by a constant orthogonal matrix. In the scalar case, this  
means that W(s) is determined to within a 1 multiplier.  
Example 1-1  
Example of Spectral Factorization  
Suppose:  
ω2 + 1  
Φ( jω) = ---------------  
ω2 + 4  
s + 1  
-----------  
Then Equation 1-3 is satisfied by W(s) =  
minimum phase.  
, which is stable and  
s + 2  
s 1  
s + 2  
s 3 s 1  
s + 2 s + 2  
sTs + 1  
-----------  
----------- -----------  
,
-----------  
e
Also, Equation 1-3 is satisfied by  
and  
and  
, and  
s + 2  
so forth, but none of these is minimum phase.  
bst( )and mulhank( )both require execution within the program of  
a spectral factorization; the actual algorithm for achieving the spectral  
factorization depends on a Riccati equation. The concepts of a spectrum  
and spectral factor also underpin aspects of wtbalance( ).  
Xmath Model Reduction Module  
1-14  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
Chapter 1  
Introduction  
Low Order Controller Design Through Order Reduction  
The Model Reduction Module is particularly suitable for achieving low  
order controller design for a high order plant. This section explains some of  
the broad issues involved.  
Most modern controller design methods, for example, LQG and H, yield  
controllers of order roughly comparable with that of the plant. It follows  
that, to obtain a low order controller using such methods, one must either  
follow a high order controller design by a controller reduction step,  
or reduce an initially given high order plant model, and then design a  
controller using the resulting low order plant, with the understanding that  
the controller will actually be used on the high order plant. Refer to  
Figure 1-2.  
High Order Plant  
High Order Controller  
Plant  
Reduction  
Controller  
Reduction  
Low Order Plant  
Low Order Controller  
Figure 1-2. Low Order Controller Design for a High Order Plant  
Generally speaking, in any design procedure, it is better to postpone  
approximation to a late step of the procedure: if approximation is done  
early, the subsequent steps of the design procedure may have unpredictable  
effects on the approximation errors. Hence, the scheme based on high order  
Controller reduction should aim to preserve closed-loop properties as far  
as possible. Hence the controller reduction procedures advocated in this  
module reflect the plant in some way. This leads to the frequency weighted  
reduction schemes of wtbalance( )and fracred( ), as described in  
Chapter 4, Frequency-Weighted Error Reduction. Plant reduction logically  
should also seek to preserve closed-loop properties, and thus should involve  
the controller. With the controller unknown however, this is impossible.  
Nevertheless, it can be argued, on the basis of the high loop gain property  
within the closed-loop bandwidth that is typical of many systems, that  
© National Instruments Corporation  
1-15  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 1  
Introduction  
multiplicative reduction, as described in Chapter 4, Frequency-Weighted  
Error Reduction, is a sound approach. Chapter 3, Multiplicative Error  
Reduction, and Chapter 4, Frequency-Weighted Error Reduction, develop  
these arguments more fully.  
Xmath Model Reduction Module  
1-16  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
2
Additive Error Reduction  
This chapter describes additive error reduction including discussions of  
truncation of, reduction by, and perturbation of balanced realizations.  
Introduction  
Additive error reduction focuses on errors of the form,  
G( jω) Gr( jω)  
where G is the originally given transfer function, or model, and Gr is the  
reduced one. Of course, in discrete-time, one works instead with:  
G(ejω) Gr(ejω)  
As is argued in later chapters, if one is reducing a plant that will sit inside  
a closed loop, or if one is reducing a controller, that again is sitting in a  
closed loop, focus on additive error model reduction may not be  
appropriate. It is, however, extremely appropriate in considering reducing  
the transfer function of a filter. One pertinent application comes specifically  
from digital filtering: a great many design algorithms lead to a finite  
impulse response (FIR) filter which can have a very large number of  
coefficients when poles are close to the unit circle. Model reduction  
provides a means to replace an FIR design by a much lower order infinite  
impulse response (IIR) design, with close matching of the transfer function  
at all frequencies.  
© National Instruments Corporation  
2-1  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
         
Chapter 2  
Additive Error Reduction  
Truncation of Balanced Realizations  
A group of functions can be used to achieve a reduction through truncation  
of a balanced realization. This means that if the original system is  
·
x1  
x1  
x2  
B1  
B2  
A11 A12  
A21 A22  
=
+
u
·
x2  
y =  
x + Du  
C1 C2  
(2-1)  
and the realization is internally balanced, then a truncation is provided by  
·
x1 = A11x1 + B1u  
y = C1x1 + Du  
The functions in question are:  
• balmoore( )  
balance( )(refer to the Xmath Help)  
• truncate( )  
• redschur( )  
One only can speak of internally balanced realizations for systems which  
are stable; if the aim is to reduce a transfer function matrix G(s) which  
contains unstable poles, one must additively decompose it into a stable part  
and unstable part, reduce the stable part, and then add the unstable part back  
in. The function stable( ), described in Chapter 5, Utilities, can be used  
to decompose G(s). Thus:  
G(s)  
=
=
=
Gs(s) + Gu(s)(Gs(s) stable, Gu(s) unstable)  
found by algorithm (reduction of Gs(s))  
Gsr(s) + Gu(s) (reduction of G(s))  
Gsr(s)  
Gr(s)  
Xmath Model Reduction Module  
2-2  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
Chapter 2  
Additive Error Reduction  
A very attractive feature of the truncation procedure is the availability  
of an error bound. More precisely, suppose that the controllability and  
observability grammians for [Enn84] are  
Σ1  
0
P = Q = Σ =  
(2-2)  
0 Σ2  
with the diagonal entries of Σ in decreasing order, that is, σ1 ≥ σ2 ···. Then  
the key result is,  
G( jω) Gr( jω) ≤ 2trΣ2  
with G, Gr the transfer function matrices of Equation 2-1 and Equation 2-2,  
respectively. This formula shows that small singular values can, without  
great cost, be thrown away. It also is valid in discrete time, and can be  
improved upon if there are repeated Hankel singular values. Provided that  
the smallest diagonal entry of Σ1 strictly exceeds the largest diagonal entry  
of Σ2, the reduced order system is guaranteed to be stable.  
Several other points concerning the error can be made:  
The error G( jω)–Gr( jω) as a function of frequency is not flat; it is zero  
at ω = , and may take its largest value at ω = 0, so that there is in  
general no matching of DC gains of the original and reduced system.  
The actual error may be considerably less than the error bound at all  
frequencies, so that the error bound formula can be no more than an  
advance guide. However, the bound is tight when the dimension  
reduction is 1 and the reduction is of a continuous-time  
transfer-function matrix.  
With g(·) and gr(·) denoting the impulse responses for impulse  
responses of G and Gr and with Gr of degree k, the following L1 bound  
holds [GCP88]  
g gr ≤ (4 2k + 1)trΣ2  
1
This bound also will apply for the Lerror on the step response.  
It is helpful to note one situation where reduction is likely to be difficult (so  
that Σ will contain few diagonal entries which are, relatively, very small).  
Suppose G(s), strictly proper, has degree n and has (n – 1) unstable zeros.  
Then as ω runs from zero to infinity, the phase of G(s) will change by  
(2n – 1)π/2. Much of this change may occur in the passband. Suppose Gr(s)  
has degree n–1; it can have no more than (n – 2) zeros, since it is strictly  
© National Instruments Corporation  
2-3  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 2  
Additive Error Reduction  
proper. So, even if all zeros are unstable, the maximum phase shift when ω  
moves from 0 to is (2n – 3)π/2. It follows that if G(jω) remains large in  
magnitude at frequencies when the phase shift has moved past (2n – 3)π/2,  
approximation of G by Gr will necessarily be poor. Put another way, good  
approximation may depend somehow on removing roughly cancelling  
pole-zeros pairs; when there are no left half plane zeros, there can be no  
rough cancellation, and so approximation is unsatisfactory.  
As a working rule of thumb, if there are p right half plane zeros in the  
passband of a strictly proper G(s), reduction to a Gr(s) of order less than  
p + 1 is likely to involve substantial errors. For non-strictly proper G(s),  
having p right half plane zeros means that reduction to a Gr(s) of order less  
than p is likely to involve substantial errors.  
An all-pass function exemplifies the problem: there are n stable poles and  
n unstable zeros. Since all singular values are 1, the error bound formula  
indicates for a reduction to order n – 1 (when it is not just a bound, but  
exact) a maximum error of 2.  
Another situation where poor approximation can arise is when a highly  
oscillatory system is to be replaced by a system with a real pole.  
Reduction Through Balanced Realization Truncation  
This section briefly describes functions that reduce( ), balance( ),  
and truncate( )to achieve reduction.  
balmoore( )Computes an internally balanced realization of a  
system and optionally truncates the realization to form an  
approximation.  
balance( )Computes an internally balanced realization of a  
system.  
truncate( )This function truncates a system. It allows  
examination of a sequence of different reduced order models formed  
from the one balanced realization.  
redschur( )These functions in theory function almost the same  
as the two features of balmoore( ). That is, they produce a  
state-variable realization of a reduced order model, such that the  
transfer function matrix of the model could have resulted by truncating  
a balanced realization of the original full order transfer function  
matrix. However, the initially given realization of the original transfer  
function matrix is never actually balanced, which can be a numerically  
hazardous step. Moreover, the state-variable realization of the reduced  
Xmath Model Reduction Module  
2-4  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
               
Chapter 2  
Additive Error Reduction  
order model is not one in general obtainable by truncation of an  
internally-balanced realization of the full order model.  
Figure 2-1 sets out several routes to a reduced-order realization. In  
continuous time, a truncation of a balanced realization is again balanced.  
This is not the case for discrete time, but otherwise it looks the same.  
Full Order Realization  
balmoore  
(with both steps)  
balmoore  
(with first step)  
balance  
redschur  
truncate  
Balanced Realization of  
Nonbalanced  
Reduced Order Model  
(in continuous time)  
Realization of  
Reduced Order Model  
Reduced Order Model Transfer Function  
Figure 2-1. Different Approaches for Obtaining the Same Reduced Order Model  
Singular Perturbation of Balanced Realization  
Singular perturbation of a balanced realization is an attractive way to  
produce a reduced order model. Suppose G(s) is defined by,  
·
x1  
x1  
x2  
B1  
B2  
A11 A12  
A21 A22  
=
+
u
·
x2  
y =  
x + Du  
C1 C2  
© National Instruments Corporation  
2-5  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 2  
Additive Error Reduction  
with controllability and observability grammians given by,  
Σ1  
0
P = Q = Σ =  
0 Σ2  
in which the diagonal entries of Σ are in decreasing order, that is,  
σ1 ≥ σ2 ···, and such that the last diagonal entry of Σ1 exceeds  
the first diagonal entry of Σ2. It turns out that Reλi(A1)<0 and  
22  
Reλ (A11A12 A1A21)< 0, and a reduced order model Gr(s) can be  
22  
i
defined by:  
1  
1  
·
x = (A11 A12A22A21)x + (B1 + –A12A22B2)u  
1  
22  
1  
22  
y = (C1 C2A A21)x + (D C2A B2)u  
The attractive feature [LiA89] is that the same error bound holds as for  
balanced truncation. For example,  
G( jω) Gr( jω) ≤ 2trΣ2  
Although the error bounds are the same, the actual frequency pattern of  
the errors, and the actual maximum modulus, need not be the same for  
reduction to the same order. One crucial difference is that balanced  
truncation provides exact matching at ω = , but does not match at DC,  
while singular perturbation is exactly the other way round. Perfect  
matching at DC can be a substantial advantage, especially if input signals  
are known to be band-limited.  
Singular perturbation can be achieved with mreduce( ). Figure 2-1 shows  
the two alternative approaches. For both continuous-time and discrete-time  
reductions, the end result is a balanced realization.  
Hankel Norm Approximation  
In Hankel norm approximation, one relies on the fact that if one chooses an  
approximation to exactly minimize one norm (the Hankel norm) then the  
infinity norm will be approximately minimized. The Hankel norm is  
defined in the following way. Let G(s) be a (rational) stable transfer  
Xmath Model Reduction Module  
2-6  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
Chapter 2  
Additive Error Reduction  
function matrix. Consider the way the associated impulse response maps  
inputs defined over (–,0] in L2 into outputs, and focus on the output over  
[0,). Define the input as u(t) for t < 0, and set v(t) = u(–t). Define the  
output as y(t) for t > 0. Then the mapping is  
y(t) = CexpA(t + r)Bv(r)dr  
0
if G(s) = C(sI-A)–1B. The norm of the associated operator is the Hankel  
G
norm  
of G. A key result is that if σ1 σ2 ···, are the Hankel singular  
H
values of G(s), then G = σ1.  
H
To avoid minor confusion, suppose that all Hankel singular values of G are  
ˆ
distinct. Then consider approximating G by some stable G of prescribed  
ˆ
degree k much that  
is minimized. It turns out that  
G G  
H
ˆ
infGˆ of degree k G G = σk + 1(G)  
H
ˆ
and there is an algorithm available for obtaining G . Further, the  
ˆ
ˆ
optimum which is minimizing  
does a reasonable job  
G
G G  
H
ˆ
of minimizing  
, because it can be shown that  
G G  
ˆ
G G  
σj  
j = k + 1  
ˆ
where n = deg G, with this bound subject to the proviso that G and G are  
allowed to be nonzero and different at s = .  
ˆ
The bound on G G is one half that applying for balanced truncation.  
However,  
It is actual error that is important in practice (not bounds).  
The Hankel norm approximation does not give zero error at ω = ∞  
or at ω = 0. Balanced realization truncation gives zero error at ω = ,  
and singular perturbation of a balanced realization gives zero error  
at ω = 0.  
There is one further connection between optimum Hankel norm  
ˆ
approximation and Lerror. If one seeks to approximate G by a sum G + F,  
ˆ
with stable and of degree k and with F unstable, then:  
G
ˆ
infGˆ of degree k and F unstable G G F = σk + 1(G)  
© National Instruments Corporation  
2-7  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 2  
Additive Error Reduction  
ˆ
Further, the G which is optimal for Hankel norm approximation also is  
optimal for this second type of approximation.  
In Xmath Hankel norm approximation is achieved with ophank( ).  
The most comprehensive reference is [Glo84].  
balmoore( )  
[SysR,HSV,T] = balmoore(Sys,{nsr,bound})  
The balmoore( )function computes an internally-balanced realization of  
a continuous system and then optionally truncates it to provide a balance  
reduced order system using B.C. Moore’s algorithm.  
When balmoore( )is being used to reduce a system, its objective mirrors  
that of redschur( ), therefore, if the same Sysand nsrare used for both  
algorithms, the reduced order system should have the same transfer  
function (though in general the state-variable realizations will be different).  
When balmoore( )is being used to balance a system, its objective, like  
that of balance, is to generate an internally-balanced state-variable  
realization. The implementations are not identical.  
balmoore( )only can be applied on systems that have a stable A matrix,  
and are controllable and observable, (that is, minimal). Checks, which are  
rather time-consuming, are included. The computation is intrinsically not  
well-conditioned if Sysis nearly nonminimal. The first part of  
balmoore( )serves to find a transformation matrix T such that the  
controllability and observability grammians after transformation are equal,  
and diagonal, with decreasing entries down the diagonal, that is, the system  
representation is internally balanced. (The condition number of T is a  
measure of the ill-conditioning of the algorithm. If there is a problem with  
ill-conditioning, redschur( )can be used as an alternative.) If this  
common grammian is Σ, then after transformation:  
(continuous)  
Σ A+ A Σ = –BB′ Σ A + A′Σ = –CC  
(discrete)  
Σ – A Σ A = –BB′ Σ - A′ Σ A = –CC  
Σ = diag1, σ2, σ3..., σns] with σi ≤ σi + 1 > 0 with theσi the Hankel  
Singular Values of Sys. In the second part of balmoore( ), a truncation  
Xmath Model Reduction Module  
2-8  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 2  
Additive Error Reduction  
of the balanced system occurs, (assuming nsris less than the number of  
states). Thus, if the state-space representation of the balanced system is  
A11 A12  
A21 A22  
B1  
B2  
A =  
B =  
C =  
C1 C2  
with A11 possessing dimension nsr× nsr, B1 possessing nsrrows and C1  
possessing nsrcolumns, the reduced order system SysRis:  
(continuous)  
(discrete)  
·
x1(k + 1) = A11x1(k) + B1u(k)  
x1 = A11x1 + B1u  
y = C1x + Du  
y(k) = C1x1(k) + Du(k)  
The following error formula is relevant:  
(continuous)  
[C(jωI A)1] [C1(jωI A1)1(B1 + D)]  
2nsr + 1 + σnsr + 2 + ... + σns]  
(discrete)  
[C(ejωI A)1B + D] [C1(ejωI A1)1B1 + D]  
2nsr + 1 + σnsr + 2 + ... + σns]  
It is this error bound which is the basis of the determination of the order  
of the reduced system when the keyword boundis specified. If the error  
bound sought is smaller than 2σns, then no reduction is possible which is  
consistent with the error bound. If it is larger than 2trΣ, then the constant  
transfer function matrix D achieves the bound.  
For continuous systems, the actual approximation error depends on  
frequency, but is always zero at ω = . In practice it is often greatest at  
ω = 0; if the reduction of state dimension is 1, the error bound is exact, with  
the maximum error occurring at DC. The bound also is exact in the special  
case of a single-input, single-output transfer function which has poles and  
zero alternating along the negative real axis. It is far from exact when the  
poles and zeros approximately alternate along the imaginary axis (with the  
poles stable).  
© National Instruments Corporation  
2-9  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 2  
Additive Error Reduction  
The actual approximation error for discrete systems also depends on  
frequency, and can be large at ω = 0. The error bound is almost never tight,  
that is, the actual error magnitude as a function of ω almost never attains  
the error bound, so that the bound can only be a guide to the selection of the  
reduced system dimension.  
In principle, the error bound formula for both continuous and discrete  
systems can be improved (that is, made tighter or less likely to overestimate  
the actual maximum error magnitude) when singular values occur with  
multiplicity greater than one. However, because of errors arising in  
calculation, it is safer to proceed conservatively (that is, work with the error  
bound above) when using the error bound to select nsr, and examine the  
actual error achieved. If this is smaller than required, a smaller dimension  
for the reduced order system can be selected.  
mreduce( )provides an alternative reduction procedure for a balanced  
realization which achieves the same error bound, but which has zero error  
at ω = 0. For continuous systems there is generally some error at ω = ,  
because the D matrix is normally changed. (This means that normally the  
approximation of a strictly proper system through mreduce( )will not be  
strictly proper, in contrast to the situation with balmoore( ).) For discrete  
systems the D matrix is also normally changed so that, for example, a  
system which was strictly causal, or guaranteed to contain a delay (that is,  
D = 0), will be approximated by a system SysRwithout this property.  
The presentation of the Hankel singular values may suggest a logical  
dimension for the reduced order system; thus if σk » σk + 1, it may be  
sensible to choose nsr = k.  
With mreduce( )and a continuous system, the reduced order system  
SysRis internally balanced, with the grammian diag1, σ2, ...,σnsr], so  
that its Hankel Singular Values are a subset of those of the original system  
Sys. Provided σnsr > σnsr + 1, SysRalso is controllable, observable, and  
stable. This is not guaranteed if σnsr = σnsr + 1, so it is highly advisable to  
avoid this situation. Refer to the balmoore( ) section for more on the  
balmoore( )algorithm.  
With mreduce( )and discrete systems, the reduced order system SysRis  
not in general balanced (in contrast to balmoore( )), and its Hankel  
singular values are not in general a subset of those of Sys. Provided  
σnsr > σsrn + 1, the reduced order system SysRalso is controllable,  
observable and stable. This is not guaranteed if σnsr = σnsr + 1, so it is  
highly advisable to avoid this situation. For additional information about  
the balmoore( )function, refer to the Xmath Help.  
Xmath Model Reduction Module  
2-10  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 2  
Additive Error Reduction  
Related Functions  
balance(), truncate(), redschur(), mreduce()  
truncate( )  
SysR = truncate(Sys,nsr,{VD,VA})  
The truncate( )function reduces a system Sysby retaining the first  
nsrstates and throwing away the rest to form a system SysR.  
If for Sysone has,  
A11 A12  
A21 A22  
B1  
B2  
A =  
B =  
C =  
C1 C2  
the reduced order system (in both continuous-time and discrete-time cases)  
is defined by A11, B1, C1, and D. If Sysis balanced, then SysRis an  
approximation of Sysachieving a certain error bound. truncate( )may  
well be used after an initial application of balmoore( )to further reduce  
a system should a larger approximation error be tolerable. Alternatively, it  
may be used after an initial application of balance( )or redschur( ).  
If Syswas calculated from redschur( )and VA,VDwere posed as  
arguments, then SysRis calculated as in redschur( )(refer to the  
redschur( ) section).  
truncate( )should be contrasted with mreduce( ), which achieves a  
reduction through a singular perturbation calculation. If Sysis balanced,  
the same error bound formulas apply (though not necessarily the same  
errors), truncate( )always ensures exact matching at s = (in the  
continuous-time case), or exacting matching of the first impulse response  
coefficient D (in the discrete-time case), while mreduce( )ensures  
matching of DC gains for Sysand SysRin both the continuous-time and  
discrete-time case. For a additional information about the truncate( )  
function, refer to the Xmath Help.  
Related Functions  
balance(), balmoore(), redschur(), mreduce()  
© National Instruments Corporation  
2-11  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 2  
Additive Error Reduction  
redschur( )  
[SysR,HSV,slbig,srbig,VD,VA] = redschur(Sys,{nsr,bound})  
The redschur( )function uses a Schur method (from Safonov and  
Chiang) to calculate a reduced version of a continuous or discrete system  
without balancing.  
Algorithm  
The objective of redschur( )is the same as that of balmoore( )when  
the latter is being used to reduce a system; this means that if the same Sys  
and nsrare used for both algorithms, the reduced order system should have  
the same transfer function matrix. However, in contrast to balmoore( ),  
redschur( )do not initially transform Systo an internally balanced  
realization and then truncate it; nor is SysRin balanced form. The fact that  
there is no balancing offers numerical advantages, especially if Sysis  
nearly nonminimal.  
Sysshould be stable, and this is checked by the algorithm. In contrast to  
balmoore( ), minimality of Sys(that is, controllability and  
observability) is not required.  
If the Hankel singular values of Sysare ordered as σ1 ≥ σ2 ... ≥ σns 0,  
then those of SysRin the continuous-time case are σ1 ≥ σ2 ... ≥ σnsr > 0.  
A restriction of the algorithm is that σnsr > σnsr + 1 is required for both  
continuous-time and discrete-time cases. Under this restriction, SysRis  
guaranteed to be stable and minimal.  
The algorithms depend on the same algorithm, apart from the calculation  
of the controllability and observability grammians Wc and Wo of the  
original system. These are obtained as follows:  
(continuous)  
(discrete)  
Wc A+ AWc = –BB′  
WoA + AWo = CC  
Wc AWcA= BB′  
Wo AWoA = CC  
The maximum order permitted is the number of nonzero eigenvalues of  
WcWo that are larger than ε.  
Xmath Model Reduction Module  
2-12  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
Chapter 2  
Additive Error Reduction  
Next, Schur decompositions of WcWo are formed with the eigenvalues of  
WcWo in ascending and descending order. These eigenvalues are the square  
of the Hankel singular values of Sys, and if Sysis nonminimal, some can  
be zero.  
VAWcWoVA = Sasc  
VDWcWoVD = Sdes  
The matrices VA, VD are orthogonal and Sasc, Sdes are upper triangular. Next,  
submatrices are obtained as follows:  
0
Insr  
Vlbig = VA  
Vrbig = VD  
Insr  
0
and then a singular value decomposition is found:  
UebigSebigVebig = VlbigVrbig  
From these quantities, the transformation matrices used for calculating  
SysRare defined:  
1 2  
Slbig = VlbigUebig  
S
S
ebig  
1 2  
ebig  
Srbig = VrbigVebig  
and the reduced order system is:  
AR = SlbigASrbig  
BR = Slbig  
B
AR = CSrbig  
D
An error bound is available. In the continuous-time case it is as follows. Let  
G( jω) and GR( jω) be the transfer function matrices of Sysand SysR,  
respectively.  
For the continuous case:  
G( jω) GR( jω) ≤ 2nsr + 1 + σnsr + 2 + ... + σns)  
© National Instruments Corporation  
2-13  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 2  
Additive Error Reduction  
For the discrete-time case:  
G(ejω) GR(ejω) 2nsr + 1 + σnsr + 2 + ... + σns)  
When {bound}is specified, the error bound just enunciated is used to  
choose the number of states in SysRso that the bound is satisfied and nsr  
is as small as possible. If the desired error bound is smaller than 2σns,  
no reduction is made.  
In the continuous-time case, the error depends on frequency, but is always  
zero at ω = . If the reduction in dimension is 1, or the system Sysis  
single-input, single-output, with alternating poles and zeros on the real  
axis, the bound is tight. It is far from tight when the poles and zeros  
approximately alternate along the jω-axis. It is not normally tight in the  
discrete-time case, and for both continuous-time and discrete-time cases,  
it is not tight if there are repeated singular values.  
The presentation of the Hankel singular values may suggest a logical  
dimension for the reduced order system; thus if σk » σk + 1, it may be  
sensible to choose nsr = k.  
Related Functions  
ophank(), balmoore()  
ophank( )  
[SysR,SysU,HSV] = ophank(Sys,{nsr,onepass})  
The ophank( )function calculates an optimal Hankel norm reduction  
of Sys.  
Restriction  
This function has the following restriction:  
Only continuous systems are accepted; for discrete systems use  
makecontinuous( )before calling bst( ), then discretize the  
result.  
Sys=ophank(makecontinuous(SysD));  
SysD=discretize(Sys);  
Xmath Model Reduction Module  
2-14  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 2  
Additive Error Reduction  
Algorithm  
The algorithm does the following. The system Sysand the reduced order  
system SysRare stable; the system SysUhas all its poles in Re[s] > 0. If  
the transfer function matrices are G(s), Gr(s) and Gu(s) then:  
Gr(s) is a stable approximation of G(s).  
Gr(s) + Gu(s) is a more accurate, but not stable, approximation of G(s),  
and optimal in a certain sense.  
Of course, the algorithm works with state-space descriptions; that of G(s)  
can be minimal, while that of Gr(s) cannot be.  
These statements are explained in the Behaviors section. If onepassis  
specified, reduction is calculated in one pass. If onepassis not called or is  
set to 0 {onepass=0}, reduction is calculated in (number of states of  
Sys - nsr) passes. There seems to be no general rule to suggest which  
setting produces the more accurate approximation Gr. Therefore, if  
accuracy of approximation for a given order is critical, both should be tried.  
As noted previously, if an approximation involving an unstable system is  
desired, the default {onepass=1}is specified.  
Behaviors  
The following explanation deals first with the keyword {onepass}.  
Suppose that σ1, σ2,..., σns are the Hankel Singular values of S, which has  
transfer function matrix G(s). Suppose that the singular values are ordered  
so that:  
σ1 = σ2 = ... = σn > σn + 1... = σn + 1... = σn > σn + 1...  
1
1
1
2
2
> σn  
= σn  
= σn (=σns) 0  
m 1 + 1  
m 1 + 2  
m
Thus, there are n1 equal values, followed by n2 n1 equal values, followed  
by n3 n2 equal values, and so forth.  
The order nsrof Gr(s) cannot be arbitrary when there are equal Hankel  
singular values. In fact, the orders shown in Table 2-1 for the strictly stable  
Gr (all poles in Re[s]<0) and strictly unstable Gu (all poles Re[s]>0) are  
possible (and there are no other possibilities).  
© National Instruments Corporation  
2-15  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
Chapter 2  
Additive Error Reduction  
Table 2-1. Orders of G  
Number of  
Eliminated States  
(Retaining Gu)  
Number of  
Eliminated States  
(Discarding Gu)  
Order of  
Gr nsr  
Order of  
Gu nsu  
0
ns n1  
ns n2  
ns n3  
n1  
n2 n1  
n3 n2  
ns  
ns n1  
ns n2  
n1  
n2  
nm – 1  
0
ns nm – 1  
ns nm – 1  
By abuse of notation, when we say that G is reduced to a certain order, this  
corresponds to the order of Gr(s) alone; the unstable part of Gu(s) of the  
approximation is most frequently thrown away. The number of eliminated  
states (retaining Gu) refers to:  
(# of states in G) – (# of states in Gr) – (# of states in Gu)  
This number is always the multiplicity of a Hankel singular value. Thus,  
when the order of Gr is ni – 1 the number of eliminated states is ni ni – 1 or  
the multiplicity of σn + 1 = σni.  
i – 1  
For each order ni – 1 of Gr(s), it is possible to find Gr and Gu so that:  
G( jω) Gr( jω) Gu( jω) ≤ σn  
i
(Choosing i = 1 causes Gr to be of order zero; identify n0 = 0.) Actually,  
among all “approximations” of G(s) with stable part restricted to having  
degree ni – 1 and with no restriction on the degree of the unstable part, one  
can never obtain a lower bound on the approximation error than σn ; in the  
i
scalar or SISO G(s) case, the Gr(s) which achieves the previous bound is  
unique, while in the matrix or MIMO G(s) case, the Gr(s) which achieves  
the previous bound may not be unique [Glo84]. The algorithm we use to  
find Gr(s) and Gu(s) however allows no user choice, and delivers a single  
pair of transfer function matrices.  
The transfer function matrix Gr( jω) alone can be regarded as a stable  
approximation of G( jω). If the D matrix in Gr( jω) is approximately  
chosen, (and the algorithm ensures that it is), then:  
G( jω) Gr( jω) ≤ σn + σn + ... + σns  
(2-3)  
i
i + 1  
Xmath Model Reduction Module  
2-16  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 2  
Additive Error Reduction  
Thus, the penalty for not being allowed to include Gu in the approximation  
is an increase in the error bound, by σn + 1 + ... +σns. A number of  
i
theoretical developments hinge on bounding the Hankel singular values of  
Gr(s) and Gu(–s) in terms of those of G(s). With Gr(s) of order ni – 1, there  
holds:  
σk(Gr ) ≤ σk(G)k = 1, 2, …, ni 1  
The transfer function matrix Gu(s), being unstable, does not have Hankel  
singular values; however, Gu(–s) (which is stable) does have Hankel  
singular values. They satisfy:  
σk[Gu(s)] ≤ σn + k(G)  
i
In most cases, the Hankel singular values of G(s) are distinct. If,  
accordingly,  
G Gr Gu = σi  
then Gr has degree (i – 1), Gu has degree ns i and  
G Gr = σi + σi + 1 + ... + σns  
(2-4)  
Observe that the bound (Equation 2-3 or Equation 2-4), which is not  
necessarily obtained, is one half that applying for both balanced truncation  
(as implemented by balmoore( )or, effectively, by redschur( )); it  
also is one half that obtained when applying mreduceto a balanced  
realization. In general, the D matrices of G and Gr are different, that is,  
G() Gr() (in contrast to balmoore( )and redschur( )). Similarly,  
G(0) Gr(0) in general (in contrast to the result when mreduceis applied  
to a balanced realization). The price paid for obtaining a smaller error  
bound overall through Hankel norm reduction is that one no longer  
(normally) secures zero error at ω = or ω = 0.  
Two special cases should be noted. If nsr= 0 then Gr(s) is a constant only.  
This constant can be added onto Gu(s), so that G(s) is then being  
approximated by a totally unstable transfer function matrix, with error σ1;  
this type of approximation is known as Nehari approximation. The second  
special case arises when nsr= nm – 1 (or NS – 1 if the smallest Hankel  
singular value has multiplicity 1). In this case, Gu(s) becomes a constant,  
which can then be lumped in with Gr(s), so that G(s), of degree NS, is then  
© National Instruments Corporation  
2-17  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 2  
Additive Error Reduction  
being approximated by a stable Gr(s) with the actual error (as opposed to  
just the error bound) satisfying:  
G(s) Gr(s) = σns  
Note Gr is optimal, that is, there is no other Gr achieving a lower bound.  
Onepass Algorithm  
The first steps of the algorithm are to obtain the Hankel singular values of  
G(s) (by using hankelsv( )) and identify their multiplicities. (Stability of  
G(s) is checked in this process.) If the user has specified nsrand this does  
not coincide with one of 0,n1,n2, ... an error message is obtained; generally,  
all the σi are different, so the occurrence of error messages will be rare.  
The next step of the algorithm is to calculate the sum G(s) = Gr(s) + Gu(s),  
following [SCL90]. (A separate function ophred( )is called for this  
purpose.) The controllability and observability grammians P and Q are  
found in the usual way.  
AP + PA= –BB′  
QA + AQ = –CC  
and then a singular value decomposition is obtained of the  
matrixQP σ2n I:  
i
V1′  
SB  
0
= QP σ2n I  
U1 U2  
i
V2′  
0 0  
There are precisely ni ni – 1 zero singular values, this being the multiplicity  
of σn . Next, the following definitions are made:  
i
U1′  
A11 A12  
A21 A22  
=
2n A+ QAP)(V1V2)  
i
U2′  
U1′  
B1  
=
QB  
U2′  
B2  
[C1 C2] = CP[V1 V2]  
Xmath Model Reduction Module  
2-18  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 2  
Additive Error Reduction  
and finally:  
#
1  
˜
A = SB (A11 A12A22A21)  
#
1  
˜
B = SB (B1 A12A22B2)  
#
#
˜
C = C1 C2A22A21  
#
˜
D = D C2A22B2  
˜
These four matrices are the constituents of the system matrix of G(s),  
where:  
˜
G(s) = Gr(s) + Gu(s)  
Digression:  
This choice is related to the ideas of [Glo84] in the following way;  
˜
in [Glo84], the complete set is identified of G(s) satisfying  
˜
G(jω) G(jω) = σn  
i
˜
with G having a stable part of order ni – 1. The set is parameterized in  
terms of a stable transfer function matrix K(s) which has to satisfy  
C2 + K(s)B2 = 0  
I K′(jω)K( jω) ≤ 0 for all ω  
with C2, B2 being two matrices appearing in the course of the algorithm  
of [Glo84], and enjoying the property C2C2 = B2B2. The particular  
choice  
K(s) = –C2(C2C2)#B2  
in the algorithm of [Glo84] and flagged in corollary 7.3 of [Glo84] is  
equivalent to the previous construction, in the sense of yielding the  
˜
same , though the actual formulas used here and in [Glo84] for the  
Gs  
construction procedure are quite different. In a number of situations,  
including the case of scalar (SISO)G(s), this is the only choice.  
The next step of the algorithm is to call stable( ), to separate G(s) into  
˜
˜
its stable and unstable parts, call them  
and  
, stable( )will  
Gu(s)  
G(s)  
˜
˜
always assign the matrix to  
, and the final step of the algorithm is  
Gr(s)  
D
© National Instruments Corporation  
2-19  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 2  
Additive Error Reduction  
˜
to choose the D matrix of Gr(s), by splitting D between Gr(s) and Gu(s).  
This is done by using a separate function ophiter( ). Suppose Gu(s) is  
the unstable output of stable( ), and let K(s) = Gu(–s). By applying the  
multipass Hankel reduction algorithm, described further below, K(s) is  
reduced to the constant K0 (the approximation), which satisfies,  
K(s) K0 ≤ σ1(K) + ... + σ σn n (K)  
s
i
≤ σn + 1(G) + ... + σn (G)  
i
s
that is, if it is larger than,  
ns  
Gu(s) K0  
σk(G)  
k = ni + 1  
then one chooses:  
˜
Gr = Gr + K0  
˜
Gu = Gu + K0  
This ensures satisfaction of the error bound for G Gr given previously,  
because:  
˜
˜
˜
G Gr  
=
G Gr Gu + (Gu K0)  
˜
˜
= G Gr Gu + K K0  
≤ σn (G) + σn + 1(G) + ... + σn (G)  
i
i
s
Multipass Algorithm  
We now explain the multipass algorithm. For simplicity in first explaining  
the idea, suppose that the Hankel singular values at every stage or pass are  
distinct.  
1. Find a stable order ns – 1 approximation Gn – 1(s) of G(s) with:  
G( jω) Gns 1jω = σns(G)  
(This can be achieved by the algorithm already given, and there is no  
unstable part of the approximation.)  
Xmath Model Reduction Module  
2-20  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
Chapter 2  
Additive Error Reduction  
2. Find a stable order ns – 2 approximation Gns – 2 of Gns – 1(s), with  
Gns 1( jω) Gns 2( jω) = σns 1(Gns 1  
)
.
.
.
3. (Step ns–nr):  
Find a stable order nsr approximation of Gnsr + 1  
,
with Gnsr + 1( jω) Gnsr( jω) = σnsr + 1(Gnsr + 1  
)
Then, because σi(Gns 1) ≤ σi(G) for i < ns,  
σi(Gns 2) ≤ σi(Gns 1  
)
for i s i, ..., this being a property of the algorithm, there follows:  
G( jω) Gnsr( jω) ≤ σnsr + 1(Gnsr + 1) + ... + σns(G)  
ns  
σi(G)  
i = nsr + 1  
The only difference that arises when singular values have multiplicity in  
excess of 1 is that the degree reduction at a given step is greater. Thus, if  
σns(G) has multiplicity k, then G(s) is approximated by Gns k(s) of degree  
ns k.  
No separate optimization of the D matrix of Gnsr is required. The  
approximation error bound is the same as for the first algorithm. The actual  
approximation error may be different. Notice that this second algorithm  
does not calculate an unstable Gu(s) such that  
G(jω) Gnsr(jω) Gu(jω) = σnsr + 1  
Discrete-Time Systems  
No special algorithm is presented for discrete-time systems. Any stable  
discrete-time transfer-function matrix G(z) can be used to define a stable  
continuous-time transfer-function matrix by a bilinear transformation, thus  
α + s  
-----------  
H(s) = G  
α s  
when α is a positive constant.  
© National Instruments Corporation  
2-21  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 2  
Additive Error Reduction  
We use sysZ to denote G(z) and define:  
bilinsys=makepoly([-1,a]/makepoly([1,a])  
as the mapping from the z-domain to the s-domain. The specification is  
reversed because this function uses backward polynomial rotation. Hankel  
norm reduction is then applied to H(s), to generate, a stable reduced order  
approximation Hr(s) and unstable Hu(s) such that:  
H Hr Hu = σi  
H Hr = σi + σn + 1 + ... + σns  
i
Here, the sni are the Hankel singular values of both G(z) and H(s); they are  
the same:  
z 1  
z + 1  
-----------  
Gr(z) = H α  
r  
z 1  
z + 1  
-----------  
α
Gu(z) = H  
u  
We then implement the s-domain equivalent with:  
sysS=subsys(sysZ,bilinsys)  
There is no simple rule for choosing α; the choice α = 1 is probably as good  
as any. The orders of Gr and Gu are the same as those of Hr and Hu,  
respectively. The error formulas are as follows:  
G(ejω) Gr(ejω) Gu(ejω) = σn  
i
G(ejω) Gr(ejω) ≤ σn + σn + 1 + ...σns  
i
i
Impulse Response Error  
If Gr is determined by the first (single-pass) algorithm, the impulse  
response error (for t > 0) between the impulse responses of G and Gr can  
be bounded. As shown in Corollary 9.9 of [Glo84], if Gr is of degree i – 1  
and the multiplicity of the ith larger singular value σi of G is r, then:  
σj[G Gr] ≤ σiG for j = 1, 2, ..., 2i 2 + r  
≤ σj i + 1(G) for j = 2i 1 + r, ...,ns + i 1  
Xmath Model Reduction Module  
2-22  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
Chapter 2  
Additive Error Reduction  
It follows by a result of [BoD87] that the impulse response error for t > 0  
satisfies:  
ns  
g(t) gr(t) 1 2 (2i 2 + ri(G) +  
σj(G)  
i + r  
Evidently, Hankel norm approximation ensures some form of  
approximation of the impulse response too.  
Unstable System Approximation  
A transfer function G(s) with stable and unstable poles can be reduced by  
applying stable( )to separate G(s) into a stable and unstable part. The  
former is reduced and then the unstable part can be added back on. For  
additional information about the ophank( )function, refer to the Xmath  
Help.  
Related Functions  
stable(), redschur(), balmoore()  
© National Instruments Corporation  
2-23  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
3
Multiplicative Error Reduction  
This chapter describes multiplicative error reduction presenting  
two reasons to consider multiplicative rather than additive error reduction,  
one general and one specific.  
Selecting Multiplicative Error Reduction  
The general reason to use multiplicative error reduction is that many  
specifications are given using decibels; 1 db corresponds to a  
multiplicative error of about 12%. Specifications regarding phase shift also  
can be regarded as multiplicative error statements: 0.05 radians of phase  
shift is like 5% multiplicative error also.  
The more specific reason arises in considering the problem of plant  
approximation, with a high order (possibly very high order) plant being  
initially prescribed, with no controller having been designed, and with a  
requirement to provide a simpler model of the plant, possibly to allow  
controller design. Consider the arrangement of Figure 3-1, controller C(s)  
designed for G(s)j with G(s) = (I + )G(s).  
ˆ
C
G
Figure 3-1. Controller C(s) Designed for Multiplicative Error Reduction  
ˆ
The full order plant is G = (I + )G, and the reduced order model is G.  
1  
ˆ ˆ  
Since  
, this means that is the multiplicative error.  
(G G)G = ∆  
Another way one could measure the multiplicative error would be as  
ˆ ˆ  
1. In the matrix plant case, interchange of the order of the  
(G G)G  
product gives two more possibilities again.  
The following multiplicative robustness result can be found in [Vid85].  
© National Instruments Corporation  
3-1  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
         
Chapter 3  
Multiplicative Error Reduction  
Multiplicative Robustness Result  
ˆ
ˆ ˆ 1  
Suppose C stabilizes G, that = (G G)G has no jω-axis poles, and  
ˆ
that G has the same number of poles in Re[s] 0 as . If for all ω,  
G
1  
ˆ
ˆ
∆( jω) G( jω)C( jω)[I + G( jω)C( jω)] < 1  
(3-1)  
then C stabilizes G.  
This result indicates that if a controller C is designed to stabilize a nominal  
ˆ
or reduced order model G , satisfaction of Equation 3-1 ensures that the  
controller also will stabilize the true plant G.  
In reducing a model of the plant, there will be concern not just to have this  
type of stability property, but also concern to have as little error as possible  
ˆ
between the designed system (based on ) and the true system (based  
G
on G). Extrapolation of the stability result then suggests that the goal  
should be not just to have Equation 3-1, but to minimize the quantity on the  
left side of Equation 3-1, or its greatest value:  
1  
ˆ
ˆ
max{ ∆( jω) G( jω)C( jω)[I + G( jω)C( jω)]  
}
ω
However, there are difficulties. The principal one is that if we are reducing  
the plant without knowledge of the controller, we cannot calculate the  
measure because we do not know C(jω). Nevertheless, one could presume  
that, for a well designed system, GC(I + GC)1 will be close to I over the  
operating bandwidth of the system, and have smaller norm than 1 (tending  
to zero as ω→∞ in fact) outside the operating bandwidth of the system.  
This suggests that in the absence of knowledge of C, one should carry out  
multiplicative approximation by seeking to minimize:  
ˆ
ˆ
max ∆(jω)  
=
∆(jω)  
ω
This is the prime rationale for (unweighted) multiplicative reduction of a  
plant.  
Two other points should be noted. First, because frequencies well beyond  
1 will be small, it is in a sense,  
ˆ
ˆ
the closed-loop bandwidth,  
GC(I + GC)  
wasteful to seek to have (jω) small at very high frequencies. The choice  
of maxω ∆(jω) as the index is convenient, because it removes a  
requirement to make assumptions about the controller, but at the same time  
it does not allow ∆(jω) to be made even smaller in the closed-loop  
Xmath Model Reduction Module  
3-2  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 3  
Multiplicative Error Reduction  
bandwidth at the expense of being larger outside this bandwidth, which  
would be preferable.  
ˆ ˆ 1  
(G G)G  
Second, the previously used multiplicative error is  
. In the  
ˆ ˆ 1  
algorithms that follow, the error  
δ = (G G)G  
appears. It is easy to  
check that:  
∆(jω)  
1 ∆(jω)  
-------------------------------  
δ(jω)  
and  
δ(jω)  
------------------------------  
∆(jω)  
1 δ(jω)  
This means that if either bound is small, so is the other, with the bounds  
approximately equal.  
Two algorithms for multiplicative reduction are presented: bst( ),  
a mnemonic for balanced stochastic truncation, and mulhank( ).  
Roughly, they relate to one another in the same way that redschur( )  
and ophank( )relate, that is, one focuses on balanced realization  
truncation and the other on Hankel norm approximation. Some of the  
similarities and differences are as follows:  
When the errors are small, the error bound formula for bst( )is  
about one half of that for bst( ).  
With bst( ), the actual multiplicative error as a function of frequency  
goes to zero as ω→∞ (or, after using an optional transformation given  
in the algorithm description, to zero as ω→ 0); with mulhank( ), the  
error tends to be more uniform as a function of frequency.  
bst( )can handle nonsquare reduction, while mulhank( )cannot.  
Both algorithms are restricted to stable G(s); both preserve right half  
plane zeros, that is, these zeros are copied into the reduced order  
object; both have difficulties with jω-axis zeros of G(s), but these  
difficulties are not insuperable.  
bst( )  
[SysR,HSV] = bst(Sys,{nsr,left,right,bound,method})  
The bst( )function calculates a balanced stochastic truncation of Sysfor  
the multiplicative case.  
© National Instruments Corporation  
3-3  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 3  
Multiplicative Error Reduction  
The objective of the algorithm is to approximate a high-order stable transfer  
function matrix G(s) by a lower-order Gr(s) with either inv(g)(g-gr)or  
(g-gr)inv(g)minimized, under the condition that Gr is stable and of the  
prescribed order.  
Restrictions  
This function has the following restrictions:  
The user must ensure that the input system is stable and nonsingular at  
s = infinity.  
The algorithm may be problematic if the input system has a zero on the  
jω-axis.  
Only continuous systems are accepted; for discrete systems use  
makecontinuous( )before calling bst( ), then discretize the  
result.  
Sys=bst(makecontinuous(SysD));  
SysD=discretize(Sys);  
Algorithm  
The modifications described in this section allow you to circumvent the  
previous restrictions.  
The objective of the algorithm is to approximate a high order stable transfer  
function matrix G(s) by a lower order Gr(s) with, in the square G(s) case,  
(G Gr)G1  
G1(G Gr)  
(approximately) minimized,  
either  
or  
under the constraint that Gr is stable and of prescribed order nsr. In case  
G is not square but has full row rank, the algorithm seeks to minimize:  
(G Gr) (GG*)1(G Gr)  
*
Recall that X*(s) = X′(s) so that when s = jω,  
X*( jω) = X*( jω)  
When G is not square but has full column rank, the algorithm seeks to  
minimize:  
(G Gr)(G*G)1(G Gr)  
*
Xmath Model Reduction Module  
3-4  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
Chapter 3  
Multiplicative Error Reduction  
These cases are secured with the keywords rightand left, respectively.  
If the wrong option is requested for a nonsquare G(s), an error message will  
result.  
The algorithm has the property that right half plane zeros of G(s) remain as  
right half plane zeros of Gr(s). This means that if G(s) has order nsr with n+  
zeros in Re[s] > 0, Gr(s) must have degree at least n+, else, given that it has  
n+ zeros in Re[s] > 0 it would not be proper, [Gre88].  
The conceptual basis of the algorithm can best be grasped by considering  
the case of scalar G(s) of degree n. Then one can form a minimum phase,  
stable W(s) with |W(jω)|2 = |G(jω)|2 and then an all-pass function (the phase  
function) W–1(–s) G(s). This all pass function has a mixture of stable and  
unstable poles, and it encodes the phase of G(jω). Its stable part has n  
Hankel singular values σi with σi 1, and the number of σi equal to 1 is the  
same as the number of zeros of G(s) in Re[s] > 0. State-variable realizations  
of W,G and the stable part of W–1(–s)G(s) can be connected in a nice way,  
and when the stable part of W–1(–s)G(s) has a balanced realization, we say  
that the realization of G is stochastically balanced. Truncating the balanced  
realization of the stable part of W–1(–s)G(s) induces a corresponding  
truncation in the realization of G(s), and the truncated realization defines an  
approximation of G. Further, a good approximation of a transfer function  
encoding the phase of G somehow ensures a good approximation, albeit in  
a multiplicative sense, of G itself.  
Algorithm with the Keywords right and left  
The following description of the algorithm with the keyword rightis  
based on ideas of [GrA86] developed in [SaC88]. The procedure is almost  
the same when leftis specified, except the transpose of G(s) is used; the  
algorithm finds an approximation in the same manner as for right, but  
transposes the approximation to yield the desired Gr(s).  
1. The algorithm checks  
That the system is state-space, continuous, and stable  
That a correct option has been specified if the plant is nonsquare  
That D is nonsingular; if the plant is nonsquare, DD´ must be  
nonsingular  
© National Instruments Corporation  
3-5  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 3  
Multiplicative Error Reduction  
2. With G(s) = D + C(sI A)–1B and stable, with DD´ nonsingular and  
G(jω) G'(–jω) nonsingular for all ω, part of a state variable realization  
of a minimum phase stable W(s) is determined such that  
W´(–s)W(s) = G(s)G´(–s) with  
W(s) = DW + CW(sI Aw)1BW  
The state variable matrices in W(s) are obtained as follows. The  
controllability grammian P associated with G(s) is first found from  
AP + PA´ + BB´=0 then AW = A, BW = PC´ + BD´.  
When G(s) is square, the algorithm checks to see if there is a zero or  
singularity of G(s) close to the jω-axis (the zeros are given by the  
eigenvalues of A BD–1C and are computed reliably with the aid of  
schur( )). If one is found, you are warned that results may be  
unreliable. Next, a stabilizing solution Q is found for the following  
Riccati equation:  
QA + AQ + (C BWQ)′(DD′)1(C BWQ) = 0  
The singriccati( )function is used; failure of the nonsingularity  
condition on G(jω)G´(–jω) will normally result in an error message  
that no stabilizing solution exists. To obtain the best numerical results,  
singriccati( )is invoked with the keyword {method="schur"}.  
Although DW, CW are not needed for the remainder of the algorithm,  
they are simply determined in the square case by  
DW = DCW = D1(C BWQ)  
with minor modification in the nonsquare case. The real point of the  
algorithm is to compute P and Q; the matrix Q satisfies (square or  
nonsquare case).  
QA + AQ + CWCW = 0  
P, Q are the controllability and observability grammians of the transfer  
function CW(sI A)–1B. This transfer function matrix, it turns out, is  
the strictly proper, stable part of θ(s) = WT(–s)G(s), which obeys the  
matrix all-pass property θ(s)θ´(–s) = I, and is the phase matrix  
associated with G(s).  
3. Compute ordered Schur decompositions of PQ, with the eigenvalues  
of PQ is ascending and descending order. Obtain the phase matrix  
Hankel singular values, that is, the Hankel singular values of the  
Xmath Model Reduction Module  
3-6  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 3  
Multiplicative Error Reduction  
strictly proper stable part of θ(s), as the square roots of the eigenvalues  
of PQ. Call these quantities νi. The Schur decompositions are,  
VAPQVA = Sasc  
VDPQVD = Sdes  
where VA, VD are orthogonal and Sasc, Sdes are upper triangular.  
4. Define submatrices as follows, assuming the dimension of the reduced  
order system nsris known:  
0
Insr  
Vlbig = VA  
Vrbig = VD  
Insr  
0
Determine a singular value decomposition,  
UebigSebigVebig = VlbigVrbig  
and then define transformation matrices:  
1 2  
Slbig = VlbigUebig  
S
S
ebig  
1 2  
ebig  
Srbig = VrbigVebig  
The reduced order system Gr is:  
AR = SlbigASrbig  
BR = Slbig  
B
AR = CSrbig  
DR = D  
where step 4 is identical with that used in redschur( ), except  
the matrices P, Q which determine VA, VD and so forth, are the  
controllability and observability grammians of CW(sI A)–1B rather  
than of C(sI A)–1B, the controllability grammian of G(s) and the  
observability grammian of W(s).  
The error formula [WaS90] is:  
vi  
G1(G Gr) 2  
------------  
(3-2)  
1 vi  
All νi obey νi 1. One can only eliminate νi where νi < 1. Hence, if nsris  
chosen so that νnsr + 1 = 1, the algorithm produces an error message. The  
algorithm also checks that nsrdoes not exceed the dimension of a minimal  
© National Instruments Corporation  
3-7  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 3  
Multiplicative Error Reduction  
state-variable representation of G. In this case, the user is effectively asking  
for Gr = G. When the phase matrix has repeated Hankel singular values,  
they must all be included or all excluded from the model, that is,  
νnsr = νnsr + 1 is not permitted; the algorithm checks for this.  
The number of νi equal to 1 is the number of zeros in Re[s]>0 of G(s), and  
as mentioned already, these zeros remain as zeros of Gr(s).  
If erroris specified, then the error bound formula (Equation 3-2) in  
conjunction with the νi values from step 3 is used to define nsrfor step 4.  
For nonsquare G with more columns than rows, the error formula is:  
ns  
(G Gr) (G*G)1(G Gr) 1 2 2  
------------  
vi  
*
1 vi  
i = nsr + 1  
If the user is presented with the νi, the error formula provides a basis for  
intelligently choosing nsr. However, the error bound is not guaranteed to  
be tight, except when nsr = ns – 1.  
Securing Zero Error at DC  
The error G–1(G Gr) as a function of frequency is always zero at ω = .  
When the algorithm is being used to approximate a high order plant by a  
low order plant, it may be preferable to secure zero error at ω = 0. A method  
for doing this is discussed in [GrA90]; for our purposes:  
1. We need a bilinear transformation of sys = 1/z. Given G(s) we generate  
H(s) through:  
bilinsys=makepoly([b3,b4]/makepoly([b1,b2])  
sys=subsys(sys,bilinsys)  
2. Reduce with the previous algorithm:  
[sr,nsr,hsv] = bst(sys)  
3. Use the bilinear transformation s = 1/z again:  
[sr1,nsr1] = bilinear(sr,nsr,[0,1,1,0])  
The νi are the same for G(s) and H(s) = G(s–1). The error bound formula is  
the same; H is stable and H(jω)H'(–jω) of full rank for all ω including  
ω = if and only if G has the same property; right half plane zeros of G are  
still preserved by the algorithm. The error G–1(G Gr), though now zero at  
ω = 0, is in general nonzero at ω = .  
Xmath Model Reduction Module  
3-8  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 3  
Multiplicative Error Reduction  
Hankel Singular Values of Phase Matrix of Gr  
The νi, i = 1,2,...,ns have been termed above the Hankel singular values of  
the phase matrix associated with G. The corresponding quantities for Gr are  
νi, i = 1,..., nsr.  
Further Error Bounds  
The introduction to this chapter emphasized the importance of the error  
measures  
(G Gr)Gr 1 or Gr1(G Gr)  
(G Gr)G1  
for plant reduction, as opposed to  
or G1(G Gr)  
The BST algorithm ensures that in addition to (Equation 3-2), there holds  
[WaS90a].  
ns  
vi  
Gr1(G Gr) 2  
------------  
1 vi  
i = nsr + 1  
which also means that for a scalar system,  
ns  
Gr  
-----  
vi  
20log  
8.69 2  
------------ dB  
10 G  
1 v  
i  
i = nsr + 1  
and, if the bound is small:  
ns  
vi  
------------  
phase(G) phase(Gr) ≤  
radians  
1 vi  
i = nsr + 1  
Reduction of Minimum Phase, Unstable G  
For square minimum phase but not necessarily stable G, it also is possible  
to use this algorithm (with minor modification) to try to minimize (for Gr  
of a certain order) the error bound  
(G Gr)Gr 1  
© National Instruments Corporation  
3-9  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
         
Chapter 3  
Multiplicative Error Reduction  
which also can be relevant in finding a reduced order model of a plant.  
The procedure requires G again to be nonsingular at ω = , and to have no  
jω-axis poles. It is as follows:  
1. Form H = G–1. If G is described by state-variable matrices A, B, C, D,  
then H is described by A BD–1C, BD–1, –D–1C, D–1. H is square,  
stable, and of full rank on the jω-axis.  
2. Form Hr of the desired order to minimize approximately:  
H1(H Hr)  
3. Set Gr = H–1  
.
r
Observe that  
H1(H Hr) = I H1Hr  
= I GGr1  
= (Gr G)Gr1  
The reduced order Gr will have the same poles in Re[s] > 0 as G, and  
be minimum phase.  
Imaginary Axis Zeros (Including Zeros at )  
We shall now explain how to handle the reduction of G(s) which has a rank  
drop at s = or on the jω-axis. The key is to use a bilinear transformation,  
[Saf87]. Consider the bilinear map defined by  
z a  
s = ------------------  
bz + 1  
s + a  
z = --------------  
bs + 1  
˜
where 0 < a < b–1 and mapping G(s) into G(s) through:  
s a  
bs + 1  
˜
------------------  
G(s) = G  
s + a  
bs + 1  
˜
--------------  
G(s) = G  
Xmath Model Reduction Module  
3-10  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 3  
Multiplicative Error Reduction  
The values of G(s), as shown in Figure 3-2, along the jω-axis are  
˜
the same as the values of  
around a circle with diameter defined by  
G(s)  
[a j0, b–1 + j0] on the positive real axis.  
a
b1  
˜
G(s)  
G(s)  
values  
values  
˜
Figure 3-2. Bilinear Mapping from G(s) to (Gs) (Case 1)  
˜
Also, the values of G(s), as shown in Figure 3-3, along the jω-axis are  
the same as the values of G(s) around a circle with diameter defined by  
[–b–1 + j0, –a + j0].  
b–1  
-a  
˜
G(s)  
G(s)  
values  
values  
˜
Figure 3-3. Bilinear Mapping from G(s) to (Gs) (Case 2)  
We can implement an arbitrary bilinear transform using the subsys( )  
function, which substitutes a given transfer function for the s- or z-domain  
operator.  
s a  
bs + 1  
˜
------------------  
To implement G(s) = G  
use:  
gtildesys=subsys(gsys,makep([-b,1]/makep([1,-a])  
s + a  
s + 1  
˜
-----------  
To implement G(s) = G  
use:  
gsys=subsys(gtildesys,makep([b,1]/makep([1,a])  
Note The systems substituted in the previous calls to subsys invert the function  
specification because these functions use backward polynomial rotation.  
© National Instruments Corporation  
3-11  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 3  
Multiplicative Error Reduction  
Any zero (or rank reduction) on the jω-axis of G(s) becomes a zero (or rank  
˜
reduction) in Re[s] > 0 of  
, and if G(s) has a zero (or rank reduction)  
G(s)  
˜
at infinity, this is shifted to a zero (or rank reduction) of  
at the point  
G(s)  
b–1, (in Re[s] > 0). If all poles of G(s) are inside the circle of diameter  
˜
[–b–1 + j0, a + j0], all poles of  
will be in Re[s] < 0, and if G(s) has no  
G(s)  
˜
zero (or rank reduction) on this circle,  
reduction) on the jω-axis, including ω = .  
will have no zero (or rank  
G(s)  
If G(s) is nonsingular for almost all values of s, it will be nonsingular or  
have no zero or rank reduction on the circle of diameter [–b–1 + j0, –a + j0]  
for almost all choices of a,b. If a and b are chosen small enough, G(s) will  
have all its poles inside this circle and no zero or rank reduction on it, while  
˜
G(s) then will have all poles in Re[s] < 0 and no zero or rank reduction on  
the jω-axis, including s = . The steps of the algorithm, when G(s) has a  
zero on the jω-axis or at s = , are as follows:  
s a  
bs + 1  
–1  
˜
------------------  
1. For small a,b with 0 < a < b , form G(s) = G  
as shown for  
gtildesys.  
˜
˜
˜
2. Reduce G(s) to Gr(s), this being possible because G(s) is stable and  
has full rank on s = jω, including ω = ∞.  
s + a  
bs + 1  
˜
--------------  
3. Form Gr(s) = G  
as shown for gsys.  
r  
The error G1(G Gr) will be overbounded by the error  
G (G Gr) , and Gr will contain the same zeros in Re[s] 0 as G.  
1  
˜
˜
˜
If there is no zero (or rank reduction) of G(s) at the origin, one can take  
a = 0 and b–1 = bandwidth over which a good approximation of G(s) is  
needed, and at the very least b–1 sufficiently large that the poles of G(s)  
lie in the circle of diameter [–b–1 + j0, –a + j0]. If there is a zero or rank  
reduction at the origin, one can replace a = 0 by a = b. It is possible to take  
b too small, or, if there is a zero at the origin, to take a too small. The user  
will be presented with an error message that there is a jω-axis zero and/or  
that the Riccati equation solution may be in error. The basic explanation is  
˜
that as b 0, and thus a 0, the zeros of  
approach those of G(s).  
G
G(s)  
˜
Thus, for sufficiently small b, one or more zeros of (s) may be identified  
as lying on the imaginary axis. The remedy is to increase a and/or b above  
the desirable values.  
The procedure for handling jω-axis zeros or zeros at infinity will be  
deficient if the number of such zeros is the same as the order of G(s)—for  
example, if G(s) = 1/d(s), for some stable d(s). In this case, it is possible  
Xmath Model Reduction Module  
3-12  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 3  
Multiplicative Error Reduction  
again with a bilinear transformation to secure multiplicative  
approximations over a limited frequency band. Suppose that  
s
˜
--------------  
G(s) = G  
εs + 1  
˜
Create a system that corresponds to  
with:  
G(s)  
gtildesys=subs(gsys,(makep([-eps,1])/makep([1,-]))  
bilinsys=makep([eps,1])/makep([1,0])  
sys=subsys(sys,bilinsys)  
Under this transformation:  
˜
Values of G(s) along the jω-axis correspond to values of G(s) around  
a circle in the left half plane on diameter (–ε–1 + j0, 0).  
˜
Values of G(s) along the jω-axis correspond to values of G(s) around  
a circle in the right half plane on diameter (0, ε–1 + j0).  
˜
Multiplicative approximation of G(s) (along the jω-axis) corresponds to  
multiplicative approximation of G(s) around a circle in the right half plane,  
touching the jω-axis at the origin. For those points on the jω-axis near the  
circle, there will be good multiplicative approximation of G( jω). If it is  
desired to have a good approximation of G(s) over an interval [–j, j],  
then a choice such as ε–1 = 5 or 10 needs to be made. Reduction then  
proceeds as follows:  
˜
1. Form G(s).  
˜
2. Reduce G(s) through bst( ).  
˜
3. Form Gr(s) = –Gr(s ⁄ (1 εs)) with:  
gsys=subsys(gtildesys(gtildesys,  
makep([-eps,-1])/makep[-1,-0]))  
Notice that the number of zeros of G(s) in the circle of diameter  
(0, ε1 + j0)s  
sets a lower bound on the degree of Gr(s)—for such zeros become right half  
˜
plane zeros of G(s), and must be preserved by bst( ). Obviously, zeros at  
s = are never in this circle, so a procedure for reducing G(s) = 1/d(s) is  
available.  
© National Instruments Corporation  
3-13  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 3  
Multiplicative Error Reduction  
There is one potential source of failure of the algorithm. Because G(s) is  
˜
stable, G(s) certainly will be, as its poles will be in the left half plane circle  
˜
on diameter (ε = j0, 0). If Gr(s) acquires a pole outside this circle  
(but still in the left half plane of course)—and this appears possible in  
principle—Gr(s) will then acquire a pole in Re [s] > 0. Should this difficulty  
be encountered, a smaller value of ε should be used.  
Related Functions  
redschur(), mulhank()  
mulhank( )  
[SysR,HSV] = mulhank(Sys,{nsr,left,right,bound,method})  
The mulhank( )function calculates an optimal Hankel norm reduction of  
Sysfor the multiplicative case.  
Restrictions  
This function has the following restrictions:  
The user must ensure that the input system is stable and nonsingular at  
s = infinity.  
The algorithm may be problematic if the input system has a zero on the  
jω-axis.  
Only continuous systems are accepted; for discrete systems use  
makecontinuous( )before calling mulhank( ), then discretize  
the result.  
Sys=mulhank(makecontinuous(SysD));  
SysD=discretize(Sys);  
Algorithm  
The objective of the algorithm, like bst( ), is to approximate a high order  
square stable transfer function matrix G(s) by a lower order Gr(s) with  
either (G Gr)G1 or G1(G Gr) (approximately) minimized,  
under the constraint that Gr is stable and of prescribed order.  
The algorithm has the property that right half plane zeros of G(s) are  
retained as zeros of Gr(s). This means that if G(s) has order NS with N+  
zeros in Re[s] > 0, Gr(s) must have degree at least N+—else, given that it  
has N+ zeros in Re[s] > 0 it would not be proper, [GrA89].  
Xmath Model Reduction Module  
3-14  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
         
Chapter 3  
Multiplicative Error Reduction  
The conceptual basis of the algorithm can best be grasped by considering  
the case of scalar G(s) of degree n. Then one can form a minimum phase,  
stable W(s) with |W(jω)|2 = |G(jω)|2 and then an all-pass function (the phase  
function) W–1(–s) G(s). This all-pass function has a mixture of stable and  
unstable poles, and it encodes the phase of G(jω). Its stable part has  
n Hankel singular values σi with σi 1, and the number of σi equal to 1  
is the same as the number of zeros of G(s) in Re[s]>0. State-variable  
realizations of W,G and the stable part of W–1(–s)G(s) can be connected in  
a nice way. The algorithm computes an additive Hankel norm reduction of  
the stable part of W–1(–s)G(s) to cause a degree reduction equal to the  
multiplicity of the smallest σi. The matrices defining the reduced order  
object are then combined in a new way to define a multiplicative  
approximation to G(s); as it turns out, there is a close connection between  
additive reduction of the stable part of W–1(–s)G(s) and multiplicative  
reduction of G(s). The reduction procedure then can be repeated on the new  
phase function of the just found approximation to obtain a further reduction  
again in G(s).  
right and left  
A description of the algorithm for the keyword rightfollows. It is based  
on ideas of [Glo86] in part developed in [GrA86] and further developed  
in [SaC88]. The procedure is almost the same when {left}is specified,  
except the transpose of G(s) is used; the following algorithm finds an  
approximation, then transposes it to yield the desired Gr(s).  
1. The algorithm checks that G(s) is square, stable, and that the transfer  
function is nonsingular at infinity.  
2. With G(s) = D + C(sIA)–1B square and stable, with D nonsingular  
[rank(d)must equal number of rows in d] and G(jω) nonsingular for  
all finite ω, this step determines a state variable realization of a  
minimum phase stable W(s) such that,  
W´(–s)W(s) = G(s)G´(–s)  
with:  
W(s) = Dw + Cw(sI–Aw)–1Bw  
The various state variable matrices in W(s) are obtained as follows. The  
controllability grammian P associated with G(s) is first found from  
AP + PA´ + BB´ = 0, then:  
Aw = ABw = PC´+BD´Dw = D´  
The algorithm checks to see if there is a zero or singularity of G(s)  
close to the jω-axis. The zeros are determined by calculating the  
© National Instruments Corporation  
3-15  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 3  
Multiplicative Error Reduction  
eigenvalues of A – B/D * C with the aid of schur( ). If any real part  
of the eigenvaluesis less than eps, a warning is displayed.  
Next, a stabilizing solution Q is found for the following Riccati  
equation:  
QA + AQ + (C BwQ)′(DD′)1(C BwQ) = 0  
The function singriccati( )is used; failure of the nonsingularity  
condition of G(jω) will normally result in an error message. To obtain  
the best numerical results, singriccati( )is invoked with the  
keyword method="schur".  
The matrix Cw is given by C = D1(C BwQ).  
w
Notice that Q satisfies QA + AQ + CwCw = 0, so that P and Q are  
the controllability and observability grammians of  
F(s) = Cw(sI A)1B  
This strictly proper, stable transfer function matrix is the strictly  
proper, stable part (under additive decomposition) of  
θ(s)=W–T(–s)G(s), which obeys the matrix all pass property  
θ(s)θ'(–s)=I. It is the phase matrix associated with G(s).  
3. The Hankel singular values νi of F(s) = C (sI A)1B are  
w
computed, by calling hankelsv( ). The value of nsris obtained if  
not prespecified, either by prompting the user or by the error bound  
formula ([GrA89], [Gre88], [Glo86]).  
ns  
vnsr + 1 G1(G Gr)  
(1 + vj) 1  
(3-3)  
j = nsr + 1  
(with νi ≥ νi + 1  
⋅⋅⋅ being assumed). If νk = νk + 1 = ... = νk + r for some  
k, (that is, νk has multiplicity greater than unity), then νk appears once  
only in the previous error bound formula. In other words, the number  
of terms in the product is equal to the number of distinct νi less than  
νnsr. There are restrictions on nsr. nsrcannot exceed the dimension  
of a minimal realization of G(s); although νi i + 1⋅⋅⋅, nsr must obey  
nnsr > nnsr+1; and while 1 ≥ νi for all i, it is necessary that 1>νnsr + 1. (The  
number of νi equal to 1 is the number of right half plane zeros of G(s).  
They must be retained in Gr(s), so the order of Gr(s), nsr, must at least  
be equal to the number of νi equal to 1.) The software checks all these  
conditions. The minimum order permitted is the number of Hankel  
Xmath Model Reduction Module  
3-16  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
Chapter 3  
Multiplicative Error Reduction  
singular values of F(s) larger than 1– ε (refer to steps 1 through 3 of the  
Restrictions section). The maximum order permitted is the number of  
nonzero eigenvalues of WcWo larger than ε.  
4. Let r be the multiplicity of νns. The algorithm approximates  
F(s) = Cw(sI A)1B  
ˆ
by a transfer function matrix F(s) of order ns – r, using Hankel norm  
approximation. The procedure is slightly different from that used in  
ophank( ).  
Construct an SVD of QP v2nsI:  
V1′  
QP v2NSI = U  
V= [U1U2]  
Σ1 0  
Σ1 0  
V2′  
0 0  
0 0  
with Σ1 of dimension (ns r) × (ns r) and nonsingular. Also, obtain  
an orthogonal matrix T, satisfying:  
B2 + Cw2T = 0  
where B2 and Cw2 are the last r rows of B and Cw , the state variable  
1  
matrices appearing in a balanced realization of Cws(I A) B. It is  
possible to calculate T without evaluating B B, Cw as it turns out (refer  
to [AnJ]), and the algorithm does this. Now with  
ˆ
ˆ
ˆ
ˆ
1 ˆ  
F(s) = DF + CF(sI AF) BF  
ˆ
ˆ
ˆ
ˆ
Fp(s) = CF(sI AF)BF  
there holds:  
ˆ
2
AF = Σ11U1[vnsA+ QAP vnsCwTB′]V1  
BF = Σ11U1[QB + vnsCwT]  
ˆ
ˆ
CF = (CwP + vnsTB′)V′  
ˆ
DF = –vnsT  
© National Instruments Corporation  
3-17  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 3  
Multiplicative Error Reduction  
ˆ
ˆ
Note The expression Fp(s) is the strictly proper part of F(s). The matrix  
vns[F(s) F(s)]  
1  
ˆ
is all pass; this property is not always secured in the multivariable case  
when ophank( )is used to find a Hankel norm approximation of F(s).  
ˆ
ˆ
5. The algorithm constructs G and W , which satisfy,  
ˆ
ˆ
G(s) = G(s) W′(s)[F(s) F(s)]  
and,  
1  
ˆ
W(s) = (I vnsT′)(I vnsT)  
ˆ
{W(s) [F(s) F(s)]G+ (s)}  
through the state variable formulas  
1 ˆ  
ˆ
ˆ
ˆ
(G(s) = (D(I vnsT))[DCF + BWUΣ1](sI AF) BF)  
and:  
1  
ˆ
W(s) = (I vnsT′)D+ (I vnsT′)(I vnsT)  
1  
ˆ
ˆ
ˆ
CF(sI AF) [BFD+ V1C′]  
ˆ
ˆ
ˆ
Continue the reduction procedure, starting with G, W , F and  
repeating the process till Gr of the desired degree nsris obtained.  
^
ˆ
For example, in the second iteration,  
is given by:  
(s)  
G
^
^
ˆ
ˆ
ˆ
ˆ
ˆ
G(s)= G(s) W+ (s)[Fp(s) F(s)]  
(3-4)  
Consequences of Step 5 and Justification of Step 6  
A number of properties are true:  
ˆ
G(s) is of order ns r, with:  
1  
ˆ
G (G G) = vns  
(3-5)  
Xmath Model Reduction Module  
3-18  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 3  
Multiplicative Error Reduction  
ˆ
ˆ
W(s) and Gs stand in the same relation as W(s) and G(s), that is:  
ˆ
ˆ
ˆ
ˆ
W′(s)W(s) = G(s)G′(s)  
ˆ ˆ ˆ ˆ  
ˆ ˆ  
With PAF + AFP = –BFBF , there holds  
ˆ
BWˆ = PCGˆ + BGˆ DGˆ  
or  
ˆ
ˆ
ˆ
ˆ
BFD+ V1C= P(DCF + BWU1Σ1)′ + BF(I vnsT′)D′  
ˆ ˆ  
ˆ ˆ  
ˆ
ˆ
With  
there holds  
ˆ
QAF + AFQ = –CFCF  
1  
CWˆ = DGˆ (CGˆ BWˆ Q)  
or  
(I v T′)(I v T) CF = [D(I vnsT)]1  
1  
ˆ
ns  
ns  
ˆ
ˆ
ˆ
{DCF + BWU11 [BFD+ V1C′]′Q)}  
DWˆ = DGˆ  
1  
ˆ
ˆ
ˆ
is the stable strictly proper part of (W (s))G(s).  
F
ˆ
ˆ
The Hankel singular values of Fp(and F ) are the first as r Hankel  
singular values of F,  
P = Σ1 U1QV1 = V1QU1Σ11  
1  
ˆ
Q = V1PU1Σ1 = Σ1U1PV1  
ˆ
ˆ
Gs has the same zeros in Re[s] > 0 as G(s).  
These properties mean that one is immediately positioned to repeat the  
ˆ
reduction procedure on , with almost all needed quantities being on  
Gs  
hand.  
© National Instruments Corporation  
3-19  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 3  
Multiplicative Error Reduction  
Error Bounds  
The error bound formula (Equation 3-3) is a simple consequence of  
iterating (Equation 3-5). To illustrate, suppose there are three reductions  
ˆ
ˆ
ˆ
G2 G3  
, each by degree one. Then,  
G
G
1  
1  
ˆ
ˆ
G (G G3) = G (G G)  
1  
1  
ˆ ˆ  
ˆ
ˆ
+ G GG (G G2)  
ˆ ˆ 1 ˆ ˆ 1  
ˆ
ˆ
1  
+ G GG G2G2 (G2 G3)  
Also,  
ˆ 1  
G (G G) + I  
ˆ
1  
ˆ
G G  
=
1 + vns  
Similarly,  
Then:  
ˆ 1  
ˆ
ˆ 1  
ˆ
G G2 1 + vns 1, G2 G3 1 + vns 2  
1  
ˆ
G (G G3) ≤ vns + (1 + vns)vns 1 + (1 + vns 1)vns 2  
= (1 + vns)(1 + vns 1)(1 + vns 2) 1  
The error bound (Equation 3-3) is only exact when there is a single  
reduction step. Normally, this algorithm has a lower error bound than  
bst( ); in particular, if the νi are all distinct and vnsr + 1 « 1, the error  
bounds are approximately  
ns  
ns  
vi  
vi  
for mulhank( )  
for bst(  
2
i = nsr + 1  
i = nsr + 1  
Xmath Model Reduction Module  
3-20  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 3  
Multiplicative Error Reduction  
For mulhank( ), this translates for a scalar system into  
ns  
ˆ
86.9  
vi dB < 20log10 Gnsr G  
ns  
i = nsr + 1  
< 8.69  
vi dB  
i = nsr + 1  
and  
ns  
phase error <  
vi radians  
i = nsr + 1  
The bounds are double for bst( ).  
The error as a function of frequency is always zero at ω = for bst( )  
(or at ω = 0 if a transformation s s–1 is used), whereas no such particular  
property of the error holds for mulhank( ).  
Imaginary Axis Zeros (Including Zeros at )  
When G(jω) is singular (or zero) on the jω axis or at , reduction can be  
handled in the same manner as explained for bst( ).  
The key is to use a bilinear transformation [Saf87]. Consider the bilinear  
map defined by  
z a  
s = ------------------  
bz + 1  
s + a  
z = --------------  
bs + 1  
˜
where 0 < a < b–1 and mapping G(s) into G(s) through  
s a  
bs + 1  
˜
------------------  
G(s) = G  
s + a  
bs + 1  
˜
--------------  
G(s) = G  
© National Instruments Corporation  
3-21  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 3  
Multiplicative Error Reduction  
˜
The values of G(s) along the jω-axis are the same as the values of G(s)  
around a circle with diameter defined by [a j0, b–1 + j0] on the positive  
˜
real axis (refer to Figure 3-2). Also, the values of  
along the jω-axis  
G(s)  
are the same as the values of G(s) around a circle with diameter defined by  
[–b–1 + j0, –a + j0].  
We can implement an arbitrary bilinear transform using the subsys( )  
function, which substitutes a given transfer function for the s- or z-domain  
operator, as previously shown.  
s a  
bs + 1  
˜
------------------  
To implement G(s) = G  
use:  
gtildesys=subsys(gsys,makep([-b,1]/makep([1,-a])  
s + a  
bs + 1  
˜
--------------  
To implement G(s) = G  
use:  
gsys=subsys(gtildesys,makep([b,1]/makep([1,a])  
Note The systems substituted in the previous calls to subsys invert the function  
specification because these functions use backward polynomial rotation.  
Any zero (or rank reduction) on the jω-axis of G(s) becomes a zero (or rank  
˜
reduction) in Re[s] > 0 of G(s), and if G(s) has a zero (or rank reduction)  
˜
at infinity, this is shifted to a zero (or rank reduction) of G(s) at the point  
b–1, again in Re[s] > 0. If all poles of G(s) are inside the circle of diameter  
˜
[–b–1 + j0, a + j0], all poles of  
will be in Re[s] < 0, and if G(s) has no  
G(s)  
˜
zero (or rank reduction) on this circle,  
reduction) on the jω-axis, including ω = .  
will have no zero (or rank  
G(s)  
If G(s) is nonsingular for almost all values of s, it will be nonsingular or  
have no zero or rank reduction on the circle of diameter [–b–1 + j0, –a + j0]  
for almost all choices of a,b. If a and b are chosen small enough, G(s) will  
have all its poles inside this circle and no zero or rank reduction on it, while  
˜
then will have all poles in Re[s] < 0 and no zero or rank reduction on  
G(s)  
the jω-axis, including s = .  
The steps of the algorithm, when G(s) has a zero on the jω-axis or at s = ,  
are as follows:  
s a  
bs +1  
˜
1. For small a,b with 0 < a < b–1, form G(s) = G  
as shown for  
------------------  
gtildesys.  
˜
˜
˜
2. Reduce G(s) to Gr(s), this being possible because G(s) is stable and  
has full rank on s = jω, including ω = .  
s + a  
˜
--------------  
3. Form Gr(s) = G  
as shown for gsys.  
r  
bs + 1  
Xmath Model Reduction Module  
3-22  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 3  
Multiplicative Error Reduction  
G1(G Gr)  
The error  
will be overbounded by the error  
1  
˜
˜
˜
, and Gr will contain the same zeros in Re[s] 0 as G.  
G (G Gr)  
If there is no zero (or rank reduction) of G(s) at the origin, one can take  
a = 0 and b–1 = bandwidth over which a good approximation of G(s) is  
needed, and at the very least b–1 sufficiently large that the poles of G(s)  
lie in the circle of diameter [–b–1 + j0, –a + j0]. If there is a zero or rank  
reduction at the origin, one can replace a = 0 by a = b. It is possible to take  
b too small, or, if there is a zero at the origin, to take a too small. In these  
cases an error message results, saying that there is a jω-axis zero and/or that  
the Riccati equation solution may be in error. The basic explanation is that  
˜
as b 0, and thus a 0, the zeros of  
approach those of G(s). Thus,  
G(s)  
˜
for sufficiently small b, one or more zeros of  
lying on the imaginary axis. The remedy is to increase a and/or b above the  
may be identified as  
G(s)  
desirable values.  
The previous procedure for handling jω-axis zeros or zeros at infinity will  
be deficient if the number of such zeros is the same as the order of G(s); for  
example, if G(s) = 1/d(s), for some stable d(s). In this case, it is possible  
again with a bilinear transformation to secure multiplicative  
approximations over a limited frequency band. Suppose that  
s
˜
--------------  
G(s) = G  
εs + 1  
˜
Create a system that corresponds to G(s) with:  
gtildesys=subs(gsys,(makep([-eps,1])/makep([1,-]))  
bilinsys=makep([eps,1])/makep([1,0])  
sys=subsys(sys,bilinsys)  
Under this transformation:  
˜
Values of G(s) along the jω-axis correspond to values of G(s) around  
a circle in the left half plane on diameter (–ε–1 + j0, 0).  
˜
Values of G(s) along the jω-axis correspond to values of G(s) around  
a circle in the right half plane on diameter (0, ε–1 + j0).  
© National Instruments Corporation  
3-23  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 3  
Multiplicative Error Reduction  
˜
Multiplicative approximation of G(s) (along the jω-axis) corresponds to  
multiplicative approximation of G(s) around a circle in the right half plane,  
touching the jω-axis at the origin. For those points on the jω-axis near the  
circle, there will be good multiplicative approximation of G(jω). If a good  
approximation of G(s) over an interval [–j, j] it is desired, then  
ε
–1 = 5 or 10 are good choices. Reduction then proceeds as follows:  
˜
1. FormG(s).  
˜
2. Reduce G(s) through bst( ).  
˜
3. Form Gr(s) = –Gr(s ⁄ (1 εs)) with:  
gsys=subsys(gtildesys(gtildesys,  
makep([-eps,-1])/makep[-1,-0]))  
Notice that the number of zeros of G(s) in the circle of diameter (0, ε–1 + j0)  
sets a lower bound on the degree of Gr(s)—for such zeros become right half  
˜
plane zeros of G(s), and must be preserved by bst( ). Zeros at s = are  
never in this circle, so a procedure for reducing G(s) = 1/d(s) is available.  
There is one potential source of failure of the algorithm. Because G(s) is  
˜
stable, G(s) certainly will be, as its poles will be in the left half plane circle  
on diameter (ε1 = j0, 0). If  
acquires a pole outside this circle  
˜
Gr(s)  
(but still in the left half plane of course)—and this appears possible in  
principle—Gr(s) will then acquire a pole in Re [s] >0. Should this difficulty  
be encountered, a smaller value of ε should be used.  
Related Functions  
singriccati(), ophank(), bst(), hankelsv()  
Xmath Model Reduction Module  
3-24  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
4
Frequency-Weighted Error  
Reduction  
This chapter describes frequency-weighted error reduction problems. This  
includes a discussion of controller reduction and fractional representations.  
Introduction  
Frequency-weighted error reduction means that the error is measured not,  
as previously, by  
E0  
=
G(jω) Gr(jω)  
but rather by  
E1  
=
G(jω) Gr(jω)V(jω)  
(4-1)  
(4-2)  
(4-3)  
or  
or  
E2  
=
W(jω)[G(jω) Gr(jω)]  
E3  
=
W(jω)[G(jω) Gr(jω)]V(jω)  
where W,V are certain weighting matrices. Their presence reflects a desire  
that the approximation process be more accurate at certain frequencies  
(where V or W have large singular values) than at others (where they  
have small singular values). For scalar G(jω), all the indices above are  
effectively the same, with the effective weight just |V(jω)|, |W(jω)|,  
or |W(jω)V(jω)|.  
When the system G is processing signals which do not have a flat spectrum,  
and is to be approximated, there is considerable logic in using a weight. If  
the signal spectrum is Φ(jω), then taking V(jω) as a stable spectral factor  
© National Instruments Corporation  
4-1  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
               
Chapter 4  
Frequency-Weighted Error Reduction  
(so that VV* = Φ) is logical. However, a major use of weighting is in  
controller reduction, which is now described.  
Controller Reduction  
Frequency weighted error reduction becomes particularly important in  
reducing controller dimension. The LQG and Hdesign procedures lead to  
controllers which have order equal to, or roughly equal to, the order of the  
plant. Very often, controllers of much lower order will result in acceptable  
performance, and will be desired on account of their greater simplicity.  
It is almost immediately evident that (unweighted) additive approximation  
of a controller will not necessarily ensure closeness of the behavior of the  
two closed-loop systems formed from the original and reduced order  
controller together with the plant. This behavior is dependent in part on the  
plant, and so one would expect that a procedure for approximating  
controllers ought in some way to reflect the plant. This can be done several  
ways as described in the Controller Robustness Result section. The  
following result is a trivial variant of one in [Vid85] dealing with robustness  
in the face of plant variations.  
Controller Robustness Result  
Suppose that a controller C stabilizes a plant P, and that Cr is a (reduced  
order) approximation to C with the same number of unstable poles. Then  
Cr stabilizes P also provided  
[C( jω) Cr( jω)]P( jω)[I + C( jω)P( jω)]1 < 1  
or  
[I + P( jω)C( jω)]1(P( jω)[C( jω)Cr( jω)]) < 1  
An extrapolation to this thinking [AnM89] suggests that Cr will be a good  
approximation to C (from the viewpoint of some form of stability  
robustness) if  
EIS  
EIS  
4-2  
=
(C Cr)P(I + CP)1  
or  
=
(C Cr)P(I + CP)1  
Xmath Model Reduction Module  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 4  
Frequency-Weighted Error Reduction  
is minimized (and of course is less than 1). Notice that these two error  
measures are like those of Equation 4-1 and Equation 4-2. The fact that the  
plant ought to show up in a good formulation of a controller reduction  
problem is evidenced by the appearance of P in the two weights.  
It is instructive to consider the shape of the weighting matrix or function  
P(Ι + CP)–1. Consider the scalar plant case. In the pass band, |PC| is likely  
to be large, and if this is achieved by having |C| large, then |P(Ι + CP)–1|  
will be (approximately) small. Also outside the plant bandwidth,  
|P(Ι + CP)–1| will be small. This means that it will be most likely to take its  
biggest values at frequencies near the unity gain cross-over frequency. This  
means that the approximation Cr is being forced to be more accurate near  
this frequency than well away from it—a fact very much in accord with  
classical control, where one learns the importance of good loop shaping  
round this frequency.  
The above measures EIS and EOS are advanced after a consideration of  
stability, and the need for its preservation in approximating C by Cr. If one  
takes the viewpoint that the important thing to preserve is the closed-loop  
transfer function matrix, a different error measure arises. With T, Tr  
denoting the closed-loop transfer function matrices,  
T Tr = PC(I + PC) PCr(I + PCr)1  
Then, to a first order approximation in C – Cr, there holds  
T Tr ≈ (I + PC)1P(C Cr)(I + PC)1  
The natural error measure is then  
EM  
=
(I + PC)1P(C Cr)(I + PC)1  
(4-4)  
and this error measure parallels E3 in Equation 4-3. Further refinement  
again is possible. It may well be that closed-loop transfer function matrices  
should be better matched at some frequencies than others; if this weighting  
on the error in the closed-loop transfer function matrices is determined by  
the input spectrum VV* = Φ, then one really wants (T Tr)V to be small,  
so that Equation 4-4 is replaced by  
EMS  
=
(I + PC)1P(C Cr)(I + PC)1V  
© National Instruments Corporation  
4-3  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 4  
Frequency-Weighted Error Reduction  
Most of these ideas are discussed in [Enn84], [AnL89], and [AnM89].  
The function wtbalance( )implements weighted reduction, with five  
choices of error measure, namely EIS, EOS, EM, EMS, and E1 with arbitrary  
V(jω). The first four are specifically for controller reduction, whereas the  
last is not aimed specifically at this situation.  
Several features of the algorithms are:  
Only the stable part of C is really reduced; the unstable part is copied  
exactly into Cr.  
A modification of balanced realization truncation underpins the  
algorithms, namely (frequency) weighted balanced truncation,  
although to avoid numerical problems, the actual construction of  
a frequency weighted balanced realization of C is avoided.  
Frequency weighted Hankel singular values can be computed,  
and although no error bound formula is available (in contrast to the  
unweighted problem), generally speaking there is little damage done in  
reducing by a number of states equal to the number of (relatively) small  
Hankel singular values.  
The error measures themselves deserve certain comments:  
The two stability based measures EIS and EOS are derived from a  
sufficiency condition for stability, rather than a necessity and  
sufficiency condition, and so capture stability a little crudely.  
For any constant nonsingular N, the error measure EIS can be replaced  
N(C Cr)P(I + CP)1N1  
by  
and the robustness result remains  
valid. Use of an N may improve or worsen the quality of the  
approximation.  
Having T Tr small normally ensures closeness of the closed-loop  
impulse and step responses.  
In classical control especially, constraints on the loop gain can be  
imposed (Minimum value of gain in one band, maximum value of gain  
in another band, for example). None of the methods presented directly  
addresses the task of retaining satisfaction of these constraints after  
reduction of a high order acceptable controller. However, judicious use  
of a weight V can assist. Suppose that above the closed-loop bandwidth  
there is an overbound constraint on the loop gain, which is violated  
when a controller reduction is performed (but not with the original  
controller). At these frequencies, roughly PC and PCr are small, so that  
T Tr = P(C Cr). Introduction of a weight V in EMS penalizing  
frequencies in the region in question will evidently encourage PCr to  
better track PC.  
Xmath Model Reduction Module  
4-4  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 4  
Frequency-Weighted Error Reduction  
Fractional Representations  
The treatment of jω-axis or right half plane poles in the above schemes is  
crude: they are simply copied into the reduced order controller. A different  
approach comes when one uses a so-called matrix fraction description  
(MFD) to represent the controller, and controller reduction procedures  
based on these representations (only for continuous-time) are found in  
fracred( ). Consider first a scalar controller C(s) = n(s) ⁄ d(s). One  
can take a stable polynomial d(s) of the same degree as d, and then  
represent the controller as a ratio of two stable transfer functions, namely  
1  
n(s)  
d(s)  
d(s)  
d(s)  
---------  
---------  
Now n d is the numerator, and d d the denominator. Write d d as  
1 + e d. Then we have the equivalence shown in Figure 4-1.  
e
--  
d
n
--  
d
C(s)  
Figure 4-1. Controller Representation Through Stable Fractions  
Evidently, C(s) can be formed by completing the following steps:  
1. Construction of the one-input, two-output stable transfer function  
matrix  
n d  
G =  
e d  
(which has order equal to that of d or d).  
2. Interconnection through negative feedback of the second output to the  
single input.  
These observations motivate the reduction procedure:  
Reduce G to Gr; notice that G is stable. Let Gr be  
nr dr  
G =  
er dr  
© National Instruments Corporation  
4-5  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
Chapter 4  
Frequency-Weighted Error Reduction  
Form the reduced controller by interconnecting using negative  
feedback the second output of Gr to the input, that is, set  
nr  
Cr(s) = ---------------  
dr + er  
Nothing has been said as to how d should be chosen—and the end result  
of the reduction, Cr(s), depends on dr. Nor has the reduction procedure  
been specified.  
When C(s) has been designed to combine a state estimator with a  
stabilizing feedback law, it turns out that there is a natural choice for d(s).  
As for the reduction procedure, one possibility is to use a weight based  
on the spectrum of the input signals to G—and in case C(s) has been  
determined by an LQG optimal design, this spectrum turns out to be white,  
that is, independent of frequency, so that no weight (apart perhaps from  
scaling) is needed. A second possibility is to use a weight based on a  
stability robustness measure. These points are now discussed in more  
detail.  
To understand the construction of a natural fractional representation for  
C(s), suppose that P(s) = C(sI A)1B and let KR, KE be matrices such  
that A BKR and A KEC are stable. The controller  
·
ˆ
ˆ
ˆ
x = Ax + Bu KE(Cx y)  
ˆ
u = –KRx  
ˆ
generates an estimate –KRx of the feedback control –KRx . The controller  
can be represented as a series compensator  
·
ˆ
ˆ
ˆ
ˆ
x = Ax + BKRx KECx + KEy  
ˆ
u = –KRx  
(with compensator input y and output u). Allowing for connection with  
negative feedback, the compensator transfer function matrix is:  
C(s) = KR(sI A + BKR + KEC)1KE  
Xmath Model Reduction Module  
4-6  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 4  
Frequency-Weighted Error Reduction  
Matrix algebra shows that C(s) can be described through a left or right  
matrix fraction description  
C(s) = DL1(s)NL(s) = NR(s)DR1(s)  
with DL, and related values, all stable transfer function matrices.  
In particular:  
DL = I + KR(sI A + KEC)1B  
NL = KR(sI A + KEC)1KE  
NR = KR(sI A + BKR)1KE  
DR = I + C(sI A + BKR)1KE  
For matrix C(s), the left and right matrix fraction descriptions are distinct  
entities. It is the right MFD which corresponds to Figure 4-1; refer to  
Figure 4-2.  
C
-
+
-
1
s
+
+
KE  
KR  
P(s)  
--  
+
A BKR  
+
C(s)  
P(s)  
-
Figure 4-2. C(s) Implemented to Display Right MFD Representation  
© National Instruments Corporation  
4-7  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 4  
Frequency-Weighted Error Reduction  
The left MFD corresponds to the setup of Figure 4-3.  
+
-
P(s)  
Kr(sI – A +KEC)–1  
B
KE  
Figure 4-3. C(s) Implemented to Display Left MFD Representation  
The setup of Figure 4-2 suggests approximation of:  
Kr  
G(s) =  
(sI A + BKr)1KE  
C
whereas that of Figure 4-3 suggests approximation of:  
H(s) = KR(sI A + KEC)1  
B KE  
In the LQG optimal case, the signal driving KE in Figure 4-2 is white noise  
(the innovations process); this motivates the possibility of using no  
frequency dependent weighting in approximating G(s) [but observe that  
after approximating, the signal will no longer be white noise, so that  
argument is questionable]. Simple appeal to duality motivates using no  
frequency dependent weighting for H(s). These are two of the options  
offered by fracred( ).  
Two more fracred( )options depend on examining stability robustness  
(the options are duals of one another). From the stability point of view, the  
ˆ
set-up of Figure 4-3 is identical to that of Figure 4-4, with  
.
P =  
P I  
Xmath Model Reduction Module  
4-8  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 4  
Frequency-Weighted Error Reduction  
-
+
-
KR  
C
(sI – A +BKR)–1KE  
P(s) I  
+
+
ˆ
C(s)  
P(s)  
-
Figure 4-4. Redrawn; Individual Signal Paths as Vector Paths  
It is possible to verify that  
1 ˆ  
(I + PG) P = [CsI A + KEC B  
1  
ˆ
I C(sI (A + KEC)1KE)]  
ˆ
1 ˆ  
and accordingly the output weight (I + PG) P = W can be used in an  
error measure W(G Gr) . It turns out that the calculations for frequency  
weighted balanced truncation of G and subsequent construction of Cr(s) are  
exceptionally easy using this weight.  
The second fracred( )option is the dual of this. The error measure is  
(H Hr)V where:  
I KR(sI A + BKR)1B  
C(sI A + BKR)1B  
V =  
It is possible to argue heuristically the relevance of these error measures  
from a second point of view. It turns out that:  
DL NL V1 NR  
I 0  
=
W1 W2 V2 Dr  
0 I  
© National Instruments Corporation  
4-9  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 4  
Frequency-Weighted Error Reduction  
(Here, the Wi and Vi are submatrices of W,V.) Evidently,  
NR  
DR  
V = I  
and  
W
= I  
DL NL  
Some manipulation shows that trying to preserve these identities after  
approximation of DL, NL or NR, DR suggests use of the error measures  
W(G Gr)  
[LAL90].  
and (H Hr)V . For further details, refer to  
and  
[AnM89]  
In all four fracred( )options, it is possible to construct (weighted)  
Hankel singular values, and to use them as a guide to the likely quality of  
approximation. The patterns tend to be different for the four options.  
The fracred( )options are normally different in outcome from the  
wtbalance( )options. However, if the controller has been designed  
by the loop transfer recovery method and is stable, then one of the  
fracred( )options is essentially the same as one of the wtbalance( )  
options, refer to [LiA90].  
More precisely, if the LTR design is performed with input noise or process  
noise weighting tending to infinity, reduction with fracred( )and  
type="left stab", which uses the error measure (H Hr)V , leads to  
effectively the same reduction as wtbalance( )with the type="input  
stab". If the LTR design is performed with state or output weighting  
tending to infinity (in the index determining the state feedback law),  
reduction with fracred( )and type="right stab"using the error  
measure W(G Gr) leads to effectively the same reduction as  
wtbalance( )with type="output stab".  
wtbalance( )  
[SysCR,SysCLR,HSV]=wtbalance(Sys,SysC,type,{nscr,SysV})  
The wtbalance( )function calculates a frequency weighted balanced  
truncation of a system.  
wtbalance( )has two separate uses:  
Reduce the order of a controller C(s) located in a stable closed-loop,  
with the plant P(s) known. Frequency-weighted balanced truncation is  
used, with the weights involving P(s) and being calculated in a  
predominantly standard way.  
Xmath Model Reduction Module  
4-10  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 4  
Frequency-Weighted Error Reduction  
Reduce the order of a transfer function matrix C(s) through  
frequency-weighted balanced truncation, a stable frequency weight  
V(s) being prescribed.  
The syntax is more accented towards the first use. For the second use,  
the user should set S = 0, NS = 0. This results in (automatically)  
SCLR = NSCLR = 0. The user will also select the type="input  
spec".  
Let Cr(s) be the reduced order approximation of C(s) which is being  
sought. Its order is either specified in advance, or the user responds to  
a prompt after presentation of the weighted Hankel singular values.  
Then the different types concentrate on (approximately) minimizing  
certain error measures, through frequency weighted balanced  
truncation. These are shown in Table 4-1.  
Table 4-1. Types versus Error Measures  
Type  
Error Measure  
[C Cr]P[I + CP]1  
[I + PC]1P[C Cr]  
[I + PC]1P[C Cr][I + PC]1  
"input stab"  
"output stab"  
"match"  
[I + PC]1P[C Cr][I + PC]1V  
"match spec"  
"input spec"  
[C Cr]V  
These error measures have certain interpretations, as shown in Table 4-2.  
In case C(s) is not a compensator in a closed-loop and the error measure  
V(jω)[C(jω) Cr(jω)]  
is of interest, you can work with type="input spec"and C', V' in lieu  
of C and V.  
There is no restriction on the stability of C(s) [or indeed of P(s)] in the  
algorithm, though if C(s) is a controller the closed-loop must be stabilizing.  
Also, V(s) must be stable. Hence all weights (on the left or right of  
C(jω) – Cr(jω) in the error measures) will be stable. The algorithm,  
however, treats unstable C(s) in a special way, by reducing only the stable  
part of C(s) (under additive decomposition) and copying the unstable part  
into Cr(s).  
© National Instruments Corporation  
4-11  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 4  
Frequency-Weighted Error Reduction  
This rather crude approach to the handling of the unstable part of a  
controller is avoided in fracred( ), which provides an alternative to  
wtbalance( )for controller reduction, at least for an important family  
of controllers.  
Table 4-2. Error Measure Interpretation for wtbalance  
Type  
Error Measure Interpretations  
"input stab"  
"output stab"  
"match"  
A stability robustness argument, based on breaking the loop at the controller  
output, indicates that if C is stabilizing for P and the error measure is less  
than 1, then Cr is stabilizing for P. The smaller the error measure is, the  
greater the stability robustness.  
A similar stability robustness argument, but based on breaking the loop  
at the controller input, indicates that if C is stabilizing for P and the error  
measure is less that 1, then Cr is stabilizing for P. The smaller the error  
measure is, the greater the stability robustness.  
If T = PC(I + PC)–1 and Tr = PCr(I + PCr)–1 are the two closed-loop transfer  
function matrices, then T Tr to first order in C Cr is given by  
(I + PC)–1P[CrC][I + PC]–1, so that the error measure looks at matching  
of the closed-loop transfer function matrix.  
"match spec"  
It may be important to match closed-loop transfer function matrices more  
at certain frequencies than others; frequency weighting is achieved by  
introducing V(s). Frequencies corresponding to larger values of |V(jω)| or  
V(jω)V (jω) will be the frequencies at which T(jω) and Tr(jω) should have  
*
smaller error.  
"input spec"  
This is the one error measure that is not associated with a plant, or  
closed-loop of some kind. It simply allows the user to emphasize certain  
frequencies in the reduction procedures.  
Algorithm  
The major steps of the algorithm are as follows:  
1. Check dimension, syntax, stability of SysV, closed-loop stability, and  
decomposition of C(s) into the sum of a stable part (poles in Re[s] < 0)  
and unstable part (poles in Re[s] 0); stable( )is used for this  
purpose.  
2. Compute input (right) weight and/or output (left) weight as appropriate  
for the specified type.  
Xmath Model Reduction Module  
4-12  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 4  
Frequency-Weighted Error Reduction  
3. Compute weighted Hankel Singular Values σi (described in more  
detail later). If the order of Cr(s) is not specified a priori, it must be  
input at this time. Certain values may be flagged as unacceptable for  
various reasons. In particular nscr cannot be chosen so that  
σnscr = σnscr + 1.  
4. Execute reduction step on stable part of C(s), based on a modification  
of redschur( )to accommodate frequency weighting, and yielding  
stable part of Cr(s).  
5. Compute Cr(s) by adding unstable part of C(s) to stable part of Cr(s).  
6. Check closed-loop stability with Cr(s) introduced in place of C(s),  
at least in case C(s) is a compensator.  
More details of steps 3 and 4, will be given for the case when there is an  
input weight only. The case when there is an output weight only is almost  
the same, and the case when both weights are present is very similar, refer  
to [Enn84a] for a treatment. Let  
C(s) = Dc + Cc(sI Ac)1Bc  
WS(s) = Dw + Cw(sI Aw)1Bw  
be a stable transfer function matrix to be reduced and its stable weight.  
Thus, W(s) may be P(I + CP)–1, corresponding to "input stab", and will  
thus have been calculated in step 2; or it maybe an independently specified  
stable V(s). Then  
1  
sI Ac BcCw  
BcDw  
Bw  
Cs(s)W(s) = DcDw +  
Cc DcCw  
0
sI Aw  
The controllability grammian P satisfying  
Ac′  
CwBC Aw  
0
BcDw  
Bw  
Ac BcCw  
DwBc Bw  
P
+
P +  
= 0  
0
Aw  
is written as  
Pcc Pcw  
Pcw Pww  
P =  
© National Instruments Corporation  
4-13  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 4  
Frequency-Weighted Error Reduction  
and the observability grammian Q, defined in the obvious way, is written as  
Qcc Qcw  
Qcw Qww  
Q =  
It is trivial to verify that QccAc + AcQcc = –CC so that Qcc is the  
c
c
observability gramian of Cs(s) alone, as well as a submatrix of Q.  
The weighted Hankel singular values of Cs(s) are the square roots of the  
eigenvalues of PccQcc. They differ from the usual or unweighted Hankel  
singular values because Pcc is not the controllability gramian of Cs(s) but  
rather a weighted controllability gramian. The usual controllability  
gramian can be regarded as E[x x] when Cs(s) is excited by white noise.  
c
c
The weighted controllability gramian is still  
excited by colored noise, that is, the output of the shaping filter W(s), which  
, but now C (s) is  
E[xcxc]  
s
is excited by white noise.  
Small weighted Hankel singular values are a pointer to the possibility  
of eliminating states from Cs(s) without incurring a large error in  
[C(jω) Cr(jω)]W(jω) . No error bound formula is known, however.  
The actual reduction procedure is virtually the same as that of  
redschur( ), except that Pcc is used. Thus Schur decompositions of  
PccQcc are formed with the eigenvalues in ascending and descending order  
VAPccQccVA = Sasc  
VDPccQccVD = Sdes  
The maximum order permitted is the number of nonzero eigenvalues of  
PccQcc that are larger than ε.  
The matrices VA, VD are orthogonal and Sasc and Sdes are upper triangular.  
Next, submatrices are obtained as follows:  
0
Inscr  
Vlbig = VA  
Vrbig = VD  
Inscr  
0
and then a singular value decomposition is formed:  
UebigSebigVebig = VlbigVrbig  
Xmath Model Reduction Module  
4-14  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 4  
Frequency-Weighted Error Reduction  
From these quantities the transformation matrices used for calculating  
Csr(s), the stable part of Cr(s), are defined  
1 2  
Slbig = VlbigVebig  
S
ebig  
1 2  
Srbig = VrbigVebig  
S
ebig  
and then  
ACR = SlbigACSrbig BCR = SlbigBC  
ACR = CCSrbig BCR = DC  
Just as in unweighted balanced truncation, the reduced order transfer  
function matrix is guaranteed stable, the same is guaranteed to be true in  
weighted balanced truncation when either a left (output) weight or a right  
(input) weight is used. It is suspected to be true when both input and output  
weights are present. The overall algorithm is not, however, at risk in this  
case, since it is stability of the closed-loop system which is the key issue of  
concern, (except for type="input spec", but here there is only a single  
weight, and so the theory guarantees preservation of stability).  
Related Functions  
balance(), redschur(), stable(), fracred()  
fracred( )  
[SysCR,HSV] = fracred(Sys,Kr,Ke,type,{nscr,Qyy})  
The fracred( )function uses fractional representations to calculate a  
reduction of a continuous-time compensator comprising a state estimator  
with state feedback law.  
Restrictions  
1. The closed-loop system (SCLR,NSCLR) is calculated from  
sysol=scr*sys  
# open loop system  
syscl=feedback(sysol)  
# closed loop system  
2. Initial state values, state names, and input and output names are not  
considered by fracred( ).  
© National Instruments Corporation  
4-15  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 4  
Frequency-Weighted Error Reduction  
3. Only continuous systems are accepted; for discrete systems use  
makecontinuous( )before calling bst( ), then discretize the  
result.  
Sys=fracred(makecontinuous(SysD));  
SysD=discretize(Sys);  
Defining and Reducing a Controller  
Suppose P(s) = C(sI A)–1B and A BKR and A KEC are stable (where  
KR is a stabilizing state feedback gain and KE a stabilizing observer gain).  
A controller for the plant P(s) can be defined by  
·
ˆ
ˆ
ˆ
x = Ax + Bu KE(Cx y)  
ˆ
u = –KRx  
(with u the plant input and y the plant output). The associated series  
compensator under unity negative feedback is  
C(s) = KR(sI A + BKR + KEC)1KE  
and this may be written as a left or right MFD as follows:  
C(s) = [I + KR(sI A + KEC)1B]1KR(sI A + KEC)1KE  
(4-5)  
(4-6)  
C(s) = KR(sI A + BKR)1KE[I + C(sI A + BKR)1KE]1  
The reduction procedures "right perf"and "left perf"have similar  
rationales. We shall describe "right perf", refer to [AnM89] and  
[LiA86]. The first rationale involves observing that to reduce C(s), one  
might as well reduce its numerator and denominator simultaneously, and  
then form a new fraction Cr(s) of lower order than C(s).  
This amounts to reducing  
(sI A + BKR)1KE  
(4-7)  
KR  
E(s) =  
C
Xmath Model Reduction Module  
4-16  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 4  
Frequency-Weighted Error Reduction  
to, for example,  
KR  
ˆ
Er(s) =  
(sI A)KE  
C
through, for example, balanced truncation, and then defining:  
1  
1  
1  
ˆ
ˆ
C(s) = KR(sI A) KE[I + C(sI A) KE]  
1  
ˆ
= KR(sI A + KEC) KE  
For the second rationale, consider Figure 4-5.  
C
-
+
-
1
--  
+
+
KE  
KR  
P(s)  
s
+
X
A BKR  
Figure 4-5. Internal Structure of Controller  
Recognize that the controller C(s) (shown within the hazy rectangle in  
Figure 4-5) can be constructed by implementing  
KR(sI A + BKR)1KE  
and  
C(sI A + BKR)1KE  
and then applying an interconnection rule (connect the output of the second  
transfer function matrix back to the input at point X in Figure 4-5).  
© National Instruments Corporation  
4-17  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 4  
Frequency-Weighted Error Reduction  
Controller reduction proceeds by implementing the same connection rule  
but on reduced versions of the two transfer function matrices.  
When KE has been defined through Kalman filtering considerations, the  
spectrum of the signal driving KE in Figure 4-5 is white, with intensity Qyy.  
It follows that to reflect in the multiple input case the different intensities  
on the different scalar inputs, it is advisable to introduce at some stage a  
weight Q1yy2 into the reduction process.  
Algorithm  
After preliminary checks, the algorithm steps are:  
1. Form the observability and weighted (through Qyy) controllability  
grammians of E(s) in Equation 4-7 by  
P(A BKR)′ + (A BKR)P = –KEQyyKE′  
(4-8)  
(4-9)  
Q(A BKR) + (A BKR)′Q = –KRKR CC  
2. Compute the square roots of the eigenvalues of PQ (Hankel singular  
values of the fractional representation of Equation 4-5). The maximum  
order permitted is the number of nonzero eigenvalues of PQ that are  
larger than ε.  
3. Introduce the order of the reduced-order controller, possibly by  
displaying the Hankel singular values (HSVs) to the user. Broadly  
speaking, one can throw away small HSVs but not large ones.  
4. Using redschur( )-type calculations, find a state-variable  
description of Er(s). This means that Er(s) is the transfer function  
matrix of a truncation of a balanced realization of E(s), but the  
redschur( )type calculations avoid the possibly numerically  
difficult step of balancing the initially known realization of E(s).  
Suppose that:  
ˆ
A = Slbig(A BKR)Srbig, KE = SlbigKE  
5. Define the reduced order controller Cr(s) by  
ACR = Slbig(A BKR KEC)Srbig  
(4-10)  
Cr(s) = CCR(sI ACR)1BCR  
Xmath Model Reduction Module  
4-18  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
         
Chapter 4  
Frequency-Weighted Error Reduction  
6. Check the stability of the closed-loop system with Cr(s). When the  
type="left perf"is specified, one works with  
E(s) = KR(sI A + KEC)1  
(4-11)  
B KE  
which is formed from the numerator and denominator of the MFD  
in Equation 4-5. The grammian equations (Equation 4-8 and  
Equation 4-9) are replaced by  
P(A KEC)′ + (A KEC)P = – BBKEKE′  
Q(A KEC) + (A KEC)′Q = –KRKR  
redschur( )-type calculations are used to reduce E(s) and Equation 4-10  
again yields the reduced-order controller. Notice that the HSVs obtained  
from Equation 4-10 or the left MFD (Equation 4-5) of C(s) will in general  
be quite different from those coming from the right MFD (Equation 4-6). It  
may be possible to reduce much more with the left MFD than with the right  
MFD (or vice-versa) before closed-loop stability is lost.  
As noted in the fracred( )input listing, type="left stab"and  
"right stab"focus on a stability robustness measure, in conjunction  
with Equation 4-5 and Equation 4-6, respectively. Leaving aside for the  
moment the explanation, the key differences in the algorithm computations  
lie solely in the calculation of the grammians P and Q. For type="left  
stab", these are given by  
P(A BKR)′ + (A BKR)P = –BB′  
Q(A KEC) + (A KEC)′Q = –KRKR  
and for "right stab",  
P(A BKR)′ + (A BKR)P = –KEKE′  
(4-12)  
(4-13)  
Q(A KEC) + (A KEC)′Q = –CC  
© National Instruments Corporation  
4-19  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 4  
Frequency-Weighted Error Reduction  
Additional Background  
A discussion of the stability robustness measure can be found in [AnM89]  
and [LAL90]. The idea can be understood with reference to the transfer  
functions E(s) and Er(s) used in discussing type="right perf". It is  
possible to argue (through block diagram manipulation) that  
C(s) stabilizes P(s) when E(s) stabilizes (as a series compensator) with  
ˆ
unity negative feedback  
.
P(s) =  
P(s) I  
Er(s) also will stabilize [P(s)I], and then Cr(s) will stabilize P(s),  
provided  
C( jωI A + KEC)1B I C(jωI A + KEC)1KE  
< 1 (4-14)  
[E( jω) Ee(jω)]  
Accordingly, it makes sense to try to reduce E by frequency-weighted  
balanced truncation. When this is done, the controllability grammian for  
E(s) remains unaltered, while the observability grammian is altered. (Hence  
Equation 4-5, at least with Qyy = I, and Equation 4-12 are the same while  
Equation 4-6 and Equation 4-13 are quite different.) The calculations  
leading to Equation 4-13 are set out in [LAL90].  
The argument for type="left perf"is dual. Another insight into  
Equation 4-14 is provided by relations set out in [NJB84]. There, it is  
established (in a somewhat broader context) that  
{C(jωI A + KEC)1B  
I C( jωI A + KEC)1KE}  
Kr(sI A + BKR)1KE  
I + C( jωI A + BKR)1KE  
×
= I  
The left matrix is the weighting matrix in Equation 4-14; the right matrix is  
the numerator of C(jω) stacked on the denominator, or alternatively  
0
E(jω) +  
I
This formula then suggests the desirability of retaining the weight in the  
approximation of E(jω) by Er(jω).  
Xmath Model Reduction Module  
4-20  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 4  
Frequency-Weighted Error Reduction  
The four schemes all produce different HSVs; it follows that it may be  
prudent to try all four schemes for a particular controller reduction. Recall  
again that their relative sizes are only a guide as to what can be thrown away  
without incurring much error. There is no simple rule to indicate which of  
the four schemes will be the most effective at controller reduction.  
Two rough rules can, however, be formulated.  
Problems with instability through reduction to too low a controller  
order are more likely with "left perf"and "right perf"than  
"left stab"or "right stab".  
If the controller has been designed using the loop transfer recovery  
idea, "left stab"will probably be attractive if the input noise  
covariance is very large, and "right stab"will probably be  
attractive if the output weighting in the performance index is very  
large, [LiA90]. The reduced controllers will then actually be very  
similar to those obtained using wtbalance( )with the option  
"input stab"in the first case and "output stab"in the second  
case.  
One example gives the HSVs summarized in Table 4-3 for an eighth order  
controller.  
Table 4-3. HSVs for an Eighth Order Controller  
1
2
3
4
5
6
7
8
right perf  
left perf  
right stab  
left stab  
.0339  
.0164  
4.8742  
.7278  
1.317  
.0128  
3.8457  
.1123  
1.1269  
.0102  
3.7813  
.0783  
1.0862  
.0040  
1.2255  
.0242  
.9638  
.0037  
1.1750  
.0181  
.5846  
.0000  
.5055  
.0107  
.5646  
.0000  
.0413  
.0099  
.3144  
4.9075  
3.3081  
1.3914  
The most attractive candidate for reducing to second order is right stab.  
This is because the HSVs being discarded (columns 3 to 8) are smaller  
relative to those being retained (columns 1 and 2) for right stabthan for  
the other three candidates.  
Note The relative values count, not the absolute values.  
Related Functions  
redschur(), wtbalance()  
© National Instruments Corporation  
4-21  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
5
Utilities  
This chapter describes three utility functions: hankelsv( ), stable( ),  
and compare( ).  
The background to hankelsv( ), which calculates Hankel singular  
values, was presented in Chapter 1, Introduction. Hankel singular values  
are also calculated in other functions, sometimes by other procedures.  
A comparison of the procedures is given in the Hankel Singular Values  
section. The function compare( )serves to facilitate the comparisons  
of an unreduced and a reduced system, from various points of views.  
The function stable( )is used to separate (additively) a system into its  
stable and unstable parts, that is, given G(s), the function determines Gs(s)  
and Gu(s), the first with all poles in Re[s] < 0, the second with all poles in  
Re[s] 0, such that  
G(s) = Gs(s) + Gu(s)  
The function is used within some of the other functions of the Model  
Reduction Module. It should also be used when reduction of an unstable  
G(s) is contemplated. The normal reduction functions, for example,  
balmoore( )or redschur( ), require stability of the transfer function  
matrix G(s) being reduced. If G(s) is unstable, stable( )should be used  
to generate Gs(s) and Gu(s); reduction of Gs(s) should be performed, and  
then Gu(s) added to the outcome using the +operator, to yield the desired  
reduction of Gs(s).  
hankelsv( )  
[HSV,Wc,Wo] = hankelsv(Sys,{noplot})  
The hankelsv( )function computes the Hankel Singular Values of a  
stable system (continuous or discrete) and displays them in a bar plot.  
© National Instruments Corporation  
5-1  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
         
Chapter 5  
Utilities  
The gramian matrices are defined by solving the equations (in continuous  
time)  
AWc + WcA= –BB′  
WoA + AW0 = –CC  
and, in discrete time  
Wc AWcA= BB′  
Wo AWoA = CC  
The computations are effected with lyapunov( )and stability is checked,  
which is time-consuming. The Hankel singular values are the square roots  
of the eigenvalues of the product.  
Related Functions  
lyapunov(), dlyapunov()  
stable( )  
[SysS,SysU] = stable(Sys,{tol})  
The stable( )function decomposes Sysinto its stable (SysS) and  
unstable(SysU) parts, such that Sys=SysS+SysU.  
Continuous systems have unstable poles if real parts > –tol.  
Discrete systems have unstable poles if magnitudes > 1-tol.  
The direct term (D matrix) is included in SysS.  
If Syshas poles clustered near -tol(or 1-tol), then SysSand SysU  
might be ill-conditioned. To avoid this problem choose tolto a value  
that is not close to the majority of poles.  
Algorithm  
The algorithm begins by transforming the A matrix to Schur form, and  
counting the number of stable and unstable eigenvalues, together with  
those for which classification is doubtful. Stable eigenvalues are those  
in either of the following:  
Re[s] < 0 (continuous time)  
|z| < 1 (discrete time)  
Xmath Model Reduction Module  
5-2  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
         
Chapter 5  
Utilities  
Doubtful ones are those for which the real part of the eigenvalue has  
magnitude less than or equal to tolfor continuous-time, or eigenvalue  
magnitude within the following range for discrete time:  
1 tol, 1 + tol  
A warning is given if doubtful eigenvalues exist.  
The algorithm then computes a real ordered Schur decomposition of A  
so that after transformation  
AS ASU  
A =  
0 AU  
where the eigenvalues of AS and AU are respectively stable and unstable.  
A matrix X satisfying –ASX + XAU + ASU = 0 is then determined by calling  
the algorithm sylvester( ). The eigenvalue properties of AS and AU  
guarantee that X exists. If doubtful eigenvalues are present, they are  
assigned to the unstable part of Sys. In this circumstance you get the  
message,  
The system has poles near, or upon, the jw-axis  
for continuous systems, and the following for discrete systems:  
The system has poles near the unit circle.  
Note If A has eigenvalues clustered near -tol(1–tolin discrete-time), then X is likely  
to be ill-conditioned and consequently SysSand SysUwill also be ill-conditioned. (For  
example, the B matrix of SysScould contain very small values, while the C matrix could  
contain large values. In this case, SysSwould be very weakly controllable and very  
strongly observable. This will cause problems when gramians and Hankel singular values  
are calculated.) To avoid this problem, change tolto a value that is not close to the  
majority of eigenvalues.  
A further transformation of A is constructed using X:  
I X  
I X  
A →  
A
0 I  
0 I  
AS  
0
=
0 AU  
© National Instruments Corporation  
5-3  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 5  
Utilities  
After this last transformation, and with  
BS  
B =  
C = [CSCU]  
BU  
it follows that  
SysS = [AS,AS;CSD]  
and  
SysU = [AU,BU;CU0]  
By combining the transformation yielding the real ordered Schur form for  
A with the transformation defined using X, the overall transformation T is  
readily identified. In case all eigenvalues of A are stable or all are unstable,  
this is flagged, and T = I.  
stable( )can be combined with a reduction algorithm such as  
redschur( )or balmoore( )to reduce the order of a system with some  
unstable and some stable poles. One uses stable( )to separate the stable  
and unstable parts, and then, for example, reduces the stable part with  
redschur( ); the reduced stable part is added to the original unstable part  
to provide the desired system reduction.  
Related Functions  
sylvester(), schur(), redschur(), balmoore()  
compare( )  
[respdiff] = compare(Sys,SysRed,FTvec,{Fmin,Fmax,npts,radians,type})  
The compare( )function provides a number of different graphical tests  
which can be used to compare two state-space system implementations.  
compare( )can be used as a tool for evaluating a reduced-order system  
by comparing it with the original full-order system from which it was  
obtained. However, it can be used for more general comparisons as well,  
such as examining the results of different discretization or identification  
Xmath Model Reduction Module  
5-4  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
6
Tutorial  
This chapter illustrates a number of the MRM functions and their  
underlying ideas. A plant and full-order controller are defined, and then  
the effects of various reduction algorithms are examined. The data for this  
example is stored in the file mr_disc.xmdin the Xmath demos directory.  
To follow the example, start Xmath, and then select File»Load from the  
Xmath Commands menu, or enter the load command with the file  
specification appropriate to your operating system from the Xmath  
Commands area. For example:  
load "$XMATH/demos/mr_disc"  
Plant and Full-Order Controller  
The plant in question comprises four spinning disks, connected by a  
flexible shaft. A motor applies torque to the third disk, and the output  
variable of interest is the angular displacement of the first disk. The plant  
transfer function, which is nonminimum phase and has a double pole at the  
origin, is as follows:  
s22ζ0ω0s + ω02 s2ζ1ω1s + ω12  
s + a  
a
----------------------------------- -------------------------------- -----------  
ω02  
ω12  
1
----------------------------------------------------------------------------------------------------------------------------  
G(s) =  
4s2 s22ζ2ω2s + ω22 s22ζ3ω3s + ω23 s22ζ4ω4s + ω42  
----------------------------------- ----------------------------------- -----------------------------------  
ω22  
ω32  
ω42  
with:  
ζ =0.02 ω =1  
0
0
ζ =-0.4 ω =5.65  
1
1
ζ =0.02 ω =0.765  
2
2
ζ =0.02 ω =1.41  
3
3
ζ =0.02 ω =1.85  
4
4
a=4.84  
© National Instruments Corporation  
6-1  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
       
Chapter 6  
Tutorial  
A minimal realization in modal coordinates is C(sI A)–1B where:  
0 1 0.015 0.765  
0.028 1.410  
0.04 1.85  
A = diag  
,
,
,
0 0 0.765 0.015 1.410 0.028 1.85 0.04  
0.026  
0.251  
0.033  
0.996  
0.105  
0.261  
0.886  
4.017  
0.145  
0.009  
B =  
C=  
0.001  
0.043  
0.002  
3.604  
0.280  
0.026  
The specifications seek high loop gain at low frequencies (for performance)  
and low loop gain at high frequencies (to guarantee stability in the presence  
of unstructured uncertainty). More specifically, the loop gain has to lie  
outside the shaded region shown in Figure 6-1.  
40 dB/decade  
Frequency (rad/sec)  
0.3  
0.07  
40 dB/decade  
Figure 6-1. Loop Gain Constraints  
Xmath Model Reduction Module  
6-2  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
With a state weighting matrix,  
Q = 1e-3*diag([2,2,80,80,8,8,3,3]);  
R = 1;  
(and unity control weighting), a state-feedback control-gain is determined  
through a linear-quadratic performance index minimization as:  
[Kr,ev] = regulator(sys,Q,R);  
A B × Kr is stable. Next, with an input noise variance matrix Q = WtBBWt  
where,  
Wt = DIAG([0.346, 0.346, 0.024, 0.0240.042, 0.0420.042, 0.042])  
ˆ
and measurement noise covariance matrix R =1, an estimation gain Ke  
(so that A KeC is stable) is determined:  
Qhat = Wt*b*b'*Wt;  
Rhat = 1;  
[Ke,ev] = estimator(sys,Qhat,Rhat,{skipChks});  
The keyword skipChkscircumvents syntax checking in most functions.  
It is used here because we know that Qhatdoes not fulfill positive  
semidefiniteness due to numerics).  
sysc=lqgcomp(sys,Kr,Ke);  
poles(sysc)  
ans (a column vector) =  
-0.296674 + 0.292246 j  
-0.296674 - 0.292246 j  
-0.15095 + 0.765357 j  
-0.15095 - 0.765357 j  
-0.239151 + 1.415  
-0.239151 - 1.415  
j
j
-0.129808 + 1.84093 j  
-0.129808 - 1.84093 j  
The compensator itself is open-loop stable. A brief explanation of how Q  
and Wtare chosen is as follows. First, Qis chosen to ensure that the loop  
gain Kr(jωI A)1B (which would be relevant were the state measurable)  
meets the constraints as far as possible. However, it is not possible to obtain  
a 40 dB per decade roll-off at high frequencies, as LQ design virtually  
always yields a 20 dB per decade roll-off. Second, a loop transfer recovery  
ˆ
Q
approach to the choice of as ρBBfor some large ρ is modified through  
the introduction of the diagonal matrix Wt. The larger entries of Wt, because  
of the modal coordinate system, in effect promote better loop transfer  
© National Instruments Corporation  
6-3  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
recovery at low frequencies; there is consequently a faster roll-off of the  
loop gain at high frequencies than for Kr(jωI A)1B , and this is desired.  
Figure 6-2 displays the (magnitudes of the) plant transfer function, the  
compensator transfer function and the loop gain, as well as the constraints;  
evidently the compensated plant meets the constraints.  
You can enter the following commands to create a plot equivalent to  
Figure 6-2:  
sysol=sys*sysc;  
svals=svplot(sys,w,{radians});  
svalsc=svplot(sysc,w,{radians});  
svalsol=svplot(sysol,w,{radians});  
plot(svals,{x_log,!grid,!ylab,  
line_width=2,hold})  
plot(svalsc,{keep})  
plot(svalsol,{keep})  
f2=plot(wc,constr,{keep,  
legend=["plant","compensator",  
"compensated plant","constraint"]})  
plot({!hold})  
Figure 6-2. Frequency Response for Plant, Compensator, and Compensated Plant  
Xmath Model Reduction Module  
6-4  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
Controller Reduction  
This section contrasts the effect of unweighted and weighted controller  
reduction. Unweighted reduction is at first examined, through  
redschur( )(using balance( )or balmoore( )will give similar  
results). The Hankel singular values of the controller transfer function are  
6.264×10–2 4.901×10–2 2.581×10–2 2.474×10–2  
1.545×10–2 1.335×10–2 9.467×10–3 9.466×10–3  
A reduction to order 2 is attempted. The ending Hankel singular values, that  
is, σ3, σ4, ..., σ8, have a sum that is not particularly small with respect to σ1  
and σ2; this is an indication that problems may arise in the reduction.  
[syscr,hsv] = redschur(sysc,2);  
svalsRol = svplot(sys*syscr,w,{radians});  
plot(svalsol, {keep})  
f3=plot(wc, constr,{keep,!grid,  
legend=["reduced","original","constrained"],  
title="Open-Loop Gain Using redschur()"})  
Figure 6-3. Open-Loop Gain Using redschur  
© National Instruments Corporation  
6-5  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 6  
Tutorial  
Figures 6-3, 6-4, and 6-5 display the outcome of the reduction. The loop  
gain is shown in Figure 6-3. The error near the unity gain crossover  
frequency may not look large, but it is considerably larger than that  
obtained through frequency weighted reduction methods, as described  
later.  
Figure 6-3 also shows the inability to suppress all three plant resonances,  
in contrast to the full-order controller. Two are such as to cause violation  
of the specifications. The closed-loop gains differ by some 4 to 5 dB  
between the full-order and reduced-order controller, in the vicinity of  
0.1 radians per second. The step response has overshoot of 50% as opposed  
to 40% and the ripple persists for longer.  
We use the compare( )function (refer to the compare( ) section of  
Chapter 5, Utilities) to reproduce Figures 6-4 and 6-5. Calculate the  
full-order closed-loop system, then the closed-loop system with the  
reduced-order compensator:  
syscl = feedback(sysol);  
sysolr=sys*syscr;  
sysclr=feedback(sysolr);  
Xmath Model Reduction Module  
6-6  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Chapter 6  
Tutorial  
Generate Figure 6-4:  
compare(syscl,sysclr,w,{radians,type=5})  
f4=plot({keep,legend=["original","reduced"]})  
Figure 6-4. Closed-Loop Gain with redschur  
© National Instruments Corporation  
6-7  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
Generate Figure 6-5:  
tvec=0:(140/99):140;  
compare(syscl,sysclr,tvec,{type=7})  
f5=plot({keep,legend=["original","reduced"]})  
Figure 6-5. Step Response with redschur  
Xmath Model Reduction Module  
6-8  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
ophank( )  
ophank( )is next used to reduce the controller with the results shown in  
Figures 6-6, 6-7, and 6-8.  
Generate Figure 6-6:  
[syscr,sysu,hsv]=ophank(sysc,2);  
svalsrol = svplot(sys*syscr,w,{radians});  
plot(svalsol, {keep})  
f6=plot(wc, constr, {keep,!grid,  
title="Open-loop gain using ophank()"})  
Figure 6-6. Open-Loop Gain Using ophank  
© National Instruments Corporation  
6-9  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 6  
Tutorial  
Generate Figure 6-7:  
syscl = feedback(sysol);  
sysolr=sys*syscr;  
sysclr=feedback(sysolr);  
compare(syscl,sysclr,w,{radians,type=5})  
f7=plot({keep,legend=["original","reduced"]})  
Figure 6-7. Closed-Loop Gain with ophank  
Xmath Model Reduction Module  
6-10  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
Generate Figure 6-8:  
tvec=0:(140/99):140;  
compare(syscl,sysclr,tvec,{type=7})  
f8=plot({keep,legend=["original","reduced"]})  
Figure 6-8. Step Response with ophank  
The open-loop gain, closed-loop gain and step response are all inferior to  
those obtained with redschur( ). This emphasizes the point that one  
cannot automatically assume that, because the error bound formula for  
ophank( )is more attractive than that for redschur( ), the error itself  
will be better for ophank( ).  
© National Instruments Corporation  
6-11  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
wtbalance  
The next command examined is wtbalancewith the option "match".  
[syscr,ysclr,hsv] = wtbalance(sys,sysc,"match",2)  
Recall that this command should promote matching of closed-loop transfer  
functions. The weighted Hankel singular values are:  
1.486 4.513 × 10–1 8.420 × 10–2 5.869 × 1–2  
1.999 × 10–2 1.382 × 10–2 7.198 × 10–3 6.336 × 10–3  
The relative magnitudes suggest that reduction to order 2 will produce less  
of an approximation error here (in the closed-loop transfer function) than a  
reduction to this order through redschur( )or ophank( )(where the  
implicit criterion is the unweighted error in approximating the controller  
transfer function). Examination of Figures 6-9, 6-10, and 6-11 reveals that  
far better approximation is now obtained.  
Violation of the specification is to be observed in the open-loop gain.  
Notice though that:  
The error measure for wtbalancedoes not reflect the open-loop gain;  
it reflects the closed-loop gain.  
While the error in dB looks large, as an absolute value it is not  
extremely so; wtbalanceworks with additive, not multiplicative  
error.  
Hence, it cannot be concluded that the algorithm is not working. Use of the  
option "match spec"with wtbalancemight be conjectured as a device  
for reducing the violation of the specification: one could introduce a weight  
V(jw) emphasizing frequencies from 0.1 radians per second to 5 radians per  
second.  
For example,  
(s + 0.1)(s + 10)  
V(jω) = ----------------------------------------  
(s + 1)(s + 1.4)  
This would tend to force the closed-loop transfer functions derived from  
the full-order and reduced controller to match better over this range;  
because their absolute value is small there, they are approximately equal  
to the open-loop gains which, accordingly, may be close. The flaw in this  
reasoning is that a second-order controller, with four independent  
parameters only, can only do so much, and the totality of designer demands  
cannot be fully met.  
Xmath Model Reduction Module  
6-12  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
The following function calls produce Figure 6-9:  
svalsrol = svplot(sys*syscr,w,{radians})  
plot(svalsol, {keep})  
f9=plot(wc, constr, {keep,!grid,  
legend=["reduced","original","constrained"],  
title="Open-Loop Gain Using wtbalance()"})  
Figure 6-9. Open-Loop Gain with wtbalance  
© National Instruments Corporation  
6-13  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
Generate Figure 6-10:  
syscl = feedback(sysol);  
sysolr=sys*syscr;  
sysclr=feedback(sysolr);  
compare(syscl,sysclr,w,{radians,type=5})  
f10=plot({keep,legend=["original","reduced"]})  
Figure 6-10. Closed-Loop Gain with wtbalance  
Xmath Model Reduction Module  
6-14  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
Generate Figure 6-11:  
tvec=0:(140/99):140;  
compare(syscl,sysclr,tvec,{type=7})  
f11=plot({keep,legend=["original","reduced"]})  
Figure 6-11. Step Response with wtbalance  
Figures 6-9, 6-10, and 6-11 are obtained for wtbalancewith the option  
"input spec". Evidently, there is little difference between this and the  
result with the option "match". One notices marginally better matching in  
the region of interest (0.1 to 5 rad per second) at the expense of matching  
at other frequencies. The weighted Hankel singular values again indicate  
that it is reasonable to seek a second order controller.  
© National Instruments Corporation  
6-15  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
Generate Figure 6-12:  
vtf=poly([-0.1,-10])/poly([-1,-1.4])  
[,sysv]=check(vtf,{ss,convert});  
svalsv = svplot(sysv,w,{radians});  
Figure 6-12. Frequency Response of the Weight V(jω)  
Xmath Model Reduction Module  
6-16  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
Generate Figure 6-13:  
[syscr,sysclr,hsv] = wtbalance(sys,sysc,  
"input spec",2,sysv)  
svalsrol = svplot(sys*syscr,w,{radians})  
plot(svalsol, {keep})  
f13=plot(wc,constr,{keep, !grid,  
legend=["reduced","original","constrained"],  
title="Open-Loop Gain with wtbal(), \"input spec\""})  
Figure 6-13. Open-Loop Gain from wtbalance with "input spec"  
© National Instruments Corporation  
6-17  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
Generate Figure 6-14:  
syscl = feedback(sysol);  
sysolr=sys*syscr;  
sysclr=feedback(sysolr);  
compare(syscl,sysclr,w,{radians,type=5})  
f14=plot({keep,legend=["original","reduced"]})  
Figure 6-14. System Singular Values of wtbalance with "input spec"  
Xmath Model Reduction Module  
6-18  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
Generate Figure 6-15:  
tvec=0:(140/99):140;  
compare(syscl,sysclr,tvec,{type=7})  
f15=plot({keep,legend=["original","reduced"]})  
Figure 6-15. Step Response of wtbalance with "input spec"  
© National Instruments Corporation  
6-19  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
fracred  
fracred, the next command examined, has four options—"right  
stab", "left stab", "right perf", and "left perf".  
The options "left stab", "right perf", and "left perf"all  
produce instability. Given the relative magnitudes of the Hankel singular  
values, this is perhaps not surprising. Figures 6-16, 6-17, and 6-18  
illustrate the results using "right stab".  
Generate Figure 6-16:  
svalsrol = svplot(sys*syscr,w,{radians})  
plot(svalsol, {keep})  
f16=plot(wc,constr,{keep,!grid,  
legend=["reduced","original","constrained"],  
title="Open-Loop Gain Using fracred()"})  
Figure 6-16. Open-Loop Gain Using fracred  
Xmath Model Reduction Module  
6-20  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
   
Chapter 6  
Tutorial  
Generate Figure 6-17:  
syscl = feedback(sysol);  
sysolr=sys*syscr;  
sysclr=feedback(sysolr);  
compare(syscl,sysclr,w,{radians,type=5})  
f17=plot({keep,legend=["original","reduced"]})  
Figure 6-17. Closed-Loop Response with fracred  
© National Instruments Corporation  
6-21  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
Generate Figure 6-18:  
tvec=0:(140/99):140;  
compare(syscl,sysclr,tvec,{type=7})  
f18=plot({keep,legend=["original","reduced"]})  
Figure 6-18. Step Response with fracred  
The end result is comparable to that from wtbalance( )with option  
"match".  
We can create a table to examine the values of the Hankel singular values  
based on different decompositions approaches.  
set precision 3 # Optional:  
set format fixed # we set a smaller precision here so we  
could fit  
# the table in the manual.  
[syscr, hsvrs] = fracred(sys, Kr, Ke, "right stab",2);  
[syscr, hsvls] = fracred(sys, Kr, Ke, "left stab",2);  
[syscr, hsvrp] = fracred(sys, Kr, Ke, "right perf",2);  
[syscr, hsvlp] = fracred(sys, Kr, Ke, "left perf",2);  
Xmath Model Reduction Module  
6-22  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Chapter 6  
Tutorial  
hsvtable = [...  
"right stab:", string(hsvrs');  
"left stab:", string(hsvls');  
"right perf:", string(hsvrp');  
"left perf:", string(hsvlp')]?  
hsvtable (a rectangular matrix of strings) =  
right stab:3.308 0.728 0.112 0.078 0.024 0.018 0.011 0.010  
left stab:1.403 1.331 1.133 1.092 0.965 0.549 0.526 0.313  
right perf:0.034 0.016 0.013 0.010 0.004 0.004 0.000 0.000  
left perf:4.907 4.874 3.846 3.781 1.225 1.175 0.505 0.041  
© National Instruments Corporation  
6-23  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
A
Bibliography  
[AnJ]  
BDO Anderson and B. James, “Algorithm for multiplicative approximation of a stable  
linear system,” in preparation.  
[AnL89]  
[AnM89]  
[BoD87]  
[Enn84]  
BDO Anderson and Y. Liu, “Controller reduction: Concepts and approaches,” IEEE  
Transactions on Automatic Control, Vol. 34, 1989, pp. 802–812.  
BDO Anderson and J. B. Moore, Optimal Control: Linear Quadratic Methods,  
Prentice-Hall Inc., Englewood Cliffs, NJ, 1989.  
S. P. Boyd and J. Doyle, “Comparison of peak and RMS gains for discrete-time systems,”  
System Control Letters, Vol. 9, No. 1, pp. 1–6, 1987.  
D. F. Enns, “Model reduction with balanced realizations: an error bound and a frequency  
weighted generalization,” Proceedings for the 23rd IEEE Conference on Decision and  
Control, Las Vegas, November 1984, pp. 127–132.  
[Enn84a]  
[GCP88]  
D.F. Enns, “Model reduction for control systems design,” PhD Thesis, Dept of  
Aeronautics and Astronautics, Stanford University, CA, USA, 1984.  
K. Glover, R. F. Curtain, and J. R. Partington, “Realisation and approximation of linear  
infinite-dimensional systems with error bounds,” SIAM J Controls and Optimization,  
Vol. 26, 1988, pp. 863–898.  
[Glo84]  
[Glo86]  
[GrA86]  
[GrA89]  
K. Glover, “All optimal Hankel norm approximations of linear multivariable systems and  
their Lerror bounds,” Int J Controls, Vol. 39, 1984, pp. 1115–1193.  
K. Glover, “Multiplicative approximation of linear multivariable systems with Lerror  
bounds,” Proceedings for American Controls Conference, Seattle, 1986, pp. 1705–1709.  
M. Green and BDO Anderson, “The approximation of power spectra by phase matching,”  
Proceedings for 25th CDC, 1986, pp. 1085–1090.  
M. Green and BDO Anderson, “Model reduction by phase matching,” Mathematics of  
Control, Signals, and Systems, Vol. 2, 1989, pp. 221–263.  
© National Instruments Corporation  
A-1  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
     
Appendix A  
[GrA90]  
[Gre88]  
[Gre88a]  
[HiP90]  
Bibliography  
M. Green and BDO Anderson, “Generalized balanced stochastic truncation,”  
Proceedings for 29th CDC, 1990.  
M. Green, “Balanced stochastic realization,” Linear Algebra and Applications, Vol. 98,  
1988, pp. 211–247.  
M. Green, “A relative error bound for balanced stochastic truncation,” IEEE Transactions  
on Automatic Control, Vol. 33, 1988, pp. 961–965.  
D. Hinrichsen and A J Pritchard, “An improved error estimate for reduced-order models  
of discrete-time systems,” IEEE Transactions on Automatic Control, Vol. 35, 1990,  
pp. 317–320.  
[LAL90]  
[Lau80]  
Y. Liu, BDO Anderson, and U-L Ly, “Coprime factorization controller reduction with  
Bezout identity induced frequency weighting,” Automatica, Vol. 26, No. 2, 1990,  
pp. 233–249.  
A. J. Laub, “On computing ‘balancing’ transformations,” Proceedings on Joint American  
Controls Conference, San Francisco, CA, 1980, Section FA8-E.  
[LHPW87] A. J. Laub, M. T. Heath, C. C. Paige, and R. C. Ward, “Computation of system balancing  
transformations and other applications of simultaneous diagonalizing algorithms,” IEEE  
Transactions on Automatic Control, Vol. AC-32, 1987, pp. 115–122.  
[LiA86]  
[LiA89]  
[LiA90]  
[Moo81]  
Y. Liu and BDO Anderson, “Controller reduction via stable factorization and balancing,”  
Int. J. Control, Vol. 44, 1986, pp. 507–531.  
Y. Liu and BDO Anderson, “Singular perturbation approximation of balanced systems,”  
International Journal of Control, Vol. 50, 1989, pp. 1379–1405.  
Y. Liu and BDO Anderson, “Frequency weighted controller reduction methods and loop  
transfer recovery,” Automatica, Vol. 26, No. 3, pp. 487–489.  
B.C. Moore, “Principal component analysis in linear systems: Controllability,  
observability and model reduction,” IEEE Transactions on Automatic Control,  
Vol. AC-26, No. 1, 1981, pp. 17–32.  
[NJB84]  
[PeS82]  
C. N. Nett, C. A. Jacobson, and M. J. Balas, “A connection between state-space  
and doubly coprime fractional representations,” IEEE Transactions on Automatic  
Control, Vol. AC-29, 1984, pp. 831–832.  
L. Pernebo, and L. M. Silverman, “Model reduction via balanced state space  
representations,” IEEE Transactions on Automatic Control, Vol. AC-27, No. 2, 1982,  
pp. 382–387.  
Xmath Model Reduction Module  
A-2  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Appendix A  
Bibliography  
[SaC88]  
[Saf87]  
[SCL90]  
M. G. Safonov and R. Y. Chiang, “Model reduction for robust control: a Schur  
relative-error method,” Proceedings for the American Controls Conference, 1988,  
pp. 1685–1690.  
M. G. Safonov, “Imaginary-axis zeros in multivariable H• optimal control,” Modeling,  
Robustness, and Sensitivity Reduction in Control, (Ed. R. F. Curtain), Springer Verlag,  
Berlin, 1987.  
M. G. Safonov, R. Y. Chiang, and DJN Limebeer, “Optimal Hankel model reduction  
for nonminimal systems,” IEEE Transactions on Automatic Control, Vol. 35 No. 4,  
pp. 496–502, 1990.  
[Vid85]  
M. Vidyasagar, Control Systems Synthesis: A Factorization Approach, MIT Press,  
Cambridge, MA, 1985.  
[WaS90]  
W. Wang and M. G. Safonov, “A tighter relative error bound for balanced stochastic  
truncation,” Systems and Control Letters, Vol. 14, 1990, pp. 307–317.  
[WaS90a] W. Wang and M. G. Safonov, “Comparison between continuous and discrete-time model  
truncation,” Proceedings for the 29th CDC, 1990.  
[BBK88]  
S. Boyd, V. Balakrishnan, and P. Kabamba, “A bisection method for computing the  
Lnorm of a transfer matrix and related problems,” Mathematical Controls, Signals,  
and Systems, Vol. 2, No. 3, pp. 207–219, 1989.  
[BeP79]  
[BoB90]  
A. Berman and R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences.  
Computer Science and Applied Mathematics series, Academic Press, 1979.  
S. Boyd and V. Balakrishnan. “A regularity result for the singular values of a transfer  
matrix and a quadratically convergent algorithm for computing its Lnorm.” Systems  
Control Letters Vol. 15, pp. 1–7, 1990.  
[BoB91]  
[BH69]  
S. Boyd and C. Barratt, Linear Controller Design: Limits of Performance, Prentice-Hall,  
1991.  
A. E. Bryson and Y. C. Ho, Applied Optimal Control, p. 149, Blaisdell Publishing Co.,  
1969.  
[DoS79]  
[DoS81]  
J. C. Doyle and G. Stein. “Robustness with Observers,” IEEE Transactions on Automatic  
Control, August 1979.  
J. C. Doyle and G. Stein. “Multivariable Feedback Design: Concepts for a  
Classical/Modern Synthesis.” IEEE Transactions on Automatic Control, Vol. AC-26,  
February 1981, pp 4–16.  
© National Instruments Corporation  
A-3  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Appendix A  
[Doy82]  
[DWS82]  
Bibliography  
J. C. Doyle. “Analysis of Feedback Systems with Structured Uncertainties.” IEEE  
Proceedings, November 1982.  
J. C. Doyle, J. E. Wall, and G. Stein. “Performance and Robustness Analysis for Structure  
Uncertainties,” Proceedings IEEE Conference on Decision and Control, pp. 629–636,  
1982.  
[FaT88]  
[FaT86]  
[Fr87]  
M. K. Fan and A. L. Tits, “m-form Numerical Range and the Computation of the  
Structured Singular Value.” IEEE Transactions on Automatic Control, Vol. 33,  
pp. 284–289, March 1988.  
M. K. Fan and A. L. Tits. “Characterization and Efficient Computation of the Structured  
Singular Value,” IEEE Transactions on Automatic Control, Vol. AC-31, pp. 734–743,  
August 1986.  
B. Francis, A Course in LControl Theory, Springer-Verlag, Berlin-New York, 1987.  
[FPGM87] D. S. Flamm, S. Boyd, G. Stein, and S. K. Mitter, “Tutorial Workshop on LControl  
Theory,” pre-conference workshop, Proceedings 26th IEEE Conference on Decision and  
Control, December 1988.  
[GD88]  
K. Glover and J. C. Doyle, “State-space formulae for all stabilizing controllers that  
satisfy an Lnorm bound and relations to risk sensitivity,” Systems and Control Letters,  
Vol. 11, pp. 167–172, 1988.  
[DGKF89] J. C. Doyle, K. Glover, P. K. Khargonekar, and B. Francis, “State-space solutions to  
standard H2 and Lcontrol problems,” IEEE Transactions on Automatic Control,  
Vol. AC-34, No. 8, pp. 831–847, August 1989.  
[Gu80]  
N. K. Gupta, “Frequency Shaping of Cost Functionals: An extension of LQG Design  
Methods,” AIAA Journal of Guidance and Control, Vol. 3, No. 6, December 1980.  
[ONR84]  
ONR/Honeywell Workshop on Advances in Multivariable Control, Lecture Notes,  
Minneapolis, MN, 1984.  
[Osb60]  
[Saf82]  
E. E. Osborne, “On Preconditioning of Matrices,” JACM, 7:338–345, 1960.  
M.G. Safonov, “Stability Margins of Diagonally Perturbed Multivariable Feedback  
Systems,” IEEE Proceedings, 129-D:251–256, November 1982.  
[SD83]  
[SD84]  
M. G. Safonov and J. C. Doyle. “Optimal Scaling for Multivariable Stability Margin  
Singular Value Computation,” Proceedings of MECO/EES 1983, Symposium, 1983.  
M. G. Safonov and J. C. Doyle, “Minimizing Conservativeness of Robust Singular  
Values,” Multivariable Control, pp. 197–207, S. G. Tzafestas, ed. D. Reidel Publishing  
Company, 1984.  
Xmath Model Reduction Module  
A-4  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Appendix A  
Bibliography  
[SLH81]  
[SA88]  
[Za81]  
M. G. Safonov, A. J. Laub, and G. L. Hartmann, “Feedback Properties of Multivariable  
Systems: The Role and Use of the Return Difference Matrix,” IEEE Transactions on  
Automatic Control, Vol. AC-26, February 1981.  
G. Stein and M. Athans. “The LQG/LTR Procedure for Multivariable Control Design,”  
IEEE Transactions on Automatic Control, Vol. AC-32, No. 2, February 1987, pp.  
105–114.  
G. Zames, “Feedback and optimal sensitivity: model reference transformations,  
multiplicative semi-norms, and approximate inverses,” IEEE Transactions on Automatic  
Control, Vol. AC-26, pp. 301–320, 1981.  
[KS72]  
H. Kwakernaak and R. Sivan, Linear Optimal Control Systems, Wiley, 1972.  
© National Instruments Corporation  
A-5  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
B
Technical Support and  
Professional Services  
Visit the following sections of the National Instruments Web site at  
ni.comfor technical support and professional services:  
Support—Online technical support resources at ni.com/support  
include the following:  
Self-Help Resources—For immediate answers and solutions,  
visit the award-winning National Instruments Web site for  
software drivers and updates, a searchable KnowledgeBase,  
product manuals, step-by-step troubleshooting wizards, thousands  
of example programs, tutorials, application notes, instrument  
drivers, and so on.  
Free Technical Support—All registered users receive free Basic  
Service, which includes access to hundreds of Application  
Engineers worldwide in the NI Developer Exchange at  
ni.com/exchange. National Instruments Application Engineers  
make sure every question receives an answer.  
Training and Certification—Visit ni.com/trainingfor  
self-paced training, eLearning virtual classrooms, interactive CDs,  
and Certification program information. You also can register for  
instructor-led, hands-on courses at locations around the world.  
System Integration—If you have time constraints, limited in-house  
technical resources, or other project challenges, NI Alliance Program  
members can help. To learn more, call your local NI office or visit  
ni.com/alliance.  
If you searched ni.comand could not find the answers you need, contact  
your local office or NI corporate headquarters. Phone numbers for our  
worldwide offices are listed at the front of this manual. You also can visit  
the Worldwide Offices section of ni.com/niglobalto access the branch  
office Web sites, which provide up-to-date contact information, support  
phone numbers, email addresses, and current events.  
© National Instruments Corporation  
B-1  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
                   
Index  
controller robustness, 4-2  
conventions used in the manual, iv  
Symbols  
*, 1-6  
´, 1-6  
diagnostic tools (NI resources), B-1  
documentation  
A
additive error, reduction, 2-1  
algorithm  
balanced stochastic truncation (bst), 3-4  
fractional representation reduction, 4-18  
Hankel multi-pass, 2-20  
optimal Hankel norm reduction, 2-15  
stable, 5-2  
E
equality bounds, tight, 1-7  
error bound, 2-7  
all-pass transfer function, 1-6, 2-4  
for balanced stochastic truncation, 3-8  
for balanced truncation, 2-2  
for impulse responses, 2-3  
for multiplicative Hankel reduction, 3-16  
for stochastic truncation, 3-9  
error formulas, ophank, 2-22  
error reduction  
B
balance, 1-5, 2-4  
algorithm, 1-11  
balanced realization  
definition, 1-10  
additive, 2-1  
internally balanced, 3-9  
singular perturbation, 2-5  
truncation, 2-2, 2-4  
balanced stochastic truncation, 3-3  
See also bst  
examples (NI resources), B-1  
balmoore, 1-5, 2-4  
algorithm, 1-10  
bst, 1-5, 1-14, 3-3  
for Hankel norm approximation, 2-7  
fracred, 1-5, 4-15  
reduction, 4-18  
frequency weighted error reduction, 1-1, 4-1  
controller reduction, 4-2  
C
compare, 1-5, 5-4  
controller reduction, 4-2  
with fractional representations, 4-5  
© National Instruments Corporation  
I-1  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  
 
Index  
G
N
grammians  
National Instruments support and  
services, B-1  
nomenclature, 1-2  
description of, 1-7  
H
Hankel matrix, 1-9  
ophank, 1-5, 2-14  
Hankel norm approximation, 2-6  
Hankel singular values, 1-8, 3-9, 5-1  
hankelsv, 1-5, 5-1  
discrete-time systems, 2-21  
error formulas, 2-22  
impulse response error, 2-22  
help, technical support, B-1  
onepass, 2-18  
I
instrument drivers (NI resources), B-1  
internal balancing, 1-10  
P
Padé approximation, 1-5  
perturbation  
K
of balanced realization, 2-5  
singular, 1-11, 2-6  
KnowledgeBase, B-1  
phase function, 3-5, 3-15  
phase matrix, 3-6, 3-16  
L
lmax(A), 1-6  
lyapunov, 1-8  
M
reduced order system, 2-3  
Reli(A), 1-6  
MATRIXx Help, 1-3  
minimality requirements, 1-5  
model reduction, schur, 2-5  
Moore-Penrose pseudo-inverse, 1-6  
mreduce, 1-5, 1-13, 2-10  
mulhank, 1-5, 1-14, 3-14  
S
singular perturbation, 1-11  
skipChks, 6-3  
software (NI resources), B-1  
spectral factorization, 1-13  
stability requirements, 1-5  
Xmath Model Reduction Module  
I-2  
ni.com  
Download from Www.Somanuals.com. All Manuals Search And Download.  
Index  
stable, 1-5, 5-2  
sup, 1-6  
unstable zeros, 2-3  
T
Web resources, B-1  
wtbalance, 1-5, 4-10  
technical support, B-1  
tight equality bounds, 1-7  
training and certification (NI resources), B-1  
transfer function, allpass, 1-6  
troubleshooting (NI resources), B-1  
truncate, 1-5, 2-4, 2-11  
© National Instruments Corporation  
I-3  
Xmath Model Reduction Module  
Download from Www.Somanuals.com. All Manuals Search And Download.  

Miele Ventilation Hood 09 968 470 User Manual
Nilfisk ALTO Pressure Washer Multi User Pressure Washer System User Manual
Nilfisk ALTO Vacuum Cleaner UZ 934 User Manual
Nokia Cell Phone 2275 User Manual
NuTone Fan PFMB42 User Manual
Olympus Binoculars 10 42 EXWP I User Manual
Onkyo Speaker System HTP 103E User Manual
Panasonic Computer Drive CF VFS721 User Manual
Parasound Portable Speaker PAL 80 User Manual
Parker Hannifin Video Game Console P N88023735 01A User Manual