Robust Multi-objective Optimization Applied to Engineering Systems Design 1

In engineering systems design, theoretical deterministic solutions can be hardly applied directly to real-world scenarios. Basically, this is due to manufacturing limitations and environmental conditions under which the real system will operate. Therefore, a small variation in the design variables vector can result in a meaningful change on the theoretical optimal design as represented by the minimization of the corresponding vector of objective functions. In this context, it is important to develop methodologies that are able to produce solutions (even suboptimal) that are less sensitive to perturbations in the design variable vector and, consequently, leading to a robust optimal design. In this contribution, first the proposed approach is tested on various mathematical functions. Then, the methodology is applied to the design of two representative engineering systems through multi-objective optimization using the Firefly Colony Algorithm in association with the Effective Mean Concept is presented.The results obtained demonstrate that the design strategy conveyed represents an interesting alternative approach to obtain robust design for a number of engineering applications.

Robust Multi-objective Optimization Applied to Engineering Systems Design

INTRODUCTION
Frequently, during engineering system design, the model, the design variable vector, andthe parameter vector are consideredfree of errors, i.e., they do not contain uncertainties.However, more realistically, small variations in the design variable vector can cause significant variations in the vector of Latin American Journal of Solids and Structures 13 (2016) 1802-1822 objective functions.As mentioned by Ritto et al. (2008), the process of modeling engineering systems introduces two types of uncertainties: i) uncertainties related to the parameters of the model, such as geometrical and constitutive parameters (data uncertainties), and ii) uncertainties due to the proposed model.In this case, regarding the existence of uncertainties in real-world engineering systems, some aspects should be highlighted (Leidemer, 2009):  Even if the real optimum is found, it will never be possible to implement it in practice, because there are uncertainties associated to the manufacture process or even due to the demand of a high degree of accuracy, which can be extremely difficult to achieve or economically unaffordable;  The formulation of optimization problems is inherently static, but the reality is essentially dynamic.In this context, the design problems can be associated with parameter vector fluctuation, such as temperature, wind speed, humidity.There might also be waste of some components that imply alteration in system behavior.
Consequently, the system to be optimized can be very sensitive to small changesin the design variable vector, and thus, small variations in this vector can cause significant changes in the vector of objective functions (Leidemer, 2009).In this context, it is important to determine a methodology that produces solutions less sensitive to small variations in the design variables vector.Solutions with this characteristic are called robust solutions and the procedure to find these solutions is named Robust Optimization (Taguchi, 1984).
Real-world problems involve the simultaneous optimization of two or more (often conflicting) objectives, known as multi-objective optimization problem (MOOP).The solution of such problems is different from that of a single-objective optimization problem, i.e., multi-objective optimization problems normally have not only one but a set of solutions, which may all be equally good (Deb, 2001).
In the literature, several methods for solving MOOP can be found.These methods follow a preference-based approach, in which a relative preference vector is used to scalarize multiple objectives.Since classical searching and optimization methods use a point-by-point approach, through which the solution is successively modified, the outcome of this classical optimization method is a single optimized solution.However, Evolutionary Algorithms (EA) can find multiple optimal solutions in one single simulation run due to their population-based search approach.Thus, EA are ideally suited for multi-objective optimization problems.A detailed account of multi-objective optimization using EA and applications using genetic algorithms are widely found in the literature (Deb, 2001).
Other authors have also contributed to evolutionary algorithms.In this sense we can cite the development of evolutionary optimization methodologies based on genetic algorithms (Castro, 2001;Deb, 2001;Djamaluddin et al. 2015); fish swarm (Lobato and Steffen Jr.); bee colony (Ghashochi-Bargh and Sadr, 2014) and firefly colony (Yang, 2008;Yang, 2009;Pfeifer and Lobato, 2010;Lobato et al. 2011), to cite few.The Firefly Colony Algorithm (FCA)is based on strategies that seek to mimic the behavior of fireflies observed in the nature to update a population of candidates to solve optimization problems.These systems have the capacity to notice and modify their environment in order to seek for diversity and convergence.In addition, this capacity makes possible the communi-Latin American Journal of Solids and Structures 13 (2016) 1802-1822 cation among the individuals of the population.Each individual is able to capture the changes in the environment as generated by local interactions (Yang, 2008).
Therefore, in this contribution a systematic methodology to solve robust multi-objective optimization problems is proposed.This optimization strategy is based on the Effective Mean Concept -EMC (DebandGupta, 2006) in association with the Multi-objective Optimization Firefly Colony Algorithm -MOFCA (Lobato et al., 2011).The organization of this article is as follows.Section 2 introduces briefly the EMC, the MOFCA and robust multi-objective optimization.For illustration purposes, section 3 presents applications involving mathematical functions and two engineering case studies.Finally, the conclusions are presented in the last section.

