A New Approach for Robust Design Optimization Based on the Concepts of Fuzzy Logic and Preference Function

1.Malek-Ashtar University of Technology – Department of Mechanical and Aerospace Engineering – Shahinshahr/Isfahan – Iran Correspondence author: Ali Babaei | Malek-Ashtar University of Technology – Department of Mechanical and Aerospace engineering | Ferdowsi Street – Malek-Ashtar University of Technology | CEP: 83.145-115 – Shahinshahr/Isfahan – Iran | E-mail: arbabaei@aut.ac.ir Received: Jun. 4, 2017 | Accepted: Sept. 19, 2017 Section Editor: Paulo Greco ABSTRACT: In this study, an efficient methodology is proposed for robust design optimization by using preference function and fuzzy logic concepts. In this method, the experience of experts is used as an important source of information during the design optimization process. The case study in this research is wing design optimization of Boeing 747. Optimization problem has two objective functions (wing weight and wing drag) so that they are transformed into new forms of objective functions based on fuzzy preference functions. Design constraints include transformation of fuel tank volume and lift coefficient into new constraints based on fuzzy preference function. The considered uncertainties are cruise velocity and altitude, which Monte Carlo simulation method is used for modeling them. The non-dominated sorting genetic algorithm is used as the optimization algorithm that can generate set of solutions as Pareto frontier. Ultimate distance concept is used for selecting the best solution among Pareto frontier. The results of the probabilistic analysis show that the obtained configuration is less sensitive to uncertainties.


INTRODUCTION
Engineering design cycle is a process of the subsystems formulation to meet human needs.To increase design performance, it is necessary to use the design optimization methodologies in the systems design (Hao et al. 2015).Complex systems usually have multiple conflicting disciplines, which difficult its compromise.Moreover, sequential or single objective optimization will result in low-efficiency solutions or even non-optimal solutions (Huang et al. 2008).Multi-objective optimization is a method that optimizes a group of objective functions simultaneously (Huang et al. 2006).
In today's world, systems should also be robust in the face of uncertainties.In other words, the performance of the optimal solutions must have fewer changes because of the variations of design variables and operational conditions.So a measure of robustness should be introduced during optimization and design (Gaspar-Cunha and Covas 2008).Taguchi introduces robust design at first to enhance product quality so that the quality is not sensitive to variations (Wan et al. 2011).The purpose of this method is to minimize the deviation of system response due to uncertainties.In classical design methods for considering uncertainty in a system design, the designers apply drastic tolerances and large safety factors (Jun et al. 2011).Assigning the values of these parameters is based on past experience of designer and has these drawbacks: a) specifying the values of safety factor for new systems and materials without any experience is difficult; b) in the design process, it is not easy to measure robustness and xx/xx 02/22 reliability, so the designer cannot compare the system from standpoint of resistant and optimality; c) using this method in the optimization problem limits design space and therefore the new design is very conservative.So with attention to these drawbacks, it is necessary to use robust design optimization for considering uncertainty in system design (Roshanian et al. 2012).
Fuzzy logic is a methodology based on the experience of human and has been developed to deal with vague and uncertain systems.Fuzzy theory is a systematic process to convert the experience of human into nonlinear mapping and it is an important aspect of this methodology.This technique is used as a modeling method for complex systems.The main core of a fuzzy system is a set of "If-then" rules obtained from the experience of experts.Some of the advantages that make this theory as an efficient technique in engineering applications are proper simplicity and speed, no need for any complex calculations, finding acceptable answers in a short period of time, and using the experience of experts (Wang 1997).Jaeger et al. (2013) proposed a new method for considering uncertainties in the aircraft conceptual design phase.This method is based on the uncertainties distribution so that their standard deviations are variable.Uncertainties have been considered in design and model variables as probabilistic setting and the uncertainties update from the historical database at each step of the optimization process.Shah et al. (2015) presented a robust optimization algorithm for airfoil design under mixed uncertainties (inherent and epistemic) using by the multi-fidelity approach.This method uses stochastic expansions to generate surrogate models.To reduce the computational cost, the high-fidelity CFD model is used.Probabilistic modeling and interval analysis methods are combined to uncertainty modeling.Dashilewicz et al. (2011) studied the disciplinary uncertainty effects in multi-objective optimization for aircraft conceptual design.By analyzing the Pareto frontier, the decision makers can judge between the expected performance and resistance so they can determine regions of design space that are proper or dangerous.Zhang et al. (2012) studied an effective computational approach for aerodynamic robust optimization under aleatory and epistemic uncertainties using by stochastic expansions that are based on non-intrusive polynomial method.The stochastic surfaces are used as surrogates in the optimization process.The stochastic surface is a function of the design variables and uncertain variables.In this paper, two stochastic optimization formulations have been investigated: a) optimization under pure aleatory uncertainty; and b) optimization under mixed uncertainty.Mach number is the aleatory uncertain variable, and a k factor multiplied by the turbulent eddy viscosity coefficient is the epistemic uncertain variable.Messac and Ismail-Yahaya (2002) developed a physical programming-based robust design optimization (PP-RDO) method.This technique is based on physical programming and robust design optimization methods that allow the designer to express robustness wishes in physical meaningful terms.Du and Wang (2016) discussed a method with dynamic surrogates models based on fuzzy clustering analysis to improve the efficiency of uncertainty analysis and reduce the effect of error propagation on uncertainty models.Fuzzy clustering analysis applied to the sample points before the generation of surrogate models.The results show that the amount of uncertainty analysis can be decreased effectively, while the optimized performance can satisfy the reliability and robustness.Bashiri and Hosseininezhad (2009) proposed a methodology for optimizing multi-response surface in robust design that optimizes mean and variance by applying fuzzy set theory.A fuzzy programming method is expressed to solve the problem by applying the degree of satisfaction from each objective.Park et al. (2011) proposed the reliability-based design optimization (RBDO) and possibility-based design optimization (PBDO) methods to handle uncertainties in multidisciplinary design optimization when low fidelity analysis tools are used.The RBDO method uses probability density function and PBDO method uses fuzzy input for uncertain parameters.Zhao and Xue (2010) studied a new parametric design approach with non-deterministic neural network and fuzzy relationships for considering uncertainties.In this paper, probabilistic uncertainties are presented by neural network relationships and possibility uncertainties are presented by fuzzy relationships.This new parametric design method has been used in the design of a solid oxide fuel cell system.
In this study, a robust design optimization methodology based on the concepts of fuzzy logic and preference function is proposed.Preference function is a function that shows the satisfactory of system response.In this paper, the preference functions of objective functions and constraints are formed using fuzzy logic.Using the experience and knowledge of experts (designers) xx/xx 04/22

