Document downloaded from: This paper must be cited as: The final publication is available at Copyright http://dx.doi.org/10.1016/j.eswa.2012.01.111 http://hdl.handle.net/10251/50835 Elsevier Reynoso Meza, G.; Sanchís Saez, J.; Blasco Ferragud, FX.; Herrero Durá, JM. (2012). Multiobjective evolutionary algorithms for multivariable PI controller design. Expert Systems with Applications. 39(9):7895-7907. doi:10.1016/j.eswa.2012.01.111,. Multiobjective Evolutionary Algorithms for Multivariable PI Controller Design Gilberto Reynoso-Meza∗, Javier Sanchis, Xavier Blasco, Juan Herrero Grupo de Control Predictivo y Optimización Heuŕıstica (CPOH) Instituto Universitario de Automática e Informática Industrial Universitat Politècnica de València Camino de Vera 14, 46022 - Valencia, Spain http://ctl-predictivo.upv.es Abstract A multiobjective optimisation engineering design (MOED) methodology for PI controller tuning in multivariable processes is presented. The MOED pro- cedure is a natural approach for facing multiobjective problems where several requirements and specifications need to be fulfilled. An algorithm based on the differential evolution technique and spherical pruning is used for this purpose. To evaluate the methodology, a multivariable control benchmark is used. The obtained results validate the MOED procedure as a practical and useful technique for parametric controller tuning in multivariable processes. Keywords: multiobjective optimisation, controller tuning, pid tuning, multiobjective evolutionary optimisation, decision making. 1. Introduction PI and PID controllers currently represent a reliable digital control solu- tion because of their simplicity and efficacy [3]. They are often used in indus- trial applications and there is ongoing research on new techniques for robust tuning in single-input single-output (SISO) systems, as well as multiple-input multiple-output (MIMO) systems. MIMO systems are very common in in- ∗Corresponding author. Tel.:+34963877007;fax:+34963879579 Email address: gilreyme@posgrado.upv.es (Gilberto Reynoso-Meza) URL: http://ctl-predictivo.upv.es/ (Gilberto Reynoso-Meza) Preprint submitted to Expert Systems with Applications April 7, 2011 dustrial processes, and their complexity relies on the dynamic interaction between inputs and outputs. New PI-PID controller tuning techniques mainly search for a trade-off so- lution among several control and operational requirements. Some approaches state the design problem as an analytical/numerical optimisation procedure [14, 51, 15, 2, 36], or as an evolutionary optimisation statement [24, 22, 23]. In both cases, a variety of specifications with several requirements and spec- ifications must be faced. Such problems involving multiple objectives are known as multiobjective problems (MOP). In an MOP, the designer (control engineer) has to deal with a list of re- quirements and searches for a solution with a desired trade-off (preferences) among objectives. A traditional approach to handle preferences in an MOP is to translate it into a single-objective problem using weighting factors. More elaborate methods have been developed [30], such as goal programming, lex- icographic methods, physical programming [55], and recently, global physical programming [33, 45]. Multiobjective optimisation (MOO) can handle these issues in a simpler manner because of its simultaneous optimisation approach. In MOO, all of the objectives and constraints are significant from the designer’s point of view. Consequently, each is optimised to obtain a set of optimal non- dominated solutions. In this set of solutions, no solution is better than the others in every objective - but each solution offers different balances between design objectives. As a result, the decision maker (DM) can obtain a better insight into the trade-off for different solutions and can analyse the tendencies. This approach produces more information for selecting the most preferable solution that meets the DM’s preferences. The difficulty involved in the PI-PID tuning process based on optimisa- tion increases when: • MIMO systems are considered instead of SISO systems. • The number of engineering requirements (objectives) increases. • MOO is required instead of single objective optimisation. • Constrained problems are treated instead of unconstrained problems. It is, therefore, worthwhile searching for new algorithms and strategies to tackle constrained MOO for PI-PID tuning in multivariable processes. 2 Therefore, this paper proposes an MOP statement for constrained MIMO PI tuning that demonstrates its viability in an easy and intuitive way. This is fulfilled by defining the MOO statement with well-known performance indexes and a graphical visualisation of the Pareto front. This is a very im- portant issue since the DM requires a useful and interpretable approximation for the decision making stage. The remainder of this paper is organised as follows: in Section 2 a review of MOO is presented; in Section 3 a multiobjective optimisation engineer- ing design (MOED) methodology for multivariable PI controller tuning is explained. In Section 4 the MOED methodology is evaluated in a multivari- able benchmark process. Finally, some concluding remarks are given. 2. Multiobjective optimisation review An MOP, without loss of generality,1 can be stated as follows: min θ∈ℜn J(θ) = [J1(θ), . . . , Jm(θ)] ∈ ℜ m (1) where θ ∈ ℜn is defined as the decision vector and J(θ) as the objective vector (see Figure 1). A unique solution does not generally exist for an MOP because no solution is better than the others for all the objectives. Let ΘP be defined as the Pareto set, or set of solutions of the MOP, and JP be defined as the Pareto front or the projection of ΘP in the objective space. Each point in the Pareto front is said to be a non-dominated solution (see Figure 2). Definition 1. (Dominance relation): given a solution θ1 with objective vec- tor J(θ1) dominates a second solution θ2 with objective vector J(θ2) if and only if: {∀i ∈ [1, 2, . . . m], Ji(θ 1) ≤ Ji(θ 2)} ∧ {∃q ∈ [1, 2, . . . m] : Jq(θ 1) < Jq(θ 2)} which is denoted as θ1 ≺ θ2. 1A maximisation problem can be converted to a minimisation problem taking into account that max Ji(θ) = min(−Ji(θ)) is applied. 3 q 1 J2 q2 q3 J1 q 2 q 3 q 4 q1 E s p a c i o d e v a r i a b l e s d e d e c i s i ó n D E s p a c i o d e f u n c i o n e s o b j e t i v o J( )q 1 J( )q 2 J( )q 3 J( )q 4J(·) Á r e a d e d o m i n a n c i a Figure 1: Pareto set (left) and Pareto front (right). Objective vector J(θ4) is dominated by J(θ2) . Two useful vectors can be defined: the ideal solution Jmin and the nadir solution Jmax: Jideal = Jmin = [ min J(θ)∈J∗ P J1(θ), . . . , min J(θ)∈J∗ P Jm(θ) ] Jnadir = Jmax = [ max J(θ)∈J∗ P J1(θ), . . . , max J(θ)∈J∗ P Jm(θ) ] (2) MOO techniques search for a discrete approximation Θ∗ P of the Pareto set ΘP capable of generating a good quality description J ∗ P of the Pareto front JP (see Figure 3). In this way, the DM has a set of solutions for a given problem and more flexibility for choosing a particular or desired solution. There are several widely used algorithms for calculating this Pareto front approximation (normal boundary intersection method [7, 35], normal constraint method [1, 44, 32, 31], and the successive Pareto front optimisation [43]). Recently, multiobjective evolutionary algorithms (MOEAs) have been used due to their flexibility in dealing with non-convex and highly constrained functions [6, 5]. For this reason, MOEAs are considered in this work. A general framework is required to successfully incorporate the MOO approach into any engineering process. A multiobjective optimisation engi- neering design (MOED) methodology is shown in Figure 4. It consists in 4 Figure 2: Dominance concept. A given objective vector A dominates the objective vectors with a better or equal cost value in all objectives (with, at least, one of them being better). Two important points are defined: the ideal solution and the nadir solution (see Equation 2). 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 J1 J 2 True Pareto front (continous) J P Solutions describing a Pareto front approximation J P * Dominated Solutions (+) Figure 3: Pareto front concept (example of two objectives). Points are a possible Pareto front approximation obtained by a particular optimisation algorithm. 5 three main steps: MOP definition: at this stage the following are defined: the design concept (how to tackle the problem at hand); the engineering requirements (what it is important to optimise); and the constraints (which solutions are not practical/allowed). The design concept implies the existence of a parametric model that defines the parameter values (the decision space) that leads to a particular design alternative and its performance [34]. MOO process: at this stage, the MOO statement, as well as the MOEA, are defined. It is important to select an MOEA that assures reasonable diversity, spread, and convergence to the Pareto front and is an efficient constraint handling mechanism. Decision making stage: finally, with the calculated approximation J∗ P , the DM can analyse the trade-off along the Pareto front. The DM will select the best vector solution according to his/her needs. A reli- able tool or methodology is required for this final step, since it is not a trivial task to perform an analysis on m-dimensional Pareto fronts. Figure 4: Multiobjective optimisation engineering methodology. The MOED methodology for PI tuning the multivariable process is then defined. 6 3. Multiobjective optimisation engineering design applied to mul- tivariable PI controller tuning MIMO systems are common in industry. Their complexity is due to their coupling effects between inputs and outputs. Consider a NxN multivariable process modeled by the following transfer matrix: Gp(s) =    Gp11(s) . . . Gp1N(s) ... ... ... GpN1(s) . . . GpNN(s)    (3) The selected controller design concept must fulfill a set of requirements, in accordance with the given process. Common choices for controlling MIMO system are: decoupled PI-PID controllers [22]; centralised PI-PID controllers [23]; state space feedback techniques [13, 40]; or predictive control [12, 47, 28]. The selection of one technique over another depends on the desired balance between complexity and trade-off between design specifications. 3.1. MOP definition In accordance with Figure 4, the first step is to define the design technique (leading to the decision space definition), the operational constraints, and the objective space (optimisation objectives). Using MOEAs in the MOO process gives a greater flexibility to use any type of parametric controller and define any type of performance objective. In this work, a set of decoupled PI controllers is proposed to tackle the control problem in a MIMO system. PI controllers are a simple but successful solution, and they can be improved with complementary techniques (see [3]). Equation (4) shows the structure of the chosen PI controller: Gc(s) = kc ( 1 + 1 Tis ) E(s) (4) where kc is the proportional gain, Ti the integral time (secs), and E(s) the error signal. The decoupled PI controller Gc(s) design has N SISO controllers: Gc(s) =    Gc1(s) . . . 0 ... ... ... 0 . . . GcN(s)    (5) 7 Therefore, the decision space is defined as: θ = [kc1, Ti1, . . . , kcN, TiN] ∈ ℜ 2N (6) The non-convex optimisation developed by [2] will be used as guideline for the SISO PI controllers. This optimisation procedure is analytical and model oriented and does not require any time domain function computations (simulations). It defines a given value of the maximum sensitivity func- tion as a design constraint Ms = max ∣ ∣ ∣ 1 1+Gc(ω)Gp(ω) ∣ ∣ ∣ and/or the maximum complementary sensitivity function Mp = max ∣ ∣ ∣ Gc(ω) 1+Gc(ω)Gp(ω) ∣ ∣ ∣ . A numerical non-convex optimisation is then used, by increasing as much as possible the integral gain ki = kc/Ti subject to the values of Ms and Mp, to obtain a desired trade-off between load rejection and robustness. The previous tuning procedure can be adapted for MOEAs and defining as engineering control objectives ki, Ms and Mp. Such objectives give the DM some insight regarding the trade-off for robustness, load rejection, and set point response as in [2]. To apply this tuning procedure in a multivariable process, an index to measure the overall MIMO system stability is required. Here, the closed loop log modulus (Lcm) will be used as a robustness indicator. This index leads to the well-known largest log modulus (BLT) tuning criteria for diagonal PID controllers in MIMO processes [29]. The criteria is defined as: Lcm = 20 log ∣ ∣ ∣ ∣ W(s) 1 + W(s) ∣ ∣ ∣ ∣ (7) where W(s) = −1 + det (I + Gp(s)Gc(s)). Therefore, the MOP at hand is to find a trade-off solution θ, that is: J(θ) = [−ki1, Ms1, Mp1, . . . , −kiN, MsN, MpN, Lcm] ∈ ℜ 3N+1 (8) The objective vector as defined by Equation (8) does not guarantee to give the DM a useful Pareto front with a good degree of flexibility to select a reliable and practical solution. It is well-known that certain practical limits to Ms, Mp and Lcm values are needed to guarantee a minimum of stability margin. Therefore, the MOP statement must consider the following practical limits: 8 kc1 + ν1 · kc1/Ti1 ≤ Ku1 ... kcN + νN · kcN/TiN ≤ KuN 1.2 ≤ Ms1,...,N ≤ 2.0 1 ≤ Mp1,...,N ≤ 1.5 0 ≤ Lcm ≤ 2N (9) Where ν is the maximum value between the time delay process and 1. Constraint kc +ν ·kc/Ti ≤ Ku is used to bound the maximum allowed control action effort to the ultimate gain Ku. Constraints 1.2 ≤ Ms and 1 ≤ Mp are used to avoid controllers with a sluggish performance, while constraints Ms ≤ 2.0 and Mp ≤ 1.5 guarantee a minimum of stability margin [2]. The empirical rule of keep Lcm ≤ 2N [29] is accepted. 3.2. The MOO process As constraints are considered in the MOP, a constraint handling mecha- nism is used. According to the practical and empirical limits defined for J∗P by Equation (9), any unfeasible solution is punished. In [8], a penalty func- tion without penalty parameter is proposed. Such penalty function enforces the following criteria: 1. Any feasible solution is preferred to any infeasible solution. 2. Between two feasible solutions, the solution with the better objective function value is preferred. 3. Between two infeasible solutions, the solution with the smaller con- straint violation is preferred. Following these ideas, the MOO vector objective takes the form: min θ∈ℜ2N J(θ) =        J(θ) ∈ ℜ3N+1 if 7 ∑ k=1 φk(θ) = 0 offset + ( 7 ∑ k=1 φk(θ) ) · R ∈ ℜ3N+1 otherwise (10) 9 where: offset = max (Jmax) · R φ1(θ) = max{0, kc1 + ν1kc1 Ti1 − Ku1, . . . , kcN + νNkcN TiN − KuN} φ2(θ) = max{0, 1.2 − Ms1 . . . , 1.2 − MsN} φ3(θ) = max{0, 1.0 − Mp1 . . . , 1.0 − MpN} (11) φ4(θ) = max{0, Ms1 − 2.0, . . . , MsN − 2.0} φ5(θ) = max{0, Mp1 − 1.5, . . . , MpN − 1.5} φ6(θ) = max{0, Lcm − 2N} φ7(θ) = max{0, −Lcm} (12) and R is a vector of 1s with dimension 1x7. An extensive list of MOEAs are available for MOO. Some examples are NSGA-II [10], MOGA [11], ǫv-MOGA [19], paǫ-MyDE [17], sp-MODE [37], among others. They have been used for SISO PID tuning [20, 41, 42, 38, 50], as well for MIMO tuning [53]. In this case, the differential evolution (DE) algorithm [49, 48] was selected as the evolutionary technique. Multiobjective optimisation algorithms based on DE have shown a remarkable performance in a variety of multiobjective optimisation problems [21]. The DE algorithm has two main operators: mutation and crossover. Mutation: At generation k for each target (parent) vector θi|k, a mutant vector vi|k is generated according to Equation (13): v i|k = θ r1|k + F(θ r2|k − θ r3|k) (13) Where r1 6= r2 6= r3 6= i and F is known as the scaling factor. Crossover: For each target vector θi|k and its mutant vector v i|k, a trial (child) vector ui|k = [u i 1|k, u i 2|k, . . . , u i n|k] is created as follows: uij|k = { vij|k if rand(0, 1) ≤ Cr θij|k otherwise (14) where j ∈ 1, 2, 3 . . . n and Cr is the crossover probability rate. 10 For single objective optimisation, a child is selected over its parent (for the next generation) if it has a better cost value. In MOO, a child is selected over his parent if child ≺ parent. The DE by itself is only capable of evolving its population towards the Pareto front, and cannot improve diversity or spread along the Pareto front. As a result, it is necessary to use a mechanism to improve diversity. The use of an external archive is very common and widely accepted in MOEAs. It consists in using an evolutive population plus an external archive where quality solutions are stored and usually used in the evolutionary strat- egy itself. Several techniques have been used to maintain diversity and spread solutions in this archive. In [27], a relaxed form of Pareto dominance, known as ǫ−dominance, is proposed. The main idea is to use an archiving strategy, where instead of using the classical relation for dominance (Definition 1), the following ideas are used: • A solution dominates the solutions that are less fit for all the objectives. • A solution dominates the solutions inside a distance that is less than a parameter ǫ (ǫ-dominance concept). Algorithms based on this approach include the ǫ−MOEA [9], ǫ−MyDE [46], ǫv−MOGA [18, 19], paǫ−MyDE [17], and paǫ−ODEMO [16]. Pruning techniques are commonly based on using some kind of measurement to select individuals in less crowded areas. The crowding distance from [10] or a based measurement are used in algorithms such NSGAII [10] or GDE3 [26]. Nevertheless, it is worthwhile to look for new techniques to deal with any geometrical characteristic of the Pareto front, such as concavity, convexity, disconnected segments, or mixed characteristics. In this work, a pruning technique is employed, based on spherical relations in the objective space [39]. The technique shows a good flexibility in dealing with diverse geometries in m-dimensional Pareto fronts and achieves a well-spread set of solutions. The basic idea of the spherical relations is to analyse the proposed solu- tions in the current Pareto front approximation by using normalised spherical coordinates from a reference solution (see Figure 5a). The spherical pruning works over a non-dominated set of solutions. This can be interpreted as if the designer stands at the ideal solution (or any desired solution) with a given direction in the objective space. The DM will then be searching for the closest non-dominated solution (Figure 5b). 11 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 J 1 J 2 J 3 (a) Spherical relations (b) Set of solutions  → (c) Dominance on  → (d) Spherical pruning on  Figure 5: Spherical relations (up) and spherical pruning (down) for MOO problems (‖ · ‖2 case). Extreme solutions are always preserved. 12 The spherical relations are compatible with any MOO algorithm or evolu- tionary technique and will be merged together with the DE algorithm. Such an approach, presented as sp-MODE in previous works, has been shown to be effective in addressing control engineering problems [37, 42, 40]. 3.3. Decision making stage An m-dimensional Pareto front J∗P is difficult to analyse without an ef- fective visualization tool. If J∗P is not clearly displayed to the DM, it will be complicated to perform a practical analysis in the Pareto front approxima- tion, or select a solution with a desired trade-off. The graphical visualisation is not a trivial task when the number of objectives is more than three. The level diagram (LD) tool [4] is used because it is flexible in performing a useful analysis on the obtained Pareto front J∗P . LDs are based on the classification of the J∗P approximation obtained. Each objective Jq(θ) is normalised with respect to its minimum and maximum values. That is: Ĵq(θ) = Jq(θ) − J min q Jmaxq − J min q , q ∈ [1, . . . , m]. (15) On each normalised objective vector Ĵ(θ) a p-norm ‖x‖p is applied to evaluate the distance to an ideal solution Jideal ≈ Jmin 2. The LD tool displays a two-dimensional graph for every objective and every decision vari- able. The ordered pairs ( Jq(θ), ‖Ĵ(θ)‖p ) in each objective sub-graph and ( θl, ‖Ĵ(θ)‖p ) , l ∈ {1, 2, . . . , n} in each decision variable sub-graph are plot- ted. Therefore, a given solution will have the same y-value in every graphic. This correspondence will help to evaluate general tendencies along the Pareto front and compare solutions in accordance with the selected norm. A deeper explanation of the LD tool capabilities can be found in [4]. A MIMO benchmark will next be considered to validate the MOED for the multivariable PI controller tuning defined in this work. 4. Procedure validation To show the applicability of the MOED proposal for multivariable PI tuning, the well-known distillation column model defined by Wood and Berry will be used [52]: 2Since Jideal is not always available 13 Gp(s) =   Gp11(s) Gp12(s) Gp21(s) Gp22(s)   =   12.8e−s 16.7s+1 −18.9e−3s 21s+1 6.6e−7s 10.9s+1 −19.4e−3s 14.4s+1   (16) As mentioned earlier, any kind of parametric controller can be tuned with the MOED methodology, but for comparison purposes two PI controllers will be used: Gc(s) =     kc1 ( 1 + 1 Ti1s ) 0 0 kc2 ( 1 + 1 Ti2s )     (17) 4.1. Engineering design process Given equations (16) and (17), the MOP at hand is to find a trade-off solution θ = [kc1, Ti1 kc2, Ti2] for the design objectives: J(θ) = [−kc1/Ti1, Ms1, Mp1, −kc2/Ti2, Ms2, Mp2, Lcm] (18) subject to: kc1 + kc1/Ti1 ≤ Ku1 ≈ 2.0 |kc2 + 3kc2/Ti2| ≤ |Ku2| ≈ | − 0.42| 1.2 ≤ Ms1,2 ≤ 2.0 (19) 1 ≤ Mp1,2 ≤ 1.5 0 ≤ Lcm ≤ 4 4.2. Multiobjective optimisation process The MOO objective vector shown in Table 1 is in accordance with Equa- tion (10). The optimisation process is performed with three different MOEAs: • A DE algorithm without archiving strategy (NA); namely, a child will be selected over his parent if child ≺ parent. Parameter values F = 0.5, Cr = 0.8 are used (which are standard parameters in accordance with [48]) and an initial population of 50 random decision vectors. 14 • A DE algorithm with spherical pruning (sp-MODE). Parameter values F = 0.5, Cr = 0.8, a population of 50 solutions, and a spherical grid resolution of 5 are used. • The gamultiobj algorithm provided by MatLab c© is used to calculate a Pareto front for reference. This algorithm uses a controlled elitist genetic algorithm (a variant of NSGA-II [10]). Diversity is maintained by controlling the elite members of the population as the algorithm progresses by using a crowding distance index. Default parameters are used and the BLT solution [29] is used in its initial population. The maximum allowable function evaluations (FEs) for each method is bound to 6000, and 25 independent runs will be evaluated to analyse their performance. Each execution from the sp-MODE and the NA strategy will be compared with the Pareto front J∗P |GA built with the executions of the gamultiobj algorithm. To evaluate the performance of each MOEA, the Iǫ binary indicator [54, 25] is used. The indicator indicates the factor Iǫ(A, B) by which an approx- imation set A is worse than another set B with respect to all the objectives. Using a comparison method (see Table 2) CIǫ,E(A, B) = E(Iǫ(A, B), Iǫ(B, A)) = {false, true} the Eps binary indicator is a compatible and complete oper- ator 3 and this is useful to determine if two Pareto fronts are incomparable, equal, or if one is better than the other [54]. The optimisation experiments were carried on an a standard PC, with a Pentium(R) processor running at 3.40 GHz and 2 GB RAM. The results after 25 independent trials with each proposal are shown in Table 3 (performance indicators) and Table 4 (non-dominated solutions attained). As evidenced by the given results, the sp-MODE algorithm represents a viable approach for generating the Pareto front. The sp-MODE algorithm outperforms the gamultiobj algorithm, since Iǫ(sp − MODE, GA) < 1. Be- sides, the sp-MODE algorithm has a better improvement over J∗P |GA than the NA-strategy (Iǫ(sp − MODE, GA) < Iǫ(NA, GA)). 3Given a binary relation on approximation sets (·), the comparison method is com- patible if CI,E(A, B) → A · B ∨ B · A. However, the comparison method is complete if A · B ∨ A · B → CI,E(A, B). 15 min θ∈ℜ4 J(θ) =        J(θ) ∈ ℜ7 if 7 ∑ k=1 φk(θ) = 0 offset + 7 ∑ k=1 φk(θ) · R ∈ ℜ 7 otherwise J(θ) = [J1(θ), J2(θ), J3(θ), J4(θ), J5(θ), J6(θ), J7(θ)] θ = [kc1, Ti1, kc2, Ti2] J1(θ) = −ki1 = −kc1/Ti1 Kc1 ∈ [0.001, Ku1] J2(θ) = Ms1 = max ∣ ∣ ∣ 1 1+Gc1(ω)Gp11(ω) ∣ ∣ ∣ Kc2 ∈ [Ku2, −0.001] J3(θ) = Mp1 = max ∣ ∣ ∣ Gc1(ω) 1+Gc1(ω)Gp11(ω) ∣ ∣ ∣ Ti1,i2 ∈ [0.001, 40] J4(θ) = −ki2 = −kc2/Ti2 J5(θ) = Ms2 = max ∣ ∣ ∣ 1 1+Gc2(ω)Gp22(ω) ∣ ∣ ∣ J6(θ) = Mp2 = max ∣ ∣ ∣ Gc2(ω) 1+Gc2(ω)Gp22(ω) ∣ ∣ ∣ J7(θ) = Lcm = 20 log ∣ ∣ ∣ W (S) 1+W (S) ∣ ∣ ∣ , W(s) = −1 det (I + Gp(s)Gc(s)) offset = max(Jnadir) ∗ R φ1(θ) = max (0, kc1 + kc1/Ti1 − 2.1, |kc2 + 3kc2/Ti2| − 0.42) φ2(θ) = max (0, 1.2 − Ms1, 1.2 − Ms2) φ3(θ) = max (0, 1.0 − Mp1, 1.0 − Mp2) φ4(θ) = max (0, Ms1 − 2.0, Ms2 − 2.0) φ5(θ) = max (0, Mp1 − 1.5, Mp2 − 1.5) φ6(θ) = max{0, Lcm − 4.0} φ7(θ) = max{0, −Lcm} Table 1: MOO statement for the multivariable PID controller approach. 16 Iǫ(A, B) < 1 → Every J(θ b) ∈ B is strictly dominated by at least one J(θa) ∈ A. Iǫ(A, B) ≤ 1 ∧ Iǫ(B, A) > 1 → Every J(θ b) ∈ B is weakly dominated by at least one J(θa) ∈ A and A 6= B. Iǫ(A, B) ≤ 1 → Every J(θ b) ∈ B is weakly dominated by at least one J(θa) ∈ A. Iǫ(A, B) = 1 ∧ Iǫ(B, A) = 1 → A = B. Iǫ(A, B) > 1 ∧ Iǫ(B, A) > 1 → Neither A weakly dominates B nor B weakly dominates A. Table 2: Comparison methods using the Eps indicator. Eps Indicator Iǫ(NA, GA) Iǫ(sp − MODE, GA) Worst 2.34E-001 1.01E-001 Pareto Best 7.74E-002 4.69E-002 Set Median 1.07E-001 7.55E-003 Reference Mean 1.19E-001 7.53E-002 (676 Sol) Std 3.32E-002 1.44E-002 Table 3: Performance achieved by MOEAs. NA sp-MODE Worst 0 109 Best 48 243 Median 43 153 Mean 4.03E+001 1.56E+002 Std 9.15E+000 3.09E+001 Table 4: Number of solutions achieved by MOEAs. 17 kc1 Ti1 kc2 Ti2 J1(θ) J2(θ) J3(θ) J4(θ) J5(θ) J6(θ) J7(θ) BLT [29] 0.3750 8.2900 -0.0750 23.6000 -0.0452 1.2953 1.1081 -0.0032 1.2513 0.9998 3.8599 WIB [22] 0.8485 326.3462 -0.0132 1.9130 -0.0026 1.6663 1.0178 -0.0069 2.0569 1.7259 0.6024 min ‖Ĵ(θ)‖2 0.4245 15.6135 -0.0397 7.0977 -0.0272 1.3090 1.0014 -0.0056 1.3090 1.1047 1.5054 min ‖Ĵ(θ)‖1 0.3351 34.1079 -0.0476 8.9239 -0.0098 1.2263 1.0000 -0.0053 1.2496 1.0427 0.7488 min ‖Ĵ(θ)‖∞ 0.7415 11.2697 -0.0431 5.4571 -0.0657 1.6220 1.0809 -0.0079 1.5097 1.2914 2.1913 J7(θ) = 2.9922 0.7687 6.9516 -0.0408 5.1598 -0.1106 1.6989 1.2144 -0.0079 1.5316 1.3170 2.9922 J7(θ) = 3.4956 0.8458 12.4453 -0.0858 17.6735 -0.0680 1.7434 1.1414 -0.0049 1.3092 1.0000 3.4956 J7(θ) = 3.995 0.92489 8.7357 -0.0783 5.8147 -0.1059 1.8880 1.2790 -0.0135 1.6731 1.4436 3.9950 Table 5: Controllers selected for further evaluation. 4.3. Decision making stage To validate the MOED method approach as a competitive and practical solution for controller tuning, the Pareto front with 153 solutions (median value in Table 4) is selected and used for controller evaluation.4 In Figure 6 the Pareto set and Pareto front using the LD tool are shown respectively. The controller selection procedure lies on the DM preferences and desired specifications. To illustrate the tradeoff achieved by different solutions, six controllers Gc(s) were selected from the Pareto front for further evaluation (see Table 5). The controllers with the lowest ‖Ĵ(θ)‖1, ‖Ĵ(θ)‖2 and ‖Ĵ(θ)‖∞ norm are selected. An overall trade-off between objectives is expected using these controllers. The remaining controllers are selected according to DM preferences. Let’s assume, for example, that the DM is interested in controllers over the A-line in objective J7(θ) (see Figure 6b) and decides to perform a further analysis on three controller from such a geometric locus. The controller resulting from the BLT tuning [29] (oriented to MIMO- stability using the Ziegler-Nichols procedure), as well as the controller pro- posed in [22] (WIB) that minimises the integral of the absolute error for a specific test, are finally included. These controllers will be tested in the multivariable model in three dif- ferent experiments: 1. Set point change in controlled Variable 1; consequently, the perfor- mance to reject the disturbance in controlled Variable 2 is evaluated. 2. Set point change in controlled Variable 2; consequently, the perfor- mance to reject disturbance in controlled Variable 1 is evaluated. 4The best solution attained could be used for this analysis, but this will not be entirely realistic, since it is not always possible to run an optimisation algorithm several times. 18 0 0.2 0.4 0.6 0.8 1 0.6 0.7 0.8 0.9 1 1.1 θ 1 : Proportional Gain 1 | |J (θ )| | ∞ 0 5 10 15 20 25 30 35 40 0.6 0.7 0.8 0.9 1 1.1 θ 2 : Integral Time 1 (secs) −0.12 −0.1 −0.08 −0.06 −0.04 −0.02 0 0.6 0.7 0.8 0.9 1 1.1 θ 3 : Proportional Gain 2 | |J (θ )| | ∞ 0 5 10 15 20 25 30 35 40 0.6 0.7 0.8 0.9 1 1.1 θ 4 : Integral Time 2 (secs) (a) Θ∗P −0.15 −0.1 −0.05 0 0.6 0.8 1 J 1 : Integral Gain (−−−) | |J (θ )| | ∞ −0.015 −0.01 −0.005 0 0.6 0.8 1 J 4 : Integral Gain (−−−) 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 0.6 0.8 1 J 2 : Sensitivity Function | |J (θ )| | ∞ 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 0.6 0.8 1 J 5 : Sensitivity Function 1 1.1 1.2 1.3 1.4 1.5 0.6 0.8 1 J 3 : Complementary Sensitivity Function | |J (θ )| | ∞ 1 1.1 1.2 1.3 1.4 1.5 0.6 0.8 1 J 6 : Complementary Sensitivity Function 0 0.5 1 1.5 2 2.5 3 3.5 4 0.6 0.8 1 J 7 : Biggest Log Modulus | |J (θ )| | ∞ A (b) J∗P Figure 6: Pareto set and Pareto front used for analysis. ‖Ĵ(θ)‖∞ norm is used. A-line indicates controllers that match a hypothetical preference of the DM. 19 3. Simultaneous set point change in both controlled variables. In all cases, the integral of the absolute error (IAE), the integral of the absolute derivative of the control action (IADU), the settling time (ST) at ±2%, the rise time (RT) from 10% to 90%, the maximum deviation (MD), and the overshoot (OS) will be evaluated. In Table 6 and in Figures 7,8 and 9 the obtained results for each controller are shown. Some expected behaviours are noted: • For controllers in the A-line (see J7 at Figure 6b) the greater the Lcm, the greater the control action and the worse are the trade-offs. That is evident since such controllers are incapable of performing well in all the experiments. Notice how these controllers become more oscillating as J7(θ) increases (Figures 7b, 8b, and 9b). • Controller WIB obtains the best value in IAE for Experiment 3; this was expected since this controller was tuned to minimise IAE for the same experiment. Notice how this outstanding performance has a lower trade-off when single set-point change in controlled variable 1 is applied (Figures 7a and 8a). • Controllers with min ‖Ĵ(θ)‖2, min ‖Ĵ(θ)‖1 and min ‖Ĵ(θ)‖∞ have a balanced trade-off between objectives, achieving good overall perfor- mance (Figures 7b, 8b, and 9b). It is important to remark that there are no bad controllers, just controllers with different trade-offs between objectives. As we can see, performances dif- fer. This analysis could assist in scheduling strategies where more than one controller is used. As a final remark, it can be noticed that operational aspects such saturation, initial states, and operational ranges are not con- sidered. MOEA flexibility allows the use of time function computations to incorporate operational aspects and re-define the MOO statement with more meaningful objectives. 5. Conclusions In this work, an MOED methodology for multivariable PI controller tun- ing has been presented. The obtained results validate the methodology as a practical approach. Thanks to the visualisation capabilities of the LD tool, 20 it is easier to perform the controller selection procedure. As the simulations reveal, the MOO approach is validated as a useful tool for control purposes. With this approach, most of the optimisation procedure uses classical control techniques supported with well-known performance objectives. The continuous use of these objectives by the control engineer community ensures practical bounds and quick interpretations for selecting suitable controllers. The Pareto front enables us to have a better insight into the objective trade- off and how it changes between solutions. The MOP definition for the Wood and Berry distillation column will allow further comparisons of MOEA performance. This MOP provides a useful multiobjective constrained problem for controller tuning in the multivariable process, and will help focus these algorithms for a specific class of engineering design problem. Acknowledgment This work was partially supported by the FPI-2010/19 grant from the Universitat Politècnica de València and project DPI2008-02133/DPI from the Spanish Ministry of Science and Innovation. [1] A. Messac, A. Ismail-Yahaya, C. M., 2003. The normalized normal con- straint method for generating the pareto frontier. Structural and multi- disciplinary optimisation (25), 86 – 98. [2] Astrom, K., Panagopoulos, H., Hagglund, T., 1998. Design of pi con- trollers based on non-convex optimisation. Automatica 34 (5), 585 – 601. [3] Åström, K. J., Hägglund, T., 2005. Advanced PID Control. ISA - The Instrumentation, Systems, and Automation Society, Research Triangle Park, NC 27709. [4] Blasco, X., Herrero, J., Sanchis, J., Mart́ınez, M., 2008. A new graphical visualization of n-dimensional pareto front for decision-making in mul- tiobjective optimisation. Information Sciences 178 (20), 3908 – 3924. [5] Coello, C. A. C., Lamont, G. B., 2004. Applications of Multi-Objective evolutionary algorithms, adavances in natural computation vol. 1 Edi- tion. World Scientific Publishing. 21 IAE IADU ST RT MD OS BLT 4.54E+002 8.71E-001 2.29E+001 3.68E+000 – 10.38% Unit Step WIB 2.48E+003 2.02E+000 +2.00E+002 1.39E+000 – 7.24% Reference min ‖Ĵ(θ)‖2 5.78E+002 8.99E-001 5.81E+001 4.02E+000 – 0.24% Y1 min ‖Ĵ(θ)‖1 1.57E+003 6.86E-001 1.40E+002 5.53E+001 – 0.00% min ‖Ĵ(θ)‖∞ 3.34E+002 1.76E+000 3.27E+001 1.50E+000 – 13.78% Y1 J7(θ) = 2.9922 3.32E+002 1.97E+000 2.18E+001 1.35E+000 – 24.43% J7(θ) = 3.4956 3.18E+002 2.28E+000 2.90E+001 1.27E+000 – 29.11% J7(θ) = 3.995 3.16E+002 2.74E+000 1.99E+001 1.11E+000 – 30.04% BLT 1.65E+003 1.02E-001 1.38E+002 — 6.70E-001 — Unit Step WIB 1.01E+003 7.02E-002 6.78E+001 — 8.07E-001 — Reference min ‖Ĵ(θ)‖2 9.54E+002 5.34E-002 5.12E+001 — 6.48E-001 — Y1 min ‖Ĵ(θ)‖1 9.92E+002 5.43E-002 7.00E+001 — 5.36E-001 — min ‖Ĵ(θ)‖∞ 8.20E+002 6.97E-002 5.93E+001 — 8.47E-001 — Y2 J7(θ) = 2.9922 8.59E+002 7.35E-002 6.35E+001 — 9.27E-001 — J7(θ) = 3.4956 1.10E+003 1.63E-001 8.40E+001 — 8.94E-001 — J7(θ) = 3.995 6.44E+002 1.72E-001 5.27E+001 — 9.79E-001 — BLT 3.38E+002 1.92E-001 4.58E+001 — 1.82E-001 — Unit Step WIB 4.22E+003 1.58E-001 +2.00E+002 — 1.49E-001 — Reference min ‖Ĵ(θ)‖2 5.63E+002 1.53E-001 7.27E+001 — 1.36E-001 — Y2 min ‖Ĵ(θ)‖1 1.56E+003 1.53E-001 1.57E+002 — 1.89E-001 — min ‖Ĵ(θ)‖∞ 2.32E+002 1.54E-001 4.03E+001 — 9.63E-002 — J7(θ) = 2.9922 1.40E+002 1.57E-001 2.86E+001 — 8.68E-002 — Y1 J7(θ) = 3.4956 2.25E+002 2.62E-001 3.85E+001 — 1.39E-001 — J7(θ) = 3.995 1.47E+002 2.42E-001 2.61E+001 — 1.44E-001 — BLT 3.26E+003 1.63E-001 1.73E+002 8.58E+001 — 0.00% Unit Step WIB 1.80E+003 1.24E-001 6.85E+001 2.05E+001 — 6.19% Reference min ‖Ĵ(θ)‖2 1.85E+003 1.04E-001 6.78E+001 3.59E+001 — 0.00% Y2 min ‖Ĵ(θ)‖1 1.94E+003 1.04E-001 9.48E+001 3.72E+001 — 0.00% min ‖Ĵ(θ)‖∞ 1.38E+003 1.10E-001 3.84E+001 2.34E+001 — 1.33% J7(θ) = 2.9922 1.43E+003 1.14E-001 5.91E+001 2.35E+001 — 2.11% Y2 J7(θ) = 3.4956 2.13E+003 1.84E-001 1.11E+002 5.19E+001 — 0.00% J7(θ) = 3.995 9.00E+002 1.80E-001 4.75E+001 6.62E+000 — 2.57% BLT 5.75E+002 1.01E+000 3.23E+001 2.98E+000 – 26.57% Unit Step WIB 2.35E+002 1.95E+000 8.33E+000 1.39E+000 – 8.04% Reference min ‖Ĵ(θ)‖2 5.41E+002 9.57E-001 5.03E+001 3.28E+000 – 13.84% Y1,Y2 min ‖Ĵ(θ)‖1 6.35E+002 7.30E-001 8.34E+001 4.68E+000 – 7.36% min ‖Ĵ(θ)‖∞ 3.74E+002 1.78E+000 2.55E+001 1.50E+000 – 19.39% Y1 J7(θ) = 2.9922 3.67E+002 1.99E+000 1.30E+001 1.35E+000 – 29.46% J7(θ) = 3.4956 3.94E+002 2.33E+000 2.50E+001 1.27E+000 – 26.38% J7(θ) = 3.995 3.86E+002 2.74E+000 2.53E+001 1.11E+000 – 35.27% BLT 1.81E+003 2.37E-001 2.01E+001 4.77E+000 – 29.61% Unit Step WIB 7.97E+002 6.91E-002 1.73E+001 2.23E+000 – 6.91% Reference min ‖Ĵ(θ)‖2 1.10E+003 1.19E-001 1.83E+001 4.59E+000 – 11.94% Y1,Y2 min ‖Ĵ(θ)‖1 1.07E+003 1.23E-001 1.80E+001 5.28E+000 – 12.28% min ‖Ĵ(θ)‖∞ 1.01E+003 1.48E-001 1.72E+001 3.61E+000 – 37.50% Y2 J7(θ) = 2.9922 1.04E+003 1.51E-001 1.67E+001 3.49E+000 – 43.74% J7(θ) = 3.4956 1.42E+003 3.31E-001 1.61E+001 3.89E+000 – 55.55% J7(θ) = 3.995 1.15E+003 3.32E-001 2.97E+001 3.68E+000 – 79.98% Table 6: Controller performance in experimental setup. In bold appears the best value and in italics the worst value in each experiment. 22 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Y 1 0 50 100 150 200 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time U 1 0 50 100 150 200 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Y 2 0 50 100 150 200 0 0.01 0.02 0.03 0.04 0.05 0.06 Time U 2 BLT WIB min 2−norm min 1−norm (a) 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Y 1 0 50 100 150 200 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time U 1 0 50 100 150 200 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Y 2 0 50 100 150 200 0 0.02 0.04 0.06 0.08 0.1 0.12 Time U 2 min inf−norm J7=2.9922 J7=3.4956 J7=3.995 0 5 10 15 20 0.9 1 1.1 1.2 1.3 (b) Figure 7: Performance in Experiment 1. 23 0 100 200 300 400 500 0 0.05 0.1 0.15 0.2 Y 1 0 100 200 300 400 500 −0.15 −0.1 −0.05 0 Time U 1 0 100 200 300 400 500 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Y 2 0 100 200 300 400 500 −0.12 −0.1 −0.08 −0.06 −0.04 −0.02 0 Time U 2 BLT WIB min 2−norm min 1−norm (a) 0 100 200 300 400 500 −0.05 0 0.05 0.1 0.15 Y 1 0 100 200 300 400 500 −0.2 −0.15 −0.1 −0.05 0 Time U 1 0 100 200 300 400 500 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Y 2 0 100 200 300 400 500 −0.14 −0.12 −0.1 −0.08 −0.06 −0.04 Time U 2 min inf−norm J7=2.9922 J7=3.4956 J7=3.995 0 10 20 30 40 0 0.05 0.1 (b) Figure 8: Performance in Experiment 2. 24 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Y 1 0 50 100 150 200 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time U 1 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Y 2 0 50 100 150 200 −0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 Time U 2 BLT WIB min 2−norm min 1−norm 0 5 10 15 20 0 0.5 1 1.5 (a) 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Y 1 0 50 100 150 200 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time U 1 0 50 100 150 200 0 0.5 1 1.5 2 Y 2 0 50 100 150 200 −0.15 −0.1 −0.05 0 0.05 Time U 2 min inf−norm J7=2.9922 J7=3.4956 J7=3.995 0 5 10 15 20 0.9 1 1.1 1.2 1.3 1.4 (b) Figure 9: Performance in Experiment 3. 25 [6] Coello, C. A. C., Veldhuizen, D. V., G.B. Lamont, e., 2002. Evolutionary algorithms for solving multi-objective problems. Kluwer Academic Press. [7] Das, I., Dennis, J., 1998. Normal-boundary intersection: a new method for generating the Pareto surface in non-linear multicriteria optimisation problems. SIAM J. Optim. 8, 631 – 657. [8] Deb, K., 2000. An efficient constraint handling method for genetic al- gorithms. Computer Methods in Applied Mechanics and Engineering 186 (2-4), 311 – 338. [9] Deb, K., Mohan, M., Mishra, S., 2003. Towards a quick computation of well spread Pareto-optimal solutions. In: Fonseca, C., Fleming, P., Zitzler, E., Deb, K., editors., L. T. (Eds.), Evolutionary Multi-criterion optimisation: Second International Conference, EMO. pp. 222 – 236. [10] Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multiobjective genetic algorithm: nsga − ii. IEEE Transactions on Evolutionary Computation Vol. 6 (2), 124 – 141. [11] Fonseca, C., Fleming, P., 1993. Genetic algorithms for multiobjective optimisation: formulation, discussion an generalization. Proceedings of the Fifth International Conference on Genetic Algorithms., 416 – 423. [12] Garćıa-Nieto, S., Mart́ınez, M., Blasco, X., Sanchis, J., 2008. Nonlinear predictive control based on local model networks for air management in diesel engines. Control Engineering Practice 16 (12), 1399 – 1413. [13] Garcia-Alvarado, M., Ruiz-Lopez, I., 2010. A design method for robust and quadratic optimal mimo linear controllers. Chemical Engineering Science 65 (11), 3431 – 3438. [14] Ge, M., Chiu, M.-S., Wang, Q.-G., 2002. Robust pid controller design via lmi approach. Journal of Process Control (12), 3 – 13. [15] Goncalves, E. N., Palhares, R. M., Takahashi, R. H., 2008. A novel approach for h2/h∞ robust pid synthesis for uncertain systems. Journal of process control (18), 19 – 26. [16] Gong, W., Cai, Z., Zhu, L., 2009. An efficient multiobjective differential evolution algorithm for engineering design. Structural and Multidisci- plinary Optimisation 38, 137 – 157, 10.1007/s00158-008-0269-9. 26 [17] Hernández-Dı́az, A. G., Snatana-Quintero, L. V., Coello, C. A. C., Molina, J., 2007. Pareto-adaptive ǫ-dominance. Evolutionary Compu- tation (4), 493 – 517. [18] Herrero, J., Blasco, X., Mart́ınez, M., Ramos, C., 2005. Nonlinear robust identification using multiobjective evolutionary algorithms. In: Mira, J., Álvarez, J. (Eds.), Artificial Intelligence and Knowledge Engineering Ap- plications: A Bioinspired Approach. Vol. LNCS 3562. Springer-Verlag, pp. 231 – 241. [19] Herrero, J., Mart́ınez, M., Sanchis, J., Blasco, X., 2007. Well-distributed Pareto front by using the epsilon-moga evolutionary algorithm. In: et al., F. S. (Ed.), Computational and Ambient Intelligence. Vol. LNCS 4507. Springer-Verlag, pp. 292 – 299. [20] Herreros, A., Baeyens, E., Perán, J. R., 2002. Design of pid-type con- trollers using multiobjective genetic algorithms. ISA Transactions 41 (4), 457 – 472. [21] Huang, V., Qin, A., Deb, K., Zitzler, E., Suganthan, P., Liang, J., Preuss, M., Huband, S., 2007. Problem definitions for perfor- mance assessment on multi-objective optimisation algorithms. Tech. rep., Nanyang Technological University, Singapore. [22] Iruthayarajan, M. W., Baskar, S., 2009. Evolutionary algorithms based design of multivariable pid controller. Expert Systems with Applications 36 (5), 9159 – 9167. [23] Iruthayarajan, M. W., Bskar, S., 2010. Covariance matrix adaptation evolution strategy based design of centralized pid controller. Expert Sys- tems with Applications 37 (8), 5775 – 5781. [24] Kim, T.-H., Maruta, I., Sugie, T., 2008. Robust pid controller tun- ing based on the constrained particle swarm optimisation. Automatica 44 (4), 1104 – 1110. [25] Knowles, J., Thiele, L., Zitzler., E., Feb. 2006. A tutorial on the per- formance assessment of stochastic multiobjective optimizers. Tech. Rep. TIK report No. 214, Computer Engineering and Networks Laboratory. ETH Zurich. 27 [26] Kukkonen, S., Lampinen., J., 2005. Gde3: The third evolution step on generalized differential evolution. Vol. 1. IEEE Congress on Evolutionary Computation (CEC 2005), pp. 443 – 450. [27] Laumanns, M., Thiele, L., Deb, K., Zitzler, E., 2002. Combining conver- gence and diversity in evolutionary multiobjective optimisation. Evolu- tionary Computation (3), 263 – 282. [28] Lauŕı, D., Salcedo, J., Garćıa-Nieto, S., Mart́ınez, M., 2010. Model pre- dictive control relevant identification: multiple input multiple output against multiple input single output. IET Control Theory & Applica- tions 4 (9), 1756–1766. [29] Luyben, W. L., 1986. Simple method for tuning siso controllers in mul- tivariable systems. Industrial and Engineering Chemistry Process De- sign (25), 654 – 660. [30] Marler, R., Arora, J., 2004. Survey of multi-objective optimisation meth- ods for engineering. Structural and Multidisciplinary Optimisation (26), 369 – 395. [31] Mart́ınez, M., Herrero, J., Sanchis, J., Blasco, X., Garćıa-Nieto, S., 2009. Applied pareto multi-objective optimisation by stochastic solvers. Engineering Applications of Artificial Intelligence 22, 455 – 465. [32] Mart́ınez, M., Nieto, S. G., Sanchis, J., Blasco, X., 2009. Genetic al- gorithms optimisation for normalized normal constraint method under Pareto construction. Advances in Engineering Software 40, 260 – 267. [33] Mart́ınez, M., Sanchis, J., Blasco, X., 2006. Multiobjective controller de- sign handling human preferences. Engineering Applications of Artificial Intelligence 19, 927 – 938. [34] Mattson, C. A., Messac, A., 2005. Pareto frontier based concept selection under uncertainty, with visualization. Optimisation and Engineering 6, 85 – 115, 10.1023/B:OPTE.0000048538.35456.45. [35] Miettinen, K. M., 1998. Nonlinear multiobjective optimisation. Kluwer Academic Publishers. 28 [36] Panagopoulos, H., Astrom, K., Hagglund, T., Jan. 2002. Design of pid controllers based on constrained optimisation. Control Theory and Ap- plications, IEE Proceedings - 149 (1), 32 – 40. [37] Reynoso-Meza, G., 2009. Design, coding and implementation of a mul- tiobjective optimisation algorithm based on differential evolution with spherical pruning: applications for system identification and controller tuning. Master’s thesis, Universidad Politécnica de Valencia. URL http://personales.alumno.upv.es/gilreyme [38] Reynoso-Meza, G., Blasco, X., Sanchis, J., 2009. Diseño multiobjetivo de controladores pid para el benchmark de control 2008-2009. Revista Iberoamericana de Automática e Informática Industrial 6 (4), 93 – 103. [39] Reynoso-Meza, G., Blasco, X., Sanchis, J., Herrero, J., (Submitted). Spherical relations to achieve diversity in multiobjective evolutionary algorithms. Evolutionary Computation. [40] Reynoso-Meza, G., Garćıa-Nieto, S., Sanchis, J., Blasco, X., (Submit- ted). Controller tuning using multiobjective optimisation algorithms: a global tuning framework. Control Engineering Practice. [41] Reynoso-Meza, G., Sanchis, J., Blasco, X., 30 November - 2 Decem- ber 2009. Multiobjective design of a digital controller for the throttle control benchmark. In: Proceedings of the IFAC Control Workshop on Engine and Powertrain Control, Simulation and Modeling. IFP Rueil- Malmaison. [42] Reynoso-Meza, G., Sanchis, J., Blasco, X., Mart́ınez, M., 2010. Mul- tiobjective design of continuous controllers using differential evolution and spherical pruning. In: Chio, C. D., Cagnoni, S., Cotta, C., Eber, M., Ekárt, A., I.Esparcia-Alcaráz, A., Goh, C.-K., J.Merelo, J., Neri, F., Preuss, M., Togelius, J., N.Yannakakis, G. (Eds.), Applications of Evolutionary Computation, Part I. Vol. LNCS 6024. Springer-Verlag, pp. 532 – 541. [43] Ruzika1, S., Wiecek, M., 2009. Successive approach to compute the bounded Pareto front of practical multiobjective optimisation problems. SIAM J. Optim. 20, 915 – 934. 29 [44] Sanchis, J., Mart́ınez, M., Blasco, X., Salcedo, J. V., 2008. A new per- spective on multiobjective optimisation by enhanced normalized normal constraint method. Structural and Multidisciplinary Optimisation (36), 537 – 546. [45] Sanchis, J., Mart́ınez, M. A., Blasco, X., Reynoso-Meza, G., 2010. Modelling preferences in multiobjective engineering design. En- gineering Applications of Artificial Intelligence, Article in press. doi:10.1016/j.engappai.2010.07.005. [46] Santana-Quintero, L., Coello, C. C., 2003. An algorithm based on difer- ential evolution for multiobjective problems. In: Dagli, C., Buczak, A., Enke, D., Embrechts, M., editors, O. E. (Eds.), Smart systems engineer- ing design: neural networks, evolutionary programming and artificial Life. Vol. 15. pp. 221 – 220. [47] Song, Z., Kusiak, A., Feb. 2009. Optimisation of temporal processes: A model predictive control approach. Evolutionary Computation, IEEE Transactions on 13 (1), 169 –179. [48] Storn, R., 2008. Sci: Differential evolution research: Trends and open questions. Vol. LNCS 143. Springer, Heidelberg, pp. 1 – 31. [49] Storn, R., Price, K., 1997. Differential evolution: A simple and effi- cient heuristic for global optimisation over continuous spaces. Journal of Global Optimisation 11, 341 – 359. [50] Tavakoli, S., Griffin, I., Fleming, P. J., September 2007. Multi-objective optimisation approach to the pi tuning problem. In: Proceedings of the IEEE Congress on Evolutionary Computation (CEC2007). pp. 3165 – 3171. [51] Toscano, R., 2005. A simple robust pi/pid controller design via numerical optimisation approach. Journal of Process Control (15), 81 – 88. [52] Wood, R. K., Berry, M. W., 1973. Terminal composition control of a binary distillation column. Chemical Engineering Science 28 (9), 1707 – 1717. 30 [53] Xue, Y., Li, D., Gao, F., 2010. Multi-objective optimisation and selec- tion for the pi control of alstom gasifier problem. Control Engineering Practice 18 (1), 67 – 76. [54] Zitzler, E., Thiele, L., Laumanns, M., Fonseca, C., da Fonseca, V., 2003. Performance assessment of multiobjective optimizers: an analysis and review. Evolutionary Computation, IEEE Transactions on 7 (2), 117 – 132. [55] A. Messac, S.M. Gupta, and B. Akbulut. Linear physical programming: A new approach to multiple objective optimization. Trans. on Opera- tional Research, 8:39–49, 1996. 31