MULTI-OBJECTIVE OPTIMIZATION FIREFLY COLONY ALGORITHM
In this section, the mean effective concept and the robust multi-objective optimization problem will be defined.In addition, the subroutine of the MOFCA algorithm will be briefly discussed.

Robust Multi-Objective Optimization
Traditionally, during the solution of optimization problems it is commonly considered that mathematical models, variables, and parameters are sufficiently reliable, i.e., there are no errors of modeling and estimation (DebandGupta, 2006).However, as previously mentioned, systems to be optimized are generally sensitive to small changes in the design variables leading to significant changes in the vector of objective functions.In this case, to minimize this effect in the solution of optimization problems, the concept of robust optimization should be used.Robust Optimization is defined as an approach that produces a solution that is not significantly sensitive to small changes in the design variables (Taguchi, 1984).In this context, robustness characterizes an important design level to be achieved, when exposed to a given condition of uncertainty.In addition, this approach is used for modeling optimization problems under uncertainty, where the modeler aims at finding decisions that are optimal for the worst-case realization of the uncertainties within a given set of values (Taguchi, 1984).
The introduction of robustness in the mono and multi-objective contexts requires the consideration of new restrictions and/or new objectives (relationshipbetween the mean and the standard deviation of the vector of objective functions) and probability distribution functions for the design variables and/or objectives.As mentioned by Ritto et al. (2008) and also considered by Soize (2005), Paenk et al. (2006), and Sampaio and Soize (2007), probability tools are required to model the uncertainties, i.e., random variables are associated with the uncertain parameters or matrices; then, probability density functions are constructed.
As an alternative to these classical formulations, Deb and Gupta (2006) extended the Effective Mean Concept originally proposed for mono objective problems to the multi-objective context.In this approach, no additional restriction is inserted into the original problem.Thus, the problem is rewritten as a mean vector of original objectives.The authors applied the present methodology to constrained and unconstrained test problems having two and three objectives.In this case, simulation results were presented by using an evolutionary multi-objective optimization algorithm as applied to an engineering design problem.Structures 13 (2016) 1802-1822 In the present contribution, a robustness measure, which does not use the hypothesis on uncertainty parameters,will be used.The robustness measure and the solution of a robust multi-objective optimization problem are defined as (DebandGupta, 2006):

Latin American Journal of Solids and
an integrable function.The EMCof f , in relation to a neighborhood d of variable x , is a function ( , ) eff f x d , given by: where n is the number of design variables.

Definition -A point *
x is defined as robust multi-objective solution if it is a Pareto-optimal solution for the following multi-objective optimization problem, defined byconsidering the neighborhood d of the variable x :   1 2 min ( , ), ( , ), ... , ( , )   1 subject to ( , ) 0, 1,..., ; where the symbol ( ) B x d indicates the hyper-volume of the neighborhood and m is the number of objectives.
To estimate the integral defined by Equation( 1), a sample is created randomly by using the Latin Hypercube method (Viana, 2008).Moreover, this procedure increases the computational cost due to additional objective function evaluations that are necessary to solve the optimization problem (DebandGupta, 2006).The Equations (2) and (3) considersan optimization problem with m objectives and k inequality constraints.

