Continuous & Discreate Control Systems

Continuous and Discrete Control Systems John Dorsey

Chapter 1

Preliminaries 1.1

Overview

The opening chapter of the book doesn't require the use of Matlab, but we can take a preliminary look at Matlab and some of its basic commands. Matlab is, to me, just a versatile, on-line calculator, with some some limited graphical capabilities that are not user friendly. I use it in preference to a hand held calculator, but if I want to do some serious number crunching I turn to a real programming language. Matlab is an interpreted language. That is, it runs along executing each instruction as it comes to it. In normal use there is no compilation of the code. The latest releases do include a compiler, however. In normal use, without compilation, `for' loops and `while' loops don't run extremely fast. However, if we put a semicolon behind each command inside such a loop things speed up considerably. It is then possible to run simple iterative routines in a couple of seconds on the current generation of desktop computers. This is plenty fast enough for the kinds of programs we generally write with Matlab. On average, the specialized programs you nd in subsequent chapters have fewer than thirty statements. Thus, Matlab serves as a `quick and dirty' super calculator. 1.2

The basic Commands

To illustrate the basic commands I will use what I refer to as `in line Matlab dialogues.' When I type something, Matlab then responds. Hence the dialogue. To multiply to numbers together we type EDU>2 * 4 ans = 8 EDU>

We can assign names to numbers and then compute with them as follows.

2

CHAPTER 1.

EDU>a = 2 a = 2 EDU>b = 4 b = 4 EDU>c = a*b c = 8 EDU>

We can form arrays in the same way. EDU>a = [1 2] a = 1

2

EDU>b = [ 1 ; 2 ] b = 1 2 EDU>c = a*b c = 5 EDU>

Note that the semicolon separates rows in forming the column vector b. For the two arrays a and b, if we can also type

PRELIMINARIES

1.2.

THE BASIC COMMANDS

3

EDU>c = b*a c = 1 2

2 4

EDU>

Matlab also contains all the usual mathematical functions. For instance, EDU>y = cos(pi) y = -1 EDU>x = atan(0.5) x = 4.636476090008061e-01 EDU>

Note that the angle x is in radians. If we want to plot the function cos t over two periods we can do so using program chap1cda, repeated below for convenience. t1 = 0 t2 = 4*pi x = linspace(t1,t2) y = cos(x) plot(x,y) print -deps cosx.eps

The resulting plot is shown in Figure 1.1 The Matlab command linspace creates 100 equally spaced points between 0 and 4 . If we want more, or less, than 100 points we can type linspace(x,y,npts),

where npts is the number of points desired. We could also generate this plot using the program chap1cdb, repeated below for convenience. t1 = 0 t2 = 4*pi npts = 200 t = t1

4

CHAPTER 1.

PRELIMINARIES

1

0.8

0.6

0.4

0.2

0

-0.2

-0.4

-0.6

-0.8

-1

0

2

4

6

8

10

12

14

Figure 1.1: Plot of cos x Generated by chap1cda.m dt = (t2 - t1)/npts for i = 1:npts y(i) = cos(t); x(i) = t; t = t + dt; end plot(x,y) print -deps cosx.eps

The commands inside the `for' loop are terminated with a semicolon. This suppresses Matlab's print function and drastically speeds up the execution of the program. We have given an idea of how to use Matlab. To see all the possible Matlab commands pull down the Help menu in the Matlab command window and select Help Window. To see all the commands in the Control Toolbox type help control

In the Matlab command window. Matlab will respond with EDU>help control Control System Toolbox. Version 4.0 Student Edition 31-Dec-1996 Creation of LTI models. ss - Create a state-space model.

1.2.

THE BASIC COMMANDS

zpk tf dss filt set ltiprops

-

Create a zero/pole/gain model. Create a transfer function model. Specify a descriptor state-space model. Specify a digital filter. Set/modify properties of LTI models. Detailed help for available LTI properties.

Data extraction. ssdata - Extract state-space matrices. zpkdata - Extract zero/pole/gain data. tfdata - Extract numerator(s) and denominator(s). dssdata - Descriptor version of SSDATA. get - Access values of LTI model properties. Model characteristics. class - Model type ('ss', 'zpk', or 'tf'). size - Input/output dimensions. isempty - True for empty LTI models. isct - True for continuous-time models. isdt - True for discrete-time models. isproper - True for proper LTI models. issiso - True for single-input/single-output systems. isa - Test if LTI model is of given type. Conversions. ss zpk tf c2d d2c d2d

-

Conversion to state space. Conversion to zero/pole/gain. Conversion to transfer function. Continuous to discrete conversion. Discrete to continuous conversion. Resample discrete system or add input delay(s).

Overloaded arithmetic operations. + and - Add and subtract LTI systems (parallel connection). * - Multiplication of LTI systems (series connection). \ - Left divide -- sys1\sys2 means inv(sys1)*sys2. / - Right divide -- sys1/sys2 means sys1*inv(sys2). ' - Pertransposition. .' - Transposition of input/output map. [..] - Horizontal/vertical concatenation of LTI systems. inv - Inverse of an LTI system. Model dynamics. pole, eig - System poles. tzero - System transmission zeros.

5

6

CHAPTER 1.

pzmap dcgain norm covar damp esort dsort pade

-

PRELIMINARIES

Pole-zero map. D.C. (low frequency) gain. Norms of LTI systems. Covariance of response to white noise. Natural frequency and damping of system poles. Sort continuous poles by real part. Sort discrete poles by magnitude. Pade approximation of time delays.

State-space models. rss,drss - Random stable state-space models. ss2ss - State coordinate transformation. canon - State-space canonical forms. ctrb, obsv - Controllability and observability matrices. gram - Controllability and observability gramians. ssbal - Diagonal balancing of state-space realizations. balreal - Gramian-based input/output balancing. modred - Model state reduction. minreal - Minimal realization and pole/zero cancellation. augstate - Augment output by appending states. Time response. step impulse initial lsim ltiview gensig stepfun -

Step response. Impulse response. Response of state-space system with given initial state. Response to arbitrary inputs. Response analysis GUI. Generate input signal for LSIM. Generate unit-step input.

