FPE7e Appendices

Appendix WA
A Review of Complex
Variables
This appendix is a brief summary of some results on complex variables
theory, with emphasis on the facts needed in control theory. For a comprehensive study of basic complex variables theory, see standard textbooks such
as Brown and Churchill (1996) or Marsden and Hoffman (1998).
WA.1 Denition of a Complex Number

The complex numbers are distinguished from purely real numbers in that
they also contain the imaginary operator, which we shall denote j. By
definition,
j2 = 1 or j = 1.
(WA.1)
A complex number may be defined as
A = + j,
(WA.2)
where is the real part and is the imaginary part, denoted, respectively, as
= Re(A),
= Im(A).
(WA.3)
Note that the imaginary part of A is itself a real number.

Graphically, we may represent the complex number A in two ways. In
the Cartesian coordinate system (Fig. WA.1a), A is represented by a single
point in the complex plane. In the polar coordinate system, A is represented
by a vector with length r and an angle ; the angle is measured in radians
counterclockwise from the positive real axis (Fig. WA.1b). In polar form the
complex number A is denoted by
A = |A| arg A = r = re j ,
0 2 ,
(WA.4)
where rcalled the magnitude, modulus, or absolute value of Ais the

length of the vector representing A, namely,

r = |A| = 2 + 2 ,
(WA.5)
and where is given by
tan =
or
= arg(A) = tan1
(WA.6)
(WA.7)
Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 1 #1
2 Appendix WA A Review of Complex Variables

Im (s)
Figure WA.1
The complex number A
represented in
(a) Cartesian and
(b) polar coordinates
Im (s)
A
jx
jx
r
u
Re (s)
(a)
Im (s)
Re (s)
(b)
Im (s)
Im (s)
A1 + A2
A1 A2
A2
A2
A2
A1 u2
A1
u2
u1 + u2
A1
u1
u1
Re (s)
Re (s)
u1 - u2
Re (s)
A1
A2
(a)
(b)
(c)
Figure WA.2
Arithmetic of complex numbers: (a) addition; (b) multiplication; (c) division
Care must be taken to compute the correct value of the angle, depending
on the sign of the real and imaginary parts (that is, one must find the quadrant
in which the complex number lies).
The conjugate of A is defined as
A = j
(WA.8)
Therefore,
(A ) = A,
A1
(WA.9)
A2 ,
(A1 A2 ) =

A
A1
= 1 ,
A2
A2
(WA.10)
(WA.11)
(A1 A2 ) = A1 A2 ,
Re(A) =
A + A
,
2
(WA.12)
Im(A) =
A A
,
2j
AA = (|A|)2 .
Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 2 #2
(WA.13)
(WA.14)
WA.2 Algebraic Manipulations
WA.2 Algebraic Manipulations

WA.2.1 Complex Addition
If we let
A1 = 1 + j1
and
A2 = 2 + j2 ,
(WA.15)
then
A1 + A2 = (1 + j1 ) + (2 + j2 ) = (1 + 2 ) + j(1 + 2 ). (WA.16)
Because each complex number is represented by a vector extending from
the origin, we can add or subtract complex numbers graphically. The sum is
obtained by adding the two vectors. This we do by constructing a parallelogram and finding its diagonal, as shown in Fig. WA.2a. Alternatively, we
could start at the tail of one vector, draw a vector parallel to the other vector,
and then connect the origin to the new arrowhead.
Complex subtraction is very similar to complex addition.
WA.2.2 Complex Multiplication

For two complex numbers defined according to Eq. (WA.15),
A1 A2 = (1 + j1 )(2 + j2 )
= (1 2 1 2 ) + j(1 2 + 1 2 ).
(WA.17)
The product of two complex numbers may be obtained graphically using

polar representations, as shown in Fig. WA.2b.
WA.2.3 Complex Division

The division of two complex numbers is carried out by rationalization. This
means that both the numerator and denominator in the ratio are multiplied
by the conjugate of the denominator:
A1 A2
A1
=
A2
A2 A2
=
(1 2 + 1 2 ) + j(1 2 1 2 )
.
22 + 22
(WA.18)
From Eq. (WA.4) it follows that

A1 =
1 j
e ,
r
r = 0.
(WA.19)
Also, if A1 = r1 e j1 and A2 = r2 e j2 , then

A1 A2 = r1 r2 e j(1 +2 ) ,
(WA.20)
where |A1 A2 | = r1 r2 and arg(A1 A2 ) = 1 + 2 , and

r1
A1
= e j(1 2 ) ,
A2
r2
r2 = 0,
Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 3 #3
(WA.21)

A1 r1
A
where A
= r2 and arg A1 = 1 2 . The division of complex numbers
2
2
may be carried out graphically in polar coordinates as shown in Fig. WA.2c.
EXAMPLE WA.1
Frequency Response of First-Order System

1 , where
Find the magnitude and phase of the transfer function G(s) = s +
1
s = j.
Solution. Substituting s = j and rationalizing, we obtain
1
+ 1 j
+ 1 + j + 1 j
+ 1 j
=
.
( + 1)2 + 2
G( j) =
Therefore, the magnitude and phase are

1
( + 1)2 + 2
|G( j)| =
=
,
( + 1)2 + 2
( + 1)2 + 2

Im(G( j))
arg(G( j)) = tan1

= tan1
Re(G( j))
+1
WA.3 Graphical Evaluation of Magnitude and

Phase
Consider the transfer function
m
(s + zi )
G(s) = ni=1
.
i=1 (s + pi )
(WA.22)
The value of the transfer function for sinusoidal inputs is found by replacing
s with j. The gain and phase are given by G( j) and may be determined
analytically or by a graphical procedure. Consider the pole-zero configuration for such a G(s) and a point s0 = j0 on the imaginary axis, as shown
in Fig. WA.3. Also consider the vectors drawn from the poles and the zero
to s0 . The magnitude of the transfer function evaluated at s0 = j0 is simply
the ratio of the distance from the zero to the product of all the distances from
the poles:
r1
|G( j0 )| =
.
(WA.23)
r2 r3 r4
The phase is given by the sum of the angles from the zero minus the
sum of the angles from the poles:
arg G( j0 ) = G( j0 ) = 1 (2 + 3 + 4 ).
(WA.24)
This may be explained as follows. The term s + z1 is a vector addition of its

two components. We may determine this equivalently as s (z1 ), which
amounts to translation of the vector s + z1 starting at z1 , as shown in
Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 4 #4
WA.4 Differentiation and Integration
Im (s)
Figure WA.3
Graphical
determination of
magnitude and phase
s0
r4
r2
r3
u2
r1
u1
u4
Re (s)
u3
Im (s)
Figure WA.4
Illustration of graphical
computation of s + z1
s0 + z1
- z1
s0
s0 + z1
z1
Re (s)
Fig. WA.4. This means that a vector drawn from the zero location to s0 is
equivalent to s+z1 . The same reasoning applies to the poles. We reflect p1 , p2 ,
and p3 about the origin to obtain the pole locations. Then the vectors drawn
from p1 , p2 , and p3 to s0 are the same as the vectors in the denominator
represented in polar coordinates. Note that this method may also be used to
evaluate s0 at places in the complex plane besides the imaginary axis.
WA.4 Differentiation and Integration

The usual rules apply to complex differentiation. Let G(s) be differentiable
with respect to s. Then the derivative at s0 is defined as
G(s) G(s0 )
,
(WA.25)
s s0
provided that the limit exists. For conditions on the existence of the
derivative, see Brown and Churchill (1996).
The standard rules also apply to integration, except that the constant of
integration c is a complex constant:

G(s)ds = Re[G(s)]ds + j Im[G(s)]ds + c.
(WA.26)
G (s0 ) = lim
ss0
Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 5 #5
WA.5 Eulers Relations

Let us now derive an important relationship involving the complex exponential. If we define
A = cos + j sin ,
(WA.27)
where is in radians, then
dA
= sin + j cos = j2 sin + j cos
d
= j(cos + j sin ) = jA.
(WA.28)
We collect the terms involving A to obtain

dA
= jd .
A
Integrating both sides of Eq. (WA.29) yields
ln A = j + c,
(WA.29)
(WA.30)
where c is a constant of integration. If we let = 0 in Eq. (WA.30), we find

that c = 0 or
A = e j = cos + j sin .
(WA.31)
Similarly,
Eulers relations
A = ej = cos j sin .
(WA.32)
From Eqs. (WA.31) and (WA.32) it follows that

e j + ej
,
2
e j ej
sin =
.
2j
cos =
(WA.33)
(WA.34)
WA.6 Analytic Functions

Let us assume that G is a complex-valued function defined in the complex
plane. Let s0 be in the domain of G, which is assumed to be finite within
some disk centered at s0 . Thus, G(s) is defined not only at s0 but also at all
points in the disk centered at s0 . The function G is said to be analytic if its
derivative exists at s0 and at each point in the neighborhood of s0 .
WA.7 Cauchys Theorem

A contour is a piecewise-smooth arc that consists of a number of smooth arcs
joined together. A simple closed contour is a contour that does not intersect
itself and ends on itself. Let C be a closed contour as shown in Fig. WA.5a,
and let G be analytic inside and on C. Cauchys theorem states that
G(s)ds = 0.
(WA.35)
C
Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 6 #6
WA.8 Singularities and Residues

Im (s)
Figure WA.5
Contours in the s-plane:
(a) a closed contour;
(b) two different paths
between A1 and A2
Im (s)
C1
A2
A1
A2
C2
A1
Re (s)
Re (s)
(a)
(b)
There is a corollary to this theorem: Let C1 and C2 be two paths connecting

the points A1 and A2 as in Fig. WA.5b. Then,

G(s)ds =
G(s)ds.
(WA.36)
C1
C2
WA.8 Singularities and Residues

If a function G(s) is not analytic at s0 but is analytic at some point in every
neighborhood of s0 , it is said to be a singularity. A singular point is said to
be an isolated singularity if G(s) is analytic everywhere else in the neighborhood of s0 except at s0 . Let G(s) be a rational function (that is, a ratio
of polynomials). If the numerator and denominator are both analytic, then
G(s) will be analytic except at the locations of the poles (that is, at roots
of the denominator). All singularities of rational algebraic functions are the
pole locations.
Let G(s) be analytic except at s0 . Then we may write G(s) in its Laurent
series expansion form:
G(s) =
An
A1
+ ... +
+ B0 + B1 (s s0 ) + . . . . (WA.37)
(s s0 )n
(s s0 )
The coefficient A1 is called the residue of G(s) at s0 and may be

evaluated as
1
A1 = Res[G(s); s0 ] =
G(s) ds,
(WA.38)
2 j C
where C denotes a closed arc within an analytic region centered at s0 that
contains no other singularity, as shown in Fig. WA.6. When s0 is not repeated
with n = 1, we have
A1 = Res[G(s); s0 ] = (s s0 )G(s)|s=s0 .
This is the familiar cover-up method of computing residues.
Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 7 #7
(WA.39)

Im (s)
Figure WA.6
Contour around an
isolated singularity
s0
Re (s)
WA.9 Residue Theorem

If the contour C contains l singularities, then Eq. (WA.39) may be generalized
to yield Cauchys residue theorem:
1
2 j
G(s) ds =
Res[G(s); si ].
(WA.40)
i=1
WA.10 The Argument Principle

Before stating the argument principle, we need a preliminary result from
which the principle follows readily.
Number of Poles and Zeros

Let G(s) be an analytic function inside and on a closed contour C, except
for a finite number of poles inside C. Then, for C described in the positive
sense (clockwise direction),

1
G (s)
ds = N P
(WA.41)
2 j
G(s)
or
1
2 j
d(ln G) = N P,
(WA.42)
where N and P are the total number of zeros and poles of G inside C,
respectively. A pole or zero of multiplicity k is counted k times.
Proof Let s0 be a zero of G with multiplicity k. Then, in some
neighborhood of that point, we may write G(s) as
G(s) = (s s0 )k f (s),
Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 8 #8
(WA.43)
WA.10 The Argument Principle
where f (s) is analytic and f (s0 ) = 0. If we differentiate Eq. (WA.43), we

obtain
G (s) = k(s s0 )k1 f (s) + (s s0 )k f (s).
(WA.44)
Equation (WA.44) may be rewritten as
k
f (s)
G (s)
=
.
+
G(s)
s s0
f (s)
(WA.45)
Therefore, G (s)/G(s) has a pole at s = s0 with residue K. This analysis may

be repeated for every zero. Hence, the sum of the residues of G (s)/G(s) is
the number of zeros of G(s) inside C. If s0 is a pole with multiplicity l, we
may write
h(s) = (s s0 )l G(s),
(WA.46)
where h(s) is analytic and h(s0 ) = 0. Then Eq. (WA.46) may be rewritten
as
h(s)
G(s) =
.
(WA.47)
(s s0 )l
Differentiating Eq. (WA.47) we obtain
G (s) =
h (s)
lh(s)
,
l
(s s0 )
(s s0 )l+1
(WA.48)
so that
G (s)
l
h (s)
=
.
(WA.49)
+
G(s)
s s0
h(s)
This analysis may be repeated for every pole. The result is that the sum of
the residues of G (s)/G(s) at all the poles of G(s) is P.
The Argument Principle

Using Eq. (WA.38), we get
1
2 j
d[ln G(s)] = N P,
(WA.50)
where d[ln G(s)] was substituted for G (s)/G(s). If we write G(s) in polar
form, then
d[ln G(s)] =
d{ln |G(s)| + j arg[ln G(s)]}

s=s2
2
= ln |G(s)||s=s
s=s1 + j arg G(s)|s=s1 .
(WA.51)
Because is a closed contour, the first term is zero, but the second term is
2 times the net encirclements of the origin:
1
d[ln G(s)] = N P.
(WA.52)
2 j
Intuitively, the argument principle may be stated as follows: We let G(s)
be a rational function that is analytic except at possibly a finite number of
points. We select an arbitrary contour in the s-plane so that G(s) is analytic
at every point on the contour (the contour does not pass through any of the
singularities). The corresponding mapping into the G(s)-plane may encircle
Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 9 #9

the origin. The number of times it does so is determined by the difference
between the number of zeros and the number of poles of G(s) encircled
by the s-plane contour. The direction of this encirclement is determined by
which is greater, N (clockwise) or P (counterclockwise). For example, if the
contour encircles a single zero, the mapping will encircle the origin once in
the clockwise direction. Similarly, if the contour encloses only a single
pole, the mapping will encircle the origin, this time in the counterclockwise
direction. If the contour encircles no singularities, or if the contour encloses
an equal number of poles and zeros, there will be no encirclement of the
origin. A contour evaluation of G(s) will encircle the origin if there is a
nonzero net difference between the encircled singularities. The mapping
is conformal as well, which means that the magnitude and sense of the
angles between smooth arcs is preserved. Chapter 6 provides a more detailed
intuitive treatment of the argument principle and its application to feedback
control in the form of the Nyquist stability theorem.
WA.11 Bilinear Transformation

A bilinear transformation is of the form
as + b
w=
,
(WA.53)
cs + d
where a, b, c, d are complex constants, and it is assumed that ad bc = 0.
The bilinear transformation always transforms circles in the w-plane into
circles in the s-plane. This can be shown in several ways. If we solve for s,
we obtain
dw + b
s=
.
(WA.54)
cw a
The equation for a circle in the w-plane is of the form
|w |
= R.
|w |
(WA.55)
If we substitute for w in terms of s in Eq. (WA.53), we get

|s |
= R ,
|s |
where
=
d b
,
a c
=
d b
,
a c

a c
R,
R =
a c
(WA.56)
(WA.57)
which is the equation for a circle in the s-plane. For alternative proofs
the reader is referred to Brown and Churchill (1996) and Marsden and
Hoffman (1998).
Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 10 #10
Appendix WB
Summary of Matrix Theory
In the text, we assume you are already somewhat familiar with matrix theory and with the solution of linear systems of equations. However, for the
purposes of review we present here a brief summary of matrix theory with
an emphasis on the results needed in control theory. For further study, see
Strang (2006) and Gantmacher (1959).
WB.1 Matrix Denitions

An array of numbers arranged in rows and columns is referred to as a matrix.
If A is a matrix with m rows and n columns, an m n (read m by n) matrix,
it is denoted by
a11 a12 a1n

a21 a22 a2n
(WB.1)
A= .
..
.. ,
.
.
.
.
am1 am2 amn
where the entries aij are its elements. If m = n, then the matrix is square;
otherwise it is rectangular. Sometimes a matrix is simply denoted by A =
[aij ]. If m = 1 or n = 1, then the matrix reduces to a row vector or a column
vector, respectively. A submatrix of A is the matrix with certain rows and
columns removed.
WB.2 Elementary Operations on Matrices

If A and B are matrices of the same dimension, then their sum is defined by
C = A + B,
(WB.2)
cij = aij + bij .
(WB.3)
where
Commutative law for

addition
Associative law for
addition
That is, the addition is done element by element. It is easy to verify the
following properties of matrices:
A + B = B + A,
(WB.4)
(A + B) + C = A + (B + C).
(WB.5)
Two matrices can be multiplied if they are compatible. Let A = m n and

B = n p. Then the m p matrix
C = AB
(WB.6)
11
Z02_FRAN6598_07_SE_WB 2014/5/6 0:26 page 11 #1
12 Appendix WB Summary of Matrix Theory

is the product of the two matrices, where
cij =
n

aik bkj .
(WB.7)
k=1
Associative law for

multiplication
Matrix multiplication satisfies the associative law

A(BC) = (AB)C,
(WB.8)
but not the commutative law; that is, in general,

AB = BA.
(WB.9)
WB.3 Trace
The trace of a square matrix is the sum of its diagonal elements:
trace A =
n

aii .
(WB.10)
i=1
WB.4 Transpose
The n m matrix obtained by interchanging the rows and columns of A is
called the transpose of matrix A:
a11 a21 . . . am1

a12 a22 . . . am2
AT = .
..
..
..
.
.
a1n a2n . . . amn
A matrix is said to be symmetric if
AT = A.
Transposition
(WB.11)
It is easy to show that

(AB)T = BT AT ,
(WB.12)
(ABC) = C B A ,
(WB.13)
(A + B) = A + B .
(WB.14)
T
T
WB.5 Determinant and Matrix Inverse

The determinant of a square matrix is defined by Laplaces expansion
det A =
n

aij ij
for any
i = 1, 2, . . . , n,
(WB.15)
j=1
where ij is called the cofactor and

ij = (1)i+j det Mij ,
Z02_FRAN6598_07_SE_WB 2014/5/6 0:26 page 12 #2
(WB.16)
WB.6 Properties of the Determinant
13
where the scalar det Mij is called a minor. Mij is the same as the matrix A
except that its ith row and jth column have been removed. Note that Mij is
always an (n 1) (n 1) matrix, and that the minors and cofactors are
identical except possibly for a sign.
The adjugate of a matrix is the transpose of the matrix of its cofactors:
adj A = [ij ]T .
(WB.17)
A adj A = (det A)I,
(WB.18)
It can be shown that

Identity matrix
where I is called the identity matrix:
1 0 ... ... 0
0 1
0 ... 0
I= .
.
.. ;
..
..
..
.
.
0 ... ... 0 1
that is, I has ones along the diagonal and zeros elsewhere. If det A = 0, then
the inverse of a matrix A is defined by
A1 =
and has the property
adj A
,
det A
AA1 = A1 A = I.
(WB.19)
(WB.20)
Note that a matrix has an inversethat is, it is nonsingularif its

determinant is nonzero.
The inverse of the product of two matrices is the product of the inverse
of the matrices in reverse order:
Inversion
and
(AB)1 = B1 A1
(WB.21)
(ABC)1 = C1 B1 A1 .
(WB.22)
WB.6 Properties of the Determinant

When dealing with determinants of matrices, the following elementary (row
or column) operations are useful:
1. If any row (or column) of A is multiplied by a scalar , the resulting
has the determinant
matrix A
= det A.
det A
(WB.23)
Hence
det(A) = n det A.
(WB.24)
then
2. If any two rows (or columns) of A are interchanged to obtain A,
= det A.
det A
Z02_FRAN6598_07_SE_WB 2014/5/6 0:26 page 13 #3
(WB.25)
3. If a multiple of a row (or column) of A is added to another to obtain A,

then
= det A.
det A
(WB.26)
4. It is also easy to show that
det A = det AT
(WB.27)
det AB = det A det B.
(WB.28)
and
Applying Eq. (WB.28) to Eq. (WB.20), we have
det A det A1 = 1.
(WB.29)
If A and B are square matrices, then the determinant of the block

triangular matrix is the product of the determinants of the diagonal
blocks:

A C
det
= det A det B.
(WB.30)
0 B
If A is nonsingular, then

A B
det
= det A det(D CA1 B).
C D
(WB.31)
Using this identity, we can write the transfer function of a scalar system
in a compact form:

(sI A) B
det
C
D
.
(WB.32)
G(s) = C(sI A)1 B + J =
det(sI A)
WB.7 Inverse of Block Triangular Matrices

If A and B are square invertible matrices, then

1 1
A C
A
A1 CB1
.
=
0 B
0
B1
(WB.33)
WB.8 Special Matrices

Diagonal matrix
Some matrices have special structures and are given names. We have already
defined the identity matrix, which has a special form. A diagonal matrix has
(possibly) nonzero elements along the main diagonal and zeros elsewhere:
0
a11
a22
a
33
(WB.34)
A=
.
.
..
0
ann
Z02_FRAN6598_07_SE_WB 2014/5/6 0:26 page 14 #4
WB.10 Characteristic Polynomial

Upper triangular matrix
15
A matrix is said to be (upper) triangular if all the elements below the main
diagonal are zeros:
a11 a12
a1n
0 a22
..
..
0
.
A= .
(WB.35)
.
..
.
.
..
..
0
.
0
0 0 ann
The determinant of a diagonal or triangular matrix is simply the product of
its diagonal elements.
A matrix is said to be in the (upper) companion form if it has the
structure
an
a1 a2
1
0
0
1
0
0
(WB.36)
Ac =
.
..
..
..
.
.
.
0
1
0
Note that all the information is contained in the first row. Variants of this form
are the lower, left, or right companion matrices. A Vandermonde matrix
has the following structure:
1 a1 a12 a1n1
1 a2 a2 an1
2
2
A= . .
.
(WB.37)
..
..
.. ..
.
.
1 an an2 ann1
WB.9 Rank
The rank of a matrix is the number of its linearly independent rows or
columns. If the rank of A is r, then all (r + 1) (r + 1) submatrices of A
are singular, and there is at least one r r submatrix that is nonsingular. It
is also true that
row rank of A = column rank of A.
(WB.38)
WB.10 Characteristic Polynomial