Firefly Colony Algorithm
The FCA is based on the interesting characteristics of the fireflies' bioluminescence.The fireflies are insects notorious for their light emission.Although biology does not have a complete knowledge to determine all the utilities that the firefly luminescence can bring to, at least three functions have been identified (Yang, 2008): i) as a communication tool and appeal to potential partners in the reproduction, ii) as a bait to lure prey for the firefly, iii) as a warning mechanism for potential predators reminding them that fireflies have a bitter taste.
In the FCA, the attractiveness between two fireflies is determined by thelight intensity emitted and by the distance between them.The intensity is a function of the objective function (Yang, 2008).An effective exploration of the design space is obtained by the movement of fireflies according to the attractiveness of other swarm members with higher intensity of light emitted and a random step vector to avoid a premature convergence.Thus, at the k-th iteration, the movement of the i-th firefly towards the firefly that is more attractive is defined by the following equation (Yang, 2008): Latin American Journal of Solids and Structures 13 (2016) 1802-1822 where r is the distance between two fireflies; w and g are parameters defined as the maximum at- tractiveness (when r =0) and the coefficient of light absorption, respectively.The second factor of the right hand side of Equation ( 4) adds the attractiveness between the fireflies, and the third one, regulated by k , adds randomness to the process.

Multi-Objective Optimization Firefly Colony Algorithm
Due to the success obtained by the FCA in different science and engineering applications (Apostolopoulos and Vlachos, 2011;Lukasik and Zak, 2009;Yang, 2009;Pfeifer and Lobato, 2010), Lobato et al. (2011) proposed the Multi-objective Optimization Firefly Colony Algorithm (MOFCA).This approach is based on the classical FCA associated with Fast Non-Dominated Sorting and has the following structure:  An initial population with Nfirefly fireflies is randomly generated;  All dominated solutions are removed from the population through the operator Fast Non-Dominated Sorting (Deb et al., 2000).In this way, the population is sorted into nondominated fronts (sets of vectors that are non-dominated with respect to each other);  Then, FCA is applied to generate the new population of fireflies (potential candidates to solve the multi-objective optimization problem);  If the number of individuals of the population is larger than the number early defined by the user, it is truncated according to a criterion named the Crowding Distance (Deb et al., 2000).
These steps are repeated until a given stopping criterion is reached.The operators used in the MOFCA are described below.

Fast Non-Dominated Sorting
Fast Non-Dominated Sorting operator was proposed by Deb et al. (2000) in order to sort a population of size N according to the level of non-domination.Each solution must be compared with every other solution in the population to find if it is dominated.At this point, all individuals in the first non-dominated front are found.In order to find the individuals in the next front, the solutions of the first front are temporarily discounted and the above procedure is repeated.The procedure is repeated to find the subsequent fronts.

Crowding Distance Operator
This operator describes the density of solutions surrounding a vector.In order to compute the Crowding Distance for a set of population members, the vectors are sorted according to their objective function values for each objective function.An infinite Crowding Distance (or an arbitrarily large number for practical purposes) is assigned to the vectors with the smallest or largest values.For all other vectors, the Crowding Distance ( i x dist ) is calculated according to (Deb et al., 2000): Latin American Journal of Solids and Structures 13 (2016) 1802-1822 where j f corresponds to the j-th objective function and m is the number of objectives.This operator is important to avoid many points close in Pareto's Curve and to promote the diversity in terms of the space of objectives (Deb et al, 2000).

Treatment of Constraints
In this work, the treatment of constraints is performed through the Static Penalization Method, proposed by Castro (2001).This approach consists in assigning limit values to each objective to play the role of penalization parameters.According to Castro (2001), it is guaranteed that any nondominated solution dominate any solution that violates at least one of the constraints.In the same way, any solution that violates only one constraint will dominate any solution that presents two constraint violations, and so on.For a constrained problem the vector containing the objective functions to be accounted for, is given by: where ( ) f x is the vector of objective functions, p r is the vector of penalty parameters that depends on the type of problem considered, and viol n is the number of violated constraints.More details about this optimization strategy can be found in Lobato et al. (2011).
In the present paper, the methodology proposed by Deband Gupta ( 2006) considering the EM-Cis used to insert robustness in the MOFCA.In this case, the original objective function vector is transformed by considering Eq. ( 2).In addition, to evaluate Eq. ( 1) the Latin Hypercube method is used.It is important to observe that in real-world engineering problems, it is not possible to perform the analytical calculation of effective mean.Thus, the integration is performed numerically using the Trapezoidal Method.
The user has to inform to the routine the vector of objective functions, the vector of constraints, the design space, the MOFCA parameters, the perturbation i d added to the design variable vector, and the sample size sample N .To compute the Effective Mean, the user has to inform the value of i d for each design variable i x , which will be calculated by the Latin Hypercube Method.

