Additional Lecture Notes for the course SC4090
Discretetime systems analysis Ton van den Boom
October 2, 2006
2
Contents 1 Introduction 1.1 Discretetime versus continuoustime . . 1.2 Sampling and Interpolation . . . . . . . 1.3 Discretetime systems . . . . . . . . . . . 1.3.1 Difference Equations . . . . . . . 1.3.2 Discretetime statespace models 1.3.3 Linear discretetime models . . .
. . . . . .
5 5 7 12 12 13 14
2 Operational methods for discretetime linear systems 2.1 Discretetime linear system operators . . . . . . . . . . . . . . . . . . . . . 2.2 Transformation of discretetime systems . . . . . . . . . . . . . . . . . . .
17 17 21
3 System Properties and Solution techniques 3.1 Singularity input functions . . . . . . . . . . . . . . . . 3.2 Classical solution of linear difference equations . . . . . 3.2.1 Homogeneous solution of the difference equation 3.2.2 Particular solution of the difference equation . . 3.3 Discretetime Convolution . . . . . . . . . . . . . . . .
. . . . .
23 23 24 25 25 27
. . . . . . . . .
31 31 31 33 34 34 35 36 37 38
5 The Discretetime transfer function 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 singleinput singleoutput systems . . . . . . . . . . . . . . . . . . . . . . .
41 41 41
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . .
. . . . . .
. . . . .
. . . . . .
. . . . .
. . . . . .
. . . . .
. . . . . .
. . . . .
. . . . . .
. . . . .
4 Solution of discretetime linear state equations 4.1 State variable responses of discretetime linear systems . . . . . . 4.1.1 Homogeneous State responses . . . . . . . . . . . . . . . . 4.1.2 The forced response of discretetime linear systems . . . . 4.1.3 The system output response of discretetime linear systems 4.1.4 The discretetime transition matrix . . . . . . . . . . . . . 4.1.5 System Eigenvalues and Eigenvectors . . . . . . . . . . . . 4.1.6 Stability of discretetime linear systems . . . . . . . . . . . 4.1.7 Transformation of state variables . . . . . . . . . . . . . . 4.2 The response of linear DT systems to the impulse response . . . .
3
. . . . . .
. . . . . . . . . . . . . .
. . . . . .
. . . . . . . . . . . . . .
. . . . . .
. . . . . . . . . . . . . .
. . . . . .
. . . . . . . . . . . . . .
4
CONTENTS 5.3 5.4
relationship to the transfer function . . . . . . . . . System poles and zeros . . . . . . . . . . . . . . . . 5.4.1 System poles and the homogeneous response 5.4.2 System stability . . . . . . . . . . . . . . . . 5.4.3 State space formulated systems . . . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
42 43 44 45 47
Chapter 1 Introduction 1.1
Discretetime versus continuoustime
Engineers and Physical scientists have for many years utilized the concept of a system to facilitate the study of the interaction between forces and matter. A system is a mathematical abstraction that is devised to serve as a model for a dynamic phenomenon. It represents the dynamic phenomenon in terms of mathematical relations among three sets of variables known as the input, the output, and the state. The input represents, in the form of a set of time functions or sequences, the external forces that are acting upon the dynamic phenomenon. In similar form, the output represents the measures of the directly observable behavior of the phenomenon. Input and output bear a causeeffect relation to each other; however, depending on the nature of the phenomenon, this relation may be strong or weak. A basic characteristic of any dynamic phenomenon is that the behavior at any time is traceable not only to the presently applied forces but also to those applied in the past. We may say that a dynamic phenomenon possesses a ”memory” in which the effect of past applied forces is stored. In formulating a system model, the state of the system represents, as a vector function of time, the instantaneous content of the ”cells” of this memory. Knowledge of the state at any time t, plus knowledge of the forces subsequently applied is sufficient to determine the output (and state) at any time t ≥ t0 . As an example, a set of moving particles can be represented by a system in which the state describes the instantaneous position and momentum of each particle. Knowledge of position and momentum, together with knowledge of the external forces acting on the particles (i.e., the system input) is sufficient to determine the position and momentum at any future time. A system is, of course, not limited to modeling only physical dynamic phenomena; the concept is equally applicable to abstract dynamic phenomena such as those encountered in economics or other social sciences. If the time space is continuous, the system is known as a continuoustime system. However, if the input and state vectors are defined only for discrete instants of time k, where k ranges 5
6
CHAPTER 1. INTRODUCTION
over the integers, the time space is discrete and the system is referred to as a discretetime system. We shall denote a continuoustime function at time t by f (t). Similarly, a discretetime function at time k shall be denoted by f (k). We shall make no distinction between scalar and vector functions. This will usually become clear from the context. Where no ambiguity can arise, a function may be represented simply by f .
Figure 1.1: Illustration of discretetime functions and quantized functions. (a) discretetime function, (b) quantized function, (c) quantized, discretetime function
It is important to distinguish between functions whose argument is discrete (i.e., functions of a discrete variable) and those that in themselves vary over a discrete set. Functions of the latter type will be referred to as quantized functions and the systems in which they appear will be called quantizeddata systems. Thus the function illustrated in Fig. 1.1.a is a discretetime function and the one shown in Fig.1.1.b is a quantized function. Note that in Fig. 1.1 f (k) can range only over discrete values (which need not necessarily be uniformly spaced). A discretetime function may, of course, be quantized as well; this is shown in Fig. 1.1.c. In certain continuoustime systems, some state variables are allowed to change only at discrete instants of time tk , where k ranges over the integers and where the spacing between successive instants may be arbitrary or uniform. Such systems are in effect a kind of hybrid between a continuous time and a discretetime system. They are encountered whenever
1.2. SAMPLING AND INTERPOLATION
7
a continuoustime function is sampled at discrete instants of time. Although they are frequently analyzed most easily by treating them as discretetime systems, they differ from discretetime systems in that special consideration may have to be given to the instants of time at which the sampling occurs. They have been given a special name and are known as sampleddata systems. We shall be concerned here primarily with the analysis of discretetime systems. Our interest in these systems is motivated by a desire to predict the performance of the physical devices for which this kind of system is an appropriate model. This is, however, not the only motive for studying discretetime systems. A second, nearly as important reason is that there are many continuoustime systems that are more easily analyzed when a discretetime model is fitted to them. A common example of this is the simulation of a continuoustime system by a digital computer. Further, a considerable body of mathematical theory has been developed for the analysis of discretetime systems. Much of this is valuable for gaining insight into the theory of continuoustime systems as well as discretetime systems. We note that a continuoustime function can always be viewed as the limit of a time sequence whose spacing between successive terms is approaching zero.
1.2
Sampling and Interpolation
Many physical signals, such as electrical voltages produced by a sound or image recording instrument or a measuring device, are essentially continuoustime signals. Computers and related devices operate on a discretetime axis. Continuoustime signals that are to be processed by such devices therefore first need be converted to discretetime signals. One way of doing this is sampling. Definition 1 Sampling. Let Tcon = IR be the continuoustime axis, and let Tdis = ZZ be the discretetime axis. Let x(t) be a continuoustime signal with t ∈ Tcon . Then, sampling the continuoustime signal x(t), t ∈ Tcon on the discretetime axis Tdis with sampling period T , results in the discretetime signal x∗ (k) defined by x∗ (k) = x(k T )
for all
k ∈ Tdis
(1.1)
A device that performs the sampling operation is called a sampler. Example 2 (Sampled real harmonic) Let the continuoustime signal x, given by x(t) = 1/2[1 − cos(2πt)], t ∈ IR be sampled on the uniformly spaced discretetime axis ZZ. This results in the sampled signal x∗ given by x∗ (k) = 1/2[1 − cos(2πkT )], k ∈ ZZ
8
CHAPTER 1. INTRODUCTION
original signal x
sampled signal x∗
Figure 1.2: Sampling. Left: a continuoustime signal. Right: its sampled version.
The converse problem of sampling presents itself when a discretetime device, such as a computer, produces signals that need drive a physical instrument requiring a continuoustime signal as input. Suppose that a discretetime signal x∗ is defined on the discretetime axis Tdis and that we wish to construct from x∗ a continuoustime signal x defined on the continuoustime axis Tcon ⊃ Tdis . There obviously are many ways to do this. We introduce a particular class of conversions from discretetime to continuoustime signals, for which we reserve the term interpolation. This type of conversion has the property that the continuoustime signal x agrees with the discretetime signal x∗ at the sampling times. Definition 3 Interpolation. Let Tcon = IR be the continuoustime axis, and let Tdis = ZZ be the discretetime axis. Let x∗ (k) be a discretetime signal with k ∈ Tdis with sampling period T . Then, any continuoustime signal x(t), t ∈ IR is called an interpolation of x∗ on Tcon , if x(k T ) = x∗ (k)
for all
k ∈ Tdis
(1.2)
Another way of saying that x is an interpolation of x∗ is the statement that sampling the continuoustime signal x generated by interpolating the discretetime signal x∗ on Tdis , reproduces the discretetime signal x∗ . Clearly, there is no unique interpolation for a given discretetime signal x∗ . Suppose that x∗ is defined on the uniformly sampled discretetime axis ZZ. Then, one possible interpolation method is step interpolation as illustrated in Fig. 1.3(a). Interpolation is done using an interpolating function. An interpolating function is any function i: IR → IR such that
1 for t = 0 i(t) = 0 for t = k, where k 6= 0, k ∈ ZZ arbitrary elsewehere
t ∈ IR
9
1.2. SAMPLING AND INTERPOLATION
Figure 1.3: Interpolation. Top: step interpolation. Bottom: linear interpolation.
If x∗ (k), k ∈ ZZ is a discretetime signal defined on the time axis Tdis , and i an interpolating function, the continuoustime signal x given by x(t) =
X
x∗ (n) i(t/T − n),
t ∈ IR
n∈Z Z
is an interpolation of x∗ . The reason is that by setting t = k T , with k ∈ ZZ, it follows that x(k T ) =
X
x∗ (n) i(kT /T − n) = x∗ (k),
k ∈ ZZ
n∈Z Z
Step interpolation is achieved with the step interpolating function istep (t) =
(
1 for 0 ≤ t < 1 0 otherwise
t ∈ IR
while linear interpolation is obtained with the linear interpolating function ilin (t) =
(
1 − t for t < 1 0 otherwise
t ∈ IR
Another interpolation function is the sinc interpolating function isinc (t) = sinc(π t) Here, sinc: IR → IR is the function defined by sin(t) for t 6= 0 sinc(t) = t 1 for t = 0
t ∈ IR
Graphs of these interpolating functions are given in Fig. 1.4.
10
CHAPTER 1. INTRODUCTION
Figure 1.4: Interpolating functions. Top left: step interpolating function. Top right: linear interpolating function. Bottom: sinc interpolating function.
Example 4 (Interpolation) Let a sampled signal x∗ with sampling period T be given by πk cos ∗ x (k) = 8 0
!
for − 4 ≤ k < 4
k ∈ ZZ
otherwise
A plot of x∗ is given in Fig. 1.5(a). Interpolation with the step interpolating function results in the staircaselike continuoustime signal of Fig. 1.5(b). Interpolation with the linear interpolating function leads to the signal of Fig. 1.5(c), which is obtained by connecting the sampled values of the original discretetime signal by straight lines. Interpolation with the sinc interpolating function, finally, yields the signal of Fig. 1.5(d), which is smooth within the interval [−4T, 4T ] but shows , “ripple” outside it.
Remark: Are sampling and interpolation inverses of each other? Interpolation of a discretetime signal followed by sampling on the original discretetime axis results in exact reconstruction of the original discretetime signal. On the other hand, sampling a continuoustime signal followed by interpolation generally does not reproduce the original continuoustime signal. Evidently sampling is an inverse operation to interpolation, but interpolation is not inverse to sampling. Figure 1.6 illustrates this.
1.2. SAMPLING AND INTERPOLATION
11
Figure 1.5: Interpolation. (a) A discretetime signal. (b) Step interpolation (c) Linear interpolation. (d) Sinc interpolation.
Figure 1.6: Sampling and interpolation. Top: sampling is the inverse of interpolation. Bottom: interpolation is not the inverse of sampling
12
CHAPTER 1. INTRODUCTION
Remark: AnalogtoDigital and digitaltoanalog conversion (a) A/D conversion. Computers and other digital equipment do not only operate on a discretetime axis but are also limited to finite signal ranges. Thus, apart from sampling, conversion of realvalued continuoustime signals to input for digital equipment also involves quantization. The combined process of sampling and quantization is called analogtodigital (A/D) conversion. (b) D/A conversion. The inverse process of converting a quantized discretetime signal to a realvalued continuoustime signal is called digitalto analog (D/A) conversion. Devices that perform D/A conversion by step interpolation are known as zeroorder hold circuits. They can function in ”real time” meaning that given the sampled signal up to, and including, time t, the continuoustime signal up to and including that same time may be generated. Devices that perform D/A conversion by linear interpolation are known as firstorder hold circuits. They cannot precisely function in real time because two successive sampled signal values have to be received before it is known how the continuoustime signal goes. Firstorder hold circuits therefore introduce a delay, equal to the sampling period T . Sinc interpolation, finally, cannot be implemented in real time at all because all sampled signal values need be received before any point of the continuoustime signal (except at the sampling times) may be computed. Sinc interpolation, however, has great theoretical importance.
1.3
Discretetime systems
1.3.1
Difference Equations
In continuoustime the relationships between different model variables are described with the help of differential equations (see Rowell & Wormley [4]). In discretetime this can be done using difference equations. There are two different ways of describing these difference equations. One is to directly relate inputs u to outputs y in one difference equation. In principle it looks like this: g(y(k + n), y(k + n − 1), . . . , y(k), u(k + m), u(k + m − 1), . . . , u(k)) = 0
(1.3)
where g(·, ·, . . . , ·) is an arbitrary, vectorvalued, nonlinear function. The other way is to write the difference equation as a system of firstorder difference equations by introducing a number of internal variables. If we denote these internal variables by x1 (k), . . . , xn (k) and introduce the vector notation
x1 (k) .. x(t) = . xn (k)
(1.4)
1.3. DISCRETETIME SYSTEMS
13
we can, in principle, write a system of firstorder difference equations as x(k + 1) = f (x(k), u(k))
(1.5)
In (1.5), f (x, u) is a vector function with n components:
f1 (x, u) .. f (x, u) = . fn (x, u)
(1.6)
The functions fi (x, u) are in turn functions of n + m variables, the components of the x and u vectors. Without vector notation, (1.5) becomes x1 (k + 1) = f1 (x1 (k), . . . , xn (k), u1(k), . . . , um (k)) x2 (k + 1) = f2 (x1 (k), . . . , xn (k), u1(k), . . . , um (k)) .. .
(1.7)
xn (k + 1) = fn (x1 (k), . . . , xn (k), u1 (k), . . . , um(k)) The outputs of the model can then be calculated from the internal variables xi (k) and the inputs ui(k): y(k) = h(x(k), u(k))
(1.8)
which written in longhand means y1 (k) = h1 (x1 (k), . . . , xn (k), u1(k), . . . , um (k)) y2 (k) = h2 (x1 (k), . . . , xn (k), u1(k), . . . , um (k)) .. .
(1.9)
yp (k) = hp (x1 (k), . . . , xn (k), u1 (k), . . . , um (k))
1.3.2
Discretetime statespace models
For a dynamic system the output depends on all earlier input values. This leads to the fact that it is not enough to know u(k) for k ≥ k0 in order to be able to calculate the output y(k) for k ≥ k0 . We need information about the system. By the state of the system at time k0 we mean an amount of information such that with this state and the knowledge of u(k), k ≥ k0 , we can calculate y(k), k ≥ k0 . This definition is well in line with the everyday meaning of the word “state.” It is also obvious from the definition of state that this concept will play a major role in the simulation of the model. The state is exactly the information that has to be stored and updated at the simulation in order to be able to calculate the output. Consider a general system of firstorder difference equations (1.5) with the output given by (1.8): x(k + 1) = f (x(k), u(k)) y(k) = h(x(k), u(k))
(1.10) (1.11)
14
CHAPTER 1. INTRODUCTION
For this system the vector x(k0 ) is a state at time k0 . The difference equation (1.10)(1.11) with x(k0 ) = x0 has a unique solution for k ≥ k0 . The solution can easily be found using successive substitution: Assume that we know x(k) at time k0 and u(k) for k ≥ k0 . We can then according to (1.10) calculate x(k + 1). From this value we can use u(k + 1) to calculate x(k + 2), again using (1.10). We can continue and calculate x(k) for all k > k0 . The output y(k), k ≥ k0 , can then also be computed according to (1.11). We have thus established that the variables x1 (k), . . . , xn (k) or, in other words, the vector
x1 (k) x(k) = ... xn (k)
in the internal model description (1.10)(1.11) is a state for the model. Herein lies the importance of this model description for simulation. The model (1.10)(1.11) is therefore called a statespace model, the vector x(k) the state vector, and its components xi (k) state variables. The dimension of x(k), that is, n, is called the model order. Statespace models will be our standard model for dynamic discretetime systems. In conclusion, we have the following model: Discretetime state space models x(k + 1) = f (x(k), u(k)) y(k) = h(x(k), u(k))
k = 0, 1, 2, . . .
(1.12) (1.13)
u(k): input at time k , an mdimensional column vector. y(k): output at time k , a pdimensional column vector. x(k): state at time k , an ndimensional column vector. The model is said to be nth order. For a given initial value x(k0 ) = x0 , (1.12)(1.13) always has a unique solution.
1.3.3
Linear discretetime models
The model (1.12)(1.13) is said to be linear if f (x, u) and h(x, u) are linear functions of x and u: f (x, u) = Ax + Bu h(x, u) = Cx + Du Here the matrices have the following dimensions A:n×n
B :n×m
C :p×n
D :p×m
(1.14) (1.15)
1.3. DISCRETETIME SYSTEMS
15
If these matrices are independent of time the model (1.14)(1.15) is said to be linear and time invariant. The model (1.3) is said to be linear if g(·, ·, . . . , ·) is a linear function in y and u: g(y(k+n), y(k+n−1), . . . , y(k), u(k+m), u(k+m−1), . . . , u(k)) = a0 y(k+n)+an−1 y(k+n−1)+. . .+a0 y(k) −bm u(k+m)−bm−1 u(k+m−1)−. . .−b0 u(k) = 0 or a0 y(k+n)+an−1 y(k+n−1)+. . .+a0 y(k) = bm u(k+m)+bm−1 u(k+m−1)+. . .+b0 u(k) (1.16) Example 5 (Savings account) Suppose that y(k) represents the balance of a savings account at the beginning of day k , and u(k) the amount deposited during day k. If interest is computed and added daily at a rate of α · 100%, the balance at the beginning of day k + 1 is given by y(k + 1) = (1 + α)y(k) + u(k), k = 0, 1, 2, · · ·
(1.17)
This describes the savings account as a discretetime system on the time axis ZZ+ . If the interest rate α does not change with time, the system is timeinvariant; otherwise, it is timevarying. Example 6 (Firstorder exponential smoother) In signal processing sometimes a procedure called exponential smoothing is used to remove unwanted fluctuations from observed time series, such as the Dow Jones index, the daily temperature, or measurements of a human’s blood pressure. Let u(k), k ∈ ZZ be a discretetime signal, then, as we observe the successive values of the signal u we may form from these values another signal y, k ∈ ZZ, which is a smoothed version of the signal u, according to y(k + 1) = a y(k) + (1 − a) u(k + 1)
(1.18)
Here, a is a constant such that 0 < a < 1. At each time k + 1, the output y(k + 1) is formed as a weighted average of the new input u(k + 1) at time k + 1 and the output y(k) at preceding time instant k. The closer the constant a is to 1, the more the preceding output value is weighted and the “smoother” the output is. Figure 1.7 illustrates the smoothing effect.
16
CHAPTER 1. INTRODUCTION
2
u(k) −−−>
1
0
−1
−2
0
5
10
15
20
25 k −−−>
30
35
40
45
50
0
5
10
15
20
25 k −−−>
30
35
40
45
50
1
y(k) −−−>
0.5 0 −0.5 −1 −1.5
Figure 1.7: Exponential smoothing (a = 0.8)
Chapter 2 Operational methods for discretetime linear systems In this chapter we will consider operational methods for discretetime linear systems. For the continuoustime case see Rowell & Wormley [4].
2.1
Discretetime linear system operators
Consider a discretetime signal f (k), for k ∈ ZZ. All linear timeinvariant systems may be represented by an interconnection of the following primitive operators: 1. The constant scaling operator: The scaling operator multiplies the input function by a constant factor. It is denoted by the value of the constant, either numerically or symbolically, for example 2.0{}, or α{}. 2. The forward shift operator: The shift operator, designated Z{}, shifts the input signal in one step ahead in time: y(k) = Z{f (k)} = f (k + 1) 3. The backward shift operator: The shift operator, designated Z −1 {}, shifts the input signal in one step back in time: y(k) = Z −1 {f (k)} = f (k − 1) 4. The difference operator: The difference operator, written ∆{}, generates the increment of the input f (k): y(k) = ∆{f (k)} = f (k) − f (k − 1) 17
18CHAPTER 2. OPERATIONAL METHODS FOR DISCRETETIME LINEAR SYSTEMS 5. The summation operator: The summation operator, written ∆−1 {}, generates the summation of the input f (k): y(k) = ∆−1 {f (k)} =
k X
f (m)
m=0
where it is assumed that at time k < 0, the output y(k) = 0. If an initial condition y(−1) is specified, y(k) = ∆−1 {f (k)} + y(−1) and a separate summing block is included at the output of the summation block to account for the initial condition. 6. The identity operator: The identity operator leaves the value f (k) unchanged, that is: y(k) = I{f (k)} = f (k) 7. The null operator: The null operator N{} identically produces an output of zero for any input, that is, y(k) = N{f (k)} = 0 Remark 1: Note that by substituting f (k − 1) = Z −1 {f (k)} we find that ∆{} = I{} − Z −1 {} Remark 2: The discretetime operators Z, Z −1 , ∆, ∆−1 , I and N are linear operators. Therefore all properties for linear operators will also hold for these discretetime operators.
Block diagram of a discretetime system The block diagram of the statespace description of a linear discretetime system is equivalent to the one for continuoustime systems if we substitute the backward shift block Z −1 for the integral block. (In the book of Rowell & Wormley [4] it is the block S −1 as in Fig. 2.1.) We obtain the block diagram of figure 2.1. The block diagram of general secondorder discretetime state system "
x1 (k + 1) x2 (k + 1)
#
=
y(k) = is shown in figure 2.2.
"
a11 a12 a21 a22
h
i
c1 c2
#" "
x1 (k) x2 (k)
x1 (k) x2 (k)
#
#
+
"
b1 b2
+ d u(k)
#
u(k)
2.1. DISCRETETIME LINEAR SYSTEM OPERATORS
19
Figure 2.1: Vector block diagram for a linear discretetime system described by state space system dynamics
Figure 2.2: Block diagram for a state equationbased secondorder discretetime system
20CHAPTER 2. OPERATIONAL METHODS FOR DISCRETETIME LINEAR SYSTEMS
Inputoutput linear discretetime system models For systems that have only one input and one output, it is frequently convenient to work with the classical system inputoutput description in Eq. (1.16) consisting of a single nth order difference equation relating the output to the system input: an y(k + n) + an−1 y(k + n − 1) + . . . + a0 y(k) = bm u(k + m) + bm−1 u(k + m − 1) + . . . + b0 u(k)
(2.1)
The constant coefficients ai and bi are defined by the system parameters and for causal systems there holds: m ≤ n. Equation (2.1) may be written in operational form using polynomial operators as: h
i
h
i
an Z n + an−1 Z n−1 + . . . + a0 {y} = bm Z m + bm−1 Z m−1 + . . . + b0 {u} h
(2.2)
i−1
If the inverse operator an Z n + an−1 Z n−1 + . . . + a0 exists, it may be applied to each side to produce an explicit operational expression for the output variable: bm Z m + bm−1 Z m−1 + . . . + b0 y(k) = {u(k)} an Z n + an−1 Z n−1 + . . . + a0
(2.3)
The operational description of the system dynamics has been reduced to the form of a single linear operator, the dynamic transfer operator H{}, y(k) = H{u(k)}
(2.4)
where H{} =
bm Z m + bm−1 Z m−1 + . . . + b0 {} an Z n + an−1 Z n−1 + . . . + a0
(2.5)
In (1.14)(1.15), the state space description was given for a linear discretetime system: x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) + Du(k)
(2.6) (2.7)
The matrix transfer operator for this system can be computed by: H{} =
C adj(Z I − A) B + det[Z I − A] D det[Z I − A]
(2.8)
21
2.2. TRANSFORMATION OF DISCRETETIME SYSTEMS
2.2
Transformation of discretetime systems
State realization of a polynomial difference system The polynomial difference system (2.3), described by h
i
h
i
bm Z m + bm−1 Z m−1 + . . . + b0 y(k) = an Z n + an−1 Z n−1 + . . . + a0 u(k)
has a state realization x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) + Du(k) such that
A=
C=
h
0 0 .. .
1
··· .. . .. .
0 .. .
0 0 .. .
0 0 .. .
0 0 ··· 1 0 0 0 ··· 0 1 −a0 /an −a1 /an · · · −an−2 /an −an−1 /an
B=
0 0 .. . 0 0 1/an
b0 −a0 bn /an b1 −a1 bn /an · · · bn−2 −an−2 bn /an bn−1 −an−1 bn /an
i
D = bn /an
Polynomial realization of a state system The state difference system (2.6)(2.7), described by x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) + Du(k) has a polynomial realization h
h
det[Z I − A]] y(k) = C adj(ZI − A)B + det[ZI − A]D] u(k)
(2.9)
and a transfer operator H = C (ZI − A)−1 B + D C adj(ZI − A)B + det[ZI − A]D = det[Z I − A]
(2.10) (2.11)
22CHAPTER 2. OPERATIONAL METHODS FOR DISCRETETIME LINEAR SYSTEMS Example 7 (Polynomial to statespace) Consider the difference system 16 y(k + 3) − 20 y(k + 2) + 8 y(k + 1) − y(k) = 5 u(k + 2) − 7 u(k + 1) + 2 u(k) the state space realization is given by
y(k) =
h
x1 (k + 1) 0 1 0 x1 (k + 1) 0 0 1 x2 (k + 1) = 0 x2 (k + 1) + 0 u(k) x3 (k + 1) 1/16 −8/16 20/16 x3 (k + 1) 1/16
x1 (k) 2 −7 5 x2 (k) x3 (k) i
Example 8 (Statespace to polynomial) Consider the state system "
x1 (k + 1) x2 (k + 1)
#
=
y(k) =
" h
−0.5 1.5 −1 2 i
1 1
"
#"
x1 (k) x2 (k)
x1 (k) x2 (k) #
#
+
"
2 0
#
u(k)
+ 2 u(k)
The determinant of [Z I − A] is det[Z I − A] = Z 2 − 1.5 Z + 0.5 and the adjoint of [Z I − A] is "
Z + 0.5 −1.5 adj 1 Z −2
#
=
"
Z −2 1.5 −1 Z + 0.5
#
Further we compute C adj(ZI − A)B + det[ZI − A]D = " #" # h i Z −2 2 1.5 1 1 + (Z 2 − 1.5 Z + 0.5) 2 = −1 Z + 0.5 0 = (2 Z − 6) + (2Z 2 − 3 Z + 1) = 2 Z2 − Z − 5
So the related difference equation is given by equation (2.9): (Z 2 − 1.5 Z + 0.5) y(k) = (2 Z 2 − Z − 5) u(k) or y(k + 2) − 1.5 y(k + 1) + 0.5 y(k) = 2 u(k + 2) − u(k + 1) − 5 u(k)
Chapter 3 System Properties and Solution techniques In this chapter we discuss system properties and solution techniques for discretetime systems. (For the continuous case, see Rowell & Wormley [4], chapter 8).
3.1
Singularity input functions
In Rowell & Wormley [4], chapter 8, some singularity input functions for continuoustime systems were defined. In this section we redefine some of the functions for discretetime: The Unit Step Function: The discretetime unit step function us (k), k ∈ ZZ is defined as: us (k) =
(
0 1
for k < 0 for k ≥ 0
(3.1)
The Unit Impulse Function: The discretetime unit Impulse function δ(k), k ∈ ZZ is defined as: δ(k) =
(
1 0
for k = 0 for k = 6 0
(3.2)
Note that for discretetime the definition of a function δT (k) is not appropriate. The discretetime impulse function (or discretetime Dirac delta function) has the property that ∞ X
δ(m) = 1
m=−∞
23
24
CHAPTER 3. SYSTEM PROPERTIES AND SOLUTION TECHNIQUES
The unit ramp function: The discretetime unit ramp function ur (k), k ∈ ZZ is defined to be a linearly increasing function of time with a unity increment: ur (k) =
(
0 k
for k ≤ 0 for k > 0
(3.3)
Note that ∆{ur (k + 1)} = ur (k + 1) − ur (k) = us (k)
(3.4)
∆{us (k)} = us (k) − us (k − 1) = δ(k)
(3.5)
and
and also in the reverse direction: ∆−1 {δ(k)} =
k X
δ(k) = us (k)
(3.6)
m=0
and ∆−1 {us (k)} =
k X
us (k) = ur (k + 1)
(3.7)
m=0
3.2
Classical solution of linear difference equations
In this section we briefly review the classical method for solving a linear nthorder ordinary difference equation with constant coefficients, given by y(k + n) + an−1 y(k + n − 1) + . . . + a0 y(k) = bm u(k + m) + bm−1 u(k + m − 1) + . . . + b0 u(k)
(3.8)
where in general m ≤ n. We define the forcing function f (k) = bm u(k + m) + bm−1 u(k + m − 1) + . . . + b0 u(k)
(3.9)
which is known, because u(k) is known for all k ∈ ZZ. Eq. (3.8) now becomes y(k + n) + an−1 y(k + n − 1) + . . . + a0 y(k) = f (k)
(3.10)
The task is to find a unique function y(k) for k ≥ k0 that satisfies the difference equation (3.8) given the forcing function f (k) and a set of initial conditions y(k0), y(k0 +1), . . . , y(k0 + n − 1). The general solution to Eq. (3.8) may be derived as the sum of two solution components y(k) = yh (k) + yp (k)
(3.11)
where yh is the solution of the homogeneous equation (so for f (k) = 0) and yp (k) is a particular solution that satisfies (3.10) for the specific f (k) (but arbitrary initial conditions).
3.2. CLASSICAL SOLUTION OF LINEAR DIFFERENCE EQUATIONS
3.2.1
25
Homogeneous solution of the difference equation
The homogeneous difference equation (an 6= 0) is given by an y(k + n) + an−1 y(k + n − 1) + . . . + a0 y(k) = 0
(3.12)
The standard method of solving difference equations assumes there exists a solution of the form y(k) = C λk , where λ and C are both constants. Substitution in the homogeneous equation gives:
C an λn + an−1 λn−1 + . . . + a0 λk = 0
(3.13)
For any nontrivial solution, C is nonzero and λk is never zero, we require that an λn + an−1 λn−1 + . . . + a0 = 0
(3.14)
For an nthorder system with an 6= 0, there are n (complex) roots of the characteristic polynomial and n possible solution terms Ci λki (i = 1, . . . , n), each with its associated constant. The homogeneous solution will be the sum of all such terms: yh (k) = C1 λk1 + C2 λk2 + . . . + Cn λkn =
n X
Ci λki
(3.15) (3.16)
i=1
The n coefficients Ci are arbitrary constants and must be found from the n initial conditions. If the characteristic polynomial has repeated roots, that is λi = λj for i 6= j, there are not n linearly independent terms in the general homogeneous response (3.16). In general if a root λ of multiplicity m occurs, the m components in the general solution are C1 λk , C2 k λk , . . . , Cm k m−1 λm If one or more roots are equal to zero (λi = 0), the solution C λk = 0 has no meaning. In that case we add a term C δ(k − k0 ) to the general solution. If the root λ = 0 has multiplicity m, the m components in the general solution are C1 δ(k − k0 ) , C2 δ(k − k0 − 1) , . . . , Cm δ(k − k0 − m + 1)
3.2.2
Particular solution of the difference equation
Also for discretetime systems, the method of undetermined coefficients can be used (see [4], page 255257), where we use table 3.1 to find the particular solutions for some specific forcing functions.
26
CHAPTER 3. SYSTEM PROPERTIES AND SOLUTION TECHNIQUES Terms in u(k)
Assumed form for yp (k)
Test value
α α k n , (n = 1, 2, 3, . . .) α λk α cos(ωk) α sin(ωk)
β1 βn k n + βn−1 k n−1 + . . . + β1 k + β0 β λk β1 cos(ωk) + β2 sin(ωk) β1 cos(ωk) + β2 sin(ωk)
0 0 λ jω jω
Figure 3.1: Definition of Particular yp (k) using the method of undetermined coefficients
Example 9 (Solution of a difference equation) Consider the system described by the difference equation 2 y(k + 3) + y(k + 2) = 7 u(k + 1) − u(k) Find the system response to a ramp input u(k) = k and initial conditions given by y(0) = 2, y(1) = −1 and y(2) = 2. Solution: homogeneous solution: The characteristic equation is: 2 λ3 + λ2 = 0 which has a root λ1 = −0.5 and a double root λ2,3 = 0. The general solution of the homogeneous equation is therefore yh (k) = C1 (−0.5)k + C2 δ(k) + C3 δ(k − 1) particular solution: From table 3.1 we find that for u(k) = k, the particular solution is selected: yp = β1 k + β0 Testing gives: 2 (β1 (k + 3) + β0 ) + (β1 (k + 2) + β0 ) = 7 (k + 1) − k = 6 k + 7 We find β1 = 2 and β0 = −3 and so the particular solution is yp (k) = 2 k − 3 complete solution: The complete solution will have the form: y(k) = yh (k) + yp (k) = C1 (−0.5)k + C2 δ(k) + C3 δ(k − 1) + 2 k − 3
27
3.3. DISCRETETIME CONVOLUTION Now the initial conditions are evaluated in k = 0, k = 1 and k = 2. C1 (−0.5)0 + C2 δ(0) + C3 δ(0 − 1) + 2 0 − 3 C1 + C2 − 3 = 2 C1 (−0.5)1 + C2 δ(1) + C3 δ(1 − 1) + 2 1 − 3 −C1 0.5 + C3 − 1 = −1 C1 (−0.5)2 + C2 δ(2) + C3 δ(2 − 1) + 2 2 − 3 C1 0.25 + 1 = 2
y(0) = = y(1) = = y(2) = =
We find a solution for C1 = 4, C2 = 1, C3 = 2 and so the final solution becomes: y(k) = 4 (−0.5)k + δ(k) + 2 δ(k − 1) + 2 k − 3
3.3
Discretetime Convolution
In this section we derive the computational form of the discretetime system H{u(k)}, defined in [4], chapter 7, and section 2.1 of these lecture notes, that is based on a system’s response to an impulse input. We assume that the system is initially at rest, that is, all initial conditions are zero at time t = 0, and examine the discretetime domain forced response y(k), k ∈ ZZ to a discretetime waveform u(k). To start with, we assume that the system response to δ(k) is a known function and is designated h(k) as shown in figure 3.2.
δ(k − n) 6
y(k) 6
s
s s s
s
s
s
s
s
s
n
s
k
s s s
s
s
n

k
Figure 3.2: System response to a delayed unit pulse Then if the system is linear and timeinvariant, the response to a delayed unit pulse, occurring at time n is simply a delayed version of the pulse response: u(k) = δ(k − n)
gives
y(k) = h(k − n)
(3.17)
Multiplication with a constant α gives: u(k) = α δ(k − n)
gives
y(k) = α h(k − n)
(3.18)
28
CHAPTER 3. SYSTEM PROPERTIES AND SOLUTION TECHNIQUES
The input signal u(k) may be considered to be the sum of nonoverlapping delayed pulses pn (k): ∞ X
u(k) =
pn (k)
(3.19)
n=−∞
where pn (k) =
(
u(k) for k = n 0 for k 6= n
(3.20)
Each component pn (k) may be written in terms of a delayed unit pulse δ(k), that is pn (k) = u(n)δ(k − n) From equation (3.18) it follows that for a delayed version of the pulse response, multiplied by a constant u(n) we obtain the output y(k) = u(n) h(k − n). For an input u(k) =
∞ X
pn (k)
n=−∞
we can use the principle of superposition and the output can be written as a sum of all responses to pn (k), so: y(k) =
∞ X
u(n) h(k − n)
(3.21)
n=−∞
This sum is denoted as the convolution sum for discretetime systems. For physical systems, the pulse response h(k) is zero for time k < 0, and future components of the input do not contribute to the sum. So the convolution sum becomes: y(k) =
k X
u(n) h(k − n)
(3.22)
n=−∞
The convolution operation is often denoted by the symbol ∗ : y(k) = u(k) ∗ h(k) =
k X
u(n) h(k − n)
(3.23)
n=−∞
Example 10 (Discretetime convolution) The firstorder exponential smoother of example 6 is described by the difference equation y(k + 1) = a y(k) + (1 − a) u(k + 1) By repeated substitution it easily follows that if y(k) = 0 for k < n0 , the output of the system is given by y(k) = (1 − a)
k X
n=n0
ak−n u(n)
k ≥ n0
3.3. DISCRETETIME CONVOLUTION
29
assuming that the input u(k) is such that the sum converges, and let n0 approach to −∞. Then the response of the system takes the form y(k) = (1 − a)
k X
ak−n u(n)
k ∈ ZZ
n=−∞
Define the function h such that h(k) =
(
0 for k < 0 k (1 − a)a for k ≥ 0
= (1 − a)ak us (k) we see that on the infinite time axis the system may be represented as the convolution sum y(k) =
∞ X
n=−∞
h(k − n) u(k)
30
CHAPTER 3. SYSTEM PROPERTIES AND SOLUTION TECHNIQUES
Chapter 4 Solution of discretetime linear state equations In this chapter we will discuss the general solution of discretetime linear state equations. (For the continuoustime case, see Rowell & Wormley [4], chapter 10).
4.1
State variable responses of discretetime linear systems
In this chapter we examine the responses of linear timeinvariant discretetime models in the standard state equation form x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) + Du(k)
(4.1) (4.2)
The solution proceeds in two steps: First the statevariable response x(k) is determined by solving the set of firstorder state equations Eq. (4.1), and then the state response is substituted into the algebraic output equations, Eq. (4.2) in order to compute y(k).
4.1.1
Homogeneous State responses
The state variable response of a system described by Eq. (4.1) with zero input and an arbitrary set of initial conditions x(0) is the solution of the set of n homogeneous firstorder difference equations x(k + 1) = Ax(k)
(4.3)
We can find the solution by successive substitution. Note that x(1) = A x(0) and x(2) = A x(1), so by substitution we obtain x(2) = A2 x(0). Proceeding in this way we derive: x(k) = Ak x(0)
(4.4) 31
32
CHAPTER 4. SOLUTION OF DISCRETETIME LINEAR STATE EQUATIONS
The solution is often written as: x(k) = Φ(k) x(0)
(4.5)
where Φ(k) = Ak is defined to be the state transition matrix. Example 11 (Firstorder exponential smoother) In example 6, the firstorder exponential smoother was presented. The state difference representation and output equations are x(k + 1) = a x(k) + a(1 − a) u(k) y(k) = x(k) + (1 − a) u(k),
k ∈ ZZ
The homogeneous state difference equation is given by x(k + 1) = ax(k) It follows that the 1 × 1 state transition matrix is given by Φ(k) = A A · · · A = a · a · · · a = ak , k > 0, k ∈ ZZ. Example 12 (Secondorder exponential smoother) Consider now the secondorder smoother described by the difference equation y(k + 2) − a1 y(k + 1) − a0 y(k) = b2 u(k + 2) + b1 u(k + 1),
k ∈ ZZ
Since the order of the system is N = 2, the dimension of the state is also 2. According to the transformation formulae in section 2.2, the system may be represented in state form as x(k + 1) = y(k) =
"
a1 1 a0 0
h
1 0
i
#
x(k) +
"
b1 + a1 b2 a0 b2
+ b2 u(k),
#
u(k),
k ∈ ZZ
The homogeneous state difference equation is x(k + 1) =
"
a1 1 a0 0
#
x(k) = Ax(k)
(4.6)
Then the 2 × 2 state transition matrix is given by Φ(k) = A · A · · · A = Ak For instance, let us take numerical values a0 = 0, and a1 = 21 , then A=
"
1 2
1 0 0
#
;
(4.7)
33
4.1. STATE VARIABLE RESPONSES OF DISCRETETIME LINEAR SYSTEMS and it can be seen that Φ(k) =
"
( 12 )k 2( 12 )k 0 0
#
;
k ≥ 0, k ∈ ZZ
(4.8)
Consider homogeneous system (4.6) with system matrix (4.7) and initial condition x(0) =
"
16 4
#
Then for k = 3 we find: x(3) = Φ(3) x(0) = Φ(k) =
4.1.2
"
( 21 )3 2( 12 )3 0 0
#"
16 4
#
=
"
1/8 1/4 0 0
#"
16 4
#
=
"
#
3 (4.9) 0
The forced response of discretetime linear systems
The solution of the inhomogeneous discretetime state equation x(k + 1) = Ax(k) + Bu(k)
(4.10) (4.11)
follows easily by induction x(k + 2) = Ax(k + 1) + Bu(k + 1)
= A Ax(k) + Bu(k) + Bu(k + 1) = A2 x(k) + ABu(k) + Bu(k + 1) Proceeding in this way we derive x(k) = Ak x(0) +
k−1 X
Ak−m−1 B u(m)
(4.12)
m=0
for k ≥ 1. For k ≤ 0, the summation does not contribute to the solution, and we therefore multiply this term by us (k − 1), which is zero for k ≤ 0. We obtain: x(k) = Ak x(0) +
k−1 X
Ak−m−1 B u(m) us (k − 1)
m=0
Example 13 (Secondorder exponential smoother) Consider again the second order exponential smoother of example 12. Since A(k) =
"
a1 1 a0 0
#
;
B(k) =
"
b1 + a1 b2 a0 b2
#
yields x(k) = Ak x(0) +
k−1 X
m=0
Ak−m−1 B u(m) us (k − 1)
(4.13)
34
CHAPTER 4. SOLUTION OF DISCRETETIME LINEAR STATE EQUATIONS
Let us adopt the numerical values a0 = 0, a1 = 12 , b1 = 1 and b2 = 0 resulting in A as in (4.7) and B, C, and D as follows B=
"
1 0
#
;
C=
h
i
1
0 ;
D=0
(4.14)
With (4.8) we may express the forced response of this system as x(k) =
4.1.3
"
( 21 )k 2( 12 )k 0 0
#
x(0) +
k−1 X
m=0
"
( 21 )k−m−1 0
#
u(m) us (k − 1)
The system output response of discretetime linear systems
The output response of a discretetime linear system is easily derived by substitution of (4.13) in (4.2): y(k) = Cx(k) + Du(k) k
= C A x(0) + = CAk x(0) +
k−1 X
k−m−1
A
!
B u(m) us (k − 1) + Du(k)
m=0 k−1 X
CAk−m−1 B u(m) us(k − 1) + Du(k)
m=0
So the forced output response is given by y(k) = CAk x(0) +
k−1 X
CAk−m−1 B u(m) us(k − 1) + Du(k)
m=0
Example 14 (Secondorder exponential smoother) Based on example 13, we may obtain for the secondorder exponential smoother the following output response k−1 X 1 1 y(k) = ( )k x1 (0) + ( )k−m−1 u(m) us(k − 1) 2 m=0 2
k−1 X 1 1 = ( )k x1 (0) + ( )m u(k − m − 1) us(k − 1) 2 m=0 2
4.1.4
The discretetime transition matrix
The discretetime transition matrix Φ(k) has the following properties: 1. Φ(0) = I. 2. Φ(−k) = Φ−1 (k).
4.1. STATE VARIABLE RESPONSES OF DISCRETETIME LINEAR SYSTEMS
35
3. Φ(k1 )Φ(k2 ) = Φ(k1 + k2 ), and so x(k2 ) = Φ(k2 )x(0) = Φ(k2 )Φ(−k1 )x(k1 ) = Φ(k2 − k1 )x(k1 ) or x(k2 ) = Φ(k2 − k1 )x(k1 ) 4. If A is a diagonal matrix, then Ak is also a diagonal matrix and each element is the kth power of the corresponding diagonal element of the A matrix, that is, akii .
4.1.5
System Eigenvalues and Eigenvectors
Let A be the system matrix of system (4.1)(4.2). The values λi satisfying the equation λi mi = A mi
for mi 6= 0
(4.15)
are known as the eigenvalues, or characteristic values, of A. The corresponding column vectors m are defined as eigenvectors, or characteristic vectors. Equation (4.15) can be rewritten as: (λi I − A) mi = 0
(4.16)
The condition for a nontrivial solution of such a set of linear equations is that det(λi I − A) = 0
(4.17)
which is defined as the characteristic polynomial of the A matrix. Eq. (4.17) may be written as λn + an−1 λn−1 + . . . + a1 λ + a0 = 0
(4.18)
or, in factored form in terms of its roots λ1 , . . . , λn , (λ − λ1 )(λ − λ2 ) · · · (λ − λn ) = 0 Define M= and
Λ=
h
m1 m2 · · · mn
λ1
0
0 .. .
λ2
0
···
··· ..
.
0 .. .
λn
i
(4.19)
36
CHAPTER 4. SOLUTION OF DISCRETETIME LINEAR STATE EQUATIONS
then
k
Λ =
λk1
0
0 .. .
λk2
0
···
··· ..
.
0 .. .
λkn
and Φ(k) = Ak = (M Λ M −1 )k = M Λk M −1
4.1.6
Stability of discretetime linear systems
For asymptotic stability, the homogeneous response of the state vector x(k) returns to the origin for arbitrary initial conditions x(0) at time k → ∞, or lim x(k) = lim Φ(k) x(0) = lim M Λk M −1 x(0) = 0
k→∞
k→∞
k→∞
for any x(0). All the elements are a linear combination of the modal components λki , therefore, the stability of a system response depends on all components decaying to zero with time. If λ > 1, the component will grow exponentially with time and the sum is by definition unstable. The requirements for system stability may therefore be summarized: A linear discretetime system described by the state equation x(k + 1) = A x(k) + b u(k) is asymptotically stable if and only if all eigenvalues have magnitude smaller than one. Three other separate conditions should be considered: 1. If one or more eigenvalues, or pair of conjugate eigenvalues, has a magnitude larger than one, there is at least one corresponding modal component that increases exponentially without bound from any initial condition, violating the definition of stability. 2. Any pair of conjugate eigenvalues that have magnitude equal to one, λi,i+1 = e±jω , generates an undamped oscillatory component in the state response. The magnitude of the homogeneous system response neither decays nor grows but continues to oscillate for all time at a frequency ω. Such a system is defined to be marginally stable. 3. An eigenvalue λ = 1 generates a model exponent λk = 1k = 1 that is a constant. The system response neither decays or grows, and again the system is defined to be marginally stable.
4.1. STATE VARIABLE RESPONSES OF DISCRETETIME LINEAR SYSTEMS
37
Example 15 (Secondorder exponential smoother) Consider the secondorder exponential smoother of example 12. With the numerical values a0 = 0, a1 = 12 , b1 = 1 and b2 = 0, the system matrices are given by 4.7 and 4.14. The characteristic polynomial of the matrix A is det(λI − A) = det
"
λ− 0
1 2
−1 λ
#
1 = λ(λ − ) 2
(4.20)
As a result, the eigenvalues of A are λ1 = 0 and λ2 = 12 . Thus since all its eigenvalues are less than one, this discretetime system is stable.
4.1.7
Transformation of state variables
Consider the discretetime system x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) + Du(k)
(4.21) (4.22)
and let consider the eigenvalue decomposition of the A matrix A = M Λ M −1 Then this system can be transformed to the diagonal form x′ (k + 1) = A′ x′ (k) + B ′ u(k) y(k) = C ′ x′ (k) + D ′ u(k)
(4.23) (4.24)
by choosing x′ (k) A′ B′ C′ D′
= = = = =
M −1 x(k) M −1 AM = Λ M −1 B CM D
Example 16 (Secondorder exponential smoother) Consider the secondorder exponential smoother of example 12 and example 15. It can be found that the corresponding eigenvectors are v1 =
"
1 − 12
#
;
v2 =
"
1 0
#
;
(4.25)
It follows that the modal transformation matrix V , its inverse V −1 and the diagonal matrix Λ are V =
"
1 1 − 12 0
#
;
V
−1
=
"
0 −2 1 2
#
;
Λ=
"
0 0 0 12
#
;
(4.26)
38
CHAPTER 4. SOLUTION OF DISCRETETIME LINEAR STATE EQUATIONS
Thus, after modal transformation the system is represented as in Eqs. (4.23) and (4.24) where ′
A =V B ′ = V −1 B =
4.2
"
0 1
#
−1
AV = Λ =
C ′ = CV =
;
"
0 0 0 21
h
1
#
;
i
1 ;
The response of linear DT systems to the impulse response
Let u(k) = δ(k) where δ(k) is the impulse response, defined in Eq. (3.2). Following Eq. (4.15) we find y(k) = CAk x(0) +
k−1 X
CAk−m−1 B δ(m) us (k − 1) + Dδ(k)
(4.27)
m=0
Note that the term k−1 X
CAk−m−1 B δ(m) us (k − 1)
m=0
only gives a contribution for m = 0, so: k−1 X
CAk−m−1 B δ(m) us (k − 1) = CAk−1 B us (k − 1)
m=0
Eq. (4.27) now becomes y(k) = CAk x(0) + CAk−1 B us (k − 1) + Dδ(k)
(4.28)
which is the impulse response of a linear DT system. Remark: Note that a state transformation does not influence the inputoutput behaviour of the system. We therefore may use the diagonal form of the state space description of the system. The impulse response becomes: y(k) = C ′ A′k x(0) + C ′ A′k−1 B ′ us (k − 1) + D ′ δ(k) = CMΛk x(0) + CMΛk−1 M −1 B us (k − 1) + Dδ(k) which can easily be computed.
4.2. THE RESPONSE OF LINEAR DT SYSTEMS TO THE IMPULSE RESPONSE 39 Example 17 (Secondorder exponential smoother) Consider again the second order exponential smoother of example 12 and example 15. The state transition matrix of the transformed system is ′
k
Φ (k) = Λ =
"
0 0 0 ( 21 )k
#
; k ≥ k0
(4.29)
It follows for the impulse response of the transformed system, which is also the impulse response of the untransformed system, h(k) = C ′ Λk−1 B ′ us (k − 1) + Dδ(k) " #" # h i ∆(k − 1) 0 0 11 = 1 0 ( 12 )k−1 1 = ( )k−1us (k − 1). 2
40
CHAPTER 4. SOLUTION OF DISCRETETIME LINEAR STATE EQUATIONS
Chapter 5 The Discretetime transfer function In this chapter we will discuss the transfer function of discretetime linear state equations. (For the continuoustime case, see Rowell & Wormley [4], chapter 12).
5.1
Introduction
The concept of the transfer function is developed here in terms of the particular solution component of the total system response when the system is excited by a given exponential input waveform of the form u(k) = U(z)z k
(5.1)
where z = ρ ej ψ is a complex variable with magnitude ρ and phase ψ and the amplitude U(z) is in general complex. Then the input u(k) = U(z) ρk ej ψ k = U(z) ρk ( cos ψk + j sin ψk )
(5.2)
is itself complex and represent a broad class of input functions of engineering interest including growing and decaying real exponential waveforms as well as growing and decaying sinusoidal waveforms.
5.2
singleinput singleoutput systems
Consider a system, described by a singleinput singleoutput differential equation an y(k + n) + an−1 y(k + n − 1) + . . . + a0 y(k) = bm u(k + m) + bm−1 u(k + m − 1) + . . . + b0 u(k)
(5.3)
where m ≤ n and all coefficients are real constants. For an exponential input signal of the form u(k) = U(z)z k
(5.4) 41
42
CHAPTER 5. THE DISCRETETIME TRANSFER FUNCTION
we may assume that the particular solution yp (k) is also exponential in form, that is, yp (k) = Y (z)z k
(5.5)
where y(z) is a complex amplitude to be determined. Substitution of (5.4), (5.5) in (5.3) gives an Y (z)z k+n + an−1 Y (z)z k+n−1 + . . . + a0 Y (z)z k = bm U(z)z k+m + bm−1 U(z)z k+m−1 + . . . + b0 U(z)z k
(5.6)
or
an z n + an−1 z n−1 + . . . + a0 ) Y (z)z k = (bm z m + bm−1 z m−1 + . . . + b0 ) U(z)z k
(5.7)
The transfer function H(z) is defined to be the ration of the response amplitude Y (z) to the input amplitude U(z) and is Y (z) bm z m + bm−1 z m−1 + . . . + b0 = H(z) = U(z) an z n + an−1 z n−1 + . . . + a0
(5.8)
The transfer function Eq. (5.8) is an algebraic rational function of the variable z.
Example 18 (Transfer function of discretetime system) Consider a discretetime system described by the third order difference equation 100 y(k + 3) − 180 y(k + 2) + 121 y(k + 1) − 41 y(k) = 100 u(k + 3) − 10 u(k + 2) + 48 u(k + 1) − 34 u(k) The transfer function is now given by H(z) =
5.3
100 z 3 − 10 z 2 + 48 z − 34 100 z 3 − 180 z 2 + 121 z − 41
(5.9)
relationship to the transfer function
The use of Z{} as the difference operator in chapter 2 and z as the exponent in the exponential input function creates a similarity in appearance. When consideration is considered to linear timeinvariant systems, these two system representations are frequently used interchangeably. The similarity results directly from the difference operator relationship for an exponential waveform Z n {U(z)z k } ≡ U(z) z k+n
(5.10)
43
5.4. SYSTEM POLES AND ZEROS
5.4
System poles and zeros
It is often convenient to factor the polynomials in the numerator and denominator of Eq. (5.8) and to write the transfer function in terms of those factors: H(z) =
N(z) (z − w1 )(z − w2 ) . . . (z − wm ) =K D(z) (z − p1 )(z − p2 ) . . . (z − pn )
(5.11)
where the numerator and denominator polynomials, N(z) and D(z), have real coefficients defined by the systems’s difference equation and K = bm /an . As written in Eq. (5.11), the wi ’s are the roots of the equation N(z) = 0 and are defined to be the system zeros, and the pi ’s are the roots of the equation D(z) = 0 and are defined to be the system poles. Note that lim H(z) = 0
and
z→wi
lim H(z) = ∞
z→pi
Example 19 (Secondorder exponential smoother) Consider difference equation of the secondorder exponential smoother as given in example 12. y(k + 2) − a1 y(k + 1) − a0 y(k) = b2 u(k + 2) + b1 (k + 1),
k ∈ ZZ
The corresponding transfer function is found using Eq. (5.8): H(z) =
Y (z) b2 z 2 + b1 z = 2 U(z) z + a1 z + a0
The zeros are given by the roots of N(z) = b2 z 2 + b1 z = 0 so for b2 6= 0 we find w1 = 0
w2 = −b1 /b2
and for b2 = 0 we only have one zero w1 = 0. The poles are given by the roots of D(z) = z 2 + a1 z + a0 = 0 so we find p1 = −a1 /2 +
q
a21 /4 − a0
p2 = −a1 /2 −
q
a21 /4 − a0
For the setting a0 = 0.64, a1 = −1.6, b1 = −0.5 and b2 = 1 we find: w1 = 0
,
w2 = 0.5
,
p1 = p2 = 0.8
44
CHAPTER 5. THE DISCRETETIME TRANSFER FUNCTION
Figure 5.1: Polezero plot of a third order discretetime system
Example 20 (Poles and zeros of thirdorder system) Consider the thirdorder system in example 18, with transfer function H(z) =
100 z 3 − 10 z 2 + 48 z − 34 100 z 3 − 180 z 2 + 121 z − 41
(5.12)
The poles of the system are finding the roots of the equation D(z) = 100 z 3 − 180 z 2 + 121 z − 41 = 100(z − 1)(z − 0.4 − j 0.5)(z − 0.4 + j 0.5) = 0 We find that the poles are equal to p1 = 1, p2 = 0.4 + j 0.5 and p3 = 0.4 − j 0.5. The zeros of the system are finding the roots of the equation N(z) = 100 z 3 − 10 z 2 + 48 z − 34 = 100(z − 0.5)(z + 0.2 − j 0.8)(z + 0.2 + j 0.8) = 0 We find that the zeros are equal to p1 = 0.5, p2 = −0.2 + j 0.8 and p3 = −0.2 − j 0.8. The polezero plot of this system is given in Fig. 5.1
5.4.1
System poles and the homogeneous response
Also in the discretetime case there holds The transfer function poles are the roots of the characteristic equation and also the eigenvalues of the system A matrix (as discussed in the previous chapter).
5.4. SYSTEM POLES AND ZEROS
45
So let the transfer function denominator polynomial be given by D(z) = an z n + an−1 z n−1 + . . . + a0 and let the roots of D(z) = 0 be given by pi , i = 1, . . . , n. Then the homogeneous response of the system is given by yh (k) =
n X
Ci pki
i=1
The location of the poles in the zplane therefore determines the n components in the homogeneous response as follows: 1. A real pole inside the unit circle (p < 1) defines an exponentially decaying component C pk in the homogeneous response. The rate is determined by the pole location; poles further away from the origin corresponds to components that decay rapidly; poles close to the origin correspond to slowly decaying components. Negative real poles alternate in sign. 2. A pole in pi = 1 defines a component that is constant in amplitude and defined by the initial conditions. 3. A real pole outside the unit circle (p > 1) corresponds to an exponentially increasing component C pk in the homogeneous response, thus defining the system to be unstable. 4. A complex conjugate pole pair p = ρ e±j ψ inside the unit circle combine to generate a response component that is decaying sinusoid of the form A ρk sin(ψ k + φ) where A and φ are determined by the initial conditions. The rate of decaying is specified by ρ; The frequency of oscillation is determined by ψ. 5. A complex conjugate pole pair on the unit circle p = e±j ψ generates an oscillatory component of the form A sin(ψ k + φ) where A and φ are determined by the initial conditions. 6. A complex pole pair outside the unit circle generates a exponentially increasing oscillatory component.
5.4.2
System stability
A nth order linear discretetime system is asymptotically stable only if all components in the homogeneous response from a finite set of initial conditions decay to zero as time increases, or lim
k→∞
n X
Ci pki = 0
i=1
where pi are the system poles. If any poles has magnitude larger than one, there is a component in the output without bound, causing the system to be unstable.
46
CHAPTER 5. THE DISCRETETIME TRANSFER FUNCTION
p = −0.5
p = 0.2 + j 0.8 J ] J J J
J J J
Q k Q Q
p = −0.9
p = −1
OCC
J J
p=1
p = ej 0.4 6
C C
J
p=0
?
p = 1.2
HH H H j H
p = 0.5
p = 0.9
Figure 5.2: The specification of the form of components of the homogeneous response from the system pole locations on the polezero plot
47
5.4. SYSTEM POLES AND ZEROS In order for a linear discretetime system to be stable, all its poles must have a magnitude strictly smaller than one, that is they must all lie inside the unit circle. An “unstable” pole, lying outside the unit circle, generates a component in the system homogeneous response that increases without bound from any finite initial conditions. A system having one or more poles lying on the unit circle has nondecaying oscillatory components in the homogeneous response and is defined to be marginally stable.
5.4.3
State space formulated systems
Consider the nth order linear system describe by the set of n state equations and one output equation: x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) + Du(k)
(5.13) (5.14)
The transfer function is given by Y (z) H(z) = U(z) = [C (z I − A)−1 B + D] C adj(z I − A)−1 B + D det(z I − A) = det(z I − A) Example 21 (Transfer function of a thirdorder system) Consider the thirdorder system in example 18. The statespace representation is found using the results of section 2.2: 0 1 0 0 0 1 A= 0 B= 0 0.41 −1.21 1.8 0.01 C=
h
7 −73 170
i
D=1
We compute the adjoint of (z I − A): z 2 − 1.8 z + 1.21 0.41 0.41 z z − 1.8 z 2 − 1.8 z −1.21 z + 0.41 adj(z I − A) = 2 1 z z
and the determinant of (z I − A):
det(z I − A) = z 3 − 1.8 z 2 + 1.21 z − 0.41 and so: C adj(z I − A)−1 B + D det(z I − A) det(z I − A) 3 2 z − 0.1 z + 0.48 z − 0.34 = 3 z − 1.80 z 2 + 1.21 z − 0.41
H(z) =
48
CHAPTER 5. THE DISCRETETIME TRANSFER FUNCTION
Bibliography [1] H. Freeman. Discretetime systems. John Wiley & Sons, New York, USA, 1965. [2] H. Kwakernaak and R. Sivan. Modern Signals and Systems. Prentice Hall, New Jersey, USA, 1991. [3] C. Phillips and J.M. Parr. Signals, Systems, and Tranforms. Prentice Hall, New Jersey, USA, 1999. [4] D. Rowell and D.N. Wormley. System Dynamics, an introduction. Prentice Hall, New Jersey, USA, 1997.
49