UNCERTAINTY SIMULATION AND MODELING OF OBJECTIVE FUNCTIONS AND CONSTRAINTS
There are many approaches to modeling uncertainty, from which one of the most famous is MCS.This method was developed in 1940 by Stan Ulam and John Von Neumann and is a computational algorithm that performs repetitive sampling and simulation.Unlike the many probability methods, this technique needs to a little information about the statistic and probability.If the sufficient number of sample points to be used, the results of MCS are completely accurate.So this method is a standard reference for evaluating other methods (Yao et al. 2011).The flowchart of MCS is shown in Fig. 2. To do simulations for different sample points, the modeling of objective functions and constraints should be performed at first.After simulation, the mean and standard deviation of objective functions and constraints are produced.These values are inputs to the fuzzy system to form preference functions (see Fig. 1).

FUZZY PREFERENCE FUNCTION Brief Description of Fuzzy Logic
Fuzzy systems are knowledge-based systems.The core of a fuzzy system is a knowledge base that has been consisted of fuzzy If-Then rules.A fuzzy If-Then rule is an If-Then statement in which some words are characterized by continuous membership functions.For example, the following statement is a fuzzy If-Then rule: If the speed of a car is high, Then apply less force to the accelerator where the words "high" and "less" are characterized by the membership functions (Wang 1997).The general configuration of a fuzzy system is shown in Fig. 3.It is composed of four main parts: fuzzifier, fuzzy rule base set, fuzzy inference engine, and defuzzifier.Fuzzy sets consist of a number of so-called membership functions to describe the fuzziness of variables.The fuzzifier component converts crisp input values into fuzzy sets.Fuzzy sets enter the inference engine as inputs.The decisions are made upon the fuzzy If-Then linguistic rule base set.The inference engine as a mathematical mechanism is applied to manipulate fuzzy sets and aggregate the rules to make the decisions.The outputs of the inference engine are fuzzy, but the system output must be crisp.Therefore, the defuzzifier part maps the fuzzy sets to crisp outputs (Oroumieh et al. 2013).

Preference Function Creation
In this block, the mean value (μ) and standard deviation (σ) of objective functions and constraints are transformed into the new objective functions and new constraints as preference functions.The reason for using this function is that experience of the designer (or experts) can be used during design optimization (by expressing his or her preferences relative to objective functions and constraints) and increases the ability to achieve more practical solutions.To create the preference function, the designer defines satisfaction degrees for any objective function (and constraint) in the horizontal axis.In other words, the designer classifies horizontal axis to different regions in terms of satisfaction of objective function (and constraint) (see Fig. 4).Then the vertical axis, which indicates the preference function, is divided into several regions same as the horizontal axis.It is worth noting that the preference function value is between zero and one, and maximization of preference function (the greatest satisfaction degree) is the purpose of optimization.In this study, the relationship between horizontal and vertical axis is created using fuzzy logic.On the other hand, the preference functions are generated using fuzzy logic.Fuzzy logic is used because we can also benefit from the experience of experts in the creation of fuzzy rules (in addition to determining the satisfaction degrees for any objective function and constraint).In other words, we can use designer experiences in design optimization using the fuzzy preference function definition by two methods:

•
Expressing the satisfaction degrees for objective functions and constraints.

•
Creating of fuzzy rules.

OPTIMIZATION ALGORITHM AND THE FINAL SELECTION METHOD
Genetic algorithms are stochastic optimization algorithms suitable for optimization of complex problems with unknown search space.These algorithms are programming techniques that use genetic evolution as a pattern of problems solution.This is a global optimization algorithm based on the principle of survival of fittest, natural selection mechanism, and reproduction.Because of the genetic algorithm's ability in solving complex problems, it is selected as the optimization algorithm.Since we are dealing with more than one objective function, an improved version of GA called non-dominated sorting genetic algorithm is used in this study.This algorithm finds set of optimal solutions (Pareto frontiers) by adding an essential operator to general single objective GA. Figure 5 shows the flowchart of this algorithm (Babaei et al. 2015).
Since the solutions in the Pareto frontier don't have the preference to each other completely, the concept of utopian distance is used in this study to select the final solution.This distance defines as follows (Eq.1):

Preference function
where: d ut is utopian distance; i is the number of objective function; and f ut is the utopian value obtained from a single objective Genetic algorithms are stochastic optimization algorithms suitable for optimization of complex problems with unknown search space.These algorithms are programming techniques that use genetic evolution as a pattern of problems solution.This is a global optimization algorithm based on the principle of survival of fittest, natural selection mechanism, and reproduction.Because of the genetic algorithm's ability in solving complex problems, it is selected as the optimization algorithm.Since we are dealing with more than one objective function, an improved version of GA called non-dominated sorting genetic algorithm is used in this study.This algorithm finds set of optimal solutions (Pareto frontiers) by adding an essential operator to general single objective GA. Figure 5 shows the flowchart of this algorithm (Babaei et al. 2015).
Since the solutions in the Pareto frontier don't have the preference to each other completely, the concept of utopian distance is used in this study to select the final solution.This distance defines as follows (Eq.1): where: dut is utopian distance; i is the number of objective function; and fut is the utopian value obtained from a single objective optimization for each objective function.Among the Pareto (1) xx/xx 06/22 optimization for each objective function.Among the Pareto frontier, the solution that has the lower value of d ut is the best solution because the standpoint of objective functions is closest to the utopian values.

Definition of Objective Functions and Design Variables
In this study, two design optimizations have been done for wing optimization, which includes: deterministic optimization (DO), and robust optimization (RO).In the first case, the problem can be formulated as follows (Eq.2): where: W W is wing weihgt, Dr cr is wing cruise drag, V is fuel tank volume, and C L is lift coefficient.
In the second case, the problem can be formulated as follows (Eq.3): where: F Ww , F Dr cr , and F C L are the preference functions of wing weight, wing cruise drag, fuel tank volume, and lift coefficient, respectively.
Because the main purpose of this paper is to present a new method for robust design and for simpler presentation of the application of the proposed method, only three variables are considered.Wing root chord (C r ), wingspan (b), and sweepback angle (Λ 0.5 ) are considered as design variables for both optimizations.Numerical ranges of these variables have been shown in Table 1.

UNCERTAINTY SIMULATION
Aircraft design usually is accomplished based on a certain mission.Often every aerial vehicle is designed for a specific mission with certain cruise altitude and speed.It is very economical that one aircraft can be optimum or near-optimum for other missions.In other words, it is suitable that one aircraft is designed for a range of cruise altitude or speed.This topic could also be the case for transport jet.Therefore, cruise altitude and flight speed are considered as uncertainties in this research.A normal distribution is used to generate sampling points (X = Xm + (∆x/3) randn (1,N)).In this study, the values of h m , V m , Δh m and ΔV are considered 25000 ft, 762 ft/s, 10000 ft and 182 ft/s, respectively.The values of h m and V m are selected based on cruise condition of Boeing 747, and the values of Δh and ΔV are selected based on designer experiences.

MODELING OF OBJECTIVE FUNCTIONS AND CONSTRAINT
As already noted, wing weight and drag force of Boeing 747 are considered as objective functions while fuel tank volume and lift coefficient are as constraints.The considered condition for wing design optimization is cruise phase.It is necessary to extract equations of wing weight, drag force, fuel tank volume and lift coefficient.Like all modeling, modeling of this research has been done with consideration of assumptions.It is worth noting that some parameters are considered constant (such as taper ratio, Mach number, dynamic pressure, wing incidence angle, cruise angle of attack, fuselage diameter, etc.) are related to Boeing 747.In robust optimization, the parameters like Mach number and dynamic pressure are considered variable in modeling process (unlike the deterministic optimization).