APPLICATIONS
For illustration purposes, in this section, applications involving mathematical functions and the design of two representative engineering systems (the optimum design of an industrial robot and the optimization of a flexible rotor)are presented.
The MOFCA parameters considered for all the applicationsare the following (Lobato et al., 2011):  Absorption ( w )and attractiveness coefficients ( g )=1;  Sample size to estimate the Effective Mean ( sample N )= 100.Regarding the mathematical problems, ten simulations for different initial seeds of the random number generator rand aretaken into account to evaluate the convergence and diversity of the solutions obtained by using the MOFCA.To compute the convergence metric (CM ) and diversity metric ( DM )for the cases tests where the analytic solution is known, the following equations are used (Deb, 2001): where i d is the smaller Euclidean distance between the firefly i x and the analytic Pareto's Curve,
ThePareto global and local curves correspond to i x =0 (i=3, 4 and 5) and are obtained from the global minimum and local minimum of ( ) 2 h x , respectively.According toDeb and Gupta (2006), the global and local curves are given by: Latin American Journal of Solids and Structures 13 (2016) 1802-1822 The objectives obtained through the calculation of the EMC are given by the following expressions: where Figure 1 presents the nominal Pareto's Curve (global and local), i.e., without robustness,and the solutions obtained by using MOFCA considering the problem MP1.Note that MOFCA was able to find the global Pareto's Curve.Table 1 presents the metrics computed in each run for the problem MP1.Note thatthe MOFCA converges to the global optimum, except for the run with initial seed equal to zero.In addition, the solutions found by using MOFCA are well distributedon the Pareto's Curve, i.e., there is a good diversity of solutions.
In order to evaluatethe effect of d on the curve, the following perturbation vector   In Table 2 it is possible to observe that MOFCA converges to the global optimum for all the simulations performed.

Mathematical Problem MP2
Mathematically, the problem MP2 can be formulated as follows (Deband Gupta, 2006): where In this problem, the nominal Pareto's Curve is obtained by the activation of the constraint equations, i.e., when 2 1 0.1 0.2 x x = -and considering i x =0 for i=2, 3, 4 and 5.As observed in the previous case, the objectives obtained by EMC are expressed by: where Latin American Journal of Solids and Structures 13 (2016) 1802-1822 ( )

Industrial Robot
This real-world engineering application was first proposed by Eschenauer et al. (1990) and considers the optimization of an industrial robot, in which a hydraulic spring mechanism is used for the balancing of the armAB (see Figure 4).
The design variables associated with the balancing of the robot's arm are the following: localization of spring mechanism in relation to base ( 1 X and 2 X , in meters), pressure of hydraulic accumulator ( 3 X , in MPa), and cylinder diameter ( 4 X , in meters).The robot's performance is a function of the static torque in steering the motor shaft, the work executed by the motor, and the static force in the joint A. The model of the industrial robot presented in Fig. 1 was obtainedbyupdatingof following meta-model: where i b are the approximation coefficients (i = 1, …, 14), i x (i = 1, …, 4) are the design variables vector, and is the vector of response of the model, representing measures to evaluate the performance of the industrial robot.For this purpose, the design vector variables were coded as (Eschenauer et al., 1990): According to Eschenauer et al. (1990), a satisfactory performance of the robot depends on the following criteria: static torque on the shaft of the driving motor ( 1 y ), work performed by the driving motor ( 2 y ), and the static force at the joint A ( 3 y ).
Table 6 presents the experiments considered by Eschenauer et al. (1990), which are now used in this application.
Table 7 presents the meta-models obtained using the Sequential Quadratic Programming Method.

Robust Optimum Design
The formulation for the multi-objective optimization problem is given by (Eschenauer et al., 1990): 4  1   2 1 2 3 4  2   3 1 2 3 4    It can be observed that for the solutions in the interval 9500 2 f £ £ 10000 there is a high sensitivity of the objectives 1 f and 3 f , thus characterizing a zone of small robustness.On the other hand, from the identification of regions with small robustness, it is also possible to identify zones with low sensitivity of the objectives, i.e., the robustness zones.

