Lecture notes for EE250 (Control Systems Analysis) Ramprasad Potluri Department of Electrical Engineering, IIT Kanpur January 2, 2012
EE250 (Control Systems Analysis) IITK
January 2, 2012
ii of 218
Lecture Notes
Ramprasad Potluri
Contents Preface
ix
Notation
xi
I
Preliminaries
1
1
Introduction to the course 1.1 About the title of the course . . . . . . . 1.2 Definition of negative feedback control 1.3 Examples of control systems . . . . . . . 1.4 Videos of control systems . . . . . . . . 1.5 Appendix: History of control systems .
. . . . .
3 3 3 4 5 6
2
Dynamic systems, transfer functions, and block diagrams 2.1 Dynamic systems . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Transfer functions, block diagrams . . . . . . . . . . . . . . . . . 2.3 Appendix: Detailed model of a DC motor . . . . . . . . . . . . .
7 7 8 10
3
Transfer functions and block diagrams revisted, state-space tions, an example of design in state space 3.1 Transfer functions revisited . . . . . . . . . . . . . . . . . . 3.2 Block diagrams revisited . . . . . . . . . . . . . . . . . . . . 3.3 State-space equations . . . . . . . . . . . . . . . . . . . . . . 3.4 An example of design in the state-space domain . . . . . .
13 13 14 15 17
4
5
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
equa. . . .
. . . .
. . . .
An example of design in the classical domain, block diagram manipulations 4.1 An example of design in the classical control domain . . . . . . 4.2 Rules of block diagram algebra . . . . . . . . . . . . . . . . . . . 4.3 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Appendix: A little more on state-space theory . . . . . . . . . . . 4.5 Appendix: Similarities and differences between classical and LTI state-space control . . . . . . . . . . . . . . . . . . . . . . . . . . . Signal flow graphs 5.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 SFG Algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Appendix: Usage of these terms . . . . . . . . . . . . . . . . . . . iii
21 21 22 23 25 25 27 27 29 29
EE250 (Control Systems Analysis) IITK
Lecture Notes
6
Mason’s gain formula, properties desired from a control system 6.1 Mason’s gain formula . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 A note on loop gain . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Properties desired from a control system . . . . . . . . . . . . . .
31 31 32 33
7
Virtues of control systems: Disturbance rejection
35
8
Virtues of control systems: Improvement in performance and overcoming deadzone 37 8.1 Improvement in speed of response . . . . . . . . . . . . . . . . . 37 8.2 Command feedforward . . . . . . . . . . . . . . . . . . . . . . . 39
II 9
Bode Plots
41
Frequency response 9.1 Definition of Frequency Response . 9.2 Problem statement . . . . . . . . . 9.3 Solution . . . . . . . . . . . . . . . 9.4 Conclusion . . . . . . . . . . . . . . 9.5 Caution . . . . . . . . . . . . . . . . 9.6 Appendix . . . . . . . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
43 43 43 43 45 45 46
10 Bode plots (Part 1) 10.1 Uses of frequency response in control systems . . . . . . . . . . 10.2 Convention in constructing Bode plots . . . . . . . . . . . . . . . 10.3 Types of Bode plots . . . . . . . . . . . . . . . . . . . . . . . . . . 10.4 Appendix: Angles in MATLAB . . . . . . . . . . . . . . . . . . . 10.5 Appendix: Angles in GNU Octave . . . . . . . . . . . . . . . . . 10.6 Appendix: Examples of departure of MATLAB from the convention . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47 47 48 50 50 51
11 Bode plots (Part 2) 11.1 First method for sketching an ABMP . . . . . . . . . . . . . . . . 11.2 Second method for sketching an ABMP . . . . . . . . . . . . . .
55 55 59
12 Bode plots (Part 3) 12.1 Comparison of the two methods for constructing ABMPs . . . 12.2 Selecting an ω to vertically fix an ABMP in the second method 12.3 Examples of construction of ABPs . . . . . . . . . . . . . . . . 12.4 BP and ABP of the second order TF . . . . . . . . . . . . . . . . 12.5 Appendix: MATLAB codes . . . . . . . . . . . . . . . . . . . .
. . . . .
61 61 61 62 62 63
13 Bode plots (Part 4) 13.1 Notes on the ABMP of second order TF 13.2 Complete BP of second order TF . . . . 13.3 Filter-related terminology used in BPs . 13.4 Minimum phase TF . . . . . . . . . . . . 13.5 Appendix: MATLAB codes . . . . . . .
. . . . .
65 65 65 67 67 67
January 2, 2012
iv of 218
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
51
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
III
Lecture Notes
Stability Theory
69
14 Stability theory (Part 1) 14.1 Minimum phase TFs (concluded) . . . . . . . . . . . . . 14.2 Why practical rational TFs are proper . . . . . . . . . . 14.3 BIBO stability of LTI systems** . . . . . . . . . . . . . . 14.4 Discussion of the theorem and its proof** . . . . . . . . 14.5 Appendix: Some literature that discusses BIBO stability
. . . . .
71 71 71 71 73 74
15 Stability theory (Part 2) 15.1 A question of interest . . . . . . . . . . . . . . . . . . . . . . . . . 15.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15.3 An interesting example . . . . . . . . . . . . . . . . . . . . . . . .
75 75 77 78
16 Nyquist Stability Criterion (Part 1)* 16.1 A quick recap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . 16.3 Cauchy’s theorem (Principle of argument) . . . . . . . . . . . . .
79 79 79 80
17 Nyquist Stability Criterion (Part 2)* 17.1 Application of Cauchy’s theorem in Nyquist’s stability criterion 17.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83 83 84
18 Nyquist Stability Criterion (Part 3) 18.1 Example 1: G (s) = K/s . . . . . . . . . . . . . . . 18.2 Convention followed in sketching Nyquist plots 18.3 Example: G (s) = K (s + 1) . . . . . . . . . . . . 18.4 Example: G (s) = K (s + 10)2 s3 . . . . . . . . . 18.5 Appendix: A useful fact . . . . . . . . . . . . . .
. . . . .
87 87 88 90 90 90
19 Nyquist Stability Criterion (Part 4) 19.1 Gain margin, phase margin . . . . . . . . . . . . . . . . . . . . . 19.2 When GM and PM are useful . . . . . . . . . . . . . . . . . . . . 19.3 Appendix: Dead-time element in MATLAB . . . . . . . . . . . .
91 91 92 93
20 Nyquist Stability Criterion (Part 5) 20.1 Demonstration of quick construction of Nyquist plots 20.2 Critical gain . . . . . . . . . . . . . . . . . . . . . . . . 20.3 Conditionally stable systems . . . . . . . . . . . . . . 20.4 Using − K1 + j0 instead of −1 + j0 . . . . . . . . . . . . 20.5 Nyquist plots in MATLAB and GNU Octave . . . . . 20.6 Rotation of a vector instead of encirclement . . . . . . 20.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . .
95 95 95 96 96 96 97 99
. . . . .
. . . . .
. . . . .
. . . . .
. . . . . . .
. . . . .
. . . . .
. . . . . . .
. . . . .
. . . . .
. . . . . . .
. . . . .
. . . . .
. . . . . . .
. . . . .
. . . . .
. . . . . . .
. . . . . . .
IV Design using bode plots-based loop-shaping for minimum-phase plants 101 21 NST (concluded), Loop-shaping (Part 1) 103 21.1 Nyquist stability theory (conclusion) . . . . . . . . . . . . . . . . 103 21.2 Specifications for control system design . . . . . . . . . . . . . . 103 January 2, 2012
v of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK 21.3 21.4 21.5 21.6
Lecture Notes
Loop gain and shaping . . . . . . . . . . . . . . . . . BMP-BPP correspondence for minimum-phase TFs Appendix: Bode’s theorem . . . . . . . . . . . . . . . Appendix: NST for SFG defintion of loop gain . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. 104 . 105 . 105 . 105
22 Loop-shaping (Part 2) 22.1 Recap of specifications for control system design 22.2 Gist of loop-shaping-based design . . . . . . . . 22.3 Stability . . . . . . . . . . . . . . . . . . . . . . . . 22.4 Robustness to plant parameter variations . . . . 22.5 Relation between OL and CL BMP . . . . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
107 107 107 108 108 108
23 Loop-shaping (Part 3) 23.1 Satisfaction of desired performance . . . 23.1.1 Steady-state error . . . . . . . . . 23.1.2 Settling time . . . . . . . . . . . . 23.1.3 Overshoot . . . . . . . . . . . . . 23.2 Respecting constraints on control input 23.3 Appendix: Additional loop-shaping tips
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
111 111 112 112 113 113 113
24 Loop-shaping (Part 4) 24.1 Robustness to disturbances . . . . . . . . 24.2 Robustness to noises . . . . . . . . . . . . 24.3 Summary of loop-shaping recipe . . . . . 24.4 MATLAB-based example of loop-shaping 24.4.1 Problem . . . . . . . . . . . . . . . 24.4.2 Solution . . . . . . . . . . . . . . . 24.4.3 How to use MATLAB . . . . . . . 24.4.4 Results . . . . . . . . . . . . . . . . 24.4.5 MATLAB code . . . . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
115 115 115 116 116 116 117 117 118 118
25 Loop-shaping (Part 5): Information from lead/lag controllers 25.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25.2 Lead and lag controllers . . . . . . . . . . . . . . . . . . . . . . . 25.3 How DD-Φmax relationship for lead/lag controllers is useful in loop-shaping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25.4 Appendix: Derivation of DD versus PM relationship . . . . . . .
123 123 123
26 Loop-shaping (Part 6): Relation between M p and PM 26.1 Relation between M p and PM for second order system . . . . . 26.2 Applying the second order relation to higher order . . . . . . . .
131 131 133
125 126
27 Loop-shaping (Part 7): A solved example 139 27.1 Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 27.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 28 Metrics for steady-state accuracy of unity-feedback around minimum-phase systems 28.1 Error constants . . . . . . . . . . . . . . . . . . . . 28.2 How to read K p , Kv , Ka off the BMPs . . . . . . . . 28.3 System type . . . . . . . . . . . . . . . . . . . . . . January 2, 2012
vi of 218
systems built 145 . . . . . . . . 146 . . . . . . . . 147 . . . . . . . . 147 Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
28.4 Appendix: Various definitions of the Laplace transform . . . . . 149 28.5 Appendix: Unit impulse function . . . . . . . . . . . . . . . . . . 150 28.6 Appendix: Care needed in using the one-sided LT . . . . . . . . 151 29 An example of application of loop-shaping in automobile control
153
V
155
Digital implementation of controllers
30 State space realization of transfer function 157 30.1 Example of state space realization . . . . . . . . . . . . . . . . . . 157 30.2 Numerical integration of state space model . . . . . . . . . . . . 161
VI
Pole-placement design using root locus
163
31 Unit step response of second order TF 31.1 State-space realization (concluded): Determining the of a transfer function to an arbitrary input . . . . . . . 31.2 Unit step response of standard second order TF . . . . 31.3 Appendix: MATLAB codes . . . . . . . . . . . . . . .
165 response . . . . . . 165 . . . . . . 167 . . . . . . 167
32 Second order TF with additional zero 32.1 Analysis of formulae developed for unit step response dard second order system . . . . . . . . . . . . . . . . . 32.2 The case of an additonal zero in a second order TF . . . 32.3 Appendix: MATLAB codes . . . . . . . . . . . . . . . .
171 of stan. . . . . 171 . . . . . 171 . . . . . 174
33 Second order TF with additional pole; recipe for pole-placement 33.1 Additional pole in second order TF . . . . . . . . . . . . . . . . 33.2 Recipe for placement of CL poles in a second order TF . . . . . 33.3 Root locus: What it is . . . . . . . . . . . . . . . . . . . . . . . . 33.4 Appendix: MATLAB codes . . . . . . . . . . . . . . . . . . . .
179 . 179 . 180 . 181 . 181
34 Root locus (Part 1) 183 34.1 Rules of sketching root locus . . . . . . . . . . . . . . . . . . . . . 183 34.2 Appendix: An alternative way to see Rule 2 . . . . . . . . . . . . 183 35 Root locus (Part 2)
185
36 Root locus (Part 3) 187 36.1 Roots of multiplicity greater than 1 . . . . . . . . . . . . . . . . . 187 36.2 Angles of departure and angles of arrival . . . . . . . . . . . . . 189 37 Root contours
191
38 Design using root locus 195 38.1 Appendix: TA09 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195 January 2, 2012
vii of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
VII
Lecture Notes
Addressing nonlinearities
197
39 Stability of nonlinear control systems: Circle criterion
199
40 Dealing with plant input constraints: An anti-windup scheme
203
VIII
205
PID Control
41 PID Controllers
207
42 Ziegler-Nichols tuning of PID controllers 42.1 What is controller tuning? . . . . . . . 42.2 What the two ZNT methods do . . . . 42.3 First method . . . . . . . . . . . . . . . 42.4 Second method . . . . . . . . . . . . . 42.5 Appendix: Further reading . . . . . .
January 2, 2012
viii of 218
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
209 209 209 210 212 212
Ramprasad Potluri
Preface These lecture notes adapt from the many good books available in the market the minimum set of topics that I have found necessary to practically build and operate a simple control system based on a 39-hour lecture course. These notes evolved as I taught over the past 5 years the course numbered EE250, titled Control Systems Analysis, which is a core course in the Department of Electrical Engineering, IIT Kanpur. These notes, along with some tutorial sessions on MATLAB and GNU Octave, and with background acquired in a basic electrical engineering course on motors and a basic electrical engineering course on signals and systems, are adequate to design and operate, for example, a permanent magnet dc motor control system. Indeed the set of experiments that we introduced in a course numbered EE380 in the Control Systems Laboratory in the Department of Electrical Engineering at IIT Kanpur in 2009, and that the students performed in 2010 and 2011 also, are entirely organized around dc motor control. The specific background needed from the course in signals and systems is mainly only a comfort with the Laplace transform and the convolution integral. The specific background needed from the course on motors is only a comfort with Maxwell’s laws and how these laws explain how a motor works. Students traditionally find certain topics difficult to understand. Those of these topics that I consider are not mandatory for a well-rounded introduction to feedback control, but are needed to see the complete story, are marked by 2 stars (**). I deemphasize these topics on the examinations. An example of this topic is the theorem-proof discussion of stability. On the other hand, certain topics that the student finds difficult, but are needed for a well-rounded introduction to controls, are marked by a single star (*). An example of such a topic is Nyquist stability criterion. The lectures on these topics go at a slower pace. The appendices of these notes shed more light on their respective lectures. These appendices usually contain material that I consider slightly advanced for the course or of only special interest, and are not discussed in the lectures. Each lecture note also contains some problems for practice. These problems are named Problem 1, Problem 2, etc. These problems are not discussed in lectures, and are in addition to the problems that the students see in the tutorials. Where needed, the appropriate literature is cited. While many good textbooks are available that may be relevant to this course, the textbooks that I cite are mainly those that are available in the Indian market. The topics of this course are organized as shown in Table 1. The last column of the table may or may not be replaced by state-space methods. The student may visit this table occasionally to navigate through the course. ix
How the location of poles and zeros in the s-plane determines the nature of response of a system, BIBO stability of LTI systems, Cauchy’s theorem of argument, Nyquist stability criterion
Theory of stability
Loop-shaping, performance specifications, lead and lag compensators, error constants, system type, examples of design using semilog graph papers
Analysis and design using Bode plots
State-space realization of a transfer function, Euler’s approximationbased discretization of a statespace equation, digital implementation, use of state-space realization to determine the response of a TF to arbitrary input
Digital implementation of controllers
PID controllers, two ZieglerNichols tuning methods
Tuning of PID controllers
Rules for sketching root loci, simple examples of design using root loci, standard second order model, MATLAB’s sisotool and rlocus, effect of additional zero and additional pole
Analysis and design using root locus
circle criterion, anti-windup
Addressing the nonlinearities in practical hardware
Table 1: How the topics of this course are organized. We will progress through the course from left to right along this table. Background
Negative feedback, LTI dynamic systems, Laplace transforms, transfer functions (TFs), state-space equations, block diagram algebra, signal flow graph algebra, Mason’s gain formula, virtues of a negative feedback control system, frequency response, Bode plots
Ramprasad Potluri
x of 218
January 2, 2012
Lecture Notes EE250 (Control Systems Analysis) IITK
Notation , < :=
“Is, by definition, equal to”. For example, ???????
x|
Transpose of the vector x.
a.k.a.
The set of all real numbers. “Denotes”. For example, y1 := y reads “y1 denotes y”, x˙ := dx/dt reads “x˙ denotes dx/dt”. Also known as.
xi
EE250 (Control Systems Analysis) IITK
January 2, 2012
xii of 218
Lecture Notes
Ramprasad Potluri
Part I
Preliminaries
1
Lecture 1
Introduction to the course This lecture begins with an overview of the syllabus and the preface.
1.1
About the title of the course
The title of this course is “control system analysis”. However, we will also learn “control system design”. Analysis involves evaluating how something works, and how well. Design involves building something to work as well as we desire. A control system is a mechanism by which a given plant works, in spite of disturbances, as desired by the designer of the control system. Under the title of control system, we are interested in this course in negative feedback control systems. Exist other kinds of control systems; see Figure 1.1. The terms control and feedback control are usually synonymous with negative feedback control. We define this last concept in the next section.
1.2
Definition of negative feedback control
To understand what a negative feedback control system is, we need to first understand what feedback and negative feedback are. Figure 1.2 states the goal of a control problem. The goal is that the output of a plant must follow a desired value zref as closely as possible. For example, if z is the actual speed of a motor’s shaft in rad/s, then zref is a number that represents the desired speed in rad/s. zref will need to be the input to, and z will be the output from, any system that will achieve the above goal. Therefore, we call the direction from zref to z as the forward direction, and that from z to zref as the backward direction. Figure 1.3 shows how we can achieve the said goal using the concept of negative feedback control. We use a sensor that outputs a quantity y that is proportional to z, we take the negative of y, and we sum the quantity −y with yref . The last is the output of a prefilter. The prefilter and the sensor are chosen such that yref and y are dimensionally consistent and comparable in magnitude. In the example of the motor, yref and y are voltages, both of which may be in the interval [0, 3] V. As y is taken from the output and fed towards the 3
EE250 (Control Systems Analysis) IITK
Lecture Notes
Control Systems (CS)
Open-loop CS (E.g., Household fan)
Closed-loop CS (also known as feedback CS)
Automatic CS
Negative feedback CS
Human-in-the-loop CS (E.g., driver and automobile)
Positive feedback CS
Figure 1.1: A classification of control systems. In this course, we are interested in negative feedback control systems.
input, we say that y is fed back. The operation of taking −y is called negative feedback. The act of using the error e = yref − y to apply an input u to the plant so that e is reduced is called negative feedback control. A control system built using the concept of negative feedback control is called negative feedback control system. is called summator or sumThe element represented by the shaded circle mer or summing element. If its inputs are a and b, then its output is the sum a + b. given that the summator performs the operation of summation, the system of Figure 1.3 may be called a positive feedback control system if y were to be fed back without a minus sign, that is, if e = yref + y were the output of the summator. Feedback action can be contrasted with feedforward action. This distinction is explained in lectures 6 and 7.
1.3
Examples of control systems
The applications of control systems are numerous. Listed below are some examples of where control systems are used in practice. 1. The read head of hard disc drive needs to be positioned at a desired track, or needs to remain at the track in spite of any mechanical shocks that the PC/laptop may receive. zref
u
Plant
z
Figure 1.2: Plant and goal. We have a plant whose input is u and output is z. We want z to follow, as closely as possible, a quantity zref that is dimensionally consistent with z.
January 2, 2012
4 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
2. The write head of an inkjet printer needs to be positioned at a desired location with a tolerance of a few micrometers. 3. The speed of the shaft of a motor needs to be maintained at a given value in spite of voltage fluctuations. 4. In automobiles: antilock braking systems (ABS), cruise control systems, active suspension, auxiliary steering, etc. 5. Aircraft autopilot. 6. Rockets need to be controlled so that they track a desired trajectory on their ascent in spite of admissible wind gusts and variations in atmospheric pressure. 7. The article [1] describes the many ways in which control theory is used in digital imaging systems manufactured by Xerox.
1.4
Videos of control systems
www.youtube.com has many nice videos of control systems. Listed below are some of these. All these videos existed on January 02, 2012. 1. Magnetic levitation. The concept of magnetic levitation is used in maglev, as well as in magnetic bearings. The idea of using a control system is to keep the magnet floating above/below another magnet at a desired distance. http://www.youtube.com/watch?v=kXodf7WKiFs. 2. Sojourner. The 6 wheels need to be driven and steered such that the rover moves along the desired path with the desired speed. http://www.youtube.com/watch?v=zZWOGcdC_PI&feature=fvsr. 3. Inverted pendulum. Keep a pole upright on a movable cart. Simulates pole balanced in the hand, and rocket launching. http://www.youtube.com/watch?v=Ci_y14y3DU4. zref
Prefilter
yref +
e = yref − y
− y
Controller
u
Plant
z
Sensor
Figure 1.3: A block diagram representing the standard architecture of a negative feedback control system (NFCS). This NFCS is one way to make z follow zref as closely as possible. Underlying this control system are the concepts of negative feedback and negative feedback control. zref is called reference input or desired output or command input, z is called actual output or regulated output, y is called sensor output, e is called error, u is called control input. The block diagrams seen in control system theory satisfy two conditions. First, each block conveys signals unidirectionally. Second, each block is non-loading, that is, it does not load the preceding block.
January 2, 2012
5 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
4. Two-wheeled inverted pendulum. This setup is a precursor to Segway. http://www.youtube.com/watch?v=DDBkg1PfvFY. 5. Segway. This setup is an advanced inverted pendulum. http://www.youtube.com/watch?v=BD-y8F1Fa7g&feature=fvst. 6. Ball on wheel. This setup is a cute one, though I have not yet imagined an application for it. http://www.youtube.com/watch?v=F9CuHVAp5k0. 7. Dual-motor ball-beam (DMBB). We developed this setup in the Networked Control Systems Laboratory1 . One application for this setup is as an inexpensive testbed for dual-manipulator coordination. If a certain analogue of a control algorithm that coordinates two manipulators to handle a load together will not run successfully on the DMBB testbed, then most likely it will not work on the actual manipulators. Details on dual-manipulator coordination can be found in [2, 3, 4, 5]. http://www.youtube.com/watch?v=2KhvPT6935k and http://www.youtube.com/watch?v=-tB6_CCJu_8.
1.5
Appendix: History of control systems
Most textbooks devote space in their first chapters to a discussion of the history of control systems. Besides textbooks, I found the discussion in [6, 7] to be especially interesting and different. Prof. Chris Bissell of the Open University of UK has written many easy-to-read articles on the history of control. These articles may be available for free download on the web. The explanation provided in Section 1.2 for negative feedback is consistent with Harold Black’s statement that may explain why he called his amplifier the negative feedback amplifier: “Then came the morning of Tuesday, August 2, 1927, when the concept of the negative feedback amplifier came to me in a flash while I was crossing the Hudson River on the Lackawanna Ferry, on my way to work. . . . I suddenly realized that if I fed the amplifier output back to the input, in reverse phase, and kept the device from oscillating (singing, as we called it then), I could have exactly what I wanted: a means of canceling out the distortion in the output.” [8]. The section in [9] titled Two cultures of feedback and stability describes the origins of the term negative feedback in more detail, and states the modern meaning of negative feedback. The explanation of Section 1.2 is consistent with this meaning. Harold Black, Hendrik Bode, and Harry Nyquist of the Bell Laboratories are credited by many textbooks for having laid the foundations of Control theory. Black is credited for having invented the negative feedback amplifier, while Bode plots and Nyquist stability theory are stuff that we will see in this course. The article [9] describes the historical context in which these men made their contributions.
1 Western
Labs 217 B, IIT Kanpur.
January 2, 2012
6 of 218
Ramprasad Potluri
Lecture 2
Dynamic systems, transfer functions, and block diagrams In this course, we will work with dynamic systems described only by linear time-invariant (LTI) differential equations. Classical control theory uses the Laplace transform to represent such systems through transfer functions. We explain these two concepts in this lecture.
2.1
Dynamic systems
A dynamic system is one whose state at the next time instant is dependant on the state at the current instant, and on the stimulus at the current instant. Traditionally, differential and difference equations, as well as algebraic equations, are used to describe dynamic systems. There are two types of differential/difference equations that describe a dynamic system. Equations such as m x¨ = F, the relation between displacement x and force F causing this displacedi ment, and L dt = v, the relation between the EMF v across an inductor and the current i that this EMF drives, are dynamic equations. Equations such as x˙ = v, the relation between displacement x and velocity v), x (1) = x (0) + v(0)∆t (relation between current position x (0), current velocity v(0), next position x (1), and time step ∆t), and q˙ = i (relation between charge q and current i), are kinematic equations. We see that not every differential equation is a dynamic equation. Only those differential equations that include the “causes” or “stimuli” are dynamic equations. Figure 2.1 summarizes the discussion of this section. Note that the terms dynamic equations and kinematic equations are not be confused with dynamics and kinematics, which are branches of mechanics1 . 1 From Oxford English Dictionary: mechanics, n. b. The branch of applied mathematics that deals with the motion and equilibrium of bodies and the action of forces, and includes kinematics, dynamics, and statics. Now often distinguished as classical mechanics (as opposed to quantum mechanics).
7
EE250 (Control Systems Analysis) IITK
Lecture Notes
Dynamic system
Algebraic equations
Differential/Difference equations
Dynamic equations
Kinematic equations
Figure 2.1: The types of equations that describe a dynamic system [].
V (s) +
− E(s)
1 sL a + RΣ
I (s)
Kt
T (s) +
TL (s) −
1 Js + B
ω (s)
Kb
Figure 2.2: The block diagram of a PMDC motor.
2.2
Transfer functions, block diagrams
We will illustrate the development of transfer functions (TF) using the example of a DC motor. A dc motor is modeled by the equations listed in Table 2.1. Four of these equations describe a permanent magnet dc (PMDC) motor. Table 2.2 shows these equations and the corresponding TFs obtained using the Laplace transform. Using these TFs, we can draw the block diagram of Figure 2.2. The advantage of working with a block diagram over working with equations is that a picture speaks a thousand words. Problem 1 Identify the algebraic and dynamic equations in the tables.
Problem 2 How will the motor block diagram change if you wish to include the kinematic equation θ˙ = ω?
Problem 3 Table 2.2 shows the four elemental TFs of a PMDC motor. Work out the composite TF of this motor from the input V (t) to the output ω (t).
January 2, 2012
8 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Table 2.1: The equations that describe how a dc motor works. Equation
Explanation
di + Ra i + V (t) = L a dt E(t)
The voltage applied to the armature pushes a current through the armature. A part of this voltage drops across the resistance R and inductance L. The balance of this voltage overcomes the back emf E(t) that is induced by the rotation of the rotor.
TM (t) = K f Φ(t)i (t)
In a dc motor, the interaction between the current i flowing through the conductors of the motor’s armature and the magnetic flux Φ produces torque TM .
TM (t) = K M i (t)
In a PMDC motor Φ is constant, resulting in K f Φ being a constant denoted K M .
TM (t) − TL (t) − Bω = J dω dt
The torque TM accelerates the rotor after overcoming any load or disturbance torques TL that may act on the motor shaft and the torque Bω due to viscous friction. B is known as the coefficient of viscous friction.
E ( t ) = Ke Φ ( t ) ω ( t )
In a dc motor, the armature’s motion through the field induces an EMF in the conductors of the armature.
E(t) = KE ω (t)
In a PMDC motor, as Φ is constant, Ke Φ is a constant denoted KE .
Table 2.2: The transfer functions of a PMDC motor. Time domain equation
Equation in s domain
Transfer function
di V (t) − E(t) = L a dt + Ra i
V (s) − E(s) = L a sI (s) − L a i (0) + R a I ( s )
I (s) V (s)− E(s)
TM (t) = K M i (t)
TM (s) = K M I (s)
TM I (s)
TM (t) − TL (t) − Bω = J dω dt
TM (s) − TL (s) − BΩ(s) = JsΩ(s) − Jω (0)
E(t) = KE ω (t)
E(s) = KE Ω(s)
January 2, 2012
9 of 218
=
= KM
Ω(s) TM (s)− TL (s) E(s) Ω(s)
1 sL a + R a
=
1 sJ + B
= KE
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
i S
+
V −
N
Figure 2.3: Illustration for a conductor in a magnetic field. The applied EMF V drives a current i through the conductor. The conductor turns about its axis due to the interaction of i and the field. A nice dynamic picture that explains the principle of working of DC motors can be found at http://www.sciencejoywagon.com/physicszone/lesson/ otherpub/wfendt/electricmotor.htm
2.3
Appendix: Detailed model of a DC motor
When a conductor is placed in a magnetic field, and a current is sent through the conductor, then the conductor experiences a thrust/force called Lorentz force. The direction of this force relative to the directions of the field and current is given by Fleming’s left hand rule which says that if the first finger, second finger, and thumb were to be held mutually perpendicular, and the direction of the field (from north pole to south pole is) indicated by the first finger, and the direction of the current is indicated by the second finger, then the direction of the thrust is indicated by the thumb: field current thrust
first finger second finger thumb
This phenomenon is utilized in motors: the current carrying conductor is mounted on an axis and the thrust that the conductor experiences works as a turning force (torque) that makes it rotate around the axis (Figure 2.3). If the field is variable (as, for example, when using an electromagnet), then we can write the torque as follows: TM (t) = K f Φi (t) where K f is a constant of proportionality, i (t) is the current in the conductor and Φ is the flux. If a permanent magnet is used to provide the field, then the field will be constant, and the motor’s torque can be written as TM (t) = K M i (t)
(2.1)
Another phenomenon comes into play. When the conductor turns in the field, it cuts through the field, and an EMF is induced in it that is, by Faraday’s January 2, 2012
10 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
law, directly proportional to the rate of change of flux linkage (which in turn is proportional to the flux Φ and the angular velocity of the conductor ω): E(t) = Ke Φω (t) If Φ is a constant, as when produced by permanent magnets, then E(t) = KE ω (t)
(2.2)
By Lenz’s law, the direction of this induced EMF is such that it tends to oppose the flux change causing it. So, the current that eventually flows through the conductor is the result of the net EMF V − E, where V is the voltage applied to the ends of the conductor as shown in Figure 2.3. This current causes a voltage drop across the inductance L a and the resistance R a of the conductor: V (t) − E(t) = L a
di + Ra i dt
(2.3)
In the present lecture, we are interested in DC motors. Equations (2.1), (2.2), (2.3) describe the operations of DC motors. A DC motor has not just one turn of wire for a conductor, but many turns of wires, all distributed evenly around a cylindrical core. This set of wires is called armature. The core is mounted on an axis. The entire rotating component (core plus conductors) of the motor is called rotor. The L a and R a in Equation (2.3) are in fact the net inductance and the net resistance of the armature, not of a single conductor. In order for the rotor to spin, the torque T generated by any motor (AC or DC) has to overcome any externally applied disturbance torques (TL ), the moment of inertia (J, assumed constant, although can be variable too) of the motor (which matters when there is a change in angular velocity ω), and the friction in the bearings as well as due to wind (summed up in B, the coefficient of viscous friction), summarized by the fundamental torque equation: T = TL + J
dω + Bω dt
In the case of the DC motor that we have been discussing, this equation can be written as follows: dω TM (t) − TL (t) = J + Bω (2.4) dt In case a load is mounted on the rotor shaft, then the moment of inertia and the coefficient of viscous friction of the motor rotor-load combination will b be different from J and B mentioned above and will be equal to some b J and B which are called respectively equivalent moment of inertia (of the rotor-load combination) referred to the rotor shaft and equivalent coefficient of viscous friction (of the rotor-load combination) referred to the rotor shaft.
January 2, 2012
11 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
12 of 218
Lecture Notes
Ramprasad Potluri
Lecture 3
Transfer functions and block diagrams revisted, state-space equations, an example of design in state space In Lecture 2, we saw the concepts of transfer functions (TFs) and block diagrams. In this lecture, we make a few observations about these concepts. We see block diagram algebra in a later lecture. We also give a brief overview of state-space theory. In this theory, ordinary differential equations (ODEs) are written as state-space equations.
3.1
Transfer functions revisited
Classical control theory depends on the concept of transfer function (TF). The TF of an LTI system is defined as the ratio of the Laplace transform of the output of this system to that of the input to this system. Note that this definition ignores the initial conditions (ICs). In classical control theory, the only properties of a control system that are of interest are 1. the stability of the control system, 2. the system’s speed of response, measured in terms of properties of the system’s response to a unit step input such as rise time tr and settling time ts , 3. the system’s tendency to oscillate, measured by the peak overshoot M p in the system’s response to a unit step input, 4. how closely the output tracks the input once the transient phase of the response is over, measured by the steady-state error ess to a unit step input, and 5. the robustness of the above properties. 13
EE250 (Control Systems Analysis) IITK
Lecture Notes
For LTI systems modeled by ODEs, the satisfaction of these properties by the system is independent of ICs []. As classical control theory restricts itself to such systems, there is no loss in ignoring ICs. On the other hand, for nonlinear systems, these properties are in general IC-dependant. Classical control theory is inadequate here. [10, Page 42] provides additional interesting observations about TFs. Problem 4 Draw a block diagram of a PMDC motor to include the initial conditions. Do you see any advantages in such an inclusion?
3.2
Block diagrams revisited
Transfer functions are defined only for LTI systems represented by ODEs of the form an
d n −1 y dy dn y + a n −1 n −1 + · · · + a 1 + a 0 y = n dt dt dt dm u du d m −1 u bm m + bm−1 m−1 + · · · + b1 + b0 u, dt dt dt
(3.1)
and by LTI algebraic equations. An example is the model of the PMDC motor seen in Lecture 2. In the block diagrams of such systems, only two blocks appear. These are the summator and the gain block shown in Figure 3.1. For Equation 3.1, the gain from u(t) to y(t) is the same as the TF from u(t) to y(t), and is the ratio of Laplace transform of the output to the Laplace transform of the input as Y (s) bm sm + bm−1 sm−1 + b1 s + b0 = U (s) a n s n + a n −1 s n −1 + · · · a 1 s + a 0 The Laplace transform handles the product of two time-domain signals as follows [11, Page 837]:
L{ f 1 (t) f 2 (t)} =
Y1 (s) +
1 2πj
Z c+ j∞ c− j∞
F1 (ξ ) F2 (s − ξ )dξ =
Y (s) = Y1 (s) + Y2 (s)
1 F (s) ∗ F2 (s). 2πj 1
U (s)
Gain
Y (s)
G (s)
+ Y2 (s)
Figure 3.1: Only two blocks appear in the block diagrams of systems modeled by LTI ODEs and by LTI algebraic equations. These are the summator and the gain block. Figure 2.2 shows a sample block diagram.
January 2, 2012
14 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Thus, we may wonder why we did not apply it to the nonlinear element K
TM (t) = K f Φ(t)i (t), obtain the s-domain expression TM (s) = 2πjf Φ(s) ∗ I (s), and somehow incorporate this term in the block diagram of Figure 2.2. Note that, once included in the block diagram, we may wonder how to work with such a term. As we see in the remainder of the course, even working with systems whose block diagrams have only summators and gain blocks is sufficiently challenging. Thus, classical control theory restricts itself to these systems. Note that apart from LTI ODEs such as that of Equation 3.1, classical control theory also works with LTI ODEs containing time delay elements. The TFs of these elements too are gain blocks. Problem 5 Prove, starting from the definition of linearity, that the equation TM (t) = K f Φ(t)i (t), representing a system with inputs Φ and i and output TM , is a nonlinear one.
Problem 6 Consider the following two examples: 1. In a thermal power plant, sensing of temperature x (t) in the path of superheated steam is done at a distance from where it is actually needed to be sensed. This introduces a delay in sensing. Thus, as a first approximation, we can say that what the controller receives is a delayed version, y(t) of x (t). Of course, more accurately speaking, even the actual temperature sensed downstream may be slightly different from even the delayed version of x (t). 2. In a distributed control system, where the sensors, actuators, and controllers are interconnected by a “control network” such as Controller Area Network (CAN), even though the sensor may have sensed the variable of interest x (t) immediately, there may be a delay introduced by the network in this sensed value reaching the controller. Thus, the controller receives y(t) which is a delayed version of x (t). In such cases, y(t) is a time-shifted version of x (t) as shown in Figure 3.2. This behavior is summarized by the equation y(t) = x (t − td )1(t − td ). Applying the Laplace transform on both sides of this equation gives Y (s) = X (s)e−td s . Now consider the TF Y (s) s + a −td s = e . U (s) s+b Write the corresponding ODE.
3.3
State-space equations
We call state-space equations of LTI systems LTI state-space equations, and the control theory that uses these equations LTI state-space control theory. January 2, 2012
15 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
x (t)
t
y(t)
td
t
Figure 3.2: Illustration for time delay. y(t) is a delayed version of x (t).
As an example to illustrate the state-space framework, consider d2 y dy + a1 + a0 = u. dt dt2
(3.2)
This equation is a linear time invariant (LTI) second order ordinary differential equation in one variable (y). This equation models a certain plant. Denoting ˙ we write this differential equation equivalently as y1 := y and y2 := y˙ 1 = y, y˙ 1 = y2 , y˙ 2 + a1 y2 + a0 y1 = u, which is a system of first order differential equations in two variables (y1 and y2 ). This system is written as " # " #" # " # " # h i y y˙ 1 0 1 y1 0 1 = + u, y = 1 0 . y˙ 2 − a0 − a1 y2 1 y2 The vector z = [y1 y2 ]| is called the state vector of the plant that this second order ODE models. Thus, z˙ = Az + Bu, y = Cz, with
" A=
0
1
− a0
− a1
#
" , B=
0
#
1
, C=
(3.3)
h
1
0
i
,
is the state-space equation (a.k.a. state-space model) of this plant. Rewriting the original ODE in the state-space form helps as follows. Note that A contains the crucial information about the natural response of the original second order ODE. If we construct u as u = −kx + yref , January 2, 2012
16 of 218
(3.4) Ramprasad Potluri
EE250 (Control Systems Analysis) IITK we have
ˆ + Byref , y = Cz, z˙ = Az
Lecture Notes
(3.5)
where Aˆ = ( A − Bk ). Thus, the state feedback control of Equation (3.4) helps us replace the inputoutput behavior given by Equation (3.3) with the input-output behavior given by Equation (3.5). As A is 2 × 2 and B is 2 × 1, k is 1 × 2. With k = [k0 k1 ], we have " # 0 1 Aˆ = . −( a0 + k0 ) −( a1 + k1 ) ˆ the state-space model of Equation (3.5) is equivalent to the LTI With this A, ODE dy d2 y + ( a1 + k 1 ) + ( a0 + k 0 ) = u dt dt2 This k can be constructed algebraically without solving any ODEs or statespace equations.
3.4
An example of design in the state-space domain
In the LTI ODE of Equation (3.2), let a0 = 1, a1 = 0.05. For a unit step input, that is, for u(t) = 1(t), where 1(t) is a unit step function, a.k.a., Heaviside function, ( 0, t < 0 1( t ) = 1, t > 0, the output of Equation (3.2) is shown in Figure 3.3. This response is highly oscillatory, and it overshoots almost 100% the value of 1 that it is required to settle to. Hence, this response is undesirable, and needs to be changed. Given that we cannot modify the plant, we decide to use state feedback to create a new system around this plant such that the input-output behavior of this new system is desirable. From prior experience with the solution of the equation d2 y dy + 2ζωn + ωn 2 = ωn 2 u, dt dt2
(3.6)
which is more general than Equation (3.2), we know that, for a given ωn , the larger the value of ζ, the less the overshoot and oscillations in the unit step response √of Equation (3.6). In the present example, ωn = 1, and ζ = 0.025. ζ = 1/ 2 with ωn = 1 give us an adequately damped response, that is, one which has an overshoot of not more than about 10 − 15% and has only about 3 – 5 overshoots and undershoots before it settles to within about 95% of the desired value of 1, as shown in Figure 3.4. ζ is called the damping ratio and ωn is called the natural frequency of Equation (3.2). √ 2 and We use this prior experience to decide that we want a + k = 1 1 √ a0 + k0 = 1. Thus, we obtain k0 = 0 and k1 = 2 − 1. In the initial stages of developing this experience, we may need to go frequently from the design in the state-space domain to solving the resulting state-space equation in the time domain, and back, to see how well the design January 2, 2012
17 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Unit step response of d2y/dt2 + 0.05 dy/dt + y = u 2 1.8 1.6 1.4
y(t)
1.2 1 0.8 0.6 0.4 0.2 0
0
50
100
150
200
250
Time (sec)
Figure 3.3: The unit step response of the plant of Equation (3.2) with a0 = 1 and a1 = 0.05. This response is poorly damped as seen from the fact that it has large overshoot and is very oscillatory.
performs in practice. With experience, however, the design can be performed without needing to see how the resulting system performs in the time domain. Problem 7 Develop the state-space model of the PMDC motor.
Problem 8 Paste the following code at the command prompt of either MATLAB or GNU Octave, and hit the enter key. num = 1; den = [1, 0.005, 1]; sys = tf(num,den); step(sys); Without exiting the software, paste the following code at the command prompt and hit the enter key. title(’Unit step response of d^2y/dt^2 + 0.05 dy/dt + y = u’); Observe the results. Modify the code to obtain the responses shown in figures 3.3 and 3.4.
January 2, 2012
18 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Unit step response of d2y/dt2 + 1.414 dy/dt + y = u 1.4
1.2
1
y(t)
0.8
0.6
0.4
0.2
0
0
1
2
3
4
5 Time (sec)
6
7
8
9
10
Figure √3.4: The unit step response of the plant of Equation (3.2) with a0 = 1 and a1 = 2. This response is well-damped as seen from the fact that its overshoot is low and it has very few overshoots and undershoots.
January 2, 2012
19 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Acknowledgment: The discussions with Mr. Manavaalan Gunasekaran, Ph.D. student, helped formulate the ideas presented in Section 3.2.
January 2, 2012
20 of 218
Ramprasad Potluri
Lecture 4
An example of design in the classical domain, block diagram manipulations 4.1
An example of design in the classical control domain
Let us see how classical control theory would solve the design problem of Section 3.4. Figure 4.1 shows the block diagram of a control system that we would want to build around the given plant so that the closed-loop behavior is to our satisfaction. We see that if we select C (s) as equal to c1 s + c0 , then the closed-loop TF is Y (s) c1 s + c0 = 2 Yref (s) s + ( a1 + c1 ) s + ( a0 + c0 )
(4.1)
We can plot the unit step response of this TF in GNU Octave or MATLAB using the code c0=30; c1=5; a0=1; a1=0.05; num=[c1,c0]; den=[1,a1+c1,a0+c0]; step(tf(num,den)); Through trial and error using this code, we determine that for the given value of a0 = 1 and a1 = 0.05, the values c0 = 30 and c1 = 10 provide us a response that is satisfactory in the sense described in Section 3.4. This response is shown in Figure 4.2. Classical control theory says that if the zero of this TF s = −c0 /c1 is deep in the left half plane in comparison to the poles of this TF, then the effect of this zero will be negligible on the response of this TF to an input, and the nature of the response will be determined by only the poles of this TF. Let us see if the zero is much deeper than the poles with the above-selected values of c1 and c0 . The following code provides us the values of these zero and poles. c0=30; c1=5; a0=1; a1=0.05; num=[c1,c0]; den=[1,a1+c1,a0+c0]; -c0/c1, roots(den) 21
EE250 (Control Systems Analysis) IITK Yref (s) +
C (s)
−
U (s)
Lecture Notes 1 s2 + a1 s + a0
Y (s)
Figure 4.1: Block diagram of the proposed control system that we wish to design in the classical control domain.
We see that the zero is at s = −3 and the poles are at −5.0250 + j2.3978 and −5.0250 − j2.3978. We see that even though the poles are deeper than the zeros, the response is still to our satisfaction. We explore the answer to this behavior in the lectures on root locus. Step Response 1.4
1.2
Amplitude
1
0.8
0.6
0.4
0.2
0
0
0.2
0.4
0.6 0.8 Time (sec)
1
1.2
1.4
Figure 4.2: The unit step response of the closed-loop system of Equation (4.1) with a1 = 0.05, a1 = 1, c0 = 30, and c1 = 10. This response is well-damped as seen from the fact that its overshoot is low and it has very few overshoots and undershoots.
4.2
Rules of block diagram algebra
In our course, as we work with block diagrams, we will frequently have the need to quickly determine the TF between a certain input and a certain output. January 2, 2012
22 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
As mentioned in Lecture 2, it is easier to work with block diagrams than with the underlying equations. Figure 4.3 lists the few simple rules that enable us to manipulate block diagrams.
4.3
Examples
Example 1 An example of block diagram manipulation.
Example 2 January 2, 2012
23 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Figure 4.3: Rules of block diagram algebra.
January 2, 2012
24 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Example 3
We can quickly understand the block diagram algebra through just a few more examples. Note that this algebra was possible because we have restricted our attention to systems that are LTI.
4.4
Appendix: A little more on state-space theory
State-space theory works with differential equations that are of the general form x˙ = f ( x, u, t). This equation is called state-space equation. Here, x˙ , dx/dt, with t the time. x and f are real-valued vectors with n components each, and u is a real-valued vector of m components. In mathematics, this last statement is compactly written as x, f ∈
4.5
Appendix: Similarities and differences between classical and LTI state-space control
Classical control theory and LTI state-space control theory work with LTI ODEs and design controllers without solving the underlying differential equations. This thought is illustrated in more detail in Figure 4.4. The difference is that, while classical control theory designs controllers using the Laplace transform and transfer functions, state-space control theory designs in the state space. Further, while classical control uses graphical techniques for analysis and design, state-space control theory uses linear algebra. Additionally, while classical control can handle time delay in the ODEs, LTI state-space theory cannot [].
January 2, 2012
25 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
LTI ODE
TF or
Analysis
LTI SS
&
form
Design
Lecture Notes
LTI ODE
Simulation
No
simulation successful?
Yes
Exit
Figure 4.4: Illustration for the similarities between classical control theory and LTI state-space control theory.
January 2, 2012
26 of 218
Ramprasad Potluri
Lecture 5
Signal flow graphs Manipulation of block diagrams as demonstrated in Lecture 4 can be tedious in the case of complex block diagrams. Example 4 Example of a complex block diagram.
The difficulty in manipulating this block diagram arises from the fact that we may be forced to move the summators through the branches that are to their right. See Rules 10 in table 1-3 of Lecture 4 about transfer of a summator into branches.
For a block diagram that looks complicated, we can determine the inputoutput gain by first drawing a signal flow graph (SFG) that is equivalent to this block diagram, and then applying a formula that has come to be known as Mason’s gain formula (MGF). We see SFGs in this lecture and MGF in the next one.
5.1
Definitions
Consider the SFG of Figure 5.1 Here are definitions of the various parts of this SFG. Nodes represent signals. For example, x1 , x2 , x3 , x4 . Transmittances / gains represent signal multipliers. For example, a, b, c, d. Branches are directed line segments, each of which joins two nodes. Arrow on a branch shows direction of signal flow. E.g. x1 to x2 , x2 to x3 , x4 to x3 , x3 to x2 , x3 to x3 . Input node / source has only outgoing branches. For example, x1 , x4 . 27
EE250 (Control Systems Analysis) IITK
Lecture Notes
Figure 5.1: A sample signal flow graph.
Output node / sink has only incoming branches. For example, x3 . Mixed node has both incoming and outgoing branches. For example, x2 , x3 . A mixed node can be converted into an output node by adding an outgoing branch with a gain equal to 1. A path is the traversal of connected branches in the direction of the arrows. An open path is a path in which no node is crossed more than once. For example, x1 − x2 − x3 , x4 − x3 , x1 − x2 , x4 − x3 − x2 . A closed path is a path that ends at same node where it began and doesn’t cross any other node more than once. For example, x2 − x3 − x2 , x3 − x2 − x3 . A path that is neither open nor closed crosses some node more than once but ends at a different node from which it began. For example, x1 − x2 − x3 − x2 , x4 − x3 − x2 − x3 , etc. A loop is a closed path. For example, x2 − x3 − x2 , x3 − x2 − x3 . Loop gain is the product of branch gains of a loop. For example, bc. Nontouching loops do not possess common node. We don’t have any in this example. Forward path is a path from input node to output node that does not cross any node more than once. For example, x1 − x2 − x3 − x3 , x4 − x3 − x3 . Forward path gain is the product of branch gains of a forward path. For example, ab, d, dc. Example 5 Draw an SFG that is equivalent to the following block diagram.
Answer:
From Example 5, we see that a node may either represent a signal or the sum of two or more signals. In this example, we have the nodes corresponding to the signals C and R, the node corresponding to the sum E1 of signals N and GE, and the node corresponding to the sum E of signals R and − HC. January 2, 2012
28 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
5.2
Lecture Notes
SFG Algebra
Similarly as in block diagrams, we can list about 13 rules for manipulations of SFGs. In fact, as block diagrams are same as SFGs, the 13 rules we saw for block diagrams also apply to SFGs. However, we don’t need to build that list here all over, as we will not use SFGs in our course — we will only use block diagrams. However, we will try to develop a little familiarity with SFG algebra through an example. Example 6 Let us work out the following five rules that we choose purely arbitrarily.
5.3
Appendix: Usage of these terms
An SFG conveys exactly the same information as a block diagram. In fact, a search on the Web for “signal flow graph” shows that even block diagrams are sometimes known as SFG in the European literature. Our usage of the term SFG is in keeping with the US literature. The difference between block diagrams and SFGs is only in the pictorial representation and the domains where each is used. Block diagrams are primarily used in the control systems community. SFGs are used in the communications community.
January 2, 2012
29 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
30 of 218
Lecture Notes
Ramprasad Potluri
Lecture 6
Mason’s gain formula, properties desired from a control system Samuel J. Mason of MIT published a paper each in the years 1953 and 1956 that resulted in “Mason’s gain formula”. Mason developed his formula for signal flow graphs (SFG).
6.1
Mason’s gain formula
The gain from the input R to the output Y in an SFG is given by the following formula: Y 1 N = ∑ Pk ∆k R ∆ k =1 Here, • ∆ is called “graph determinant” of the given SFG, and is given by the following expression:
∆ = 1 − ∑ loop gains of all loops taken one at a time
+ ∑ products of loop gains of all non-touching loops taken two at a time
− ∑ products of loop gains of all non-touching loops taken three at a time + ... • N is the number of forward paths from R to Y. • Pk is the kth forward path gain. • ∆k is the ∆ for that part of the SFG that doesn’t touch the kth forward path. 31
EE250 (Control Systems Analysis) IITK
Lecture Notes
Figure 6.1: A sample signal flow graph for which to find the gain.
Example 7 Find x5 /x1 for the SFG of Figure 6.1. 1 x5 = x1 ∆
N
∑ Pk ∆k
k =1
with N=2 ∆ = 1 − ( G23 G32 + G23 G34 G42 + G23 G34 G45 G52 + G23 G32 G44 + G23 G35 G52 G44 ) P1 P2 ∆1 ∆2
+ ( G23 G32 × G44 + G23 G35 G52 × G44 ) = G12 G23 G34 G45 = G12 G23 G35 = 1−0 = 1 = 1 − G44
6.2
A note on loop gain
In classical control theory, we use the term loop gain frequently. This term is defined slightly differently for SFGs and block diagrams. Figure 21.4 reproduces the block diagram and the corresponding SFG from Example 5. For an SFG, the loop gain for a given loop is defined as the product of branch gains of the loop. By this definition, in Figure 21.4, the loop gain is − G (s) G2 (s) H (s), and the gain from R(s) to C (s) is G (s) G2 (s) (1 + G (s) G2 (s) H (s)) . For a block diagram the loop gain for a given loop is defined as the product of the gains around the loop ignoring the signs at the summator. By this definition, in Figure 21.4, the loop gain is G (s) G2 (s) H (s), and the gain from R(s) to C (s) is G (s) G2 (s) (1 + G (s) G2 (s) H (s)) . We note that this last expression has accounted for the signs at the summator. January 2, 2012
32 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Figure 6.2: A block diagram and the corresponding signal flow graph.
We often need to be able to quickly write the gain of a single-loop block diagram from an input X to an output Y. We can write this gain easily as Gain of the block diagram from X to Y =
Gain of the forward path from X to Y 1 + loop gain of this block diagram
Applying this formula, we can easily write C (s) G1 (s) G2 (s) = , R(s) G1 (s) G2 (s) H (s) C (s) G2 (s) = . N (s) G1 (s) G2 (s) H (s)
6.3
Properties desired from a control system
A control system (CS) needs to have the following properties: • Stability: This is the most fundamental property desired of any CS. • Input amplitude constraints u(t) should remain within certain bounds so that the approximation that the plant is linear remains valid. r (t)
K (s)
u(t)
G (s)
y(t)
−
• Disturbance rejection • Noise filtering • Low sensitivity to parameter variations • Good performance — both in transient and in steady state.
January 2, 2012
33 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
34 of 218
Lecture Notes
Ramprasad Potluri
Lecture 7
Virtues of control systems: Disturbance rejection We may decide that the final CS may be a feedforward CS or a combined feedforward feedback CS. Example 8 Disturbance rejection One of the main purposes of a CS is to obtain the desired behavior from the plant/object in the face of disturbances that may be acting on the object. A negative feedback CS helps in disturbance rejection as follows: Plant
η
+
yd
+
C
P
−
y
+
Here, C and P are gains for amplification factors of the controller, plant, sensor respectively. The goal is to minimize the effect of η on y. Let us see if the feedback scheme shown achieves this goal.
Y (s) = η (s) + PCE(s)
=⇒ Y (s) = η (s) + PC (Yd (s) − Y (s)) =⇒ Y (s)(1 + PC ) = η (s) + PCYd (s) η (s) PC =⇒ Y (s) = + Y (s) 1 + PC 1 + PC d if PC 1, then we see that η appears in y with a strength that is PC times less than that with which yd appears as y.
Example 9 Disturbance feedforward. 35
EE250 (Control Systems Analysis) IITK
Lecture Notes
In Example 8, we see that η will need to be inside the control loop for it to be attenuated. That is, it will need to be sent around the loop, and it will need to create the error e between yd and y before its effect can be minimized. Is it possible to reject the effect of η without needing to first create an error? In those cases, where η is measurable (by a sensor), we can reduce its effect as follows: Gff
+ yd
C
η Plant P
+
y
+
From the block diagram, we have Y (s) = η (s) + PC (Yd (s) + Gff η (s)) That is, Y (s) = (1 + PCGff )η (s) + PCYd (s) We want Y (s) = Yd (s). We may choose Gff such that PCGff = −1, that is, Gff = −1/( PC ). This is useful when the disturbance η may be sensed almost as quickly as it hits the plant.
Where is the feedforward action in the figure of Example 9? In fact, which is the “forward” direction? In that example, it was probably difficult to see that the forward direction is from the input side, and that the backward direction is from the output to the input. But if adopt the convention that all inputs come from left, and all outputs exit from the right, thus Inputs
Outputs System
agree that the forward path is from input to output, and redraw the block diagram in Example 9 thus η
Gff Plant
+ yd +
C
P
+
y
+
then, it is clear that the disturbance input is being injected along an alternative forward path into that system. On the other hand, in Example 8, the output was being fed back into the input channel. January 2, 2012
36 of 218
Ramprasad Potluri
Lecture 8
Virtues of control systems: Improvement in performance and overcoming deadzone In this lecture, we discuss how a control system can be employed to achieve the goal of the output of a plant following a reference signal. The task of following a reference signal can be broken into two components. One, when the reference signal changes to a new constant value yref (∞), the output needs to also eventually reach a new constant value y(∞) which is as close as desired to yref (∞). Second, the while the output of the plant needs
8.1
Improvement in speed of response
We saw in the last lecture that, when we wanted a motor’s speed to quickly reach a unit speed, a closed-loop configuration was more useful than an openloop configuration. We saw this as follows. We make the following assumptions: • The motor is powered from a power amplifier, and a tachogenerator is mounted on the motor’s shaft to measure the speed of the motor. The TF from the power amplifier’s input to the tachogenerator’s output is τs1+1 . • The signals rOL and yOL (and rOL and yOL ) are dimensionally consistent and of comparable magnitudes. • rOL (t) = rCL (t) = 1(t). We see that yOL and yCL behave as shown in the following figure: ) (t) rOL (t)y=1 (utOL OL τs + 1
yOL (t) = (1 − e−t/τ )1(t)
yrueCL (`(t1)(pt()tt)) CL CL + − K τs6 +1
yCL (t) = K 1+ K
37
τ 1 − e−t/ 1+K 1(t)
EE250 (Control Systems Analysis) IITK
Lecture Notes
Apparently, the CL system provides faster settling time at the small cost of a steady-state error ess , limt→∞ e(t) equal to 1+1 K . But how is the CL system achieving this? Let us compare the control effort in the OL and CL cases. While, uOL (t) = 1(t), uCL (t) can be obtained as follows: UCL (s) K = RCL (s) 1 + τsK+1 Given that rCL (t) = 1(t), we have uCL (t) = L−1 {
K (τs + 1) 1 } τs + 1 + K s
This gives uCL (t) =
τ K 1 + Ke−t/ K+1 1(t) 1+K
We see that, at t = 0+ , uCL = K, and at t = ∞, uCL = the following figure,
K 1+ K .
As we can see from
K
u CL C
y
K 1+ K
L
1 uOL yOL
t in order to provide the improved settling time, the CL system is applying a greater control effort: uCL (t) > uOL (t),
for most t.
Thus, if we want a motor to settle as quickly as possible to 1000 rpm, we drive it, using uCL (t), for example, towards 10000 rpm, and as it is tending towards 10000 rpm, we ask it to go to 8000 rpm, 7000, 6000,...,1000. The idea is that no matter what inertia the motor has, we can always supply sufficient power to help it overcome the inertia quickly. But, in practice, it is not possible to provide enough power to overcome any inertia, mostly because the power source or the motor have their own limitations. The power source may be able to apply only so much voltage, or only so much current may be permitted to pass through the armature of the motor. That is, in practice, there are always constraints on the inputs to the system, and the input cannot exceed certain values. The feedback control system needs to respect these “input saturation constraints’. January 2, 2012
38 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
8.2
Lecture Notes
Command feedforward
There is also what is known as command feedforward: For improving performance through improving the shape of response. Y (s) (C + C2 ) P = 1 Yd (s) 1 + PC1
yC y` p pd`P + + − 21? 6
In the above scheme, we can use command feedforward to alter the zeros of (C +C ) P the TF from yd (t) to y(t). For this, we can modify the numerator of 11+ PC2 to 1 obtain certain desirable properties of the response y to input yd (for example, rise time). In the chapter on performance analysis of control systems, we will see what effect the zeros of a TF have on its response. Let us compare through an example the command feed-forward (CFF) structure with the negative feedback (NFB) structure. We saw that the NFB structure helps achieve a good speed of response. In the following figure, yruCL e(−`(1tp()t()tt)) CL CL + K τs 6 +1 yCL (t) =
τ K 1 − e−t/ 1+K 1(t) 1+K
Thus, while the plant by itself has a speed of response governed by its time constant τ, the closed-loop structure has a time constant τ/(1 + K ), and thus is faster by 1 + K times. However, the steady-state value of yCL (t) is yss = K/(1 + K ), and is not equal to 1. On the other hand, the CFF scheme below
has the TF
K + C2 1 Y = Yd τ s + 1+τ K
and, for yd (t) = 1(t), this gives y(t) =
τ K + C2 1 − e−t/ 1+K 1(t) K+1
This scheme, while allowing us to have the same speed of response as in the case of the NFB configuration, allows us, for C2 = 1, to have yss = 1. That is, this scheme also improves the steady-state error. This is not the only advantage of the CFF structure. Note that, real-world plants do not respond if their stimuli are not strong enough. That is, these plants are said to have “dead zone” nonlinearities. For example, a DC motor’s shaft will not turn unless the voltage applied to its armature is larger than a certain minimum value. When the value is below this minimum, the power supplied to the motor is not enough to help the shaft overcome the friction in the bearings, and to overcome the losses in the electric and magnetic circuit of the motor. The adjacent figure illustrates two nonlinearities that we find in January 2, 2012
39 of 218
Ramprasad Potluri
Lecture Notes
Output
EE250 (Control Systems Analysis) IITK
Li n
ea
rr eg
io
n
Saturation
Dead zone 0
Input
a power amplifier-DC motor combination − → Power-amp − → DC motor . The power amplifier (such as an H-bridge circuit) and the DC motor may have both a dead zone nonlinearity and a saturation nonlinearity. A saturation nonlinearity is the situation where raising the magnitude of the stimulus beyond a point does not help anymore. In the given NFB scheme, the magnitude of the voltage applied to the power-amplifier depends on the error. If the error is small, then this voltage will also be small, and may lie in the dead zone of the amplifier. In the CFF scheme, on the other hand, a certain minimum voltage is applied to the power amplifier thanks to feeding forward of the command input. This feature can be utilized to stay out of the dead zone of the amplifier.
January 2, 2012
40 of 218
Ramprasad Potluri
Part II
Bode Plots
41
Lecture 9
Frequency response 9.1
Definition of Frequency Response
Frequency response (FR) is the steady-state (SS) response of a system to a sinusoidal input. SS is when the transient response of the system has died out. So, by definition, in the development of the concept of FR, we will talk about systems where transient responses die out.
9.2
Problem statement
We will show that the SS value of y(t) — ySS (t) — in the following figure r (t)= A sin ωt
y(t)
R(s)
Y (s)
−−−−−−−→ G (s) −−→ is given by yss (t) = | G ( jω )| A sin (ωt + ∠G ( jω )) G ( jω ) is called “sinusoidal TF” of the system G ( jω ) = G (s)|s= jω
9.3
(9.1)
Solution
We study systems modeled by LTI differential equations of the form of Equation (3.1) including a dead time component. As these differential equations represent practical systems, they have real coefficients. A corollary of the fundamental theorem of algebra is that every polynomial with real coefficients can be expressed as the product of polynomials with real coefficients of degrees 1 or 2. The above two paragraphs mean that the TFs that we study in this course are of the form G (s) = a
n
n
d
d
2 2 2 ∏i=1 0 (s + zi ) ∏k= 0 ( s + 2ζ k ωnk s + ωnk ) 2 2 2 ∏i=1 0 (s + pi ) ∏k= 0 ( s + 2ξ k ωdk s + ωdk )
43
e−td s
(9.2)
EE250 (Control Systems Analysis) IITK
Lecture Notes
where a, zi , pi , ωnk , ωdk , td , ζ k , ξ k are all real. In the TF of Equation (9.2), the part that does not include the exponential term e−td s is said to be rational as it is the ratio of two terms. This part is said to be proper if the order of the numerator is not more than the order of the denominator, that is, if n1 + 2n2 ≤ d1 + 2d2 and strictly proper if the order of the numerator is less than the order of the denominator, i.e., if n1 + 2n2 < d1 + 2d2 As we are interested in the SS response of G (s), G (s) will need to be stable. Else, the SS value will not exist. We will discuss stability in the required detail in a later lecture. The term e−td s represents “dead time” or “transportation lag”. We evaluate yss (t) for r (t) = A sin ωt under the assumption that G (s) is stable and that e−td s is absent.
Steps 1. R(s) =
Aω . s2 + ω 2
2. Y (s) = G (s) R(s) can be written as follows: Y (s) = G (s) ·
Aω (s + jω )(s − jω )
(9.3)
3. Split Y (s) into partial fractions: Y (s) = Ψ( G (s)) +
γ2 γ1 + s + jω s − jω
(9.4)
where Ψ( G (s)) represents all the partial fractions that G (s) generates. 4. Evaluate γ1 & γ2 using (9.3) & (9.4). G (s)
Aω γ1 γ2 = Ψ( G (s)) + + (s + jω )(s − jω ) s + jω s − jω
γ1 and γ2 can be evaluated as follows: Aω γ1 = (s + jω ) G (s) (s + jω )(s − jω ) s=− jω Aω γ2 = (s − jω ) G (s) (s + jω )(s − jω ) s=+ jω
Therefore, γ1 = January 2, 2012
G (− jω ) Aω , −2jω
γ2 =
44 of 218
G (+ jω ) Aω . +2jω Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
5. Are G (− jω ) & G (+ jω ) related? In particular is G (− jω ) = G (+ jω )∗ ? We can see that this is indeed so as all the constant in G (s) (i.e., a, zi , pi , ωnk , ωdk , td , ζ k , ξ k ) are all real numbers. If they were not real, then in general G (− jω ) would not be the complex conjugate of G ( jω ). Since we are dealing with real systems, these constants are real in our case. 6. Now, let us take the inverse LT of (9.4): y(t) = L−1 {Ψ( G (s))} + γ1 e− jωt + γ2 e+ jωt
(9.5)
Here L−1 {Ψ( G (s))} contains all the transient terms. 7. As we have assumed that G (s) is stable, then L−1 {Ψ( G (s))} dies out eventually. Then, yss (t) = γ1 e− jωt + γ2 e+ jωt
=
G (− jω ) Aω − jωt G ( jω ) Aω jωt e + e −2jω −2jω
| G ( jω )|e− jφ A − jωt | G ( jω )|e jφ A jωt e + e −2j 2j | G ( jω )| j(ωt+φ) = A e − e− j(ωt+φ) 2j ⇒ yss (t) = | G ( jω )| A sin(ωt + ∠G ( jω )) =
Problem 9 Evaluate yss (t) assuming that e−td s is also present.
9.4
Conclusion
We have proved that the S.S. response of G (s) to an input r (t) = A sin ωt is yss (t) = | G ( jω )| A sin(ωt + ∠G ( jω )) So, the term G ( jω ) = | G ( jω )|∠G ( jω ) = | G ( jω )|e j∠G( jω ) fully characterizes the S.S. response of the stable TF G (s) to a sine input. So, the output sine from G (s) is r (t) amplified | G ( jω )| times and shifted ahead by a phase of ∠G ( jω ).
9.5
Caution
While G ( jω ) can be written as in Equation (9.1), and though G (s) =
Y (s) , R(s)
Y (s) Y ( jω ) G ( jω ) 6= = R(s) s= jω R( jω ) January 2, 2012
45 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
For, while G ( jω ) can be evaluated as done above, Y ( jω )/R( jω ) works out to zero because R( jω ) = R(s)|s= jω =
9.6
Aω Aω = =∞ ( jω )2 + ω 2 −ω 2 + ω 2
Appendix
FR is developed in most books. My derivation is based on [10]. But, I have taken care to avoid Ogata’s errors. Note that, as per Section 9.1, FR exists only for stable systems. We develop the machinery of Bode plots on the basis of FR. However, Bode plots are drawn for unstable systems too. I have not seen that the standard textbooks explain this anomaly. Indeed, they do not even mention it.
January 2, 2012
46 of 218
Ramprasad Potluri
Lecture 10
Bode plots (Part 1) 10.1
Uses of frequency response in control systems
Frequency response (FR) G ( jω ) of a system may be sketched with | G ( jω )| against ω and ∠G ( jω ) against ω as shown in the Figure 10.1. One value of an FR is that it allows us to compare the speed of response of two systems. For example, from the Figure 10.2, we can say that system 2 with TF G2 (s) is faster than system 1 that has the TF G1 (s). This statement can be understood as follows. In order to amplify a signal, the system needs to work, needs to put in energy. To amplify a fast sinusoid (i.e., one that has a high frequency), the system needs to work fast, to put in energy fast. If the system can process a fast sinusoid, then it can process any other fast signal whose fastest sinusoidal component is this fast sinusoid. System 2 above is able to amplify with the same constant amplification factor as system 1 more frequencies than system 1 is able to. The speed of response of a system is measured in the time domain through such quantities as time constant, rise time, settling time. For example, in Lecture 7, we saw that the OL system had a slower response than the CL system. In the real world, all physical quantities evolve in the time domain, and we tend to think of speed of response of a system in the time domain. It may not seem very intuitive to think in the frequency domain. But, the frequency domain does have conveniences that we will exploit in this course. For example, the frequency-domain quantity “bandwidth” is roughly the inverse of the time-domain quantity “time-constant” of a system. Thus, in order to make a system faster, we need to increase its bandwidth. Here is how we can determine the frequency response of a system in practice. If it is an electrical circuit, input a sine voltage of a fixed amplitude from a function generator. Sweep the frequencies from 0 to some large value, say, 20 kHz or upto the frequency that the function generator will support. Take readings of the output sine’s amplitude. Also, read the output sine’s phase w.r.t. to the input by observing both on an oscilloscope. Overall, take about 10 - 15 readings of the amplitudes and the corresponding phases. Plot the amplitude versus frequency, and the phase versus frequency. This is the frequency response plot. Catch: the readings must be logarithmically spaced. Otherwise, the plot 47
EE250 (Control Systems Analysis) IITK
Lecture Notes
ω, [rad/s]
0
0
linear scale
G ( jω ) [rad] 6
log scale
ω, [rad/s]
G ( jω ) [rad] linear scale
linear scale
6
linear scale 20 log10 | G ( jω )| [dB]
linear scale
| G ( jω )|
ω, [rad/s]
0
log scale
linear scale
ω, [rad/s]
ω, [rad/s]
10 log scale
linear scale 20 log10 | G2 ( jω )| [dB]
linear scale 20 log10 | G1 ( jω )| [dB]
Figure 10.1: Sample sketches of | G ( jω )| and ∠G ( jω ) versus ω.
ω, [rad/s]
100 log scale
Figure 10.2: Magnitude versus frequency plots with different bandwidths. will not be useful or even clear. How to logarithmically space readings? This brings us to the topic of Bode plots (BPs). Bode plots are convenient tools to display the FR of a system.
10.2
Convention in constructing Bode plots
Consider the general form of the TF of Equation (9.2). The FR of the system whose TF G (s) is is given by the sinusoidal TF n
G ( jω ) = a
n2 2 2 ∏i=1 0 ( jω + zi ) ∏k= 0 ( jω ) + 2ζ k ωnk jω + ωnk d
d
2 2 2 ∏i=1 0 ( jω + pi ) ∏k= 0 (( jω ) + 2ξ k ωdk jω + ωdk )
e−td jω
The magnitude of G ( jω ) is q
2 (ωnk 2 − ω 2 ) + (2ζ k ωnk ω )2 q | G ( jω )| = a 2 d p d2 (ωdk 2 − ω 2 ) + (2ξ k ωdk ω )2 ∏i=1 0 ω 2 + pi 2 ∏k= 0 n ∏i=1 0
January 2, 2012
p
2 ω 2 + zi 2 ∏nk= 0
48 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
The phase of G ( jω ) is
n1
G ( jω ) =
∑
n2
jω + zi +
i =0
∑
.
ωnk 2 − ω 2 + j2ζ k ωnk ω − td ω −
k =0 d1
∑
d2
jω + pi −
i =0
∑
.
ωdk 2 − ω 2 + j2ξ k ωdk ω
k =0
While | G ( jω )| is straightforward to evaluate, ∠G ( jω ) needs care. Here is the convention that control systems theory follows in evaluating ∠G ( jω ): 1. ∠ (ωk + jω ) can be in the interval [0, 2π ] rad. This restriction does not apply to the exponential term e− jωtd ; here the angle is −ωtd in radians. 2. Angles of complex numbers are measured in the counter clockwise sense. Example 10 Determine ∠G ( jω ) given that G (s) =
(s + z1 )(−s − z2 )e−td s (s + p1 )(s − p2 )(−s + p2 )
Note that all the coefficients in the TF are positive1 .
∠G ( jω ) = ∠ (z1 + jω ) + ∠ (−z2 − jω ) + (−td ω )− − ∠ ( p1 + jω ) − ∠ (− p2 + jω ) − ∠ ( p2 − jω ) Let us calculate this angle for some fixed value ω = ω0 . Angle
6
6
6
6
Value of angle π 4
(z1 + jω0 )
π+
(−z2 − jω0 ) 6
1 We
Location of complex number
π 6
π 3
( p1 + jω0 )
(− p2 + jω0 )
π−
π 4
( p3 − jω0 )
2π −
π 4
talk of positiveness or negativeness of only real numbers, and not of complex numbers.
January 2, 2012
49 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Then,
∠G ( jω ) =
π π π π π − ( t d ω0 ) − − π − − 2π − rad + π+ 4 6 3 4 4
The method demonstrated in this example to determine the angle of a complex number can prove convenient in determining the .phase of sinusoidal TFs. k ω0 j ωω0
jω ω0
+ 1 by noting that as ω goes from 0 to ∞, the complex number 1 + goes from 1 + j0 to 1 + j∞, and consequently, the angle of this complex number goes from 0 to 90◦ . Therefore, the plot of ∠G ( jω ) versus ω lies between 0◦ and −90◦ . For example, we can find the phase of G ( jω ) =
10.3
Types of Bode plots
We will use the following tools in our course: • Bode plot (BP) • Bode magnitude plot (BMP) • Bode phase plot (BPP) • Asymptotic Bode plot (ABP) • Asymptotic Bode magnitude plot (ABMP) • Asymptotic Bode phase plot (ABPP) The ABP is a straight line approximation to BP. The (A)BP has two components: (A)BMP and (A)BPP.
10.4
Appendix: Angles in MATLAB
In MATLAB, to determine the angle of a complex number x + jy, we have several functions. We discuss here the suitability of these functions for the purpose of drawing Bode plots. Suppose we wish to determine the angle of the complex number -10-j10. atan x = -10, y = -10, atan(x+i*y) will not work. 270 atan(abs(y/x))*180/pi will work. That is, we have to figure out the quadrant in which x + jy lies2 . In this case, x + jy lies in the third quadrant. atan2 Takes into account the quadrant. However, the solution, similar to MATHEMATICA’s arctan, is in the range [−π, +π ]. So, for example, atan2(-10,-10)*180/pi returns −135, whereas for the purposes of plotting BPPs, we would like to obtain the answer as +225. 2 On the other hand, in MATHEMATICA, ArcTan(x,y) will take into account the quadrant in which x + jy lies.
January 2, 2012
50 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
angle Takes into account the quadrant in which the complex number lies. For example, angle(-10-j*10)*180/pi gives −135. That is, angle returns the angle in the range [−π, +π ]. bode This function does not have the drawback of atan2 or angle. For example, the commands sys = tf([10,10,0,0],1); [mag,phase]=bode(sys,1) will return phase = 225.
10.5
Appendix: Angles in GNU Octave
GNU Octave’s atan, atan2, angle seem to be behaving the same way as their counterparts in MATLAB. The bode function in these two software behave differently from each other. While MATLAB seems to have implemented its bode function to calculate angles as explained in Section 10.2, GNU Octave’s bode seems to be following a different convention to calculate angles. In this connection, two observations about GNU Octave version 3.2.4 are 1. The commands sys = tf([10,10,0,0],1); [mag,phase]=bode(sys,1) returns an error message that the number of poles is less than the number of zeros. That is, GNU Octave’s bode seems to work only with proper transfer functions. 2. sys = tf(1,[10,10,0,0]); [mag,phase]=bode(sys,1) returns phase = 135. This command returns phase = -225 in MATLAB, which is in keeping with the convention of Section 10.2. My guess is that when x + jy is in the first and second quadrants, GNU Octave 3.2.4’s bode calculates its angle in keeping with the convention of Section 10.2, while if the complex number is in the third and the fourth quadrants, a different convention is used.
10.6
Appendix: Examples of departure of MATLAB from the convention
Problem 10 If you were to design a software that would draw the Bode plots as per the convention of Section 10.2, sketch how the phase plot in response to the following commands would look: sys=tf([1],conv([1 1],[1 1 1])); bode(sys) sys=tf([1],conv([1 1],[1 0 1])); bode(sys) Figures 10.3 and 10.4 show the result generated by the first line of commands in GNU Octave 3.2.4 and MATLAB respectively. Figures 10.5 and 10.6 show the result of the second line of commands generated by GNU Octave 3.2.4 and MATLAB respectively.
January 2, 2012
51 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
|[Y/U](jw)|, u=u1, y=y1 0
-20
Gain in dB
-40
-60
-80
-100
-120 10-1
100
101
102
101
102
Frequency in rad/sec
phase([Y/U](jw)), u=u1, y=y1 0
-50
Phase in deg
-100
-150
-200
-250
10-1
100 Frequency in rad/sec
Figure 10.3: Result of sys=tf([1],conv([1 1],[1 1 1])); bode(sys) from GNU Octave version 3.2.4. Bode Diagram
Magnitude (dB)
0
−50
−100
Phase (deg)
−150 0
−90
−180
−270 −2 10
−1
10
0
10 Frequency (rad/sec)
1
10
2
10
Figure 10.4: Result of sys=tf([1],conv([1 1],[1 1 1])); bode(sys) from MATLAB. January 2, 2012
52 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
|[Y/U](jw)|, u=u1, y=y1
50
Gain in dB
0
-50
-100
10-1
100
101
102
101
102
Frequency in rad/sec
phase([Y/U](jw)), u=u1, y=y1 140 120 100
Phase in deg
80 60 40 20 0 -20 -40 10-1
100 Frequency in rad/sec
Figure 10.5: Result of sys=tf([1],conv([1 1],[1 0 1])); bode(sys) from GNU Octave version 3.2.4. Bode Diagram 200 150 Magnitude (dB)
100 50 0 −50 −100 −150
Phase (deg)
−200 −180
−270
−360
−450 −2 10
−1
10
0
10 Frequency (rad/sec)
1
10
2
10
Figure 10.6: Result of sys=tf([1],conv([1 1],[1 0 1])); bode(sys) from MATLAB. January 2, 2012
53 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
54 of 218
Lecture Notes
Ramprasad Potluri
Lecture 11
Bode plots (Part 2) 11.1
First method for sketching an ABMP k s + ω0 .
As an example, let us sketch the BP and ABP of G (s) =
Step 1 First convert G (s) into what we shall call in this course the “Bode form”: G (s) =
k/ω0 s ω0 + 1
This form is convenient for sketching BPs and ABPs. Step 2 Write the sinusoidal TF corresponding to G (s). G ( jω ) =
k/ω0 jω0 ω0
+1
; | G ( jω )| = r
k/ω0 2 ; ∠G ( jω ) = − 1 + ωω0
1+j
ω ω0
Step 3 To construct the BP, populate Table 11.1. In the table, dB-gain is as follows: s 2 k ω +1 + −20 log10 dB-gain = 20 log10 | G ( jω )| = 20 log10 ω0 ω0 | {z } | {z } C1
C2
Construct the BP using the data populated in Table 11.1.
Table 11.1: Table to populate to build the Bode plot of G (s) = k (s + ω0 ) . ω, [rad/s]
ω0 100
ω0 10
ω0 2
dB-gain, [dB]
∠G ( jω ), [rad]
55
ω0
2ω0
10ω0
100ω0
Lecture Notes
linear scale 20 log10 |G(jω)| [dB]
EE250 (Control Systems Analysis) IITK
0 log10
³
ω, [rad/s], if log scale ´ ω ω0 , [decades], if linear scale
linear scale 6 G(jω) [rad] or degrees
log10 ω, [decades], if linear scale
log10
³
ω, [rad/s], if log scale ´ ω ω0 , [decades], if linear scale
log10 ω, [decades], if linear scale
Figure 11.1: Setting up the axes on the graph paper on which we sketch the Bode plot.
Step 4 To construct the ABMP, note that taking the logarithm of a product results in a sum, and that a sum is easier to sketch than a product. Thus, we can sketch the ABMP of G ( jω ) as the sum of the ABMP corresponding to C1 with that corresponding to C2. Let us consider C2. We note the following about C2: s ( 2 ω 0, ω ω0 1 (i.e., when ω / ω0 /10) −20 log10 +1 = ω0 −20 log10 ωω0 , ωω0 1 (i.e., when ω ' 10ω0 ) Step 5 Set up the axes as shown in Figure 11.1. We can use either a linear-linear graph paper or a semilog graph paper. A 2-cycle semilog graph paper is shown in Figure 11.2 and a 5-cycle paper is provided in a larger size at the end of this book. Step 6 Sketch the straight line approximations corresponding to C1 and to Step 4’s approximation of C2 as shown in Figure 11.3. Here is a brief explanation for how section C2-2 is constructed in Figure 11.3. C2-2 is the plot of −20 log10 ω/ω0 versus log10 ω/ω0 . Sketching this section is equivalent to sketching y = −20x versus x for x ≥ 0 as shown in January 2, 2012
56 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
20 log10 |G(jω)|, [dB]
Figure 11.2: A 2-cycle semilog grid. The ordinate is graduated linearly. The abscissa is graduated logarithmically meaning that within each cycle on the abscicca from the left, the marks corresponding to 1, 2, 3, . . . , 10 are respectively at log10 1 = 0, log10 2 = 0.30103, log10 3 = 0.47712, . . . , 1 units.
C1
20 C2-1 0
ω − 10
ω0 −10ω0
−20 −40
−20 dB/dec
C2-2
Figure 11.3: A sketch of the straight line approximations to C1 and C2.
January 2, 2012
57 of 218
Ramprasad Potluri
Lecture Notes
y = −20x
EE250 (Control Systems Analysis) IITK
0
1
3
2
x
4
−20 −40 −60
ABMP
C1 of k
s+
ω0
20 log10 | G ( jω )|, [dB]
Figure 11.4: Sketching section C2-2 is equivalent to sketching y = −20x versus x for x ≥ 0.
20 C2-1 0
ω − 10
ω0 −10ω0
−20 −40
−20 dB/dec
C2-2
Figure 11.5: The final asymptotic Bode magnitude plot of k (s + ω0 )
Figure 11.4. On a graph of 20 log10 | G ( jω )| versus log10 ω/ω0 , the 0 on the abscissa is at ω = ω0 . Therefore, beginning at ω = ω0 , the ABMP of −20 log10 ω/ω0 is a straight line at a slope of −20 dB/dec. Step 7 Add C1 and C2 point by point to obtain the ABMP of k (s + ω0 ) as shown in Figure 11.5. The maximum difference between the BMP and the ABMP is 3 dB and occurs at the corner frequency as shown in Figure 11.6. Step 8 Draw the BPP using the data filled into the table of Step 3. The BPP must be presented below the BP with the grids aligned vertically as shown in Figure 11.6. Step 9 Draw a straight line approximation to the BPP. There can be multiple choices for this straight line approximation. See Figure 11.6 for the most commonly used choice. January 2, 2012
58 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Bode Diagram 0 −5
Magnitude (dB)
−10
3 dB
−15
−20 dB/dec
−20 −25 −30 −35
Phase (deg)
−40 0
−45 deg/dec
o
5.7 −45
−90 −2 10
−1
10
0
10 Frequency (rad/sec)
1
2
10
10
Figure 11.6: Bode plot and asymptotic bode plot of G (s) = 1 (s + 1) . The straight line plot in the thick solid line is the asymptotic plot, while the smooth line plot in the thin solid line is the Bode plot. The maximum difference between these plots is at the corner frequencies. The maximum error between the magnitude plots is 3 dB, while that between the phase plots is 5.7◦ .
Examples of straight line approximations to BPPs are in [12] and [13].
11.2
Second method for sketching an ABMP
The procedure shown above constructs ABMPs as sum of elemental ABMPs. In practice, however, a much simpler method exists to quickly sketch ABMPs, and works as follows. Step 1 The shape of the ABMP is visualized first in terms of the slopes of its various sections, and the corner frequencies. This step fixes the horizontal location of the ABMP. Step 2 The vertical position of the ABMP can be fixed by working out the value of the dB-gain at any one frequency. Thus, to sketch the ABMP of k (s + ω0 ) , we note that the ABMP of this TF has the shape 0 dB/dec, −20 dB/dec. The corner frequency is at ω = ω0 . One January 2, 2012
59 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
point on this ABMP is 20 log10 k/ω0 at ω ω0 , say ω = ω0 /10. This point fixes the ABMP in the vertical direction. As another example, to sketch the ABMP of k s(s + ω0 ) , we note that the ABMP of this TF has the shape −20 dB/dec, −40 dB/dec. The corner frequency is at ω = ω0 . One point on this ABMP is 20 log10 10k/ω0 2 at ω ω0 , say ω = ω0 /10. This point fixes the ABMP in the vertical direction.
January 2, 2012
60 of 218
Ramprasad Potluri
Lecture 12
Bode plots (Part 3) 12.1
Comparison of the two methods for constructing ABMPs
Working out a few problems with the first method helps develop a comfort with the second method. While the first method is simplified by using the Bode form of the TF, the second method does not need this form. In practice, whenever there is any confusion in using the second method, we can employ the first method. While the elemental ABMPs are constructed through an analytical approximation and do not need the BMPs to be constructed first, elemental ABPPs are not as easy to construct analytically. Therefore, the elemental BPPs are constructed first, and only then their approximations are constructed as ABPPs. We memorize the elemental ABPPs and sum them up when a composite ABPP is needed. While ABMPs can be constructed easily using the second method, ABPPs are more easy to construct by summing the elemental ABPPs. Elemental BPs are those corresponding to the various elemental TFs in the composite TF of (9.2).
12.2
Selecting an ω to vertically fix an ABMP in the second method
To fix the vertical location of the ABMP, working out the dB-gain in the lowfrequency section of the ABMP is less error-prone than working out the dBgain near one of the corner frequencies. For example, .p to fix the ABMP of k (s + ω0 ) , we use the expression for magnitude k ω 2 + ω0 2 . If we select ω = ω0 as the frequency at which to .√ work out the magnitude, then we obtain the point ω0 , 20 log10 k 2 ω0 , which lies on the BMP but not on the ABMP. Indeed, the ABMP at this ω = ω0 lies 3 dB above this point. We could instead use the point 61
EE250 (Control Systems Analysis) IITK
Lecture Notes
.√ 1.01ω0 , which is only 0.04 dB below the point on ω0 /10 , 20 log10 k the ABMP ω0 /10, 20 log10 (k /ω0 ) . Indeed, to fix the ABMP, it is convenient to use not just points from the lowfrequency section, but any points that are at least one decade away from their closest corner frequency.
12.3
Examples of construction of ABPs
Consider the TF k (s (s + ω0 )) . Consider the TF k s2 (s + ω0 ) . Consider the TF k ((s + 10)(s + 1000)) . Consider the TF k (s+ ω0 )n . Consider the TF ke−s (s(s + 10)(s + 1000)) . Most of the time of the lecture was spent on this section. I will fill this section soon.
12.4
BP and ABP of the second order TF
We consider the second order TF ωn 2 s2 + 2ζωn s + ωn 2 for 0 < ζ < 1. The ωn 2 in the numerator helps normalize the magnitude of this TF at low frequencies to 1. ζ is known as damping ratio. The magnitude of the sinusoidal version of this TF is
| G ( jω )| = q
ωn 2
( ωn
2 − ω 2 )2
+ (2ζωn ω )
2
= r
1−
1 2
+
ω2 1− ωn 2
ω2 ωn 2
2ζ ωωn
2
(12.1)
and the phase is
∠G ( jω ) = −
.
2
ωn − ω
2
+ j2ζωn ω = −
+ j2ζ
ω ωn
We see that, as r = ωωn goes from 0 to ∞, the complex number (1 − r2 ) + j2ζr moves from the positive real axis to the first quadrant, to the positive imaginary axis, to the second quadrant, to the negative real axis. So, the phase of this number goes from 0 to 180◦ . So, the phase of our 2nd order TF goes from 0◦ to −180◦ . The phase goes from 0◦ to −180◦ as ω goes from 0 to ∞. To obtain a straight line approximation to the BMP, we note that ω 0, ω0 1 (i.e., when ω / ω0 /10) −20 log10 | G ( jω )| = 2 ω −20 log , ωω0 1 (i.e., when ω ' 10ω0 ) 10 ω0 This means that the ABMP has the slopes 0 dB/dec, −40 dB/dec, and that the corner frequency is at ω = ωn . While the ABMP of the second order TF is convenient, the difference between the BMP and the ABMP can be significant for values of ζ much smaller than 1 as Figure 12.1 shows. Problem 11 January 2, 2012
62 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
20
ζ=0 10
ζ=0.15713
dB gain of 2nd order TF
ζ=0.31427 ζ=0.4714 ζ=0.62854 ζ=0.78567 ζ=0.94281
0
ζ=1.0999 ζ=1.2571
−10
ζ=1.4142
−20
−30
−40 −1 10
0
1
10 r = ω/ωn
10
Figure 12.1: Bode magnitude plots of the second order transfer func tion ωn 2 s2 + 2ζωn s + ωn 2 . This figure is generated using the m-file order2tf.m.
Write a MATLAB code to plot the BMP of the second order TF for various values of ζ, and to mark the locus of the maxima of these BMPs.
12.5
Appendix: MATLAB codes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % order2tf.m: m-file to evaluate the value of the dB-gain of % the TF % wn^2 % -------------------% s^2 + 2*zeta*wn + wn^2 % % at various frequencies and for various values of zeta. % % This task can be implemented using a FOR loop. % There may be other ways too. This code works nicely in % GNU Octave version 3.2.4 and in MATLAB. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% January 2, 2012
63 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
clear all % Clears workspace of the values of all variables % that have been evaluated. close all % Closes all figure windows. % Specify values of the damping ratio for which Bode plots % will be drawn: nzeta = 10; zeta = linspace(sqrt(2),0,nzeta); % Specify logarithmically-spaced frequencies at which dB-gain % & phase will be evaluated: nr = 150; r = logspace(-3,+1,nr)’; % r = w/wn for i = 1:nzeta for j = 1:nr dBgain(j,i) = -10*log10((1-r(j)^2)^2 + (2*zeta(i)*r(j))^2); % What this is meant to generate is a matix of dB gains % where each column corresponds to a given value of zeta and % each row corresponds to a given value of r. end end % Plot the Bode diagram for the first value of zeta: a = [’\zeta=’,num2str(zeta(1))] semilogx(r,dBgain(:,1)), grid, gtext(a); axis([0.1,10,-40,+20]); xlabel(’r = \omega/\omega_n’), ylabel(’dB gain of 2nd order TF’); % Once the first diagram has been plotted, plot all % subsequent diagrams on the same figure: hold for i = 2:nzeta a = [’\zeta=’,num2str(zeta(i))] semilogx(r,dBgain(:,i)); gtext(a); end % Next, show on the same plot, the locus of the peaks of each % Bode magnitude plot: zeta1 = linspace(sqrt(2),0,5*nzeta); r = sqrt(1-2*zeta1.^2); dBmax = -20*log10(2.*zeta1.*sqrt(1-zeta1.^2)); semilogx(r,dBmax,’-.r’) print -depsc order2tf.eps % -dbitmap, -depsc, -deps are available that are relevant % to us.
January 2, 2012
64 of 218
Ramprasad Potluri
Lecture 13
Bode plots (Part 4) 13.1
Notes on the ABMP of second order TF
The following points are of interest about | G ( jω )| of (12.1). 1. The maximum of | G ( jω )| occurs at ω = ωd := ωn 2. dB-gain at ω = ωd is 20 log10
√1 2ζ
1− ζ 2
p
1 − 2ζ 2 .
.
√ √ 3. At ζ = 1/ 2, dB-gain equals zero at ω = ωd . For ζ < 1/ 2, dB-gain at ωd is greater than zero. That is, a hump begins to appear in the BMP for √ ζ < 1/ 2. √ √ 4. At ζ = √ 1/ 2, | G ( jωn )| = 1/(2ζ ) = 1/ 2.√Therefore, dB-gain at ω = ωn for ζ = 1/ 2 is −3 dB. Therefore, for ζ = 1/ 2, ω = ωn is the 3-dB frequency.
13.2
Complete BP of second order TF
The complete BP of the second order √ TF is shown in Figure 13.1 for values of the damping ratio ζ going from 2 to nearly 0 and is built using the m-file order2bp.m. Problem 12 Write a GNU Octave or MATLAB code to plot the BP of the second order TF for various values of ζ.
Problem 13 Run the m-file order2bp.m in MATLAB and modify it suitably if necessary to plot a figure that reproduces the salient information provided by Figure 13.1.
65
EE250 (Control Systems Analysis) IITK
Lecture Notes
40
20
dB-gain
0
-20
-40
-60
-80
-100 10-2
10-1
100
101
102
101
102
ω, rad/s
0 -20 -40
degrees
-60 -80 -100 -120 -140 -160 -180 10-2
10-1
100 ω, rad/s
Figure 13.1: √The BP of the second order TF for values of the damping ratio ζ going from 2 to close to 0. The smaller the value of ζ, the more pronounced the hump in the BMP, and the more abrupt the transition from 0◦ to −180◦ .
January 2, 2012
66 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
13.3
Lecture Notes
Filter-related terminology used in BPs
Review of ideal low pass, band pass, high pass filters. Review of concepts of bandwidth and cut off frequency for ideal filters. That we are talking about frequency response of control systems means that we can treat control systems as filters. In practical filters, as it is not possible to have sharp demarcations of passband and stopband, the passband is defined as “the frequency range where the filter’s amplitude response is greater than some arbitrarily √ chosen minimum value. Quite often, this minimum value is chosen as 1/ 2 of the maximum amplitude response” [14, pages 178-179]. Certain frequency responses (e.g., those whose BMP has a non-zero first slope, or those whose BMP has a non-zero last slope) are not even non-ideal versions of the ideal filters. So, for these, we cannot define bandwidth and cut-off frequency. 2
ωn For example, for the TF s2 +2ζω 2 , the cutoff frequency is where the BMP n s + ωn √ falls to 1/ 2 of its low-frequency value (of 1).
13.4
Minimum phase TF
Example 11 Construct the Bode phase plots of the following TFs as sums of elemental sτ1 1−sτ1 1+sτ1 1−sτ1 phase plots: G1 (s) = 11+ +sτ2 , G2 ( s ) = 1+sτ2 , G3 ( s ) = 1−sτ2 , G4 (s ) = 1−sτ2 . Note that the BPP of G1 (s) is in the smallest interval of phases. In this sense, G1 (s) is a minimum-phase TF.
13.5
Appendix: MATLAB codes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % order2bp: m-file to plot the complete Bode plot of the TF % % wn^2 % -------------------% s^2 + 2*zeta*wn + wn^2 % % for various values of zeta. % % This code works nicely in GNU Octave version 3.2.4, and % in MATLAB. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all % Clears workspace of the values of all variables % that have been evaluated. close all % Closes all figure windows. January 2, 2012
67 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
% Specify values of the damping ratio for which Bode plots % will be drawn: nzeta = 10; zeta = linspace(sqrt(2),0.001,nzeta) wn = 1; w = logspace(log10(wn/100),log10(100*wn),500); for i = 1:nzeta [mag(i,:),ph(i,:),w] = bode(tf(wn^2,[1,2*zeta(:,i)*wn,wn^2]),w); end subplot(2,1,1), semilogx(w,20*log10(mag(:,:)));grid; ylabel(’dB-gain’); xlabel(’\omega, rad/s’); subplot(2,1,2), semilogx(w,ph(:,:));grid, axis([1/100,100,-180,0]); ylabel(’degrees’); xlabel(’\omega, rad/s’); print -depsc order2bp.eps
January 2, 2012
68 of 218
Ramprasad Potluri
Part III
Stability Theory
69
Lecture 14
Stability theory (Part 1) 14.1
Minimum phase TFs (concluded)
We define minimum-phase TFs as follows. Definition 1 (Minimum-phase transfer function) Of all stable TFs, and those that have the right-most poles only at the origin, that have the same BMP, the one that has the phase plot spread over the smallest interval of phases is said to be minimumphase TF. By Definition 1, of the TFs in Example 11 G1 (s) and − G1 (s) are minimumphase TFs, while G2 (s) and G1 (s)e−td s are non-minimum phase TFs. Some authors include unstable TFs too into the class of non-minimum phase TFs. Definition 1 is in agreement with [15, pages 332-333], and allows us to include the TFs that are not stable, for example, (s + 1) s and (s + 1) s2 , into the class of minimum phase TFs.
14.2
Why practical rational TFs are proper
If they were not, then the magnitude of their frequency response would go to infinity as ω → ∞. Real systems cannot provide such arbitrarily high gains; their BMPs level off or drop after a certain frequency, meaning their TFs are proper.
14.3
BIBO stability of LTI systems**
Definition 2 A signal r (t) is bounded if it is bounded from above and below by a non-negative finite number M, meaning − M ≤ r (t) ≤ M, 0 < M < ∞, for all t ∈ [0, ∞). As signals could also be of bounded energy and bounded power, whereas Definition 2 implies only bounds on the magnitude of the signal, to avoid confusion, we use the term bounded in magnitude (BIM) for a bounded signal. 71
EE250 (Control Systems Analysis) IITK
Lecture Notes
Definition 3 An LTI system is said to be bounded-input bounded-output (BIBO) stable if and only if every bounded-in-magnitude (BIM) input r (t) results in a BIM output y(t). Note that, for a time-domain scalar function r (t), the magnitude is simply the absolute value of r (t), denoted |r (t)| and equal to |r (t)| = sgn(r (t))r (t). Here, sgn(·) is the signum function that is defined in Equation (14.2). Here, we develop a necessary and sufficient condition for the BIBO stability of a causal LTI system that will helps us evaluate the BIBO stability of a given TF G (s). From a basic course on signals and systems, we know that the response y(t) of a system to an arbitrary input r (t) (for t ≥ 0) can be expressed in terms of its impulse response g(t) as y(t) =
Z ∞ 0
r (τ ) g(t − τ )dτ
This is summarized in the following figure: causal
δ(t)
−−−→
LTI
g(t)
−−−→
r (t)=
R∞ 0
r (τ )δ(t−τ )dτ
causal
−−−−−−−−−−−−−−→
=⇒
y(t)=
R∞
r (τ ) g(t−τ )dτ
0 −−−−−−−−−−−−−− →
LTI plant
plant
As δ(t − τ ) 6= 0 only at t − τ = 0, we see that the impulse response g(t − τ ) 1 is non-zero R ∞ only beginning t − τ ≥ 0 , that is only for τ ≤ t. Therefore, the integral 0 r (τ ) g(t − τ )dτ does not need to be evaluated for all τ ≥ 0, but only for 0 ≤ τ ≤ t. Therefore, Z ∞ 0
r (τ ) g(t − τ )dτ =
Z t 0
r (τ ) g(t − τ )dτ
We shall use this result in order to develop the required condition. Theorem 1 In a causal LTI system, forRa BIM input r (t), the corresponding response ∞ y(t) of the system is BIM if and only if 0 | g(τ )| dτ < ∞. Proof: Given that |r (t)| ≤ M < ∞, we have Z t Z t Z t |y(t)| = r (τ ) g(t − τ )dτ ≤ |r (τ ) g(t − τ )| dτ = |r (τ )| | g(t − τ )| dτ 0
0
Therefore,
|y(t)| ≤ M
Z t 0
0
| g(t − τ )| dτ = M
Z t 0
| g(θ )| dθ
Thus, for y(t) to be bounded for a bounded r (t), for all t ∈ [0, ∞), it is Rt sufficient that 0 | g(τ )| dτ be bounded for all t ∈ [0, ∞). That is, For all t ∈ [0, ∞), for a BIM r (t), y(t) is BIM if 1 Here,
R∞ 0
| g(τ )| dτ < ∞.
we are using the assumption of causality of the system. Real-world systems are causal.
January 2, 2012
72 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
R ∞ Next, to show that for y(t) to be BIM given r (t) is BIM it is necessary that only need to show an example where y(t) is not BIM for 0 | g ( τ )| dτ < ∞, we R∞ a BIM r (t) because 0 | g(τ )| dτ 6< ∞. Consider, for example, the bounded input r (τ ) = sgn ( g(t − τ )) where the signum function sgn is defined as follows: +1, sgn ( x ) = 0, −1,
x>0 x=0 x<0
(14.1)
Then, y(t) =
Z t 0
r (τ ) g(t − τ )dτ =
Z t 0
| g(t − τ )| dτ =
Z t 0
| g(θ )| dθ
Rt Thus, in this case, we see that if 0 | g(τ )| dτ 6< ∞, then y(t) 6< ∞. Thus, we have a situation where, even though the input to the system is BIM, the output is not. Thus, for y(t) to be BIM for a BIM r (t), for all t ∈ [0, ∞), it is necessary Rt that 0 | g(τ )| dτ be bounded for all t ∈ [0, ∞). That is, For all t ∈ [0, ∞), for a BIM r (t), y(t) is BIM only if
R∞ 0
| g(τ )| dτ < ∞. 2
14.4
Discussion of the theorem and its proof**
The proof of Theorem 1 is not entirely convincing. In particular, the proof of the statement For all t ∈ [0, ∞), for a BIM r (t), y(t) is BIM only if
R∞ 0
| g(τ )| dτ < ∞.
is not convincing. Put another way, this statement says For all t ∈ [0, ∞), for a BIM r (t), if
R∞ 0
| g(τ )| dτ 6< ∞, then y(t) is not BIM.
We may claim that we have notR really proved this statement. We may claim ∞ that we have only proved that “if 0 | g(τ )| dτ 6< ∞, then y(t) is not BIM” only R∞ for this particular r (t). For a different r (t), maybe even if 0 | g(τ )| dτ 6< ∞, y(t) is BIM. A careful survey of some literature shows that Theorem 1 is stronger than what the literature claims. Theorem 2 summarizes what the literature says. Theorem 2 A causal LTI system is BIBO stable if and only if that is, its impulse response is absolutely integrable. January 2, 2012
73 of 218
R∞ 0
| g(τ )| dτ < ∞,
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Proof: Given that |r (t)| ≤ M < ∞, we have Z t Z t Z t |y(t)| = r (τ ) g(t − τ )dτ ≤ |r (τ ) g(t − τ )| dτ = |r (τ )| | g(t − τ )| dτ 0
0
Therefore,
|y(t)| ≤ M
Z t 0
0
| g(t − τ )| dτ = M
Z t
| g(θ )| dθ
0
Thus, for y(t) to be bounded for a bounded r (t), for all t ∈ [0, ∞), it is Rt sufficient that 0 | g(τ )| dτ be bounded for all t ∈ [0, ∞). That is, For all t ∈ [0, ∞), for a BIM r (t), y(t) is BIM if
R∞ 0
| g(τ )| dτ < ∞.
R∞ Next, to show that if 0 | g(τ )| dτ 6< ∞, then the sytem is not BIBO stable, we need to show that exist some BIM signals for which the output is not BIM. Consider, for example, the bounded input r (τ ) = sgn ( g(t − τ )) where the signum function sgn is defined as follows: +1, x > 0 sgn ( x ) = 0, x=0 −1, x < 0
(14.2)
Then, y(t) =
Z t 0
r (τ ) g(t − τ )dτ =
Z t 0
| g(t − τ )| dτ =
Z t 0
| g(θ )| dθ
Rt Thus, in this case, we see that if 0 | g(τ )| dτ 6< ∞, then y(t) 6< ∞. Thus, we have a situation where, even though the input to the system is BIM, the output Rt is not. Thus, we have shown that when 0 | g(τ )| dτ 6< ∞, there may be BIM inputs for which the output is not BIM. 2
14.5
Appendix: Some literature that discusses BIBO stability
Before I realized the fine difference between theorems 1 and 2, I investigated [16, page 217], [17, pages 175–176], [18, pages 113 – 114], [14, page 77], and [19, page 38]. [17, pages 175–176] guesses that [19] may have originated Theorem 2. Problem 14 By the way, [19, page 54] says that Z ∞ 0−
dτg1 (τ ) g2 (t − τ ) =
Z ∞ 0−
dτg1 (t − τ ) g2 (τ )
Prove or disprove this statement.
January 2, 2012
74 of 218
Ramprasad Potluri
Lecture 15
Stability theory (Part 2) 15.1
A question of interest
Question 1 Given that in our course we are studying systems modeled by the TF (9.2) R ∞ what is the effect of the poles and zeros of G (s) on whether or not the condition 0 | g ( t )| dτ < ∞ is satisfied for |r ( t )| ≤ M < ∞? That is, how do the poles and zeros of G (s) affect its BIBO stability? We may wonder why Question 1 does not ask about the effect of e−td s on the BIBO stability of G (s). The following example answers our question. Example 12 Problem 6 of Section 3.2 stated that y(t) being a delayed-in-time version of x (t) can be summarized by the equation y ( t ) = x ( t − t d )1( t − t d ), and taking LT on both sides of this equation gives Y (s) = X (s)e−td s . So, e−td s represents a time-delay. What is the effect of e−td s on the BIBO stability of G (s)? Note that we can redraw the block diagram r (t)
y(t)
R(s)
Y (s)
−−−→ G (s) −−−→ as follows: n
n
r (t)
2 2 2 ∏i=1 0 (s + zi ) ∏k= 0 ( s + 2ζ k ωnk s + ωnk )
R(s)
d ∏i=1 0 (s +
−−−→
2 2 pi ) ∏dk= 0 (s
+ 2ξ k ωdk s + ωdk 2 )
z(t)
y(t)
Z (s)
Y (s)
−−−→ e−td s −−−→
The time delay element only delays z(t). It does not change the nature of z(t) in any other manner. Thus, if z(t) is BIM, then y(t) is also BIM, and if z(t) is not BIM, y(t) is not BIM. That is, e−td s has no effect on the BIBO stability of G (s).
Problem 15 Can we not express the fact that y(t) is a delayed-in-time version of x (t) as either y(t) = x (t − td ) or y(t) = x (t)1(t − td )? 75
EE250 (Control Systems Analysis) IITK
Lecture Notes
Table 15.1: Case of unrepeated poles. Impulse response of each term of G (s) depending on the location of the poles of this term in the s-plane. A1 s + B1 (s + a1 )2 + b12 g1 (t) = Ae− a1 t cos(b1 t + φ1 ) G1 (s) =
Poles in left half (LH) s-plane Poles on jω-axis
R∞ 0
0
R∞
t
0 g1 ( t )
| g1 (τ )|dτ < ∞
R∞
t
0
0
t
g2 ( t )
d1 = 0
0
| g1 (τ )|dτ 6< ∞
g2 ( t )
d1 > 0
0
| g1 (τ )|dτ < ∞
a1 = 0
R∞
g2 ( t ) = c 1 e − d 1 t
g1 ( t )
a1 > 0
c1 s + d1
G2 (s) =
| g1 (τ )|dτ 6< ∞
0
t
g1 ( t )
Poles in right half (RH) s-plane
a1 < 0
R∞ 0
g2 ( t )
d1 < 0 0
| g1 (τ )|dτ 6< ∞
t
R∞ 0
| g1 (τ )|dτ 6< ∞
0
Example 12 showed why we may safely ignore any time delay elements in G (s) in studying its BIBO stability. To study Question 1, we considered 2 cases: case of G (s) containing unrepeated poles, and case of G (s) containing repeated poles. In each of these two cases, G (s) can be written as follows (ignoring the td part): Case of unrepeated poles: We assume that the number of distinct real poles is q and that the number of pairs of distinct complex conjugate poles is r. G (s) =
A1 s + B1
( s + a1
)2
+ b1
2
+···+
Ar s + Br
( s + ar
)2
+ br
2
+
cq c1 +···+ s + d1 s + dq
Table 15.1 shows the form of the impulse response of each term of G (s) depending on the location of the poles of this term in the s-plane: Case of repeated poles: For convenience, we shall assume that we have one real pole and one pair of complex conjugate poles, and that the real pole is repeated k times, while the second is complex conjugate pole is repeated r times. Then, G (s) =
A1k A B1r B11 + · · · + 11 + r +···+ k 2 2 s + α [(s + β) + b ] ( s + β )2 + b2 (s + α) | {z } | {z } H1 (s)
January 2, 2012
H2 (s)
76 of 218
Ramprasad Potluri
t
EE250 (Control Systems Analysis) IITK
Lecture Notes
Table 15.2: Case of repeated poles. Impulse response of each term of H1 (s) and H2 (s) of G (s) depending on the location of the poles of this term in the s-plane. Response due to real roots of multiplicity k ≥ 2
β>0
α>0 Poles in LH splane
D t i −1
k
h1 ( t ) =
i e−αt ∑ (i − 1) !
i =1
Exponentially decaying response
R∞ 0
|h1 (τ )|dτ < ∞
D t i −1
i ∑ (i − 1) !
i =1
Growing response
|h1 (τ )|dτ 6< ∞
R∞ 0
k
D t i −1
i e−αt ∑ (i − 1) !
i =1
Exponentially growing response 0
∑
|h2 (τ )|dτ 6< ∞
β<0
h1 ( t ) =
R∞
Ei ti−1 cos(bt + φi ) − βt e ( i − 1) ! i =1 r
h2 ( t ) =
Growing response
α<0 Poles in RH splane
|h2 (τ )|dτ < ∞
β=0 k
0
∑
Exponentially decaying response 0
h1 ( t ) =
R∞
Ei ti−1 cos(bt + φi ) − βt e ( i − 1) ! i =1 r
h2 ( t ) =
R∞
α=0 Poles on jω-axis
Response due to complex conjugate root pairs of multiplicity r ≥ 2
|h1 (τ )|dτ 6< ∞
Ei ti−1 cos(bt + φi ) − βt e ( i − 1) ! i =1 r
h2 ( t ) =
∑
Exponentially growing response
R∞ 0
|h2 (τ )|dτ 6< ∞
Table 15.2 shows the form of the impulse response of H1 (s) and H2 (s) of G (s) depending on the location of the poles of this term in the s-plane:
15.2
Conclusions
R∞ 1. When a TF has poles in the left half s-plane (LHP), 0 | g(τ )|dτ < ∞, and the TF is BIBO stable. 2. When a TF has poles in the right half s-plane (RHP) or has repeated poles R∞ on the jω-axis, then 0 | g(τ )|dτ 6< ∞, and the TF is unstable. R∞ 3. When a TF has unrepeated poles on the jω-axis, then 0 | g(τ )|dτ 6< ∞, and the TF is not BIBO stable. However, the TF’s response to signals excluding only the following is BIM: 3.1. In the case of an unrepeated complex conjugate pole pair, if r (t) is a sinusoid whose frequency coincides with that of the pole pair, then y(t) is unbounded. January 2, 2012
77 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
For example, if G (s) = 1/ s2 + ω0 2 and the input is u(t) = sin (ω0 t), then, using the information from Table 15.2, the output is E1 cos (ω0 t + φ1 ) + E2 t cos (ω0 t + φ2 ). 3.2. In the case of a pole at the origin, if r (t) is of the form ti 1(t), where i ≥ 0 and integer, then y(t) is unbounded. That is, a system with unrepeated poles on the jω-axis is BIBO unstable for a narrow class of signals. For this reason, a system with unrepeated poles on the jω-axis is called marginally stable. Marginal stability is the boundary between stability and instability. We do not attach any more meaning to it in this course.
15.3
An interesting example
Consider the TF G (s) = 1/(s − 1). In the light of the preceding discussion in this lecture, we see that this TF is unstable. Which of the following is a correct statement? 1. This TF is unstable because the response of this TF to no BIM input is BIM. 2. This TF is unstable because the response of this TF to not every BIM input is BIM. If we select the first one as the correct statement, then consider the BIM input whose TF is R(s) = (ss+−11)2 . It seems that the response to this BIM input is BIM1 . Let us examine this more carefully. Indeed, in this case, the response to this R(s) is Y (s) = 1/(s + 1)2 , that is, y(t) = te−t . So, the first choice is not the correct one. Is the second choice a correct one? Does the second choice follow from the definition of BIBO stability that we reproduce below from the previous lecture? An LTI system is said to be bounded-input bounded-output (BIBO) stable if and only if every bounded-in-magnitude (BIM) input r (t) results in a BIM output y(t). Yes, it does. This fact is an additional pointer to the fact that Theorem 2 is all that we need to prove.
1 This
input was constructed by Rohit Gupta, a Y8 student
January 2, 2012
78 of 218
Ramprasad Potluri
Lecture 16
Nyquist Stability Criterion (Part 1)* “Nyquist’s original analysis was fairly complicated, but his main result had a lasting value. It was later pointed out by MacColl [20] that the result can in fact be easily obtained from a theorem of Cauchy or from the Argument Principle of complex analysis” [21].
16.1
A quick recap
We saw in the previous lecture that the stability of a system described by a TF of the form (9.2) is determined by the location of the poles. In particular, if the poles are not in the open LHP, the system is not stable. Except in the case of unrepeated poles on the jω axis, if the poles are anywhere other than in the open LHP, the system is unstable. The systems that have unrepeated poles on the jω axis and all the remaining poles in the open LHP are said to be marginally stable as their response is bounded for all signals except those in a narrow class.
16.2
Problem statement
When we want to close the loop around a controller-plant-sensor combination, we are concerned if closing the loop may result in an unstable system. The Nyquist stability criterion (NSC), part of Nyquist stability theory, provides us a way to discuss the stability of the CL system based on the knowledge of the location of the open loop poles. Harry Nyquist developed this criterion in 1932. However, textbooks present a version of the NSC that seems to be MacColl’s version, apparently not realizing that credit is due to him. As MacColl apparently used Cauchy’s princple of argument, we shall review this principle first. But, we will present MacColl’s version of Nyquist’s stability criterion. 79
EE250 (Control Systems Analysis) IITK
16.3
Lecture Notes
Cauchy’s theorem (Principle of argument)
1. Suppose we have the TF Q(s) = (s + z)/(s + p). The pole is s = − p and the zero is s = −z. 2. Suppose some closed contour1 Γs in the s-plane encloses − p and −z, as shown in the following figure. jω
s-plane Γs
−p
−z
σ
Enclosing − p and −z means that the contour does not pass through − p and −z. 3. Take a point s on Γs and consider the angles θz and θ p made by the vectors s + z and s + p. jω s
s-plane Γs
θp
−p
θz
−z
σ
Then,
∠Q(s) = ∠(s + z) − ∠(s + p) = θz − θ p 4. Suppose we traverse Γs in the clockwise direction. For one clockwise traversal of Γs , the overall change in ∠Q(s) is ∆Γs ∠Q(s) = ∆Γs θz − ∆Γs θ p = −360◦ − (−360◦ ) = 0◦ 5. On the other hand, if Γs encloses only the zero and not the pole, as shown below, 1 In
this context, by closed we mean that the contour starts and ends at the same point.
January 2, 2012
80 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes jω
s
s-plane Γs
θp
−p
θz
−z
σ
then for one clockwise traversal of Γs , the overall change in ∠Q(s) is ∆Γs ∠Q(s) = ∆Γs θz − ∆Γs θ p = −360◦ − (0◦ ) = −360◦ 6. In general, if the s-plane contour Γs encloses Z zeros and P poles of the TF G (s) =
∏in=0 (s + zi ) ∏id=0 (s + pi )
then for one clockwise traversal of Γs , the change in the angle of G (s) is ∆Γs ∠G (s) = ( Z − P)(−2π ) 7. Another way of stating this is as follows. Each point on Γs is mapped by G (s) to a point on a new complex plane. This collection of points is a closed contour ΓG . Note that this mapping may be one-to-one or many-to-one, but not one-to-many. This is because G (s) is of the form of (9.2). We call this new complex plane the G (s)-plane, and label its axes Re{ G (s)} and Im{ G (s)}) to emphasize that on this complex plane, we draw only the map from Γs through G (s). s-plane
jω
G (s)
Re { G (s)} ΓG
Γs
−p
−z 0
σ
0
Re { G (s)}
G (s)-plane
For one complete clockwise traversal of Γs , ΓG encircles the origin of the G (s)plane N = Z − P times in the clockwise direction. For example, in the above figure, N = 2. Problem 16 Given the manner in which Γs has been constructed, can it be a contour that is not closed?
January 2, 2012
81 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
82 of 218
Lecture Notes
Ramprasad Potluri
Lecture 17
Nyquist Stability Criterion (Part 2)* 17.1
Application of Cauchy’s theorem in Nyquist’s stability criterion
Consider the unity feedback closed-loop system GCL (s) formed around a plant with transfer function G (s). GCL (s) =
G (s) 1 + G (s)
The stability of GCL (s) is decided by the location of its poles. We wish to determine if GCL (s) has any poles in the closed RHP. NSC helps determine if there are any poles of GCL (s) in the open RHP using only knowledge of G (s). The following two theorems are useful in the discussion. Theorem 3 The set of poles of GCL (s) equals the set of zeros of 1 + G (s). 2
Proof: Theorem 4 The set of poles of 1 + G (s) equals the set of poles of G (s).
Proof: 2 Theorem 3 allows us to convert the problem of examining the poles of GCL (s) into a problem of examining the zeros of 1 + G (s). NSC involves the following procedure. Construct an s-plane contour that encloses the entire closed RHP (i.e., the RHP including the jω-axis) while passing infinitesimally to the right of any poles of 1 + G (s) that are on the jω-axis. Construct the (1 + G (s))-plane contour corresponding to this s-plane contour. For one clockwise traversal of the s-plane contour, the (1 + G (s))-plane contour encircles the origin N = Z − P times in the clockwise direction. Here, Z is 83
EE250 (Control Systems Analysis) IITK
Lecture Notes
the number of zeros, and P is the number of poles, of 1 + G (s) enclosed by the s-plane contour. Alternatively, this means that the G (s)-plane contour encircles the point −1 + j0 in the clockwise direction N = Z − P times. Therefore, by knowing P and N, we can determine Z. From Theorem 4 we know P as the number of poles of G (s) enclosed by the s-plane contour.
17.2
Discussion
Theorems 3 and 4 may seem obvious and needing no proof. Then, why state them and invoke them? NSC was published in 1932. More than half a decade later since then, I find that the literature seems to have taken the NSC so much for granted that it is stated rather imprecisely. Occassionally, this imprecision may lead us into difficulties as is illustrated by Example 13. Example 13 In the previous lecture, I made one of the following statements (I don’t re call which) about the closed-loop system GCL (s) = G (s) (1 + G (s) H (s)) : 1. The poles of GCL (s) are the zeros of 1 + GH, 2. The zeros of 1 + GH are the poles of GCL (s). My intention was to claim that, to discuss the stability of GCL (s), it is adequate to examine the zeros of 1 + GH. If the first statement were true, then it would be appropriate for our purpose. However, Anurag Dash (Y9115) showed that for G = 1/(s + 5) and H = (s + 5)/(s − 5) the closed-loop poles are −5 and +4, while the only zero of 1 + GH is +4. So, the first statement is not true. Siddharth Sharma (Y9584) showed that for G = (s − 2)/(s − 4) and H = s/(s − 2), there are no closed-loop poles, while the only zero of 1 + GH is +2. So, neither the first nor the second statement is true. Yet, as seen in Table 17.1, certain famous textbooks do claim to the effect that for closed-loop stability, all roots of 1 + GH must lie in the LHP.
While I do some investigation to settle the applicability of NSC to GCL (s) of Example 13, we shall work with the case of H (s) = 1. The applicability of NSC for this case is clear from Section 17.1.
January 2, 2012
84 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Table 17.1: Examples of the imprecision found in the literature in discussing Nyquist stability criterion. Statement and source
Remarks
Consider the function F (s) = 1 + G0 (s)C (s). “. . . the zeros of F (s) are the closed-loop poles in a unity feedback control system” [22, page 140].
Denote by Z the set of zeros of F (s), and by P the set of the closed-loop poles of the unity feedback control system (that is, G0 C/(1 + G0 C )). Then, this statement only means that Z ⊆ P , It does not mean that P ⊆ Z or that P = Z . To apply NSC, it would have been better to have P ⊆ Z . Theorem 3 shows that P = Z , which is adequate.
“Also, the poles of F (s) are the open-loop poles of plant and controller” [22, page 140].
Denote by P F the set of poles of F (s), and by POL the set of open-loop poles. Then, this statement only means that P F ⊆ POL . Theorem 4 showed a stronger result, namely that P F = POL . However, this statement is adequate for use with NSC.
“The closed-loop transfer function is G/(1 + GH ). For stability, all roots of the characteristic equation 1 + GH = 0 must lie in the left-half s plane.”[10, page 462].
The tacit assumption here is that the set of the roots (that is, zeros) of 1 + GH is a superset of the set of poles of the closed-loop transfer function. As we see in Example 13, this may not always be true.
“the closed-loop transfer function of a SISO system is M(s) = G (s) (1 + G (s) H (s)) , where G (s) H (s) can assume the form (9.2). Since the characteristic equation is obtained by setting the denominator polynomial of M(s) to zero, the roots of the characteristic equation are also the zeros of 1 + G (s) H (s). Or, the characteristic equation roots must satisfy ∆(s) = 1 + G (s) H (s) = 0.” [23, page 557 – 558], [24, page 427]
The assumption is that the set of the poles of the closed-loop transfer function equals the set of the zeros of 1 + GH. As we see in Example 13, this may not always be true.
“we consider the closed-loop transfer function Y (s)/R(s) = G/(1 + GH ) . . . The characteristic equation of the closed-loop system is obtained by setting the denominator of Y/R to zero, which is the same as setting the numerator of 1 + GH to zero. Thus, the roots of the characteristic equation must satisfy 1 + GH = 0.” [16, page 375]
This version is the same as the one in the previous row.
January 2, 2012
85 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
86 of 218
Lecture Notes
Ramprasad Potluri
Lecture 18
Nyquist Stability Criterion (Part 3) We apply the Nyquist stability criterion (NSC) to three sample unity feedback closed-loop systems. The TF of the closed-loop system is GCL (s) = G (s) (1 + G (s)) . For G (s) we consider K/s, K/(s + 1), and K (s + 10)2 s3 .
18.1
Example 1: G (s) = K/s
Here are the steps in applying NSC. 1. Mark on the s-plane the poles of G (s) that are in the closed RHP. In this case, G (s) has only one pole in the closed RHP, and this pole is at the origin. 2. Draw the s-plane contour as shown in Figure 18.1. Note that, for the contour to not pass through the poles on the jω-axis, it makes slight detours around these poles at an infinitesimally small distance away from them. In the figure, the section of the s-plane contour near the origin has taken a detour to the right of the pole. In practice, it can also take a detour to the left of the pole. 3. Use the Bode plot (BP) of G (s) to map section C1 to the G (s)-plane. Note that several textbooks map the jω-axis to the G (s)-plane using a few points taken from this axis. The reader may verify that this way of mapping is tedious in comparison to the BP way of mapping. 4. The section at infinity of the s-plane contour maps via G (s) into an infinitesimally small area around the origin of the G (s)-plane. As this area does not affect the count of the encirclements, we say that “the section at infinity of the s-plane contour of G (s) maps into the origin of the G (s)-plane.” 5. The section of the s-plane contour around the origin maps into the section at infinity of the G (s)-plane contour (same as NP of G (s)). We can sketch this contour with just 1 – 2 points. 87
EE250 (Control Systems Analysis) IITK
Lecture Notes
jω
s-plane C2
C1 C6 C5 C4
jφ
ρe
jθ
Re
σ ρ→0 R → ∞ Γs C3
Figure 18.1: The s-plane contour for G (s) = K/s. This contour has six sections of interest, C1 — C6. Note that, by construction, the s-plane contours for any G (s) are symmetric about the real axis. In this figure, the section of the s-plane contour near the origin has taken a detour to the right of the pole, making P = 0. If it had taken a detour to the left of the pole, then P = 1. To estimate the location of the image of C6 under G (s), we select a convenient point, marked ◦ by •, which, in this case is ρe j45 .
6. Here is another fact that simplifies sketching NPs. As the s-plane contour is symmetric about the σ-axis (i.e., Re{s}-axis), then the NP too is symmetric about the Re{ G (s)}-axis. The reason for this fact is that for the TFs of the form (9.2), G (s) = ( G (s∗ ))∗ , where x ∗ denotes the complex conjugate of x. 7. The complete G (s)-plane contour is shown in Figure 18.2.
18.2
Convention followed in sketching Nyquist plots
Multiple poles are marked no matter what the multiplicity. Sections of the s-plane contour are labeled C1, C2, · · · . The corresponding sections on the NP are identified by the labels C10 , C20 , · · · . Points on the s-plane contour are labeled 1, 2, · · · . The corresponding points on the NP are labeled 10 , 20 , · · · . The section of the NP that maps from the jω-axis section of the s-plane contour is called polar plot. However, by “polar plot” some books refer to the entire NP. The sections of the NP that map from the sections of the s-plane contour that are above the σ-axis (i.e., the Re{s}-axis) are shown solid. Their mirror images about the real axis are shown dashed. January 2, 2012
88 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Im { G (s)} G (s)-plane C50 C40
−1 + j0
C30 C20
− jφ Ke ρ
Re { G (s)} ΓG
C10
C60
Figure 18.2: The Nyquist plot (NP) (also know as “G (s)-plane contour”) for G (s) = K/s. This contour has six sections of interest, C10 — C60 . Note that the G (s)-plane contours are symmetric about the real axis. In this figure, the sections C20 and C30 are at the origin. The • on section C60 is the image of the • on section C6 in Figure 18.1. In our procedure, the section C10 is constructed first. This section is the polar plot of G (s), and for convenience is shown thicker than the other sections. We note that the NP of K/s does not encircle −1 + j0. Therefore, N = 0, and, consequently, Z = N + P = 0. This means that the unity feedback closed-loop system built around K/s does not have any open RHP poles for an K ≥ 0.
January 2, 2012
89 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
18.3
Example: G (s) = K (s + 1)
18.4
Example: G (s) = K (s + 10)2 s3
Lecture Notes
Problem 17 Sketch the NP of the TFs K s ; ; s2 ( s2 + 2) s + 1
18.5
Appendix: A useful fact
A necessary condition for a polynomial an s2 + an−1 sn−1 + . . . + a1 s + a0 to have roots only in the LHP is that ai > 0, i = 0, 1, . . . , n.
January 2, 2012
90 of 218
Ramprasad Potluri
Lecture 19
Nyquist Stability Criterion (Part 4) 19.1
Gain margin, phase margin
The concept of relative stability is used to express how close a control system is to instability. To evaluate the relative stability of the unity-feedback CL system, the terms gain margin (GM) and phase margin (PM) are used frequently with respect to Nyquist plots (NPs) that have simple forms such as those shown in Figure 19.1. Note that the transfer functions (TFs) corresponding to these NPs usually have simple forms such as K ∏m j =1 ( s + z i ) s ∏in=1 (s + pi )
e−td s
where K, td , pi and zi are real. We define GM and PM using Figure 19.2. The GM of this TF is defined as the amount by which G ( jω ) can be amplified at the phase crossover frequency ω p before G ( jω p ) passes through −1 + j0. From this definition, it follows that GM = 1 G ( jω p ) . The PM of this TF is the angle through which the polar plot needs to be rotated in the clockwise direction before G ( jω g ) passes through −1 + j0. ω g is known as the gain crossover frequency of this TF. From this definition, it follows that PM = G ( jω g ) − (−180◦ ). It has been found for many systems that GM 1 and PM > 30◦ help keep the CL system stable in the face of parameter variations. The GM is also expressed in logarithmic terms as 20 log10 ( GM), in dB, as shown in Figure 19.2. 91
EE250 (Control Systems Analysis) IITK
Lecture Notes
Figure 19.1: Nyquist plots for which concepts of GM and PM are relevant. Note that the transfer functions (TFs) corresponding to these NPs usually have simple forms such as
19.2
K ∏m j =1 ( s + z i ) − t d s e , s ∏in=1 (s+ pi )
where K, td , pi and zi are real.
When GM and PM are useful
Consider the minimum-phase TF [11] 85(s + 1)(s2 + 2s + 43.25) + 2s + 82)(s2 + 2s + 101)
s2 ( s2
The polar plot of this TF is shown in Figure 19.3. The following GNU Octave code generates this figure. plot((-1:0.01:1),sqrt(1-(-1:0.01:1).^2),’-r’) hold, plot((-1:0.01:1),-sqrt(1-(-1:0.01:1).^2),’-r’) nyquist(tf(conv(85*[1,1],[1,2,43.25]),conv([1,2,82,0,0],[1,2,101]))) axis([-2,1,-1.5,1.5],"square") We see that the multiple intersections of the polar plot with the unit circle make it difficult to determine the PM of this transfer function. When the NP has multiple crossings of the unit circle or of the −180◦ axis, then we may not be able to convey much information by talking in terms of GM and PM. Other examples of situations where it is better to view the Nyquist plot as a whole instead of discussing relative stability through GM and PM are the January 2, 2012
92 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
G ( jω p )
20 log10 | G ( jω )|
ωp 0
−1 + j0
ωg
φ
ωp
0 dB 20 log10 G ( jω p )
ωg
−90◦ 6
G ( jω g )
−180◦ 6
−270◦
Figure 19.2: Illustration for the concepts of PM and GM. The left subfigure shows the NP of a certain minimum phase TF, while the right figure shows the approximate BP of the same TF. We see from the NP that, for the stability of this closed-loop system, φ = −180◦ − ∠G ( jω g ) > 0 and 1 G jω p > 1 are desirable. φ is the PM. From the BP, we see that, while the PM is still expressed as −180 ◦ − ∠G ( jω g ), the GM is expressed in dB. Specifically, the GM is 0 − 20 log10 G ( jω p ) dB. In this figure, the GM in absolute terms works out to a value greater than 1, and to a positive value in dB terms. following TFs: K ( s + 1) ; s s ( − 1) 10
19.3
K (s + 10)2 . s3
Appendix: Dead-time element in MATLAB
How to specify TF of a First-Order Plus Dead Time Model in MATLAB is explained at http://www.mathworks.com/ products/control/demos.html?file=/products/demos/shipping/control/ January 2, 2012
93 of 218
Ramprasad Potluri
G ( jω )
EE250 (Control Systems Analysis) IITK
Lecture Notes
Nyquist plot from u1 to y1, w (rad/s) in [1.000000e-01, 1.000000e+03] 1.5
+w -w
1
Im( G(jw) )
0.5
0
-0.5
-1
-1.5 -2
-1.5
-1
-0.5
0
0.5
1
Re( G(jw) )
Figure 19.3: An example of a complicated Nyquist plot. The red circle is the unit circle centered at the origin. Of the Nyquist plot, only the polar plot (shown by the thin solid blue line) and its mirror image (shown by the ‘+’ line) are shown. The multiple intersections of the polar plot with the unit circle make it difficult to determine the PM of this transfer function. This TF is 85(s+1)(s2 +2s+43.25) . s2 (s2 +2s+82)(s2 +2s+101)
GSSpecifyingDelays.html#2.
January 2, 2012
94 of 218
Ramprasad Potluri
Lecture 20
Nyquist Stability Criterion (Part 5) 20.1
Demonstration Nyquist plots
of
quick
construction
of
Consider some sample transfer functions for G (s) of Figure 20.1 s−1 s−1 (s + 10)2 (s + 10)2 ; ; ; . s(s − 10) s(s + 10) s2 s3
20.2
Critical gain
The point −1 + j0 is known as critical point, and the value of K for which the Nyquist plot of the system of Figure 20.1 passes through the critical point is known as critical gain. Example 14 shows how to evaluate the critical gain. Example 14 Consider the TF
s −1 . s(s−10)
Yref (s) +
E(s)
KG (s)
Y (s)
− Figure 20.1: Unity feedback closed-loop system. The parameter K is to be designed through analysis such that this system is stable.
95
EE250 (Control Systems Analysis) IITK
20.3
Lecture Notes
Conditionally stable systems
In practice, many plants are such that, the unity feedback closed-loop system formed around them goes closer to instability as the controller gain K of Figure 20.1 increases. A PMDC motor is one such plant. K (s+10)2
. The unity feedback CL However, there are plants with TFs such as s3 systems corresponding to these TFs are examples of conditionally stable systems. They become more stable as K increases. In the example of
K (s+10)2 , s3
we see that (
N=
2 0
K is small K is large.
Thus, as the number of open-loop poles P = 0, by NSC, the number of closedloop poles in the open RHP is ( 2 K is small Z= 0 K is large.
20.4
Using − K1 + j0 instead of −1 + j0
From the examples, which we saw thus far, of applying NSC, we see that the position of polar plot of KG (s) relative to −1 + j0 is the same as the position of the polar plot of G (s) relative to − K1 + j0. This observation helps simplify using NSC. We need to only sketch the polar plot of G (s), and not that of KG (s).
20.5
Nyquist plots in MATLAB and GNU Octave
Consider some sample transfer functions for G (s) of Figure 20.1 s−1 (s + 10)2 (s + 10)2 s−1 ; ; ; . s(s − 10) s(s + 10) s2 s3 The following code generates the Nyquist plots of these TFs in MATLAB and GNU Octave. Figures 20.2 and 20.3 show the output of this code. subplot(2,2,1), nyquist(tf([1,-1],[1,-10,0])), xlim([-0.2,0.05]), ylim([-1,1]); title(’(s-1)/(s(s-10))’); subplot(2,2,2), nyquist(tf([1,-1],[1,10,0])), xlim([-0.1,0.2]), ylim([-0.5,0.5]); title(’(s-1)/(s(s+10))’); subplot(2,2,3), nyquist(tf([1,20,100],[1,0,0])), title(’(s+10)^2/s^2’); subplot(2,2,4), nyquist(tf([1,20,100],[1,0,0,0])), xlim([-1.5,0.1]), ylim([-0.5,0.5]); title(’(s+10)^2/s^3’); January 2, 2012
96 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
(s−1)/(s(s−10))
(s−1)/(s(s+10)) 0.5 Imaginary Axis
Imaginary Axis
1 0.5 0
0
−0.5 −1 −0.2
−0.15
−0.1 −0.05 Real Axis
0
−0.5 −0.1
0.05
0
(s+10)2/s2
0.2
(s+10)2/s3 0.5 Imaginary Axis
50 Imaginary Axis
0.1 Real Axis
0
−50 −600
−400
−200 Real Axis
0
200
0
−0.5 −1.5
−1
−0.5 Real Axis
0
Figure 20.2: Nyquist plots of sample values of G (s) of Figure 20.1 generated by MATLAB.
We see that in the name of Nyquist plots, we only have the polar plots and their mirror images about the real axis. This picture may not seem enough to talk about encirclements of the critical point. We may wonder how to use this picture. To use this picture, we need to perform one additional step. We need to sketch the sections at infinity of the Nyquist plots.
20.6
Rotation of a vector instead of encirclement
Instead of visualizing as the encirclement of the origin by the 1 + G (s) contour, we can visualize as the turning of the vector 1 + G (s), with its tail at the origin and head on the 1 + G (s) contour, around the origin. Similarly, instead of visualizing as the encirclement of the origin by the G (s) contour, we can visualize as the turning around −1 + j0 of a vector that has its tail at −1 + j0 and head on the G (s) contour. Such a visualization may help where we are unable to visualize encirclement. January 2, 2012
97 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
(s-1)/(s(s-10))
(s-1)/(s(s+10))
1
+w -w
+w -w
0.4
0.5
Im( G(jw) )
Im( G(jw) )
0.2
0
0
-0.2 -0.5
-0.4 -1 -0.2
-0.15
-0.1
-0.05
0
0.05
-0.1
-0.05
0
0.05
0.1
Re( G(jw) )
Re( G(jw) )
(s+10)2/s2
(s+10)2/s3
+w -w
20
0.15
0.2
+w -w
0.4
15
0.2
5
Im( G(jw) )
Im( G(jw) )
10
0
0
-5 -0.2
-10 -15
-0.4 -20 -100
-80
-60
-40
-20
0
Re( G(jw) )
-1.4
-1.2
-1
-0.8
-0.6
-0.4
-0.2
0
Re( G(jw) )
Figure 20.3: Nyquist plots of sample values of G (s) of Figure 20.1 generated by GNU Octave.
January 2, 2012
98 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
20.7
Lecture Notes
Discussion
For those TFs that have one or more poles at the origin, the section at infinity of the s-plane contour (this section is also known as the D-section) maps into an infinitisimally small section of a circle around the origin of the Nyquist plane. This section does not affect the count of the encirclements. So, we can safely assume that the D-section maps into the origin of the Nyquist plane in this case. When the OL TF does not have a pole at the origin, the D-section does not map into the origin any more. The Nyquist plane is also known as the G (s)-plane.
January 2, 2012
99 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
100 of 218
Lecture Notes
Ramprasad Potluri
Part IV
Design using bode plots-based loop-shaping for minimum-phase plants
101
Lecture 21
NST (concluded), Loop-shaping (Part 1) 21.1
Nyquist stability theory (conclusion)
Examples 15 and 16 wrap up our discussion of Nyquist stability theory. Example 15 Sketch the NP of K/(s2 ) and discuss the stability of the unity feedback closed-loop system formed around this TF.
Example 16 Sketch the NP of K/(s2 + 1) and discuss the stability of the unity feedback closed-loop system formed around this TF.
When the NP passes through the critical point, then N is unknown.
21.2
Specifications for control system design
The closed-loop (CL) system should 1. be stable, 2. be robust to plant parameter variations, disturbances (usually mechanical such as wind gusts on vehicle, tyre-road friction, etc), noises (usually electromagnetic interference (EMI)), 3. demonstrate desired performance (such as fast rise time (tr ), limited overshoot (M p ), fast settling time (ts ), and small steady-state error (ess )), and 4. apply limited control input to the plant. 103
EE250 (Control Systems Analysis) IITK
Lecture Notes
+ R(s) +
E(s)
+
GOL (s)
−
D (s) Y (s)
+
N (s) +
Figure 21.1: A unity feedback closed-loop system formed around the openloop transfer function GOL (s).
R(s) +
Cascade configuration R(s) +
E(s)
C (s)
U (s)
G (s)
Y (s)
− E(s)
C1 (s)
−
U (s)
+ +
G (s)
Y (s)
C2 (s)
Minor loop feedback
Figure 21.2: Two sample configurations of plant and controller.
21.3
Loop gain and shaping
In classical control theory, the design technique called loop-shaping provides a collection of graphical tools by which we can design the CL system’s properties as desired by modifying the open-loop (OL) system’s transfer function (TF). This modification is called “shaping” because this TF is modified by shaping the corresponding Nyquist/Bode/Nichols plot of this TF. We already saw an instance of discussing the CL system’s properties via the OL system’s TF in the form of Nyquist stability criterion (NSC). Consider the CL system of Figure 21.1. GOL (s) is called the loop gain of the control system. See Section 6.2 for details. Given that the plant gain G (s) is fixed, the loop gain can be as desired if we design a controller for the plant. The plant-controller configuration can be of many types. Two configurations are shown in Figure 21.2, ignoring disturbances and noises. The configurations we can have are probably only limited by our creativity. In this course, we will consider only cascade configuration. However, the techniques developed may be applied to other configurations through a suitable modification of the block diagram. January 2, 2012
104 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Slope: −20 dB/dec
1 decade
0◦
1
−90◦ 2
−180◦
3
−270◦
Figure 21.3: Illustration for Example 17.
21.4
BMP-BPP correspondence for minimum-phase TFs
For minimum phase TFs, we can uniquely construct the BPP given the BMP. A theorem by Hendrick Bode says that for a minimum-phase TF, if a section of the BMP persists with a slope of n for at least one decade, then in the middle of that section, the BPP has a value of n × 90◦ [25, pages 291-293]. Here n = . . . , −2, −1, 0, +1, +2, . . ., where −2 ≡ −40 dB/dec, −1 ≡ −20 dB/dec, etc. Example 17 Consider the BMP and BPPs shown in Figure 21.3 for a certain minimumphase TF. Which of the three BPPs are valid candidates?
21.5
Appendix: Bode’s theorem
21.6
Appendix: NST for SFG defintion of loop gain
Applying the discussion of Section 6.2 to the SFG equivalent shown in Figure 21.4 of the block diagram of Figure 21.1, the loop gain is bc. Thus, if we used the signal flow graph definition, we should have called − GOL (s) the loop gain. Then NSC would have been modified to examine the January 2, 2012
105 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Figure 21.4: The signal flow graph equivalent of the block diagram of Figure 21.1.
encirclement of the point +1 + j0 by the Nyquist plot of − GOL (s). Indeed, “the critical point in Nyquist’s paper is +1 instead of −1 as is commonly used today” [21].
January 2, 2012
106 of 218
Ramprasad Potluri
Lecture 22
Loop-shaping (Part 2) 22.1
Recap of specifications for control system design
The following properties are desired of a CL system 1. Stability, 2. Robustness to plant parameter variations, 3. Robustness to disturbances, 4. Robustness to sensor noise, 5. Demonstration of desired performance (transient as well as steady-state), and 6. Respect for constraints on control input to the plant.
22.2
Gist of loop-shaping-based design
As we will consider only cascade configuration in our course, GOL (s) = C (s) G (s), where C (s) is the controller TF, and G (s) is the plant TF. C (s) needs to be designed, and G (s) is fixed. Loop-shaping using Bode plots for minimum-phase TFs involves the following simple steps: 1. Draw the ABMP of the plant (G) on a semilog graph paper. 2. Draw the ABMP of the desired OL TF (GOL ) on the same semilog graph paper 3. Subtract the first ABMP from the second. This gives the ABMP of the controller. 4. Write the TF of the controller. 107
EE250 (Control Systems Analysis) IITK
22.3
Lecture Notes
Stability
For minimum-phase TFs with simple-looking Nyquist plots (NPs), we saw that gain margin (GM) and phase margin (PM) were adequate measures of how far the NP was from encircling the −1 + j0 point. Thus, PM and GM are adequate measures of stability. For minimum-phase TFs that do not have second order elements, taking care of PM takes care of CL stability. Problem 18 Cook up a minimum-phase TF that includes a second order component such that the unity-feedback closed-loop system built around this TF has BP and NP that show adequate PM but not adequate GM.
22.4
Robustness to plant parameter variations
Plant parameters variations may bring the CL system closer to instability. Provide adequate PM so that the CL system remains stable in spite of coming closer to instability. How much PM is adequate? A value of 60◦ has been found to be adequate for many systems occurring in practice. A PM greater than this makes the system relatively more stable but also sluggish. How to obtain a PM less than 90◦ ? Though the statement of Section 21.4 talks in terms of phases which are integral multiples of 90◦ , Bode’s theorem says more. For example, if a section of the BMP persists with a slope of n for a certain fraction of a decade, then in the middle of that the BPP has a value that is proportional to this fraction of 90◦ , as illustrated in Figure 22.1. Loop-shaping Recipe 1 From Bode’s theorem discussed in Section 21.4, it follows that, to have adequate PM, make the BMP of GOL (s) have a −20 dB/dec section that is centered at ω g and has appropriate width. Doing so provides the CL system both stability and robustness to plant parameter variations.
22.5
Relation between OL and CL BMP
The closed-loop (CL) TF from input r (t) to output y(t) is GCL (s) =
GOL (s) 1 + GOL (s)
The GOL ( jω ) designed by loop-shaping usually has a high low-frequency gain and a low high-frequency gain, while being minimum-phase. It follows that ( 1, ∀ω ω g | GCL ( jω )| = | GOL ( jω )|, ∀ω ω g January 2, 2012
108 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Slope: −20 dB/dec 1 decade
1 2 0◦
−90◦ −180◦ −270◦
Figure 22.1: If a section of the BMP persists with a slope of −20 dB/dec for less than a decade, then in the middle of that section, the phase is appropriately away from −90◦ . This figure shows the ABMPs in contrast to Figure 21.3, which showed the BMP. Note that if a section of the BMP has to persist for a decade with a slope of, for example, −20 dB/dec, then the same section on the ABMP needs to persist for more than a decade with a slope of −20 dB/dec. ABMP 1 of this figure corresponds to the BMP of Figure 21.3. Note that the PM shown by phase plot 1 is +90◦ , while that shown by phase plot 2 is about +50◦ .
dB-gain
For a minimum-phase TF, in the vicinity of the gain crossover frequency ω g , where | GOL ( jω g )| = 1, | GCL ( jω g )| depends heavily on the phase margin (PM)[25, pages 295-296], as illustrated by Example 18. The form of 20 log10 | GCL | for two values of PM (45◦ and 90◦ ) is shown in the following figure (for the moment, ignore the Low ω, Mid ω, and Hi ω bands).
0 dB
GOL
GCL
PM = 45◦ ωg PM = 90◦
Low ω
January 2, 2012
Mid ω
109 of 218
ω (log scale)
Hi ω
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Example 18 Evaluate GCL ( jω g ) for PM = 90◦ and PM = 45◦ . We have GOL ( jω ) GCL ( jω ) = 1 + GOL ( jω ) Note that GOL ( jω g ) = 1. Thus, we only need to evaluate 1 + GOL ( jω g ) . The following figure shows 1 + GOL ( jω ) graphically for the case of PM = 90◦ (left subfigure) and PM = 45◦ (right subfigure).
G (s)-plane
G (s)-plane 1
)
L
) ωg
1
)
√ 1 + GOL ( jω g ) = 2 √ ⇒ GCL ( jω g ) = 1/ 2 = 0.7
0
( jω g G OL
(j
L
GO
GO
g
1+
( jω
1+
GOL ( jω g )
0
1 + GOL ( jω g ) ≈ 0.75 ⇒ GCL ( jω g ) = 1/0.75 = 1.33 More accurately, 1 + GOL ( jω g ) = r 2 2 √1 + 1 − √1 ≈ 0.77 2
2
Historically, through calculations such as done in Example 18, | GCL ( jω )| were found around ω g , and the BMPs of GCL (s) were sketched.
January 2, 2012
110 of 218
Ramprasad Potluri
Lecture 23
Loop-shaping (Part 3) 23.1
Satisfaction of desired performance
The Performance of a control system is expressed through how well the control system tracks the test signals 1(t), t1(t), t2 1(t), etc. Which of these signals is used to measure the performance depends on which signals the control system most frequently encounters in practice. The control system does not actually encounter any of these test signals in practice. However, for a given plant, one or more of these test signals serve to compare the performance of one control system built around this plant with the peformance of another control system built around this plant. These test signals are usually idealized versions of the signals that the control system will encounter in practice. The unit step response presents the most pretty and clear picture of the performance of a control system. So, this is most often shown in books. Once we become comfortable with the unit step response, talking about the performance to a ramp or parabolic input is not very difficult. The unit step response of a general control system is shown in Figure 23.1. The parameters of this response that are of interest to us are • Percent overshoot M p defined as Mp ,
y(t p ) − y(∞) × 100% y(∞)
• Settling time ts is the time taken by the system to enter the x% tube with the intention of remaining in it. Here, x% of y(∞) is the radius of the tube. The tube is centered about the steady-state value y(∞) of the unit step response. • Steady state error ess = 1 − y(∞) is a measure of the steady-state accuracy of the control system. For example, in a position control system, if, after the transients die out, the position achieved is different from what was needed to be achieved, we say that this control system has a non-zero ess . • Rise time tr . How do loop-shaping techniques help meet these performance specs? The following subsections answer this question. 111
EE250 (Control Systems Analysis) IITK
Lecture Notes
1 y(∞) x% tube
ess
0
tr
ts
tp
Diameter = 2x%
y(t)
t
Figure 23.1: The unit step response of a control system. First order response and second order responses are special cases of this response.
23.1.1
Steady-state error
Consider the following equation: E(s) =
1 1 1 R(s) − D (s) − N (s) 1 + GOL (s) 1 + GOL (s) 1 + GOL (s)
When R(s) = 1/s (that is, r (t) is unit step), we can write the following using final value theorem: ess , lim e(t) = lim sE(s) = lim s t→∞
s →0
s →0
1 1 1 = 1 + GOL (s) s 1 + GOL (0)
Therefore, GOL (0) ≈ 1/ess . Therefore, to reduce ess , increase the LF value of | GOL |. We can define klf , 20 log10 | GOL (0)|. Then, klf can define the LF forbidden zone seen in the earlier figure.
23.1.2
Settling time
Consider the unit impulse response of the first order TF τs1+1 . This has the equation 1 y(t) = e−t/τ 1(t) τ where, τ is the time constant of this response. This response reaches x% of its initial value of 1 in time t x% . Therefore, t x% = τ ln
100 x
For a x% tube centered about the y = 0 line, t x% is the settling time ts . The following table summarizes the settling time ts for standard values of the tube:
January 2, 2012
x
1
2
5
10
ts
4.6τ
3.9τ
3τ
2.3τ
112 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
From the table, we can say that transients in a first order system die out in approximately 5τ, that is, in five time constants. Therefore, for a first order system, ts ≈ 5τ. For a first order system, the TF τs1+1 shows that ωb = 1/τ. Therefore, ωb ≈ 5/ts . We can carry this approximation to higher order systems too, and say that transients in any system die out — if they die out — in approximately five times the largest time constant τmax . Then, assuming that ωb ≈ 1/τmax , the bandwidth of any system can be defined roughly as being related to the system’s ts as ωb ≈ 5/ts . The inequality ω g ≤ ωb ≤ 2ω g , implies that ω g ≈ ωb . Therefore, in order to provide the desired ts , we need to set ω g ≈ 5/ts .
23.1.3
Overshoot
The rule of thumb is, “the closer the system is to instability, the more oscillatory its step response becomes, and the more will be its overshoot.” Of course, this statement may not be applicable to a first order system. But, most practical systems are not purely first order. If they are modeled as of first order, then that may be because their higher order dynamics have been neglected. To limit the value of M p , provide suitable PM. We will see soon the relation between the desired maximum M p and the appropriate PM. The measures we adopted in subsections 22.3 and 22.4 and are also adequate to limit M p .
23.2
Respecting constraints on control input
We have the following equation for the cascade configuration of the controller: U (s) =
C (s) C (s) C (s) R(s) − D (s) − N (s) 1 + GOL (s) 1 + GOL (s) 1 + GOL (s)
This equation shows that, to keep the initial spike in u(t) within limits, keep 20 log10 |C ( jω )| small at the appropriate frequencies. What are these appropriate frequencies? Section 23.3 provides an answer to this question.
23.3
Appendix: Additional loop-shaping tips
This appendix starts where subsection 23.2 left off. These additional loopshaping tips use the following facts: 1. The final value theorem relates the low-frequency section of the loop-gain to the area of the time-domain response that is near t = ∞ lim y(t) = lim sY (s)
t→∞
s →0
2. The initial value theorem relates the high-frequency section of the loop-gain to the area of the step response that is near t = 0+ lim y(t) = lim sY (s)
t →0+
January 2, 2012
s→∞
113 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
3. From the previous two points, I am conjecturing that, as we travel right to left on the ω axis, we travel left to right on the t axis of the time-domain response. We summarize this idea through the following rule of thumb: lim y(t) =
t→0++
lim sY (s).
s→∞−−
Here, t = 0+ is the time instant immediately to the right of t = 0, while t = 0 + +, t = 0 + ++, etc are time instants progressively further to the right of t = 0+. s = ∞−, s = ∞ − −, etc carry similar meaning.
How the loop gain affects M p of x (t) and magnitude of u(0+). Here u x = u.
t ux t
20 log10 | Gdes ( jω )|
x (t) xd
ωu ωl
ωg
ω
Finally, the additional tips are as follows: 1. To reduce |u(0+)| reduce | Gdes ( jω )| immediately after ωu , as shown in the above figure, while also keeping the loop gain at higher frequencies small as shown. This can be seen using initial value theorem as follows: u(0+) =
C (∞−) 1 + C (∞−) G (∞−)
assuming R(s) = 1/s. So, making |C (∞−)| small will make |u(0+)| small. Note that, the choice of ωu and ωl itself involves trial and error in practice. 2. To reduce M p (that occurs at t = 0 + +) of x reduce | Gdes ( jω ) in the region much below high frequencies, but above the low frequency region, as shown in the above figure. This can be seen using the above-mentioned conjecture as follows: lim y(t) ≈
t→0++
C (s) G (s) s→∞−− 1 + C ( s ) G ( s ) lim
assuming R(s) = 1/s. So, making |C (∞ − −)| small will make M p small. 3. To reduce ts of x, we have two options: 3.1. Increase ω g . 3.2. Decrease M p as mentioned in item 2 above.
January 2, 2012
114 of 218
Ramprasad Potluri
Lecture 24
Loop-shaping (Part 4) 24.1
Robustness to disturbances
Any undesired phenomenon acting on the plant can be modeled as disturbance irrespective of whether it is electrical, mechanical or chemical in nature. Disturbances are usually relatively slow phenomena. Disturbance needs to be processed quickly. So, the system needs to be quicker than the disturbance. So, have the bandwidth of the CL system much larger than disturbance frequencies. So, design so that disturbance lies in the low frequency (LF) region of the CL system’s passband. So, designate the frequencies in which the disturbances lie as the “Low ω” frequencies as shown in Figure 24.1. “Processing” the disturbance involves suppressing its effect on the output y(t). This can be done as follows. Consider the following the equation: Y (s) =
GOL (s) 1 GOL (s) R(s) + D (s) − N (s) 1 + GOL (s) 1 + GOL (s) 1 + GOL (s)
We see that in order to reject disturbances, we need to make | GOL ( jω )| large in the Low ω region. How large? Keep the low frequency (LF) loop gain above a certain LF forbidden zone that lies in the interval [0, ωd ].
24.2
Robustness to noises
Noises are usually electrical in nature, such as the high frequency EMI induced by transistors switching in a AC-DC or DC-DC converter circuit. They are also of much higher frequencies than disturbances. We see from the last equation above that keeping the noises within the bandwidth of the CL system will not be the right thing to do, because the noise n(t) will then appear unmitigated in the output y(t). So, the strategy adopted is to select ωb smaller than the noise frequencies. That is, keep the noises outside the passband of the CL system. The inequality ω g ≤ ωb ≤ 2ω g , implies that ω g ≈ ωb . This means that | GOL ( jω )| needs to be small in the Hi ω region. How small? Keep the high frequency (HF) loop gain below a certain HF forbidden zone that lies in the interval [ωn , ∞). 115
EE250 (Control Systems Analysis) IITK
Lecture Notes
khf
ωu ωd
ωn
ωg ωl −20 dB/dec
HF forbidden zone
LF forbidden zone
klf
Figure 24.1: Summary of loop-shaping. The width of the −20 dB/dec section centered at ω g accounts for the CL properties of stability, robustness to disturbance rejection. This width, and, consequently, these properties are related to the phase margin (PM). The height of the low frequency (LF) forbidden zone accounts for steady-state accuracy of tracking an input signal, and for rejecting disturbances. The right side boundary of this zone is at a frequency ωd that is sufficiently small in comparison to ωb (the closed-loop bandwidth), and, consequently, to ω g . This way, disturbances are tackled quickly. The high frequency (HF) forbidden zone accounts for sensor noise. The lower level of this zone is at a value of khf dB; the amplification of the sensor noise does not rise above this level.
24.3
Summary of loop-shaping recipe
The loop-shaping technique we developed in these lectures is valid for minimum-phase GOL (s). The loop-shaping tips developed thus far are summarized by Figure 24.1 [25, pages 322]. Note that ω g being in the center of the −20 dB/dec section implies that √ ω g = ωl ωu because ωl and ωu are logarithmically spaced. Note that, if the −20 dB/dec section centered at ω g is at least one decade wide on the BMP, then the PM will be 90◦ . If the section is one decade wide on the ABMP, then the PM will be ≈ 60◦ . Problem 19 Why are we not using GM in our design?
24.4
MATLAB-based example of loop-shaping
24.4.1
Problem
Design the TF GOL (s) such that the unity-feedback closed-loop system’s unit step response satisfies the following specs: M p ≤ 30%; ess ≤ 0.05; ts = 5 s for a 2% tube. January 2, 2012
116 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
24.4.2
Lecture Notes
Solution
1 ess ≤ 0.05 implies that klf = 20 log10 | GOL (0)| & 20 log10 = 26 dB. ess 5 ≈ 1 rad/s. ts ◦ M p ≤ 30% can be achieved with a PM ≈ 60 . This translates into a −20 dB/dec section that is approximately one decade wide centered at ω g . Next, we will use MATLAB to generate a BMP that will have these three features. ts = 5 s for a 2% tube can be achieved with ω g ≈
24.4.3
How to use MATLAB
1. In MATLAB, we cd to the directory that contains the files a.m and cleanup.m that are shown below. 2. We type the word magshape at the command prompt. This launches the GUI tool MAGSHAPE. Note that MAGSHAPE presents us a semilog grid where the x-axis represents frequencies in rad/s, and y-axis represents dB-gain. 3. In MAGSHAPE — 3.1. We declare a filter named test. 3.2. We select the add point radio button and add points on the GUI’s grid approximately along the desired BMP that we have in mind. 3.3. We specify the filter order as 3. We could have specified other filter orders too. We want to have a filter of the smallest order that will match our desired BMP. 3.4. We press the fit data button. This generates a filter of order 3 that best approximates the filter that we specified through the points. 3.5. If the generated filter does not fit the specified filter properly, then select the move point or delete point radio button, move or delete the points, and fit data once again. Repeat this cycle until MAGSHAPE generates a filter that whose BMP resembles what we want. 4. The filter generated by MAGSHAPE is now in MATLAB’s workspace. Run the m-file a.m that I wrote. a.m does the following: 4.1. Converts the filter generated by MAGSHAPE into a TF form using the function ltitf provided by MATLAB. 4.2. Plots this filter’s BP. 4.3. Plots the unit step response of the CL system corresponding to this filter. 5. If the step response does not satisfy the problem’s specs, then repeat the entire exercise beginning with step 3. January 2, 2012
117 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
24.4.4
Lecture Notes
Results
A few trials of the above algorithm gave us the filter shown in the MAGSHAPE window of Figure 24.2. Figure 24.3 shows the Bode plot of GOL (s). Figure 24.4 shows the unit step response of the unity-feedback CL system built around GOL (s). The step response approximately satisfies the three specs given by the problem.
24.4.5
MATLAB code
% a.m: File to capture results from MAGSHAPE % and plot various plots. This file is to % be used after using MAGSHAPE. % % What I did in MAGSHAPE: I declared a filter % named "test". Having selected the radio % button corresponding to "test", I added % many points that I wanted my Gol (which is % the desired loop gain) to fit, and gave 3 % as the filter order (I could have given % other numbers too). Then, I clicked on % "fit data" button. This generated the % filter named "test" that I am using in % the following. %------------------------------------------close all % Close any open figure windows % (except MAGSHAPE window). [num,den]=ltitf(test); % % % % % % % % % %
I found that num or den may contain some terms that are about 0.0001 or less times less in absolute value compared to the largest-in-absolute-value term of num (den). These terms caused Gol (earlier obtained as Gol=tf(num,den)) to contain high-order nondominating terms that made it a non-proper TF. So, I wrote a little function "cleanup" that helps remove the non-dominating high order terms from the TF of "test".
[n,d]=cleanup(num,den); gol=tf(n,d); gcl=feedback(gol,1,-1); % Plot the OL and CL BPs, % and the CL step response. % January 2, 2012
118 of 218
Ramprasad Potluri
Lecture Notes
Figure 24.2: Snapshot of MAGSHAPE showing the BMP of GOL that gives the desired CL response.
EE250 (Control Systems Analysis) IITK
January 2, 2012
119 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
t = linspace(0,18,200); % % The above line was absent in the a.m % mailed to the students immediately after the % lecture. I addded it later because I found % that the step reponse plotted by MATLAB was % jagged. Jagged responses are not natural. % Adding the above line gives a smooth response. % step(gcl,t), grid print -depsc step.eps % figure, % bode(gcl), grid, title(’BP of Gcl’) figure, bode(gol), grid, title(’BP of Gol’) print -depsc OLbode.eps % cleanup.m function [n,d]=cleanup(num,den) numerr = 0.0001; denerr = 0.0001; % In future, we can specify the values of % numerr \& denerr based on various factors % such as the orders of the numerator and % denominator, and maybe others that I am % not seeing at the moment. for i=1:size(num,2) if abs(num(i))/abs(max(den)) < numerr num(i) = 0; end end for j=1:size(den,2) if abs(den(j))/abs(max(den)) < denerr den(j) = 0; end end n = num; d = den; Problem 20 Construct an ordinary differential equation of the smallest order possible so that the solution of this differential equation to a unit step forcing function has the form shown in Figure 23.1. Assume M p ≈ 30%, ess ≈ 0.04, ts ≈ 1 January 2, 2012
120 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
BP of Gol 50
Magnitude (dB)
0
−50
−100
−150 0
Phase (deg)
−45
−90
−135
−180 −3 10
−2
10
−1
10
0
1
10 10 Frequency (rad/sec)
2
10
3
10
4
10
Figure 24.3: BP of GOL . The PM is what we expected based on the width of the section of slope −20 dB/dec centered at ω g .
s. The permitted tolerance for your solution is 5%. Multiple overshoots and undershoots are not mandatory.
January 2, 2012
121 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Step Response 1.4
1.2
Amplitude
1
0.8
0.6
0.4
0.2
0
0
2
4
6
8 10 Time (sec)
12
14
16
18
Figure 24.4: Unit step response of the CL system. The step response approximately satisfies the three specs given by the problem.
January 2, 2012
122 of 218
Ramprasad Potluri
Lecture 25
Loop-shaping (Part 5): Information from lead/lag controllers 25.1
Background
We are interested in designing a controller C (s) for the cascade combination shown in Figure 21.2. Figure 25.1 shows the steps, which are involved in the design and deployment of a controller. In the flow chart, Step7 involves tuning. Tuning is required because the controller may have been designed for an idealized model of the plant. The model may have ignored nonlinearties and certain dynamics of the plant. So, the controller may not work right unless it is tuned. The fewer the number of parameters to be tuned in practice, the easier the tuning, as we can track the effect of each parameter on the control system. The simplest controller is a proportional controller. But, a proportional controller is limited in what it can help the closed-loop system achieve. When a controller is more than a proportional controller, we call it a compensator. A transfer function (TF) that does not contain terms in s is called a static TF. A TF that has terms in s is called a dynamic TF. A compensator contains terms in s; so it is dynamic. Some of the simplest controllers are shown in Table 25.1 In this lecture, we study lead and lag controllers. We look at the PID controller and its special cases (that is, PI and PD) after the second midsemester exam.
25.2
Lead and lag controllers
The lead controller provides a positive phase, while the lag controller provides a negative phase as shown in Figure 25.2. We can develop the relationship between the maximum positive phase Φmax (or the minimum negative phase Φmin ) provided by a lead (or lag) controller and the corresponding decade distance (DD) between the corner fre123
EE250 (Control Systems Analysis) IITK
Lecture Notes
Start Exit Step1 Mathematically model the plant as LTI ODE
Step2 Analysis of the plant model (Is it minimum phase? Is it proper or strictly proper? What is its order? What are its poles and zeros?)
Step3 Formulate design goals in terms of which signals to track, with what accuracy, what disturbances to reject & how much, how much to reject sensor noises & at what frequencies?
YES
Review step 3
Mechanism behaves well?
Step8 Deploy the control system in the full mechanism whose component it is YES
Review the steps that lead here
Step4 Convert design goals into loopshaping spec-s such as PM, ω g , forbidden zones
CL system behaves as desired?
Step7 Tune the controller so that CL system behaves as desired
Step5 Design & simulate controller until satisfactory on paper
Step6 Program the controller in µC / DSP
Figure 25.1: Flow chart of the design and deployment of a controller.
January 2, 2012
124 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Table 25.1: Some of the simplest controllers. Proportional
kp
Proportional plus derivative (PD)
k p + kd s
Proportional plus integral (PI)
k p + k I /s
Proportional plus integral plus derivative (PID)
k p + k d s + k I /s Ts + 1 , α<1 K αTs + 1 Ts + 1 , α>1 K αTs + 1
Lead Lag
−20 dB/dec
+20 dB/dec 0◦
ωmax
90◦
Φmin
Φmax 0◦
1 T
−90◦
1 αT
1 αT
ωmax
1 T
Figure 25.2: Bode plots of lead and lag controllers. The lead controller provides a positive phase, while the lag controller provides a negative phase.
quencies of the controller through the steps shown in Figure 25.3. See Section 25.4 for details. Table 25.2 tabulates the required DD (∆) between the corner frequencies of the lead (or lag) controller and the corresponding Φmax (or Φmin ) that the controller needs to provide for certain commonly occuring values of Φmax (or Φmin ).
25.3
How DD-Φmax relationship for lead/lag controllers is useful in loop-shaping
We can use the insights obtained so far for the lead/lag controller in our loopshaping techniques. Specifically, I claim that if the desired minimum-phase √ loop gain is as shown in Figure 24.1, with ω g = ωl ωu , then the desired phase margin (PM) in radians can be provided by choosing the DD ∆ between ωl and ωu as 1 + sin PM decades, (25.1) ∆ = log10 1 − sin PM with PM expressed in degrees. Thus, for a given PM, we can choose the required DD between ωl and ωu from the above table by replacing Φmax with PM. January 2, 2012
125 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK Determine ωmax as a function of T and α
Lecture Notes
Determine Φmax as a function of α using ωmax
Determine α as a function of Φmax
Determine ∆ as a function of Φmax
Determine ∆ as a function of α
Figure 25.3: Steps in the development of the relation between the decade distance (DD, ∆) and Φmax . See Section 25.4 for details.
Table 25.2: The required DD (∆) between the corner frequencies of the lead (or lag) controller and the corresponding Φmax (or Φmin ) that the controller needs to provide for certain commonly occuring values of Φmax (or Φmin ). Note that the DD is the same on the BMP and the AMBP. Φmax (Φmin ) needed
90◦ (−90◦ )
60◦ (−60◦ )
45◦ (−45◦ )
30◦ (−30◦ )
DD between corner
∞
1.1439
0.76555
0.47712
≈ 1.144
≈ 0.766
≈ 0.477
freq-s of BP in decades
A plot of ∆ versus PM can be generated using the code shown in Figure 25.4 and is shown in Figure 25.5. This code uses Equation 25.1. Example 19 Design using loop-shaping techniques on the 5-cycle semilog paper provided a controller of the minimum order possible to control the speed of a moω (s) tor with the transfer function (TF) V (s) = s100 +20 for the following time domain specifications: steady-state error to a step input ess ≤ 2%, settling time ts ≈ 0.5 s, peak overshoot M p ≤ 20%. Label your axes appropriately. Show how the specifications are reflected in your loop-shaping construction. Mark ω g suitably. The solution is shown in Figure 25.6.
25.4
Appendix: Derivation of DD versus PM relationship
January 2, 2012
126 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
% ddpm.m: Plots DD verus PM. % Works in GNU Octave and Matlab. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all, close all, pm = linspace(0,pi/2,100); D = log10((1+sin(pm))./(1-sin(pm))); plot(pm*180/pi,D); grid; xlabel(’PM (in degrees)’); ylabel(’Decade Distance (in decades)’); Figure 25.4: The code that generates Figure 25.5. 4.5
4
Decade Distance (in decades)
3.5
3
2.5
2
1.5
1
0.5
0
0
10
20
30
40 50 PM (in degrees)
60
70
80
90
Figure 25.5: Decade distance (DD) versus phase margin (PM). We note that the DD is fairly proportional to the PM upto about 60◦ with a relationship of DD = PM/60.
Derivation of DD versus PM relationship Lead compensator C (s)
January C (2, jω2012 ) =
Ts + 1 , α<1 αTs + 1 K s + T1 1 α s + αT
= K =
Φ dΦ ⇒ dω
Lag compensator
K
jωT + 1 jωαT + 1
dy/dx
C ( jω )
Φ dΦ ⇒ dω
dΦ
Ts + 1 , α>1 αTs + 1 K s + T1 1 α s + αT
= K =
127 of 218
= arctan Tω − arctan αTω T αT = − 1 + T2 ω2 1 + α2 T 2 ω 2
d arctan y
C (s)
= K
jωT + 1 Potluri Ramprasad jωαT + 1
= arctan Tω − arctan αTω T αT = − 1 + T2 ω2 1 + α2 T 2 ω 2
EE250 (Control Systems Analysis) IITK
Lecture Notes
Figure 25.6: The solution of Example 19. January 2, 2012
128 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Derivation of DD versus PM relationship (concluded) Lead compensator Plugging the expression for ωmax in that for Φ, we obtain Φmax
= arctan Tωmax − arctan αTωmax √ 1 = arctan √ − arctan α α √ 1 √ − α α = arctan[ √ ] 1 + √1α α
⇒ Φmax ⇒ Φmax
Lag compensator
⇒ Φmin ⇒ Φmin ⇒ Φmin
1−α = arctan √ 2 α
= arctan Tωmin − arctan αTωmin √ 1 = arctan √ − arctan α α 1−α = arctan √ 2 α ⇒ tan Φmin =
1−α √ 2 α
1−α √ 2 α
⇒ tan Φmax =
Using the relation sin θ = √ tan θ 2 , 1+tan θ we can write Φmax alternatively as follows: sin Φmax =
⇒α=
∆ = log10
= log10 = log10 E.g., for Φmax = √ ∆ = log10
3 √2 1− 23
sin Φmin =
1 − sin Φmax 1 + sin Φmax
The DD (∆) between
1+
1−α 1+α 1 − sin Φmin ⇒α= 1 + sin Φmin
1−α 1+α
1 T
and
1 αT
is
= log10
and
1 αT
is
1 1 − log10 T αT = log10 α > 0, α > 1 1 − sin Φmin = log10 1 + sin Φmin √
∆ needs to be 1+0.877 1−0.877
1 T
∆ = log10
1 1 − log10 αT T 1 > 0, α < 1 α 1 + sin Φmax 1 − sin Φmax
+60◦ ,
The DD (∆) between
E.g., Φmin =
−60◦
: ∆ = log10
⇒ ∆ = 1.1439 decades
3 √2 1− 23
1+
⇒ ∆ = 1.1439 decades.
January 2, 2012
129 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
130 of 218
Lecture Notes
Ramprasad Potluri
Lecture 26
Loop-shaping (Part 6): Relation between M p and PM In this lecture, for the closed-loop (CL) system of Figure 21.2, we see the relation between the desired overshoot M p in the system’s unit step response and the phase margin (PM) that would provide this M p . For this purpose, we develop this relationship for a second order system, and assume this relationship as being approximately valid for systems of higher order.
26.1
Relation between M p and PM for second order system
Consider the standard second order TF G (s) =
ωn 2 s(s + 2ζωn )
(26.1)
Sample Bode and polar plots that correspond to this standard second order TF, for ωn = 1 rad/s and ζ = 0.5, plotted using the m-file ord2.m included in this lecture note are shown in figures 26.1 and 26.2 respectively. We show these plots so that the PM of this TF is clear. The CL system is, of course, a unity-feedback one. Let us work out the PM of this TF. For this, we first determine the gain crossover frequency ω g by setting the magnitude of G ( jω ) = 1 at ω g . ωn 2 G ( jω g ) = −ω g 2 + j2ζωn ω g = 1 This gives ω g = ωn
rq
As
∠G ( jω ) = −
4ζ 4 + 1 − 2ζ 2
π ω − arctan 2 2ζωn 131
EE250 (Control Systems Analysis) IITK
Lecture Notes
Bode Diagram
Magnitude (dB)
50
0
−50
Phase (deg)
−100 −90
−135 PM
−180 −2 10
−1
0
10
1
10 Frequency (rad/sec)
10
2
10
Figure 26.1: Bode plot of the standar open-loop second order system G (s) = ωn 2 . The item of interest in this plot is the phase margin of this open-loop s(s+2ζωn ) system.
90
4
120
60 3 2
150
30
1
180
0
210
330
240
300 270
Figure 26.2: Polar plot of the standar open-loop second order system G (s) = ωn 2 . The item of interest in this plot is the phase margin of this open-loop s(s+2ζωn ) system.
January 2, 2012
132 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK we have
∠G ( jω g ) = −
Lecture Notes
ωg π − arctan 2 2ζωn
The famililar expression PM = π + ∠G ( jω g ). gives us, using the fact that
π 2
− arctan x = arctan 1x ,
PM = arctan qp
2ζ
(26.2)
4ζ 4 + 1 − 2ζ 2
Let us see how this PM is related to the peak overshoot M p (expressed as M p = y(t p )−y(∞) y(∞)
× 100%) of the unit step response of the unity-feedback CL system G (s) =
ωn 2 . + 2ζωn s + ωn 2
s2
that corresponds to the above OL TF. Problem 21 Show that, for the second order unity-feedback CL system built around the TF of (26.1), M p is given by the expression: Mp = e
√−πζ
1− ζ 2
× 100%
(26.3)
A plot of M p versus ζ using Equation (26.3) and a plot of PM versus ζ using Equation (26.2) is shown in Figure 26.3, and a plot of PM versus M p based on the plot of Figure 26.3 is shown in Figure 26.4, using the m-file pmmp.m that is included in this lecture note. Also, plotted in Figure 26.3 is the linear approximation (denoted PM∼) PM ≈ 100 × ζ,
for 0 ≤ ζ . 0.707
to the PM vs. ζ curve. In the figure, we see that PM ≈ 100ζ is a reasonable approximation to PM given by Equation (26.2) only for 0 ≤ ζ . 0.707. This approximation provides us the PM-to-M p approximate correspondence of Table 26.1 and Figure 26.4 for a second order system. Finally, we are interested in the DD-to-M p correspondence. This correspondence is shown in Figure 26.5 by putting together information from figures 26.3 and 26.4.
26.2
Applying the second order relation to higher order
Henceforth in our loop-shaping methodology, we will use the approximate correspondence of Figure 26.5 as a guideline to choosing the width of the −20 January 2, 2012
133 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
100 PM (in °) Mp (in %)
90
PM∼ (in °)
80
PM, Mp, PM∼
70
60
50
40
30
20
10
0
0
0.1
0.2
0.3
0.4
0.5 ξ
0.6
0.7
0.8
0.9
1
Figure 26.3: Relationship between M p and ζ and PM and ζ. Note that the approximation shown by PM∼ is only valid for 0 ≤ ζ . 0.707.
100 PM PM∼
90
80
PM, PM∼ (in degrees)
70
60
50
40
30
20
10
0
0
10
20
30
40
50 M_p (in \%)
60
70
80
90
100
Figure 26.4: Relationship between M p and PM. Note that the PM∼ versus M p relationship is only valid for PM between 0 and 70◦ .
January 2, 2012
134 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
2 ∆ PM/100 1.8
DD (in decades) & PM (in degrees)
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0
10
20
30
40 50 60 Peak overshoot M (in %)
70
80
90
100
p
Figure 26.5: Plot of DD versus M p constructed from figures 26.3 and 26.4. Also reproduced is the information from Figure 26.4 of PM versus M p .
Table 26.1: The approximate correspondence of PM to M p of the standard second order system. PM
25◦
40◦
60◦
70◦
Mp
45%
25%
10%
5%
dB/dec section centered about ω g . This figure is constructed from figures 26.3 and 26.4. Note that this approximation seemed to have worked in the MATLAB-based example given in Section 24.4. There we have a PM of about 60◦ giving an M p of about 10%. See figures 24.3 and 24.4. Note that this approximation is only an approximation. It may work very well if the system is close to second order. It may not work very well for higher order systems. For example, in Section 23.3 we mention that the section of the Bode magnitude plot of Gdes ( jω ) in the region much below high frequencies, but above the low frequency region, also has a bearing on M p . Problem 22 Verify that the design of Example 19 meets its specifications, and see what information you could have applied to that design from this lecture.
% ord2.m: m-file to plot Bode and polar plot of January 2, 2012
135 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
% the 2nd second order system open-loop system. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; close all; clc; % Define the standard OL 2nd order system’s % coefficients: wn = 1; zeta = 0.5; s=tf(’s’); g=(wn^2)/(s*(s+2*zeta*wn)); % Plot the BP of the OL 2nd order system: bode(g); grid; % Plot the polar plot of the OL 2nd order system: n=1000; w = logspace(-0.6,2,n); [mag,ph] = bode(g,w); % % The mag and phase outputs of bode command in MATLAB % even when the system is Single Input Single Output, % are 3-D arrays. The polar command of Matlab is unable % to handle these arrays. So we convert them into 1-D % arrays as follows: % for i=1:n mag1(i) = mag(:,:,i) ; ph1(i) = ph(:,:,i)*pi/180; % converted degrees to radians end % figure; polar(ph1,mag1) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % pmmp.m: m-file to see the relationship % between PM of the second order system % $\frac{{\omega_n}^2}{s(s+2\xi\omega_n)}$, % and the peak overshoot $M_p$ of the unit % step response of the corresponding % closed-loop system. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; close all; xi = linspace(0,1,40); % Declare an array of values of xi at which to % evaluate pm = atan(2*xi./sqrt(sqrt(4*xi.^4+1)-2*xi.^2)); January 2, 2012
136 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
% This gives PM in radians. Convert this into % degrees: pm = pm*180/pi; % % % % % % %
Note the construct "./". The dot preceding the arithmetic operations *, /, and ^ allows to perform elementwise multiplication, division, and raising to a power respectively. For example, [1, 2].*[2, 3] gives [2, 6]. [1,2]./[2,5] gives [0.5,0.4], and [1,2].^2 gives [1,4].
mp = exp(-pi*xi./sqrt(1-xi.^2))*100; % This gives the peak overshoot in percent % Peak overshoot is defined as % $M_p = \frac{y(t_p)-y(\infty)}{y(\infty)}\times 100\%$. pmapprox = 100*xi; plot(xi,pm,xi,mp,xi,pmapprox) xlabel(’\xi’); ylabel(’PM, M_p, PM\sim’) legend(’PM (in \circ)’,’M_p (in %)’,’PM\sim (in \circ)’) grid print -depsc pm-zeta-mp.eps figure plot(mp,pm,mp,pmapprox) xlabel(’M_p (in \%)’); ylabel(’PM, PM\sim (in degrees)’); legend(’PM’,’PM\sim’); grid print -depsc mppm.eps Problem 23 If necessary, edit ord2.m and pmmp.m so that these files run successfully in GNU Octave, giving the same results that they gave in MATLAB.
January 2, 2012
137 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
138 of 218
Lecture Notes
Ramprasad Potluri
Lecture 27
Loop-shaping (Part 7): A solved example 27.1
Problem
Consider the control system of Figure 27.1. We wish to design a compensator D (s) that satisfies the following design specifications: (a) Kv = 100. (b) PM ≈ 60◦ . (c) Sinusoidal inputs of up to 1 rad/sec to be reproduced with ≤ 2% error. (d) Sinusoidal inputs with a frequency of greater than 100 rad/sec to be attenuated at the output to ≤ 5% of their input value. 1. Determine K1 and K2 of Figure 27.2. 2. Write the numerical values of ω1 and ω2 of Figure 27.2. 3. What is the decade distance needed between the corner frequencies ωl and ωh for the desired Bode plot of Figure 27.2? 4. Is this DD the distance on the BMP or on the ABMP or both? 5. On the semilog grid provided, draw the ABMPs of the desired D (s) G (s) and of G (s). Your figure must contain all the necessary labels. 6. On the semilog grid provided, show the ABMP of the resulting D (s). Write the TF of D (s). 7. For the resulting CL system, given that ωB ∈ [ωmin , ωmax ], where ωB is the bandwidth, what are the values of ωmin and ωmax ?
139
EE250 (Control Systems Analysis) IITK Yref (s) +
−
D (s)
Lecture Notes
U (s) G (s) =
10 s s +1 10
Y (s)
Figure 27.1: Control system to be designed.
20 log10 (K1 )
ωg ω1
ωl
ωh −20 dB/dec
ω2
20 log10 (K2 )
Figure 27.2: Parameters for loop-shaping for the problem.
27.2
Solution
We first translate the specifications into actionable information. (a) Kv = 100. At this point in the course, we have not yet seen the concept of Kv . We will do so soon. In the meanwhile, it suffices to know that for a unity-feedback control system built around the TF GOL (s) that contains a single integrator, Kv is the frequency where the left-most asymptote of the ABMP of GOL (s) has a value of 0 dB. Kv has the units of rad/s. In this problem, Kv = 100 implies that the left most asymptote of the BMP of D (s) G (s) has a slope of −20 dB/decade and passes through the point { ω = 100 rad/s, dB-gain = 0 dB }. (b) PM ≈ 60◦ . This specification can give us the decade distance (DD) between the ωl and ωh of Figure 27.2. Table 25.2 says that a DD of 1.144 decades between ωl and ωh with a −20 dB/decade decline between them implies PM ≈ 60◦ . We have from this specification: log10 ωh − log10 ωl = 1.144 January 2, 2012
140 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
(c) Sinusoidal inputs of up to 1 rad/sec to be reproduced with ≤ 2% error. E(s) 1 = R(s) 1 + D (s) G (s) 1 ⇒ Ge ( jω ) = 1 + D ( jω ) G ( jω )
Ge (s) ,
Let us define e(t) , r (t) − y(t). For an input sinusoid r (t) = A sin ωt, the steady-state value of e(t) is ess (t) = | Ge ( jω )| A sin(ωt + ∠Ge ( jω )). It is given that ess (t) ≤ 2% when r (t) is upto 1 rad/s in frequency. Is this a specification for the magnitude and phase of ess (t), or is it only for the phase of ess (t)? Note that, the frequency responses of even ideal filters are defined as having phases changing linearly with frequency. So, most likely, we are not being asked to constrain the phase as well as magnitude of a practical filter within 2%. This guess is verified by looking at specification (d). So, the 2% specification can only be on the magnitude of ess (t). So, we shall constrain the magnitude of ess (t) within 2%. This can be done as follows:
| Ge ( jω )| =
1 ≤ 2%, |1 + D ( jω ) G ( jω )|
∀ω ≤ 1 rad/s. ⇒ | D ( jω ) G ( jω )| ' 50, ∀ω ≤ 1 rad/s. So, K1 ≥ 50 and ω1 = 1 rad/sec. (d) Sinusoidal inputs with a frequency of greater than 100 rad/sec to be attenuated at the output to ≤ 5% of their input value. Y (s) D (s) G (s) = R(s) 1 + D (s) G (s) D ( jω ) G ( jω ) ⇒ GCL ( jω ) = 1 + D ( jω ) G ( jω )
GCL (s) ,
Want | GCL ( jω )| ≤ 5% (wanting to keep ∠GCL ( jω ) ≤ 5% is unrealistic). This can be achieved thus:
| D ( jω ) G ( jω )| / 5%, ∀ω ≥ 100 rad/s So, K2 / 0.05 and ω2 = 100 rad/s. Now that we have actionable information, we answer the questions. 1. Determine K1 and K2 of Figure 27.2. K1 = 50, K2 = 0.05 2. Write the numerical values of ω1 and ω2 of Figure 27.2. ω1 = 1 rad/sec, ω2 = 100 rad/sec January 2, 2012
141 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
3. What is the decade distance needed between the corner frequencies ωl and ωh for the desired Bode plot of Figure 27.2? 1.144 decades. 4. Is this DD the distance on the BMP or on the ABMP or both? Both. 5. On the semilog grid provided, draw the ABMPs of the desired D (s) G (s) and of G (s). Your figure must contain all the necessary labels. See graph of Figure 27.3. 5.1. From the values of K1 , K2 and ω1 , ω2 we get the forbidden regions. The BMP of D (s) G (s) cannot lie in this region. 5.2. Draw the left-most asymptote of the BMP of D (s) G (s). For this, using Kv = 100, draw a straight line at a slope of −20 dB/decade and passing through 100 rad/sec at 0 dB. 5.3. As no data is given to fix ω g , we choose ω g to be approximately equidis√ tant, on logarithmic scale, from ω1 and ω2 . That is, ω g ≈ ω1 ω2 . This gives ω g = 10 rad/sec. 5.4. We draw a 1.144 decade-wide −20 dB/decade straight line centered at ωg . 5.5. We connect the left end of this −20 dB/decade section to the asymptote of step 5.2 via a straight line while being careful that this line does not intersect the LF forbidden region. This line will need to have a slope that is integer multiple of −20 dB/decade. We get this slope as −60 dB /dec in our case. 5.6. From the right end of the −20 dB/decline (of step 5), we draw a straight line at a slope of −40 dB/dec (integral multiple of 20), that passes below the HF forbidden zone. The ABMP of the desired D (s) G (s) is shown in the graph of Figure 27.3. Also shown on the same graph is the ABMP of G (s). 6. On the semilog grid provided, show the ABMP of the resulting D (s). Write the TF of D (s). See graph of Figure 27.3. 6.1. Subtract the ABMP of G (s) from that of D (s) G (s) to obtain the ABMP of D (s) as shown on the graph. 6.2. The TF of D(s) is s 2 s 10 +1 +1 D (s) = 2.5 2 10 . s s +1 +1 0.7 40 January 2, 2012
142 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
7. For the resulting CL system, given that ωb ∈ [ωmin , ωmax ], where ωb is the bandwidth, what are the values of ωmin and ωmax ? Since ω g ≤ ωb ≤ 2ω g , we have ωmin = ω g = 10 rad/s, and ωmax = 2ω g = 20 rad/s.
s 2 s 10 +1 +1 Figure 27.3: Solution to the problem. D (s) = 2.5 2 10 . s s +1 +1 0.7 40
January 2, 2012
143 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
144 of 218
Lecture Notes
Ramprasad Potluri
Lecture 28
Metrics for steady-state accuracy of unity-feedback systems built around minimum-phase systems Consider the closed-loop system of Figure 21.1. Assuming that d(t) and n(t) are absent, the error in tracking a reference input R(s) is given by E(s) = R(s) (1 + G (s)). How accurately a control system tracks, for example, a step or a ramp input is of interest in the design. More generally, test signals of the form r (t) = tn 1(t) n! may be used to evaluate the performance of a control system. Note that L− { tm m!} = 1 sm+1 . So, R(s) = 1 sm+1 . See Section 28.4 for the definition of L− {·}. The steady-state error ess is a measure of this accuracy, and is defined as ess , lim e(t). t→∞
We can use the final value theorem (FVT) to calculate ess as follows: ess = lim e(t) = lim sE(s) = lim s t→∞
s →0
s →0
1 R(s) 1 + G (s)
Note that the FVT works only if ess to a given input r (t) exists and is finite. Equivalently, FVT requires all poles of sE(s) to be in the open left half s-plane. Problem 24 In practice, where is a step or a ramp used as a test signal for a control system?
The TFs that we study in this course are of the form shown in (9.2). Assume that G (s) of (9.2) is minimum-phase. This means that all the constants in G (s) are non-negative, and that G (s) is rational. In the remainder of this lecture, we 145
EE250 (Control Systems Analysis) IITK
Lecture Notes
study the concepts of error constants and system type, which are used to discuss the steady-state accuracy of minimum-phase systems.
28.1
Error constants
Example 20 Evaluate ess for the case of a minimum-phase G (s) having 2 poles at the origin and the CL system tracking a parabolic input. Note that, though G (s) has 2 poles at the origin, if we wish to apply 2 FVT, then G (s) is not permitted 2 to be, for example, of the form 1/s . For m + 1 R(s) = 1 s , G (s) = 1 s would result in sE(s) having a pair of complex conjugate poles on the jω-axis, disqualifying this sE(s) for FVT. In the rest of this example, we assume that G (s) is such that sE(s) is amenable to FVT. ess = lim s s →0
1 1 1 1 = lim 2 = 3 2 1 + G (s) s s →0 s + s G ( s ) lims→0 s2 G (s)
As G (s) has exactly two poles at the origin and is minimum-phase, s2 G (s) will have all poles and zeros in the open left half s-plane. Therefore, lims→0 s2 G (s) is a positive finite constant.
We summarize in Table 28.1 calculations similar to those in Example 20 for the cases of r (t) being a step, ramp, and parabola, and G (s) having up to 2 poles at the origin, while being such that sE(s) has all poles in the open LHP. Problem 25 Show a G (s), which is minimum-phase, but the sE(s) corresponding to which has some poles that are not in the open left half s-plane.
In Table 28.1, K p , Kv , Ka are called error constants or error coefficients (ECs). K p is called position (or step) constant or position (or step) error constant, Kv is called velocity (or ramp) (error) constant, and Ka is called acceleration (or parabolic) (error) constant. We may wonder why control engineers would want to talk about steadystate error as the inverse of some quantity. We provide an answer to this question here. ECs are attributes of the loop gain of a control system when this loop gain is minimum-phase. They relate the loop gain of a unity feedback control system to the steady-state error of the control system under various inputs such as step, ramp, parabola, etc. ECs are especially obvious on the Bode magnitude plot (BMP) of the loop gain. Indeed, their importance stems from the fact that they can be read easily off the BMP []. Control system design specifications on the steady state error are provided through ECs, particularly when designing using Bode plots[]. Given the frequency response of the open loop system, we can fit a TF to this frequency response, and an important parameter in this TF is an EC. Therefore, an EC is a finite quantity. January 2, 2012
146 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Table 28.1: Calculations of the steady-state error (ess ) in tracking step, ramp, and parabolic inputs by a unity-feedback closed-loop system whose loop gain has upto 2 poles at the origin. Number of poles of G (s) at origin 0 lims→0 s 1+G1 (s) 1s = 1 . Here, 1+lims→0 G (s)
0 < lim G (s) < ∞ s →0
R(s)
1 s
∵ G (s) is minimum-phase. Denote it by K p . Then, ess = 1+1K p .
1 s2
1 s3
28.2
Cannot use FVT to evalute ess (∵ sE(s) has 1 pole at origin). Determine y(t). Then, e(t) = r (t) − y(t) gives ess = ∞. Cannot use FVT to evalute ess . Determine y(t). Then, e(t) = r (t) − y(t) gives ess = ∞.
1
2
Here, lim G (s) = ∞.
Here, lim G (s) = ∞.
∴ ess = 0.
∴ ess = 0.
s →0
s →0
lims→0 s 1+G1 (s) s12 = 1 . lims→0 sG (s)
Here, 0 < lim sG (s) <
Here, lim sG (s) = ∞.
∞. Denote it by Kv . Then, ess = K1v .
∴ ess = 0.
s →0
Cannot use FVT to evalute ess . Determine y(t). Then, e(t) = r (t) − y(t) gives ess = ∞.
s →0
lims→0 s 1+G1 (s) s13 = 1 . Here, lims→0 s2 G (s) 0 < lim s2 G (s) < s →0 ∞. Denote it by Ka . Then, ess = K1a .
How to read K p , Kv , Ka off the BMPs
The ECs can be read off BMPs easily as shown in the Figure 28.1 and as justified by Table 28.2.
28.3
System type
As we saw, ess depends on the input that the CL system is trying to track, as well as the number of integrators (that is, poles at the origin) present in G (s). The latter number is assigned the special name “system type”. Definition 4 [23, page 371] The unity-feedback CL system built around G (s) is type i (i = 0, 1, 2, . . .) if G (s) contains i integrators, i.e., if G (s) contains the element 1/si . Thus, we have type 0, type 1, type 2, etc systems. The definition of system type is much confused in literature[]. For example, which system is type i, G (s) or its CL version? Books also define the system as type 0 if its position error constant is finite, as type 1 if its velocity error constant is finite, etc. In my view, these definitions are roundabout versions of Definition 4. January 2, 2012
147 of 218
Ramprasad Potluri
OL TF has no poles or zeros at the origin
OL TF has exactly one pole at the origin
OL TF has exactly two poles at the origin
Table 28.2: Explanation for why the error constants can be read off the Bode magnitude plots.
n
d
2 2 ∏i=1 0 ( pi ) ∏k= 0 ( ωdk )
d
2 2 ∏i=1 0 (zi ) ∏k= 0 ( ωnk )
n
!
Information from the BMP of the minimum-phase version of G (s) of (9.2)
On the BMP of G (s), the equation dB-gain = 20 log10
that is, at ω = Kv .
ω=
n
d
n
d
2 2 ∏i=1 0 (zi ) ∏k= 0 ( ωnk )
n
,
!
2 2 1 ∏i=1 0 (zi ) ∏k= 0 ( ωnk ) ω 2 ∏d1 ( pi ) ∏d2 (ωdk 2 ) i =2 k =0
n
2 2 ∏i=1 1 ( pi ) ∏k= 0 ( ωdk )
On the BMP of G (s), the equation dB-gain = 20 log10
n
2 2 ∏i=1 0 (zi ) ∏k= 0 ( ωnk )
n
,
describes the left-most asymptote. This asymptote intersects the ω-axis at
!
describes the left-most asymptote. This asymptote intersects the ω-axis at
n n2 2 1 ∏i=1 0 (zi ) ∏k= 0 ( ωnk ) ω ∏d1 ( pi ) ∏d2 (ωdk 2 ) i =1 k =0
describes the left-most asymptote. This observation helps in reading K p off the BMP as equal to G (0).
dB-gain = 20 log10
ECs calculated for the minimum-phase version of G (s) of (9.2) using Table 28.1
n
On the BMP of G (s), the equation
n
d
2 2 ∏i=1 0 ( pi ) ∏k= 0 ( ωdk )
d
2 2 ∏i=1 0 (zi ) ∏k= 0 ( ωnk )
If a unit step isinput to the CL system, then we say ess = 1 (1 + K p ) , where K p = G (0) =
n n2 2 ∏ 1 ( zi ) ∏ k = 0 ( ωnk ) = di=0 d
2 2 ∏i=1 1 ( pi ) ∏k= 0 ( ωdk )
If a unit ramp is input to the CL system, then we say ess = 1 /Kv , where Kv = lim sG (s) s →0
n n2 2 ∏ 1 ( zi ) ∏ k = 0 ( ωnk ) = di=0 d
2 2 ∏i=1 2 ( pi ) ∏k= 0 ( ωdk )
If a unit parabola is input to the CL system, then we say ess = 1 /Ka , where Ka = lim s2 G (s) s →0
ω2 =
d d2 2 ∏i=1 2 ( pi ) ∏k= 0 ( ωdk ) √ Ka .
that is, at ω 2 = Ka , that is, at ω =
Ramprasad Potluri
148 of 218
January 2, 2012
Lecture Notes EE250 (Control Systems Analysis) IITK
EE250 (Control Systems Analysis) IITK
Lecture Notes
20 log10 K p
−20
−40 Kv
ω
√
ω
Ka
Figure 28.1: How to read K p , Kv , Ka off the BMPs.
28.4
Appendix: Various definitions of the Laplace transform
To what extent do we use the Laplace transform (LT) in our course? For us, the LT is only a tool, not the focus of our course. But, like any tool, we need to learn this one just enough to use it wisely. Firstly, the purpose of the LT in classical control theory is only mathematical rigor. Before the second world war, classical control used to be done using the Heaviside operator p instead of the LT []. Though it worked, the Heaviside operator was not supported by rigorous mathematics [26]. There are two kinds of LT we encountered in our earlier courses: one-sided LT and two-sided LT. For the function f (t), the various pairs of the definitions of the one-sided LT and the corresponding inverse LT are shown in Table 28.3. Here, c is any real number greater than the abscicca of convergence, α, of the LT. Given that s = σ + jω, the smallest value α of σ for which the LT exists is called the absicca of convergence of the LT [17]. In our course, we will not need to use the inverse LT. Instead, we will use tables of LT, where needed. The two-sided LT is as follows:
L{ f (t)} ,
Z ∞ −∞
f (t)e−st dt
But, we will not use this in our course because we will work with signals that are zero always but abruptly acquire non-zero value at some point in time, and this value of time can be arbitrarily chosen as t = 0. So, one-sided LT is appropriate in this case. Even when the signals may have non-zero values before t = 0, we assume that we know their effect upto either t = 0− (immediately before t = 0), or t = 0, or t = 0+ (immediately after t = 0) in the form of initial conditions of differential equations at these time instants. So, one-sided LT is appropriate in this case too. Distinction between the three definitions of one-sided LT is needed when an impulse is involved at the origin [25, page 62]. So, in the next section, we will recap our knowledge of the unit impulse function. January 2, 2012
149 of 218
Ramprasad Potluri
ω
EE250 (Control Systems Analysis) IITK
Lecture Notes
Table 28.3: Pairs of definitions of the one-sided Laplace Transform. Definition F− (s) = L− { f (t)} , 1 L− − { F− ( s )} ,
Z c+ j∞ 1
2πj
c− j∞
1 2πj
Z c+ j∞ c− j∞
1 L− + { F+ ( s )} ,
28.5
2πj
c− j∞
0−
Z ∞
[17] f (t)e−st dt
0
f (t)e−st dt
F (s)est ds = f (t) , t > 0
F+ (s) = L+ { f (t)} , Z c+ j∞ 1
Z ∞
F− (s)est ds = f (t) , t > 0−
F (s) = L{ f (t)} ,
L−1 { F (s)} ,
Source
Z ∞ 0+
[17] f (t)e−st dt
F+ (s)est ds = f (t) , t > 0+
Appendix: Unit impulse function
A unit impulse function (a.k.a Dirac delta function), denoted δ(t), has the following properties among others that are familiar to us and that we may use in this course: 1. δ(τ − t) = 0, ∀τ 6= t [27]. 2. δ(t) =
d1(t) dt ,
where 1(t) is a unit step function, a.k.a., Heaviside function, ( 1( t ) =
0 1
t<0 t>0
The Heaviside function has a jump discontinuity at t = 0. Its value at t = 0 is usually taken to be 1/2. Sometimes, its value is taken equal to some constant c, where 0 < c < 1 [27]. R +∞ R +∞ R 0+ R 0+ 3. −∞ δ(τ − t)dτ = −∞ δ(t − τ )dτ = 0− δ(t − τ )dτ = 0− δ(τ )dτ = 1. 4.
R +∞
−∞ δ ( τ − t ) f (τ ) dτ = f ( t ), if f ( τ ) is continuous at τ = t. “In other words, the impulse is so short and so intense that no value of f matters except over the short range where the δ occurs” [25, page 53].
We can work this property out as follows. The product f (τ )δ(τ − t) is illustrated in Figure 28.2. January 2, 2012
150 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
δ 0 t
τ
t
τ
f
f ×δ
0
0 t
τ
Figure 28.2: Graphical illustration of the product f (τ )δ(τ − t).
We see that this product is equal to zero everywhere, except at τ = t where it equals the unit impulse multiplied by f (t). Similarly, Z +∞ −∞
δ(t − τ ) f (τ )dτ =
Z +∞ −∞
δ(t − τ ) f (t)dτ = f (t)
Z +∞ −∞
δ(t − τ )dτ = f (t)
The unit impulse function is not a function in the classical sense. For a classical function r (t) that is defined as zero everywhere except at finitely many R +∞ points, the integral −∞ r (t)dt evaluates to zero. Functions such as δ are called generalized functions and are studied in the theory of distributions or measure theory (see, for example, [27]).
28.6
Appendix: Care needed in using the one-sided LT
“It is the transform of the unit-impulse function that led us to choose the L− transform rather than the L+ transform” [25, page 64]. Much confusion exists in literature about the definition of the one-sided LT and the delta function. The book [17] was probably the first to notice this confusion. The paper [28] presents a survey of this confusion, and also lists the literature that does not have this confusion. The book [25] is cited in [28] as treating the one-sided LT and the delta function in a consistent manner.
January 2, 2012
151 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
152 of 218
Lecture Notes
Ramprasad Potluri
Lecture 29
An example of application of loop-shaping in automobile control
153
EE250 (Control Systems Analysis) IITK
January 2, 2012
154 of 218
Lecture Notes
Ramprasad Potluri
Part V
Digital implementation of controllers
155
Lecture 30
State space realization of transfer function 30.1
Example of state space realization
Consider the following TF G (s) =
a s + a0 Y (s) = 2 1 U (s) s + b1 s + b0
We wish to find a state space model whose input and output are u(t) and y(t) respectively, and whose TF is G (s). The following is one way to accomplish this goal: 1. Introduce an intermediate variable X (s) thus: G (s) =
Y (s) Y (s) X (s) = U (s) X (s) U (s)
2. Let X (s) 1 = 2 U (s) s + b1 s + b0
and
Y (s) = a1 s + a0 X (s)
3. From the last step above, we can write the following: s2 X (s) = U (s) − b1 sX (s) − b0 X (s) Y (s) = a1 sX (s) + a0 X (s) 4. Using this last step, we can construct the simulation diagram of Figure 30.1. Note that here, the block containing the 1/s represents integration operation. Books show either of the two equivalent blocks shown in Figure 30.2. The 1/s is more appropriate R when the block diagram is in the Laplace Transform domain, while the is more appropriate when the block diagram is in time domain. 157
EE250 (Control Systems Analysis) IITK
Lecture Notes
a1 s2 X ( s )
U (s) +
−
1 s
sX (s)
X (s)
1 s
+ a0
Y (s)
+
− b1 b0
Figure 30.1: Simulation diagram for Step 4.
1 s
=
R
Figure 30.2: The two equivalent blocks that represent an integrator. The 1/s is more appropriate when the block diagram is in the Laplace Transform domain, R while the is more appropriate when the block diagram is in time domain.
a1 U (s) u(t) + −
s2 X ( s ) x˙ 2 −
1 s
sX (s) x2 = x˙ 1
1 s
X (s) x1
+ a0
+
Y (s) y(t)
b1 b0
Figure 30.3: Writing in the time-domain quantities into the simulation diagram of Figure 30.1.
January 2, 2012
158 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes a1
U (s) u(t) + −
s2 X ( s )
−
x˙ 1
sX (s)
1 s
x1 = x˙ 2
X (s)
1 s
+ a0
x2
Y (s)
+
y(t)
b1 b0
Figure 30.4: An alternative ssignment of state variables to the outputs of the integrators of the simulation diagram of Figure 30.1.
5. In the block diagram, write the time-domain quantities u(t) and y(t) as shown in Figure 30.3. Also, assign a state variable to the output of each integrator block. In the above case, we started the assignment at the output of the right-most integrator block, going left. 6. Write the differential equations based on the time-domain quantities in the simulation diagram x˙ 1 = x2 x˙ 2 = −b0 x1 − b1 x2 + u y = a0 x1 + a1 x2 This system of first order differential equations gives us the state-space (SS) model # " # " # " x1 0 0 1 u with x = x+ x˙ = 1 x2 −b0 −b1 Problem 26 Assign a state variable to the output of each integrator block starting from the left-most integrator, going right as shown in Figure 30.4. Verify that the SS model is as follows: " x˙ =
−b1 0
−b0 1
#
" x+
1 0
#
" u,
with x =
x1
#
x2
In Example 26 and the disucssion preceding it, we illustrated two separate SS realizations for a given TF. In practice, to each TF correspond an infinite number of SS realizations. The dimension of a SS model is the dimension of the space from which the state vector (also known as state) x comes. In our example, x comes from a 2-D space. The simulation diagram approach shown here is not the only simulation diagram approach for obtaining the SS realization. Also, a simulation diagram approach is not the only approach for obtaining a SS realization. January 2, 2012
159 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK Summator
Lecture Notes
Integral controller Power amplifier
Desired output
− +
− +
DC motor
− +
Signal inversion Tachogenerator or potentiometer
− +
Figure 30.5: An example of an op-amp circuit for controlling a DC motor.
The above diagram is called a simulation diagram because it can be built using hardware and simulates the behavior of the TF G (s) and of the corresponding differential equation. Simulation diagrams can be implemented for LTI systems using only operational amplifiers, capacitors, resistors, and, of course, power supplies. Indeed, every LTI state-space equation can be implemented using an op-amp circuit [29, page 16]. A simulation diagram represents an analog computer. Analog computers for LTI systems implement the operations of integration, summation, and amplification. Analog computers used to be used for simulation before the days of simulation on digital computers. An op-amp circuit can be built for controlling the speed or position of the shaft of a DC motor as shown in Figure 30.5. In this figure, each of the op-amp-based blocks has been built in a inverting feedback configuration. We need to be aware that this configuration introduces a phase shift of 180◦ . An alternative is to build the op-amp-based blocks using a non-inverting feedback configuration. The integral controller has the transfer 1 . The above figure does not show the power supply to each opfunction sCR amp. Why did we build the circuit using op-amps? For example, a capacitor too has the transfer function of an integrator. Then why not use a simple capacitor in place of an op-amp based implementation of an integrator? Note that, as mentioned in Lecture 1, each block in a block diagram needs to be non-loading. Only then are our manipulations with block diagrams, the way we have been doing them so far, valid. An op-amp has an almost infinite input impedance and an almost zero output impedance. The first means that an op-amp will not load the block that preceds it in a circuit. The second means that it keeps a minimal share of the signal (voltage) that it outputs. On the other hand, a January 2, 2012
160 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
simple capacitor, if used as replacement for an op-amp based integrator, would load the preceding block. Op-amps, in either inverting or non-inverting feedback configuration, can be used to build arbitrary rational TFs.
30.2
Numerical integration of state space model
A SS model has more than one use. We can design in the SS framework a control system for a given plant. In this sense, the SS modeling framework is an alternative to the TF modeling framework. Another use of the SS model is illustrated in the next section. In preparation for that section, let us see how we can numerically integrate the differential equation that a SS model represents. Suppose we have the following SS model: x˙ = Ax + Bu y = Cx + Du Example 21 Which of the following two definitions of a derivative do we use in our differential equations, and why? x˙ , lim
∆t→0
x (t + ∆t) − x (t) ∆t
or
x˙ , lim
∆t→0
x (t) − x (t − ∆t) ∆t
The first one is sometimes called right derivative in order to contrast with the left derivative which is the second one. The concept of a left derivative is used when talking about the differentiability of a function, where we say that for the function to be differentiable at a point the left derivative needs to be equal to the right derivative at that point. The left and right derivatives are thus mathematical artefacts that are used mainly to talk about differentiability of functions. For all other purposes, the derivative is defined only in the first manner, and has the following two well-known physical interpretations: (a) the slope of the function at each point, the slope being calculated by looking ahead from the point; (b) as estimating the value of the variable at the next time instant, given its value at the current instant and its rate of change at the current instant: y(t + ∆t) = y(t) + v(t)∆t. In fact, the second interpretation is actually the origin of the derivative operation [30, page 112].
The derivative as a forward-looking element is the meaning and purpose with which it is used in control systems. We will see evidence of this statement when we talk about derivative control. Digital implementation of a controller is necessarily constrained by the speed of the processor. So, ∆t cannot be made arbitrarily small. This means January 2, 2012
161 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
that in a digital implementation of the derivative action, using the left derivative may not give the same result as does using the right derivative. Now let us return to the SS model. Euler’s approximation of the derivative is x (t + ∆t) − x (t) x˙ ≈ , for small ∆t. ∆t Using Euler’s approximation, we can write the state space model as x (t + ∆t) − x (t) ≈ Ax (t) + Bu(t), ∆t y(t) = Cx (t) + Du(t)
for small ∆t.
which can be rewritten as x (t + ∆t) ≈ ( A∆t + I ) x (t) + B∆tu(t), y(t) = Cx (t) + Du(t) We denote by x (k) the state at the k-th instant. The k-th instant is the time k∆t starting from t = 0. Then the state at the instant (k + 1)∆t is x (k + 1). Using this notation, the SS model can be written as follows: x (k + 1) = ( A∆t + I ) x (k ) + B∆tu(k ), y(k) = Cx (k) + Du(k) In the literature, this system of equations is also written as xk+1 = ( A∆t + I ) xk + B∆tuk , yk = Cxk + Duk This system of equations is called discrete-time state-space model, and is convenient for numerical integration of the continuous-time version. The numerical integration can be made sufficiently accurate by reducing ∆t. In situations where there is a lower bound on ∆t, there are other methods of discretizing a continuous-time SS model that are more accurate than the one resulting from Euler’s approximation. See, for example, the chapter on ordinary differential equations in [31]. Apart from numerical integration of the continuous-time SS model, another use of the discrete-time version is that it can be used for the digital implementation of a controller. We design a controller, for example, using Bode plots. We have the TF of the controller. We convert this TF into a SS model. We discretize it and program it in a microcontroller/microprocessor/DSP. We see this process in more detail in EE380. In the next section, we see how a discrete-time version of a SS model can be used to numerically simulate the response of a control system to an arbitrary input.
January 2, 2012
162 of 218
Ramprasad Potluri
Part VI
Pole-placement design using root locus
163
Lecture 31
Unit step response of second order TF 31.1
State-space realization (concluded): Determining the response of a transfer function to an arbitrary input
MATLAB provides us the step function to plot the unit step response of an arbitrary TF. However, if we wish to plot the reponse of the same TF to an arbitrary input, then MATLAB does not yet have any ready-made functions. We have the following way out. 1. Convert the TF into a SS model. 2. Discretize this model, for example, using Euler’s approximation. 3. Numerically integrate this model using the chosen input. As an example that illustrates the above three steps, consider the control system we designed in Section 27.2. Example 22 The plant TF was G (s) =
10 s s +1 10
and the controller TF was s 2 s 10 +1 +1 D (s) = 2.5 2 10 s s +1 +1 0.7 40 We wish to see if the following two design specs are satisfied by our design: 1. Sinusoidal inputs of up to 1 rad/sec to be reproduced with ≤ 2% error. 165
EE250 (Control Systems Analysis) IITK
Lecture Notes
2. Sinusoidal inputs with a frequency of greater than 100 rad/sec to be attenuated at the output to ≤ 5% of their input value. The m-file simul.m — that uses the m-file cleanup.m (that we used earlier in Lecture 24) — implements the three steps that we listed above. We see from figures 31.1 and 31.2 that both the specs are satisfied.
Problem 27 Peform the minimum needed modifications in the m-file simul.m so that this file executes successfully in GNU Octave.
Problem 28 Which of the two methods of realization that we saw in this lecture does the MATLAB function tf2ss perform?
1.5 y u 1
0.5
0
−0.5
−1
−1.5
0
1
2
3
4
5
6
7
8
9
10
Figure 31.1: Results for Example 22: Response y of the control system designed in Lecture 27 to a sinusoidal input u with ω = 1 rad/s, T = 0.001 s, N = 10000.
January 2, 2012
166 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
1 y u 0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figure 31.2: Results for Example 22: Response y of the control system designed in Lecture 27 to a sinusoidal input u with ω = 100 rad/s, T = 0.001 s, N = 1000.
31.2
Unit step response of standard second order TF
31.3
Appendix: MATLAB codes
% simul.m: generates the response of the CL system % designed in Lecture 27 to a sinusoidal input. % This file runs in MATLAB. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all, close all s = tf(’s’); % Define G(s) and D(s) from Lecture 27: G = 10/s/(s/10+1); D = 10*(s/2.5+1)^2*(s/10+1)/(s/0.7+1)^2/(s/40+1); % Form the closed-loop (CL) TF: Gcl = D*G/(1+D*G); % Extract the numerator and denominator of this % CL TF for use with tf2ss: [num,den] = tfdata(Gcl,’v’); January 2, 2012
167 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
2.5
2
y(t)
1.5
1
0.5
0
tsa 0
0.5
1
tse 1.5
2
2.5
3
t (in s)
Figure 31.3: The unit step response of the standard second order TF. tsa and tse are respectively the actual settling time of the response and the settling time of the envelope of the response.
January 2, 2012
168 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
[n,d] = cleanup(num,den); % Use tf2ss to convert the TF model into a state % space (SS) model: [A,B,C,D] = tf2ss(n,d); % % % w
Define the step size T and the number of time steps N over which we wish to compute the sinusoidal response: = 1; T = 0.001; N = 10000;
% Define an identify matrix whose order is that of A. I = eye(size(A,1),size(A,2)); % Declare and initialize variables x and y: x(:,1) = zeros(size(A,1),1); y(:,1) = zeros(size(C,1),1); % Compute the sinusoidal response using the discret % version of the SS model (we discretize using % Euler’s approximation): for k = 1:N u(k) = sin(w*k*T); x(:,k+1) = (A*T + I)*x(:,k) + B*T*u(k); y(:,k) = C*x(:,k) + D*u(k); end % Plot y versus t and u versus t on one plot: t = (1:N)*T; plot(t,y,’-’,t,u,’--’); legend(’y’,’u’); grid print -depsc sineresp.eps % ord2sys.m: Plot the step response and envelopes of the 2nd % order system seen in Lecture 7. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all clear all zeta = 0.1; omegad = 20; tau = 1/(zeta*omegad); timelen = 6*tau; numsteps= 1000; t = linspace(0.001,timelen,numsteps); y = 1 - exp(-zeta*omegad*t).*sin(omegad*t+acos(zeta))/sqrt(1-zeta^2);
January 2, 2012
169 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
% Equations of the upper and lower envelopes: yeu = 1 + exp(-zeta*omegad*t)/sqrt(1-zeta^2); yel = 1 - exp(-zeta*omegad*t)/sqrt(1-zeta^2); % Specify the x% tube: x = 5; xl = (1-x/100)*ones(1,numsteps); xu = (1+x/100)*ones(1,numsteps); plot(t,y,t,yeu,t,yel,t,xl,’--’,t,xu,’--’), grid xlabel(’t’); ylabel(’y(t)’);
January 2, 2012
170 of 218
Ramprasad Potluri
Lecture 32
Second order TF with additional zero 32.1
Analysis of formulae developed for unit step response of standard second order system
Figure 32.1 shows the plots.
32.2
The case of an additonal zero in a second order TF
The m-file addzero.m plots the responses due to the additional zeros in the second-order transfer function. This file calculates the derivative dy/dt using the five point numerical differentiation method given by the expression1 f 0 (x) ≈ 1 See
− f ( x + 2h) + 8 f ( x + h) − 8 f ( x − h) + f ( x − 2h) 12h
wikipedia.org.
171
tr*omegan
EE250 (Control Systems Analysis) IITK 100 50 0 tp*omegan
Lecture Notes
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5 xi
0.6
0.7
0.8
0.9
1
100 50 0
Mp
1 0.5
ts*omegan
0 4000 2000
ess*omegan
0 2 1 0
Figure 32.1: Plot of the formulae developed for the unit step response of the standard second order system.
Note that addzero.m implements this formula without using a FOR loop. Problem 29 Figure 32.2 shows the results from executing addzero.m. Identify the various curves in this figure.
January 2, 2012
172 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
8
6
4
2
0
−2
−4 0
0.5
1
1.5
2
2.5
3
Figure 32.2: Results from executing addzero.m. Can you identify each curve?
January 2, 2012
173 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Problem 30 In the example discussed in Section 4.1, we see that even though the poles are deeper in the left half s-plane than the zero, the response is to our satisfaction. We intended to revisit this example in our lectures on root locus. Here, we revisit this example. Are there an contradications or errors in that section?
32.3
Appendix: MATLAB codes
% allinone.m: m-file to plot tr, tp, Mp, % ts, essr vs zeta. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% close all clear all zeta = linspace(0.001,0.999,1000); %Rise time: tromn = (pi-acos(zeta))./sqrt(1-zeta.^2); %tromn is the product of tr and omega_n. % Peak time: tpomn = pi./sqrt(1-zeta.^2); % Peak overshoot: Mp = exp(-pi*zeta./sqrt(1-zeta.^2)); % Settling time of envelope (x% tube): x = 2; tsomn = (1./zeta).*log(100./(x.*sqrt(1-zeta.^2))); % Exact expression for ess*omega_n|ramp: essomnr = 2*zeta; subplot(5,1,1), plot(zeta,tromn), ylabel(’t_r*omega_n’), grid January 2, 2012
174 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
subplot(5,1,2), plot(zeta,tpomn), ylabel(’t_p*omega_n’), grid subplot(5,1,3), plot(zeta,Mp), ylabel(’Mp’), grid subplot(5,1,4), plot(zeta,tsomn), ylabel(’t_s*omega_n|env’), grid, subplot(5,1,5), plot(zeta,essomnr), xlabel(’zeta’), ylabel(’ess*omega_n|ramp’), grid % addzero.m: plots the responses due to additional zeros in the % second-order transfer function. % The TF with the additional zero is the following: % % wn^2*(s+alpha) % Gaz = ---------------------------% alpha*(s^2+2*zeta*wn*s+wn^2) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all close all alpha = -5; zeta = 0.2; wn = 10; tau = 1/(zeta*wn); steps = 2500; t = linspace(0,6*tau,steps); % 6*tau because the transient will die out in about 4-5 tau. % Evaluate y(t) and its envelopes yel(t) & yeh(t): yel = 1 - exp(-zeta*wn*t)/sqrt(1-zeta^2); yeh = 1 + exp(-zeta*wn*t)/sqrt(1-zeta^2); ysin = sin(wn*sqrt(1-zeta^2).*t+acos(zeta)); y = 1 - exp(-zeta*wn*t)/sqrt(1-zeta^2).*ysin; % % % % %
Numerically differentiate y using the ‘‘five point method’’ method described at http://en.wikipedia.org/wiki/Numerical_differentiation I want to realize this formula using matrix manipulations alone, and not using a FOR loop. Here is how I am doing it:
yleftleft = [0,0,y(1,1:size(y,2)-2)]; yleft = [0,y(1,1:size(y,2)-1)]; yright = [y(1,2:size(y,2)),0]; yrightright = [y(1,3:size(y,2)),0,0];
January 2, 2012
175 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
h = t(1,2)-t(1,1); dybydt = (yleftleft - 8*yleft + 8*yright - yrightright)/(12*h); % % % % % % % % % %
The last point of dybydt comes out to be pretty large in comparison to the previous points. My guess is that the first point should also turn out to be large for the same reason that the last point is large (my guess turned out to be incorrect; don’t know why). As our dybydt was computed by prefixing two zeros and suffixing two zeros, the correct dybydt will be the one obtained by dropping the first 2 and the last 2 points. The slight gap at the origin resulting from this will not be so visible if we take sufficiently large number of steps. The following commands implements these thoughts:
tnet = t(1,3:size(t,2)-2); ynet = y(1,3:size(y,2)-2); yelnet = yel(1,3:size(yel,2)-2); yehnet = yeh(1,3:size(yeh,2)-2); dybydtnet = dybydt(1,3:size(dybydt,2)-2); yaznet = ynet+dybydtnet/alpha; % Specify the x% tube: x = 5; xl = (1-x/100)*ones(1,steps-4); xu = (1+x/100)*ones(1,steps-4); realpartofcomplexpole = -zeta*wn realzero = -alpha % Plot the y(t), yaz(t), envelopes of y(t), and x% tube: plot(tnet,ynet,tnet,dybydtnet,tnet,yaznet,tnet,yelnet,tnet,yehnet,tnet,xl,tnet,xu), grid legend(’ynet’,’dybydtnet’,’yaznet’,’yelnet’,’yehnet’,’xl’,’xu’) % % % % %
The default grid created by GNU Octave may not be dense enough for our purpose. We can have control over the location of the xticks and yticks as shown by the following code. If these commands do not work in MATLAB, then they can be commented out and analogous commands can be written for MATLAB.
% set(gca(),"xtick", (linspace(min(tnet),max(tnet),10))); % all = [ynet,yelnet,yehnet,yaznet]; % set(gca(),"ytick", linspace(min(all),max(all),10)); % How can we achieve a similar control over the xticks and yticks % in MATLAB? % Finally, I tried to increase the thickness of the lines using % the following command: % January 2, 2012
176 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK % % % % %
Lecture Notes
set(gca,"linewidth",1) This command is returning the message: warning: set: invalid property ‘linewidth’ Maybe "linewidth" has not yet been implemented in Octave.
print -depsc addzero.eps
January 2, 2012
177 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
178 of 218
Lecture Notes
Ramprasad Potluri
Lecture 33
Second order TF with additional pole; recipe for pole-placement 33.1
Additional pole in second order TF
179
EE250 (Control Systems Analysis) IITK
Lecture Notes
The m-file addpole.m helps us verify the assertions made in this section.
33.2
Recipe for placement of CL poles in a second order TF
January 2, 2012
180 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK +
G (s)
Lecture Notes
+
−
K
G (s)
−
K
Figure 33.1: Two block diagrams for which exist software support for the determination of root loci.
33.3
Root locus: What it is
For the CL systems shown in Fig 37.1, the root locus is the “set of values of s for which 1 + KG (s) = 0 is satisfied as the real parameter K varies from 0 to ∞” [11, page 295]. Alternatively, the root locus is the “graph of all possible roots of 1 + KG (s) = 0 relative to parameter K . . .” [11, page 290]. The MATLAB and GNU Octave function rlocus helps plot root loci. MATLAB’s rlocus function works with the left hand block diagram of Figure 37.1, while GNU Octave’s rlocus function works with the right hand block diagram. Example 23 Let G (s) = 1/(s(s + 2)). Determine the root locus.
33.4
Appendix: MATLAB codes
% addpole.m: Helps see the effect an additional pole % in the standard second order TF. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all close all wn = 10; alpha = 1; zeta = 0.2; Gap = tf(wn^2*alpha,conv([1,alpha],[1,2*zeta*wn,wn^2])); realpartofcomplexpole = -zeta*wn realpole = -alpha [gap,tap] = step(Gap); [g,t] = step(tf(wn^2,[1,2*zeta*wn,wn^2])); plot(tap,gap,’*’,t,g), grid on, legend(’gap’,’g’)
January 2, 2012
181 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
182 of 218
Lecture Notes
Ramprasad Potluri
Lecture 34
Root locus (Part 1) For G (s) =
∏im=1 (s + zi ) , where m ≤ n, ∏in=1 (s + pi )
the equation that describes the root locus is the characteristic equation 1+K
∏im=1 (s + zi ) =0 ∏in=1 (s + pi )
(34.1)
stated the magnitude criterion (MC) and phase criterion (PC). Note that all the s that satisfy the PC also satisfy the MC in the sense that for each s that satisfies the PC, we can determine a real non-negative K from the MC. On the other hand, not all points s that satisfy the MC satisfy the PC. G (s) shown above has m finite zeros and n − m zeros at infinity. This is because the TF approaches zero as s approaches infinity.
34.1
Rules of sketching root locus
Root Locus Rule 1 The RL has n branches. Root Locus Rule 2 The RL starts at K = 0 at the poles and ends at K = ∞ on the m zeros, and n − m of its branches end at ∞. Root Locus Rule 3 On the real axis, the RL exists only on those segments to the right of which the sum of the number of poles and the number of zeros is odd.
34.2
Appendix: An alternative way to see Rule 2
Chirag Gupta (Y9186) showed the following way of arriving at Rule 2. Rewrite (34.1) as ∏in=1 (s + pi ) = −K ∏im=1 (s + zi ) We see that when K = 0, the roots of (34.1) are at the open-loop poles. 183
EE250 (Control Systems Analysis) IITK Rewrite (34.1) as
Lecture Notes
1 ∏im=1 (s + zi ) =− n ( s + p ) K ∏ i =1 i
We see that when K = ∞, the roots of (34.1) are at the open-loop zeros. This way seems to be more convincing than simply saying that Rule 2 follows from (34.1).
January 2, 2012
184 of 218
Ramprasad Potluri
Lecture 35
Root locus (Part 2) Root Locus Rule 4 The RL is symmetric about the real axis. Root Locus Rule 5 As K → ∞, n − m branches of the RL go asymptotically along rays that emanate from the centroid s = −σ = −
∑in=1 pi −∑im=1 zi . n−m
Root Locus Rule 6 The angles that these rays form with the real axis are ±180◦ (2q+1) , q = 0, 1, 2, . . .. n−m
185
EE250 (Control Systems Analysis) IITK
January 2, 2012
186 of 218
Lecture Notes
Ramprasad Potluri
Lecture 36
Root locus (Part 3) 36.1
Roots of multiplicity greater than 1
Consider the root locus shown in Figure 36.1. The point where the root locus breaks away from the real axis is called the break away point, and where it breaks in to the real axis is called the break in point. While this is the meaning which much of the literature associates with these two terms, a root locus can also break away or break in at points not on the real axis as shown in Figure 36.2. Note that where the root locus breaks away, or where it breaks in, the characteristic equation 1 + KG (s) = 0 has roots of multiplicity greater than 1, for the appropriate value of K. Therefore, we are interested in finding such roots. Root Locus Rule 7 The roots of multiplicity greater than one occur on the RL where s satisfies dG/ds = 0, or, equivalently, dK/ds = 0. To determine the multiplicity of a root of the characteristic equation 1 + KG (s) = 0, we can utilize the following observation. Suppose s = −r is a root of 1 + KG (s) of multiplicity 2. Then, 1 + KG (s) can be written as (s + r )2 M0 (s), where M0 (s) does not contain the factor s + r. d{1+KG (s)} ds
1 = K dG ds ; dG/ds can be written as ( s + r ) M1 ( s ). Here, M1 ( s ) does not contain the factor s + r. d2 {1+KG (s)} ds2
2
= K ddsG2 ; d2 G/ds2 can be written as (s + r )0 M1 (s).
Thus, we see that both 1 + KG (s) and dG/ds evaluate to zero at s = −r, and that d2 G/ds2 does not evaluate to zero at s = −r. This reasoning shows that if s = −r is a root of multiplicity q of 1 + KG (s) d2 G for some value of K ∈ (0, ∞). Then, at s = −r, dG ds has q − 1 roots, ds2 has q − 2 ( q −1)
q
roots, · · · , d (q−1G) has 1 root, and ddsGq has no roots. ds We can use the above result to determine the roots of the characteristic equation of multiplicity greater than 1 as follows. Step 1 Determine the roots of dG/ds by setting dG/ds = 0. 187
EE250 (Control Systems Analysis) IITK
Lecture Notes
Root Locus 2
1.5 Break away point 1
Imaginary Axis
Break in point 0.5
0
−0.5
−1
−1.5
−2 −9
−8
−7
−6
−5
−4 Real Axis
−3
−2
−1
0
1
Figure 36.1: A root locus with break away and break in points on the real axis.
Step 2 Plug each of these roots into 1 + KG (s) = 0 and identify the root that results in a K that is real and non-negative. Suppose s = sˆ is such a root. Step 3 Evaluate di G/dsi at s = sˆ beginning i = 2 until we find an i for which this derivative does not go to zero at s = sˆ. Let this value of i be i.ˆ Step 4 The characteristic equation has iˆ roots at s = sˆ. An alternative to the algorithm shown in the above steps is to work with dK/ds instead of with dG/ds. The justification for this alternative is as follows. 0 0 Suppose G (s) is of the form N (s)/D (s). Then dG/ds = N DD−2ND , where N 0 = dN/ds and D 0 = dD/ds. Step 1 above is equivalent to solving N 0 D − ND 0 = 0. On the other hand, the characteristic equation can be written as K = −1/G = − D/N. Setting dK/ds = 0 gives N 0 D − ND 0 = 0. Thus, both dG/ds = 0 and dK/ds = 0 are equivalent to N 0 D − ND 0 = 0. Problem 31 1. In the above algorithm, can we also use di K/dsi instead of di G/dsi , for i = ˆ 2, . . . , i? 2. In the first algorithm, in working with dG/ds and higher derivates, our assumption is that K is independent of s. On the other hand, in the alternative algorithm, we are working with dK/ds, thus seemingly contradicting the assumption. Explain. January 2, 2012
188 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
36.2
Lecture Notes
Angles of departure and angles of arrival
Root Locus Rule 8 To determine the angle of departure from a pole, use the PC along with a point s = sˆ that is infinitesimally close to the pole the angle of departure from which needs to be calculated. To determine the angle of arrival at a zero, use the PC along with a point s = sˆ that is infinitesimally close to the zero the angle of arrival at which needs to be calculated. We can use the idea presented in Example 24 to determine the angle of departure of the RL of 1 + s(sK+2) = 0 from the two roots that are at s = −1. Example 24 Demonstrate that the locus of the roots of the equation 1 + s(sK+2) = 0 beginning from the roots s = −1 is the same as the locus of the roots of the equation 1 + (s+K1)2 = 0 for K ∈ [0, ∞). Consider the characteristic equation 1 + the equation
s2
K s ( s +2)
= 0, which is equivalent to
+ 2s + K = 0, which is equivalent to s2 + 2s + 1 + K − 1 = 0,
Root Locus of 1/(s(s+2)(s2+2s+5)) 8
6 Break away points
Imaginary Axis
4
2
0
−2
−4
−6
−8 −8
−6
−4
−2
0
2
4
6
Real Axis
Figure 36.2: A root locus with break away points not only on the real axis.
January 2, 2012
189 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK which is equivalent to 1 +
K −1 ( s +1)2
Lecture Notes
= 0.
= 0 and 1 + (sK+−11)2 = 0 have the same RL — even when K is varied only on a subset of [0, ∞). ˆ Therefore, the RL of 1 + s(sK+2) = 0, K ∈ [1, ∞), and the RL of 1 + (s+K1)2 = 0, Kˆ (= K − 1) ∈ [0, ∞), are the same. The above method can be applied to any other characteristic equation too. ∏m (s+z ) If 1 + K ∏ni=1(s+ pi ) = 0 has the roots pˆ 1 , pˆ 2 , . . . , pˆ n for some value K0 of K, then Therefore, both the equations 1 +
i =1
K s ( s +2)
i
∏m (s+z )
∏m (s+z )
the RL of 1 + K ∏ni=1(s+ pi ) = 0, K ∈ [K0 , ∞), and the RL of 1 + Kˆ ∏ni=1(s+ pˆi ) = 0, i i i =1 i =1 where Kˆ (= K − K0 ) ∈ [0, ∞), are the same. Using this idea, we see that the angle of departure of the root locus from the root s = −1 can be calculated from the phase criterion as
−2/sˆ + 1 = ±180◦ (2q + 1) where, s = sˆ is a point that is infinitesimally close to s = −1. We can set q = 0 as this RL has only two branches that emanate from s = −1, and ±180◦ is enough to calculate them. Thus, /sˆ + 1 = ±90◦ are the angles of departure of the RL from s = −1. This conclusion leads to the fact that a RL that departs from s = −1 perpendicular to the real axis is the only valid RL of 1 + s(sK+2) = 0. This completes the example started in Section 34.1.
January 2, 2012
190 of 218
Ramprasad Potluri
Lecture 37
Root contours Consider, for example, the equation s3 + s2 + βs + α = 0
(37.1)
We wish to study how the roots of this equation behave with change in α ∈ (0, ∞) and β ∈ (0, ∞). We could write (37.1) as 1+
βs =0 s3 + s2 + α
(37.2)
and determine the roots of s3 + s2 + α for various values of α, thus generating the locus of the roots of s3 + s2 + α, as shown in the upper subfigure of Figure 37.1. Then, keeping α fixed at a certain value, we could generate the root locus of (37.2). The upper subfigure of Figure 37.2 shows a collection of root loci of (37.2) generated in this manner. Such a collection is called root contour of (37.2). We could write (37.1) also as 1+
s3
α =0 + s2 + βs
(37.3)
and determine the roots of s3 + s2 + βs for various values of β, thus generating the locus of the roots of s3 + s2 + βs, as shown in the lower subfigure of Figure 37.1. Then, keeping β fixed at a certain value, we could generate the root locus of (37.3). The lower subfigure of Figure 37.2 shows a collection of root loci of (37.3) generated in this manner. This collection is the root contour of (37.3). Problem 32 What are some possible applications of root contours?
Problem 33 See Figure 37.2. It seems that the root contour of the upper subfigure shows that even though the roots of (37.1) are unstable for large values of α, they are stabilized by a corresponding increment in the value of β. On the other hand, 191
EE250 (Control Systems Analysis) IITK
Lecture Notes
rlocus(tf(1,[1,1,0,0])) 4
3
Imaginary Axis
2
1
0
−1
−2
−3
−4 −5
−4
−3
−2
−1
0
1
2
0
0.2
Real Axis rlocus(tf([1,0],[1,1,0,0])) 0.8
0.6
Imaginary Axis
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8 −1.2
−1
−0.8
−0.6 −0.4 Real Axis
−0.2
Figure 37.1: The locus of the roots of s3 + s2 + α (upper subfigure), and the locus of the roots of s3 + s2 + βs (lower subfigure).
January 2, 2012
192 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
rlocus(tf(1,[1,1,0,0])); hold; for alpha = 1:11, rlocus(tf([1,0],[1,1,0,alpha])); end; 15
10
Imaginary Axis
5
0
−5
−10
−15 −5
−4
−3
−2
−1
0
1
2
Real Axis rlocus(tf([1,0],[1,1,0,0])); hold; for beta = 0.1:21, rlocus(tf(1,[1,1,beta,0])); end; 10 8 6
Imaginary Axis
4 2 0 −2 −4 −6 −8 −10 −10
−8
−6
−4
−2
0 Real Axis
2
4
6
8
10
Figure 37.2: The root contour (upper subfigure) of (37.1) rewritten as (37.2), and the root contour (lower subfigure) of (37.1) rewritten as (37.3).
January 2, 2012
193 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
the root locus of the lower subfigure seems to show that even for large values of β, a slight increment in the value of α could make the complex conjugate pole pair unstable. There probably is some error in this observation as both the root contours eventually represent the same equation (37.1). Find this error.
January 2, 2012
194 of 218
Ramprasad Potluri
Lecture 38
Design using root locus Consider the control system of Figure 38.1. Discuss ways to improve its speed of response while maintaining its overshoot at a value determined by the damping ratio of ζ = 0.5.
38.1
Appendix: TA09
Solve the problems of TA06 using root locus based design. Demonstration of the use of MATLAB’s SISOTOOL for solving problem 1 of TA06 using root locus.
Yref (s) +
−
K
U (s)
1 s(s + 1)(s/5 + 1)
Y (s)
Figure 38.1: A control system whose speed of response we wish to improve while maintaining the damping ratio at ζ = 0.5.
195
EE250 (Control Systems Analysis) IITK
January 2, 2012
196 of 218
Lecture Notes
Ramprasad Potluri
Part VII
Addressing nonlinearities
197
Lecture 39
Stability of nonlinear control systems: Circle criterion Assume that G (s) in Figure 39.1 has no poles in the RHP and that the nonlinearity f is sector-bounded as defined in Figure 39.2. A sufficient condition for input-output (from r1 , r2 to e1 , e2 , y1 , y2 ) stability of the system in Figure 39.1 is that the polar plot (i.e., G ( jω ), 0 ≤ ω ≤ ∞) does not enter or enclose the circle that lies in the interval [−1/k1 , −1/k2 ] with center on the real axis in the Nyquist plane (see Figure 39.3). When the loop has nonlinearities in the forward path, instead of in the feedback path as in Figure 39.1, we can use the block diagram transformations shown in Figure 39.4 and then apply the circle criterion.
r1
e1
G (s)
y1
− y2
f (·)
e2
r2
Figure 39.1: Block diagram for the circle criterion[32].
199
EE250 (Control Systems Analysis) IITK
Lecture Notes
k2 y
f (y) k1 y y
0
Figure 39.2: A sector-bounded nonlinearity[32]. For such a nonlinearity, the following holds: f (0) = 0, k1 ≤ f (y)/y ≤ k2 for y 6= 0.
− k11
− k12
Im 0
Re G ( jω )
Figure 39.3: Illustration of the circle criterion. The circle’s center is on the real axis and radius equals 1/k1 −2 1/k2 .
January 2, 2012
200 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Figure 39.4: Block diagram transformations for the circle criterion[32].
January 2, 2012
201 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
202 of 218
Lecture Notes
Ramprasad Potluri
Lecture 40
Dealing with plant input constraints: An anti-windup scheme The problem we are interested in this lecture is as follows. In designing controllers in this course, we have assumed that the plant behaves more or less linearly. However, all real-world plants are nonlinear to lesser or greater extent. In some cases, the performance we see with our linear design and simulation may not be realizable when we deploy the controller on the actual plant. In such cases, we want to have techniques that help us design controllers for nonlinear plants. Such techniques are the subject of a course on nonlinear control. In our course, we adopt the simpler path of designing controllers assuming the plants are linear, and then modifying the control architecture to do the best we can with the linear controller. To simplify the discussion, we focus our attention on only a saturation nonlinearity. Figure 40.1 shows one possible solution. C (s) is a biproper transfer function. While [22] recommends selecting c∞ such that C (s) = c∞ + C (s), where C (s) is a strictly proper transfer function, the results of the simulink file test1.mdl show seem to suggest that c∞ can be any constant.
r (t) +
−
e(t) +
−
c∞
uˆ (t)
u(t) Lim
Lim
G (s)
y(t)
[C (s)]−1 − [c∞ ]−1
Figure 40.1: This anti-wind-up scheme is taken from [22]. The controller C (s) needs to be a biproper transfer function.
203
EE250 (Control Systems Analysis) IITK
January 2, 2012
204 of 218
Lecture Notes
Ramprasad Potluri
Part VIII
PID Control
205
Lecture 41
PID Controllers
207
EE250 (Control Systems Analysis) IITK
January 2, 2012
208 of 218
Lecture Notes
Ramprasad Potluri
Lecture 42
Ziegler-Nichols tuning of PID controllers 42.1
What is controller tuning?
Figure 42.1 shows PID control of a plant. “If a mathematical model of the plant can be derived, then it is possible to apply various design techniques for determining parameters of the controller that will meet the transient and steady-state specifications of the closed-loop (CL) system. However, if the plant is so complicated that its mathematical model cannot be easily obtained, then an analytical approach to the design of a PID controller is not possible. Then we must resort to experimental approaches the tuning of PID controllers.”[33]. “The process of selecting the controller parameters to meet given performance specifications is known as controller tuning.”[33]. In practice, tuning is used in the following three situations: 1. We have designed a controller and analyzed its performance on paper and using simulation tools such as Matlab. We wish to deploy this controller on the actual plant. We find that this controller does not give us the kind of behavior in practice that it gave in simulation. The design has only brought us roughly near what we want. Tuning of the parameters of the controller is needed to obtain the behavior that we saw in simulation. 2. The plant model may not be known; so we may not have designed a controller. We may use on-line estimation and tuning to stabilize the CL system. 3. We know that the plant belongs to a class of plants, and that for these plants a certain kind of controller will work. So, we bypass the design stage and go straight to tuning for a plant in this class. This is the situation that is addressed by the two Ziegler-Nichols tuning (ZNT) techniques.
42.2
What the two ZNT methods do
The design specifications for both the ZNT methods is quarter amplitude decay (QAD) [34]. What this means is that the tuned controller will impart the CL 209
EE250 (Control Systems Analysis) IITK r +
−
T s 1 u + D kp 1 + Plant TI s τs + 1
Lecture Notes y
Figure 42.1: The CL system. The ZNT methods help tune k p , TI , TD . The time constant τ has two purposes. Firstly, it makes implementing the derivative controller possible. Secondly, it helps create a low pass filter for high frequency noise. τ is selected to be between 0.1TD and 0.2TD [22, page 161].
unit step response a second overshoot whose ratio to the first overshoot is 25% [11, page 240, Figure 4.14]1 . The ZNT rules “suggest a set of k p , TI , and TD that will give a stable operation of the system. However, the resulting system may exhibit a large maximum overshoot in the step response, which is unacceptable. In such a case we need a series of fine tunings until an acceptable result is obtained. In fact, the ZNT tuning rules give an educated guess for the parameter values and provide a starting point for fine tuning, rather than giving the final settings for k p , TI and TD .” [33]. The first method is a time-domain method, while the second is a frequency domain method[34]. In the remainder of this section, we will describe these methods based on [33]. Though [22, page 162] only says that the second method is valid only for open-loop stable plants, we can see in the following that, as an “S” curve can be obtained for only such plants, the first method too is valid only for such plants.
42.3
First method
This method applies to plants whose step response looks approximately like an “S”. Step 1: Obtain the unit step response of the open-loop (OL) plant as shown in Figure 42.2. Step 2: Determine T and L as shown, by drawing a tangent as shown to the step response through the point of inflexion. Alternatively draw a tangent to this curve of maximum slope. Step 3: The approximate transfer function (TF) of the OL plant is ω (s) Ke− Ls = U (s) Ts + 1 Step 4: Tune the PID controller parameters according to Table 42.1. 1 Interestingly,
the fifth edition of this book does not have a chapter on ZNT.
January 2, 2012
210 of 218
Ramprasad Potluri
Unit step response of OL system
EE250 (Control Systems Analysis) IITK
Lecture Notes
K
Point of inflexion
t
0 L
T
Figure 42.2: Unit step response of the plant is needed for the first ZNT method. This response is also known as the reaction curve of the plant. As the first ZNT method uses these curves, the first method is also known as reaction curve method. See [22, page 167] for an alternative way of obtaining the reaction curve.
Table 42.1: ZNT rules for the First ZNT method [33, 34]. Type of controller P PI PID
January 2, 2012
kp
Ti
TD
T LK 0.9T LK 1.2T LK
∞
0
L 0.3
0
2L
0.5L
211 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Table 42.2: ZNT rules for the Second ZNT method [33, 34].
42.4
Type of controller
kp
Ti
TD
P PI PID
0.5k cr 0.45k cr 0.6k cr
∞ (1/1.2) Pcr 0.5Pcr
0 0 0.125Pcr
Second method
This method is also called the “ultimate gain method”. This method applies to plants that exhibit sustained oscillations in CL under proportional control for some value of k p > 0. Step 1: Form the CL system with k p > 0, TI = ∞, and TD = 0. Step 2: Apply a step input to the CL system and observe its response. Step 3: With the step input on, increase the value of k p from 0 to some value k cr (called the “critical gain” or “ultimate gain”) at which the CL system will exhibit sustained oscillations. Step 4: Determine the period Pcr of these oscillations. Step 5: Tune the PID controller parameters according to Table 42.2. Example 25 Take the plant transfer function as 1/(s + 1)3 , but assuming that you do not know what it is, use the ultimate gain method to tune the three controllers. The results are shown in Figure 42.3
Problem 34 Consider the first of the two simulink files that I gave you, goodwin_ex_ultimate_gain.mdl and goodwin_ex_reaction_curve.mdl. Theoretically determine the value of the critical gain, the value of the critical period, and the values of the closed-loop poles corresponding to this critical period.
42.5
Appendix: Further reading
Chapter 6 of [22] and the works by Karl Astrom are recommended for further reading on PID control, including tuning methods.
January 2, 2012
212 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
Output from loop to identify Kc 2 1 0
0
2
4
6
8
10
12
14
16
18
20
14
16
18
20
14
16
18
20
14
16
18
20
Output from loop with P−controller 2 1 0
0
2
4
6
8
10
12
Output from loop with PI−controller 2 1 0
0
2
4
6
8
10
12
Output from loop with PID controller 2 1 0
0
2
4
6
8
10 Time
12
Figure 42.3: Results from tuning using the ultimate gain method.
January 2, 2012
213 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
214 of 218
Lecture Notes
Ramprasad Potluri
Bibliography [1] Lalit mestha [people in control]. Control Systems Magazine, IEEE, 30(6):38 –41, December 2010. [2] Eric Paljug and Xiaoping Yun. Experimental study of two robot arms manipulating large objects. IEEE Transactions on Control Systems Technology, 3(2):177–188, June 1995. [3] Dong Sun, James K. Mills, and Yunhui Liu. Position control of robot manipulators manipulating a flexible payload. The International Journal of Robotics Research, 18(3):319 – 332, March 1999. [4] H. Osumi and T. Arai. Cooperative control between two positioncontrolled manipulators. In Proceedings of the IEEE International Conference on Robotics and Automation, volume 2, pages 1509–1514, 1994. [5] M. Uchiyama, N. Lwasawa, and K. Hakomori. Hybrid position/force control for coordination of a two-arm robot. In Proceedings of IEEE International Conference on Robotics and Automation, volume 4, pages 1242– 1247, 1987. [6] Tamer Basar. Preface. In Tamer Basar, editor, Control Theory: Twenty-Five Seminal Papers (1932-1981). Wiley-IEEE Press, 2000. This chapter is available at http://as.wiley.com/WileyCDA/WileyTitle/productCd-0780360214.html. [7] Christopher Bissell. A history of automatic control. In Shimon Y. Nof, editor, Springer handbook of automation, Springer handbook series (LXXVI). Springer Verlag, Heidelberg, Germany, 2009. This chapter is available on the web through a search using the author’s name and the title of the book. [8] Harold S. Black. Inventing the negative feedback amplifier. IEEE Spectrum, pages 55 – 60, December 1977. [9] David A. Mindell. Opening black’s box: Rethinking feedback’s myth of origin. Technology and Culture, 41(3):405 – 434, July 2000. http://www.jstor.org/stable/25147536. [10] Katsuhiko Ogata. Modern Control Engineering. Prentice-Hall, Inc., 2nd edition, 1990. Later editions available in India. [11] Gene F. Franklin, J. David Powell, and Abbas Emami-Naeini. Feedback Control of Dynamic Systems. Pearson Education, fourth edition, 2002. 215
EE250 (Control Systems Analysis) IITK
Lecture Notes
[12] William G. Vogt and Mohamed I. Younis. Straight-line approximations for phase-angle plots. IEEE Tranasactions on Education, pages 249 – 252, December 1968. [13] Jr. Thomas F. Schubert. A quantitative comparison of three bode straightline phase approximations for second-order, underdamped systems. IEEE Tranasactions on Education, 40(2):135 – 138, May 1997. [14] Rodger E. Ziemer, William H. Tranter, and D. R. Fannin. Signals and Systems: Continuous and Discrete. Prentice Hall, fourth edition, 1998. [15] Issac M. Horowitz. Synthesis of Feedback Systems. Academic Press Inc., 1963. [16] Madan Gopal. Control Systems — Principles and Design. Tata McGraw-Hill, 3rd edition, 2008. [17] Thomas Kailath. Linear Systems. Prentice-Hall, Inc., 1980. [18] Alan V. Oppenheim and Alan S. Willsky. Signals and Systems. Prentice Hall of India, 2nd edition, 1997. [19] Hubert M. James, Nathaniel B. Nichols, and Ralph S. Phillips, editors. Theory of servomechanisms. McGraw-Hill Book Company, 1947. [20] L.A. MacColl. Fundamental theory of servomechanisms. Princeton, NJ, USA, 1945.
Van Nostrand,
[21] K.J. Astrom. Regeneration Theory: H. Nyquist. In Tamer Basar, editor, Control Theory: Twenty-Five Seminal Papers (1932-1981). Wiley-IEEE Press, 2000. [22] Graham C. Goodwin, Stefan F. Graebe, and Mario E. Salgado. Control System Design. Prentice-Hall, 2001. [23] Benjamin C. Kuo. Automatic Control Systems. Prentice-Hall, 7th edition, 1995. [24] Farid Golnaraghi and Benjamin C. Kuo. Automatic Control System. John Wiley & Sons, Inc., 2010. [25] Gene F. Franklin, J. David Powell, and Abbas Emami-Naeini. Feedback Control of Dynamic Systems. Pearson Education, fifth edition, 2006. [26] Paul J. Nahin. Behind the Laplace transform. IEEE Spectrum, pages 60 – 60, March 1991. [27] Ram P. Kanwal. Generalized functions: theory and technique. Academic press, Inc. (London) Ltd., 1983. [28] Kent H. Lundberg, Haynes R. Miller, and David L. Trumper. Initial conditions, generalized functions, and the Laplace transform: Troubles at the origin. IEEE Control Systems Magazine, pages 22 – 35, February 2007. [29] Chi-Tsong Chen. Linear system theory and design. Oxford University Press, Inc., 3rd edition, 1999. January 2, 2012
216 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
Lecture Notes
[30] Feynman, Leighton, and Sands. The Feynmann Lectures on Physics. Narosa, c India, 2003. Copyright 1963 by Addison-Wesley Publishing Company, Inc. [31] Cleve Moler. Numerical Computing with MATLAB. Society for Industrial and Applied Mathematics, 2004. The electronic version of this book can be downloaded from http://www.mathworks.com/moler/chapters.html. [32] Torkel Glad and Lennart Ljung. Control Theory: Multivariable and Nonlinear methods. Taylor & Francis, 2000. [33] Katsuhiko Ogata. Modern Control Engineering. Pearson Education, fourth edition, 2002. Available in India. [34] Tore Hägglund and Karl J. Åström. Automatic tuning of PID controllers. In William S. Levine, editor, The Control Handbook, pages 817 – 826. CRC Press LLC, 1996.
January 2, 2012
217 of 218
Ramprasad Potluri
EE250 (Control Systems Analysis) IITK
January 2, 2012
218 of 218
Lecture Notes
Ramprasad Potluri