Frequency response. bode - Bode plot of the frequency response. sigma - Singular value frequency plot. nyquist - Nyquist plot. nichols - Nichols chart. ltiview - Response analysis GUI. evalfr - Evaluate frequency response at given frequency. freqresp - Frequency response over a frequency grid. margin - Gain and phase margins. System interconnections. append - Group LTI systems by appending inputs and outputs. parallel - Generalized parallel connection (see also overloaded +). series - Generalized series connection (see also overloaded *). feedback - Feedback connection of two systems.

1.2.

THE BASIC COMMANDS

star connect

7

- Redheffer star product (LFT interconnections). - Derive state-space model from block diagram description.

Classical design tools. rlocus - Evans root locus. rlocfind - Interactive root locus gain determination. acker - SISO pole placement. place - MIMO pole placement. estim - Form estimator given estimator gain. reg - Form regulator given state-feedback and estimator gains. LQG design tools. lqr,dlqr - Linear-quadratic (LQ) state-feedback regulator. lqry - LQ regulator with output weighting. lqrd - Discrete LQ regulator for continuous plant. kalman - Kalman estimator. kalmd - Discrete Kalman estimator for continuous plant. lqgreg - Form LQG regulator given LQ gain and Kalman estimator. Matrix equation lyap dlyap care dare -

solvers. Solve continuous Lyapunov equations. Solve discrete Lyapunov equations. Solve continuous algebraic Riccati equations. Solve discrete algebraic Riccati equations.

Demonstrations. ctrldemo - Introduction to the Control System Toolbox. jetdemo - Classical design of jet transport yaw damper. diskdemo - Digital design of hard-disk-drive controller. milldemo - SISO and MIMO LQG control of steel rolling mill. kalmdemo - Kalman filter design and simulation. EDU>

All the commands that we have discussed, plus many, many more are listed. Note that there are also ve demonstration programs that can be run.

Chapter 2

The Laplace Transform 2.1

Important Commands

The Matlab commands of most use in this chapter are 1. conv 2. roots 3. residue The conv command is a quick means of multiplying together polynomials. For instance, suppose in the Matlab command window we type p = [1 2 1]

This creates a row vector that stands for the polynomial ( ) = s2 + 2s + 1:

p s

If we next create the polynomial ( ) = s3 + 2s2 + 3s + 1

q s

by typing q = [1 2 3 1]

we can then nd the multiplication of these two polynomials by typing r = conv(p,q)

What this would like if you were in the Matlab command window is: 1

2

CHAPTER 2.

THE LAPLACE TRANSFORM

EDU>p = [1 2 1] p = 1

2

1

EDU>q = [1 2 3 1] q = 1

2

3

1

EDU>r = conv(p,q) r = 1

4

8

9

5

1

EDU>

Thus

( ) = p(s)q (s) = s5 + 4s4 + 8s3 + 9s2 + 5s + 1:

r s

The roots command can be used to nd the roots of a polynomial. For instance, suppose we type: a = [1 2] b = [1 4] c = [1 8] p = conv{a,b} p = conv(p,c) roots(p)

What we have done is create the polynomials s

+2

;

s

+4

;

and

s

+ 8;

and then multiply them together to create the polynomial p(s). We then use the Matlab command roots to nd roots nd the roots ( which we already know). If we typed these commands in the Matlab command window we would get: EDU>a = [1 2] a = 1

2

2.1.

3

IMPORTANT COMMANDS

EDU>b = [1,4] b = 1

4

EDU>c = [1,8] c = 1

8

EDU>p = conv(a,b) p = 1

6

8

EDU>p= conv(p,c) p = 1

14

56

64

EDU>roots(p) ans = -7.99999999999998 -4.00000000000002 -2.00000000000000 EDU>

Note that we put commas between the entries in `b.' They aren't required, but they sometimes are helpful. Note also that we have some round-o error when we compute the roots. The residue command enables us to quickly nd the residues in a partial fraction expansion. For instance, suppose we have Y

(s) =

10(s + 4) s(s + 10)(s + 20)

4

CHAPTER 2.

THE LAPLACE TRANSFORM

We rst create the two polynomials ( ) = 10(s + 4)

a s

() =

b s

( + 10)(s + 20);

s s

and then use the Matlab command [R,P,k} = residue[a,b]

The vector R contains the residues, the vector P the corresponding pole and the vector K the direct terms. The direct terms would occur if we had ( )=

F s

C1

s

+ p1

+

C2

s

+ p2

+ :::+

C s

n

+ pn

+ K (s):

We won't often encounter a situation where the term K (s) is present, but Matlab can handle it. For the Y (s) considered here, the Matlab entries and responses are as follows EDU>a = [1 4] a = 1

4

EDU>a = 10*a a = 10

40

EDU>c = [1 0] c = 1

0

EDU>d = [1 10] d = 1

10

EDU>e = [1 20] e =

2.1.

5

IMPORTANT COMMANDS

1

20

EDU>b = conv(c,d) b = 1

10

0

EDU>b = conv(b,e) b = 1

30

200

0

EDU>[R,P,K] = residue(a,b) R = -0.80000000000000 0.60000000000000 0.20000000000000 P = -20 -10 0 K = [ ] EDU>

Thus, the partial fraction expansion is 0:2 0:6 10(s + 4) = + s(s + 10)(s + 20) s s + 10

0:8 : s + 20

Chapter 3

The Transfer Function 3.1

Useful Commands

Several useful Matlab commands are 1. mesh 2. view 3. print 4. tf The command mesh generates a three dimensional plot similar to those in Figures 3.4 and 3.5. The following program was used to generate Figure 3.4. xt = 3 xb = -6 yt = 3 yb = -3 dx = (xt -xb)/100 dy = (yt -yb)/100 x=xb y=yb for l = 1:100 for k = 1:100 s = x + i*y; a = sqrt( (x + 2)^2 + y^2); b = sqrt(x^2 + y^2); Z(l,k) = 20*log( b/a );

1

2

CHAPTER 3.

THE TRANSFER FUNCTION

if Z(l,k) > 10000 Z(l,k) = 10000; end X(k) = x; x = x + dx; end x = xb; Y(l) = y; y = y + dy; end mesh(X,Y,Z) view(37,20) print tent.eps clear X clear Y clear Z

This program is typical of the kind of Matlab programs I write: crude, short, and eective. There are a lot of sophisticated Matlab programs that come in handy once in a while, but from my perspective the most eective way to use Matlab is to write short programs to solve speci c problems, problems that almost invariably fail to t nicely into some `canned' routine. Note in this problem that I set up the range of the variables x and y , namely xb, xt, yb, and yt, and then de ned the increments dx and dy by which x and y would change in the nested `for' loop. The nested `for' loops in this program are not very eÆcient because Matlab is an interpreted language, that is it is not compiled and then executed, although that is changing because the latest edition comes with a compiler. However, for routine use it is still an interpreted language. This proves to be inconsequntial because the current generation of desktop computers are so fast that even ineÆcient programs run in a second or two. Therefore, it is not worth spending an hour writing an eÆcient program when you can write a crude one in ve minutes that will execute in two seconds or less. This is my view of Matlab. With rare exceptions, the programs I write are programs that I will use to solve a speci c problem and then discard. But they are very, very useful, because they let me do analysis, like simple optimizations, that I would not even attempt without a tool like Matlab. Returning to the code, once the nested `for' loops have calculated all the values of the functions, and values of x, y, and z have been stored

