Evolutionary algorithms based design of multivariable PID controller Expert Systems with Applications 36 (2009) 9159–9167 Contents lists available at ScienceDirect Expert Systems with Applications j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / e s w a Evolutionary algorithms based design of multivariable PID controller M. Willjuice Iruthayarajan *, S. Baskar Department of EEE, Thiagarajar College of Engineering, Madurai 625 015, Tamilnadu, India a r t i c l e i n f o Keywords: PID Control Evolutionary algorithm MIMO system On-line tuning Off-line tuning 0957-4174/$ - see front matter � 2008 Elsevier Ltd. A doi:10.1016/j.eswa.2008.12.033 * Corresponding author. Tel.: +91 94434 87093. E-mail address: willjuice@tce.edu (M.W. Iruthayar a b s t r a c t In this paper, performance comparison of evolutionary algorithms (EAs) such as real coded genetic algo- rithm (RGA), modified particle swarm optimization (MPSO), covariance matrix adaptation evolution strategy (CMAES) and differential evolution (DE) on optimal design of multivariable PID controller design is considered. Decoupled multivariable PI and PID controller structure for Binary distillation column plant described by Wood and Berry, having 2 inputs and 2 outputs is taken. EAs simulations are carried with minimization of IAE as objective using two types of stopping criteria, namely, maximum number of func- tional evaluations (Fevalmax) and Fevalmax along with tolerance of PID parameters and IAE. To compare the performances of various EAs, statistical measures like best, mean, standard deviation of results and average computation time, over 20 independent trials are considered. Results obtained by various EAs are compared with previously reported results using BLT and GA with multi-crossover approach. Results clearly indicate the better performance of CMAES and MPSO designed PI/PID controller on multivariable system. Simulations also reveal that all the four algorithms considered are suitable for off-line tuning of PID controller. However, only CMAES and MPSO algorithms are suitable for on-line tuning of PID due to their better consistency and minimum computation time. � 2008 Elsevier Ltd. All rights reserved. 1. Introduction Proportional-integral-derivative (PID) control offers the sim- plest and yet most efficient solution to many real-world control problems. Three-term functionality of PID controller covers treat- ment of both transient and steady state responses. The popularity of PID control has grown tremendously, since the invention of PID control in 1910 and the Ziegler–Nichol’s straight forward tuning method in 1942. With the advances in digital technology, the sci- ence of automatic control now offers a wide spectrum of choices for control schemes such as adaptive control (Astrom & Witten- mark, 1995), neural network control (Fukuda & Shibata, 1992) and fuzzy logic control (Lee, 1990). However more than 90% of industrial controllers are still implemented based around PID con- trol algorithms, as no other controllers match the simplicity, clear functionality, applicability and ease of use offered by the PID con- trollers (Ang, Chang, & Li, 2005). Several approaches have been reported in literature for tuning the parameters of PID controllers. Ziegler–Nichols and Cohen– Coon are the most commonly used conventional methods for tuning PID controllers and neural network, fuzzy based approach, neuro-fuzzy approach and evolutionary computation techniques are the recent methods (Astrom & Hagglund, 1995). ll rights reserved. ajan). Many researches have already reported the optimal design of PID controller parameters using various EAs such as GA (Chen, Cheng, & Lee, 1995), MPSO (Gaing, 2004; Ghoshal, 2004; Mukher- jee & Ghoshal, 2007; Wang, Zhang, & Wang, 2006) and DE (Bingul, 2004) for SISO system. In general, EAs are robust search and optimization methodology, able to cope with ill-defined problem domain such as multimodality, discontinuity, time-variance, randomness and noise. GA approach for tuning of PID controllers for multi-input multi-output (MIMO) process is also reported (Chang, 2007; Zuo, 1995). In Chang (2007), decoupled multivariable PI controller tuning using GA with multi-parent crossover approach was presented. Simple three-parent differential crossover and uniform mutation operators have been employed. The better performance of three- parent crossover RGA over BLT and traditional two-parent cross- over based RGA was demonstrated in the paper. Recently, several modifications are carried out in crossover and mutation mechanisms of RGA such as SBX crossover, PCX crossover and non-uniform polynomial mutation to improve the perfor- mance of RGA. Self-adaptive simulated binary crossover (SBX) based RGA was successfully applied to various engineering optimi- zation problems (Deb, 2001). SBX crossover is self-adaptive in nat- ure which creates children solutions in proportion to the difference in parent solutions. The near parent solutions are monotonically more likely to be chosen as offspring than solutions distant from parents. mailto:willjuice@tce.edu http://www.sciencedirect.com/science/journal/09574174 http://www.elsevier.com/locate/eswa 9160 M.W. Iruthayarajan, S. Baskar / Expert Systems with Applications 36 (2009) 9159–9167 Another EA, namely, covariance matrix adaptation evolution strategy (CMAES) with the ability of learning of correlations be- tween parameters and the use of the correlations to accelerate the convergence of the algorithm is recently proposed. Due to the learning process, the CMAES algorithm performs the search independent of the coordinate system, reliably adapts topologies of arbitrary functions, and significantly improves convergence rate especially on non-separable and/or badly scaled objective func- tions. CMAES algorithm has been successfully applied in varieties of engineering optimization problems (Baskar, Alphones, Sugan- than, Ngo, & Zheng, 2005).This algorithm outperforms all other similar classes of learning algorithms on the benchmark multi- modal functions (Kern et al., 2004). Covariance matrix adaptation evolution strategy algorithm and also recent modifications in other EAs were not applied for the tun- ing of PID controllers. Also, all the reported papers for EA based PID controller design, have considered one or two algorithms for the purpose of comparison. This paper focuses mainly on the performance evaluation of various EAs such as Self-adaptive RGA, MPSO, DE and CMAES on optimum design of multivariable PI and PID controllers for binary distillation column plant described by Wood and Berry (Chang, 2007). The essence of the paper lies in the determination of suit- able EA method for the tuning of PID controller for MIMO system. The remaining part of the paper is organized as follows. Section 2 introduces PID controller structure for SISO and MIMO systems. Section 3 describes the various EAs methods. Section 4 introduces the MIMO system considered for PID controller tuning. Section 5 presents the implementation of EA based multivariable PID con- troller design. Section 6 reveals the simulation results. Finally, con- clusions are given in Section 7. 2. PID controller structure A standard PID controller structure is also known as the ‘‘three- term” controller, whose transfer function is generally written in the ideal form in (1) or in the parallel form in (2) GðsÞ¼ K P 1 þ 1 T Is þ T Ds � � ð1Þ GðsÞ¼ K P þ K I s þ K Ds ð2Þ where KP is the proportional gain, TI is the integral time constant, TD is the derivative time constant, K I ¼ K P=T I is the integral gain and K D ¼ K P T D is the derivative gain. The ‘‘three-term” functionalities are highlighted below. � The proportional term – providing an overall control action pro- portional to the error signal through the all pass gain factor. � The integral term – reducing steady state errors through low fre- quency compensation by an integrator. � The derivative term – improving transient response through high frequency compensation by a differentiator. For optimum performance, KP, KI (or TI) and KD (or TD) are tuned by EAs by minimizing the performance measures such as IAE, ISE and ITAE. Yd + - Multivariable PID controllers K(s) E Fig. 1. A multivariable P 2.1. PID controller for MIMO system Consider a multivariable PID control structure as in Fig. 1, where, desired output vector: Yd = [yd1, yd2, . . ., ydn] T; Actual output vector: Y = [y1, y2, . . ., yn] T; Error vector: E = Yd � Y = [yd1 � y1, yd2 � y2, . . ., ydn � yn]T = [e1, e2, . . ., en] T; Control input vector: U = [u1, u2, . . ., un] T; n � n Multivariable processes: GðsÞ¼ g11ðsÞ � � � g1nðsÞ .. . . . . .. . gn1ðsÞ � � � gnnðsÞ 2 664 3 775 ð3Þ n � n Multivariable PID controller: KðsÞ¼ k11ðsÞ � � � k1nðsÞ .. . . . . .. . kn1ðsÞ � � � knnðsÞ 2 664 3 775 ð4Þ In this work, decoupled multivariable PID controller is consid- ered. So K(s) becomes KðsÞ¼ k1ðsÞ � � � 0 .. . . . . .. . 0 � � � knðsÞ 2 664 3 775 ð5Þ The form of ki(s) is either in (1) or (2). In this work, ‘‘parallel for- m” of PID controller in (2) is used and can be rewritten as kiðsÞ¼ kPi þ kIi s þ kDi s; ð6Þ For convenience, let hi ¼ ½kPi ; kIi ; kDi�, represents the gains vector of ith diagonal sub PID controller in K(s) .For multivariable PI con- troller, ki(s) in (6) can be rewritten as kiðsÞ¼ kPi þ kIi s ð7Þ hi ¼ ½kPi ; kIi� represents the gains vector of the ith diagonal sub PI controller in K(s). 2.2. Performance index In the design of PID controller, the performance criterion or objective function is first defined based on the desired specifica- tions such as time-domain specifications, frequency domain specifications and time-integral performance. The commonly used time-integral performance indexes are integral of the square error (ISE), integral of the absolute value of the error (IAE) and integral of the time-weighted absolute error (ITAE). Minimization of IAE as given in (8) is considered as the objective of this paper IAE ¼ Z 1 0 ðje1ðtÞjþ je2ðtÞjþ � � �þ jenðtÞjÞdt ð8Þ Multivariable Processes G(s) YU ID control system. M.W. Iruthayarajan, S. Baskar / Expert Systems with Applications 36 (2009) 9159–9167 9161 3. Evolutionary algorithms Evolutionary algorithms differ from the traditional optimization techniques in that EAs make use of a population of solutions, not a single point solution. EAs are inherently parallel, because they simultaneously evaluate many points in the parameter space (search space). Considering many points in the search space, EA has a reduced chance of converging to the local optimum and would be more likely to converge to the global optimum. An iter- ation of EA involves a competitive selection that weeds out poor solutions and offspring generation mechanism. Several evolution- ary search algorithms like GA, EP, MPSO, DE and CMAES were developed independently. These algorithms differ in selection, off- spring generation and replacement mechanisms. For solving global functional optimization problems, RGA, MPSO DE and CMAES algo- rithms are normally used and hence these algorithms are em- ployed in this paper. Over the years, several variants of these algorithms were proposed. Some of the recent variants of these algorithms are briefly explained in this section. 3.1. Real coded genetic algorithm with SBX crossover In general, a genetic algorithm has five components as follows: 1. A genetic representation of solutions to the problem. 2. A way to create an initial population of solutions. 3. An evaluation function rating solution in terms of its fitness. 4. Parent selection mechanism and genetic operators that alter the genetic composition of children during reproduction. 5. Values for the parameter of genetic algorithm. Real-number encoding is best used for function optimization problems. It has been widely confirmed that real-number encoding performs better than binary or gray encoding for constrained opti- mization. Owing to the adaptive capability, SBX crossover and polynomial mutation operators are employed. Tournament selec- tion is used as selection mechanism in order to avoid premature convergence. Simulated binary crossover (SBX) and polynomial mutation are briefly explained below. 3.1.1. Simulated binary crossover In SBX crossover (Deb, 2001), two children solutions are created from two parents as follows: Choose a random number ui, e [0, 1] and calculate bqi as given in (9) bqi ¼ ð2uiÞ 1 gcþ1; ui 6 0:5 1 2ð1�uiÞ � � 1 gcþ1 ; otherwise 8< : ð9Þ A spread factor bqi is defined as the ratio of the absolute differ- ence in offspring values to that of the parents. gc is the crossover index. Then compute the offspring xð1;tþ1Þi & x ð2;tþ1Þ i as xð1;tþ1Þi ¼ 0:5 ð1 þ bqiÞx ð1;tÞ i þð1 � bqiÞx ð2;tÞ i h i xð2;tþ1Þi ¼ 0:5 ð1 � bqiÞx ð1;tÞ i þð1 þ bqiÞx ð2;tÞ i h i : ð10Þ 3.1.2. Polynomial mutation Newly generated offspring undergoes polynomial mutation operation. Like in the SBX operator, the probability distribution can also be a polynomial function, instead of a normal distribution. The new offspring yð1;tþ1Þi is determined as follows: yð1;tþ1Þi ¼ x ð1;tþ1Þ i þðx U i � x L i Þdi: ð11Þ xUi and x L i are the upper and lower limit values.where the parameter di is calculated from the polynomial probability distribution PðdÞ¼ 0:5ðgm þ 1Þð1 �jdjÞ gm di ¼ ð2riÞ 1=ðgmþ1Þ � 1; if ri < 0:5 1 �½2ð1 � riÞ� 1=ðgmþ1Þ; if ri P 0:5 ( ð12Þ where gm is the mutation index. Newly generated individuals replace their parents and form the parents for the next generation and the above procedure is re- peated until a maximum number of function evaluations are completed. 3.2. Modified particle swarm optimization (MPSO) Particle swarm optimization provides a population-based search procedure in which individuals called particles, change their positions (states) with time. In a PSO system, particles fly around in a multidimensional search space. During flight, each particle ad- justs its position according to its own experience, and the experi- ence of neighboring particles, making use of the best position encountered by itself and its neighbors. The swarm direction of a particle is defined by the set of particles neighboring the particle and its history of experience. Let Xi and Vi represent ith particle position and its correspond- ing velocity in a search space. The best previous position of ith par- ticle is recorded and represented as Pbesti. The best particle among all the particles in the group is represented as Gbest. The updated velocity of individual i is given in (13) (Shi & Eberhart, 1998) V kþ1i ¼ x kV ki þ c1 rand1 Pbest k i � X k i � � þ c2 rand2 Gbest k � Xki � � ð13Þ where V ki Velocity of ith individual at iteration k xk Inertia weight at iteration k c1, c2 Acceleration factors rand1, rand2 Uniform random numbers between 0 and 1 Xki Position of ith individual at iteration k Pbestki Best position of ith individual at iteration k Gbestk Best position of the group until iteration k Each individual moves from the current position to the next one with the modified velocity as Xkþ1i ¼ X k i þ V kþ1 i ð14Þ The above procedure is repeated until a maximum number of func- tional evaluations have been reached. 3.3. Differential evolution Differential Evolution is a simple population-based, stochastic parallel search evolutionary algorithm for global optimization (Storm & Price, 1997). The initial population is chosen randomly and should cover the entire parameter space. At each generation, DE employs both mutation and crossover (recombination) to pro- duce one trial vector Ui,J+1 for each target vector Xi,J. Then, a selec- tion phase takes place, where each trial vector is compared with the corresponding target vector; the better one will enter the pop- ulation of the next generation. For each target vector Xi,J, a mutant vector is generated using Rand/Rand strategy as Ui;Jþ1 ¼ Xr3;J þ FðXr1;J � Xr2;JÞ ð15Þ 9162 M.W. Iruthayarajan, S. Baskar / Expert Systems with Applications 36 (2009) 9159–9167 where Ui;Jþ1 ith mutated individual for the next iteration X Population set F Mutation constant [0, 2] J Current iteration J+1 Next iteration Xr1,J. . . Xr3,J Randomly selected individuals from the population of the current iteration. To increase the diversity of the perturbed parameter vectors cross- over is performed after mutation V i;Jþ1 ¼ Xi;jð1 � CRÞþ Ui;Jþ1CR ð16Þ where, CR = crossover probability constant from interval [0, 1] The parents for the next iteration are selected as follows: Xi;Jþ1 ¼ V i;Jþ1 if fðV i;Jþ1Þ > fðXi;Jþ1Þ Xi;Jþ1 if fðXi;Jþ1Þ > fðV i;Jþ1Þ � ð17Þ where f(Vi,J+1) is the fitness function value of the ith individual of the population to which the mutation and crossover operators are ap- plied and f(Xi,J+1) is the fitness function value of the ith individual in the original population. The loss of the best individuals in the fol- lowing iteration is avoided by this selection mechanism, as the worst individuals are replaced by the best individuals. This process continues until the maximum function evaluation is reached. 3.4. Covariance matrix adaptation evolution strategy (CMAES) Covariance matrix adaptation evolution strategy is a class of continuous EA that generates new population members by sam- pling from a probability distribution that is constructed during the optimization process. One of the key concepts of this algorithm involves the learning of correlations between parameters and the use of the correlations to accelerate the convergence of the algo- rithm. The adaptation mechanism of CMAES consists of two parts, (1) the adaptation of the covariance matrix C and (2) the adapta- tion of the global step size. The covariance matrix C is adapted by the evolution path and difference vectors between the l best individuals in the current and previous generation. The detailed CMAES algorithm is presented in Kern et al. (2004). 3.4.1. CMAES algorithm Step 1: Generate an initial random solution. Step 2: The offspring at g + 1 generation xgþ1k are sampled from a Gaussian distribution using covariance matrix and global step size at generation g xðgþ1Þk ¼ zk; zk ¼ N hxi ðgÞ l ; r ðgÞ2 CðgÞ � � k ¼ 1; :::; k ð18Þ where hxiðgÞl ¼ Pl i¼1x ðgÞ i with l being the selected best individuals from the population. The parameters cc, ccov, cr and d required for further computa- tions are by default given in terms of the number of decision vari- ables (n) and l as follows: cc ¼ 4 n þ 4 ; cr ¼ 10 n þ 20 ; d ¼ max 1; 3l n þ 10 � � þ cr; ccov ¼ 1 l 2 ðn þ ffiffiffiffiffi 2Þ p 2 þ 1 � 1l � � min 1; 2l � 1 ðn þ 2Þ2 þ l ! ð19Þ The parameters cr and ccov control independently the adaptation time scales for the global step size and the covariance matrix. Note that if l � n, d is large and the change in r is negligible compared to that of C. The initial values are Pð0Þr ¼ P ð0Þ c ¼ 0 and Cð0Þ ¼ I. Step 3: The evolution path Pðgþ1Þc is computed as follows: Pðgþ1Þc ¼ð1 � ccÞ � P ðgÞ c þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ccð2 � ccÞ p � ffiffiffiffi l p rðgÞ hxiðgþ1Þl �hxi g l � � ð20Þ Cðgþ1Þ ¼ ð1 � ccovÞ � CðgÞ þ ccov � 1 l Pðgþ1Þc P ðgþ1Þ c � �T þ 1 � 1 l � � 1 l Xl i¼1 1 rðgÞ2 xðgþ1Þi �hxi ðgÞ l � � xðgþ1Þi �hxi ðgÞ l � �T� ð21Þ The strategy parameter ccov 2 ½0; 1� determines the rate of change of the covariance matrix C. Step 4: Adaptation of global step size r(g+1) is based on a conju- gate evolution path Pðgþ1Þr Pðgþ1Þr ¼ð1 � crÞ � P ðgÞ r þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi crð2 � crÞ p � BðgÞðDðgÞÞ�1ðBðgÞÞ�1 � ffiffiffiffi l p rðgÞ hxiðgþ2Þl �hxi g l � � ð22Þ the matrices B(g) and D(g) are obtained through a principal compo- nent analysis: CðgÞ ¼ BðgÞðDðgÞÞ2ðBðgÞÞT ð23Þ where the columns of B(g) are the normalized eigen vectors of C(g) and D(g) is the diagonal matrix of the square roots of the given eigen values of C(g). The global step size r(g+1) is determined by rðgþ1Þ ¼ rðgÞ exp cr d kPðgþ1Þr k EðkNð0; IÞkÞ ! � 1 ! ð24Þ Step 5: Repeat Steps 2–4 until the stopping criteria are satisfied. 4. MIMO system Most of the industrial processes belong to the category of MIMO system, which requires MIMO control techniques to improve per- formance, even though they are naturally more difficult to exploit than SISO system. Binary distillation column plant described by Wood and Berry (Chang, 2007; Wang, Zou, Lee, & Qiang, 1997) is considered. The transfer function of the above process is given in (25) GðsÞ¼ 12:8e�s 1þ16:7s �18:9e�3s 1þ21s 6:6e�7s 1þ10:9s �19:4e�3s 1þ14:4s " # ð25Þ The transfer function concerned with multivariable process has first order dynamics and significant time delays and it has a strong interaction between inputs and outputs. In this paper, multivari- able controller with PI and PID structures are used for optimizing the IAE performance for set point regulation using EAs. 5. EA implementation of multivariable PID controller System described by (25), has two inputs and two outputs. The decoupled multivariable PID controller K(s) for this system is given in (26) KðsÞ¼ k1ðsÞ 0 0 k2ðsÞ � ð26Þ In order to obtain the optimum performance, the parameters of K(s) i.e., ½h1; h2� ¼ ½kP1 ; kI1 ; kD1 ; kP2 ; kI2 ; kD2� are optimized by optimi- zation algorithms. The chromosome/particle representation is gi- ven in Fig. 2. First three elements are the parameters of k1(s) and next three elements are the parameters of k2(s). For multivariable PI control- ler, the structure of the chromosome is given in Fig. 3. First two 1P k 1I k 1D k 2P k 2I k 2D k Fig. 2. Chromosome/particle of multivariable PID controller. 1P k 1I k 2P k 2I k Fig. 3. Chromosome/particle of multivariable PI controller. Table 1 Parameter selection. Evolutionary algorithms Parameter RGA Pc = 0.8 Pm = 1/n gc = 5 gm = 20 MPSO C1 = 1 C2 = 1 Vmax = 0.1 x is linearly decreased from 0.9 to 0.2 over the iterations DE CR = 0.5 F = 0.8 CMAES Parameters are fixed automatically by the algorithm Table 2 Optimum parameters of multivariable PI controller – without tolerance. Method Optimum parameters of multivariable PI controller IAE kP1 kI1 kP2 kI2 BLT* 0.3750 0.0452 �0.0750 �0.0032 23.5568 RGA-multi-crossover* 0.9971 0.0031 �0.0141 �0.0071 10.5778 RGA-SBX crossover 0.8433 0.0026 �0.0127 �0.0069 10.4395 MPSO 0.8485 0.0026 �0.0132 �0.0069 10.4378 DE 0.8485 0.0026 �0.0132 �0.0069 10.4378 M.W. Iruthayarajan, S. Baskar / Expert Systems with Applications 36 (2009) 9159–9167 9163 elements are the parameters of k1(s) and the next two elements are the parameters of k2(s). All the elements of chromosomes/particles of population are randomly initialized within the search space specified by their lower and upper bounds of individual parameters as given in (Chang, 2007). The inequality conditions for the parameter ranges of PI and PID controllers are given in (27) and (28), respectively �1 6 K P1 6 1; �1 6 K I1 6 1; �1 6 K P2 6 1; �1 6 K I2 6 1 ð27Þ �1 6 K P1 6 1; �1 6 K I1 6 1; �1 6 K D1 6 1; �1 6 K P2 6 1; �1 6 K I2 6 1; �1 6 K D2 6 1: ð28Þ CMAES 0.8485 0.0026 �0.0132 �0.0069 10.4378 * Data taken from Chang (2007). Table 3 Optimum parameters of multivariable PI controller – with tolerance. Method Optimum parameters of multivariable PI controller IAE kP1 kI1 kP2 kI2 BLT* 0.3750 0.0452 �0.0750 �0.0032 23.5568 RGA-multi-crossover* 0.9971 0.0031 �0.0141 �0.0071 10.5778 RGA-SBX crossover 0.8622 0.0026 �0.0135 �0.0069 10.44 MPSO 0.8485 0.0026 �0.0132 �0.0069 10.4378 DE 0.8485 0.0026 �0.0132 �0.0069 10.4379 CMAES 0.8485 0.0026 �0.0132 �0.0069 10.4378 * Data taken from Chang (2007). Table 4 Statistical performance of EAs of multivariable PI controller – without tolerance. Method Best value Mean value Standard deviation Average computation time (s) RGA-SBX crossover 10.4395 10.9366 1.6638 137.8698 MPSO 10.4378 10.8157 1.6570 135.2763 DE 10.4378 10.4379 4.5272E-5 137.9376 CMAES 10.4378 10.4378 0 142.5264 Table 5 Statistical performance of EAs of multivariable PI controller – tolerance. Method Best value Mean value Standard deviation Average computation time (s) Average functional evaluations RGA-SBX crossover 10.44 10.4961 0.0810 134.8504 5840 MPSO 10.4378 10.5394 0.3410 110.3138 4805 DE 10.4379 11.2228 1.5740 107.3965 4676 CMAES 10.4378 10.4378 0 84.8255 3572 6. Simulation results In this work, two experiments namely design of multivariable PI and PID controller for binary distillation column plant described by Wood and Berry, using various EAs such as RGA, MPSO, DE and CMAES are conducted. For simulating Binary distillation column plant, MATLAB-SIMULINK software is employed. Simulations are carried out using Core 2 Duo Processor 2.2 GHz, 2GB RAM PC. IAE is determined for step response over 150 min time period. EAs sim- ulations are carried out using two types of stopping criteria namely, Fevalmax, Fevalmax with tolerance of design variables and tolerance of objective function. In this work, the Fevalmax is set at 6000 functional evaluations and tolerance is fixed as 10�5 for the last 20 generations. Owing to the randomness of the EAs, many trials with independent population initializations should be made to acquire a useful conclusion of the performance of the algorithm. Hence, best, mean, standard deviation of IAE measure and average computation in 20 independent trials of various EAs are reported and compared with the already reported results (Chang, 2007). 6.1. Parameter tuning Optimal parameter combinations for different EA methods are experimentally determined by conducting a series of experiments with different parameter settings before conducting actual runs to collect the results. The parameters actually used in the simulations are summarized in the Table 1. 6.2. Tuning of multivariable PI controller Best PI parameters and the corresponding IAE values for the 20 trials of multivariable PI controllers using different EAs with and without tolerance in stopping criteria are reported in Tables 2 and 3, respectively. For the purpose of comparison, already re- ported values obtained by conventional BLT method and GA with multi-crossover approach are directly taken from (Chang, 2007) and given in the Tables. 9164 M.W. Iruthayarajan, S. Baskar / Expert Systems with Applications 36 (2009) 9159–9167 Almost all the EAs are giving equal performance with respect to the best IAE. The statistical performances such as best, mean, stan- dard deviation of IAE and average computation time (ACT) of 20 trials using various EAs with and without considering tolerance are given in Tables 4 and 5, respectively. Without considering tol- erance in stopping criteria, CMAES and DE are the same with re- spect to mean value but CMAES is computationally slightly expensive due to complex mathematical manipulations. While 0 20 40 60 80 1 10 1 10 2 10 3 Generation IAE Fig. 4. Convergence characteristics of 0 50 100 150 200 0.7 0.8 0.9 1 Iteration K p1 1 0 50 100 150 200 -0.2 0 0.2 0.4 0.6 Iteration K p2 2 Fig. 5. Convergence characteristics considering tolerance, CMAES algorithm gives better performance in consistently achieving good results as compared to all other algorithms and also it requires less number of average functional evaluations (AFeval). MPSO algorithm gives better performance than DE and RGA. Fig. 4 shows the convergence characteristics of various EAs con- sidering tolerance. Higher value of IAE during the initial genera- tions/iterations indicates unstable PID control settings. For 00 120 140 160 180 200 / Iteration RGA-SBX MPSO DE CMAES EAs – multivariable PI controller. 0 50 100 150 200 -0.5 0 0.5 1 Iteration K i1 1 0 50 100 150 200 -0.15 -0.1 -0.05 0 0.05 Iteration K i2 2 of PID parameters using MPSO. 0 50 100 150 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time in min O ut pu t re sp on se y 1 BLT RGA-MC PI-CMAES Fig. 6. Output response y1. 0 50 100 150 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time in min O ut pu t re sp on se y 2 BLT RGA-MC PI-CMAES Fig. 7. Output response y2. M.W. Iruthayarajan, S. Baskar / Expert Systems with Applications 36 (2009) 9159–9167 9165 convenience, IAE for the unstable response is limited to 1000. Other than CMAES algorithm, the convergence characteristics of all algorithms are smooth. Due to the self-learning behavior of CMAES algorithm, convergence characteristics show large varia- tions during the initial search. Fig. 5 shows the convergence char- acteristics of multivariable PI parameters obtained by MPSO algorithm with tolerance in stopping criteria. Output responses y1 and y2 of the MIMO system with multivar- iable PI controller using best PI parameter obtained out of the 20 trials by using CMAES/MPSO algorithm is shown in Figs. 6 and 7, respectively. For comparison purpose, output responses using RGA with multi-crossover approach are also given in the same fig- ure. PI controller simulation results show the better performance of CMAES and MPSO algorithms as compared to the previously re- ported results obtained using GA with multi-crossover approach (Chang, 2007) and conventional BLT method. Time-responses spec- ifications for the designed PI controllers such as peak overshoot in%, rise time (min) and settling time for ±5% tolerance (min) for output responses y1 and y2 are summarized in Table 6. Peak over- shoot of system response with multivariable PI controller designed by CMAES/MPSO is approximately 50% lesser than GA with multi- crossover approach with a slight increase in rise time and settling time. 6.3. Tuning of multivariable PID controller Best PID parameters and the corresponding IAE values for the 20 trials of multivariable PID controllers using various EAs with and without considering tolerance in stopping criteria are reported in Tables 7 and 8, respectively. The CMAES and MPSO algorithms are almost giving equal performance with respect to the best IAE. The statistical performances such as best, mean, standard deviation of IAE, ACT and AFeval of 20 trials using various EAs with and with- out considering tolerance are given in Tables 9 and 10, respec- tively. CMAES, MPSO and RGA algorithms give better performance in consistently achieving good results as compared to DE. But AFeval required and ACT for RGA is larger than CMAES and MPSO. For clarity, transient portion (25 min) of output responses y1 and y2 of the MIMO system with multivariable PID controller using best PID parameter obtained using CMAES/MPSO algorithm is shown in Figs. 8 and 9, respectively. For comparison purpose, out- put responses using RGA with multi-crossover approach and mul- tivariable PI controller designed by CMAES are also given in the same figure. PID controller simulation results show the better per- formance of CMAES and MPSO algorithms as compared to the pre- Table 6 Time-response specifications of output responses – multivariable PI controller. Method y1 %Mp Rise time (min) Settling tim BLT 27.4365 4.9 32.2 RGA with multi-crossover 21.9866 2.5 7.1 CMAES 10.8390 2.8 7.3 Table 7 Optimum parameters of multivariable PID controller – without tolerance. Method Optimum PID parameters kP1 kI1 kD1 RGA-SBX crossover 1.0 0.0025 0.3892 PID-MPSO 1.0 0.0025 0.3872 DE 0.9945 0.0026 0.4021 CMAES 1.0 0.0025 0.3872 viously reported results obtained using GA with multi-crossover approach (Chang, 2007) and multivariable PI controller designed y2 e 5% (min) %Mp Rise time (min) Settling time 5% (min) 30.9972 9.2 88.2 17.8161 8.5 14.9 9.2414 8.9 10.6 IAE kP2 kI2 kD2 �0.0317 �0.0072 �0.0885 9.6859 �0.0332 �0.0073 �0.0909 9.6824 �0.0289 �0.0071 �0.0709 9.7108 �0.0332 �0.0073 �0.0909 9.6824 Table 8 Optimum parameters of multivariable PID controller – with tolerance. Method Optimum PID parameters IAE kP1 kI1 kD1 kP2 kI2 kD2 RGA-SBX crossover 1.0 0.0025 0.3860 �0.0358 �0.0073 �0.0991 9.6871 PID-MPSO 1.0 0.0025 0.3872 �0.0332 �0.0073 �0.0909 9.6824 DE 0.9825 0.0026 0.4239 �0.0563 �0.0082 �0.1563 9.9251 CMAES 1.0 0.0025 0.3872 �0.0332 �0.0073 �0.0909 9.6824 Table 9 Statistical performance of EAs of multivariable PID controller – without tolerance. Method Best value Mean value Standard deviation Average computation time RGA-SBX crossover 9.6859 10.2370 0.6968 136.6054 MPSO 9.6824 10.0619 0.4304 136.6904 DE 9.7479 9.8855 0.1135 138.9478 CMAES 9.6824 10.4451 0.7116 146.4358 Table 10 Statistical performance of EAs of multivariable PID controller – with tolerance. Method Best value Mean value Standard deviation Average computation time (s) Average functional evaluations RGA-SBX crossover 9.6871 10.0954 0.3128 139.4493 5943 MPSO 9.6824 10.1992 0.6073 133.0369 5705 DE 9.9251 46.4502 85.6846 71.8283 3075 CMAES 9.6824 10.3186 0.7618 138.6396 5677 0 5 10 15 20 25 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time in min O ut pu t re sp on se y 2 RGA-MC PI-CMAES PID-CMAES Fig. 9. Transient portion of output response y2. 0 5 10 15 20 25 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time in min O ut pu t re sp on se y 1 RGA-MC PI-CMAES PID-CMAES Fig. 8. Transient portion of output response y1. Table 11 Time-response specifications of output responses – multivariable PID controller. Method y1 y2 %Mp Rise time (min) Settling time 5% (min) %Mp Rise time (min) Settling time 5% (min) RGA-SBX 0.9602 2.8 12 11.9161 8.3 10.9 DE 2.1241 5.1 12.3 20.7696 8.1 11.7 CMAES/MPSO 0.8059 2.9 11.9 10.4336 8.4 10.6 9166 M.W. Iruthayarajan, S. Baskar / Expert Systems with Applications 36 (2009) 9159–9167 by CMAES. Time responses of y1 and y2 using best PID parameters using various EAs are summarized in Table 11. Results show the better performance of CMAES/MPSO in terms of very less peak overshoot and settling time of response y1 and y2 with slight in- crease in rise time. 7. Conclusions In this paper, performance evaluation of EAs such as RGA with SBX crossover, MPSO, DE and CMAES on the optimal design of mul- tivariable PI and PID controller for the binary distillation column plant is conducted. EAs simulations are carried out using two types of stopping criteria, namely, Fevalmax and Fevalmax along with tolerance of PID parameters and IAE. Multivariable PI/PID control- lers are designed by minimizing IAE and the results are compared with those of the already reported in literature namely, BLT and GA with multi-crossover approach. Simulation results are summarized as follows: (i) The better performance of evolutionary designed multivari- able PI controller over the already reported results. Also, multivariable PID controllers designed for the same system by various EAs are better than multivariable PI controller. (ii) Without tolerance, the best IAE obtained in 20 independent trials by all EAs is almost equal. This reveals that all EAs are equally applicable to off-line PID controller tuning. (iii) Considering the tolerance of PID parameters and IAE, CMAES and MPSO algorithms are more suitable for on-line tuning of PID controller due to their better consistency and minimum computation time. Also, MPSO is much more suitable for on- line tuning of PID controller due to the reduced computation time. M.W. Iruthayarajan, S. Baskar / Expert Systems with Applications 36 (2009) 9159–9167 9167 Acknowledgments The authors gratefully acknowledge the Management of the Thiagarajar College of Engineering, Madurai 625 015, Tamilnadu, India, for their continued support, encouragement, and the exten- sive facilities provided to carry out this research work. They also gratefully acknowledge the support of Dr. M. Chidambaram, Direc- tor, NIT, Trichy. References Ang, K. H., Chang, G., & Li, Yun (2005). PID control system analysis, design and technology. IEEE Transaction on Control System Technology, 13(4), 559–577. Astrom, K. J., & Wittenmark, B. (1995). Adaptive control (2nd ed.). Addison Wesley. Astrom, K. J., & Hagglund, T. (1995). PID controllers: theory, design, and tuning (2nd ed.). Instrument society of America. Baskar, S., Alphones, A., Suganthan, P. N., Ngo, N. Q., & Zheng, R. T. (2005). Design of optimal length low-dispersion FBG filter using covariance matrix adapted evolution. IEEE Photonics Technology Letters, 17(10), 2119–2121. Bingul, Z. (2004). A new PID tuning technique using differential evolution for unstable and integrating processes with time delay. ICONIP, Proceedings Lecture Notes in Computer Science, 3316, 254–260. Chang, W. D. (2007). A multi-crossover genetic approach to multivariable PID controllers tuning. Expert Systems with Applications, 33, 620–626. Chen, B. S., Cheng, Y. M., & Lee, C. H. (1995). A genetic approach to mixed H2/H1 Optimal PID control. IEEE Control Systems, 15(5), 51–60. Deb, K. (2001). Multiobjective optimization using evolutionary algorithms. Chichester, UK: Wiley. Fukuda, T., & Shibata, T. (1992). Theory and application of neural networks for industrial control systems. IEEE Transactions on Industrial Electronics, 39(6), 472–489. Gaing, Z. L. (2004). A particle swarm optimization approach for optimum design of PID controller in AVR system. IEEE Transactions on Energy Conversion, 19(2), 384–391. Ghoshal, S. P. (2004). Optimizations of PID gains by particle swarm optimizations in fuzzy based automatic generation control. Electric Power Systems Research, 72, 203–212. Kern, S., Müller, S. D., Hansen, N., Büche, D., Ocenasek, J., & Koumoutsakos, P. (2004). Learning probability distributions in continuous evolutionary algorithms – A comparative review. Natural Computation, 3(1), 77–112. Lee, C. C. (1990). Fuzzy logic in control systems: Fuzzy logic controller – Part I and II. IEEE Transactions on Systems Man and Cybernetics, 20(2), 404–435. Mukherjee, V., & Ghoshal, S. P. (2007). Intelligent particle swarm optimized fuzzy PID controller for AVR system. Electric Power Systems Research, 77(12), 1689–1698. Shi, Y., & Eberhart, R. C. (1998). A modified particle swarm optimizer. Proceedings of IEEE International Conference on Evolutionary Computation, 69–73. Anchorage, AK. Storm, R., & Price, K. (1997). Differential evolution – a simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization, 11, 341–359. Wang, J. S., Zhang, Y., & Wang, W. (2006). Optimal design of PI/PD controller for non-minimum phase system. Transactions of the Institute of Measurement and Control, 28(1), 27–35. Wang, Q. G., Zou, B., Lee, T. H., & Qiang, B. (1997). Auto-tuning of multivariable PID controllers from decentralized relay feedback. Automatica, 33(3), 319–330. Zuo, W. (1995). Multivariable adaptive control for a space station using genetic algorithms. IEE Proceedings – Control Theory and Applications, 142(2), 81–87. Evolutionary algorithms based design of multivariable PID controller Introduction PID controller structure PID controller for MIMO system Performance index Evolutionary algorithms Real coded genetic algorithm with SBX crossover Simulated binary crossover Polynomial mutation Modified particle swarm optimization (MPSO) Differential evolution Covariance matrix adaptation evolution strategy (CMAES) CMAES algorithm MIMO system EA implementation of multivariable PID controller Simulation results Parameter tuning Tuning of multivariable PI controller Tuning of multivariable PID controller Conclusions Acknowledgments References