The characteristic polynomial of a matrix A is defined by
a(s) det(sI A)
= sn + a1 sn1 + + an1 s + an ,
(WB.39)
where the roots of the polynomial are referred to as eigenvalues of A. We

can write
a(s) = (s 1 )(s 2 ) (s n ),
(WB.40)
Z02_FRAN6598_07_SE_WB 2014/5/6 0:26 page 15 #5

where {i } are the eigenvalues of A. The characteristic polynomial of a
companion matrix (for example, Eq. (WB.36)) is
a(s) = det(sI Ac )
= sn + a1 sn1 + + an1 s + an .
(WB.41)
WB.11 CayleyHamilton Theorem

The CayleyHamilton theorem states that every square matrix A satisfies
its characteristic polynomial. This means that if A is an n n matrix with
characteristic equation a(s), then
a(A) An + a1 An1 + + an1 A + an I = 0
(WB.42)
WB.12 Eigenvalues and Eigenvectors

Any scalar and nonzero vector v that satisfy
Av = v,
(WB.43)
are referred to as the eigenvalue and the associated (right) eigenvector of

the matrix A [because v appears to the right of A in Eq. (WB.43)]. By
rearranging terms in Eq. (WB.43), we get
(I A)v = 0.
(WB.44)
det(I A) = 0,
(WB.45)
Because v is nonzero,
so is an eigenvalue of the matrix A as defined in Eq. (WB.43). The normalization of the eigenvectors is arbitrary; that is, if v is an eigenvector, so
is v. The eigenvectors are usually normalized to have unit length; that is,
v2 = vT v = 1.
If wT is a nonzero row vector such that
wT A = wT ,
(WB.46)
then w is called a left eigenvector of A [because wT appears to the left of

A in Eq. (WB.46)]. Note that we can write
AT w = w
(WB.47)
so that w is simply a right eigenvector of AT .
WB.13 Similarity Transformations

Consider the arbitrary nonsingular matrix T such that
= T1 AT.
A
(WB.48)
The matrix operation shown in Eq. (WB.48) is referred to as a similarity

transformation. If A has a full set of eigenvectors, then we can choose T
will be diagonal.
to be the set of eigenvectors and A
Z02_FRAN6598_07_SE_WB 2014/5/6 0:26 page 16 #6
WB.14 Matrix Exponential
17
Consider the set of equations in state-variable form:

x = Ax + Bu.
(WB.49)
T = x,
(WB.50)
T = AT + Bu,
(WB.51)
If we let
then Eq. (WB.49) becomes
and premultiplying both sides by T1 , we get

= T1 AT + T1 Bu
+ Bu,
= A
(WB.52)
= T1 AT,
A
B = T1 B.
(WB.53)
where
is
The characteristic polynomial of A
= det(sI T1 AT)
det(sI A)
= det(sT1 T T1 AT)
= det[T1 (sI A)T]
= det T1 det(sI A) det T.
(WB.54)
Using Eq. (WB.29), Eq. (WB.54) becomes

= det(sI A).
det(sI A)
(WB.55)
and A both have the same characteristic

From Eq. (WB.55) we can see that A
polynomial, giving us the important result that a similarity transformation
does not change the eigenvalues of a matrix. From Eq. (WB.50) a new state
made up of a linear combination from the old state has the same eigenvalues
as the old set.
WB.14 Matrix Exponential

Let A be a square matrix. The matrix exponential of A is defined as the
series
1
A3 t 3
eAt = I + At + A2 t 2 +
+
(WB.56)
2!
3!
It can be shown that the series converges. If A is an n n matrix, then eAt
is also an n n matrix and can be differentiated:
d At
e = AeAt .
dt
Z02_FRAN6598_07_SE_WB 2014/5/6 0:26 page 17 #7
(WB.57)

Other properties of the matrix exponential are
eAt1 eAt2 = eA(t1 +t2 )
(WB.58)
eA eB = eB eA .
(WB.59)
and, in general,
(In the exceptional case where A and B commutethat is, AB = BAthen
eA eB = eB eA .)
WB.15 Fundamental Subspaces

The range space of A, denoted by R(A) and also called the column space
of A, is defined by the set of vectors
x = Ay
(WB.60)
for some vector y. The null space of A, denoted by N (A), is defined by the
set of vectors x such that
Ax = 0.
(WB.61)
If x N (A) and y R(AT ), then yT x = 0; that is, every vector in the null
space of A is orthogonal to every vector in the range space of AT .
WB.16 Singular-Value Decomposition

The singular-value decomposition (SVD) is one of the most useful tools
in linear algebra and has been widely used in control theory during the last
few decades. Let A be an m n matrix. Then there always exist matrices U,
S, and V such that
A = USVT .
(WB.62)
Here U and V are orthogonal matrices; that is,
UUT = I, VVT = I.
(WB.63)
S is a quasidiagonal matrix with singular values as its diagonal elements;

that is,

0
S=
,
(WB.64)
0 0
where is a diagonal matrix of nonzero singular values in descending order:
1 2 r > 0.
(WB.65)
The unique diagonal elements of S are called the singular values. The
maximum singular value is denoted by (A), and the minimum singular
value is denoted by (A). The rank of the matrix is the same as the number
of nonzero singular values. The columns of U and V,
U = [ u1
u2
. . . um ],
V = [ 1
...
n ],
Z02_FRAN6598_07_SE_WB 2014/5/6 0:26 page 18 #8
(WB.66)
WB.18 Matrix Identity
19
are called the left and right singular vectors, respectively. SVD provides
complete information about the fundamental subspaces associated with a
matrix:
N (A) = span[ r+1
r+2
...
n ],
R(A) = span[ u1
u2
. . . ur ],
R(A ) = span[ 1
...
N (AT ) = span[ ur+1
ur+2
r ],
. . . um ].
(WB.67)
Here N denotes the null space and R, the range space, respectively.
The norm of the matrix A, denoted by A2 , is given by
A2 = (A).
(WB.68)
If A is a function of , then the infinity norm of A, A , is given by

A(j) = max (A).
(WB.69)
WB.17 Positive Denite Matrices

A matrix A is said to be positive semidefinite if
xT Ax 0
for all x.
(WB.70)
The matrix is said to be positive definite if equality holds in Eq. (WB.70)

only for x = 0. A symmetric matrix is positive definite if and only if all of
its eigenvalues are positive. It is positive semidefinite if and only if all of its
eigenvalues are nonnegative.
An alternate method for determining positive definiteness is to test the
minors of the matrix. A matrix is positive definite if all the leading principal
minors are positive, and positive semidefinite if they are all nonnegative.
WB.18 Matrix Identity

If A is an n m matrix and B is an m n matrix, then
det[In AB] = det[Im BA],
where In and Im are identity matrices of size n and m, respectively.
Z02_FRAN6598_07_SE_WB 2014/5/6 0:26 page 19 #9
Appendix WC
Controllability and
Observability
Controllability and observability are important structural properties of
dynamic systems. First identified and studied by Kalman (1960) and later
by Kalman et al. (1961), these properties have continued to be examined
during the last five decades. We will discuss only a few of the known results
for linear constant systems with one input and one output. In the text, we
discuss these concepts in connection with control law and estimator designs.
For example, in Section 7.4 we suggest that, if the square matrix given by
C = [B AB A2 B . . . An1 B]
(WC.1)
is nonsingular, by transformation of the state we can convert the given

description into control canonical form. We can then construct a control law
that will give the closed-loop system an arbitrary characteristic equation.
WC.1 Controllability
We begin our formal discussion of controllability with the first of four
definitions.
Definition WC.1 The system (A, B) is controllable if, for any given nthorder polynomial c (s), there exists a (unique) control law u = Kx such
that the characteristic polynomial of A BK is c (s).
From the results of Ackermanns formula (see Appendix WD), we have
the following mathematical test for controllability: (A, B) is a controllable
pair if and only if the rank of C is n. Definition WC.1 based on pole placement
is a frequency-domain concept. Controllability can be equivalently defined
in the time domain.
Definition WC.2 The system (A, B) is controllable if there exists a (piecewise continuous) control signal u(t) that will take the state of the system
from any initial state x0 to any desired final state xf in a finite time interval.
We will now show that the system is controllable by this definition if and
only if C is full rank. We first assume that the system is controllable but that
rank[B AB A2 B . . . An1 B] < n.
(WC.2)
We can then find a vector v such that

v[B AB A2 B . . . An1 B] = 0
20
Z03_FRAN6598_07_SE_WC 2014/5/8 0:19 page 20 #1
(WC.3)
21
or
vB = vAB = vA2 B = . . . = vAn1 B = 0.
(WC.4)
The Cayley-Hamilton theorem states that A satisfies its own characteristic

equation, namely,
An = a1 An1 + a2 An2 + . . . + an I.
(WC.5)
Therefore,
vAn B = a1 vAn1 B + a2 vAn2 B + . . . + an vB = 0.
(WC.6)
By induction, vA B = 0 for k = 0, 1, 2, . . ., or vA B = 0 for

m = 0, 1, 2, . . ., and thus

1
veAt B = v I + At + A2 t 2 + . . . B = 0,
(WC.7)
2!
n+k
for all t. However, the zero initial-condition response (x0 = 0) is

t
x(t) =
veA(t ) Bu( ) d
0
= eAt
eA Bu( ) d .
(WC.8)
Using Eq. (WC.7), Eq. (WC.8) becomes

t
vx(t) =
veA(t ) Bu( ) d = 0
(WC.9)
for all u(t) and t > 0. This implies that all points reachable from the origin are
orthogonal to v. This restricts the reachable space and therefore contradicts
the second definition of controllability. Thus if C is singular, (A, B) is not
controllable by Definition WC.2.
Next we assume that C is full rank but (A, B) is uncontrollable by
Definition WC.2. This means that there exists a nonzero vector v such that
tf
(tf )
v
eA
Bu( ) d = 0,
(WC.10)
0
because the whole state-space is not reachable. But Eq. (WC.10) implies
that
(tf )
veA
B = 0,
0 tf .
(WC.11)
If we set = tf , we see that vB = 0. Also, differentiating Eq. (WC.11) and

letting = tf gives vAB = 0. Continuing this process, we find that
vB = vAB = vA2 B = . . . = vAn1 B = 0,
(WC.12)
which contradicts the assumption that C is full rank.

We have now shown that the system is controllable by Definition WC.2
if and only if the rank of C is n, which is exactly the same condition we
found for pole assignment.
Z03_FRAN6598_07_SE_WC 2014/5/8 0:19 page 21 #2
22 Appendix WC Controllability and Observability

Figure WC.1
Block diagram of a
system with a diagonal
matrix
b1
1
s + n1
c1
b2
1
s + n2
c2
1
s + nn
cn
bn
y
+
Our final definition comes closest to the structural character of controllability.

Definition WC.3 The system (A, B) is controllable if every mode of A is
connected to the control input.
Because of the generality of the modal structure of systems, we will treat
only the case of systems for which A can be transformed to diagonal form.
(The double-integration plant does not qualify.) Suppose we have a diagonal
matrix Ad and its corresponding input matrix Bd with elements bi . The
structure of such a system is shown in Fig. (WC.1). By definition for a
controllable system, the input must be connected to each mode so that the bi
are all nonzero. However, this is not enough if the poles (i ) are not distinct.
Suppose, for instance, that 1 = 2 . The first two state equations are then
x 1d = 1 x1d + b1 u,
x 2d = 1 x2d + b2 u.
(WC.13)
If we define a new state, = b2 x1d b1 x2d , the equation for is

= b2 x 1d b1 x 2d = b2 1 x1d +b2 b1 ub1 1 x2d b1 b2 u = 1 , (WC.14)
which does not include the control u; hence, is not controllable. The point
is that if any two poles are equal in a diagonal Ad system with only one input,
we effectively have a hidden mode that is not connected to the control, and
the system is not controllable (Fig. WC.2a). This is because the two state
variables move together exactly, so we cannot independently control x1d
1
s+1
u
1
s+1
1
s+1
y
u
s-1
s+1
1
s-1
1
s-1
(a)
(b)
(c)
Figure WC.2
Examples of uncontrollable systems
Z03_FRAN6598_07_SE_WC 2014/5/8 0:19 page 22 #3
23
and x2d . Therefore, even in such a simple case, we have two conditions for
controllability:
1. All eigenvalues of Ad are distinct.
2. No element of Bd is zero.
Now let us consider the controllability matrix of this diagonal system.
By direct computation,
b1 b1 1 . . . b1 n1
1
..
b2 b2 2 . . .
C= .
.
.
..
.
..
.
n1
bn bn n . . . bn n
1 2 . . . n1
1
b1
0
1
1
..
2 ...
b
2
1
2
2
.
=
(WC.15)
. .
.
.
.
.
. .
.
.
.
. .
.
.
0
bn
1 n 2n . . . n1
n
Note that the controllability matrix C is the product of two matrices and
is nonsingular if and only if both of these matrices are invertible. The first
matrix has a determinant that is the product of bi , and the second matrix
(called a Vandermonde matrix) is nonsingular if and only if the i are distinct.
Thus Definition WC.3 is equivalent to having a nonsingular C also.
Important to the subject of controllability is the PopovHautus
Rosenbrock (PHR) test (see Rosenbrock, 1970, and Kailath, 1980), which
is an alternate way to test the rank (or determinant) of C. The system (A, B)
is controllable if the system of equations
vT [sI A
B] = 0T ,
(WC.16)
has only the trivial solution vT = 0T that is, if the matrix pencil
[sI A
B] = n,
(WC.17)
is full rank for all s, or if there is no nonzero vT such that1

vT A = svT ,
(WC.18)
v B = 0.
(WC.19)
This test is equivalent to the rank-of-C test. It is easy to show that, if such
a vector v exists, C is singular. For if a nonzero v exists such that vT B = 0
then by Eqs. (WC.18) and (WC.19), we have
vT AB = svT B = 0.
(WC.20)
Then, multiplying by AB, we find that

vT A2 B = svT AB = 0,
1 vT is a left eigenvector of A.
Z03_FRAN6598_07_SE_WC 2014/5/8 0:19 page 23 #4
(WC.21)

and so on. Thus we determine that vT C = 0T has a nontrivial solution, that
C is singular, and that the system is not controllable. To show that a nontrivial
vT exists if C is singular requires more development, which we will not give
here (see Kailath, 1980).
We have given two pictures of uncontrollability. Either a mode is
physically disconnected from the input (Fig. WC.2b), or else two parallel subsystems have identical characteristic roots (Fig. WC.2a). The control
engineer should be aware of the existence of a third simple situation, illustrated in Fig. WC.2c, namely, a pole-zero cancellation. Here the problem
is that the mode at s = 1 appears to be connected to the input but is masked
by the zero at s = 1 in the preceding subsystem; the result is an uncontrollable system. This can be confirmed in several ways. First let us look at the
controllability matrix. The system matrices are

1 0
2
A=
, B=
,
1 1
1
so the controllability matrix is
C=[ B

AB ] =
2
1
2
1
(WC.22)
which is clearly singular. The controllability matrix may be computed using

the ctrb command in Matlab: [cc]=ctrb(A,B). If we compute the transfer
function from u to x2 , we find that

s1
1
1
H(s) =
=
.
(WC.23)
s+1 s1
s+1
Because the natural mode at s = 1 disappears from the inputoutput description, it is not connected to the input. Finally, if we consider the PHR
test,

s+1
0
2
[sI A B] =
,
(WC.24)
1 s 1 1
and let s = 1, then we must test the rank of

2 0 2
,
1 0 1
which is clearly less than 2. This result means, again, that the system is
uncontrollable.
Definition WC.4 The asymptotically stable system (A, B) is controllable if
the controllability Gramian, the square symmetric matrix Cg , given by the
solution to the Lyapunov equation
ACg + Cg AT + BBT = 0
(WC.25)
is nonsingular. The controllability Gramian is also the solution to the

following integral equation:

T
e A BBT e A d .
(WC.26)
Cg =
0
Z03_FRAN6598_07_SE_WC 2014/5/8 0:19 page 24 #5
WC.2 Observability
25
One physical interpretation of the controllability Gramian is that, if the

input to the system is white Gaussian noise, Cg is the covariance of the state.
The controllability Gramian (for an asymptotically stable system) can be
computed with the following command in Matlab: [cg]=gram(A,B).
In conclusion, the four definitions for controllabilitypole assignment
(Definition WC.1), state reachability (Definition WC.2), mode coupling
to the input (Definition WC.3), and controllability Gramian (Definition WC.4)are equivalent. The tests for any of these four properties are
found in terms of the rank of the controllability, controllability Gramian
matrices, or the rank of the matrix pencil [sI A B]. If C is nonsingular, we
can assign the closed-loop poles arbitrarily by state feedback, we can move
the state to any point in the state space in a finite time, and every mode is
connected to the control input.2
WC.2 Observability
So far we have discussed only controllability. The concept of observability
is parallel to that of controllability, and all of the results we have discussed
thus far may be transformed to statements about observability by invoking
the property of duality, as discussed in Section 7.7.2. The observability
definitions are analogous to those for controllability.
Definition WC.1: The system (A, C) is observable if, for any nth-order
polynomial e (s), there exists an estimator gain L such that the characteristic
equation of the state estimator error is e (s).
Definition WC.2: The system (A, C) is observable if, for any x(0), there is a
finite time such that x(0) can be determined (uniquely) from u(t) and y(t)
for 0 t .
Definition WC.3: The system (A, C) is observable if every dynamic mode
in A is connected to the output through C.
Definition WC.4: The asymptotically stable system (A, C) is observable if
the observability Gramian is nonsingular.
As we saw in the discussion for controllability, mathematical tests can
be developed for observability. The system is observable if the observability
matrix
C
CA
O=
(WC.27)
..
.
CAn1
is nonsingular. If we take the transpose of O and let CT = B and
AT = A, then we find the controllability matrix of (A, B), which is another
manifestation of duality. The observability matrix O may be computed
2 We have shown the latter for diagonal A only, but the result is true in general.
Z03_FRAN6598_07_SE_WC 2014/5/8 0:19 page 25 #6

using the obsv command in Matlab: [oo]=obsv(A,C). The system (A, C)
is observable if the following matrix pencil is full rank for all s:

sI A
rank
= n.
(WC.28)
C
The observability Gramian Og , which is a symmetric matrix, and the solution
to the integral equation

T
e A CT Ce A d ,
(WC.29)
Og =
0
as well as the Lyapunov equation

AT Og + Og A + CT C = 0,
(WC.30)
also can be computed (for an asymptotically stable system) using the gram
command in Matlab: [og]=gram(A,C). The observability Gramian has an
interpretation as the information matrix in the context of estimation.
Z03_FRAN6598_07_SE_WC 2014/5/8 0:19 page 26 #7
Appendix WD
Ackermanns Formula
for Pole Placement
Given the plant and state-variable equation
x = Ax + Bu,
(WD.1)
our objective is to find a state-feedback control law

u = Kx
(WD.2)
such that the closed-loop characteristic polynomial is

c (s) = det(sI A + BK).
(WD.3)
First we have to select c (s), which determines where the poles are to be
shifted; then we have to find K such that Eq. (WD.3) will be satisfied. Our
technique is based on transforming the plant equation into control canonical
form.
We begin by considering the effect of an arbitrary nonsingular transformation of the state, as
x = Tx,
(WD.4)
where x is the new transformed state. The equations of motion in the new
coordinates, from Eq. (WD.4), are
x = Tx = Ax + Bu = ATx + Bu,
x = T
ATx + T
x + Bu.
Bu = A
(WD.5)
(WD.6)
Now the controllability matrix for the original state,

Cx = [ B
AB
A2 B
An1 B ],
(WD.7)
provides a useful transformation matrix. We can also define the controllability matrix for the transformed state:
A
B
Cx = [ B
2 B
A
An1 B ].
(WD.8)
The two controllability matrices are related by

Cx = [ T1 B
T1 ATT1 B
] = T1 Cx
(WD.9)
and the transformation matrix

T = Cx Cx1
.
(WD.10)
From Eqs. (WD.9) and (WD.10), we can draw some important conclusions. From Eq. (WD.9), we see that if Cx is nonsingular, then for
any nonsingular T, Cx is also nonsingular. This means that a similarity
27
Z04_FRAN6598_07_SE_WD 2014/5/8 0:22 page 27 #1
28 Appendix WD Ackermanns Formula for Pole Placement

transformation on the state does not change the controllability properties
of a system. We can look at this in another way. Suppose we would like
to find a transformation to take the system (A, B) into control canonical
form. As we shall shortly see, Cx in that case is always nonsingular. From
Eq. (WD.9), we see that a nonsingular T will always exist if and only if Cx
is nonsingular. We derive the following theorem.
Theorem WD.1 We can always transform (A, B) into control canonical
form if and only if the system is controllable.
Let us take a closer look at control canonical form and treat the thirdorder case, although the results are true for any nth-order case:
a1 a2 a3
1
= Ac = 1
0
0 , B = Bc = 0 . (WD.11)
A
0
1
0
0
The controllability matrix, by direct computation, is
1 a1 a12 a2
Cx = Cc = 0
1
a1 .
0
0
1
(WD.12)
Because this matrix is upper triangular with ones along the diagonal, it is
always invertible. Also note that the last row of Cx is the unit vector with all
zeros, except for the last element, which is unity. We shall use this fact in
what we do next.
As we pointed out in Section 7.5, the design of a control law for the
state x is trivial if the equations of motion happen to be in control canonical
form. The characteristic equation is
s3 + a1 s2 + a2 s + a3 = 0,
(WD.13)
and the characteristic equation for the closed-loop system comes from
Acl = Ac Bc Kc
(WD.14)
and has the coefficients shown:

s3 + (a1 + Kc1 )s2 + (a2 + Kc2 )s + (a3 + Kc3 ) = 0.
(WD.15)
To obtain the desired closed-loop pole locations, we must make the