3.1.

3

USEFUL COMMANDS

in the approprtiate vectors, then two commands mesh and view are used. The command mesh draws the three dimensional plot. The command view determines the angle at which you view the plot. The command has two arguments, so that when we use it we tupe view(AZ,EL),

where the rst argument is the azimuth from which the plot is viewed and second argument is the elevation from which the plot is viewed. The azimuth and elevation chosen, combined with the order in which the data was stored in the vectors X, Y, and Z, yields what I consider the `natural' view of the plot, with values on the x-axis increasing from right to left and values on the y -axis increasing from front to back. If you don't view this as `natural,' type: help view

to nd your own preferred perspective. The command print tent.eps stored the gure as a postscript le called `tent.eps.' You could also use print -deps tent.eps,

to store the gure as encapsulated postscript. Indeed if you type help print

you will nd that there are more ways to store the gure than I care to enumerate here. The command tf(a,b)

creates a transfer function from the the two polynomials a(s) and b(s), where a(s) is the numerator polynomial and b(s) is the denominator polynomial. For instance, suppose ( )=

G s

10(s + 1) : ( + 2)(s + 10)

s s

Then the following Matlab `dialogue' creates this transfer function.

4

CHAPTER 3.

EDU>a = 10*[1 1] a = 10

10

EDU>c = [1 0] c = 1

0

EDU>d = [1 2] d = 1

2

EDU>e = [1 10] e = 1

10

EDU>b = conv(c,d) b = 1

2

0

EDU>b = conv(b,e) b = 1

12

20

EDU>g = tf(a,b) Transfer function: 10 s + 10

0

THE TRANSFER FUNCTION

3.1.

USEFUL COMMANDS

------------------s^3 + 12 s^2 + 20 s

5

Chapter 4

Introducing Feedback 4.1

Useful Commands

Typing help control

will generate all the commands in the Matlab control tool box. There are a great many commands. Several that are worth introducing at this juncture are: 1. pole 2. damp 3. pzmap 4. step The command \verb+pole+ finds the poles of a transfer function. For instance, consider the following Matlab dialogue. EDU>a = [5 1] a = 5

1

EDU>b =[1 13 44 32]

1

2

CHAPTER 4.

INTRODUCING FEEDBACK

b = 1

13

44

32

EDU>g = tf(a,b) Transfer function: 5 s + 1 -----------------------s^3 + 13 s^2 + 44 s + 32 EDU>pole(g) ans = -7.999999999999991e+00 -4.000000000000006e+00 -9.999999999999998e-01 EDU>

The pole are at s

= 1;

4; and

8:

Note that there is some round o error. The command damp nds the damping ratio and natural frequency of the poles of a transfer function. To see how this command works consider the following Matlab dialogue. EDU>a = 18 a = 18 EDU>b = [1 7 24 18] b = 1

7

24

18

4.1.

3

USEFUL COMMANDS

EDU>g = tf(a,b) Transfer function: 18 ----------------------s^3 + 7 s^2 + 24 s + 18 EDU>damp(g) Eigenvalue

Damping

Freq. (rad/s)

-1.00e+00 -3.00e+00 + 3.00e+00i -3.00e+00 - 3.00e+00i

1.00e+00 7.07e-01 7.07e-01

1.00e+00 4.24e+00 4.24e+00

EDU>

Matlab uses the term eigenvalue for poles of the transfer function. We will encounter the terminology eigenvalue in Chapter 18. when we begin talking about state models. We will see there that the eigenvalues of a state model correspond to the poles of a transfer function. The command pzmap plots the poles of a transfer function. For g

=

3 s

+ 7s2

18 ; + 24s + 18

typing pzmap(g)

yields Figure 4.1 The command step generates the unit step response of a transfer function. The following Matlab dialogue generates the step response of the simple transfer function and saves the reponse as an encapsulated postscript le. 5 G(s) = : s+ 5 EDU>a = 5 a = 5

4

CHAPTER 4.

INTRODUCING FEEDBACK

Pole-zero map 3

2

Imag Axis

1

0

-1

-2

-3 -3.5

-3

-2.5

-2

-1.5 Real Axis

-1

-0.5

Figure 4.1: Pole/zero Map Created bu Command pzmap(g)

0

4.1.

USEFUL COMMANDS

5

EDU>b = [1 5] b = 1

5

EDU>g = tf(a,b) Transfer function: 5 ----s + 5 EDU>step(g) EDU>print -deps stepg.eps EDU>

The response is shown in Figure 4.2 There are many more commands in the Matlab controls toolbox worth investigating. We will discuss some of these in subsequent chapters, but the student is encouraged to investigate these commands.

6

CHAPTER 4.

INTRODUCING FEEDBACK

Step Response 1 0.9 0.8 0.7

Amplitude

0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.2

0.4

0.6

0.8

1

Time (sec.)

Figure 4.2: Step Response to Command step(g)

1.2

Chapter 1

Root Locus Analysis 4.1

Useful Commands