Flexible Rotor
Rotating machines are used in a wide range of applications regarding different engineering fields, such as automotive, aerospace and power generation industries.The generator units of hydraulic power plants, jet engines, steam turbines, are some examples of applications (Tsuzuki, 2012;Cavalini Jr, 2013).Therefore, the identification of optimal design and operating conditions for such kind of machines is an interesting engineering challenge (Quitzrau, 2002).Assis(1999) applied optimization techniques to identify unknown parameters of a vertical rotating machine.A similar procedure was adopted by Steffen Jr et al. (1999) by using multi-objective optimization techniques.He et al. (2001) applied the Genetic Algorithm associated to a finite element model to detect cracks in a rotating shaft.Bueno (2007) used optimization techniques to obtain the optimal localization of piezoelectric sensors and actuators for active vibration control.Cavalini Jr(2013) identified incipient cracks in a rotating shaft by using a nonlinear phenomenon and the Differential Evolution Algorithm.

Robust Optimum Design
It is well known that critical speeds in rotating machines are associated with high vibration amplitudes, which can lead to irreversible damages on the system.Thus, it is important keep the rotor operating as always as possible from those rotation speeds.Regarding the associated minimization process, the optimal (and robust) solution should have a low sensibility to small variations on the design variables (for example, changes due to wear).In this context, the maximization between the first and second forward critical speeds of the rotating machine presented in Fig. 8    The robust optimization process was formulated to deal with 4 design variables, namely the external radius of the finite elements confined between nodes 1 to 5; 5 to 11; 11 to 15 and 15 to 18 (see Fig. 8).Table 9 shows the design space used in the fitting of the finite element model.In order to avoid assembly problems, the external radius of the shaft associated with the nodes 5 to 11 is given by the sum of the radius determined for the first shaft section (nodes 1 to 5) and the ones for the second section (nodes 5 to 11).Similar procedure was adopted for the sections associated with nodes 11 to 15 and 15 to 18.The robustness evaluation of the solution was performed considering the deviation parameters 1 d = 0.0005 m, 2 d = 0.00075 m and 3 d = 0.001 m, separately, which represent perturbations of 8.33%, 12.5%, and 16.67%, respectively, on the nominal radius of the shaft.
Figure 9 shows the deterministic and robust Pareto's Curve obtained by using the MOFCA algorithm.In this case, both objectives are sensible to small perturbations on the design variables.The first critical speed presents a variation of 2000 rev/min approximately, while the third one reaches 4800 rev/min (nominal value of 2940 rev/min).The robust region for the optimum design of the rotating machine is verified.Note that the third critical speed increases according to d (per- turbation parameter) for the case in which a single first critical speed is considered.Table 10 presents the optimal solutions determined from the deterministic and robust process considering the highly fireflies of Figure 9. 1 to 5 7.76 10 -3 6.01 10 -3 6.10 10 -3 7.89 10 -3 5 to 11 3.73 10 -3 3.47 10 -3 4.05 10 -5 2.67 10 -4 11 to 15 5.67 10 -4 4.03 10 -4 5.34 10 -4 4.20 10 -4 15 to 18 5.27 10 -3 4.90 10 -4 5.62 10 -3 4.60 10 Figures 10 and 11 shows the Campbell diagrams obtained from the nominal radius of the shaft and the one obtained at the end of the optimization process considering the robust ( 2 r ) solution, respectively.As expected, the difference between the first and third critical speeds was maximized by the associated optimization process.
Latin American Journal of Solids and Structures 13 (2016) 1802-1822 text of robust design, this characterization is of great importance, since regions with greater sensitivity have to be crowded out in relation to those with lower sensitivity.
Considering the results obtained, it is important to mention that the robustness parameter influences the robust Pareto's Curve.
Since a systematic study introducing robustness in multi-objective optimization problems (DebandGupta, 2006;Paenk et al., 2006) is not easily available, the problems studied may serve as comparison for future evaluations of other methodologies for robust multi-objective optimization.Regarding optimal robust design, the determination of robustness regions may represent a criterion for the choice of a specific point of the Pareto's Curve for a possible practical implementation.However, it is important to observe that the main disadvantage of this approach is the increase of the number of objective function evaluations, which are necessary to evaluate the integral considered in the Effective Mean Concept, independently from the optimization strategy considered.
It is worth mentioning that the operators used in the MOFCA algorithm were not developed in this paper.In addition, the authors of the present contribution did not propose the Mean Effective Concept originally.However, the use of Mean Effective Concept associated to Firefly Colony Algorithm represents a new approach, considering the literature associated with this work.The development of a multi-objective algorithm for the robust optimization using the FCA is considered as being the main contribution of the present contribution.
Further works will be dedicated to approaches related to dynamically updating the parameters and mutation strategies of the Multi-objective Optimization Firefly Colony Algorithmtogether with its parallelization to reduce the number of objective function evaluations.
f d and l d represents the Euclidean distance between the extreme solutions of the analytic Pareto's Curve and the non-dominated solutions obtained through the MOFCA, and m d is the arithmetic average of i d .In this case, the analytic Pareto's Curve was obtained by using 1000 points equally spaced in order to compute CM and DM .