coefficients of s in Eq. (WD.15) match those in
c (s) = s3 + 1 s2 + 2 s + 3 ,
(WD.16)
so
a1 + Kc1 = 1 ,
a2 + Kc2 = 2 ,
a3 + Kc3 = 3 ,
(WD.17)
or in vector form,
a + Kc = ,
(WD.18)
where a and are row vectors containing the coefficients of the characteristic
polynomials of the open-loop and closed-loop systems, respectively.
Z04_FRAN6598_07_SE_WD 2014/5/8 0:22 page 28 #2
Appendix WD Ackermanns Formula for Pole Placement
29
We now need to find a relationship between these polynomial coefficients and the matrix A. The requirement is achieved by the Cayley
Hamilton theorem, which states that a matrix satisfies its own characteristic
polynomial. For Ac , this means that
Anc + a1 An1
+ a2 An2
+ + an I = 0.
c
c
(WD.19)
Now suppose we form the polynomial c (A), which is the closed-loop

characteristic polynomial with the matrix A substituted for the complex
variable s:
c (Ac ) = Anc + 1 An1
+ 2 An2
+ + n I.
c
c
If we solve Eq. (WD.19) for
that
Anc
(WD.20)
and substitute into Eq. (WD.20), we find
c (Ac ) = (a1 + 1 )An1

+ (a2 + 2 )An2
+ + (n + n )I.
c
c
(WD.21)
But, because Ac has such a special structure, we observe that if we multiply
it by the transpose of the nth unit vector, enT = [ 0 0 1 ], we get
enT Ac = [ 0
T
0 ] = en1
,
(WD.22)
as we can see from Eq. (WD.11). If we multiply this vector again by Ac ,

getting
(enT Ac )Ac = [ 0
0 ]Ac
=[ 0
T
0 ] = en2
,
(WD.23)
and continue the process, successive unit vectors are generated until

enT An1
= 1 0 0 = e1T .
(WD.24)
c
Therefore, if we multiply Eq. (WD.21) by enT , we find that
enT c (Ac ) = (a1 + 1 )e1T + (a2 + 2 )e2T + + (an + n )enT
= [ Kc1
Kc2
Kcn ] = Kc ,
(WD.25)
where we use Eq. (WD.18), which relates Kc to a and .

We now have a compact expression for the gains of the system in control
canonical form as represented in Eq. (WD.25). However, we still need the
expression for K, which is the gain on the original state. If u = Kc x , then
u = Kc T1 x, so that
K = Kc T1 = enT c (Ac )T1
= enT c (T1 AT)T1
= enT T1 c (A).
(WD.26)
In the last step of Eq. (WD.26), we used the fact that (T1 AT)k = T1 Ak T
and that c is a polynomial, that is, a sum of the powers of Ac . From
Eq. (WD.9), we see that
T1 = Cc Cx1 .
(WD.27)
Z04_FRAN6598_07_SE_WD 2014/5/8 0:22 page 29 #3
30 Appendix WD Ackermanns Formula for Pole Placement

With this substitution, Eq. (WD.26) becomes
K = enT Cc Cx1 c (A).
Ackermanns formula
(WD.28)
Now, we use the observation made earlier for Eq. (WD.12) that the last row
of Cc , which is enT Cc , is again enT . We finally obtain Ackermanns formula:
K = enT Cx1 c (A).
(WD.29)
We note again that forming the explicit inverse of Cx is not advisable

for numerical accuracy. Thus we need to solve bT such that
enT Cx1 = bT .
(WD.30)
We solve the linear set of equations

bT Cx = enT
(WD.31)
K = bT c (A).
(WD.32)
and then compute

Ackermanns formula, Eq. (WD.29), even though elegant, is not recommended for systems with a large number of state variables. Even if it is used,
Eqs. (WD.31) and (WD.32) are recommended for better numerical accuracy.
Z04_FRAN6598_07_SE_WD 2014/5/8 0:22 page 30 #4
Appendix W2.1.4
Complex Mechanical Systems
In some cases, mechanical systems contain both translational and rotational
portions. The procedure is the same as that described in Section 2.1: sketch
the free-body diagrams, define coordinates and positive directions, determine all forces and moments acting, and apply Eqs. (2.1) and/or (2.14).
EXAMPLE W2.1
Rotational and Translational Motion: Hanging Crane

Write the equations of motion for the hanging crane pictured in Fig. W2.1
and shown schematically in Fig. W2.2. Linearize the equations about
Figure W2.1
Crane with a hanging
load
Source: Photo courtesy of
Harnischfeger Corporation,
Milwaukee, Wisconsin
31
Z05_FRAN6598_07_SE_W2.1.4 2014/5/13 9:13 page 31 #1
32 Appendix W2.1.4 Complex Mechanical Systems

= 0, which would typically be valid for the hanging crane. Also linearize
the equations for = , which represents the situation for the inverted
pendulum shown in Fig. W2.3.
Solution. A schematic diagram of the hanging crane is shown in Fig. W2.2,
while the free-body diagrams are shown in Fig. W2.4. In the case of the
pendulum, the forces are shown with bold lines, while the components of
the inertial acceleration of its center of mass are shown with dashed lines.
Because the pivot point of the pendulum is not fixed with respect to an inertial reference, the rotation of the pendulum and the motion of its mass center
must be considered. The inertial acceleration needs to be determined because
the vector a in Eq. (2.1) is given with respect to an inertial reference. The
inertial acceleration of the pendulums mass center is the vector sum of the
three dashed arrows shown in Fig. W2.4b. The derivation of the components
of an objects acceleration is called kinematics and is usually studied as a
prelude to the application of Newtons laws. The results of a kinematic study
are shown in Fig. W2.4b. The component of acceleration along the pendulum
is l 2 and is called the centripetal acceleration. It is present for any object
whose velocity is changing direction. The x -component of acceleration is a
consequence of the pendulum pivot point accelerating at the trolleys acceleration and will always have the same direction and magnitude as those of
the trolleys. The l component is a result of angular acceleration of the
pendulum and is always perpendicular to the pendulum.
These results can be confirmed by expressing the center of mass of the
pendulum as a vector from an inertial reference and then differentiating that
vector twice to obtain an inertial acceleration. Figure W2.4c shows i and j
Figure W2.2
Schematic of the crane
with hanging load
mt
I, mp
u
Figure W2.3
Inverted pendulum
I, mp
u
mt
Z05_FRAN6598_07_SE_W2.1.4 2014/5/13 9:13 page 32 #2
Appendix W2.1.4 Complex Mechanical Systems