We introduce here the commands 1. rlocus, 2. rlocfind, The simplest way to draw a root locus is to use the rlocus command in program chap5rld,m, repeated below for convenience. %Root locus drawn using rlocus(g,gain) where gain is vector of gains that %user generates. gainmax = 150 gainmin = 0 npts = 400 dgain = (gainmax - gainmin)/npts newgain = gainmin % % generate gain vector % gain = linspace(0,400,npts) % % generate transfer function % a = [1 3] b = [1 19 18 0] g = tf(a,b) % % generate root locus % rlocus(g)

1

2

CHAPTER 4.

ROOT LOCUS ANALYSIS

1.5

1

Imag Axis

0.5

0

-0.5

-1

-1.5

-15

-10

-5 Real Axis

0

5

Figure 4.1: Root Locus Generated by Use of rlocus(g) Command print -deps chap5rld.eps

The result of this program is shown in Figure 4.1 Note the stray line from the zero up to one limb of the root locus. This program is not very good, and I seldom use it. We can get around the problem of stray lines by using the Matlab program chap5rla.m, repeated below for convenience that draws the root locus shown in Figure 4.2. %Root locus drawn using rlocus(g,gain) where gain is vector of gains that %user generates. gainmax = 150 gainmin = 0 npts = 400 dgain = (gainmax - gainmin)/npts newgain = gainmin % % generate gain vector % for i = 1:npts gain(i) = newgain; newgain = newgain + dgain; end % % generate transfer function

4.1.

3

USEFUL COMMANDS

8

6

4

Imag Axis

2

0

-2

-4

-6

-8

-18

-16

-14

-12

-10 -8 Real Axis

-6

-4

-2

0

Figure 4.2: Root Locus Generated by Use of rlocus(g,gain) Command % a = [1 3] b = [1 19 18 0] g = tf(a,b) % % generate root locus % rlocus(g,gain) print -deps chap5rla.eps

The command rlocus(g,gain)

prints an x at each point of the root locus. Figure 4.2 the x's are so close together they form a thick solid line. We can spread the x's out by plotting fewer points. We can also write our own root locus program. Suppose for unity feedback we have (s2 + 12s + 72) + 12s3 + 72s2 + 120s + 100 (s + 6 j 6)(s + 6 + j 6) ; = (s + 1 j )(s + 1 + j )(s + 5 j 5)(s + 5 + j 5)

( ) =

Gp s

4 s

4

CHAPTER 4.

ROOT LOCUS ANALYSIS

and

160(s + 1)2 s(s + 24:71) Then the Matlab program chap5rl.m listed below will plot the root locus for 0 < K ( )=

Gc s

%Root Locus Program chap5rl.m gaininc = 0.25 gainmin = 0 gainmax = 160 gain = gainmin; kount = 1; while (gain <= gainmax); a = [1,2,1]; b= [1,12,72]; top = conv(a,b); c = [1,0]; d = [1,24.71]; e =[1,10,50]; f = [1,2,2]; bot = conv(c,d); bot = conv(bot,e); bot = conv(bot,f); roots(top); roots(bot); top = top * gain; top=[0,0,top]; p= bot + top; x = roots(p); a1(kount) = real ( x(1) ); b1(kount) = imag ( x(1) ); c1(kount) = real( x(2) ); d1(kount) = imag ( x(2) ); e1(kount) = real( x(3) ); f1(kount) = imag( x(3) ); g1(kount) = real( x(4) ); h1(kount) = imag( x(4) ); i1(kount) = real( x(5) ); j1(kount) = imag( x(5) ); k1(kount) = real( x(6) ); l1(kount) = imag( x(6) ); gain = gain + gaininc; kount = kount + 1; end save real1.dat a1 -ascii

<

160

4.1.

5

USEFUL COMMANDS

8

6

4

2

0

-2

-4

-6

-8 -25

-20

-15

-10

-5

0

Figure 4.3: Root Locus Generated by Program chap5rl.m save imag1.dat b1 -ascii save real2.dat c1 -ascii save imag2.dat d1 -ascii save real3.dat e1 -ascii save imag3.dat f1 -ascii save real4.dat g1 -ascii save imag4.dat h1 -ascii save real5.dat i1 -ascii save imag5.dat j1 -ascii save real6.dat k1 -ascii save imag6.dat l1 -ascii plot(a1,b1,'k.',c1,d1,'k.',e1,f1,'k.',g1,h1,'k.',i1,j1,'k.',k1,l1,'k.') print -deps chap5rl.eps

In the plot command 'k.' places a '.' at each point that is on the root locus. The command print -deps chap5rl.eps

saves the root locus plot as an encapsulated postscript le. The root locus is shown in Figure 4.3 The command rlocfind can be used to nd the gain associated with a speci c set of poles. Suppose we have displayed in the Matlab command window the root locus shown in Figure 4.2. Then the following Matlab diaglogue allows us to nd the gain associated with a particular set of poles.

6

CHAPTER 4.

ROOT LOCUS ANALYSIS

EDU>[K,poles]=rlocfind(g) Select a point in the graphics window selected_point = -6.146464646464647e+00 + 3.820512820512819e+00i K = 1.167208589859281e+02 poles = -6.160840653110633e+00 + 3.804837293076535e+00i -6.160840653110633e+00 - 3.804837293076535e+00i -6.678318693778726e+00

When we enter the command [K,poles]=rlocfind(g),

a set of crosshairs appear, the intersection of which can be controlled by the computer mouse. When we have the maneuvered the crosshairs to the desired position, clicking the mouse causes the value of K and the poles locations shown above to be displayed. In Figure 4.4 the location of the three poles corresponding to the gain 162.72 are marked by the + signs. In general, I nd the Matlab plotting routines to be diÆcult to use and inferior to the quality of graph that can be obtained with a software package like Deltagraph. Even Excel is superior to Matlab for graphing data. So the user has to make the choice whether to use Matlab to generate the plot, or to use Matlab to generate the data, and then plot the data using some other software package.

7

USEFUL COMMANDS

8

6

4

2 Imag Axis

4.1.

0

-2

-4

-6

-8

-18

-16

-14

-12

-10 -8 Real Axis

-6

-4

-2

0

Figure 4.4: Root Locus Generated by Use of rloc nd Command

Chapter 6

Specifying Performance 6.1

Useful Commands

Some useful commands that we can use with this chapter are 1. damp 2. pole 3. zpk 4. zpkdata 5. tfdata 6. ltiview The Matlab command damp nds the natural frequency all the poles of transfer function. For instance, let Tc(s)