Modeling of Wing Weight
The wing weight of an aircraft can be calculated from Eq. 4 (Sadraey 2012): xx/xx 08/22 process (unlike the deterministic optimization).

Modeling of Wing Weight
The wing weight of an aircraft can be calculated from Eq. 4 (Sadraey 2012): ) q.r  q.qt  (4) where: SW is wing area, MAC is mean aerodynamic chord, ( % A ) Zde is the maximum of thickness to chord ratio, ρmat is density of construction material, Kρ is density factor, AR is aspect ratio, nult is ultimate load factor, λ is taper ratio and g is gravity constant.
If root chord, taper ratio, and span are known, then wing area is calculated as follows (Eq.

5)
: Taper ratio for Boeing 747 is 0.28, therefore the wing area is obtained with placement of this value in Eq. 5, obtaining (Eq.6): The mean chord of the wing is calculated as follows (Eq.7): where: S W is wing area, MAC is mean aerodynamic chord, (t/c)max is the maximum of thickness to chord ratio, ρ mat is density of construction material, K ρ is density factor, AR is aspect ratio, n ult is ultimate load factor, λ is taper ratio and g is gravity constant.
If root chord, taper ratio, and span are known, then wing area is calculated as follows (Eq.5): Taper ratio for Boeing 747 is 0.28, therefore the wing area is obtained with placement of this value in Eq. 5, obtaining (Eq.6): The mean chord of the wing is calculated as follows (Eq.7): It is assumed that the airfoil of the wing is NACA2412, so (Eq.8): The density of construction material and wing density factor for the aerospace aluminum alloy are (Eq.9): It is assumed that the maximum load factor of Boeing 747 is 3 and therefore the value of ultimate load factor is considered (Eq.10) (Roskam 1997): Note that 1.5 is a safety factor and it is considered based on structural requirements.The equation of aspect ratio according to the Eq. 6 is expressed as follows (Eq.11): Finally, the main equation of wing weight is obtained with placement above equations in Eq. 4, obtaining (Eq.12): Finally, the main equation of wing weight is obtained with placement above e Eq. 4, obtaining (Eq.12): (1

Modeling of Wing Drag
The main equation of the drag is expressed as follows (Eq.13):

Modeling of Wing Drag
The main equation of the drag is expressed as follows (Eq.13): xx/xx 09/22

Modeling of Wing Drag
The main equation of the drag is expressed as follows (Eq.13): where:  ô is dynamic pressure, and CD is drag coefficient that it is obtained from In Eq. 14, CD0 is zero-lift coefficient that can be calculated as (Eq.15) (R where: f is equivalent parasite area.This parameter is obtained as follows (Eq.16 where: a and b are constant, and Swet is the wetted area. The main equation of the drag is expressed as follows (Eq.13): where:  ô is dynamic pressure, and CD is drag coefficient that it is obtained from Eq  R =  Rq +  L 0 In Eq. 14, CD0 is zero-lift coefficient that can be calculated as (Eq.15) (Ros where: f is equivalent parasite area.This parameter is obtained as follows (Eq.16): where: a and b are constant, and Swet is the wetted area.

𝐷𝐷 = 𝑞𝑞 ô𝑆𝑆𝐶𝐶 R
where:  ô is dynamic pressure, and CD is drag coefficient that it is obtained from In Eq. 14, CD0 is zero-lift coefficient that can be calculated as (Eq.15) ( where: f is equivalent parasite area.This parameter is obtained as follows (Eq.
where: a and b are constant, and Swet is the wetted area.
Since the purpose of this research is wing optimization, therefore it is evident th area changes in addition to wing area in the optimization process and finally equivalen area and zero-lift drag coefficient will change.So, to calculate zero-lift drag coeffi equation for the wetted area must be determined at first and then according to Eqs. 15 zero-lift drag coefficient can be calculated.The wetted area is shown in Fig. 6.To find the equation for the wetted area in terms of design variables, it is nec find an equation to express chord distribution along the wingspan.This equation is using a simple geometric relationship as follows (Eq.17): where: y is chord distribution, and x is along the wing span.To find the wetted area, the is the only unknown parameter that is obtained from Eq. 17. Fuselage diame meter for Boeing 747.Finally, the main equation is obtained for the wetted area (Eq.18) Since the purpose of this research is wing optimization, therefore it is evident that wetted area changes in addition to wing area in the optimization process and finally equivalent parasit area and zero-lift drag coefficient will change.So, to calculate zero-lift drag coefficient, an equation for the wetted area must be determined at first and then according to Eqs. 15 and 16 zero-lift drag coefficient can be calculated.The wetted area is shown in Fig. 6.To find the equation for the wetted area in terms of design variables, it is necessary to find an equation to express chord distribution along the wingspan.This equation is obtaine using a simple geometric relationship as follows (Eq.17): where: y is chord distribution, and x is along the wing span.To find the wetted area, the chord i  = ü † 0 is the only unknown parameter that is obtained from Eq. 17. Fuselage diameter is 6.5 meter for Boeing 747.Finally, the main equation is obtained for the wetted area (Eq.18): where: q is dynamic pressure, and C D is drag coefficient that it is obtained from Eq. 14: In Eq. 14, C D0 is zero-lift coefficient that can be calculated as (Eq.15) (Roskam 1997): where: f is equivalent parasite area.This parameter is obtained as follows (Eq.16): where: a and b are constant, and S wet is the wetted area.
Since the purpose of this research is wing optimization, therefore it is evident that wetted area changes in addition to wing area in the optimization process and finally equivalent parasite area and zero-lift drag coefficient will change.So, to calculate zero-lift drag coefficient, an equation for the wetted area must be determined at first and then according to Eqs. 15 and 16, zero-lift drag coefficient can be calculated.The wetted area is shown in Fig. 6.To find the equation for the wetted area in terms of design variables, it is necessary to find an equation to express chord distribution along the wingspan.This equation is obtained using a simple geometric relationship as follows (Eq.17): where: y is chord distribution, and x is along the wing span.To find the wetted area, the chord in x = d f /2 is the only unknown parameter that is obtained from Eq. 17. Fuselage diameter is 6.5 meter for Boeing 747.Finally, the main equation is obtained for the wetted area (Eq.18): Swet w where: C t is wing tip chord, and d f is fuselage diameter.
The equivalent parasite area (Eq.19) and zero-lift drag coefficient (Eq.20) can be calculated using Eqs.6, 15, 16 and 18. q.rtV T v (20) To calculate the drag coefficient according to Eq. 14, an equation for lift coefficient must be achieved.The lift coefficient is assumed as CL0 + CLα, where CL0 is lift coefficient at zero angles of attack, CLα is lift curve slope and α is the angle of attack.The lift curve slope is calculated as follows (Eq.21) (Roskam 1997).
In Eq.21, k and β are obtained from the Eqs.22 and 23.
To calculate the drag coefficient according to Eq. 14, an equation for lift coefficie be achieved.The lift coefficient is assumed as CL0 + CLα, where CL0 is lift coefficient angles of attack, CLα is lift curve slope and α is the angle of attack.The lift curve calculated as follows (Eq.21) (Roskam 1997).
In Eq.21, k and β are obtained from the Eqs.22 and 23.
To calculate the drag coefficient according to Eq. 14, an equation for lift c be achieved.The lift coefficient is assumed as CL0 + CLα, where CL0 is lift coe angles of attack, CLα is lift curve slope and α is the angle of attack.The lift calculated as follows (Eq.21) (Roskam 1997).
In Eq.21, k and β are obtained from the Eqs.22 and 23.
To calculate the drag coefficient according to Eq. 14, an equation for lift be achieved.The lift coefficient is assumed as CL0 + CLα, where CL0 is lift co angles of attack, CLα is lift curve slope and α is the angle of attack.The lif calculated as follows (Eq.21) (Roskam 1997).
In Eq.21, k and β are obtained from the Eqs.22 and 23. (24) To calculate the drag coefficient according to Eq. 14, an equation for lift coefficient must be achieved.The lift coefficient is assumed as C L0 + C Lα , where C L0 is lift coefficient at zero angles of attack, C Lα is lift curve slope and α is the angle of attack.The lift curve slope is calculated as follows (Eq.21) (Roskam 1997).
In Eq.21, k and β are obtained from the Eqs.22 and 23.
where: M is Mach number and C lα is airfoil lift curve slope.
The Mach number of Boeing 747 in cruise condition is 0.85, therefore, C lα | atM is calculated as follows (Eq.24).The value of C lα | atM = 0 is assumed to be equal to the value of the NACA 2412.
where: M is Mach number and Clα is airfoil lift curve slope.
The Mach number of Boeing 747 in cruise condition is 0.85, therefore,  à ± ª d% Ω calculated as follows (Eq.24).The value of  à ± ª d% Ω3q is assumed to be equal to the value of NACA 2412.
With placement of Eqs.22, 23, 24, and aspect ratio in Eq.21, the final equation (Eq. 2 is obtained for lift curve slope. The lift coefficient for zero angles of attack is not constant, similar to zero-lift dr coefficient and changes during the optimization process.This parameter is calculated as follo (Eq.26): where: M is Mach number and Clα is airfoil lift curve slope.
The Mach number of Boeing 747 in cruise condition is 0.85, therefore,  à ± ª d% Ω is calculated as follows (Eq.24).The value of  à ± ª d% Ω3q is assumed to be equal to the value of the NACA 2412.
With placement of Eqs.22, 23, 24, and aspect ratio in Eq.21, the final equation (Eq.25) is obtained for lift curve slope.
The lift coefficient for zero angles of attack is not constant, similar to zero-lift drag coefficient and changes during the optimization process.This parameter is calculated as follows (Eq.26): With placement of Eqs.22, 23, 24, and aspect ratio in Eq.21, the final equation (Eq.25) is obtained for lift curve slope.
The lift coefficient for zero angles of attack is not constant, similar to zero-lift drag coefficient and changes during the optimization process.This parameter is calculated as follows (Eq.26): where: C L0 wf is lift coefficient at zero angles of attack for combination of wing body, is horizontal tail lift curve slope, η h is horizontal tail efficiency, S h is horizontal tail area, i h is horizontal tail incidence angle, ε 0h is downwash angle at the horizontal tail, is canard lift curve slope, η c is canard efficiency, S c is canard area, i c is canard incidence angle, and ε 0c is up-wash angle at the canard.
The last term in the Eq. 26 is related to the canard that has been neglected in this study.The second term is related to horizontal tail, which is negligible compared to the first term.Therefore lift coefficient for zero angles of attack is calculated from Eq. 27: The last term in the Eq. 26 is related to the canard that has been neglected in thi The second term is related to horizontal tail, which is negligible compared to the fir Therefore lift coefficient for zero angles of attack is calculated from Eq. 27: where: iw is wing incidence angle, α0LW is wing angle of attack for zero lift and  L ± Q † body lift curve slope.
The incidence angle is 2° for Boeing 747 and α0LW is calculated using Eq.28: where: α0l is airfoil zero lift angle, is the variation of wing zero lift angle of attack twist angle, and εt is wing twist angle.
With the placement of relevant parameters in Eq. 28, the value of this param obtained -1.92.Wing body lift curve slope is obtained as follows (Eq.29): where: i w is wing incidence angle, α 0LW is wing angle of attack for zero lift and is wing body lift curve slope.
The incidence angle is 2° for Boeing 747 and α 0LW is calculated using Eq.28: where: α 0l is airfoil zero lift angle, Δα 0 /ε t is the variation of wing zero lift angle of attack to wing twist angle, and ε t is wing twist angle.
With the placement of relevant parameters in Eq. 28, the value of this parameter is obtained -1.92.Wing body lift curve slope is obtained as follows (Eq.29): where: k wf is wing-body interference factor, which is calculated from Eq. 30: The final equation for lift curve slope is (Eq.31): Finally, according to Eqs. 27 and 29, lift coefficient for zero angles of attack is calculated as follows (Eq.32): Now using Eqs.13 and 14, final equation is obtained for drag force (Eq.33): Dynamic pressure and angle of attack for the cruise condition of Boeing 747 are 8706.72kg/ms 2 and 2.5°, respectively.The main drag equation (Eq.34) is obtained by substituting Eqs. 20, 25 and 30 in Eq. 33:

Modeling of the Fuel Tank Volume
One of the wing tasks is to provide an enclosure for the fuel tank in addition to producing the lift force.Since the design variables change in the optimization process the fuel tank volume will also change, and because the aircraft mission depends on the volume of fuel so this volume should not change relative to the initial volume.
The cross-section fuel tank and its parameters are shown in Fig. 7.It is assumed that the fuel tank is between wing spars in this study.In the actual condition, the wing is the taper and the maximum of airfoil thickness changes from root to tip, so fuel tank is modeled as taper with variable thickness.Fig. 8 shows the considered fuel tank.In the actual condition, the wing is the taper and the maximum of airfoil thickness changes from root to tip, so fuel tank is modeled as taper with variable thickness.Fig. 8 shows the considered fuel tank.In Fig. 8, 2w is the length of the fuel tank, 2d is the width of the fuel tank and h is the fuel tank shell thickness.
The fuel tank volume is obtained from the following equation (Eq.35): It is clear from In Fig. 8, 2w is the length of the fuel tank, 2d is the width of the fuel tank and h is the fuel tank shell thickness.The fuel tank volume is obtained from the following equation (Eq.35): It is clear from Fig. 6 that, to calculate the fuel tank volume, the values of 2d, 2w and h must be determined in the root, mean and tip chords at first.Another important point is the position of spars along the wing because 2w and 2d are dependent on this position.The position of the front and rear spars are assumed as follows (Eq.36) (Roskam 1997): (37) From Eq. 36, the distance of two spars is 0.55 C, so the value of 2w is 0.55 C in different parts of the wing.The values of 2d and h for Boeing 747 are assumed (Eq.37) (Setayandeh 2011): Finally, the values of 2w and 2d are listed in Table 2.The fuel tank volume is 148.8 m 3 for baseline design using Eq.41.This value is considered as the first constraint in the optimization design (ΔV = -148.8= 0).
Using the described equations for the lift coefficient, the lift coefficient was calculated 0.28 for baseline design, so this value is considered as the second constraint.This constraint makes the obtained lift coefficient from optimization algorithm (as an important parameter in cruise phase) not to change, compared to the baseline design (ΔC L = C L -0.28 = 0).

MAKING PREFERENCE FUNCTION USING FUZZY LOGIC
It has been mentioned that an important advantage of this method is using the experience of experts in the design optimization in defining preference function.It was also stated that preference functions are created using fuzzy logic.For this purpose, fuzzy rules are defined by the designers (experts) and based on their experiences.In addition, simplicity and reducing the computational time due to the use of fuzzy logic are other advantages of this function.To use these advantages, preference functions are formed for all objective functions and constraints.

xx/xx 14/22
The fuzzy preference function is achieved using product inference engine, singleton fuzzifier, and center average defuzzifier, as follows (Eq.42): (42) computational time due to the use of fuzzy logic are other advantages of this functio these advantages, preference functions are formed for all objective functions and constr The fuzzy preference function is achieved using product inference engine, fuzzifier, and center average defuzzifier, as follows (Eq.42): ä k¿ë (42 In Eq. 42, y is the center of the fuzzy set, l is the number of fuzzy rules and i is the number of inputs.To form preference functions of drag and lift coefficient, because both of them change with uncertainties, two inputs (μ, σ) must be considered.Mean values and standard deviation are divided into four and three sections, respectively (the membership functions of them are shown in Fig. 9), so 12 fuzzy rules must be created to form preference function.To form preference functions of weight and fuel tank volume only one input (μ) must be considered (these parameters don't change with uncertainties).The mean values of them are divided into four sections so four fuzzy rules must be created to form preference function.The preference function (as the output of fuzzy logic) is divided into four sections.The membership function of preference function is shown in Fig. 9.

Preference Function of Drag
Since the uncertainties affect the drag, standard deviation and the mean value of drag are as inputs of the fuzzy system for making preference function of drag.This process is shown in Fig. 10.   6.If μ D is high and σ D is low then drag preference function is medium (0.7).
7. If μ D is medium and σ D is high then drag preference function is Low (0.35).
8. If μ D is medium and σ D is medium then drag preference function is medium (0.7). 9.If μ D is medium and σ D is low then drag preference function is high (1).
10.If μ D is low and σ D is high then drag preference function is low (0.35).
11.If μ D is low and σ D is medium then drag preference function is high (1).
12. If μ D is low and σ D is low then drag preference function is high (1).

Preference Function of Wing Weight
Unlike the drag, uncertainties don't affect the wing weight, so there is no need to calculate mean and standard deviation in this case.The values of normalized wing weight are as input to the fuzzy system (Fig. 13).
The membership function of normalized wing weight is shown in Fig. 14.The weight preference function is similar to Eq. 42 with a difference in the number of fuzzy rules.

xx/xx 16/22
The fuzzy rule set for wing weight is: 1.If μ W is very high then wing weight preference function is very Low (0). 2. If μ W is high then wing weight preference function is low (0.35).
3. If μ W is medium then wing weight preference function is medium (0.7). 4. If μ W is low then wing weight preference function is high (1).The lift coefficient preference function is similar to Eq. 42 and the fuzzy rule set for lift coefficient is: 1.If μ ΔCl is very high and σ ΔCl is high then lift coefficient preference function is very low (0).

Preference Function of Fuel Tank Volume
Like the wing weight, uncertainties don't affect the fuel tank volume.The process of making preference function for this constraint is similar to Fig. 13 and normalized ΔV is the only input to the fuzzy system.The membership function of normalized fuel tank volume is shown in Fig. 17 The fuzzy rule set for fuel tank volume is: 1.If μ ΔV is very high then fuel tank volume preference function is very low (0).2. If μ ΔV is high then fuel tank volume preference function is low (0.35).
3. If μ ΔV is medium then fuel tank volume preference function is medium (0.7).

OPTIMIZATION ALGORITHM
As mentioned before, the NSGA algorithm is considered as the optimizer.Since the optimization problem is a constrained problem, penalty function method is used for implementation.The specifications of NSGA for optimization in this study are in Table 3.

DESIGN OPTIMIZATION RESULTS
As already mentioned, two design optimizations have been done in this paper that includes: deterministic optimization, and robust optimization.It was also stated that the utopian distance concept is used for final solution selection so the utopian values of wing drag and weight must be calculated at first.From a single objective optimization, utopian values of drag and weight and their specifications were obtained as follows (Tables 4 and 5): Four optimal solutions (Pareto frontier) were achieved for deterministic optimization, presented in Table 6.It is worth noting that the details of deterministic optimization are available in Babaei et al. (2015).As already noted, each point that has lower distance is the best solution among Pareto frontier.So for deterministic optimization, number four has lower distance and is the final solution.The full specifications of this point are given in Table 7. it can be argued that the proposed method is a simple and efficient method for robust design optimization, which can be used from the experience of experts as a valuable source of information during the optimization.Table 11.Variance comparison of deterministic, robust and baseline configuration.

Baseline Configuration OD Configuration RD Configuration
Variance (N) 8.99 × 10 5 11.3 × 10 5 7.476 × 10 5 The importance of design cycle should be noted.In aircraft design, different disciplines and the communications between them are considered.The final design will be obtained from the compromise between these disciplines.The main purpose of this paper is to present an efficient method for robust design optimization in which wing design optimization is considered for implementation.In this paper, the wing performance is optimized alone, and more precise considerations with more analysis must be considered to replace these designs (optimal and robust) with the base design.

CONCLUSION
In this study, a robust optimization methodology was proposed and was applied for the wing of Boeing 747.As already mentioned, the main advantage of this method is the preference function definition using fuzzy logic.The designer can help optimization algorithm to find more practical design using the preparation of fuzzy rules and also determining desirable and undesirable ranges of objective functions and constraints.The preference function definition has advantages such as using of designer experience, good accuracy, and simplicity.Cruise altitude and velocity were considered as uncertainties, and MCS method was used for uncertainty modeling.Deterministic optimization and robust optimization were done in this study and their results

Figure 3 .
Figure 3. Basic configuration of a fuzzy system.

Figure 4 .
Figure 4.An example of preference function for minimization problem

Figure 7 .
Figure 7. Wing cross-section and fuel tank parameters.
Fig.6that, to calculate the fuel tank volume, the values of 2d, 2w and h must be determined in the root, mean and tip chords at first.Another important point is the position of spars along the wing because 2w and 2d are dependent on this position.The position of the front and rear spars are assumed as follows (Eq.36)(Roskam 1997):

Figure 9 .
Figure 9. Membership function of preference functions.

Figure 12 .
Figure 12.Membership function of the standard deviation of normalized drag.

Figure 13 .
Figure 13.Fuzzy preference function of weight force.

Figure 14 .
Figure 14.Membership function of normalized wing weight.

Figure 15 .
Figure 15.Membership function of the mean value of normalized lift coefficient.

Figure 16 .
Figure 16.Membership function of the standard deviation of normalized lift coefficient.

Figure 17 .
Figure 17.Membership function of normalized fuel tank volume. .
4. If μ ΔV is low then fuel tank volume preference function is high (1)

Figure 19 .
Figure 19.Probability distribution function of normalized drag force for robust design (RD), optimum design (OD) and baseline design.

Table 1 .
Limits for design parameters.

Table 2 .
Final values of fuel tank.Therefore, the final equation for the fuel tank volume is calculated as follows(Eqs.38,39, 40, and 41):

ΔC L ) Membership function
2. If μ ΔCl is very high and σ ΔCl is medium then lift coefficient preference function is very low (0).3.If μ ΔCl is very high and σ ΔCl is low then lift coefficient preference function is low (0.35). 4. If μ ΔCl is high and σ ΔCl is high then lift coefficient preference function is very low (0). 5.If μ ΔCl is high and σ ΔClv is medium then lift coefficient preference function is very low (0).6.If μ ΔCl is high and is low then lift coefficient preference function is low (0.35). 7. If μ ΔCl is medium and is high then lift coefficient preference function is low (0.35).If μ ΔCl is medium and σ ΔCl is medium then lift coefficient preference function is medium (0.7). 9.If μ ΔCl is medium and σ ΔCl is low then lift coefficient preference function is medium (0.7). 10.If μ ΔCl is low and σ ΔCl is high then lift coefficient preference function is high (1).11.If μ ΔCl is low and σ ΔCl is medium then lift coefficient preference function is high (1).12.If μ ΔCl is low and σ ΔCl is low then lift coefficient preference function is high (1).

Table 4 .
Utopian values for drag.

Table 5 .
Utopian values for weight.

Table 6 .
Results of deterministic optimization.

Table 7 .
The full specification of final solution.