Figure 1 :
Figure 1: Nominal Pareto's Curves (global e local) and the solution obtained by MOFCAfor the problem MP1.
d d 2d 2d 2d ] was considered(Deb and Gupta, 2006).Figure2(a) presents the Pareto's Curves (global, local, nominal and robust) obtained analytically through the Equations 13 and 14 for different values of d .Figure 2(b) presents the robust Pareto's Curve obtained by using MOFCA considering the problem MP1.It is observed that the convergence to the global solution is reached for the consideredvalues of d .In addition, the solutions are satisfactory as observed by different values of d for each robust Pareto's Curve.

Figure 2 :
Figure 2: Influence of parameter d on the robust Pareto's Curves for the problemMP1.

Figure 3
Figure 3(a) presents the nominal Pareto's Curve without and with the constraint and the solution obtained by using the MOFCA.It is possible to observe that satisfactory resultsare obtained for both convergence and diversity.

Figure 3 :
Figure 3: Nominal Pareto's Curve obtained by MOFCA for the problem MP2 considering rp equal to 10 6 .
[ d d d d ]  was considered to evaluate the robustness of the design, in which d is equal to 0.1 and 0.2.These values represent a disturbance of 5% and 10%, respectively, in the amplitude of the domain of definition of the design variables.Figures 5 to 7 present the vector of objective functions, where the comparison between the nominal and robust cases for d equal to 0.1 and 0.2 is performed.
is considered.A finite element model with 18 elements is used to represent the rotating machine.It is composed of a Latin American Journal of Solidsand Structures 13 (2016) 1802-1822flexible steel shaft with 552 mm length and 12 mm of diameter (nominal diameter; E = 210 GPa, r = 7800 kg/m 3 , u = 0.3), three rigid discs of steel (see Table8; r = 7800 kg/m 3 ), and two identi- cal ball bearings ( m = 10 kg, .Considering the nominal parameters, the first two forward critical speeds of the rotor are approximately: 2620 rev/min and 2940 rev/min.

Figure 8 :
Figure 8: Finite element model for the studied rotor.

Figure 9 :
Figure 9: Deterministic and robust Pareto's Curve obtained by using the MOFCA algorithm.

Table 1 :
Results of the CM and DM metrics for the problem MP1.

Table 2 :
Results of the CM and DM metrics obtained for the robust problem MP1.

Table 3
presents the CM and DM metrics for the nominalproblemMP2.

Table 3 :
Results of the CM and DMmetrics for the problem MP2.

Table 4 :
Table4presents the CM and DM metrics for different values of d .It is possible to observe that the obtained values for the metrics are satisfactory, i.e., the CM tends to small values and the DM metric values present good distribution (see Figure3(b)).Results of the CM and DM metrics obtained for the robust problemMP2.

Table 5 :
Table 5 presents the robot parameters considered in this application.Industrial robot parameters

Latin American Journal of Solids and Structures 13 (2016) 1802-1822Table 7 :
Meta-model coefficients and determination coefficient ( 2 R ) for the industrial robot application.

Table 8
presents the geometric characteristics of the discs of the rotor.

Table 8 :
Geometric characteristics of the discs of the rotating machine.

Table 9 :
Design space used in the fitting of the finite element model.

Table 10 :
Optimal solutions determined from the deterministic and robust process ( i Cr is the critical speed).