=

(s + 4

!n

and the damping ratio

250

j 3)(s + 4 + j 3)(s + 10)

The following Matlab dialogue shows how to use zpk and damp. EDU>g = zpk([],[-4+j*3,-4-j*3,-10],250) Zero/pole/gain: 250 ---------------------(s+10) (s^2 + 8s + 25) EDU>damp(g) Eigenvalue

Damping

Freq. (rad/s)

1



for

2

CHAPTER 6.

-4.00e+00 + 3.00e+00i -4.00e+00 - 3.00e+00i -1.00e+01

8.00e-01 8.00e-01 1.00e+00

SPECIFYING PERFORMANCE

5.00e+00 5.00e+00 1.00e+01

EDU>

The arguments of the command zpk are: 1. [ ] indicating there are no zeros in the transfer function 2. [-4j*3,-4-j*3,-10]+ de ning the poles of the system 3. 250 de ning the gain of the system. The command damp(g)

6.1.

USEFUL COMMANDS

3

The extracts the damping ratio and naturalfrequency of all poles of the system. The matlab command tfdata returns the zeros, poles and gain of the system. For instance, for the system we just de ned we have EDU>[Z,P,K] = zpkdata(g,'v') Z = Empty matrix: 0-by-1 P = -4.000000000000000e+00 + 3.000000000000000e+00i -4.000000000000000e+00 - 3.000000000000000e+00i -1.000000000000000e+01 K = 250 EDU>

The Matlab command tfdata works in a similar way, as demonstrated by the following Matlab dialogue. EDU>num = 250 num = 250 EDU>p1 = [1 4-j*3] p1 = Column 1 1.000000000000000e+00 Column 2 4.000000000000000e+00 - 3.000000000000000e+00i

4

CHAPTER 6.

SPECIFYING PERFORMANCE

EDU>p2 = [1 4+j*3] p2 = Column 1 1.000000000000000e+00 Column 2 4.000000000000000e+00 + 3.000000000000000e+00i EDU>p3 = [1 10] p3 = 1

10

EDU>den = conv(p1,p2) den = 1

8

25

EDU>den = conv(den,p3) den = 1

18

105

250

EDU>g1 = tf(num,den) Transfer function: 250 -------------------------s^3 + 18 s^2 + 105 s + 250 EDU>[n,d] = tfdata(g1,'v') n = 0

0

0

250

6.1.

USEFUL COMMANDS

5

d = 1

18

105

250

EDU>

One additional command of some use is ltiview. This command opens a graphical user interface (GUI) window that allows the user to select any system currently de ned in the workspace and make various plots, such as the step response, Bode magnitude and phase plot, nyquist plots, and so on. At this juncture in the book the only plots that the reader `oÆcially' knows about are the step and impulse responses. One of the features of this GUI window is that the reader can select various characterstics of the plot and have them appear on the plot. Of course, this could be done almost as easily using the command plot, and each user has to decide which display method is best. It is really a matter of taste. This command is probably best used in the early stages of a design process to compare two or three preliminary designs.

Chapter 7

Root Locus Design 7.1

Useful Commands

Most of the commands that pertain to root locus design were discussed in Chapters 5 and 6. One command of great use is the command feedback.

The following Matlab dialogue shows how it is used.

1

2

CHAPTER 7.

ROOT LOCUS DESIGN

EDU>clear all EDU>g = zpk([-3],[0 -1 -18], 78) Zero/pole/gain: 78 (s+3) -------------s (s+1) (s+18) EDU>h = 1 h = 1 EDU>tc = feedback(g,h) Zero/pole/gain: 78 (s+3) ---------------------(s+13) (s^2 + 6s + 18) EDU>step(tc,3) EDU>print -deps stepresp.eps

The feedback command lets us determine the closed loop transfer function for output feedback. Here we have chosen unity feedback with the command h = 1

We can then nd the unit step response, shown in Figure 7.1, with the command step(tc,3)

One additional command is of some use, namely ltiview. This command opens a graphical user interface (GUI) window that allows the user to select any system currently de ned in the workspace and make various plots, such as the step response, Bode magnitude and phase plots, Nyquist plots, and so on. At this juncture in the book the only plots that the reader `oÆcially' knows about are the step and impulse responses. One of the features of this GUI window is that the reader can select various characteristics of the plot and have them appear on the plot. Of course, this could be done almost as easily using the command plot, and each user has to decide which display method is best. It is really a matter of taste. This command is probably best used in the early stages of a design process to compare two or three preliminary designs.

3

USEFUL COMMANDS

Step Response 1.4

1.2

1

Amplitude

7.1.

0.8

0.6

0.4

0.2

0

0

0.5

1

1.5

2

2.5

Time (sec.)

Figure 7.1: Step Response of Closed Loop System

3

Chapter 8

Speed Control 8.1

Data Description

All the data les in this chapter are in column form with the rst column being the time and the second column being the system step respoonse. For the step responses of the d c motors the armature voltage response is also given. The motor data is contained in a folder called Motordata. The data for transfer function identi cation in problems 8,9, and 10 is contained in a folder called xferid. For problems 8,9, and 10, the inputis a unit step.

1

Chapter 9

Frequency Response 9.1

Useful Commands

Several useful Matlab commands are 1. bode, 2. axis, 3. subplot The Matlab program chap9bode1, repeated below for convenience, creates the Bode magnitude and phase plots shown in Figure 9.1. a = 100*[1 5] b = [1 19 18 0] g = tf(a,b) grid on bode(g) print -deps chap9bode1.eps

The transfer function in this case is ( )=

Gp s

100(s + 5) : s(s + 1)(s + 18)

The disadvantage of this plot is that I would prefer to have the scaling of the magnitude plot be in increments of 20 dB because in the transfer indenti cation process we look for slopes of  20, 40 and 60 db/dec. We can write our own program to compute the Bode magnitude and phase plots. An example is program bodeplot.m, listed below for convenience.

1

2

CHAPTER 9.

FREQUENCY RESPONSE

Bode Diagrams

100

Phase (deg); Magnitude (dB)

50

0

-50 -80 -100 -120 -140 -160 -180 -2 10

-1

10

0

10

1

10

2

10

Frequency (rad/sec)