P
x
x
33
lu2
lu
x
I, mp
u
mt
bx, friction
N
Reaction
P
force from
pendulum
u
u
mpg
(a)
(b)
(c)
Figure W2.4
Hanging crane: (a) free-body diagram of the trolley; (b) free-body diagram of the pendulum; (c) position
vector of the pendulum
axes that are inertially fixed and a vector r describing the position of the
pendulum center of mass. The vector can be expressed as
r = xi + l(i sin j cos ).
The first derivative of r is
i cos + j sin ).
r = x i + l (
Likewise, the second derivative of r is
i cos + j sin ) l 2 (i sin j cos ).
r = x i + l (
Note that the equation for r confirms the acceleration components shown in
Fig.W2.4b. The l 2 term is aligned along the pendulum pointing toward the
axis of rotation, and the l term is aligned perpendicular to the pendulum
pointing in the direction of a positive rotation.
Having all the forces and accelerations for the two bodies, we now
proceed to apply Eq. (2.1). In the case of the trolley, Fig.W2.4a, we see that
it is constrained by the tracks to move only in the x-direction; therefore,
application of Eq. (2.1) in this direction yields
mt x + bx = u N,
(W2.1)
where N is an unknown reaction force applied by the pendulum. Conceptually, Eq. (2.1) can be applied to the pendulum of Fig. W2.4b in the vertical
and horizontal directions, and Eq. (2.14) can be applied for rotational motion
to yield three equations in the three unknowns: N, P, and . These three equations then can be manipulated to eliminate the reaction forces N and P so
that a single equation results describing the motion of the pendulumthat
Z05_FRAN6598_07_SE_W2.1.4 2014/5/13 9:13 page 33 #3

is, a single equation in . For example, application of Eq. (2.1) for pendulum
motion in the x-direction yields
N = mp x + mp l cos mp l 2 sin .
(W2.2)
However, considerable algebra will be avoided if Eq. (2.1) is applied

perpendicular to the pendulum to yield
P sin + N cos mp g sin = mp l + mp x cos .
(W2.3)
Application of Eq. (2.14) for the rotational pendulum motion, for which the
moments are summed about the center of mass, yields
Pl sin Nl cos = I ,
(W2.4)
where I is the moment of inertia about the pendulums mass center. The
reaction forces N and P can now be eliminated by combining Eqs. (W2.3)
and (W2.4). This yields the equation
(I + mp l 2 ) + mp gl sin = mp lx cos .
(W2.5)
It is identical to a pendulum equation of motion, except that it contains a

forcing function that is proportional to the trolleys acceleration.
An equation describing the trolley motion was found in Eq. (W2.1), but
it contains the unknown reaction force N. By combining Eqs. (W2.2) and
(W2.1), N can be eliminated to yield
(mt + mp )x + bx + mp l cos mp l 2 sin = u.
(W2.6)
Equations (W2.5) and (W2.6) are the nonlinear differential equations that
describe the motion of the crane with its hanging load. For an accurate
calculation of the motion of the system, these nonlinear equations need to
be solved.
To linearize the equations for small motions about = 0, let cos
= 1,
sin
= , and 2
= 0; thus the equations are approximated by
(I + mp l 2 ) + mp gl = mp lx ,
(mt + mp )x + bx + mp l = u.
(W2.7)
Neglecting the friction term b leads to the transfer function from the
control input u to the hanging crane angle :
mp l
(s)
=
.
2
U(s)
((I + mp l )(mt + mp ) mp2 l 2 )s2 + mp gl(mt + mp )
Z05_FRAN6598_07_SE_W2.1.4 2014/5/13 9:13 page 34 #4
(W2.8)
W2.1 Additional Problems for Translational and Rotational Systems
Inverted pendulum
equations
35
, assume =
For the inverted pendulum in Fig. W2.3, where =
+ , where represents motion from the vertical upward direction. In this
case, cos
= 1, sin
= in Eqs. (W2.5) and (W2.6), and Eqs. (W2.7)
1
become
(I + mp l 2 ) mp gl = mp lx ,
(mt + mp )x + bx mp l = u.
(W2.9)
As noted in Example 2.2, a stable system will always have the same
signs on each variable, which is the case for the stable hanging crane modeled
by Eqs. (W2.7). However, the signs on and in the top Eq. (W2.9) are
opposite, thus indicating instability, which is the characteristic of the inverted
pendulum.
The transfer function, again without friction, is
mp l
(s)
=
.
(W2.10)
U(s)
((I + mp l 2 ) mp2 l 2 )s2 mp gl(mt + mp )
W2.1 Additional Problems for Translational and

Rotational Systems
Assume the driving force on the hanging crane of Fig. W2.2 is provided by a
motor mounted on the cab with one of the support wheels connected directly
to the motors armature shaft. The motor constants are Ke and Kt , and the
circuit driving the motor has a resistance Ra and negligible inductance. The
wheel has a radius r. Write the equations of motion relating the applied
motor voltage to the cab position and load angle.
Solution. The dynamics of the hanging crane are given by Eqs. (W2.5) and
(W2.6),

I + mp l 2 + mp gl sin = mp lx cos ,

mt + mp x + bx + mp l cos mp l 2 sin = u,
where x is the position of the cab, is the angle of the load, and u is the
applied force that will be produced by the motor. Our task here is to find the
force applied by the motor. Normally, the rotational dynamics of a motor is
J1 m + b1 m = Tm = Kt ia ,
where the current is found from the motor circuit, which reduces to
Ra ia = Va Ke m
for the case where the inductance is negligible. However, since the motor
is geared directly to the cab, m and x are related kinematically by
x = rm
1 The inverted pendulum is often described with the angle of the pendulum being positive for
clockwise motion. If defined that way, then reverse the sign on all terms in Eq. (W2.9) in
or .
Z05_FRAN6598_07_SE_W2.1.4 2014/5/13 9:13 page 35 #5

and we can neglect any extra inertia or damping from the motor itself compared to the inertia and damping of the large cab. Therefore we can rewrite
the two motor equations in terms of the force applied by the motor on the cab
u = Tm /r = Kt ia /r,
ia = (Va Ke m )/Ra ,
where
m = x /r.
These equations, along with

I + mp l 2 + mp gl sin = mp lx cos ,

mt + mp x + bx + mp l cos mp l 2 sin = u,
constitute the required relations.
Z05_FRAN6598_07_SE_W2.1.4 2014/5/13 9:13 page 36 #6
Appendix W3.2.3
W3.2.3 Masons Rule and the Signal-Flow Graph
Signal-ow graph
A compact alternative notation to the block diagram is given by the signalflow graph introduced by S. J. Mason (1953, 1956). As with the block
diagram, the signal-flow graph offers a visual tool for representing the causal
relationships between the components of the system. The method consists of
characterizing the system by a network of directed branches and associated
gains (transfer functions) connected at nodes. Several block diagrams and
their corresponding signal-flow graphs are shown in Fig. W3.1. The two
ways of depicting a system are equivalent, and you can use either diagram
to apply Masons rule (to be defined shortly).
In a signal-flow graph, the internal signals in the diagram, such as the
common input to several blocks or the output of a summing junction, are
called nodes. The system input point and the system output point are also
nodes; the input node has outgoing branches only, and the output node has
incoming branches only. Mason defined a path through a block diagram as
a sequence of connected blocks, the route passing from one node to another
in the direction of signal flow of the blocks without including any block more
than once. A forward path is a path from the input to output such that no
node is included more than once. If the nodes are numbered in a convenient
order, a forward path can be identified by the numbers that are included.
Any closed path that returns to its starting node without passing through any
node more than once is a loop, and a path that leads from a given variable
back to the same variable is a loop path. The path gain is the product of
component gains (transfer functions) making up the path. Similarly, the loop
gain is the path gain associated with a loopthat is, the product of gains in a
loop. If two paths have a common component, they are said to touch. Notice
particularly in this connection that the input and the output of a summing
junction are not the same and that the summing junction is a one-way device
from its inputs to its output.
Masons rule relates the graph to the algebra of the simultaneous equations it represents.1 Consider Fig. W3.1c, where the signal at each node has
been given a name and the gains are marked. Then the block diagram (or the
signal-flow graph) represents the system of equations:
X1 (s) = X3 (s) + U(s),
X2 (s) = G1 (s)X1 (s) + G2 (s)X2 (s) + G4 (s)X3 (s),
Y (s) = 1X3 (s).
1 The derivation is based on Cramers rule for solving linear equations by determinants and is
described in Masons papers.
37
Z06_FRAN6598_07_SE_W3.2.3 2014/5/6 3:01 page 37 #1
38 Appendix W3.2.3
G3
G3
1
+
G1
G2
G1
G2
G4
G4
(a)
G3
0
+
G1
G3
3
G2
0 1
G2
G1
G4
G4
(b)
G4
X1
+
G1
X2
G4
X3 Y
G3
U 1
G1
1 Y
G2
G3
G2
1
1
(c)
Figure W3.1
Block diagrams and corresponding signal-ow graphs
Masons rule
Masons rule states that the inputoutput transfer function associated with a
signal-flow graph is given by
G(s) =
1
Y (s)
=
G i i ,
U(s)

i
where
Gi = path gain of the ith forward path,
= the system determinant

=1
(all individual loop gains) + (gain products of all possible

two loops that do not touch) (gain products of all possible three
loops that do not touch) + . . . ,
Z06_FRAN6598_07_SE_W3.2.3 2014/5/6 3:01 page 38 #2
Appendix W3.2.3
39
i = ith forward path determinant

= value of for that part of the block diagram that does not touch the
ith forward path.
We will now illustrate the use of Masons rule using some examples.
EXAMPLE W3.1
Masons Rule in a Simple System

Find the transfer function for the block diagram in Fig. W3.2.
Solution. From the block diagram shown in Fig. W3.2, we have
Forward Path
1236
12346
123456
232
2342
23452
Path Gain

1
G1 = 1
(b1 )(1)
s

1
1
(b2 )(1)
G2 = 1
s
s

1
1
1
(b3 )(1)
G3 = 1
s
s
s
Loop Path Gain
l1 = a1 /s
l2 = a2 /s2
l3 = a3 /s3
and the determinants are

a
a3
a2
1
= 1 2 3 + 0,
s
s
s
1 = 1 0,
2 = 1 0,
3 = 1 0.
Applying Masons rule, we find the transfer function to be
G(s) =
Y (s)
(b1 /s) + (b2 /s2 ) + (b3 /s3 )
=
U(s)
1 + (a1 /s) + (a2 /s2 ) + (a3 /s3 )
s3
b1 s2 + b2 s + b3
.
+ a1 s 2 + a2 s + a 3
Z06_FRAN6598_07_SE_W3.2.3 2014/5/6 3:01 page 39 #3
40 Appendix W3.2.3
b1
b2
U(s)
1 +
+
2
+
+
1
s
3
x1c
1
s
1
s
x2c
5
x3c
+ 6
Y(s)
b3
- a1
- a2
- a3
Figure W3.2
Block diagram for Example W3.1
Masons rule is particularly useful for more complex systems where
there are several loops, some of which do not sum into the same point.
EXAMPLE W3.2
Masons Rule in a Complex System

Find the transfer function for the system shown in Fig. W3.3.
Solution. From the block diagram, we find that
Forward Path
12456
1236
242
454
565
236542
Path Gain
G 1 = H1 H 2 H 3
G2 = H4
Loop Path Gain
l1 = H1 H5 (does not touch l3 )
l2 = H2 H6
l3 = H3 H7 (does not touch l1 )
l4 = H4 H7 H6 H5
and the determinants are

= 1 (H1 H5 + H2 H6 + H3 H7 + H4 H7 H6 H5 ) + (H1 H5 H3 H7 ),
1 = 1 0,
2 = 1 H 2 H6 .
Therefore,
Y (s)
H1 H2 H3 + H4 H4 H2 H6
G(s) =
=
.
U(s)
1 H 1 H5 H 2 H 6 H 3 H7 H 4 H7 H6 H5 + H 1 H5 H3 H7
Masons rule is useful for solving relatively complicated block diagrams
by hand. It yields the solution in the sense that it provides an explicit input
output relationship for the system represented by the diagram. The advantage
Z06_FRAN6598_07_SE_W3.2.3 2014/5/6 3:01 page 40 #4
Appendix W3.2.3
H4
41
H6
U
1 +
+
H1
H2
+
H3
H5
H7
Figure W3.3
Block diagram for Example W3.2
compared with path-by-path block-diagram reduction is that it is systematic
and algorithmic rather than problem dependent. MATLAB and other control
systems computer-aided software allow you to specify a system in terms of
individual blocks in an overall system, and the software algorithms perform
the required block-diagram reduction; therefore, Masons rule is less important today than in the past. However, there are some derivations that rely
on the concepts embodied by the rule, so it still has a role in the control
designers toolbox.
Z06_FRAN6598_07_SE_W3.2.3 2014/5/6 3:01 page 41 #5
Appendix W.3.6.3.1
Routh Special Cases
Special case I
If only the first element in one of the rows is zero, then we can consider
a modified equation with one of the coefficients perturbed by > 0 and
applying the test by taking the limit as 0.
EXAMPLE W3.3
Rouths Test for Special Case I

Consider the polynomial
a(s) = s5 + 3s4 + 2s3 + 6s2 + 6s + 9.
Determine whether any of the roots are in the RHP.
Solution. In this example, let the coefficient of s3 be 2 + . The test follows
from there. The Routh array is
2 6
s5 : 1
s4 : 3
6 9
s3 :
3 0
s2 : 69
9
0

3 2
s: 3 23 0 0
s0 : 9
0
There are two sign changes in the first column of the array, which means
there are two poles not in the LHP.1
Special case II
Another special2 case occurs when an entire row of the Routh array is
zero. This indicates that there are complex conjugate pairs of roots that are
mirror images of each other with respect to the imaginary axis. To apply
Rouths test correctly, we follow the ensuing procedure. If the ith row is
zero, we form an auxiliary equation from the previous (nonzero) row:
a1 (s) = 1 si+1 + 2 si1 + 3 si3 + .
(W3.1)
Here {i } are the coefficients of the (i+1)th row in the array. We then replace
the ith row by the coefficients of the derivative of the auxiliary polynomial
and complete the array. However, the roots of the auxiliary polynomial in
Eq. (W3.1) are also roots of the characteristic equation, and these must be
tested separately.
1 The actual roots computed with Matlab are at 2.9043, 0.6567 1.2881j, 0.7046 0.9929j.
2 Special case II.
42
Z07_FRAN6598_07_SE_W3.6.3.1 2014/5/13 9:03 page 42 #1
Appendix W.3.6.3.1 Routh Special Cases
EXAMPLE W3.4
43
Routh Test for Special Case II

For the polynomial
a(s) = s5 + 5s4 + 11s3 + 23s2 + 28s + 12,
determine whether there are any roots on the j axis or in the RHP.
Solution. The Routh array is
s5 :
s4 :
s3 :
s2 :
s:
New s :
s0 :
1
5
6.4
3
0
6
12.
11
23
25.6
12
0
0
28
12
0
a1 (s) = 3s2 + 12
1 (s)
= 6s
dads
There are no sign changes in the first column. Hence all the roots have
negative real parts except for a pair on the imaginary axis. We may deduce
this as follows: When we replace the zero in the first column by > 0, there
are no sign changes. If we let < 0, then there are two sign changes. Thus,
if = 0, there are two poles on the imaginary axis, which are the roots of
a1 (s) = 3s2 + 12 = 0,
or
s = j2.
This agrees with the fact that the actual roots are at 3, 2j, 1, and 1,
as computed using the roots command in Matlab.
Z07_FRAN6598_07_SE_W3.6.3.1 2014/5/13 9:03 page 43 #2
Appendix W3.7
System Identication
W3.7.1 A Perspective on System Identication
In order to design controls for a dynamic system, it is necessary to have a
model that will adequately describe the systems dynamics. The information
available to the designer for this purpose is typically of three kinds.
1. Physical model: First, there is the knowledge of physics, chemistry,
biology, and the other sciences which have over the years developed equations of motion to explain the dynamic response of rigid and flexible bodies,
electric circuits and motors, fluids, chemical reactions, and many other constituents of systems to be controlled. The model based on this knowledge
is referred to as a physical model. There are many advantages to this
approach, including ease of controller development and testing. One disadvantage of this approach is that a fairly high-fidelity physical model must be
developed.
2. Black box model: It is often the case that for extremely complex
physical phenomena the laws of science are not adequate to give a satisfactory
description of the dynamic plant that we wish to control. Examples include
the force on a moving airplane caused by a control surface mounted on a wing
and the heat of combustion of a fossil fuel of uncertain composition. In these
circumstances, the designer turns to data taken from experiments directly
conducted to excite the plant and measure its response. The second approach
uses an empirical or heuristic model referred to as the black box model.
In this approach, the control engineer injects open-loop commands into the
system and records the sensor response. The process of constructing models from experimental data is called system identification. Standard system
identification techniques (for example, linear least-squares) are used to identify a dynamic input/output model. The advantage of this technique is that the
control engineer does not need to have a deep understanding of how the system physically behaves, but instead can design a controller solely based on
the derived model. There are several major disadvantages to this approach.
First, the control engineer must have access to working hardware. Another
serious disadvantage of this approach is that it does not provide insight or
physical understanding of how specific hardware modifications will affect
the controlusually hardware modifications require the control engineer to
repeat the full cycle of system identification, control design, and validation.
The advantage of this approach is that we use logic and data to model inputs
and outputs and the detailed knowledge of the physics is not required.
3. Grey box model: The third approach is the use of the combination
of physical and empirical models referred to as grey box modeling.
44
Z08_FRAN6598_07_SE_W3.6.4 2014/5/13 9:14 page 44 #1
Appendix W3.7 System Identication
Parametric model
Nonparametric model
45
In identifying models for control, our motivation is very different from

that of modeling as practiced in the sciences. In science, one seeks to develop
models of nature as it is; in control, one seeks to develop models of the plant
dynamics that will be adequate for the design of a controller that will cause
the actual dynamics to be stable and to give good performance. The initial
design of a control system typically considers a small signal analysis and
is based on models that are linear and time-invariant (LTI). This is referred
to as a control relevant model. Having accepted that the model is to be
linear, we still must choose between several alternative descriptions of linear
systems. If we examine the design methods described in the earlier chapters,
we find that the required plant models may be grouped in two categories:
parametric and nonparametric. For design via root locus or pole assignment,
we require a parametric description such as a transfer function or a statevariable description from which we can obtain the poles and zeros of the
plant. These equivalent models are completely described by the numbers
that specify the coefficients of the polynomials, the elements of the statedescription matrices, or the numbers that specify the poles and zeros. In
either case, we call these numbers the parameters of the model, and the
category of such models is a parametric description of the plant model.
In contrast to parametric models, the frequency-response methods of
Nyquist, Bode, and Nichols require the curves of amplitude and phase of
the transfer function G(j) = Y (j)/U(j) as functions of . Clearly, if
we happen to have a parametric description of the system, we can compute
the transfer function and the corresponding frequency response. However,
if we are given the frequency response or its inverse transform, the impulse
response, without parameters (perhaps obtained from experimental data),
we have all we need to design a lead, lag, notch, or other compensation
to achieve a desired bandwidth, phase margin, or other frequency response
performance objective without ever knowing what the parameters are. We
call the functional curves of G(j) a nonparametric model because in
principle there is no finite set of numbers that describes it exactly.
Because of the large data records necessary to obtain effective models
and the complexity of many of the algorithms used, the use of computer
aids is essential in identification. Developments such as Matlabs System
Identification Toolbox are enormous aids to the practical use of the system
identification techniques. For detailed discussion on system identification,
the reader is referred to Franklin, Powell, and Workman (1998).
W3.7.2 Obtaining Models from Experimental Data

There are several reasons for using experimental data to obtain a model of
the dynamic system to be controlled. In the first place, the best theoretical
model built from equations of motion is still only an approximation of reality.
Sometimes, as in the case of a very rigid spacecraft, the theoretical model
is extremely good. Other times, as with many chemical processes such as
papermaking or metalworking, the theoretical model is very approximate.
In every case, before the final control design is done, it is important and
Z08_FRAN6598_07_SE_W3.6.4 2014/5/13 9:14 page 45 #2
46 Appendix W3.7 System Identication
Our sources of
experimental data
Transient response
Frequency response
prudent to verify the theoretical model with experimental data. Second, in

situations for which the theoretical model is especially complicated or the
physics of the process is poorly understood, the only reliable information on
which to base the control design is the experimental data. Finally, the system
is sometimes subject to online changes that occur when the environment of
the system changes. Examples include when an aircraft changes altitude
or speed, a paper machine is given a different composition of fiber, or a
nonlinear system moves to a new operating point. On these occasions, we
need to retune the controller by changing the control parameters. This
requires a model for the new conditions, and experimental data are often the
most effective, if not the only, information available for the new model.
There are four kinds of experimental data for generating a model:
1. Transient response, such as comes from an impulse or a step;
2. Frequency-response data, which result from exciting the system with
sinusoidal inputs at many frequencies;
3. Stochastic steady-state information, as might come from flying an
aircraft through turbulent weather or from some other natural source of
randomness; and
4. Pseudorandom-noise data, as may be generated in a digital computer.
Each class of experimental data has its properties, advantages, and
disadvantages.
Transient-response data are quick and relatively easy to obtain. They
are also often representative of the natural signals to which the system is
subjected. Thus a model derived from such data can be reliable for designing
the control system. On the other hand, in order for the signal-to-noise ratio
to be sufficiently high, the transient response must be highly noticeable.
Consequently, the method is rarely suitable for normal operations, so the
data must be collected as part of special tests. A second disadvantage is that
the data do not come in a form suitable for standard control systems designs,
and some parts of the model, such as poles and zeros, must be computed
from the data.1 This computation can be simple in special cases or complex
in the general case.
Frequency-response data (see Chapter 6) are simple to obtain but substantially more time consuming than transient-response information. This is
especially so if the time constants of the process are large, as often occurs
in chemical processing industries. As with the transient-response data, it
is important to have a good signal-to-noise ratio, so obtaining frequencyresponse data can be very expensive. On the other hand, as we will see
in Chapter 6, frequency-response data are exactly in the right form for
frequency-response design methods; so once the data have been obtained,
the control design can proceed immediately.
1 Ziegler and Nichols (1943), building on the earlier work of Callender et al. (1936), use the
step response directly in designing the controls for certain classes of processes. See Chapter 4
for details.
Z08_FRAN6598_07_SE_W3.6.4 2014/5/13 9:14 page 46 #3

Stochastic steady-state
Pseudorandom noise
(PRBS)
47
Normal operating records from a natural stochastic environment at first

appear to be an attractive basis for modeling systems, since such records
are by definition nondisruptive and inexpensive to obtain. Unfortunately, the
quality of such data is inconsistent, tending to be worse just when the control
is best, because then the upsets are minimal and the signals are smooth. At
such times, some or even most of the system dynamics are hardly excited.
Because they contribute little to the system output, they will not be found
in the model constructed to explain the signals. The result is a model that
represents only part of the system and is sometimes unsuitable for control.
In some instances, as occurs when trying to model the dynamics of the
electroencephalogram (brain waves) of a sleeping or anesthetized person to
locate the frequency and intensity of alpha waves, normal records are the
only possibility. Usually they are the last choice for control purposes.
Finally, the pseudorandom signals that can be constructed using digital
logic have much appeal. Especially interesting for model making is the
pseudorandom binary signal (PRBS). The PRBS takes on the value +A or
A according to the output (1 or 0) of a feedback shift register. The feedback
to the register is a binary sum of various states of the register that have been
selected to make the output period (which must repeat itself in finite time)
as long as possible. For example, with a register of 20 bits, 220 1 (over
a million) steps are produced before the pattern repeats. Analysis beyond
the scope of this text has revealed that the resulting signal is almost like
a broadband random signal. Yet this signal is entirely under the control of
the engineer who can set the level (A) and the length (bits in the register)
of the signal. The data obtained from tests with a PRBS must be analyzed
by computer and both special-purpose hardware and programs for generalpurpose computers have been developed to perform this analysis.
W3.7.3 Models from Transient-Response Data

To obtain a model from transient data, we assume that a step response is
available. If the transient is a simple combination of elementary transients,
then a reasonable low-order model can be estimated using hand calculations.
For example, consider the step response shown in Fig. W3.1. The response is
monotonic and smooth. If we assume that it is given by a sum of exponentials,
we can write
y(t) = y() + Aet + Bet + Ce t + .
Figure W3.4
A step response
characteristic of many
chemical processes
y(t)
1.0
Z08_FRAN6598_07_SE_W3.6.4 2014/5/13 9:14 page 47 #4
(W3.2)

Subtracting off the final value and assuming that is the slowest pole, we
write
y y()
= Aet ,
log [y y()]
= log A t log e,
10
10
10
(W3.3)
= log10 A 0.4343t.
This is the equation of a line whose slope determines and intercept determines A. If we fit a line to the plot of log10 [y y()] (or log10 [y() y]
if A is negative), then we can estimate A and . Once these are estimated,
we plot y [y() + Aet ], which as a curve approximates Bet and on
the log plot is equivalent to log10 B 0.4345t. We repeat the process, each
time removing the slowest remaining term, until the data stop is accurate.
Then we plot the final model step response and compare it with data so we
can assess the quality of the computed model. It is possible to get a good fit
to the step response and yet be far off from the true time constants (poles)
of the system. However, the method gives a good approximation for control
of processes whose step responses look like Fig. W3.4.
EXAMPLE W3.5
Determining the Model from Time-Response Data

Find the transfer function that generates the data given in Table W3.1 and
plotted in Fig. W3.5.
Solution. Table W3.1 shows and Fig. W3.5 implies that the final value of
the data is y() = 1. We know that A is negative because y() is greater
than y(t). Therefore, the first step in the process is to plot log10 [y() y],
which is shown in Fig. W3.6. From the line (fitted by eye), the values are
log10 |A| = 0.125,
0.4343 =
Thus
1.602 1.167
0.435
=
= 1.
t
1
A = 1.33,
= 1.0.
TABLE W3.1
Step Response Data

t
y(t)
y(t)
0.1
0.1
0.2
0.3
0.4
0.5
0.000
0.005
0.034
0.085
0.140
0.215
1.0
1.5
2.0
2.5
3.0
4.0
0.510
0.700
0.817
0.890
0.932
0.975
1.000
Sinha and Kuszta (1983).
Z08_FRAN6598_07_SE_W3.6.4 2014/5/13 9:14 page 48 #5
49
1.0
Figure W3.5
Step response data in
Table W3.1
0.8
y(t)
0.6
0.4
0.2
0
-0.2
Figure W3.6
log10 [y() y]
versus t
3
Time (sec)
0.5
0.125
log10 |1 - y|
-0.5
-1.0
-1.167
-1.5
-1.602
-2.0
Time (sec)
If we now subtract 1 + Aet from the data and plot the log of the result,
we find the plot of Fig. W3.7. Here we estimate
log10 B = 0.48,
0.4343 =
0.48 (1.7)
= 2.5,
0.5
Z08_FRAN6598_07_SE_W3.6.4 2014/5/13 9:14 page 49 #6

Figure W3.7
log10 [y (1 + Aet )]
versus t
-0.48
log10 Be-bt
-0.5
-1.0
-1.5
-1.7
-2.0
0.1
0.2
0.3
0.4
0.5
Time (sec)
Figure W3.8
Model ts to the
experimental data
1.0
0.8
y(t)
0.6
0.4
A = -1.33, a = 1, B = 0.33, b = 5.8
0.2
A = -1.37, a = 1, B = 0.37, b = 4.3
0
Data
-0.2
3
Time (sec)
= 5.8,
B = 0.33.
Combining these results, we arrive at the y estimate
y (t)
= 1 1.33et + 0.33e5.8t .
(W3.4)
Equation (W3.4) is plotted as the colored line in Fig. W3.8 and shows a
reasonable fit to the data, although some error is noticeable near t = 0.
From y (t), we compute
1
1.33
0.33
Y (s) =
+
s
s + 1 s + 5.8
(s + 1)(s + 5.8) 1.33s(s + 5.8) + 0.33s(s + 1)
=
s(s + 1)(s + 5.8)
0.58s + 5.8
=
.
s(s + 1)(s + 5.8)
Z08_FRAN6598_07_SE_W3.6.4 2014/5/13 9:14 page 50 #7
51
The resulting transfer function is

G(s) =
0.58(s 10)
.
(s + 1)(s + 5.8)
Notice that this method has given us a system with a zero in the RHP,
even though the data showed no values of y that were negative. Very small
differences in the estimated value for A, all of which approximately fit the
data, can cause values of to range from 4 to 6. This illustrates the sensitivity
of pole locations to the quality of the data and emphasizes the need for a
good signal-to-noise ratio.
By using a computer to perform the plotting, we are better able to iterate
the four parameters to achieve the best overall fit. The data presentation in
Figs. W3.6 and W3.7 can be obtained directly by using a semilog plot. This
eliminates having to calculate log10 and the exponential expression to find
the values of the parameters. The equations of the lines to be fit to the data are
y(t) = Aet and y(t) = Bet , which are straight lines on a semilog plot. The
parameters A and , or B and , are iteratively selected so that the straight
line comes as close as possible to passing through the data. This process
produces the improved fit shown by the dashed black line in Fig. W3.8. The
revised parameters, A = 1.37, B = 0.37, and = 4.3 result in the transfer
function
G(s) =
0.22s + 4.3
.
(s + 1)(s + 4.3)
The RHP zero is still present, but it is now located at s

= +20 and has no
noticeable effect on the time response.
This set of data was fitted quite well by a second-order model. In many
cases, a higher-order model is required to explain the data and the modes
may not be as well separated.
Least-squares system
identication
If the transient response has oscillatory modes, then these can sometimes
be estimated by comparing them with the standard plots of Fig. 3.18. The
period will give the frequency d , and the decay from one period to the next
will afford an estimate of the damping ratio. If the response has a mixture
of modes not well separated in frequency, then more sophisticated methods
need to be used. One such is least-squares system identification, in which
a numerical optimization routine selects the best combination of system
parameters so as to minimize the fit error. The fit error is defined to be a
scalar cost function

J=
(ydata ymodel )2 , i = 1, 2, 3, , for each data point,
i
so that fit errors at all data points are taken into account in determining the
best value for the system parameters.
Z08_FRAN6598_07_SE_W3.6.4 2014/5/13 9:14 page 51 #8

W3.7.3.1 Models from Other Data
As mentioned early in Section 3.1.2, we can also generate a model using
frequency-response data, which are obtained by exciting the system with a
set of sinusoids and plotting H(j). In Chapter 6, we show how such plots can
be used directly for design. Alternatively, we can use the frequency response
to estimate the poles and zeros of a transfer function using straight-line
asymptotes on a logarithmic plot.
The construction of dynamic models from normal stochastic operating
records or from the response to a PRBS can be based either on the concept
of cross-correlation or on the least-squares fit of a discrete equivalent model,
both topics in the field of system identification. They require substantial
presentation and background that are beyond the scope of this text. An introduction to system identification can be found in Chapter 8 of Franklin et al.
(1998), and a comprehensive treatment is given in Ljng (1999). Based
largely on the work of Professor Ljng, the MATLAB Toolbox on Identification provides substantial software to perform system identification and to
verify the quality of the proposed models.
W3.7.4 Obtaining a Pole-Zero Model from

Frequency-Response Data
As we pointed out earlier, it is relatively easy to obtain the frequencyresponse of a system experimentally. Sometimes it is desirable to obtain
an approximate model, in terms of a transfer function, directly from the
frequency response. The derivation of such a model can be done to various
degrees of accuracy. The method described in this section is usually adequate
and is widely used in practice.
There are two ways to obtain a model from frequency-response data.
In the first case, we can introduce a sinusoidal input, measure the gain
(logarithm of the amplitude ratio of output to input) and the phase difference
between output and input, and accept the curves plotted from this data as
the model. Using the methods given in previous sections, we can derive the
design directly from this information. In the second case, we wish to use the
frequency data to verify a mathematical model obtained by other means. To
do so we need to extract an approximate transfer function from the plots,
again by fitting straight lines to the data, estimating break points (that is,
finding the poles and zeros), and using Fig. 6.3 to estimate the damping
ratios of complex factors from the frequency overshoot. The next example
illustrates the second case.
EXAMPLE W3.6
Transfer Function from Measured Frequency Response

Determine a transfer function from the frequency response plotted in
Fig. W3.9, where frequency f is plotted in hertz.
Solution. Drawing an asymptote to the final slope of 2 (or 40 db per
decade), we assume a break point at the frequency where the phase is 90 .
This occurs at f1
= 1.66 Hz (1 = 2 f1 = 10.4 rad/sec). We need to know
Z08_FRAN6598_07_SE_W3.6.4 2014/5/13 9:14 page 52 #9

Figure W3.9
Experimental frequency
response
53
10
Magnitude with
second-order term
removed
Magnitude, |G|
f2
f0
Slope = -2
0.1
0.01
f1
0.1 0.2
2
f (Hz)
10
20
100
20
100
305
05
f2
Phase, jG
-305
f0
-605
Phase with
second-order
term removed
-905
-1205
-1505
-1805
0.1 0.2
2
f (Hz)
10
the damping ratio in order to subtract out this second-order pole. For this,
the phase curve may be of more help. Since the phase around the break-point
frequency is symmetric, we draw a line at the slope of the phase curve at f1
to find that the phase asymptote intersects the 0 line at f0
= 0.71 Hz (or
4.46 rad/sec). This corresponds to f1 /f0
= 2.34, which in time corresponds
to
= 0.5, as seen on the normalized response curves in Fig. 6.3b. The
magnitude curve with the second-order factor taken out shows an asymptotic
amplitude gain of about 6.0 db, or a factor of 106.0/20 = 2.0. As this is a
gain rise, it occurs because of a lead compensation of the form
s/a + 1
,
s/b + 1
where b/a = 2.0. If we remove the second-order terms in the phase curve,
we obtain a phase curve with a maximum phase of about 20 , which also
Z08_FRAN6598_07_SE_W3.6.4 2014/5/13 9:14 page 53 #10

corresponds to a frequency separation of about 2. To locate the center of the
lead compensation, we must estimate the point of maximum phase based
on the lead term alone, which occurs at the geometric mean of the two
break-point frequencies. The lead center seems to occur at f2
= 0.3 Hz (or
2 = 1.88 rad/sec).
Thus we have the relations
ab(1.88)2 = 3.55,
b
= 2,
a
from which we can solve
2a2 = 3.55,
a = 1.33,
b = 2.66.
Model from measured
response
Our final model is given by

(s/1.33) + 1
.
(W3.5)
[(s/2.66) + 1][(s/10.4)2 + (s/10.4) + 1]
The actual data were plotted from
(s/2) + 1
.
G(s) =
[(s/4) + 1][(s/10)2 + (s/10) + 1]
As can be seen, we found the second-order term quite easily, but the location
of the lead compensation is off in center frequency by a factor of 4/2.66
=
1.5. However, the subtraction of the second-order term from the composite
curve was not done with great accuracy, rather, by reading the curves. Again,
as with the transient response, we conclude that by a bit of approximate
plotting we can obtain a crude model (usually within a factor of 1.4 (3 db)
in amplitude and 10 in phase) that can be used for control design.
G(s)
=
Refinements on these techniques with computer aids are rather obvious,

and an interactive program for removing standard first- and second-order
terms and accurately plotting the residual function would greatly improve the
speed and accuracy of the process. It is also common to have computer tools
that can find the parameters of an assumed model structure by minimizing
the sum of squares of the difference between the models frequency response
and the experimental frequency response.
Further Reading for System Identification:
[1] L. Ljung, Perspectives on System Identification, http://users.isy.liu
.se/en/rt/ljung/seoul2dvinew/plenary2.pdf.
[2] L. Ljung, System Identification: Theory for the User, 2nd Ed.,
Prentice-Hall, 1999.
[3] G. F. Franklin, J. D. Powell, M. L. Workman, Digital Control of
Dynamic Systems, 3rd Ed. Ellis-Kagle Press, 1998.
Z08_FRAN6598_07_SE_W3.6.4 2014/5/13 9:14 page 54 #11
55
[4] M. B. Tischler and R. K. Remple, Aircraft and Rotorcraft System Identification: Engineering Methods with Flight-Test Examples, AIAA,
2006.
[5] R. Pintelon and J. Schoukens, System Identification: A Frequency
Domain Approach, 2nd ed., Wiley-IEEE Press, 2012.
[6] System Identification Toolbox, The Mathworks.
Z08_FRAN6598_07_SE_W3.6.4 2014/5/13 9:14 page 55 #12
Appendix W3.8
Amplitude and Time Scaling
The magnitude of the values of the variables in a problem is often very different, sometimes so much so that numerical difficulties arise. This was a
serious problem years ago when equations were solved using analog computers and it was routine to scale the variables so that all had similar magnitudes.
Todays widespread use of digital computers for solving differential equations has largely eliminated the need to scale a problem unless the number
of variables is very large, because computers are now capable of accurately
handling numbers with wide variations in magnitude. Nevertheless, it is wise
to understand the principle of scaling for the few cases in which extreme variations in magnitude exist and scaling is necessary or the computer word size
is limited.
W3.8.1 Amplitude Scaling

There are two types of scaling that are sometimes carried out: amplitude
scaling and time scaling, as we have already seen in Section 3.1.4. Amplitude scaling is usually performed unwittingly by simply picking units that
make sense for the problem at hand. For the ball levitator, expressing the
motion in millimeters and the current in milliamps would keep the numbers
within a range that is easy to work with. Equations of motion are sometimes
developed in the standard SI units such as meters, kilograms, and amperes,
but when computing the motion of a rocket going into orbit, using kilometers makes more sense. The equations of motion are usually solved using
computer-aided design software, which is often capable of working in any
units. For higher-order systems, it becomes important to scale the problem
so that system variables have similar numerical variations. A method for
accomplishing the best scaling for a complex system is first to estimate the
maximum values for each system variable and then to scale the system so
that each variable varies between 1 and 1.
In general, we can perform amplitude scaling by defining the scaled
variables for each state element. If
x = Sx x,
(W3.6)
then
x = Sx x
and
x = Sx x .
(W3.7)
We then pick Sx to result in the appropriate scale change, substitute

Eqs. (W3.6) and (W3.7) into the equations of motion, and recompute the
coefficients.
56
Z09_FRAN6598_07_SE_W3.6.5 2014/5/13 9:06 page 56 #1
Appendix W3.8 Amplitude and Time Scaling
EXAMPLE W3.7
57
Scaling for the Ball Levitator

The linearized equation of motion for the ball levitator (see Example 9.2,
Chapter 9) is
x = 1667x + 47.6i,
(W3.8)
where x is in units of meters and i is in units of amperes. Scale the variables
for the ball levitator to result in units of millimeters and milliamps instead
of meters and amps.
Solution. Referring to Eq. (W3.6), we define
x = Sx x
and
i = Si i,
such that both Sx and Si have a value of 1000 in order to convert x and i
in meters and amps to x and i in millimeters and milliamps. Substituting
these relations into Eq. (W3.8) and taking note of Eq. (W3.7) yields
x = 1667x + 47.6
Sx
i .
Si
In this case, Sx = Si , so Eq. (W3.8) remains unchanged. Had we scaled the

two quantities by different amounts, there would have been a change in the
last coefficient in the equation.
W3.8.2 Time Scaling

The unit of time when using SI units or English units is seconds. Computeraided design software is usually able to compute results accurately no matter
how fast or slow the particular problem at hand. However, if a dynamic system responds in a few microseconds, or if there are characteristic frequencies
in the system on the order of several megahertz, the problem may become
ill conditioned, so that the numerical routines produce errors. This can be
particularly troublesome for high-order systems. The same holds true for
an extremely slow system. It is therefore useful to know how to change the
units of time should you encounter an ill-conditioned problem.
We define the new scaled time to be
= o t
(W3.9)
such that, if t is measured in seconds and o = 1000, then will be measured

in milliseconds. The effect of the time scaling is to change the differentiation
so that
dx
dx
dx
x =
=
= o
(W3.10)
dt
d(/o )
d
and
x =
d2x
d2x
= o2 2 .
2
dt
d
Z09_FRAN6598_07_SE_W3.6.5 2014/5/13 9:06 page 57 #2
(W3.11)
58 Appendix W3.8 Amplitude and Time Scaling
EXAMPLE W3.8
Time Scaling an Oscillator

The equation for an oscillator was derived in Example 2.5. For a case with
a very fast natural frequency n = 15,000 rad/sec (about 2 kHz), Eq. (2.23)
can be rewritten as
+ 15,0002 = 106 Tc .
Determine the time-scaled equation so that the unit of time is milliseconds.
Solution. The value of o in Eq. (W3.9) is 1000. Equation (W3.11) shows
that
d2
= 106 ,
d 2
and the time-scaled equation becomes
d2
+ 152 = Tc .
d 2
In practice, we would then solve the equation
+ 152 = Tc
(W3.12)
and label the plots in milliseconds instead of seconds.
W3.8.3 Time and Amplitude Scaling in State-Space

We have already discussed time and amplitude scaling in Chapter 3. We
now extend the ideas to the state-variable form. Time scaling with = o t
replaces Eq. (7.3) with
dx
1
1
+ Bu.
=
Ax +
Bu = Ax
d
o
o
(W3.13)
Amplitude scaling of the state corresponds to replacing x with z = D1

x x,
where Dx is a diagonal matrix of scale factors. Input scaling corresponds to
replacing u with v = D1
u u. With these substitutions,
Dx z =
1
1
ADx z +
BDu v.
o
o
(W3.14)
Then
z =
1 1
1 1
+ Bv.
D ADx z +
D BDu v = Az
o x
o x
(W3.15)
Equation (W3.15) compactly expresses the time- and amplitude-scaling

operations. Regrettably, it does not relieve the engineer of the responsibility of actually thinking of good scale factors so that scaled equations are
in good shape.
Z09_FRAN6598_07_SE_W3.6.5 2014/5/13 9:06 page 58 #3
Appendix W3.8 Amplitude and Time Scaling
EXAMPLE W3.9
59
Time Scaling an Oscillator

The equation for an oscillator was derived in Example 2.5. For a case with
a very fast natural frequency n = 15,000 rad/sec (about 2 kHz), Eq. (2.23)
can be rewritten as
+ 15,0002 = 106 Tc .
Determine the time-scaled equation so that the unit of time is milliseconds.
Solution. In state-variable form with a state vector x = [ ]T , the unscaled
matrices are

0
0
1
and
B
=
.
A=
106
15,0002 0
Applying Eq. (W3.13) results in

1
0
1000
=
A
2
15,000
0
1000
and
B =
0
103
which yields state-variable equations that are scaled.
Z09_FRAN6598_07_SE_W3.6.5 2014/5/13 9:06 page 59 #4
Appendix W4.1.4.1
The Filtered Case
Thus far the analysis has been based on the simplest open- and closed-loop
structures. A more general case includes a dynamic filter on the input and
also dynamics in the sensor. The filtered open-loop structure is shown in
Fig. W4.1 having the transfer function Tol = GDol F. In this case, the openloop controller transfer function has been simply replaced by DF and the
discussion given for the unfiltered open-loop case is easily applied to this
change.
For the filtered feedback case shown in Fig. W4.2, the changes are more
significant. In that case, the transform of the system output is given by
Y=
G
HGDcl
GDcl F
R+
W
V.
1 + GDcl H
1 + GDcl H
1 + GDcl H
(W4.1)
As is evident from this equation, the sensor dynamics, H, is part of the loop
transfer function and enters into the question of stability with Dcl H replacing
the Dcl of the unity feedback case. In fact, if F = H, then with respect to
stability, tracking, and regulation, the filtered case is identical to the unity
case with Dcl H replacing Dcl . On the other hand, the filter transfer function F
can play the role of the open-loop controller except that here the filter F would
GDcl
be called on to modify the entire loop transfer function, 1+GD
, rather than
cl H
simply GDol . Therefore the filtered closed-loop structure can realize the best
W(s)
R(s)
Filter
F(s)
Rf (s)
Controller
Dol (s)
U(s)
+
Plant
G(s)
Y(s)
Figure W4.1
Filtered open-loop system
W(s)
R(s)
Filter
F(s)
Rf (s) +
Controller
Dcl (s)
U(s)
+
Plant
G(s)
Sensor
H(s)
Y(s)
+
V(s)
+
Figure W4.2
Filtered closed-loop. R = reference, U = control, Y = output, V = sensor noise
60
Z10_FRAN6598_07_SE_W4.1.4.1 2014/5/7 2:13 page 60 #1
Appendix W4.1.4.1 The Filtered Case
61
properties of both the open-loop and the unity feedback closed-loop cases.
The controller, Dcl , can be designed to effectively regulate the system for
the disturbance W and the sensor noise V , while the filter F is designed to
improve the tracking accuracy. If the sensor dynamics H are accessible to
the designer, this term also can be designed to improve the response to the
sensor noise. The remaining issue is sensitivity.
Using the formula given in Eq. (4.18), with changes in the parameter of
interest, we can compute
SFTcl = 1.0,
(W4.2)
1
,
(W4.3)
1 + GDcl H
GDcl H
.
(W4.4)
SHTcl =
1 + GDcl H
Of these, the most interesting is the last. Notice that with respect to H, the
sensitivity approaches unity as the loop gain grows. Therefore it is particularly important that the transfer function of the sensor be not only low in
noise but also very stable in gain. Money spent on the sensor is money well
spent!
SGTcl =
EXAMPLE W4.1
If S is the sensitivity of the filtered feedback system to changes in the plant

transfer function and T is the transfer function from reference to output,
compute the sum of S + T . Show that S + T = 1 if F = H.
(a) Compute the sensitivity of the filtered feedback system shown in
Fig. W4.2 with respect to changes in the plant transfer function, G.
(b) Compute the sensitivity of the filtered feedback system shown in
Fig. W4.2 with respect to changes in the controller transfer function, Dcl .
(c) Compute the sensitivity of the filtered feedback system shown in Fig.
W4.2 with respect to changes in the filter transfer function, F.
(d) Compute the sensitivity of the filtered feedback system shown in Fig.
W4.2 with respect to changes in the sensor transfer function, H. If S is
the sensitivity of the filtered feedback system to changes in the plant
transfer function and T is the transfer function from reference to output,
compute the sum of S + T . Show that S + T = 1 if F = H.
Solution. To answer the first question, we need the answer to part (a), so
lets start there.
(a) Applying the formula for sensitivity of T to changes in G:
FDG
,
1 + DGH
1 + DGH (1 + DGH)FD FDG(DH)
S=G
FDG
(1 + DGH)2
1
=
.
1 + DGH
T=
Z10_FRAN6598_07_SE_W4.1.4.1 2014/5/7 2:13 page 61 #2
62 Appendix W4.1.4.1 The Filtered Case

Now we can do
1
FDG
+
1 + DGH
1 + DGH
1 + FDG
=
1 + DGH
= 1 if F = H.
S+T =
(W4.5)
(b) Applying the formula for sensitivity of T to changes in D:

1 + DGH (1 + DGH)FG FDG(GH)
FDG
(1 + DGH)2
1
=
.
1 + DGH
This is not surprising as D and G are in series.
(c) Applying the formula for sensitivity of T to changes in F:
STD = D
1 + DGH (1 + DGH)(DG)
FDG
(1 + DGH)2
1 + DGH
1 + DGH
= 1.
STF = F
In this case, the F term is in the open loop so it has a sensitivity of

unity.
(d) Applying the formula for sensitivity of T to changes in H:
1 + DGH (1 + DGH)0 FDG(DG)
FDG
(1 + DGH)2
DGH
=
.
(1 + DGH)
STH = H
Z10_FRAN6598_07_SE_W4.1.4.1 2014/5/7 2:13 page 62 #3
Appendix W4.2.2.1
Truxals Formula for the Error
Constants
Truxal (1955) derived a formula for the velocity constant of a Type 1 system
in terms of the closed-loop poles and zeros, which is a formula that connects
the steady-state error to the systems dynamic response. Since control design
often requires a trade-off between these two characteristics, Truxals formula
can be useful to know. Its derivation is quite direct. Suppose the closed-loop
transfer function T (s) of a Type 1 system is
(s z1 )(s z2 ) (s zm )
.
(W4.6)
T (s) = K
(s p1 )(s p2 ) (s pn )
Since the steady-state error in response to a step input in a Type 1 system is
zero, the DC gain is unity; thus
T (0) = 1.
The system error is given by
(W4.7)

Y (s)
= R(s)[1 T (s)].
(W4.8)
E(s) R(s) Y (s) = R(s) 1
R(s)
The system error due to a unit ramp input is given by
1 T (s)
.
(W4.9)
E(s) =
s2
Using the Final Value Theorem, we get
1 T (s)
.
(W4.10)
ess = lim
s0
s
Using LHpitals rule, we rewrite Eq. (W4.10) as
dT
,
(W4.11)
ess = lim
s0 ds
or
dT
1
.
(W4.12)
ess = lim
=
s0 ds
Kv
Equation (W4.12) implies that 1/Kv is related to the slope of the transfer
function at the origin, which is a result that will also be shown in Section
6.1.2. Using Eq. (W4.7), we can rewrite Eq. (W4.12) as
dT 1
ess = lim
,
(W4.13)
s0 ds T
or
d
(W4.14)
ess = lim [ln T (s)].
s0 ds
63
Z11_FRAN6598_07_SE_W4.2.2.1 2014/5/7 2:17 page 63 #1
64 Appendix W4.2.2.1 Truxals Formula for the Error Constants

Substituting Eq. (W4.6) into Eq. (W4.14), we get
m

d
(s zi )
ess = lim
ln K ni=1
,
s0 ds
i=1 (s pi )

m
m

d
= lim
ln K +
ln(s zi )
ln(s pi ) ,
s0 ds
i=1
or
(W4.15)
(W4.16)
i=1
n
m

d ln T
1 1
1
=
=
+
.
Kv
ds
s=0
pi
zi
i=1
(W4.17)
i=1
We observe from Eq. (W4.17) that Kv increases as the closed-loop

poles move away from the origin. Similar relationships exist for other error
coefficients, and these are explored in the problems.
EXAMPLE W4.2
Truxals Formula
Truxals formula
A third-order Type 1 system has closed-loop poles at 2 2j and 0.1.

The system has only one closed-loop zero. Where should the zero be if a
Kv = 10 sec1 is desired?
Solution. From Truxals formula, we have
1
1
1
1
1
=
+ ,
Kv
2 + 2j 2 2j 0.1
z
or
1
0.1 = 0.5 + 10 + ,
z
1
= 0.1 0.5 10
z
= 10.4.
Therefore, the closed-loop zero should be at z = 1/ 10.4 = 0.096.
Z11_FRAN6598_07_SE_W4.2.2.1 2014/5/7 2:17 page 64 #2
Appendix W4.5
Introduction to Digital
Control
As a result of the revolution in the cost-effectiveness of digital computers,
there has been an increasing use of digital logic in embedded applications such as controllers in feedback systems. A digital controller gives the
designer much more flexibility to make modifications to the control law after
the hardware design is fixed, because the formula for calculating the control
signal is in the software rather than the hardware. In many instances, this
means that the hardware and software designs can proceed almost independently, saving a great deal of time. Also, it is relatively easy to include binary
logic and nonlinear operations as part of the function of a digital controller as
compared to an analog controller. Special processors designed for real-time
signal processing and known as digital signal processors (DSPs) are particularly well suited for use as real-time controllers. Chapter 8 includes a more
extensive introduction to the math and concepts associated with the analysis
and design of digital controllers and digital control systems. However, in
order to be able to compare the analog designs of the next three chapters
with reasonable digital equivalents, we give here a brief introduction to the
most simple techniques for digital designs.
A digital controller differs from an analog controller in that the signals
must be sampled and quantized.1 A signal to be used in digital logic needs
to be sampled first and then the samples need to be converted by an analog-todigital converter or A/D2 into a quantized digital number. Once the digital
computer has calculated the proper next control signal value, this value
needs to be converted back into a voltage and held constant or otherwise
extrapolated by a digital-to-analog converter or D/A3 in order to be applied
to the actuator of the process. The control signal is not changed until the
next sampling period. As a result of sampling, there are strict limits on the
speed and bandwidth of a digital controller. Discrete design methods that
tend to minimize these limitations are described in Chapter 8. A reasonable
rule of thumb for selecting the sampling period is that, during the risetime of the response to a step, the input to the discrete controller should be
sampled approximately six times. By adjusting the controller for the effects
1A controller that operates on signals that are sampled but not quantized is called discrete,
while one that operates on signals that are both sampled and quantized is called digital.
2 Pronounced A to D.
3 Often spelled DAC and pronounced as one word to rhyme with quack.
65
Z12_FRAN6598_07_SE_W4.5 2014/5/13 9:10 page 65 #1
66 Appendix W4.5 Introduction to Digital Control

of sampling, the sample period can be as large as two to three times per rise
time. This corresponds to a sampling frequency that is 10 to 20 times the
systems closed-loop bandwidth. The quantization of the controller signals
introduces an equivalent extra noise into the system; to keep this interference
at an acceptable level, the A/D converter usually has an accuracy of 10 to 12
bits although inexpensive systems have been designed with only 8 bits. For
a first analysis, the effects of the quantization are usually ignored, as they
will be in this introduction. A simplified block diagram of a system with a
digital controller is shown in Fig. W4.3.
For this introduction to digital control, we will describe a simplified
technique for finding a discrete (sampled but not quantized) equivalent to a
given continuous controller. The method depends on the sampling period, Ts ,
being short enough that the reconstructed control signal is close to the signal
that the original analog controller would have produced. We also assume that
the numbers used in the digital logic have enough accurate bits so that the
quantization implied in the A/D and D/A processes can be ignored. While
there are good analysis tools to determine how well these requirements are
met, here we will test our results by simulation, following the well-known
advice that The proof of the pudding is in the eating.
Finding a discrete equivalent to a given analog controller is equivalent
to finding a recurrence equation for the samples of the control, which will
approximate the differential equation of the controller. The assumption is
that we have the transfer function of an analog controller and wish to replace
it with a discrete controller that will accept samples of the controller input
e(kTs ) from a sampler and, using past values of the control signal u(kTs ) and
present and past samples of the input e(kTs ), will compute the next control
signal to be sent to the actuator. As an example, consider a PID controller
with the transfer function
U(s) = (kP +
kI
+ kD s)E(s),
s
(W4.18)
which is equivalent to the three terms of the time-domain expression

t
u(t) = kP e(t) + kI
e( )d + kD e (t),
(W4.19)
0
= uP + uI + uD .
A/D
e(kT)
Digital controller
D(z)
u(kT)
D/A
Plant
G
(W4.20)
Clock
Sensor
H
Figure W4.3
Block diagram of a digital controller
Z12_FRAN6598_07_SE_W4.5 2014/5/13 9:10 page 66 #2
Appendix W4.5 Introduction to Digital Control
67
Based on these terms and the fact that the system is linear, the next control
sample can be computed term-by-term. The proportional term is immediate:
uP (kTs + Ts ) = kP e(kTs + Ts ).
(W4.21)
The integral term can be computed by breaking the integral into two parts
and approximating the second part, which is the integral over one sample
period, as follows.
kTs +Ts
uI (kTs + Ts ) = kI
e( )d ,
0
kTs
= kI
(W4.22)
e( )d + kI
kTs +Ts
e( )d ,
(W4.23)
kTs
= uI (kTs ) + {area under e( ) over one period},

Ts
= uI (kTs ) + kI {e(kTs + Ts ) + e(kTs )}.

2
(W4.24)
(W4.25)
In Eq. (W4.25), the area in question has been approximated by that of the
trapezoid formed by the base Ts and vertices e(kTs +Ts ) and e(kTs ) as shown
by the dashed line in Fig. W4.4.
The area also can be approximated by the rectangle of amplitude e(kTs )
and width Ts shown by the solid blue in Fig. W4.4 to give uI (kTs +
Ts ) = uI (kTs ) + kI Ts e(kTs ). These and other possibilities are considered
in Chapter 8.
In the derivative term, the roles of u and e are reversed from integration
and the consistent approximation can be written down at once from Eqs.
(W4.25) and (W4.19) as
Ts
{uD (kTs + Ts ) + uD (kTs )} = kD {e(kTs + Ts ) e(kTs )}.
2
(W4.26)
As with linear analog transfer functions, these relations are greatly simplified
and generalized by the use of transform ideas. At this time, the discrete
transform will be introduced simply as a prediction operator z much as if we
.
x
Figure W4.4
Graphical
interpretation of
numerical integration
.
x = f (x, u)
.
x (ti)
t.
x dt
0
ti
ti + 1
Z12_FRAN6598_07_SE_W4.5 2014/5/13 9:10 page 67 #3

described the Laplace transform variable, s, as a differential operator.4 Here
we define the operator z as the forward shift operator in the sense that if U(z)
is the transform of u(kTs ) then zU(z) will be the transform of u(kTs + Ts ).
With this definition, the integral term can be written as
zUI (z) = UI (z) + kI
Ts
[zE(z) + E(z)] ,
2
(W4.27)
Ts z + 1
E(z),
(W4.28)
2 z1
and from Eq. (W4.26), the derivative term becomes the inverse as
UI (z) = kI
UD (z) = kD
2 z1
E(z).
Ts z + 1
The complete discrete PID controller is thus described by

Ts z + 1
2 z1
U(z) = kP + kI
+ kD
E(z).
2 z1
Ts z + 1
(W4.29)
(W4.30)
Comparing the two discrete equivalents of integration and differentiation

with the corresponding analog terms, it is seen that the effect of the discrete
approximation in the z domain is as if everywhere in the analog transfer
function the operator s has been replaced by the composite operator T2s z1
z+1 .
This is the trapezoid rule5 of discrete equivalents.
The discrete equivalent to Dc (s) is

2 z 1
Dd (z) = Dc
.
(W4.31)
Ts z + 1
Trapezoid rule
EXAMPLE W4.3
Discrete Equivalent
Find the discrete equivalent to the analog controller having transfer function
Dc (s) =
U(s)
11s + 1
=
,
E(s)
3s + 1
(W4.32)
using the sample period Ts = 1.

Solution. The discrete operator is 2(z1)
z+1 and thus the discrete transfer
function is
U(z)
Dd (z) =
= Dc (s)|
(W4.33)
z1
E(z)
s= T2s
,
z+1

2(z1)
11 z+1 + 1

=
.
(W4.34)
3 2(z1)
+
1
z+1
4 This is defined as the z-transform in Chapter 8.
5 The formula is also called Tustins method after the English engineer who used the technique
to study the responses of nonlinear circuits.
Z12_FRAN6598_07_SE_W4.5 2014/5/13 9:10 page 68 #4
69
Clearing fractions, the discrete transfer function is

Dd (z) =
U(z)
23z 21
=
.
E(z)
7z 5
(W4.35)
Converting the discrete transfer function to a discrete difference equation

using the definition of z as the forward shift operator is done as follows.
First we cross-multiply in Eq. (W4.35) to obtain
(7z 5)U(z) = (23z 21)E(z),
(W4.36)
and interpreting z as a shift operator, this is equivalent to the difference

equation6
7u(k + 1) 5u(k) = 23e(k + 1) 21e(k),
(W4.37)
where we have replaced kTs + Ts with k + 1 to simplify the notation. To

compute the next control at time kTs + Ts , therefore, we solve the difference
equation
u(k + 1) =
5
23
21
u(k) + e(k + 1) e(k).
7
7
7
(W4.38)
Now lets apply these results to a control problem. Fortunately, Matlab

provides us with the Simulink capability to simulate both continuous and
discrete systems allowing us to compare the responses of the systems with
continuous and discrete controllers.
EXAMPLE W4.4
Equivalent Discrete Controller for Speed Control

A motor speed control is found to have the plant transfer function
Y
45
=
.
U
(s + 9)(s + 5)
(W4.39)
A PI controller designed for this system has the transfer function

Dc (s) =
s+6
U
= 1.4
.
E
s
(W4.40)
The closed-loop system has a rise time of about 0.2 sec and an overshoot
of about 20%. Design a discrete equivalent to this controller and compare
the step responses and control signals of the two systems. (a) Compare the
responses if the sample period is 0.07, which is about three samples per rise
time. (b) Compare the responses with a sample period of Ts = 0.035 sec,
which corresponds to about six samples per rise time.
6 The process is entirely similar to that used in Chapter 3 to find the ordinary differential equation
to which a rational Laplace transform corresponds.
Z12_FRAN6598_07_SE_W4.5 2014/5/13 9:10 page 69 #5

Solution. (a) Using the substitution given by Eq. (W4.31), the discrete
2 z1
equivalent for Ts = 0.07 sec is given by replacing s by s 0.07
z+1 in
Dc (s) as
2 z1
+6
0.07
z+1
Dd (z) = 1.4
,
(W4.41)
2 z1
0.07 z + 1
2(z 1) + 6 0.07(z + 1)
= 1.4
,
(W4.42)
2(z 1)
1.21z 0.79
= 1.4
.
(W4.43)
(z 1)
Based on this expression, the equation for the control is (the sample period
is suppressed)
u(k + 1) = u(k) + 1.4 [1.21e(k + 1) 0.79e(k)].
(b) For Ts = 0.035 sec, the discrete transfer function is
1.105z 0.895
Dd = 1.4
,
z1
for which the difference equation is
(W4.44)
(W4.45)
u(k + 1) = u(k) + 1.4[1.105 e(k + 1) 0.895 e(k)].

A Simulink block diagram for simulating the two systems is given in
Fig. W4.5, and plots of the step responses are given in Fig. W4.6a. The
respective control signals are plotted in Fig. W4.6b. Notice that the discrete
controller for Ts = 0.07 sec results in a substantial increase in the overshoot
in the step response, while with Ts = 0.035 sec the digital controller matches
the performance of the analog controller fairly well.
For controllers with many poles and zeros, making the continuousto-discrete substitution called for in Eq. (W4.31) can be very tedious.
PI Control s+ 6
s
Slider Kc
1.4
Step
Tau 1
9
s+9
Mux
Discrete 1.21z - 0.79

z-1
PI Control
Slider Kd
1.4
Tau 1
9
s+9
Tau 2
5
s+5
Control
Mux 1
Tau 2
5
s+5
Figure W4.5
Simulink block diagram to compare continuous and discrete controllers
Z12_FRAN6598_07_SE_W4.5 2014/5/13 9:10 page 70 #6
Output

1.4
2.5
Digital controller (Ts = 0.07 sec)
1.2
Continuous controller
2.0
Digital controller (Ts = 0.07 sec)
1.0
1.5
0.8
Discrete controller (Ts = 0.035 sec)
0.6
1.0
Discrete controller (Ts = 0.035 sec)
0.4
0.5
0.2
0
71
0.5
1.5
2 2.5 3
Time (sec)
(a)
3.5
4.5
0.5
1.5
2.5 3
Time (sec)
3.5
4.5
(b)
Figure W4.6
Comparison plots of a speed control system with continuous and discrete controllers: (a) output
responses and (b) control signals
Fortunately, Matlab provides a command that does all the work. If one has a
continuous transfer function given by Dc (s) = numD
denD represented in Matlab
as sysDc = tf(numD,denD), then the discrete equivalent with sampling period
Ts is given by
sysDd = c2d(sysDc,Ts ,' t' ).
(W4.46)
In this expression, of course, the polynomials are represented in Matlab
form. The last parameter in the c2d function given by t calls for the conversion to be done using the trapezoid (or Tustin) method. The alternatives
can be found by asking Matlab for help c2d. For example, to compute the
polynomials for Ts = 0.07 sec for Example W4.4, the commands would be
numD = [1 6];
denD = [1 - 0];
sysDc = tf(numD,denD)
sysDd = c2d( sysDc,0.07,' t' )
Z12_FRAN6598_07_SE_W4.5 2014/5/13 9:10 page 71 #7
Appendix W4.6
Sensitivity of Time Response
to Parameter Change
We have considered the effects of errors on the steady-state gain of a dynamic
system and showed how feedback control can reduce these errors. Since
many control specifications are in terms of the step response, the sensitivity of
the time response to parameter changes is sometimes very useful to explore.
For example, by looking at the sensitivity plot, we can tell if increasing a
particular parameter will increase or decrease the overshoot of the response.1
The following analysis is also a good exercise in small-signal linearization.
To consider the sensitivity of the output y(t, ) of a system having
a parameter of interest, , we compute the effect of a perturbation in
the parameter, , on the nominal response by using the Taylors series
expansion
y(t, + ) = y(t, ) +
y
+ .
(W4.47)
The first-order approximation of the parameter perturbation effect is the term

y(t) =
y
.
(W4.48)
This function can be generated from the system itself as shown by

Perkins et al., 1991. We assume that the response depends linearly on the
parameter and therefore that the overall transfer function T (s, ) is composed
of component transfer functions that can be defined to bring out the dependence on the parameter explicitly. A block diagram of the transfer function
in terms of the components Tij (s) can be expressed as shown in Fig. W4.7,
where we have labeled the parameter as and its input signal as Z. In terms
Figure W4.7
Block diagram showing
the dependence of
output Y on
parameter
T11
+
X
Y
+
T22
T12
T21
1As we will see in Chapter 5, the development of the MATLABTM root locus interface rltool
gives the designer a computer aid to this result.
72
Z13_FRAN6598_07_SE_W4.6 2014/5/8 3:23 page 72 #1
Appendix W4.6 Sensitivity of Time Response to Parameter Change
73
of this block diagram, the equations relating Y and Z to the reference input
can be written immediately.
Y = T11 R + T21 Z,
(W4.49)
Z = T12 R + T22 Z.
(W4.50)
The perturbed equations are

Y + Y = T11 R + T21 ( + )(Z + Z),
(W4.51)
Z + Z = T12 R + T22 ( + )(Z + Z).
(W4.52)
Multiplying these out and ignoring the small term Z, the expressions for
the perturbations in Y and Z are given by
Y = T21 (Z + Z),
(W4.53)
Z = T22 (Z + Z).
(W4.54)
The solutions to these equations can be best presented as a block diagram,

y
, and we notice
shown in Fig. W4.8a. The output of this figure is Y =
that the input Z is multiplied by a gain of . Therefore, if we drop the

y
block , the output will be simply
as shown in Fig. W4.8b. Finally, to
compute the sensitivity as the variation to a percent change in the parameter,

Figure W4.8
Block diagrams showing
the generation of (a) Y
y
and Z ; (b)
; and (c)
du
T21
0y
dY = 0u du
T22
dZ
u
(a)
0y
0u
T21
Z
T22
u
(b)
u
+
T21
0y
u 0u
T22
(c)
Z13_FRAN6598_07_SE_W4.6 2014/5/8 3:23 page 73 #2
74 Appendix W4.6 Sensitivity of Time Response to Parameter Change

Figure W4.9
the computation of
y
from the original
transfer function
0y
u 0u
Y
+
u
u
+
y(t, )
y
ln = , we need only shift the input Z from
the output side of the block to its input as shown in Fig. W4.8c. We are
now in a position to give the final block diagram of the system as it is to be
implemented, shown in Fig. W4.9.
In this figure, it is clear that, to compute the sensitivity of the output to
a parameter, one needs to simulate two copies of the system. The input to
the first system is the reference input of interest, and the input to the second
system is at the input to the parameter of interest of the variable Z taken
from the input to the parameter in the original system. The transfer function
from the reference input to the output sensitivity is readily computed to be
y
=
which is
ln
T12 T21
.
(1 T22 )2
(W4.55)
From this function it is clear that:

To keep the sensitivity of the output signal to a parameter change
low, it is important to have feedback with high gain around the
parameter in question.
Response sensitivity
EXAMPLE W4.5
Time-Domain Sensitivity
Compute the sensitivity of the output of the speed control example shown in
the upper portion of Fig. W4.10 with respect to the control gain, Kcl . Take
the nominal values to be Kcl = 9, = 0.01 sec, and A = 1 rad/volt-sec.
Figure W4.10
the computation of the
sensitivity of the output
of the speed control
example
Kcl
+
A
ts + 1
Kcl
A
ts + 1
0Y
0K
Z13_FRAN6598_07_SE_W4.6 2014/5/8 3:23 page 74 #3
Appendix W4.6 Sensitivity of Time Response to Parameter Change
75
Solution. The required block diagram for the computation is given in Fig.
W4.10 based on Fig. W4.9. In MATLAB, we will construct the several transnij
fer functions with Tij =
and implement Eq. (W4.55). For comparison,
dij
we compute the nominal response from Fig. W4.7 and add 10% of the sensitivity to the nominal response. The instructions to do the computation in
MATLAB are
% script to compute sensitivity for Fig. W4.10
% First input the data for the component transfer functions Tij
% and the nominal parameter, Kcl for this problem
Kcl=9; tau=0.01;
n11=0; d11=1;
n12=1; d12=1;
n22=[0 -1]; d22= [tau 1];
n21=1; d21=[tau 1];
% Now compute the numerator and denominator polynomials of the
transfer functions
% using the convolution function conv to multiply the polynomials
% and put them into system transfer function forms with the MATLAB
function tf.
% The over-all transfer function is
% Y/R =n11/d11 + (n12*n21*d22)/(d12*d21*[d22-Kcl*n22]) = sysy
% The transfer function from the reference input to the sensitivity is
% Kcl*(dy/dKcl)/R = sysdy
% Now dene the numerators and denominators of several intermediate
transfer functions
n1=Kcl*conv(n21,n12);
d1=conv(d21,d12);
n2=d22;
d2=[d22-Kcl*n22];
ny=conv(n1,n2);
dy=conv(d1,d2);
% Now put these together to form two intermediate transfer functions
sysy1 = tf(ny,dy);
sysy2 = tf(n11,d11);
% Now construct the nal transfer functions
% The overall transfer function Y/R
sysy=sysy1+sysy2;
% The sensitivity transfer function
ndy=conv(ny,n2);
ddy=conv(dy,d2);
sysdy=tf(ndy,ddy);
% Now use these to compute the step responses and
% plot the output, the sensitivity and a perturbed response
[y,t]=step(sysy);
[yd,t]=step(sysdy);
plot(t,[y yd y+.1*yd]);
Z13_FRAN6598_07_SE_W4.6 2014/5/8 3:23 page 75 #4
76 Appendix W4.6 Sensitivity of Time Response to Parameter Change

Figure W4.11
Plots of the output, the
sensitivity, and the
result of a 10% change
in the parameter value
for the speed control
example
1.0
y + dy
0.8
y
0.6
0.4
Kcl
dy
dKcl
0.2
3
4
Time (m sec)
These instructions are constructed to compute the sensitivity for any

system, given the several transfer functions. The script input is for the specific
example. Plots of the output, its sensitivity, and the result of a 10% change
in the parameter value are given in Fig. W4.11.
Z13_FRAN6598_07_SE_W4.6 2014/5/8 3:23 page 76 #5
Appendix W5.4.4
Analog and Digital
Implementations
Lead compensation can be physically realized in many ways. In analog
electronics, a common method is to use an operational amplifier, an example
of which is shown in Fig. W5.1. The transfer function of the circuit in
Fig. W5.1 is readily found by the methods of Chapter 2 to be
s+z
Dlead (s) = a
,
(W5.1)
s+p
where
p
a = , if Rf = R1 + R2 ,
z
1
z=
,
R1 C
R1 + R2
1
p=
.
R2
R1 C
If a design for Dc (s) is complete and a digital implementation is desired,
then the technique ofAppendix W4.5 can be used by first selecting a sampling
period Ts and then making substitution of T2s z1
z+1 for s. For example, consider
s+2
the lead compensation Dc (s) = s+13 . Then, since the rise time is about
0.3, a sampling period of six samples per rise time results in the selection
2 z1
of Ts = 0.05 sec. With the substitution of 0.05
z+1 for s into this transfer
function, the discrete transfer function is
40 z1
U(z)
z+1 + 2
=
z1
E(z)
40 z+1 + 13
1.55z 1.4
.
(W5.2)
1.96z 1
Clearing fractions and using the fact that operating on the time functions
zu(kTs ) = u(kTs + Ts ), we see that Eq. (W5.2) is equivalent to the formula
for the controller given by
1
1.55
1.4
u(kTs ) +
e(kTs + Ts )
e(kTs ). (W5.3)
u(kTs + Ts ) =
1.96
1.96
1.96
=
Figure W5.1
Possible circuit of a lead
compensation
Vin
R1
Rf
R2
V0
77
Z14_FRAN6598_07_SE_W5.4.4 2014/5/13 9:18 page 77 #1
78 Appendix W5.4.4 Analog and Digital Implementations

The Matlab commands to generate the discrete equivalent controller are
sysC=tf([1 2],[1 13]);
sysD=c2D(sysC,0.05,tustin)
Fig. W5.2 shows the Simulink diagram for implementing the digital
controller. The result of the simulation is contained in Fig. W5.3, which
shows the comparison of analog and digital control outputs, and Fig. W5.4,
which shows the analog and digital control outputs.
As with lead compensation, lag or notch compensation can be implemented using a digital computer and following the same procedure. However,
Step
s+2
s + 13
91
Pl
control
Slider
Kc
1
s
1
s+1
Tau 1
Tau 2
Mux
Mux 1
Mux
Mux 1
1.55z - 1.4
1.96z + 1
91
Discrete
Pl control
Slider
Kd
Output
Control
1
s
1
s+1
tau 1
tau 2
Figure W5.2
Simulink diagram for comparison of analog and digital control
Figure W5.3
Comparison of analog
and digital control
output responses
1.4
Digital controller
1.2
Amplitude
1
0.8
0.6
0.4
0.2
0
0.5
1
Time (sec)
1.5
Z14_FRAN6598_07_SE_W5.4.4 2014/5/13 9:18 page 78 #2
Appendix W5.4.4 Analog and Digital Implementations

Figure W5.4
Comparison of analog
and digital control time
histories
79
100
80
Amplitude
60
40
Digital controller
20
0
-20
-40
0.5
1
Time (sec)
1.5
Figure W5.5
Possible circuit of lag
compensation
R1
Ri
R2
they too can be implemented using analog electronics, and a circuit diagram
of a lag network is given in Fig. W5.5. The transfer function of this circuit
can be shown to be
s+z
Dc (s) = a
,
s+p
where
R2
a=
,
Ri
R1 + R2
z=
,
R1 R 2 C
1
p=
.
R1 C
Usually Ri = R2 , so the high-frequency gain is unity, or a = 1, and the
low-frequency increase in gain to enhance Kv or other error constant is set
2
by k = a pz = R1R+R
.
2
Z14_FRAN6598_07_SE_W5.4.4 2014/5/13 9:18 page 79 #3
Appendix W5.6.3
Root Locus with Time Delay
Time delays always reduce

the stability of a system
An example of a root locus

with time delay
Time delays often arise in control systems, both from delays in the process
itself and from delays in the processing of sensed signals. Chemical plants
often have processes with a time delay representing the time material takes
to be transported via pipes or other conveyer. In measuring the attitude of a
spacecraft en route to Mars, there is a significant time delay for the sensed
quantity to arrive back on Earth due to the speed of light. There is also a
small time delay in any digital control system due to the cycle time of the
computer and the fact that data is processed at discrete intervals. Time delay
always reduces the stability of a system; therefore, it is important to be able
to analyze its effect. In this section, we discuss how to use the root locus for
such analysis. Although an exact method of analyzing time delay is available
in the frequency-response methods described in Chapter 6, knowing several
different ways to analyze a design provides the control designer with more
flexibility and an ability to check the candidate solutions.
Consider the problem of designing a control system for the temperature
of the heat exchanger described in Chapter 2. The transfer function between
the control As and the measured output temperature Tm is described by two
first-order terms plus a time delay Td of 5 sec. The time delay results because
the temperature sensor is physically located downstream from the exchanger,
so that there is a delay in its reading. The transfer function is
G(s) =
e5s
,
(10s + 1)(60s + 1)
(W5.4)
where the e5s term arises from the time delay.1

The corresponding root-locus equations with respect to proportional
gain K are
1 + KG(s) = 0,
1+K
e5s
(10s + 1)(60s + 1)
= 0,
600s2 + 70s + 1 + Ke5s = 0.
(W5.5)
How would we plot the root locus corresponding to Eq. (W5.5)? Since it
is not a polynomial, we cannot proceed with the methods used in previous examples. So we reduce the given problem to one we have previously
1 Time delay is often referred to as transportation lag in the process industries.
80
Z15_FRAN6598_07_SE_W5.6.3 2014/5/14 7:48 page 80 #1
Appendix W5.6.3 Root Locus with Time Delay
Pad approximant
81
solved by approximating the nonrational function e5s with a rational function. Since we are concerned with control systems and hence typically with
low frequencies, we want an approximation that will be good for small s.2
The most common means for finding such an approximation is attributed to
H. Pad. It consists of matching the series expansion of the transcendental
function e5s with the series expansion of a rational function whose numerator is a polynomial of degree p and whose denominator is a polynomial of
degree q. The result is called a (p, q) Pad approximant3 to e5s . We will
initially compute the approximants to es , and in the final result, we will
substitute Td s for s to allow for any desired delay.
The resulting (1, 1) Pad approximant (p = q = 1) is
1 (Td s/2)
.
eTd s
=
1 + (Td s/2)
(W5.6)
If we assume p = q = 2, we have five parameters and a better match is

possible. In this case, we have the (2, 2) approximant, which has the transfer
function
1 Td s/2 + (Td s)2 /12
eTd s
.
=
1 + Td s/2 + (Td s)2 /12
(W5.7)
The comparison of these approximants can be seen from their polezero

configurations as plotted in Fig. W5.6. The locations of the poles are in the
LHP and the zeros are in the RHP at the reflections of the poles.
Figure W5.6
Poles and zeros of the
Pad approximants to
es , with superscripts
identifying the
corresponding
approximants; for
example, x 1 represents
the (1,1) approximant
Im (s)
2
2
1
1
23
22
1
21
Re (s)
21
2
22
2 The nonrational function e5s is analytic for all finite values of s and so may be approximated
by a rational function. If nonanalytic functions such as s were involved, great caution would
be needed in selecting an approximation that is valid near s = 0.
3 The (p, p) Pad approximant for a delay of T sec is most commonly used and is computed by
the Matlab command [num,den] = pade(T, P).
Z15_FRAN6598_07_SE_W5.6.3 2014/5/14 7:48 page 81 #2
82 Appendix W5.6.3 Root Locus with Time Delay
Contrasting methods of
approximating delay
In some cases, a very crude approximation is acceptable. For small

delays, the (0, 1) approximant can be used, which is simply a first-order lag
given by
1
eTd s =
.
(W5.8)
1 + Td s
To illustrate the effect of a delay and the accuracy of the different approximations, root loci for the heat exchanger are drawn in Fig. W5.7 for four
cases. Notice that, for low gains and up to the point where the loci cross the
imaginary axis, the approximate curves are very close to exact. However,
the (2, 2) Pad curve follows the exact curve much further than does the
first-order lag, and its increased accuracy would be useful if the delay were
larger. All analyses of the delay show its destabilizing effect and how it limits
the achievable response time of the system.
While the Pad approximation leads to a rational transfer function, in
theory, it is not necessary for plotting a root locus. A direct application of
the phase condition can be used to plot portions of an exact locus of a system
with time delay. The phase-angle condition does not change if the transfer
function of the process is nonrational, so we still must search for values of
s for which the phase is 180 + 360 l. If we write the transfer function as
G(s) = eTd s G(s),
the phase of G(s) is the phase of G(s)

minus for s = + j. Thus
we can formulate a root-locus problem as searching for locations where
the phase of G(s)

is 180 + Td + 360 (l 1). To plot such a locus, we
would fix and search along a horizontal line in the s-plane until we found
a point on the locus, then raise the value of , change the target angle, and
repeat. Similarly, the departure angles are modified by Td , where is
the imaginary part of the pole from which the departure is being computed.
Matlab does not provide a program to plot the root locus of systems with
Im(s)
0.6
0.4
0.2
No delay
Delay approximated by
first-order lag, Eq. W5.8
Delay approximated by
Pad (2, 2) approximant, Eq. W5.7
-5s
Delay = e
(exact)
-0.5
0.5
Re(s)
-0.2
-0.4
-0.6
Figure W5.7
Root loci for the heat exchanger with and without time delay
Z15_FRAN6598_07_SE_W5.6.3 2014/5/14 7:48 page 82 #3
Appendix W5.6.3 Root Locus with Time Delay
83
delay, so we must be satisfied here with Pad approximants. Since it is

possible to plot the frequency response (or Bode plot) of delay exactly and
easily, if the designer feels that the Pad approximant is not satisfactory,
the expedient approach is to use the frequency-response design methods
described in Chapter 6.
Z15_FRAN6598_07_SE_W5.6.3 2014/5/14 7:48 page 83 #4
Appendix W6.7.2
Digital Implementation of
Example 6.15
EXAMPLE W6.1
Lead Compensation for a DC Motor

As an example of designing a lead compensator, let us repeat the design of
compensation for the DC motor with the transfer function
G(s) =
1
s(s + 1)
that was carried out in Section 5.4.1. This also represents the model of a
satellite-tracking antenna (see Fig. 3.60). This time we wish to obtain a
steady-state error of less than 0.1 for a unit-ramp input. Furthermore, we
desire an overshoot Mp < 25%.
1. Determine the lead compensation satisfying the specifications.
2. Determine the digital version of the compensation with Ts = 0.05 sec.
3. Compare the step and ramp responses of both implementations.
Solution.
1. The steady-state error is given by

Dc
ess = lim s
R(s),
s0
1 + KDc (s)G(s)
(W6.1)
where R(s) = 1/s2 for a unit ramp, so Eq. (W6.1) reduces to

1
1
ess = lim
=
.
s0 s + KDc (s)[1/(s + 1)]
KDc (0)
Therefore, we find that KD(0), which is the steady-state gain of the
compensation, cannot be less than 10 (Kv 10) if it is to meet the
error criterion, so we pick K = 10. To relate the overshoot requirement
to PM, Fig. 6.37 shows that a PM of 45 should suffice. The frequency
response of KG(s) in Fig. W6.1 shows that the PM = 20 if no phase
lead is added by compensation. If it were possible to simply add phase
without affecting the magnitude, we would need an additional phase of
only 25 at the KG(s) crossover frequency of = 3 rad/sec. However,
maintaining the same low-frequency gain and adding a compensator
zero would increase the crossover frequency; hence, more than a 25
phase contribution will be required from the lead compensation. To be
84
Z16_FRAN6598_07_SE_W6.7.2 2014/5/13 9:26 page 84 #1
Appendix W6.7.2 Digital Implementation of Example 6.15

Figure W6.1
Frequency response for
lead-compensation
design
85
200
40
100
KDc(s)G(s)
20
10
db
Magnitude, L
KG(s)
20
Lead zero
2
1
0
0.1
0.2
1
2
v (rad/sec)
10
Lead pole
Phase (deg) jL
- 905
- 1205
- 1505
- 1805
PM = 205
- 2105
- 2405
0.1
0.2
1
2
v (rad/sec)
535
10
safe, we will design the lead compensator so that it supplies a maximum phase lead of 40 . Fig. 6.53 shows that 1/ = 5 will accomplish
that goal. We will derive the greatest benefit from the compensation if
the maximum phase lead from the compensator occurs at the crossover
frequency. With some trial and error, we determine that placing the
zero at = 2 rad/sec and the pole at = 10 rad/sec causes the maximum phase lead to be at the crossover frequency. The compensation,
therefore, is
s/2 + 1
KDc (s) = 10
.
s/10 + 1
The frequency-response characteristics of L(s) = KDc (s)G(s) in
Fig. W6.1 can be seen to yield a PM of 53 , which satisfies the design
goals.
The root locus for this design, originally given as Fig. 5.24, is
repeated here as Fig. W6.2, with the root locations marked for K = 10.
The locus is not needed for the frequency-response design procedure; it
Z16_FRAN6598_07_SE_W6.7.2 2014/5/13 9:26 page 85 #2
86 Appendix W6.7.2 Digital Implementation of Example 6.15

Figure W6.2
Root locus for leadcompensation design
Im(s)
s
+1
2
KDc(s) = K s
+1
10
1
G(s) =
s(s + 1)
- 15
-10
10
8
6
4
2
-5
K = 10
-2
Re(s)
-4
-6
-8
- 10
is presented here only for comparison with the root locus design method
presented in Chapter 5. The entire process can be expedited by the use
of Matlabs sisotool routine, which simultaneously provides the root
locus and the Bode plot through an interactive GUI interface. For this
example, the Matlab statements
G=tf(1,[1 1 0]);
Dc =tf(10*[1/2 1],[1/10 1]);
sisotool(G,Dc )
will provide the plots as shown in Fig. W6.3. It also can be used to
generate the Nyquist and time-response plots if desired.
2. To find the discrete equivalent of Dc (s), we use the trapezoidal rule
given by Eq. (W4.31). That is,
Dd (z) =
2 z1
Ts z+1 /2 + 1
,
2 z1
Ts z+1 /10 + 1
(W6.2)
which, with Ts = 0.05 sec, reduces to

4.2z 3.8
.
z 0.6
This same result can be obtained by the Matlab statement
Dd (z) =
(W6.3)
sysDc = tf([0.5 1],[0.1 1]);

sysDd = c2d(sysDc , 0.05, tustin).
Because
U(z)
= KDd (z),
E(z)
the discrete control equation that results is
u(k + 1) = 0.6u(k) + 10(4.2e(k + 1) 3.8e(k)).
(W6.4)
(W6.5)
3. The Simulink block diagram of the continuous and discrete versions

of Dc (s) controlling the DC motor is shown in Fig. W6.4. The step
Z16_FRAN6598_07_SE_W6.7.2 2014/5/13 9:26 page 86 #3
Appendix W6.7.2 Digital Implementation of Example 6.15
Figure W6.3
SISOTOOL graphical interface for Example W6.1
+
Step
Ramp
0.5 s + 1
0.1 s + 1
10
Lead
control
Slider
Kc
4.2 z - 3.8
z - 0.6
10
Discrete lead
control
Slider
Kd
1
s
1
s+1
Tau 1
Tau 2
Mux
1
s
1
s+1
tau 1
tau 2
Mux 1
Output
Figure W6.4
Simulink block diagram for transient response of lead-compensation design
Z16_FRAN6598_07_SE_W6.7.2 2014/5/13 9:26 page 87 #4
87
88 Appendix W6.7.2 Digital Implementation of Example 6.15

Figure W6.5
Lead-compensation
design: (a) step
response; (b) ramp
response
Digital controller
1.2
y
0.8
0.4
0
0.5
1.5
2.5 3
Time (sec)
3.5
4.5
1.6
1.8
(a)
2
Input ramp
1.5
y
1
0.5
0
Digital controller
0
0.2
0.4
0.6 0.8 1 1.2 1.4

Time (sec)
(b)
responses of the two controllers are plotted together in Fig. W6.5a and
are reasonably close to one another; however, the discrete controller
does exhibit slightly increased overshoot, as is often the case. Both overshoots are less than 25%, and thus meet the specifications. The ramp
responses of the two controllers, shown in Fig. W6.5b, are essentially
identical, and both meet the 0.1 specified error.
Z16_FRAN6598_07_SE_W6.7.2 2014/5/13 9:26 page 88 #5
Appendix W6.8.1
Time Delay via the Nyquist
Diagram
EXAMPLE W6.2
Nyquist Plot for a System with Time Delay

Consider the system with
KeTd s
,
s
where Td = 1 sec. Determine the range of K for which the system is stable.
Solution. Because the Bode plotting rules do not apply for the phase of a
time-delay term, we will use an analytical approach to determine the key
features of the frequency response plot. As just discussed, the magnitude
of the frequency response of the delay term is unity, and its phase is
radians. The magnitude of the frequency response of the pure integrator is
1/ with a constant phase of /2. Therefore,
1
G( j) = ej(+/2)
1
(W6.6)
= ( sin j cos ).
Using Eq. (W6.6) and substituting in different values of , we can make the
Nyquist plot, which is the spiral shown in Fig. W6.6.
Let us examine the shape of the spiral in more detail. We pick a Nyquist
path with a small detour to the right of the origin. The effect of the pole at the
origin is the large arc at infinity with a 180 sweep, as shown in Fig. W6.6.
From Eq. (W6.6), for small values of > 0, the real part of the frequency
response is close to 1 because sin
= and Re[G( j)]
= 1. Similarly,
for small values of > 0, cos = 1 and Im[G( j)]

= 1/that is,
very large negative values, as shown in Fig. W6.6. To obtain the crossover
points on the real axis, we set the imaginary part equal to zero:
cos
= 0.
(W6.7)
The solution is then

(2n + 1)
0 =
,
n = 0, 1, 2, . . . .
(W6.8)
2
After substituting Eq. (W6.8) back into Eq. (W6.6), we find that

(1)n
2
G( j0 ) =
, n = 0, 1, 2, . . . .
(2n + 1)
KG(s) =
89
Z17_FRAN6598_07_SE_W6.8.1 2014/5/13 9:33 page 89 #1
90 Appendix W6.8.1 Time Delay via the Nyquist Diagram

Figure W6.6
Nyquist plot for
Example W6.2
Im[G(s)]
v = 0-
2
-p
0.6
0.4
1
p
2
3p
0.2
-1 -0.8
-0.4 -0.2
-0.2
-0.4
-0.6
0.4
0.6
0.8
Re[G(s)]
v= +
-q
1
-p
v = 0+
So the first crossover of the negative real axis is at 2/ , corresponding

to n = 0. The first crossover of the positive real axis occurs for n = 1 and
is located at 2/3 . As we can infer from Fig. W6.6, there are an infinite
number of other crossings of the real axis. Finally, for = , the Nyquist
plot converges to the origin. Note that the Nyquist plot for < 0 is the
mirror image of the one for > 0.
The number of poles in the RHP is zero (P = 0), so for closed-loop
stability, we need Z = N = 0. Therefore, the Nyquist plot cannot be allowed
to encircle the 1/K point. It will not do so as long as
1
2
< ,
K
which means that, for stability, we must have 0 < K < /2.
Z17_FRAN6598_07_SE_W6.8.1 2014/5/13 9:33 page 90 #2
(W6.9)
Appendix W6.9.2
The Inverse Nyquist Diagram
The inverse Nyquist plot is simply the reciprocal of the Nyquist plot
described in Section 6.3 and used in Section 6.4 for the definition and discussion of stability margins. It is obtained most easily by computing the inverse
of the magnitude from the Bode plot and plotting that quantity at an angle in
the complex plane, as indicated by the phase from the Bode plot. It can be
used to find the PM and GM in the same way that the Nyquist plot was used.
When |G(j)| = 1, |G1 (j)| = 1 also, so the definition of PM is identical
on the two plots. However, when the phase is 180 or +180 , the value
of |G1 (j)| is the GM directly; no calculation of an inverse is required, as
was the case for the Nyquist plot.
The inverse Nyquist plot for the system in Fig. 6.24 (Example 6.9) is
shown in Fig. W6.7 for the case where K = 1 and the system is stable.
Note that GM = 2 and PM
= 20 . As an example of a more complex case,
Fig. W6.8 shows an inverse Nyquist plot for the sixth-order case whose
Nyquist plot was shown in Fig. 6.41 and whose Nichols chart was shown in
Fig. 6.83. Note here that GM = 1.2 and PM
= 35 . Had the two crossings of
the unit circle not occurred at the same point, the crossing with the smallest
PM would have been the appropriate one to use.
Figure W6.7
Inverse Nyquist plot for
Example 6.9
Im[G -1(s)]
2.0
1.5
1.0
0.5
0
-0.5
-1.0
Re[G -1(s)]
PM
GM
-1.5
-2.0
-2
-1
91
Z18_FRAN6598_07_SE_W6.9.2 2014/5/7 21:21 page 91 #1
92 Appendix W6.9.2 The Inverse Nyquist Diagram

Figure W6.8
Inverse Nyquist plot of
the system whose
Nyquist plot is in
Fig. 6.41
Im[G -1(s)]
2.0
1.5
1.0
0.5
PM
0
-1.0
GM
-1.5
-2.0
-2
Re[G -1(s)]
-1
-0.5
-1
Z18_FRAN6598_07_SE_W6.9.2 2014/5/7 21:21 page 92 #2
Appendix W7.8
Digital Implementation
of Example 7.31
EXAMPLE W7.1
Redesign of the DC Servo Compensator

For Example 7.31, derive an equivalent discrete controller with a sampling
period of Ts = 0.1 sec (10 times the fastest pole), and compare the continuous
and digital control outputs and control efforts. Verify the design by plotting
the step response and commenting on the comparison of the continuous and
discrete responses.
Solution. The discrete equivalent for the controller is obtained from Matlab
with the c2d command, as in
nc=94.5*conv([1 7.98],[1 2.52]); % form controller numerator
dc=conv([1 8.56 59.5348],[1 10.6]); % form controller denominator
sysDc=tf(nc,dc); % form controller system description
ts=0.1;% sampling time of 0.1 sec
sysDd=c2d(sysDc,ts,' zoh' ); % convert controller to discrete time
Discrete controller
The resulting controller has the discrete transfer function

Dd (z) =
5.9157(z + 0.766)(z + 0.4586)

.
(z 0.522 0.3903j)(z + 0.3465)
The equation for the control law (with the sample period suppressed for
clarity) is
u(k + 1) = 1.3905u(k) 0.7866u(k 1) + 0.1472u(k 2)
+ e(k) 7.2445e(k 2) + 2.0782e(k 2).
Simulink simulation
A Simulink diagram for simulating both the continuous and discrete systems is shown in Fig. W7.1. A comparison of the continuous and discrete
step responses and control signals is shown in Fig. W7.2. Better agreement
between the two responses can be obtained if the sampling period is reduced.
93
Z19_FRAN6598_07_SE_W7.8 2014/5/13 9:35 page 93 #1
94 Appendix W7.8 Digital Implementation of Example 7.31

+
Step
94.5 s2 + 992.25s + 1900.3572
s3 + 19.16 s2 + 150. 27s + 631.06
Plant
10
s3 + 10 s2 + 16 s
Mux
Discrete controller
+
5.9157z2 - 7.2445z + 2.0782

z3 - 1.3905z2 + 0.7866z - 0.1472
Mux 1
Control
Plant 1
10
s3 + 10 s2 + 16 s
Figure W7.1
1.4
Digital controller
1.2
1.0
0.8
y
Figure W7.2
Comparison of step
responses and control
signals for continuous
and discrete
controllers: (a) step
responses; (b)control
signals
0.6
0.4
0.2
0
2
3
Time (sec)
(a)
7
6
5
4
u
3
2
1
0
2.1
2.2
Digital controller
0
2
3
Time (sec)
(b)
Z19_FRAN6598_07_SE_W7.8 2014/5/13 9:35 page 94 #2
Output
Appendix W7.9
Digital Implementation
of Example 7.33
EXAMPLE W7.2
Servomechanism
Matlab c2d
For Example 7.33, derive an equivalent discrete controller with a sampling

period of Ts = 0.1 sec (20 n = 20 0.05 = 0.1 sec), and compare the
continuous and digital control outputs, as well as the control efforts. Verify
the design by plotting the step response and commenting on the comparison
of the continuous and discrete responses.
Solution. The discrete equivalent for the controller is obtained from Matlab
by using the c2d command, as in
nc=conv([1 1],[8.32 0.8]); % controller numerator
dc=conv([1 4.08],[1 0.0196]); % controller denominator
sysDc=tf(nc,dc); % form controller system description
ts=0.1; % sampling time of 0.1 sec
sysDd=c2d(sysDc,ts,' zoh' ); % convert to discrete time controller
The discrete controller has the discrete transfer function
Dd (z) =
8.32z2 15.8855z + 7.5721

8.32(z 0.9903)(z 0.9191)
=
.
2
z 1.6630z + 0.6637
(z 0.998)(z 0.6665)
The equation for the control law (with sample period suppressed for
clarity) is
u(k + 1) = 1.6630u(k) + 0.6637u(k 1) + 8.32e(k + 1)
15.8855e(k) + 7.5721e(k 1).
Simulink simulation
A Simulink diagram for simulating both the continuous and discrete systems
is shown in Fig. W7.3. A comparison of the continuous and discrete step
responses and control signals is shown in Fig. W7.4. Better agreement
between the two responses can be achieved if the sampling period is reduced.
95
Z20_FRAN6598_07_SE_W7.9 2014/5/13 9:53 page 95 #1
96 Appendix W7.9 Digital Implementation of Example 7.33

+
Plant 1
1
8.32s2 + 9.12s + 0.8
s2 + 4.0996s + 0.08
Step 1
s2 + s
Mux
Discrete controller
+
Plant 1
1
8.32z2 - 15.8848z + 7.5714

z2 - 1.6630z + 0.6637
Mux 1
Control
s2 + s
Figure W7.3
1.4
Digital controller
1.2
1.0
0.8
Figure W7.4
Comparison of step
responses and control
signals for continuous
and discrete
controllers: (a) step
responses; (b) control
signals
0.6
0.4
0.2
0
2
3
Time (sec)
(a)
10
8
6
4
2
0
22
Digital controller
0
2
3
Time (sec)
(b)
Z20_FRAN6598_07_SE_W7.9 2014/5/13 9:53 page 96 #2
Output
Appendix W7.14
Solution of State Equations
In this section, we consider the solution of state variable equations. This
material is not necessary to understand the design of pole placement but will
give a deeper insight into the method of state variables. It is instructive to
consider first the unforced, or homogenous, system, which has the form
x = A(t)x,
x(0) = x0 .
(W7.1)
If the elements of A(t) are continuous functions of time, then the above
equation has a unique solution for any initial state vector x0 . There is a
useful representation for the solution of this equation in terms of a matrix,
called the transition matrix. Let i (t, t0 ) be the solution to the special initial
condition
0
0
ith row
1
x(0) = ei = 0 .
(W7.2)
..
.
0
If x(t0 ) is the actual initial condition at t0 , then we can express it in the
decomposed form
x(t0 ) = x01 e1 + x02 e2 + + x0n en .
(W7.3)
Because Eq. (W7.1) is linear, the state x(t) also can be expressed as a sum
of the solutions to the special initial condition i , as
x(t) = x01 1 (t, t0 ) + x02 2 (t, t0 ) + + x0n n (t, t0 ),
or in matrix notation, as
x(t) =
1 (t, t0 ), 2 (t, t0 ),
n (t, t0 )
(W7.4)
x01
x02
..
.
(W7.5)
x0n
So we can define the transition matrix1 to be

(t, t0 ) = 1 (t, t0 ), 2 (t, t0 ), ,
n (t, t0 )
(W7.6)
and write the solution as

x(t) = (t, t0 )x(t0 ),
(W7.7)
1 This is also referred to as the fundamental matrix of the differential equation.
97
Z21_FRAN6598_07_SE_W7.14 2014/5/13 10:06 page 97 #1
98 Appendix W7.14 Solution of State Equations

where as the name implies, the transition matrix provides the transition
between the state at time t0 to the state at time t. Furthermore, from
Eq. (W7.7), we have
d
d
[x(t)] =
[(t, t0 )] x(t0 ),
dt
dt
(W7.8)
and from Eqs. (W7.1) and (W7.8), we have

d
x(t) = Ax(t) = A(t, t0 )x(t0 ).
dt
(W7.9)
d
[(t, t0 )] = A(t, t0 ),
dt
(W7.10)
(t, t) = I.
(W7.11)
Therefore
and also
The transition matrix can be shown to have many interesting properties.

Among them are the following:
1. (t2 , t0 )
1
2.
3.
(t, )
= (t2 , t1 )(t1 , t0 );
(W7.12)
= ( , t);
(W7.13)
d
(t, ) = (t, )A( );
dt
4. det (t, t0 ) = e t0
(W7.14)
traceA( )d
(W7.15)
The second property implies that (t, ) is always invertible. What this
means is that the solution is always unique so that, given a particular value
of state at time , we can not only compute the future states from (t, )
but also past values 1 (t, ).
For the inhomogenous case, with a forcing function input u(t), the
equation is
x (t) = A(t)x(t) + B(t)u(t),
(W7.16)
and the solution is

t
x(t) = (t, t0 )x0 +
(t, )B( )u( )d .
(W7.17)
t0
We can verify this by substituting the supposed solution, Eq. (W7.17), into
the differential equation, Eq. (W7.16), as
d
d
d
x(t) = (t, t0 )x0 +
dt
dt
dt
t
(t, )B( )u( )d .
t0
Z21_FRAN6598_07_SE_W7.14 2014/5/13 10:06 page 98 #2
(W7.18)
Appendix W7.14 Solution of State Equations
99
The second term from calculus (using the Leibnitz formula) is

d
dt
t
(t, )B( )u( )d =
t0
A(t)(t, )B( )u( )d + (t, t)B(t)u(t).

t0
(W7.19)
Using the basic relations, we have
d
x(t) = A(t)(t, t0 )x0 + A(t)
dt
t
(t, )B( )u( )d
t0
+ B(t)u(t),
(W7.20)
= A(t)x(t) + B(t)u(t),
(W7.21)
which shows that the proposed solution satisfies the system equation.
For the time-invariant case
(t, t0 ) = eA(tt0 ) = (t t0 ),
where
At
A2 t 2
= I + At +
+
2!

Ak t k
k=0
k!
(W7.22)
(W7.23)
is an invertible n n exponential matrix, and by letting t = 0, we see that

e0 = I.
(W7.24)
The state solution is now

x(t) = e
A(tt0 )
t
x0 +
eA(t ) Bu( )d ,
(W7.25)
t0
and
y(t) = Cx(t) = Ce
A(tt0 )
t
x0 + C
eA(t ) Bu( )d + Du(t).
(W7.26)
t0
Suppose x(t0 ) = x0 0, then the output is given by the convolution integral

t
y(t) =
h(t )Bu( )d ,
(W7.27)
t0
where h(t) is the impulse response. In terms of the state variables matrices,
h(t) = CeAt B + D(t).
(W7.28)
While there is no uniformly best way to compute the transition matrix,

there are several methods that can be used to compute accurate approximations to it (See Moler, 2003; Franklin, Powell, and Workman, 1998). Three
of these methods are matrix exponential series, inverse Laplace transform,
Z21_FRAN6598_07_SE_W7.14 2014/5/13 10:06 page 99 #3
100 Appendix W7.14 Solution of State Equations

and diagonalization of the system matrix. In the first technique, we use
Eq. (W7.23):
A3 t 3
A2 t 2

eAt = I + At +
+
+ ,
(W7.29)
2!
3!
and the series should be computed in a reliable fashion. For the second
method, we notice that if we define
(s) = (sI A)1 ,
(W7.30)
then we can compute (s) from the A matrix and matrix algebra. Given this
matrix, we can use the inverse Laplace Transform to compute
(t) = L1 {(s)} ,

= L1 (sI A)1 .
(W7.31)
(W7.32)
The last method we mentioned operates on the system matrix. If the system
matrix can be diagonalized, that is, if we can find a transformation matrix T
so that
= T1 AT,
(W7.33)
where A is reduced to the similar but diagonal matrix
= diag {1 , 2 , , n } ,
(W7.34)
then from the series, Eq. (W7.23), we see that we need only compute scalar
exponentials, since

eAt = T1 diag e1 t , e2 t , , en t T.
(W7.35)
Z21_FRAN6598_07_SE_W7.14 2014/5/13 10:06 page 100 #4
Appendix W8.7
Discrete State-Space Design
Methods
We have seen in previous chapters that a linear, constant-coefficient continuous system can be represented by a set of first-order matrix differential
equations of the form
x = Ax + Bu,
(W8.1)
where u is the control input to the system. The output equation can be
expressed as
y = Cx + Du.
(W8.2)
The solution to these equations (see Franklin et al., 1998) is
t
A(tt0 )
x(t) = e
x(t0 ) +
eA(t ) Bu( ) d .
(W8.3)
t0
It is possible to use Eq. (W8.3) to obtain a discrete state-space representation

of the system. Because the solution over one sample period results in a
difference equation, we can alter the notation a bit (letting t = kT + T and
t0 = kT ) to arrive at a particularly useful version of Eq. (W8.3):
kT +T
eA(kT +T ) Bu( ) d .
(W8.4)
x(kT + T ) = eAT x(kT ) +
kT
This result is not dependent on the type of hold, because u is specified

in terms of its continuous time history u( ) over the sample interval. To find
the discrete model of a continuous system where the input u(t) is the output
of a ZOH, we let u( ) be a constant throughout the sample intervalthat is,
u( ) = u(kT ), kT < kT + T .
To facilitate the solution of Eq. (W8.4) for a ZOH, we let
= kT + T ,
which converts Eq. (W8.4) to
x(kT + T ) = e
AT
x(kT ) +
d Bu(kT ).
If we let
= eAT
and
=
T
0

eA d B,
(W8.5)
101
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 101 #1
102 Appendix W8.7 Discrete State-Space Design Methods

Difference equations in
standard form
Eqs. (W8.4) and (W8.2) reduce to difference equations in standard form:

x(k + 1) = x(k) + u(k),
(W8.6)
y(k) = Cx(k) + Du(k).
(W8.7)
Here x(k + 1) is a shorthand notation for x(kT + T ), x(k) for x(kT ), and
u(k) for u(kT ). The series expansion
= eAT = I + AT +
A2 T 2
A3 T 3
+
+
2!
3!
also can be written as

= I + AT ,
(W8.8)
where
=I+
AT
A2 T 2
+
+ .
2!
3!
The integral in Eq. (W8.5) can be evaluated term by term to give

=

Ak T k+1
B
(k + 1)!
k=0

Ak T k
=
TB
(k + 1)!
k=0
= T B.
(W8.9)
We evaluate by a series in the form

AT
AT
AT
AT

I+
I +
I+
,
=I+
2
3
N 1
N
Matlab c2d
which has better numerical properties than the direct series. We then find
from Eq. (W8.9) and from Eq. (W8.8). For a discussion of various methods
of numerical determination of and , see Franklin et al. (1998) and Moler
and van Loan (1978, 2003). The evaluation of the and matrices in
practice is carried out by the c2d function in Matlab.
To compare this method of representing the plant with the discrete transfer function, we can take the z-transform of Eqs. (W8.6) and (W8.7) with
J = 0 to obtain
(zI )X(z) = U(z),
(W8.10)
Y (z) = CX(z).
(W8.11)
Therefore,
Y (z)
= G(z) = C(zI )1 .
U(z)
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 102 #2
(W8.12)
Appendix W8.7 Discrete State-Space Design Methods
EXAMPLE W8.1
103
Discrete State-Space Representation of 1/s2 Plant

Use the relation in this section to verify that the discrete model of the 1/s2
plant preceded by a ZOH is that given in the solution to Example 8.4.
Solution. The and matrices can be calculated using Eqs. (W8.8) and
(W8.9). Example 7.1 (with I = 1) showed that the values for A and B are

0 1
0
A=
,B =
.
0 0
1
Because A2 = 0 in this case, we have
A2 T 2
= I + AT +
+
2!

1 0
0 1
1
=
+
T=
0 1
0 0
0

T
= I+A
TB
2!

T 0
0 1 T2
=
+
0 T
0 0
2
T
1
0
1
Hence, using Eq. (W8.12), we obtain

Y (z)
1 0
1
G(z) =
= [1 0] z
0 1
0
U(z)

T2
z+1
=
.
2 (z 1)2

,

=
T
1
T 2 /2
T
1

.
T 2 /2
T
(W8.13)
This is the same result we obtained using Eq. (8.33) and the z-transform
tables in Example 8.4.
Note that to compute Y /U, we find that the denominator of Eq. (W8.13)
is det(zI ), which was created by the matrix inverse in Eq. (W8.12). This
determinant is the characteristic polynomial of the transfer function, and
the zeros of the determinant are the poles of the plant. We have two poles at
z = 1 in this case, corresponding to two integrations in this plants equations
of motion.
We can further explore the question of poles and zeros and the statespace description by considering again the transform formulas [Eqs. (W8.10)
and (W8.11)]. One way to interpret transfer-function poles from the perspective of the corresponding difference equation is that a pole is a value of z
such that the equation has a nontrivial solution when the forcing input is
zero. From Eq. (W8.10), this interpretation implies that the linear equations
(zI )X(z) = 0
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 103 #3

have a nontrivial solution. From matrix algebra the well-known requirement for a nontrivial solution is that det(zI ) = 0. Using the system in
Example W8.1, we get

z 0
1 T
det(zI ) = det
0 z
0 1

z 1 T
= det
0
z1
= (z 1)2 = 0,
which is the characteristic equation, as we have seen. In Matlab, the poles
of the system are found by P = eig(Phi).
Along the same line of reasoning, a system zero is a value of z such that
the system output is zero even with a nonzero state-and-input combination.
Thus, if we are able to find a nontrivial solution for X(z0 ) and U(z0 ) such
that Y (z0 ) is identically zero, then z0 is a zero of the system. In combining
Eqs. (W8.10) and (W8.11), we must satisfy the requirement that

zI
X(z)
= 0.
C
0
U(z)
Once more the condition for the existence of nontrivial solutions is that
the determinant of the square coefficient system matrix be zero. For
Example W8.1, the calculation is

z 1 T T 2 /2
T T 2 /2
det 0
z1
T = det
z1
T
1
0
0
T2
(z 1)
2
T2
T2
=
z+
2
2
T2
=
(z + 1).
2
= T2 +
Thus we have a single zero at z = 1, as we have seen from the transfer

function. In Matlab, the zeros are found by Z=tzero(Phi,Gam,C,D).
Much of the algebra for discrete state-space control design is the same as
for the continuous-time case discussed in Chapter 7. The poles of a discrete
system can be moved to desirable locations by linear state-variable feedback
u = Kx
such that
det(zI + K) = c (z),
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 104 #4
(W8.14)
Appendix W8.7 Discrete State-Space Design Methods
105
provided that the system is controllable. The system is controllable if the

controllability matrix
CX = [
2
. . . n1 ]
is full rank.
A discrete full-order estimator has the form
x (k + 1) = x(k) + u(k) + L[y(k) Cx(k)],
where x is the state estimate. The error equation,
x (k + 1) = ( LC)x(k),
can be given arbitrary dynamics e (z), provided that the system is observable,
which requires that the observability matrix
C
C2
O=
..
.
Cn1
be full rank.
As was true for the continuous-time case, if the open-loop transfer
function is
Y (z)
b(z)
G(z) =
=
,
U(z)
a(z)
then a state-space compensator can be designed such that
Y (z)
Ks (z)b(z)
=
,
R(z)
c (z)e (z)
where r is the reference input. The polynomials c (z) and e (z) are
selected by the designer using exactly the same methods discussed in
Chapter 7 for continuous systems. c (z) results in a control gain K such
that det(zI + K) = c (z), and e (z) results in an estimator gain L
such that det(zI + LC) = e (z). If the estimator is structured according to Fig. 7.48a, the system zeros (z) will be identical to the estimator
poles e (z), thus removing the estimator response from the closed-loop system response. However, if desired, we can arbitrarily select the polynomial
(z) by providing suitable feed-forward from the reference input. Refer to
Franklin et al. (1998) for details.
EXAMPLE W8.2
State-Space Design of a Digital Controller

Design a digital controller for a 1/s2 plant to meet the specifications given
in Example 8.2. Use state-space design methods, including the use of an
estimator, and structure the reference input in two ways: (a) Use the error
command shown in Fig. 7.48b, and (b) use the state command shown in
Fig. 7.15 and Fig. 7.48a.
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 105 #5

Solution. We find the state-space model of the 1/s2 plant preceded by a
ZOH using the Matlab statements
sysSSc = ss([0 1;0 0], [0; 1], [1 0], 0];
T = 1;
sysSSd = c2d(sysSSc, T);
[Phi,Gam,H] = ssdata(sysSSd);
Using discrete analysis for Example 8.4, we find that the desired z-plane
roots are at z = 0.78 0.18j. Solving the discrete pole-placement problem
involves placing the eigenvalues of K, as indicated by Eq. (W8.14).
Likewise, the solution of the continuous pole-placement problem involves
placing the eigenvalues of A BK, as indicated by Eq. (7.72). Because
these two tasks are identical, we use the same function in Matlab for the
continuous and discrete cases. Therefore, the control feedback matrix K is
found by
pc = [0.78 + 0.18*j; 0.78 - 0.18*j];
K = acker(Phi,Gam,pc);
which yields
K = [0.0808 0.3996].
To ensure that the estimator roots are substantially faster than the control
roots (so that the estimator roots will have little effect on the output), we
choose them to be at z = 0.20.2j. Therefore, the estimator feedback matrix
L is found by
pe = [0.2 + 0.2*j; 0.2 - 0.2*j];
L = acker(Phi, C, pe);
which yields

1.6
.
L=
0.68
The equations of the compensation for r = 0 (regulation to xT = [0

are then
x (k + 1) = x(k) + u(k) + L[y(k) Cx(k)],
u(k) = Kx (k).
0])
(W8.15)
(W8.16)
1. For the error command structure where the compensator is placed in

the feed-forward path, as shown in Fig. 7.48b in the book, y(k) from
Eq. (W8.15) is replaced with y(k) r, so the state description of the
plant plus the estimator (a fourth-order system whose state vector is
[x x ]T ) is
A = [Phi - Gam*K; L*C Phi - Gam*K - L*C];
B = [0; 0; -L];
C = [1 0 0 0];
D = 0;
step(A,B,C,D).
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 106 #6
Summary of State-Space Design

Figure W8.1
Step response of
Example W8.2
107
1.6
Plant output y(t)
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0
0
10
15
20
Time (sec)
25
30
Command structure from Fig. 7.48b

Command structure from Fig. 7.15
The resulting step response in Fig. W8.1 shows a response similar

to that of the step responses in Fig. 8.19 in the book.
2. For the state command structure described in Section 7.9 in the book,
we wish to command the position element of the state vector so that

1
Nx =
,
0
and the 1/s2 plant requires no steady control input for a constant output y. Therefore Nu = 0. To analyze a system with this command
structure, we need to modify matrix B from the preceding Matlab statement to properly introduce the reference input r according to Fig.7.15.
The Matlab statement
B = [Gam*K*Nx; Gam*K*Nx];
channels r into both the plant and estimator equally, thus not exciting the
estimator dynamics. The resulting step response in Fig. W8.1 shows
a substantial reduction in the overshoot with this structure. In fact,
the overshoot is now about 5%, which is expected for a second-order
system with
= 0.7. The previous designs all had considerably greater
overshoot because of the effect of the extra zero and pole.
SUMMARY OF STATE-SPACE DESIGN
The continuous state-space form of a differential equation,

x = Ax + Bu,
y = Cx + Du,
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 107 #7

has a discrete counterpart in the difference equations
x(k + 1) = x(k) + u(k),
y(k) = Cx(k) + Du(k),
where
= eAT
T

=
eA d B.
0
These matrices can be computed in Matlab by [Phi, Gam] = c2d(F,G,H,J)

and used in state-space discrete design methods.
The pole placement and estimation ideas are identical in the continuous
and discrete domains.
PROBLEMS
W8.1
In Problem 8.11 in the book, we dealt with an experiment in magnetic

levitation described by Eq. (8.62) that reduces to
x = 1000x + 20i.
Let the sampling time be 0.01 sec.
(a) Use pole placement to design a controller for the magnetic levitator so
that the closed-loop system meets the following specifications: settling
time, ts 0.25 sec, and overshoot to an initial offset in x that is less
than 20%.
(b) Plot the step response of x, x , and i to an initial displacement in x.
(c) Plot the root locus for changes in the plant gain, and mark the pole
locations of your design.
(d) Introduce a command reference input r (as discussed in Section 7.9) that
does not excite the estimate of x. Measure or compute the frequency
response from r to the system error r x and give the highest frequency for which the error amplitude is less than 20% of the command
amplitude.
W8.2
Servomechanism for Antenna Elevation Control: Suppose it is desired to

control the elevation of an antenna designed to track a satellite. A photo of
such a system is shown in Fig. W8.2, and a schematic diagram is depicted in
Fig. W8.3. The antenna and drive parts have a moment of inertia J and damping B, arising to some extent from bearing and aerodynamic friction, but
mostly from the back emf of the DC drive motor. The equation of motion is
J + B = Tc + Td ,
where
Tc = net torque from the drive motor,
(W8.17)
Td = disturbance torque due to wind.
(W8.18)
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 108 #8
Problems
Figure W8.2
Satellite-tracking
antenna (Courtesy
Space Systems/Loral)
Figure W8.3
Schematic diagram of
satellite-tracking
antenna
u
If we define
B
= a,
J
the equation simplifies to
u=
Tc
,
B
and
T
wd = d ,
B
1
+ = u + wd .
a
After using the Laplace transformation, we obtain
(s) =
1
[u(s) + wd (s)],
s(s/a + 1)
or with no disturbance,
(s)
1
=
= G2 (s).
u(s)
s(s/a + 1)
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 109 #9
109

With u(k) applied through a ZOH, the transfer function for an equivalent
discrete-time system is
G2 (z) =
z+b
(z)
=K
,
u(z)
(z 1)(z eaT )
where
1 eaT aTeaT
aT 1 + eaT
, b=
.
a
aT 1 + eaT
and write the continuous-time state equations
(a) Let a = 0.1 and x1 = ,
for the system.
(b) Let T = 1 sec, and find a state feedback gain K for the equivalent
discrete-time system that yields closed-loop polescorresponding
to the

3
following points in the s-plane: s = 1/2 j
2 . Plot the step
K=
response of the resulting design.

(c) Design an estimator: Select L so that e (z) = z2 .
(d) Using the values for K and L computed in parts (b) and (c) as the
gains for a combined estimator/controller, introduce a reference input
that will leave the state estimate undisturbed. Plot the response of the
closed-loop system due to a step change in the reference input. Also
plot the system response to a step wind-gust disturbance.
(e) Plot the root locus of the closed-loop system with respect to the plant
gain, and mark the locations of the closed-loop poles.
W8.3
Tank Fluid Temperature Control: The temperature of a tank of fluid with a

constant inflow and outflow rate is to be controlled by adjusting the temperature of the incoming fluid. The temperature of the incoming fluid is
controlled by a mixing valve that adjusts the relative amounts of hot and
cold supplies of the fluid (see Fig. W8.4). The distance between the valve
and the point of discharge into the tank creates a time delay between the
application of a temperature change at the mixing valve and the discharge
of the flow with the changed temperature into the tank. The differential
equation governing the tank temperature is
T e =
1
(qin qout ),
cM
where
Te = tank temperature,
c = specific heat of the fluid,
Figure W8.4
Tank temperature
control
Hot
Tec
Negligible heat
loss in walls
Tei
Cold
Mixing
value
Tank fluid at
temperature Te
Te
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 110 #10
Problems
111
M = fluid mass contained in the tank,

in Tei ,
qin = cm
out Te ,
qout = cm
out ),
m
= mass flow rate (m
in = m
Tei = temperature of fluid entering tank.
However, the temperature at the input to the tank at time t is equal to the control temperature d seconds in the past. This relationship may be expressed
as
Tei (t) = Tec (t d ),
where
d = delay time,
Tec = temperature of fluid immediately after the control valve and directly
controllable by the valve.
Combining constants, we obtain
T e (t) + aTe (t) = aTec (t d ),
where
.
M
The transfer function of the system is thus
a=
e ds
Te (s)
= G3 (s).
=
s/a + 1
Tec (s)
To form a discrete transfer function equivalent to G3 preceded by a ZOH,
we must compute

ed s
1 esT
G3 (z) = Z
.
s
s/a + 1
We assume that for some integer l, d = lT mT , where 0 m < 1. Then

elsT emsT
1 esT
G3 (z) = Z
s
s/a + 1

emsT
1
l
= (1 z )z Z
s(s/a + 1)

emsT
emsT
1
l
= (1 z )z Z
s
s+a

z1 1
=
Z{1(t + mT ) ea(t+mT ) 1(t + mT )}
z
zl

z
eamT z
z1 1
=
z
z 1 z eaT
zl
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 111 #11

1
= l
z

=
(1 eamT )z + eamT eaT

z eaT

z+
1 eamT
,
z eaT
zl
and
eamT eaT
.
1 eamT
The zero location varies from = at m = 0 to = 0 as m 1.
Note also that G3 (1) = 1.0 for all a, m, and l. For the specific values
d = 1.5, T = 1, a = 1, l = 2, and m = 21 , the transfer function reduces to
=
z + 0.6065
.
G3 (z) = 0.3935 2
z (z 0.3679)
(a)
(b)
(c)
(d)
(e)
W8.4
Write the discrete-time system equations in state-space form.

Design a state feedback gain that yields c (z) = z3 .
Design a state estimator with e (z) = z3 .
Plot the root locus of the system with respect to the plant gain.
Plot the step response of the system.
Consider the linear equation Ax = b, where A is an n n matrix. When b

is given, one way of solving for x is to use the discrete-time recursion
x(k + 1) = (I + cA)x(k) cb,
where c is a scalar to be chosen.
(a) Show that the solution of Ax = b is the equilibrium point x of the
discrete-time system. An equilibrium point x of a discrete-time system
x(k + 1) = f(x(k)) satisfies the relation x = f(x ).
(b) Consider the error e(k) = x(k) x . Write the linear equation that
relates the error e(k + 1) to e(k).
(c) Suppose |1 + ci (A)| < 1, i = 1, . . . , n, where i (A) denotes the ith
eigenvalue of A. Show that, starting from any initial guess x0 , the algorithm converges to x . [Hint: For any matrix B, i (I + B) = 1 + i (B).]
W8.5
The open-loop plant of a unity feedback system has the transfer function
G(s) =
1
.
s(s + 2)
Determine the transfer function of the equivalent digital plant using a sampling period of T = 1 sec, and design a proportional controller for the
discrete-time system that yields dominant closed-loop poles with a damping
ratio of 0.7.
W8.6
Write a computer program to compute and from A, B, and the sample period T . It is okay to use Matlab, but dont use C2D. Write code in
Matlab to compute the discrete matrices using the relations developed in
this chapter. Use your program to compute and when
(a)

A=
1
0
0
2

,
B=
1
1

,
T = 0.2 sec,
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 112 #12
Problems
(b)
113

1
3 2
, T = 0.2sec.
, B=
0
1
0
Consider the following discrete-time system in state-space form:

x1 (k + 1)
x1 (k)
0
0
1
u(k).
+
=
10
x2 (k + 1)
x2 (k)
0 1

A=
W8.7
Use state feedback to relocate all of the systems poles to 0.5.

W8.8

2

1 T
T /2
and =
.
0 1
T
(a) Find a transformation matrix T so that, if x = Tw, the state equations
for w will be in control canonical form.
(b) Compute the gain Kw so that, if u = Kw w, the characteristic equation
will be c (z) = z2 1.6z + 0.7.
(c) Use T from part (a) to compute Kx , which is the feedback gain
required by the state equations in x to achieve the desired characteristic
polynomial.

Let
=
W8.9
Consider a system whose plant transfer function is 1/s2 and has a piecewise
constant input of the form
u(t) = u(kT ),
kT t < (k + 1)T .
(a) Show that, if we restrict attention to the time instants kT , k = 0, 1, 2, . . .,

the resulting sampled-data system can be described by the equations

T
x1 (k + 1)
1 0 x1 (k)
=
+ 2
u(k).
T 1 x2 (k)
x2 (k + 1)
T /2
y(k) = [0
1][x1 (k) x2 (k)]T .
(b) Design a second-order estimator that will always drive the error in the
estimate of the initial state vector to zero in time 2T or less.
(c) Is it possible to estimate the initial state exactly with a first-order
estimator? Justify your answer.
W8.10 In this problem, you will show how to compute by changing states so
that the system matrix is diagonal.
(a) Using an infinite series expansion, compute eAT for

1
0
A=
.
0
2
(b) Show that if A = TAT1 for some nonsingular transformation matrix
T, then
eAT = TeAT T1 .
(c) Show that if

A=
3
2
1
0

,
there exists a T such that TAT1 = A. (Hint: Write TA = AT, assume

four unknowns for the elements of T, and solve. Next show that the
columns of T are the eigenvectors of A.)
(d) Compute eAT .
z22_fran6598_07_se_w8.7 2014/5/13 10:24 page 113 #13

FPE7e Appendices

Uploaded by

Copyright:

Available Formats

FPE7e Appendices

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

FPE7e Appendices

Uploaded by

Copyright:

Available Formats

Appendix WA

WA.1 Denition of a Complex Number

Note that the imaginary part of A is itself a real number.

where rcalled the magnitude, modulus, or absolute value of Ais the

Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 1 #1

2 Appendix WA A Review of Complex Variables

Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 2 #2

WA.2 Algebraic Manipulations

WA.2 Algebraic Manipulations

WA.2.2 Complex Multiplication

The product of two complex numbers may be obtained graphically using

WA.2.3 Complex Division

From Eq. (WA.4) it follows that

Also, if A1 = r1 e j1 and A2 = r2 e j2 , then

where |A1 A2 | = r1 r2 and arg(A1 A2 ) = 1 + 2 , and

Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 3 #3

4 Appendix WA A Review of Complex Variables

Frequency Response of First-Order System

Therefore, the magnitude and phase are

arg(G( j)) = tan1

WA.3 Graphical Evaluation of Magnitude and

This may be explained as follows. The term s + z1 is a vector addition of its

Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 4 #4

WA.4 Differentiation and Integration

WA.4 Differentiation and Integration

Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 5 #5

6 Appendix WA A Review of Complex Variables

WA.5 Eulers Relations

We collect the terms involving A to obtain

where c is a constant of integration. If we let = 0 in Eq. (WA.30), we find

From Eqs. (WA.31) and (WA.32) it follows that

WA.6 Analytic Functions

WA.7 Cauchys Theorem

Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 6 #6

WA.8 Singularities and Residues

There is a corollary to this theorem: Let C1 and C2 be two paths connecting

WA.8 Singularities and Residues

The coefficient A1 is called the residue of G(s) at s0 and may be

Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 7 #7

8 Appendix WA A Review of Complex Variables

WA.9 Residue Theorem

WA.10 The Argument Principle

Number of Poles and Zeros

Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 8 #8

WA.10 The Argument Principle

where f (s) is analytic and f (s0 )  = 0. If we differentiate Eq. (WA.43), we

Therefore, G (s)/G(s) has a pole at s = s0 with residue K. This analysis may

The Argument Principle

Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 9 #9

10 Appendix WA A Review of Complex Variables

WA.11 Bilinear Transformation

If we substitute for w in terms of s in Eq. (WA.53), we get

Z01_FRAN6598_07_SE_WA 2014/5/7 0:32 page 10 #10

WB.1 Matrix Denitions

a11 a12 a1n

WB.2 Elementary Operations on Matrices

cij = aij + bij .

Commutative law for

Two matrices can be multiplied if they are compatible. Let A = m n and

where f (s) is analytic and f (s0 ) = 0. If we differentiate Eq. (WA.43), we

Therefore, G (s)/G(s) has a pole at s = s0 with residue K. This analysis may