REAL TIME CONTROL FOR TIME DELAY SYSTEM Rihem Farkh
[email protected]Kaouther Laabidi
[email protected]Mekki Ksouri
[email protected]Research unit on Systems Analysis and Control National School Engineering of Tunis BP, 37, le Belvedere, 1002 Tunis, Tunisia KEYWORDS all values of the stabilization gains (Kp ,Ki ,Kd ) of the delay system, PI controller, stability region, Hermit- regulator for the case of first order delay system. The Biehler theorem, optimization, genetic algorithm. same technique is used to provide a complete characteri- zation of all P and PI controller that stabilize a first order delay system which considered as less complicated than ABSTRACT the PID stabilization problem (Silva et al, 2000). In order In this paper, we consider the control of time delay sys- to have good performances in closed loop, it is necessary tem by Proportional-Integral (PI) controller. By Us- to suitably choose the parameters of the regulator in the ing the Hermite-Biehler theorem, which is applicable to zone stability. quasi-polynomials, we seek the stability region of the In this work, we look for optimum regulators under controller and the computation of its optimum parame- different criteria and we present the results of this ap- ters. We have used the genetic algorithms to lead the proach when applied to the temperature control of a complexity of the optimization problem. An applica- heater air stream P T − 326. tion of the suggested approach for real-time control on This paper is structured as follows: in section 2, we a heater P T − 326 was also considered. present the theorem of Hermit-Biehler applicable to the quasi-polynomials. Section 3 is devoted to the problem formulation for first order delay system controlled via INTRODUCTION PI controller. In order to obtain optimal regulator in the Systems with delays represent a class within infinite size zone of stability, a description of the genetic algorithms largely used for the modelling and the analysis of trans- is presented in section 4. Section 5 is reserved for real port and propagation phenomena (matter, energy or in- time example to test the control method. formation)(Niculescu, 2001; Zhong, 2006). They nat- urally appear in the modelling of processes found in PRELIMINARY RESULTS FOR ANALYZING physics, mechanics, biology, physiology, economy, dy- TIME DELAY SYSTEM namics of the populations, chemistry, aeronautics and aerospace. In addition, even if the process itself does Several problems in process control engineering are re- not contain delay, the sensors, the actuators and the com- lated to the presence of delays. These delays intervene putational time implied in the development of its control in dynamic models whose characteristic equations are of law can generate considerable delays (Niculescu, 2001) . the following form (Silva et al, 2002, 2005): The latter have a considerable influence on the behaviour of closed-loop system and can generate oscillations and δ(s) = d(s) + e−L1 s n1 (s) + e−L2 s n2 (s) even instability (Dambrine, 1994). (1) PID controllers are of high interest thanks to their +... + e−Lm s nm (s) broad use in industrial circles (Cedric, 2002). Tra- Where: d(s) and n(s) are polynomials with real coeffi- ditional methods of PID parameter tuning are usually cients and Li represent time delays. These characteristic used in the case of the systems without delays (Lequin equations are recognized as quasi-polynomials. Under et al, 2003; Liu et al, 2001). An analytical approach the following assumptions: was developed in (Silva et al, 2005; Zhong, 2006) (A1 ) deg(d(s)) = n and deg(ni (s)) < n and allowed the characterization of the stability region f or i = 1, 2, ..., m (2) of delayed systems controlled via PID. Indeed, by us- (A2 ) L1 < L2 < ... < Lm ing the Hermit-Biehler theorem applicable to the quasi- polynomials (Bhattacharyya, 1995; Silva et al, 2005), the One can consider the quasi-polynomials δ ∗ (s) described authors have developed an analytical characterization of by : Proceedings 23rd European Conference on Modelling and Simulation ©ECMS Javier Otamendi, Andrzej Bargiela, José Luis Montes, Luis Miguel Doncel Pedrera (Editors) ISBN: 978-0-9553018-8-9 / ISBN: 978-0-9553018-9-6 (CD) yc (t) e(t) y(t) δ ∗ (s) = esLm δ(s) u(t) C(s) G(s) = esLm d(s) + es(Lm −L1 ) n1 (s) (3) +es(Lm −L2 ) n2 (s) + ... + nm (s) The zeros of δ ( s) are identical to those of δ ∗ (s) since esLm does not have any finite zeros in the complex plan. Figure 1: Closed-loop control of a time delay system However, the quasi-polynomial δ ∗ (s) has a principal term since the coefficient of the term containing the high- est powers of s and es is nonzero. If δ ∗ (s) does not have The PI Controller is described by the following trans- a principal term, then it has an infinity roots with positive fer function: real parts (Bhattacharyya, 1995; Silva et al, 2005). Ki C(s) = Kp + (6) The stability of the system with characteristic equa- s tion (1) is equivalent to the condition that all the zeros of Our objective is to analytically determine the region in δ ∗ (s) must be in the open left half of the complex plan. the (Kp , Ki ) parameter space for which the closed-loop We said that δ ∗ (s) is Hurwitz or is stable. The following system is stable. theorem gives a necessary and sufficient condition for the theorem 3 stability of δ ∗ (s) (Silva et al, 2000, 2001, 2002, 2005). The range of Kp value, for which a solution to PI stabi- theorem 1 Let δ ∗ (s) be given by (3), and write: lization problem for a given stable open-loop plant exists, is given by (Silva et al, 2000, 2001, 2002, 2005): δ ∗ (jω) = δr (ω) + jδi (ω) (4) r 1 T L2 where δr (ω) and δi (ω) represent respectively the real − < Kp < α21 + (7) and imaginary parts of δ ∗ (jω) . Under conditions (A1 ) K KL T2 and (A2 ), δ ∗ (s) is stable if and only if: T Where α1 the solution of the equation tan(α) = − L α 1: δr (ω) and δi (ω) have only simple, real roots and π in the interval [ 2 , π] these interlace, ′ ′ 2: δi (ω0 )δr (ω0 ) − δi (ω0 )δr (ω0 ) > 0 for some w0 in Proof The closed loop characteristic equation of the [−∞, +∞] system is given by: ′ ′ where δr (ω) and δi (ω) denote the first derivative with respect to w0 of δr (ω) and δi (ω), respectively. δ(s) = (KKi + KKp s)e−Ls + (1 + T s)s (8) A crucial stage in the application of the precedent the- orem is to verify that and have only real roots. Such a we deduce the quasi-polynomials: property can be checked while using the following theo- rem (Silva et al, 2000, 2001, 2002, 2005). δ ∗ (s) = eLs δ(s) = (KKi +KKp s)+(1+T s)seLs (9) theorem 2 Let M and N designate the highest powers of s and es which appear in δ ∗ (s) . Let η be an appro- substituting s = jw , we have: priate constant such that the coefficient of terms of high- δ ∗ (jω) = δr (ω) + jδi (ω) est degree in δr (ω) and δi (ω) do not vanish at ω = η . where: Then a necessary and sufficient condition that δr (ω) and δr (ω) = KKi − ω sin(Lω) − T w2 cos(Lω) δi (ω) have only real roots is that in each of the intervals δi (ω) = w [KKp + cos(Lω) − T w sin(Lω)] −2lπ + η < ω < 2lπ + η, l = l0 , l0 + 1, l0 + 2... (10) δr (ω) or δi (ω) have exactly 4lN + M real roots for a Clearly, the parameter Ki only affects the real part of sufficiently large l0 . δ ∗ (jω) whereas the parameter Kp affects the imaginary part. According to the first condition of theorem 1, we PI CONTROL FOR FIRST ORDER DELAY SYS- should check that the roots of δi and δr are simple. By TEM using the theorem 2, while choosing , M = 2, N = l = 1 and η = π4 , we observe that δi (ω) has simple roots for We consider the functional diagram of figure 1, in which any Kp checking (7) (Bhattacharyya, 1995; Silva et al, the transfer function of delayed system is given by (5) 2005). K the application of the second condition of theorem 2, led G(s) = e−Ls (5) us to: 1 + Ts ′ ′ E(ω0 ) = δi (ω0 )δr (ω0 ) − δi (ω0 )δr (ω0 ) > 0 Where K , T and L represent respectively the state gain, for ω0 = 0 (a value that annul δi (ω)) we obtain: KK +1 the constant time and the time delay of the plant. These E(ω0 ) = ( Lp )KKi > 0 which implies Kp > −a K 0 three parameters are supposed to be positive. since K > 0 and Ki > 0. This proofs the first inequality given by(2)in Theorem 1. which are applied to all the elements of the population We consider that z = Lω , we get: (Godelberg, 1991). Couples of parents are selected ac- cording to their functions costs. The crossover operator z T is applied with a Pc probability and generates couples of δr (z) = K Ki − (sin(z) + z cos(z)) (11) KL L children. The mutation operator is applied to the children z T with a Pm probability and generates mutant individuals It results that a(z) = KL (sin(z) + L z cos(z)) then who will be inserted in the new population. The reach- δr (z) = K [Ki − a(z)] (12) ing of a maximum number of generations is the criterion of stop for our algorithm. Figure 2 shows the basic flow for z0 = 0, we obtain : chart of AGs. The principle of regulator parameters opti- δr (z0 ) = K(Ki − a(0)) = KKi > 0 (13) generation of a random population for zj 6= z0 , j = 1, 2, 3..., we obtain: cost function δr (zj ) = K(Ki − a(zj )) (14) computation Interlacing the roots of δr (z) and δi (z) is equivalent to stop condition end δr (z0 ) > 0 (since Ki > 0), δr (z1 ) < 0, δr (z2 ) > 0... We can use the interlacing property and the fact that δi (z) has selection only real roots to establish that δr (z) possess real roots too. From the previous equations we get the following crossover inequalities: mutation δr (z0 ) > 0 Ki > 0 δ (z ) < 0 K i < a1 r 1 new population δr (z2 ) > 0 K i > a2 δr (z3 ) < 0 ⇒ K i < a3 (15) Figure 2: Steps of the genetic algorithm evolution δ (z ) > 0 K > a r 4 i 4 .. .. . . mization by the genetic algorithms is shown by figure 3. where the bounds aj for j = 1, 2, 3... are expressed by: aj = a(zj ) (16) Genetic Algorithm Now, according to these inequalities, it is clear that we Model need only odd bounds (which to say a1 ,a3 ). It has to be yc(t) e(t) u(t) y(t) strictly positive to get a feasible range for the controller C(s) G(s) parameter Ki . From (Silva et al, 2000, 2002), aj is posi- tive (for the odd values of j ) for every Kp verifying (7). Hence, the conditions giving by (15) are reduced to: 0 < Ki < min {aj } (17) j=1,3,5... Figure 3: The optimization principle by genetic algo- Once the stability domain is determined, the question is rithm what are the optimum parameters of the PI controller which guarantee the good performance of the closed- In our case, we are interest to search the optimum con- loop system? On the following, the genetic algorithm troller parameters in the area of stability using one of fol- is proposed to answer this need. lowing criterion ISE,IAE,IT AE and IT SE described by relation (23): GENETIC ALGORITHMS tP The Genetic Algorithms (AGs) are iterative algorithms of max ISE = e(t)2 global search of which the goal is to optimize a specific 0 tP function called criterion of performance or cost function max IAE = |e(t)| (Godelberg, 1991). In order to find the optimal solution 0 tP (18) of a problem by using AGs, we start by generating a pop- max IT AE = t |e(t)| ulation of individuals in a random way. The evolution 0 from one generation to the following is based on the use tmax te(t)2 P IT SE = of the three operators’ selection, crossover and mutation 0 If we want to minimize the tuning energy, the IT AE and the IAE criteria are considered. In the case where we privilege the rise time, we take the IT SE criterion. In order to guarantee the tuning energetic cost, we choose the ISE criterion ( Villain, 1996). The calculation steps of the control law are summarized by the following algorithm: 1. Introduction of the following parameters: Figure 4: PT-326 heat process trainer • maxpop individuals number by population • initial population tube, spaced at 28, 140 and 280 mm from both the heater • genmax generation number and the damper position. The system input, u(t) , is the voltage applied to the power circuit feeding the heating 2. Initialization of the generation counter (gen = 1) resistance, and the output,y(t) , is the outlet air tempera- ture, expressed by a voltage, between −10 and 10V . the 3. Initialization of the individual counter (j = 1 ) P T − 326 heating process is shown in figure 5. 4. For t = 1s to t = tmax efficiency evaluation of j th 1 Heater population individual f itness(J) = 1+J 5. Individual counter incrementing (j = j + 1 ). I II III • If j < maxpop , going back to step 4 Air output 30 • If not: application of the genetic operators (se- lection, crossover, mutation) for the founding a new population Air input Power Bridge 6. Generation counter incrementing (gen = gen + 1 ), Supply Circuit If , going back to step 3. 7. taking Kp opt and Ki opt which correspond to the u(t) y(t) best individual in the last population (individual making the best fitness). Figure 5: Schematic illustration On the following, the genetic algorithm characterized by The behavior of the P T − 326 thermal process is generation number equal to 100, Pc = 0.8 , Pm = 0.08 governed by the balance of heat energy. When the air and individual number by population equivalent to 20. temperature inside the tube is supposed to be uniform, a linear delay system model can be obtained. Thus, APPLICATIONS the transfer function between the heater input voltage This section presents an application example for the tem- and the sensor output voltage can be obtained as: y(s) K u(s) = 1+T s e −Ls perature control of an air stream heater (process trainer P T − 326), figure 4. This type of process is found in For the experiment, the damper position is set to many industrial systems such as furnaces, air condition- 30, and the temperature sensor is placed in the third ing, etc. The P T − 326 has the basic characteristics of position. The measurements acquisition is done by a large plant, with a tube through which atmospheric air ”P CI − DAS1002” card’s, the sampling time is taken is drawn by a centrifugal blower. Before being released equal to Te = 0.03 second. This choice takes account of into the atmosphere, air becomes heated by passing into the time computing and the constant time of the plant. a heater grid. Temperature control is realized by varying Firstly data is taken by operating the heater open loop the electrical power supplied to the heater grid. Air is with square input 0V/2V as shown as figure 6. The pushed through a tube by a fan blower and heated at the data is then used to obtain models to represent the plant inlet. The mass flow of air through the duct can be ad- An enlargement of the first step response up to 1000 justed by setting the opening of the throttle. There is an iterations is shown in the following figure. We define: energized electric resistance inside the tube. Heat, due to Step1: step response up to 1500 iterations the Joule effect, is released by resistance and transmitted Step3: step response from 3000 up to 4500 iterations by convection Step5: step response from 6000 up to 7500 iterations This process can be characterized as a linear time de- Step7: step response from 9000 up to 1000 iterations lay system. Time delay depends on the position of the temperature sensor element that can be inserted into the The transfer’s functions corresponding to each step re- air stream at any one of the three points throughout the sponse are estimated by conventional graphical method; 2.5 8 7 2 6 1.5 u(k), y(k) [V] 5 Ki 1 4 3 0.5 2 0 1 −0.5 0 0 2000 4000 6000 8000 10000 12000 −2 0 2 4 6 8 10 iteratons Kp Figure 6: Square input 0V /2V and outputs of P T − 326 Figure 8: PI controller stability domain for P T − 326 2.5 1. Initialization, k = 1 2 2. Output acquisition y(k) y(k) u(k) 1.5 u(k), y(k) [V] 3. Error computation e(k) = yc (k) − y(k) 1 4. Control law Computation 0.5 u(k) = u(k − 1) + Kpe(k) + (Ki Te − Kp )e(k − 1) 0 5. Application of the control law u(k) to the process −0.5 0 200 400 600 800 1000 iteratons 6. Exhausting of the sampling time, k = k + 1 then go to step 2. Figure 7: First step response for P T − 326 The followings figures show the responses of the closed loop systems in the case of PI controller designed by the and they are represented in table 1. So the P T − 326 genetic algorithms as described in table 2. 8 Table 1: system parameter Step 1 Step 3 Step 5 Step 7 mean 6 L 0.57 0.6 0.54 0.54 0.5625 K 0.587 0.586 0.5872 0.5872 0.5869 u(k),y(k) [V] 4 T 1.572 1.5671 1.572 1.572 1.572 2 y(s) is described by the following transfer function: u(s) = 0.58 1+1.57s e −0.56s 0 output control law In order to determine Kp values, we look for α in inter- val [0, π] satisfying tan(α) = 2.89α ⇒ α = 1.764. −2 0 500 1000 1500 2000 2500 3000 Kp range is given by:−2 < Kp < 8.25.The system iterations stability region, obtained in plane is presented in figure 8. From figure 8, our Kp and Ki population’s individ- Figure 9: Evolution of the output and the control law uals are choosing between [−2, 8.25] and [0, 7.45]. PI (Case: P I − ISE) controller optimum parameters supplied by genetic algo- rithm are presented by the following table. The following It is clear from these figures that the closed loop sys- tem is stable and the output y(k) tracks the step input sig- nal. From Table 3 we conclude that the minimal variance Table 2: optimum PI parameters of control law was generate by the P I − ISE controller. criterion ISE IAE IT AE IT SE Kp opt 3.67 3.33 2.94 3.95 Ki opt 4.24 4.36 4.04 4.33 Table 3: performance comparison ISE IAE IT AE IT SE algorithm describes the real-time implementation of the var(u(k)) 2.6708 2.779 2.8728 2.7234 PI control law: 8 REFERENCES Bhattacharyya, S. P. , Chapellat, H. , and Keel, L. H. (1995). 6 Robust Control:The Parametric Approach. Upper Saddle River, NJ: Prentice-Hall. [1, 2] u(k),y(k) [V] 4 Cedric, L. , (2002). Pontryagin revisit: implication sur la sta- 2 bilit des PID retards. Conference Internationale Franco- phone d’Automatique CIFA’02. [1] 0 output control law Dambrine, M. (1994). Contribution l’tude de la stabilit des systmes retards. Thse de doctorat, Universit de Lille. [1] −2 0 500 1000 1500 2000 2500 3000 iterations Godelberg, D.E.(1991). Genetic Algorithms in search, op- timization and machine learning. Addison-Wesley, Mas- Figure 10: Evolution of the output and the control law sachusetts. [3] (Case: P I − IAE) Lequin, O. , Gevers, M., Mossberg, M., Bosmans, E. and Triest , L. (2003). Iterative feedback tuning of PID parameters: 8 comparison with classical tuning rules. Control Engineering Practice,N11, pp. 1023-1033. [1] 6 Liu, G.P. and Daley, S. (2001). Optimal-tuning PID control for industrial systems. Control Engineering Practice,N 9, pp. u(k),y(k) [V] 4 1185-1194. [1] 2 Niculescu, S.-I. (2001). Delay effects on stability. Springer, London. [1] 0 control law output Silva, G. J. , Datta, A. and Bhattacharyya, S. P. (2000). Stabi- lization of Time Delay Systems. Proceedings of the Ameri- −2 0 500 1000 1500 2000 2500 3000 can Control Conference, pp. 963-970. [1, 2, 3] iterations Silva, G. J. , Datta, A. and Bhattacharyya, S. P. (2001). Sta- Figure 11: Evolution of the output and the control law bilization of First-order Systems with Time Delay using the (Case: P I − IT AE) PID controller. Proceedings of the American Control Con- ference, pp. 4650-4655. [2] 8 Silva, G. J. , Datta, A. and Bhattacharyya, S. P. (2002). New synthesis of PID controller. IEEE transactions on automatic 6 control,vol. 47,No. 2. [1, 2, 3] Silva, G. J. , Datta, A. and Bhattacharyya, S. P. (2005). PID u(k),y(k) [V] 4 controllers for time delay systems. Springer, London. [1, 2] 2 Villain, M. (1996). Systmes asservis linaires. ellipses/ dition marketing S.A. [4] 0 output control law Zhong, Q.C. (2006). Robust control of time delay system. Springer, London. [1] −2 0 500 1000 1500 2000 2500 3000 iteratons AUTHOR BIOGRAPHIES Figure 12: Evolution of the output and the control law RIHEM FARKH was born in Tunis, Tunisia, in 1982. (Case: P I − IT SE) He received the engineering degree in electrical engi- neering and the M.S degree in Control and Signal Pro- cessing from the National School Engineering of Tunis CONCLUSION in 2006 and 2007. In December 2007, he began pursu- In this work, we use the Hermit-Biehler theorem to com- ing his Ph.D. degree. His research interests include PID pute the region stability for first order delay system con- control, digital control and time-delay systems. trolled by PI regulator. Lastly, we were interested in KAOUTHER LAABIDI was born in Tunis, Tunisia. search of optimal PI for a given performance criteria She received the Master degree from the ”Ecole Suprieur (ISE, IAE, IT AE, IT SE), inside the stability region. de Sciences et de Technologie de Tunis”, in 1995 and the In regard to the complexity of the optimization problem, Ph.D degree in Genie Electic from the ”Ecole Nationale we used the genetic algorithms. The validation of these d’Ingnieurs de Tunis”, in 2005. She is currently prepar- results was tested in real time temperature control. ing the ability degree in the laboratory ACS ”Analyse et Commande des Systmes”.His research is related to the Identification and control of complex systems. MEKKI KSOURI received his B.Sc. degree in physics from the Faculty of Sciences of Tunis in 1971 and the Diplme d’Ingnieur from the Ecole Suprieure d’Electricit, Paris in 1973. He received the Diplme de Docteur Ing- nieur and the Diplme de Docteur d’Etat from the Uni- versit de Paris VI respectively in 1975 and 1977. He received the Outstanding Contribution Award for Lead- ership in the organization of the CIFA 2004 Sympo- sium, the Chevalier dans l’Ordre des Palmes Acad- miques, France (2002), Outstanding Contribution Award IEEE SMC 1998 and the Ordre National du Mrite de l’Education, Tunisia (1997). He is senior member of IEEE. His research interests include identification and control of nonlinear systems and industrial applications of advanced control. He is the author or co-author of more than 150 papers in international conferences and journals. He has also published six books.