Figure 9.1: Bode Plots Generated by Program chap9bode1.m fmin =0.01 fmax = 1000 ftop = log10(fmax) fbot = log10(fmin) npts = 100 dw = ( ftop - fbot) / npts for i = 1:npts w = fmin*( 10^(i*dw) ); s = j* w; mag(i) = 20*log10(abs(100*(s + 3)/( s * (s + 10)*(s + 100)))); phase(i) =( atan(w/3)-pi/2 -atan(w/10) -atan(w/100) )*180/pi; radfreq(i) = w; end subplot(2,1,1),semilogx(radfreq,mag) grid on axis([0.1 1000 -100 20]) subplot(2,1,2), semilogx(radfreq,phase) grid on axis([0.1 1000 -180 0]) print -deps chap9bode2.eps

The result is shown in Figure 9.2 Note we have used the subplot command here to get both

9.1.

3

USEFUL COMMANDS

20 0 -20 -40 -60 -80 -100 -1 10

0

10

1

10

2

10

3

10

0

-50

-100

-150 -1

10

0

10

1

10

2

10

3

10

Figure 9.2: Bode Plots Generated by Program bodeplot.m the magnitude and phase plots. This is a reasonable looking plot, but if you want easy control of the axes, tick marks and line widths your best bet is to generate the data using Matlab and then transfer the data to a good graphics package.

Chapter 1

The Nyquist Criterion 9.1

Useful Commands?

We discuss brie y the Matlab commands 1. nyquist 2. nichols and conclude they are worthless. The program chap10nyq.m, repeated below for convenience draws the Nyquist plot shown in Figure 9.1. As can be seen it is of very little help because of the scaling. %Program to draw Nyquist Plot for % G(s) = 100/[s (s + 10)(s + 40)] %user generates. % % generate transfer function % a = 100 b = [1 50 400 0] roots(b) g = tf(a,b) % % generate Nyquist plot % nyquist(g) print -deps chap10nyq.eps

The command nichols is no better. The program chap10nic.m, repeated below for convenience, generates the plot shown in Figure 9.2. %Program to draw Nichols Plot for % G(s) = 100/[s (s + 10)(s + 40)] %user generates. %

1

2

CHAPTER 9.

THE NYQUIST CRITERION

Chapter 11

Bode Design 11.1

Useful Commands

We have already discussed the basic command for frequency response analysis, namely bode. Three other closely related commands are 1. evalfr, 2. freqresp, 3. margin The rst two are very minor variations of bode. The third is not much more than that, but we illustrate it with the short Matlab program chap11margin.m, repeated below for convenience. This program creates the Bode magnitude and phase plots shown in Figure 11.1. a = 100000*[1 5] c = [1 0] d = [1 1] e = [1 50] f = [1 100] b = conv(c,d) b = conv(b,e) b = conv(b,f) g = tf(a,b) grid on margin(g) print -deps chap11mar.eps

As usual Matlab has chosen very inconvenient scaling for the Bode magnitude plot. This command to me is only of marginal worth, as are evalfr and freqresp.

1

2

CHAPTER 11.

BODE DESIGN

Bode Diagrams Gm=16.5 dB (Wcg=66.3); Pm=46.7 deg. (Wcp=19.0) 100

50

Phase (deg); Magnitude (dB)

0

-50

-100 -50 -100 -150 -200 -250 -300 -2 10

-1

10

0

10

1

10

2

10

3

10

Frequency (rad/sec)

Figure 11.1: Bode Plots Generated by Program chap11margin,m

Chapter 12

Robust Control There is a robust control tool box, but it is devoted to multi-input, multi-output models. A useful exercise is to write a program to nd the line of closest approach to the point 1 in the GH -plane.

1

Chapter 13

Position Control 13.1

Impulse Data

The data is stored in column form. The rst column is the time in seconds the second column the impulse response of the postioning table. The armature voltage is also provided for each turntable. The lename ptkixx

is decoded as: postioning table `k' impulse of xx miliseconds. The armature voltage lename is similar but with the suÆx `av.' 13.2

Frequency Response Data

The data is stored in column form as shown in Figure 13.1 Freq.(Rad/sec) Arm Voltage Heliopot Voltage vvv www xxx .. .. .. . . .

Mag.(dB) Phase(deg) yyy zzz .. .. . .

Figure 13.1: Data Format for Frequency Response Data

1

Chapter 14

Discrete Systems 14.1

Useful Commands

We introduce here the commands 1. c2d, 2. step, 3. zpk 4. bode We have talked about the continuous versions of the last three instructions already. The following Matlab dialogue forms the tranfer function

G(s) =

a ; s+a

and then converts it to the corresponding zero order hold digital equivalent for f = 10 hz. EDU>z = [ ] z = [ ] EDU>p = -5

1

2

CHAPTER 14.

DISCRETE SYSTEMS

p = -5 EDU>K = 5 K = 5 EDU>g = zpk(z,p,K) Zero/pole/gain: 5 ----(s+5) EDU>gz = c2d(g,0.1) Zero/pole/gain: 0.39347 ---------(z-0.6065) Sampling time: 0.1 EDU>

In getting gz we have used the command zpk to form the continuous transfer function before converting it to the discrete version. The command c2d can form several types of discrete equivalents, two of which are discussed in Chapter 15. We have used the default version of zpk which does the zero order hold equivalent. We could have typed gz = c2d(g,0.1,\zoh')

and obtained the same result. The use of the the Matlab command step for discrete systems is demonstrated by the following dialogue, where we use the function gz obtained above. EDU>step(gz,0:0.1:4)

14.1.

3

USEFUL COMMANDS

Step Response 1.4

1.2

Amplitude

1

0.8

0.6

0.4

0.2

0

0

0.5

1

1.5

2

2.5

3

3.5

4

Time (sec.)

Figure 14.1: Step Response Generated by step(gz,0:0.1:4) EDU>print -deps dstep.eps EDU>

The notation `0:0.1:4' sets the starting time (t = 0), the ending time (t = 4) and the intersample period, T = 0.1. The second command stores the time response as an encapsulated postscript le. That le is shown in Figure 14.1 The Matlab command zpk can be used with discrete systems as demonstrated by the following dialogue. EDU>gz1 = zpk([ ],-0.5,2,0.1) Zero/pole/gain: 2 ------(z+0.5) Sampling time: 0.1

4

CHAPTER 14.

DISCRETE SYSTEMS

Step Response 2.5

2

Amplitude

1.5

1

0.5

0

0

0.5

1

1.5

2

2.5

3

3.5

4

Time (sec.)

Figure 14.2: Step Response of gz1 EDU>step(gz1,4) EDU>print -deps dstepgz1.eps EDU>

Note that we did not have to specify the sampling rate since when we created gz1 we speci ed the sampling rate of 10 hz. The step response is shown in Figure 14.2. The command bode can be used for discrete systems as well. In the following Matlab dialogue we nd the frequency response of

G(z ) =

z

0:39347 : 0:6065

EDU>bode(gz) EDU>print -deps dbode.eps EDU>

The Bode phase and magnitude plots are shown in Figure 14.3 In this par-

5

USEFUL COMMANDS

Bode Diagrams

0

-5

Phase (deg); Magnitude (dB)

14.1.

-10

-15 0

-50

-100

-150

-200 0 10

1

10

Frequency (rad/sec)

Figure 14.3: Bode Phase and Magnitude Plots

2

10

6

CHAPTER 14.

DISCRETE SYSTEMS

ticular case the sampling rate is 10 Hz so that the maximum frequency that can be sampled without aliasing is

! = =T = 31:4 rad/sec: The vertical dark lines denote the maximum sampling frequency.

Chapter 15

Digital Control 15.1

Useful Commands

In Chapter 14 we discussed most of the commands of use in analyzing digital control systems. The command c2d discussed in chapter 14 has several alternative modes. Two of importance to us are the bilinear and ad hoc mappings, both of which can be done with c2d. In addition the commands 1. d2d 2. damp, and 3. filt 4. ltiview are also of some use. The command c2d rst discussed in Chapter 14, can be used to determine the bilinear and ad hoc digital equivalents of continuous transfer functions. Consider the transfer function G(s) =

10(s + 2) ; (s + 1)(s + 4)

discussed in Sections 15.2.3 and 15.2.4 of the text. Then the following Matlab dialogue creates the equivalent discrete transfer function under the bilinear (Tustin's) mapping. EDU>g = zpk(-2,[-1 -4],10) Zero/pole/gain: 10 (s+2) ----------(s+1) (s+4)

1

2

CHAPTER 15.

DIGITAL CONTROL

EDU>gz = c2d(g,0.1,'tustin') Zero/pole/gain: 0.43651 (z-0.8182) (z+1) -----------------------(z-0.9048) (z-0.6667) Sampling time: 0.1

The following dialogue generates the ad hoc version of

G(s).

EDU>gzah = c2d(g,0.1,'matched') Zero/pole/gain: 0.86538 (z-0.8187) --------------------(z-0.9048) (z-0.6703) Sampling time: 0.1

The command d2d can be used to change an existing discrete transfer function at one sampling frequency to an equivalent discrete transfer function at a new sampling frequency. For instance, we found the discrete equivalent of G(s)

=

10(s + 2) (s + 1)(s + 4)

using the bilinear mapping and a sampling frequency of 10 Hz. The following Matlab dialogue shows how to change to a sampling rate of 20 Hz. EDU>gz20 = d2d(gz,0.05) Zero/pole/gain: 0.43651 (z+0.0758) (z-0.9046) ----------------------------(z-0.9512) (z-0.8165) Sampling time: 0.05

The Matlab command damp can be used to nd the equivalent natural frequency of damping ratio of complex poles in the z -plane. The following Matlab dialogue rst forms the continuous transfer function G(s)

=

(s + 4

25 ; j 3)(s + 4 + j 3)

nds the discrete, zero order hold equivalent, and then uses damp to nd !n and  .

15.1.

3

USEFUL COMMANDS

EDU>gs = zpk([],[-4+j*3,-4-j*3],25) Zero/pole/gain: 25 --------------(s^2 + 8s + 25) EDU>gzc = c2d(gs,0.1,'zoh') Zero/pole/gain: 0.095495 (z+0.7652) ----------------------(z^2 - 1.281z + 0.4493) Sampling time: 0.1 EDU>damp(gzc) Eigenvalue

Magnitude

Equiv. Damping

Equiv. Freq. (rad/s)

6.40e-01 + 1.98e-01i 6.40e-01 - 1.98e-01i

6.70e-01 6.70e-01

8.00e-01 8.00e-01

5.00e+00 5.00e+00

The matlab command filt creates a digital lter using the digital signal processing preference for powers of z 1 . The followng Matlab dialogue creates such a lter. EDU>a = [1 2 2 1] a = 1

2

2

1

EDU>b = [1 2 4 2 ] b = 1

2

4

2

EDU>gdsp = filt(a,b) Transfer function: 1 + 2 z^-1 + 2 z^-2 + z^-3 ---------------------------1 + 2 z^-1 + 4 z^-2 + 2 z^-3

4

CHAPTER 15.

DIGITAL CONTROL

Sampling time: unspecified EDU>gdsp10 = filt(a,b,0.1) Transfer function: 1 + 2 z^-1 + 2 z^-2 + z^-3 ---------------------------1 + 2 z^-1 + 4 z^-2 + 2 z^-3 Sampling time: 0.1

The Matlab command ltiview was discussed in Chapter 7. It is equally applicable to discrete systems.

Chapter 16

Aircraft Pitch Control 16.1

Impulse Data

The data is stored in column form. The rst column is the time in seconds the second column the impulse response of the postioning table. The armature voltage is also provided for each turntable. The lename ptkixx

is decoded as: postioning table `k' impulse of xx miliseconds. The armature voltage lename is similar but with the suÆx `av.' 16.2

Frequency Response Data

The data is stored in column form as shown in Figure 16.1 Freq.(Hz) Freq.(Rad/sec) Arm Voltage vvv www xxx .. .. .. . . .

Heliopot Voltage yyy .. .

Phase(deg) zzz .. .

Figure 16.1: Data Format for Frequency Response Data

1

Chapter 17

The Transportation Lag 17.1

Useful Commands

There is not much in the Matlab tool boxes that will help us with a system that has a transportation lag. You can either use Simulink or one of the two Matlab programs ufeedbak or ufeedbakd. The program ufeedbak is an interactive program that requires you to enter the data from the keyboard. It simulates the response of a continuous system to inputs of the form ( ) = at3 + bt2 + ct + d:

u t

The program ufeedbakd does the same thing for sampled data systems. If you do not like the interactive nature of these programs, then I encourage you to modify them to your own satisfaction. The program rwtest.m that uses the data le called data is short introduction to how to read and write data in Matlab.

1

Chapter 18

The State Model 18.1

Useful Commands

Several useful Matlab commands related to the analysis of state models. Five are listed below. 1. ss, 2. ss2ss, 3. canon 4. eig 5. inv The Matlab program chap18ss, repeated below for convenience, creates the state space model called stmod A = [0 1 0;0 0 1;0 -2 -3] B = [0;0;1] C = [1 0 0] D = 0 stmod = ss(A,B,C,D)

The output generated in the Matlab command window by executing this program is A = 0 0 0

1 0 -2

0 1 -3

1

2

CHAPTER 18.

THE STATE MODEL

B = 0 0 1 C = 1

0

0

D = 0 a = x1 x2 x3

x1 0 0 0

x1 x2 x3

u1 0 0 1.00000

y1

x1 1.00000

y1

u1 0

b =

c =

d =

x2 1.00000 0 -2.00000

x3 0 1.00000 -3.00000

x2 0

x3 0

Continuous-time system.

The system stmod can then be used as an argument for many of the Matlab commands we have already discussed. For instance, the command

18.1.

3

USEFUL COMMANDS

Step Response 2.5

2

Amplitude

1.5

1

0.5

0

0

1

2

3

4

5

Time (sec.)

Figure 18.1: Step Response of State Model stmod step(stmod)

creates the step response shown in Figure 18.1 As another example, typing eig(stmod)

in the command window produces the dialogue EDU>eig(stmod) ans = 0 -1 -2 EDU>

Thus the eigenvalues of the plant matrix A are  = 0 :  = 1 ; and = 2: If we type [T,F] = eig(A)

6

4

CHAPTER 18.

THE STATE MODEL

we obtain T = 1.000000000000000e+00 0 0

-5.773502691896258e-01 5.773502691896258e-01 -5.773502691896258e-01

2.182178902359923e-01 -4.364357804719847e-01 8.728715609439694e-01

1.500000000000000e+00 3.464101615137753e+00 2.291287847477920e+00

4.999999999999999e-01 1.732050807568877e+00 2.291287847477920e+00

F = 0 0 0

0 -1 0

0 0 -2

If we then type tinv = inv(T)

we obtain tinv = 1.000000000000000e+00 0 0

If we next type cstmod = ss2ss(stmod,tinv)

we obtain a =

x1 x2 x3 x1 1.82407e-33 -6.40988e-17 1.93816e-16 x2 2.16731e-16 -1.00000 3.05716e-16 x3 0 0 -2.00000

b = x1 x2 x3 c =

u1 0.50000 1.73205 2.29129

18.1.

5

USEFUL COMMANDS

y1

x1 1.00000

y1

u1 0

d =

x2 -0.57735

x3 0.21822

Continuous-time system.

Except for the round o error, we now have the original system stmod in Jordan canonical form. We can accomplish the same thing with the canon command. That is, if we type cstmod = canon(stmod)

we obtain a = x1 x2 x3

x1 0 0 0

x1 x2 x3

u1 0.50000 1.73205 2.29129

y1

x1 1.00000

y1

u1 0

b =

c =

d =

x2 0 -1.00000 0

x3 0 0 -2.00000

x2 -0.57735

x3 0.21822

Continuous-time system.

This is the same model, sans the round o error.

Chapter 19

Observability and Controllability 19.1

Useful Commands

The following commands are of use in studying controllability and observability, and in converting continuous transfer functions to both continuous and discrete state models. 1. ctrb 2. obsv The Matlab program dialogue below forms the matrices

x_ ( ) = Ax( ) + b t

t

A and b for the system

()

u t ;

and then uses the Matlab commands ctrb and rank to from the controllability matrix and nd its rank. EDU>A = [0 1; 0 0] A = 0 0

1 0

EDU>b = [0;1] b = 0 1 EDU>co = ctrb(A,b)

1

2

CHAPTER 19.

OBSERVABILITY AND CONTROLLABILITY

co = 0 1

1 0

EDU>rank(co) ans = 2 EDU>

A

Thus, we see that even though does not have full rank the system is controllable. The Matlab dialogue below forms the matrices and for the system

x_ ( )

= y (t) = t

A Ax( ) cT x

c

t

and then forms the observability matrix and checks its rank. A = [0 1;-2 -3] A = 0 -2

1 -3

EDU>c = [0 1] c = 0

1

EDU>ob = obsv(A,c) ob = 0 -2

1 -3

EDU>rank(ob) ans =

19.1.

3

USEFUL COMMANDS

2 EDU>

The command c2d can be useful here in nding the discrete state model equivalent of a transfer function, or continuous state model. For instance, the following dialogue shows how we can turn a transfer function rst into a continuous state model and then into the step invariant discrete state model. EDU>g = zpk([],[0 -2],2) Zero/pole/gain: 2 ------s (s+2) EDU>stmod = ss(g) a = x1 x2

x1 -2.00000 1.00000

x1 x2

u1 1.41421 0

y1

x1 0

y1

u1 0

b =

c =

d =

x2 0 0

x2 1.41421

Continuous-time system. EDU>stmodd = c2d(stmod,0.1,'zoh') a =

4

CHAPTER 19.

x1 x2

x1 0.81873 0.09063

x1 x2

u1 0.12818 0.00662

y1

x1 0

y1

u1 0

b =

c =

d =

Sampling time: 0.1 Discrete-time system. EDU>

OBSERVABILITY AND CONTROLLABILITY

x2 0 1.00000

x2 1.41421

Chapter 20

The Controller/Observer I do not nd much in the Matlab control toolbox that will help you with the design of controller/observers. I have included some programs that I have written to support this activity. They are a start. They are by no means everything you would want. and that is intentional. These programs are crude and cumbersome, but with some work you can turn them into useful tools. The programs that begin with `build' compute the controller/observer matrices for some useful choices of dominant poles. The programs that begin with `sim' are simulation routines. These routines are far from perfect, but they work, and by studying them you should be able to improve them and make them your own. Good Luck.

1