Adjoint optimization methods for multi-objective aerodynamic design

Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Constraints handling in Nash/Adjoint optimization methods for multi-objective aerodynamic design q Zhili Tang a,⇑, Jacques Périaux b, Jun Dong c a

College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China International Center for Numerical Methods in Engineering (CIMNE), Edificio C1, 08034 Barcelona, Spain c AVIC Aerodynamics Research Institute, Shenyang, China b

a r t i c l e

i n f o

Article history: Received 10 January 2013 Received in revised form 10 December 2013 Accepted 17 December 2013 Available online 31 December 2013 Keywords: Constraint handling Nash equilibrium Multi-objective optimization Aerodynamic design Adjoint methods Pareto front

a b s t r a c t Game theory and its particular Nash games are gaining importance in multi-objective optimization in engineering problems over the past decade. Among ongoing advances in optimization methods and tools, many applications including mathematical modeling with constraints in different areas remain challenges in industrial design environments. This paper describes a constraint handling algorithm to extend the use of Nash games to a more realistic multi-objective aerodynamic optimization when taking into account constraints. A competitive Nash game with constraints is adapted to a cooperative game, in which Nash procedure acts as a sub-game to calculate the equilibrium point and collaborates with a player in charge of constraints. The partial gradient of each objective instead of design variables is used as elitist information exchanged symmetrically in sub-game. Existence and equivalence of the solution are analyzed and proved based on Brouwer fixed point theorem. Numerical experiments on 2D multi-objective aerodynamic optimization with lift and/or geometry constraints are implemented, and results show that constraints are satisfied. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Since John Nash published in 1950 his original contribution to non-cooperative game theory [1], introducing the fundamental concept of Nash equilibrium, the theory has been applied to numerous cases in various areas. Over the past decade, Nash equilibrium theory became an efficient tool to solve multi-objective optimization problems in aerodynamics [2,3] and other relative fields [4,5]. Among ongoing advances in optimization methods and tools, many applications including mathematical modeling with constraints in different areas remain challenges in industrial design environments. However, practical aerodynamic optimization problems involve constraints, such as drag reduction under fixed lift or lift to drag ratio maximization under fixed volume of a wing, etc. The computational design system will attempt to find the optimal solution satisfying the constraints. The Nash equilibrium is the solution of a game based on symmetric information exchanged by Nash players. Equilibrium is reached when each player, constrained by the strategy of others, cannot improve further his own criterion. Constraints associated to objective function are implemented by partial design variables in Classical Nash Game (CNG). Thus symmetric elitist information exchange in the Nash game definitively changes the constraint value enforced by each player’s decision. In

q

This work was completed with the support of National Natural Science Foundation of China under Grant No. NSFC-11272149.

⇑ Corresponding author. Tel.: +86 13913805512.

E-mail address: [email protected] (Z. Tang). 0045-7825/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2013.12.006

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

131

this paper, authors propose a new perspective on constraints handling in Nash/Adjoint optimization methods for multiobjective aerodynamic design. Literature survey Many optimization problems in science and engineering involve a number of constraints which the optimal solution must satisfy. This section presents a brief review of literatures with focussing on how various multi-objective/criterion optimization methods deal with constraints’ satisfaction. K. Deb [6] classified constraint handling methods used in classical optimization algorithms into two groups: (i) generic methods that do not exploit the mathematical structure (whether linear or nonlinear) of the constraint, and (ii) specific methods that are only used to decide if a search point is feasible or not. Generic methods, such as the penalty function method [6–10] and the constraint aggregation method [12–15] as well as hybrid EAs for constrained optimization [16], are popular, because each one of them can be easily applied to any problem without important change in the algorithm. But since these methods are generic, the performance of these methods in most cases is not satisfactory. Penalty functions decrease the fitness of captured infeasible solutions to prefer feasible solutions in the selection process. With Evolutionary Algorithms (EAs), the individual is penalized according to its constraint violation which is the sum of the violation of all constraints. A penalty term can be constructed based on constraint violation. An extended objective function is defined by introducing the penalty term into the original objective function, then to optimize the extended objective function. The penalty function approach involves a number of penalty parameters which must be set right in any problem to obtain feasible solutions. It is difficult to select an appropriate value for the penalty coefficient that adjusts the strength of the penalty. Taking into account the difficulty to determine the penalty parameter, this led researchers to devise sophisticated penalty function approaches such as multi-level penalty functions, dynamic penalty functions [7]. However, ideal control of the coefficient is problem dependent [11] and it is difficult to determine a general control scheme. Farmani and Wright [7] proposed a two-stage adaptive fitness formulation method to ensure that slightly infeasible solutions with a low objective function value remain fitted. The main advantage of this method is that it does not require any parameter tuning. Huang et al. [9] proposed a co-evolutionary differential evolution for constrained optimization. Two populations are used in this method. The first population contains a set of penalty factors and is used to evolve decision solutions, while the second population consists of decision solutions and is employed to adapt penalty factors. These two populations evolve interactively and self-adaptively. K. Deb [6] proposed a pairwise comparison used in tournament selection in EAs which does not need any penalty parameter. In this method, when comparing pairwise individuals, (1) any feasible solution is preferred to any infeasible solution; (2) between two feasible solutions, the one with better objective function value is chosen; and (3) between two infeasible solutions, the one with smaller constraint violation is chosen. Runarsson and Yao [10] proposed a stochastic-rank-based approach. Takahama and Sakai [11] proposed the a constrained method. This method can convert an algorithm for unconstrained optimization problems into an algorithm for constrained optimization problems by the a level comparison [11]. Constraint aggregation transforms the Nonlinear Programming Problem (NLP) into a Multi-objective Optimization Problem (MOP) [13]. A constrained optimization problem is viewed as a constrained satisfaction problem. The main characteristic of this kind of methods are twofold: (1) converting the original constrained optimization problem into unconstrained multiobjective optimization problem and (2) exploiting multi-objective optimization techniques to solve the converted problem. Three mechanisms taken from multi-objective optimization are frequently incorporated into constraint-handling techniques: (1) using Pareto dominance as a selection criterion [14]; (2) using Pareto ranking to assign fitness in such a way that non dominated individuals are assigned a higher fitness value; and (3) splitting the population into sub-populations that are evaluated with respect to the objective function or with respect to a single constraint of the problem. The difficulty is that solving multi-objective optimization problems is a more difficult and expensive task than solving single objective optimization problems. Hybrid methods are a combination of different algorithms and/or mechanisms. Y. Wang et al. [16] proposed a hybrid EAs and an adaptive constraint-handling technique for numerical and engineering constrained optimization problems. The hybrid EAs incorporates simplex crossover and two mutation operators (diversity mutation and improved Breeder Genetic Algorithms (BGAs) mutation) to generate the offspring population. In addition, the adaptive constraint handling technique consists of three main situations, i.e., infeasible situation, semi-feasible situation, and feasible situation. Only one situation is applied at each generation according to whether all the individuals are infeasible, there are feasible and infeasible individuals, or all the individuals are feasible. Meanwhile, one constraint-handling mechanism is designed according to the characteristic of the current situation. Specific methods, such as the cutting plane method, the reduced gradient method, and the gradient projection method, are applicable either to problems having convex feasible regions only or to problems having a few variables, because of increased computational burden in specific methods with the large number of variables. However, generating initial feasible points is difficult and computationally demanding when the feasible region is very small. Especially if the problem has equality constraints, it is almost impossible to find initial feasible points. Special representations and operators method as well as repair algorithms are others specific methods. Special representations and operators are designed to represent only feasible

132

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

solutions and operators are able to preserve the feasibility of the offspring generated. Repair algorithms aim to transform an infeasible solution into a feasible one. The interest in games that allow for constraints in the combined action space is Generalized Nash Equilibrium (GNE) [17], which is almost as old as Nash equilibrium [1]. Ichiishi [18] made a remark that the GNE concept is only useful as a mathematical tool to establish existence theorems in various applied contexts. However, the remark did not stop the mathematical operations-research (OR) researchers from applying and developing this solution concept [19] [20] [21]. The concept has also been used in the politico-economic context of environmental management [22] [23]. The GNE is a very useful and flexible modeling tool and its use is increasing steadily. However, with possibly the exception of existence theory, the study of GNEs in general form is still largely incomplete. The study of GNEPs is still in its infancy, its age notwithstanding [24]. There are still many stimulating open problems to be attacked, both on the theoretical side and on the algorithmic/numerical side. Furthermore, beyond those we listed, there are further interesting topics about which little is known and that are emerging as important. Among these we just mention two. First the problem of selecting one specific solution among the many the GNE usually has (what criteria can we use to establish that a certain solution is preferable to another, how can we compute it?). Second the equilibrium programming with equilibrium constraints (EPEC). EPECs can be viewed as GNEPs whose constraints are in turn defined by some kind of equilibrium condition. These problems naturally arise when considering multi-agents games used in modeling complex competition situations [25–29]. These problems are challenging and extremely hard to analyze. Currently, the literature on games with coupled constraints in contrast to games with fixed strategy sets is much smaller, arguably because of the challenges posed by the constraint coupling [30]. Most of GNEP are used to solve the problems in economics. The application of GNEP in aerodynamic optimization is rare. Furthermore, constraints in aerodynamics are implicit functions of design variables, and they should be evaluated by solving partial differential equations. So the gradient projection to the subspace defined by constraints is not possible to perform numerically in aerodynamic design. Nevertheless, the utilities of each player associated with that player’s own performance and constraints are not sufficient to model a constrained game and to define equilibria; one also needs to model how a player values the fact that other players meet their constraints [31].

2. Observations from survey The following observations can be done from the above survey. Observation 1: Most of constraint handling methods are focused on EAs. Observation 2: Literature concerning the constraints handing in competitive Nash algorithms for multiobjective optimization is not frequently available. No matter how popular a penalty or aggregation method can be, the symmetric elitist information exchange definitely changes the constraint value enforced by each player’s decision once they are used in Nash algorithms to deal with constraints. Here, authors propose a constraint handling method to extend CNG to a more realistic multi-objective aerodynamic optimization taking into account constraints. It can be noticed that most constraint handling methods proposed in literature have succeeded in cooperative games (e.g. Multi-Criteria EAs based on Pareto front concept), this remark inspires authors to convert a competitive Nash game with constraints into a cooperative game, in which the Nash procedure acts as a sub-game to calculate the equilibrium among objectives. The player dealing with Nash equilibrium collaborates with players in charge of constraints. The partial gradient of each objective instead of design variables is used as elitist information exchange in the present algorithm. Numerical experiments on 2D multi-objective aerodynamic optimization with lift and/or geometry constraints are implemented, and results show that constraints are satisfied at convergence. The paper is organized as follows: In the next section, the Nash game for multi-objective aerodynamic optimization is surveyed. Then, Section 4 presents a new constraint handling method in Nash/Adjoint optimization method. Next, in Section 5, we describe numerical experiments on 2D multi-objective constrained aerodynamic shape optimization of an airfoil operating at transonic regimes to validate present algorithms. Finally, concluding remarks are given in Section 6.

3. Nash equilibrium theory for multi-objective aerodynamic optimization Game theory is the formal study of competitive decision making when several players must make choices that potentially affect interests of other players. After John Nash published in 1950 his original work on non-cooperative game theory. The theory has been applied to numerous cases in various areas. Over the past decade, a Nash game has become an efficient technique to solve multi-objective optimization problems in aerodynamics [2,3]. Below the Nash game is surveyed mathematically [32] on game theory. 3.1. Definition of a Nash equilibrium Let N be the number of players. Each player m 2 f1; . . . ; Ng controls the variables xm 2 Rnm . Let x ¼ ðx1 ; . . . ; xN ÞT 2 Rn be the vector formed by all these decision variables, where n :¼ n1 þ    þ nN P N. To emphasize the mth player’s variables within the vector x, we sometimes write x ¼ ðxm ; xm ÞT , where xm subsumes all the other players’ variables.

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

133

Let J m : Rn #R be the mth player’s payoff (or loss) function. We assume that these payoff functions are at least continuous, and we further assume that the functions J m ðxÞ ¼ J m ðxm ; xm Þ are convex in the variable xm . In the classical Nash Equilibrium Problem (NEP), the variable xm belongs to a nonempty, closed and convex set Xm # Rnm ; m ¼ 1; . . . ; N. Let

X :¼ X1      XN

ð1Þ

be the Cartesian product of the strategy sets of each player. Definition 1. A vector xH 2 X is called a Nash equilibrium, or a solution of the NEP, if the block component xH m satisfies H H J m ðxH m ; xm Þ 6 J m ðxm ; xm Þ 8xm 2 Xm

ð2Þ

for all m ¼ 1; . . . ; N. A Nash equilibrium, also called strategic equilibrium, is the result of a game based on symmetric information exchange among players. Equilibrium is reached when each player, constrained by the strategy of others, cannot improve further his own criterion. 3.2. Application to multi-objective aerodynamic optimization In paper [2,3], the adjoint method is combined with Nash equilibrium theory to solve unconstrained multi-objective aerodynamic design problems. The procedure is summarized as follows. Let us consider the following two objective optimization problem

Player1 : min J 1 ðx1 ; x2 Þ Player2 : min J 2 ðx1 ; x2 Þ x1 2X 1

x2 2X 2

ð3Þ

in this game N ¼ 2, vectors x1 and x2 are sub-vectors of design variables x. x1 2 X 1 # Rn1 and x2 2 X 2 # Rn2 ; n1 þ n2 ¼ n P N; x ¼ ðx1 ; x2 Þ 2 X ¼ X 1  X 2 # Rn . Here, J 1 and J 2 can be mathematical functions or functional with can be evaluated for the solution of PDEs (Euler/NS equations in aerodynamics). The Nash Equilibrium (NE) is computed iteratively as follows. Suppose xm1 ¼ ðx1m1 ; xm1 Þ be the best design vector (which component are design variables) at the ðm  1Þst design iter2 m1 ation. At step m, Player1 optimizes sub-vector x1 starting from xm1 to achieve a better choice xm to evaluate 1 1 while using x2 m1 m m1 J 1 . Meanwhile, Player2 optimizes sub-vector x2 starting from x2 to achieve a better choice x2 while using x1 to evaluate m1 J 2 . Each player works with adjoint method to achieve xm ; i ¼ 1; 2. In this case, i from xi m1 xm Þ; 1 ¼ infn J 1 ðx1 ; x2 x1 2R

1

m1 xm ; x2 Þ 2 ¼ infn J 2 ðx1 x2 2R

2

ð4Þ

m Then, players exchange their information symmetrically to form xm ¼ ðxm 1 ; x2 Þ. A Nash equilibrium is reached when neither Player1 nor Player2 can further improve his own criterion by repeating the m above procedure. Design variables xm 1 and x2 in formulation (4) are calculated by using Adjoint methods, see references [2,3,33].

4. Constraints handling in Nash/Adjoint optimization methods for multi-objective constrained aerodynamic design Nash equilibrium theory has been introduced successfully into aerodynamic optimization to solve various unconstrained multi-objective design problems. However, practical design optimization problems involve constraints. Many research activities have been recently focused on evolutionary algorithms to handle constraints in multi-objective aerodynamic optimization. In this section, a new constraints handling method is introduced to extend the use of CNG to more realistic multiobjective aerodynamic optimizations when taking into account constraints. Existing successful constraint handling methods inspire authors to adapt the competitive Nash procedure acting as a sub game for calculating the equilibrium solution among objectives. A player looking for Nash equilibrium collaborates with a player in charge of constraints. The partial gradient of each objective is substituted for design variable and is used as elitist information symmetrically exchanged in present Nash algorithm. 4.1. Algorithms for constraints handling in Nash game In references [2,34], the coupling between an adjoint method and Nash game is described to compute the Nash equilibrium for multi-objective aerodynamic optimization without constraint. Here, a constrained Nash equilibrium algorithm is developed. The following multi-objective constrained aerodynamic optimization problem and game formulation is considered:

( 8 min J 1 ðx1 ; x2 Þ > x1 > > Player1 > > > Subject to g 1 ðx1 ; x2 Þ ¼ g 01 > < ( min J 2 ðx1 ; x2 Þ x2 > > Player2 > > > Subject to g 2 ðx1 ; x2 Þ ¼ g 02 > > : Subject to T ðx1 ; x2 Þ ¼ T 0

ð5Þ

134

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

In the above game formulation, N ¼ 2; J i and g i could be aerodynamic performance (e.g. lift, drag or lift to drag ratio, etc.), and T could be geometry representation (e.g. thickness, contained volume of a wing). Vectors x1 and x2 are sub-vectors of vector x which consists of design variables, where x1 2 X 1 # Rn1 and x2 2 X 2 # Rn2 ; n1 þ n2 ¼ n P N; x ¼ ðx1 ; x2 Þ 2 X 1  X 2 # Rn . It is noticed that constraints enforced by partial design variables in player’s decision can not be satisfied at the end of Nash game iterations, because of the symmetric elitist information exchange in the Nash game [3]. Let us suppose that xm1 ¼ ðx1m1 ; xm1 Þ is the solution at ðm  1Þth dynamic Nash iteration, then the adjoint procedure of each player ensures 2 that m1 g 1 ðxm Þ ¼ g 01 1 ; x2

0 and g 2 ðxm1 ; xm 1 2 Þ ¼ g2

m Once players exchange their information symmetrically to form xm ¼ ðxm 1 ; x2 Þ, then

m 0 g 1 ðxm 1 ; x2 Þ – g 1

m 0 and g 2 ðxm 1 ; x2 Þ – g 2

This is illustrated in Section 5.3 of paper [3] as well as in CNG case1 and CNG case2 of present paper. In this study, we propose two new Nash/Adjoint algorithms and their implementation for multi-objective optimization problems with constraints. Actually, the above Nash problem (5) is equivalent to

8 Player1 : min J 1 ðx1 ; x2 Þ > > x1 > > > > < Player2 : min J 2 ðx1 ; x2 Þ x2

> 2 > X > 2 > 1 > Hi ðg i ðx1 ; x2 Þ  g 0i Þ2 þ 12 cðT ðx1 ; x2 Þ  T 0 Þ > : min h ¼ 2 x1 ;x2

ð6Þ

i¼1

or

8 Player1 : min J 1 ðx1 ; x2 Þ > > x1 > > > > > Player2 : min J 2 ðx1 ; x2 Þ > > x2 > > > < min h1 ¼ 12 H1 ðg 1 ðx1 ; x2 Þ  g 01 Þ2 x1 ;x2 > > > > min h2 ¼ 12 H2 ðg 2 ðx1 ; x2 Þ  g 02 Þ2 > > > x1 ;x2 > > > > 0 2 > : min h3 ¼ 1 cðT ðx1 ; x2 Þ  T Þ x1 ;x2

ð7Þ

2

where Hi and c are arbitrary positive numbers which are not sensitive for satisfying the presented value of constraints. Eq. (6) uses the constraint aggregation technique, all constraints to be satisfied being considered as only one objective. Accordingly, Eq. (7) treats constraints as separate objectives. In the two above games introduced in (6) and (7), constraints are enforced by global design variables. Player 1 and Player 2 of the competitive Nash game can put their resources cooperatively with players in charge of constraints. All players make decision collaboratively following a Pareto game procedure (by ranking, sorting and niching). Doing so, Eqs. (6) and (7) can be reformulated as (8) and (9).

8 8 J 1 ðx1 ; x2 Þ > < Subplayer1 min > x1 > > > Player1 > < : Subplayer2 min J 2 ðx1 ; x2 Þ x2

> 2 > X > 2 > Hi >  ðg i ðx1 ; x2 Þ  g 0i Þ2 þ 2c ðT ðx1 ; x2 Þ  T 0 Þ h ¼ Player2 min > 2 : x1 ;x2

ð8Þ

i¼1

and

8 8 J 1 ðx1 ; x2 Þ > < Subplayer1 min > x1 > > > Player1 > > : Subplayer2 min J 2 ðx1 ; x2 Þ > > x2 > > > < Player2 min h1 ¼ 12 H1 ðg 1 ðx1 ; x2 Þ  g 01 Þ2 x1 ;x2 > > > > h2 ¼ 12 H2 ðg 2 ðx1 ; x2 Þ  g 02 Þ2 Player3 min > > > x1 ;x2 > > > > 0 2 > : Player4 min h3 ¼ 1 cðT ðx1 ; x2 Þ  T Þ x1 ;x2

ð9Þ

2

As formulated in (9), new optimization algorithms (8) or (9) are considered integrally in a cooperative form, in which Nash game acts as a competitive sub-game. In above two new Nash/Adjoint optimization algorithms, the elitist information exchanged at convergence of a Nash sub game is provided by the partial gradient, then constraints are enforced by using

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

135

global gradient of constraint at the end of Nash cycle. Descent gradient is the summation of partial gradient computed by each competitive Nash sub-player and global gradients of constraints. The new geometry is modified along the descent gradient. This new algorithm is implemented as follows: Let xm1 ¼ ðxm1 ; x2m1 Þ are best design variables at step m  1. At step m, Subplayer1 optimizes x1 from xm1 in order to 1 1 m1 m1 achieve a better value xm while using x for evaluating J . Meanwhile, Subplayer2 optimizes x from x in order to 1 2 1 2 2 m1 achieve a better value xm for evaluating J 2 . A Nash cycle for computing the partial gradient [2] of each 2 while using x1 sub-player as follows

8 < Subplayer1 : g subplayer1 ¼ @J 1 ðx1 ;x2 Þ rad @x1

ð10Þ

: Subplayer2 : g subplayer2 ¼ @J 2 ðx1 ;x2 Þ rad @x2

At the end of Nash iteration, Subplayer1 and Subplayer2 send simultaneously their partial gradient to constitute the descent T 1 @J 2 Nash gradient ½@J ;  . @x1 @x2 In order to implement constraints, gradients of constraints with respect to global design variables should been calculated as well, i.e. the following expression (11):

rh ¼

2 X

Hi ðg i ðx1 ; x2 Þ  g 0i Þrg i ðx1 ; x2 Þ þ cðT ðx1 ; x2 Þ  T 0 ÞrT ðx1 ; x2 Þ

ð11Þ

i¼1

Then, the descent gradient of the constrained multi-objective design problem (5) calculated by using Nash game is given in (12)

" @J 1 # GTrad

¼

@x1 @J 2 @x2

þ

2 X

Hi ðg i ðx1 ; x2 Þ 

i¼1

" @gi # g 0i Þ

@x1 @g i @x2

" @T # 0

þ cðT ðx1 ; x2 Þ  T Þ

@x1 @T @x2

ð12Þ

Finally, update aerodynamic shape along the descent gradient, for example the steepest descent method, say:

xm ¼ xm1  kGTrad ðxm1 Þ

ð13Þ

where k is step size. Repeat the above procedure until no one player can improve his objective function only by modifying his own design variables. This compromised solution is the equilibrium of optimization problem (5) with constraints. 4.2. Existence of solution of present Nash algorithms The above iteration procedure (13) can been summarized as

xm ¼ f ðxm1 Þ where f ðxÞ ¼ x  kGTrad ðxÞ

ð14Þ

so, f is a mapping from search space on itself, i.e. design vector x 2 S; f : S#S. If the above iteration process (14) converge, then the solution of present Nash algorithms exists. That is to say there exists a fixed point in search space under the mapping f. In order to prove the existence of the solution, we introduce the famous Brouwer fixed point theorem [35]. Lemma 1. Every continuous function f from a n-dimensional convex compact subset S of a Euclidean space to S itself, then f has a fixed point: for some xH 2 S  Rn ; f ðxH Þ ¼ xH . Among many fixed-point theorems, Brouwer’s is particularly well known, due to its use among numerous fields of mathematics. Hereby, this theorem is used to prove the existence of solution of algorithms developed in Section 4.1. Theorem 1. Suppose N criteria J i ðxÞ and constraints g i ðxÞ as well as T of Nash problem (5) to be smooth functions of the design vector x 2 S # Rn (S: admissible domain; J i and g i 2 C 1 ðSÞ; i ¼ 1; . . . ; N 6 n), then the solution of constrained Nash optimization algorithm (5)–(13) exists.

Proof. In practical aerodynamic optimization, the configuration is parameterized by Bézier splines (or NURBS etc.) for 2-D or Free-Form Deformation (FFD) technique for 3-D geometries. So, geometry constraints (e.g. thickness, contained volume of a h iT @T @T wing) T are easy to be represented as a continuous function, i.e. T 2 C 1 ðSÞ. Hence, the last term cðT ðx1 ; x2 Þ  T 0 Þ @x ; in @x 1 2 (12) is continuous. h iT h iT P @g i @g i 1 @J 2 Since J i and g i 2 C 1 , the first two terms @J and Hi ðg i ðx1 ; x2 Þ  g 0i Þ @x ; in (12) are also continuous. @x1 ; @x2 1 @x2 Therefore, the mapping f in (14) is continuous. According Brouwer fixed point theorem, there exists a fixed point xH 2 S  Rn , such that f ðxH Þ ¼ xH . This fixed point xH is the solution of algorithm (5)–(13). h The above theorem provides a sufficient condition for the existence of a solution, i.e. the iterative algorithm will find a solution once the condition of Theorem 1 is satisfied. Otherwise the existence of a solution is not certain.

136

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

Although the above theorem indicates that there exists a solution for algorithm (5)–(13), it does not answer if the solution converged to the optimization problem (5) or not. That means the solution should be the Nash equilibrium of J 1 jx1 ðx1 ; x2 Þ and J 2 jx2 ðx1 ; x2 Þ in problem (5), meanwhile it should also satisfy all of constraints. The following theorem gives the proof. Theorem 2. if the solution of constrained Nash optimization problem (5) is unique, then the solution of the algorithm is a solution of optimization problem (5). Proof. For the iteration procedure (14), if there exists a unique xH 2 S  Rn , such that f ðxH Þ ¼ xH , this implies GTrad ðxH Þ ¼ 0. In algorithm (5)–(13), Hi and c are arbitrary positive numbers. For any arbitrary positive coefficients Hi and c; GTrad ðxH Þ  0 means that each term in (12) is equal to zero. i.e.

2

 3 @J 1   6 @x1 xH 7 4  5¼0 @J 2 @x2  H

ð15Þ

x

X

2



@g i   6 @x1 H g 0i Þ4 x @g i  @x2  H x

Hi ðg i ðxH Þ 



2

@T   0 6 @x1 xH

3 7 5¼0

ð16Þ

3

cðT ðxH Þ  T Þ4  7 5¼0 @T @x2 

ð17Þ

xH

According to the definition of Nash equilibrium, Eq. (15) means that the solution xH is a Nash equilibrium of optimization problem (5). In general, geometry representations are often linear functions of design variables. So @T @xi ¼ constant for i ¼ 1; . . . ; n. Hence, Eq. (17) implies T ðxH Þ ¼ T 0 , i.e. geometry constraints are satisfied at the solution xH . Concerning aerodynamic constraints, since Hi ; i ¼ 1; . . . ; N are arbitrary positive numbers, so Eq. (16) implies  iT h  @g i  @g  ¼ 0; i ¼ 1; . . . ; N. This can be divided into two cases as follows: ðg i ðxH Þ  g 0i Þ @x  ; i 1 H @x2 H x

x

1. For 8i ¼ 1; . . . ; N, if



@g i  @xj 

xH

¼ 0; j ¼ 1; . . . ; n P N, this means xH is stationary point of function ðg i ðxÞ  g 0i Þ2 . The unique sta2

tionary point of function ðg i ðxÞ  g 0i Þ2 is zero, that is to say ðg i ðxH Þ  g 0i Þ ¼ 0, then g i ðxH Þ ¼ g 0i . I.e. aerodynamic constraints are satisfied at the solution xH .   iT h  @g i  @g i  @g  2. For 8i ¼ 1; . . . ; N, if @x  H – 0; j ¼ 1; . . . ; n P N. Then, ðg i ðxH Þ  g 0i Þ @x ¼ 0 implies g i ðxH Þ ¼ g 0i . I.e. aerodynamic  ; i 1 H @x2 H j

x

x

x

constraints are also satisfied at the solution xH . Therefore, solution xH is the Nash equilibrium of J 1 jx1 ðx1 ; x2 Þ and J 2 jx2 ðx1 ; x2 Þ in problem (5), meanwhile, all of constraints are satisfied. That is to say it is the solution of optimization problem (5). h Hence, algorithms introduced in this paper can converge to the solution of optimization problem (5). 4.3. Distributed numerical procedure of Constrained Nash/Adjoint Algorithms (CNA) Concerning the Nash optimization problem defined in (5), it can be transformed into (9) by using present constraint handling method. Therefore, the distributed computing can be implemented in the following numerical procedure.

Algorithm 1. constraints handling in Nash/Adjoint optimization method n1 old 1. (initialization) input the initial guess of geometry, e.g. xold ¼ ðxold  Rn2  Rn . 1 ; x2 Þ 2 R 2. Distributed computing for each player Load each player’s optimization task on independent processor. For the optimization strategy defined in (9), the following four processes can be run in parallel. (a) Process 1 (Player 1) Do loop for the Nash game cycle: Load each sub-player’s optimization task on an independent sub-processor, then run at its own design point. For sub-player i, adjust xi to minimize J i ; i ¼ 1; 2, simultaneously. i. Do loop for adjoint optimization iterations at design condition i.

137

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

    ii. iii.

Run the CFD solver, outputs being C l ; C d , and flow field variables. Run the adjoint solver. Compute the global gradient. Project the global gradient onto the corresponding subspace, obtain the partial gradient

@J i . @xi

End of the adjoint loop h iT 1 @J 2 Exchange the elitist information symmetrically, and construct the Nash descent gradient @J ; @x1 @x2

End of the Nash cycle End of process 1 (b) Process 2 (Player 2) Do loop for adjoint optimization iterations at design condition 1. i. Run the CFD solver at design condition 1, outputs being C l ; C d , and flow field variables. ii. Run the adjoint solver. h iT 1 @g 1 iii. Compute the global gradient of constraint at design condition 1, H1 ðg 1 ðx1 ; x2 Þ  g 01 Þ @g ; . @x1 @x2 End of the adjoint loop End of process 2 (c) Process 3 (Player 3) Do loop for adjoint optimization iterations at design condition 2. i. Run the CFD solver at design condition 2, outputs being C l ; C d , and flow field variables. ii. Run the adjoint solver. h iT 2 @g 2 iii. Compute the global gradient of constraint at design condition 2, H2 ðg 2 ðx1 ; x2 Þ  g 02 Þ @g ; . @x1 @x2

3. 4. 5. 6.

End of the adjoint loop End of process 3 (d) Process 4 (Player 4) h iT @T @T i. Compute the global gradient of geometry constraints with respect to design variables, cðT ðx1 ; x2 Þ  T 0 Þ @x ; . 1 @x2 End of process 4End of the distributed computing Compute the descent gradient defined in (12). new Update the aerodynamic shape along descent gradient as in (13), e.g. the new shape is xnew ¼ ðxnew 1 ; x2 Þ. Identify if the Nash equilibrium reached or not Reached: Stop. Otherwise: xold ¼ xnew , go to the beginning of the distributed computing.

Technical details of the above procedure, e.g. state and adjoint equations’ discretization, can be found in [34,2]. 5. Application to multi-objective constrained aerodynamic optimization In order to show the validation of present algorithm in handling various constraints (aerodynamic and/or geometry), the following 2D two-point airfoil drag reduction under the lift and/or geometry constraints is performed. The optimization problem is defined as follows.

( 8 min C d1 > > > Player1 at M 1 ¼ 0:74; a ¼ 1:1788930 > > Subject to C l1 ¼ C l1 > > < ( min C d2 0 > > Player2 at M 1 ¼ 0:76; a ¼ 0:825314 > > Subject to C ¼ C l > l 2 2 > > : Subject to Sairfoil ¼ Constant1 or=and T hickness ¼ Constant2

ð18Þ

Where T hickness is airfoil thickness, Sairfoil represents area of an airfoil, constant1 and constant2 are specified thickness and area of airfoil respectively. C d1 and C l1 are drag and lift coefficients of airfoil at M 1 ¼ 0:74 and a ¼ 1:1788930 , accordingly C d2 and C l2 are drag and lift coefficients of airfoil at M 1 ¼ 0:76 and a ¼ 0:8253140 . Here we chose the RAE2822 airfoil as baseline shape. C l1 and C l2 are lift coefficients of the baseline shape at two design points, they are equal to 0:703014 and 0:666488 respectively. In a Nash game, the territory of design variables are split. Désidéri [36] studied the effect of parameterization and split with two players. A split of territory was proposed based on the eigen system of a real symmetric matrix representing the restriction of the primary-criterion Hessian matrix to the subspace tangent to the constraint hyper-surfaces. In our study, the split of design variables is based on the physics of optimization problem since we are dealing with the validation of the algorithm. We use Cc1 and Cc2 to represent the different portion of the airfoil, where Cc1 \ Cc2 ¼ ; and Cc1 [ Cc2 ¼ Cc ; Cc

138

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

1.0

0.8

0.6

0.4

0.2

0.0

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Fig. 1. Fourteen Hicks-Henne shape functions used in airfoil parameterization.

represent airfoil surface. Hick Henne functions (HHF) are used to parameterize an airfoil. The parameterization consists of 28 shape functions to generate the airfoil geometry. HHF are split into two groups: 14 for defining the extrados of the airfoil and 14 for defining the intrados of the airfoil. Design variables are coefficients of these shape functions. Fourteen Hicks-Henne bell-shaped functions are shown on Fig. 1. In order to ensure both players having enough flexibility to modify the airfoil’s shape, we choose Cc1 as basic functions with odd numbers and Cc2 as basic functions with even numbers. With this strategy, either Cc1 or Cc2 can describe individually any smooth airfoil shape. With this split there is no bias on any objective since their importance in problem (18) is equivalent. The flow field is analyzed by 2D Euler equations, which are discretized with the Jameson’s finite-volume approach [37], including a blend of second and fourth-order artificial dissipation terms to suppress oscillations in the shock region. Then the system is solved iteratively by an implicit diagonalized factorization scheme. The adjoint equations are discretized and solved in a way similar to that of the state equations. For a detailed description and the treatment of boundary conditions, we refer to [38]. Using the classical Nash games defined in [4,2,5,34] to solve the above aerodynamic optimization problem (18), we can associate to the problem (18) the following two design optimization problems (19) and (20) without or with geometry constraints (volume, area and thickness) respectively. (1) CNG design case1 (fixed lift constraints)

8 2 > J 1 ¼ C d1 þ H21 ðC l1  C l1 Þ < Player1 min C c1

ð19Þ

2 > : Player2 min J 2 ¼ C d2 þ H22 ðC l2  C l2 Þ

Cc2

(2) CNG design case2 (fixed lift, contained area and thickness constraints)

8 2 Player1 min J 1 ¼ C d1 þ H21 ðC l1  C l1 Þ > > Cc1 > < 2

ð20Þ

Player2 min J 2 ¼ C d2 þ H22 ðC l2  C l2 Þ > Cc2 > > : Subject to Sairfoil ¼ Constant1 and T hickness ¼ Constant2 Table 1 Objectives’ and constraints’ values of design cases (Cc1 is odd, and Cc2 is even).

Objective1

Objective2

Geometry Constraints

Initial Airfoil

CNG Case1

CNG Case2

CNA Case1

CNA Case2

CNA Case3

CNA Case4

Cd Variation C1 Variation

0.0054

0.0018 66.66% 0.6665 5.165%

0.0020 62.96% 0.6702 4.639%

0.0020 62.96% 0.7031 0.043%

0.0020 62.96% 0.7027 0.014%

0.0020 62.96% 0.7011 0.242%

0.0019 64.81% 0.6990 0.541%

Cd Variation C1 Variation

0.0110

0.0040 63.63% 0.6427 3.556%

0.0054 50.91% 0.6399 3.977%

0.0035 68.18% 0.6670 0.090%

0.0040 63.63% 0.6700 0.540%

0.0044 59.99% 0.6671 0.105%

0.0047 57.27% 0.6680 0.240%

Thickness Variation Area Variation

0.1211

0.1166 3.749% 0.0762 2.194%

0.1206 0.421% 0.0779 0.071%

0.1148 5.179% 0.0753 3.306%

0.1181 2.444% 0.0779 0.010%

0.1210 0.096% 0.0793 1.788%

0.1205 0.464% 0.0779 0.079%

0.7028

0.6664

0.0779

139

1.5

1.5

1

1

0.5

0.5

5 × y/c & −C p

5 × y/c & −C p

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

0

Optimized airfoil Optimized pressure initial optimized variation Cl 0.7028 0.6665 −5.171% Initial airfoil Cd 0.0054 0.0018 −67.15% Initial pressure

−0.5

−1

Optimized airfoil Optimized pressure initial optimized variation Cl 0.6664 0.6427 Initial airfoil −3.555% Cd 0.0110 0.0040 Initial pressure −63.28%

−0.5

−1

−1.5

−2

0

−1.5

0

0.1

0.2

0.3

0.4

0.5 x/c

0.6

0.7

0.8

0.9

1

−2

0

0.1

0.2

0.3

0.4

0.5 x/c

0.6

0.7

0.8

0.9

1

Fig. 2. Pressure distribution comparisons between optimized and baseline airfoils at two design points: constraints are treated by using classical Nash game.

Fig. 3. Nash equilibrium procedure by using classical Nash game.

In above two design cases, lift constraints are performed by each player’s decision through partial design variables. So, they can not be satisfied at the end of Nash cycle because the symmetric elitist information exchange modify lift values. It can be seen from Table 1 that lift coefficients of optimized airfoils at two design points are reduced around 4% when compared with initial values. Fig. 2 shows the pressure distribution comparisons between optimized and baseline airfoils at two design points of CNG design case1. The Nash equilibrium on Fig. 3 indicates that lift coefficients at two design points are reduced significantly during the optimization. Among different kinds of constraints used in aerodynamic optimization, the following four design cases are performed with present constrained Nash/Adjoint algorithm (5)–(13).

140

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

(3) CNA design case1 (fixed lift constraints)

8 8 J 1 ¼ C d1 > < Subplayer1 min > Cc1 > > > Player1 > > : Subplayer2 min J 2 ¼ C d2 > < Cc2

ð21Þ

2 > Player2 min h1 ¼ H21 ðC l1  C l1 Þ > > Cc > > > > > : Player3 min h2 ¼ H22 ðC l2  C l2 Þ2

Cc

(4) CNA design case2 (fixed lift and contained area constraints)

8 8 J 1 ¼ C d1 > < Subplayer1 min > Cc1 > > > Player1 > > : Subplayer2 min J 2 ¼ C d2 > > > Cc2 > > < 2 H1 Player2 min h1 ¼ 2 ðC l1  C l1 Þ Cc > > > > > h2 ¼ H22 ðC l2  C l2 Þ2 > > Player3 min Cc > > > > > Player4 min h3 ¼ c1 ðSairfoil  Constant1Þ2 :

ð22Þ

2

Cc

(5) CNA design case3 (fixed lift and thickness constraints)

8 8 J 1 ¼ C d1 > < Subplayer1 min > > Cc1 > > Player1 > > : Subplayer2 min J 2 ¼ C d2 > > > Cc2 > > < 2 H1 Player2 min h1 ¼ 2 ðC l1  C l1 Þ Cc > > > > > h2 ¼ H22 ðC l2  C l2 Þ2 > > Player3 min Cc > > > > > Player4 min h3 ¼ c2 ðT hickness  Constant2Þ2 :

ð23Þ

2

Cc

(6) CNA design case4 (fixed lift, contained area and thickness constraints)

8 8 J 1 ¼ C d1 > < Subplayer1 min > > Cc1 > > Player1 > > : Subplayer2 min J 2 ¼ C d2 > > > Cc2 > > < 2 H1 Player2 min h1 ¼ 2 ðC l1  C l1 Þ Cc > > > > > h2 ¼ H22 ðC l2  C l2 Þ2 > > Player3 min Cc > > > > > Player4 min h3 ¼ c1 ðSairfoil  Constant1Þ2 þ c2 ðT hickness  Constant2Þ2 : 2

2

1.5

1.5

1

1

0.5

0.5

0 Optimized airfoil Optimized pressure initial optimized variation Cl 0.7028 0.7031 Initial airfoil 0.040% C Initial pressure −62.08% d 0.0054 0.0020

−0.5

5 × y/c & −Cp

5 × y/c & −Cp

Cc

0

−1

−1.5

−1.5

0

0.1

0.2

0.3

0.4

0.5 x/c

0.6

0.7

0.8

0.9

1

Optimized airfoil Optimized pressure Initial airfoil Initial pressure

−0.5

−1

−2

ð24Þ

−2

0

0.1

0.2

0.3

0.4

0.5 x/c

initial optimized variation Cl 0.6668 0.6670 0.027% Cd 0.0109 0.0035 −67.57%

0.6

0.7

0.8

0.9

1

Fig. 4. Pressure distribution comparisons between optimized and baseline airfoils at two design points: constraints are treated by using constrained Nash/ Adjoint algorithms.

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

141

In above four design cases, lift constraints are treated as separate objectives, and they are enforced in additional players’ decision by global design variables. The optimization results on Table 1 (the territory of design variables are split into two parts Cc1 and Cc2 , where Cc1 represents basic functions with odd numbers and Cc2 represents basic functions with even numbers) reveal that all of constraints are satisfied. The following Fig. 4 shows the pressure distribution comparisons between optimized and baseline airfoils at two design points of CNA design case1. The Nash equilibrium procedure of this design case is shown on Fig. 5, which indicates that lift coefficients at two design points are preserved at the end of optimization. Above six design cases are performed on a PC with 2.1 GHz CPU. CPU costs are given on the following Table 2, which indicates that CPU cost of present constraint handling algorithm is equivalent to that of CNG algorithm. This is due to the same number of calls to flow solver and adjoint solver in CNG and CNA algorithms. In addition, we exchange territories of subplayer1 with subplayer2: this permutation means that Cc1 is basic functions with even numbers and Cc2 is basic functions with odd numbers. Then re-perform the optimization experiments of above six design cases. Results are shown on Table 3, which are similar with that of on Table 1. This reveals that both players do have enough flexibility to modify the airfoil’s shape, since either Cc1 or Cc2 are able to describe any smooth airfoil shape individually. It also reveals that there is no bias considering in present split. In order to compare performance of Nash equilibria with that of Pareto front, the Pareto front of the constrained multi objective optimization problem (18) has also been calculated by minimizing weighted sum of both objectives with con2 2 strains, i.e. minCc J ¼ x1 C d1 þ x2 C d2 þ H21 ðC l1  C l1 Þ þ H22 ðC l2  C l2 Þ þ c21 ðSairfoil  Constant1Þ2 þ c22 ðT hickness  Constant2Þ2 by P2 varying two weights x1 and x2 , where i¼1 xi ¼ 1 and xi 2 ½0; 1 for i ¼ 1; 2. Nash equilibria are positioned on Fig. 6 (where C d1 and C d2 are drag coefficients of airfoil at two optimization conditions respectively), which indicates that NEs are almost as good as candidate solutions of the Pareto front. The present CNA approach has been used successfully in other complex optimizations problems like 2-D constrained airfoil drag reduction taking account of uncertainties on flight condition parameters [39]. As a general comment on the numerical treatments of constraints in optimization problems, other approaches (such as complex design of gas turbines) have been considered using modular design instead of integrated design. The main idea has been to subdivide the global design problems into a number of simpler sub problems with the pros and cons to delay or not the treatment of constraints during the multi objective optimization procedure. The choice of satisfying constraints depends also about its complexity. In aerodynamics, there are constraints (fixed constant lift) depending of the Mach number through non linear Partial Differential Equations and other (fixed geometry constraint independent to the Mach number). In authors’ present approach it was decided to treat lift constraints dynamically during the Nash optimization introducing the geometry constraints as an objective of an inverse problem. The flexibility of the approach considered in this paper allows the algorithm satisfying constraints only at the end of the multi objective optimization providing to Nash players search spaces including non feasible solutions. This strategy maximizes efficiency while minimizing risk of non convergence of the mathematical formulation of a Nash game. Other strategies for constraints treatment can be found [40–42].

Fig. 5. Nash equilibrium procedure by using constrained Nash/Adjoint algorithms.

142

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

Table 2 CPU costs of CNG and CNA algorithms.

CPU cost (mins)

CNG case1

CNG case2

CNA case1

CNA case2

CNA case3

CNA case4

25.1

25.1

25.2

25.2

25.2

25.2

Table 3 Objectives’ and constraints’ values of design cases (Cc1 is even, and Cc2 is odd).

Objective 1

Objective 2

Geometry Constraints

Initial Airfoil

CNG Case1

CNG Case2

CNA Case1

CNA Case2

CNA Case3

CNA Case4

Cd Variation C1 Variation

0.0054

0.0021 61.11% 0.6733 4.197%

0.0022 59.26% 0.6744 4.041%

0.0023 57.41% 0.7029 0.014%

0.0023 57.41% 0.7027 0.014%

0.0022 59.26% 0.7021 0.099%

0.0022 59.26% 0.7011 0.242%

Cd Variation C1 Variation

0.0110

0.0043 60.91% 0.6506 2.371%

0.0057 48.18% 0.6444 3.301%

0.0038 65.45% 0.6669 0.075%

0.0043 60.91% 0.6700 0.540%

0.0047 57.27% 0.6670 0.090%

0.0049 55.45% 0.6674 0.150%

Thickness Variation Area Variation

0.1211

0.1162 4.069% 0.0757 2.800%

0.1207 0.329% 0.0779 0.054%

0.1145 5.486% 0.0748 3.904%

0.1184 2.261% 0.0779 0.021%

0.1209 0.127% 0.0790 1.458%

0.1206 0.411% 0.0779 0.059%

0.7028

0.6664

0.0779

−3

6

x 10

Pareto Front NE (Γc1 is odd, Γc2 is even)

Cd2

5.5

NE (Γc1 is even, Γc2 is odd)

5 4.5 4 3.5 3 2.5 2 1.5 1.5

2

2.5

3

3.5 Cd1

4 −3

x 10

Fig. 6. Pareto front and Nash equilibria of CNA design case1.

6. Conclusion In this paper, a constraint handling method is developed to apply Nash games to realistic multi objective aerodynamic optimization with constraints. A competitive Nash game for constrained multi-objective optimization is converted into a cooperative game, in which Nash competitive game acts as a sub-game for calculating the equilibrium among objectives. Then it further collaborates with players in charge of constraints. The partial gradient of each objective instead of design variables is used as elitist information exchange in present Nash algorithm. Existence and equivalence of a solution of the present CNA algorithms are analyzed and proved in the paper. Numerical experiments on 2-D two-objective aerodynamic optimization with lift and/or geometry constraints are implemented, and numerical results show that constraints are accurately satisfied. The algorithm developed in this paper can be extended easily to 3-D applications by replacing the 2-D adjoint flow analyzer by its 3-D extension. References [1] J.F. Nash, Non-cooperative games, Ann. Math. 54 (2) (1951) 286–295. [2] Z. Tang, J.A. Désidéri, J. Périaux, Multi-criteria aerodynamic shape-design optimization and inverse problems using control theory and Nash games, J. Optim. Theory Appl. 135 (1) (2007) 599–622.

Z. Tang et al. / Comput. Methods Appl. Mech. Engrg. 271 (2014) 130–143

143

[3] Z. Tang, W. Bai, J. Dong, Distributed optimization using virtual and real game strategies for multi-criterion aerodynamic design, Sci. China Ser. E: Technol. Sci. 51 (11) (2008) 1939–1956. [4] A. Habbal, J. Petersson, M. Thellner, Multidisciplinary topology optimization solved as a Nash game, Int. J. Numer. Methods Eng. 61 (2004) 949–963. [5] S. Li, Y. Zhang, Q. Zhu, Nash-optimization enhanced distributed model predictive control applied to the Shell benchmark problem, Inf. Sci. 170 (2005) 329–349. [6] K. Deb, An efficient constraint handling method for genetic agorithms, Comput. Methods Appl. Mech. Eng. 186 (2–4) (2000) 311–338. [7] R. Farmani, JA. Wright, Self-adaptive fitness formulation for constrained optimization, IEEE Trans. Evol. Comput. 7 (5) (2005) 445–455. [8] Y.G. Woldesenbet, G.G. Yen, Constraint handling in multiobjective evolutionary optimization, IEEE Trans. Evol. Comput. 13 (3) (2009) 514–525. [9] F. Huang, L. Wang, Q. He, An effective co-evolutionary differential evolution for constrained optimization, Appl. Math. Comput. 186 (1) (2007) 340–356. [10] T.P. Runarsson, X. Yao, Stochastic ranking for constrained evolutionary optimization, IEEE Trans. Evol. Comput. 4 (3) (2000) 284–294. [11] T. Takahama, S. Sakai, Constrained optimization by applying the a constrained method to the nonlinear simplex method with mutations, IEEE Trans. Evol. Comput. 9 (5) (2005) 437–451. [12] N. Poon, J. Martins, An adaptive approach to constraint aggregation using adjoint sensitivity analysis, Struct. Multi. Optim. 34 (2007) 61–73. [13] S. Venkatraman, G.G. Yen, A generic framework for constrained optimization using genetic algorithms, IEEE Trans. Evol. Comput. 9 (4) (2005) 424–435. [14] Z. Cai, Y. Wang, A multiobjective optimization-based evolutionary algorithm for constrained optimization, IEEE Trans. Evol. Comput. 10 (6) (2006) 658–675. [15] T. Ray, K.M. Liew, Society and civilization: an optimization algorithm based on the simulation of social behavior, IEEE Trans. Evol. Comput. 7 (4) (2003) 386–396. [16] Y. Wang, Z. Cai, Y. Zhou, Z. Fan, Constrained optimization based on hybrid evolutionary algorithm and adaptive constraint-handling technique, Struct. Multi. Optim. 37 (2009) 395–413. [17] G. Debreu, A social equilibrium existence theorem, Proc. Nat. Acad. Sci. 38 (10) (1952) 886–893. [18] T. Ichiishi, Game Theory for Economic Analysis, Academic Press, New York, 1983. [19] S.M. Robinson, Shadow prices for measures of effectiveness. II. General model, Oper. Res. 41 (3) (1993) 536–548. [20] P.T. Harker, Generalized Nash games and quasivariational inequalities, Eur. J. Oper. Res. 54 (1991) 81–94. [21] J. Krawczyk, Numerical solutions to coupled-constraint (or generalised Nash) equilibrium problems, CMS 4 (2007) 183–204, http://dx.doi.org/10.1007/ s10287-006-0033-9. [22] A. Haurie, J.B. Krawczyk, Optimal charges on river effluent from lumped and distributed sources, Environ. Model. Assess. 2 (3) (1997) 93–106. [23] J.B. Krawczyk, S. Uryasev, Relaxation algorithms to find Nash equilibria with economic applications, Environ. Model. Assess. 5 (2000) 63–73. [24] F. Facchinei, C. Kanzow, Generalized Nash equilibrium problems, Ann. Oper. Res. 175 (2010) 177–211. [25] A. Ehrenmann, Equilibrium problems with equilibrium constraints and their application to electricity markets, Ph.D. Dissertation, Judge Institute of Management, The University of Cambridge, Cambridge, 2004. [26] M. Fukushima, J.S. Pang, Quasi-variational inequalities, generalized Nash equilibria, and multi-leaderfollower games, Comput. Manage. Sci. 2 (2005) 21–56. [27] X. Hu, D. Ralph, Using EPECs to model bilevel games in restructured electricity markets with locational prices, Technical Report CWPE 0619, 2006. [28] S. Leyffer, T. Munson, Solving multi-leader-follower games, Argonne National Laboratory Preprint ANL/MCS-P1243-0405, Illinois, 2005. [29] J. Zhou, W.H.K. Lam, B.G. Heydecker, The generalized Nash equilibrium model for oligopolistic transit market with elastic demand, Transp. Res. B 39 (2005) 519–544. [30] G. Arslan, MF. Demirkol, S. Yuksel, On games with coupled constraints, in: 51st IEEE Conference on Decision and Control, IEEE, 2012, USA, pp. 7151– 7156. [31] E. Altman, E. Solan, Constrained Games: the impact of the attitude to adversary’s constraints, IEEE Trans. Autom. Control 54 (10) (2009) 1–6. [32] M.J. Osborne, A. Rubinstein, A course in game theory, Technical report, Massachusetts Institute of Technology, 1997. [33] Z. Tang, J. Dong, Couplings in multi-criterion aerodynamic optimization problems using adjoint methods and game strategies, Chinese J. Aeronaut. 22 (2009) 1–8. [34] Z. Tang, Multi-objective optimization strategies using adjoint method and game theory in aerodynamics, Acta Mech. Sin. 22 (2006) 307–314. [35] H. Samelson, On the Brouwer fixed point theorem, Port. Math. 22 (1963) 189–191. [36] J.A. Désidéri, Split of territories in concurrent optimization, INRIA Research Report No. 6108, 2007. [37] A. Jameson, Aerodynamic design via control theory, J. Sci. Comput. 3 (3) (1988) 233–260. [38] Z. Tang, The research on optimum aerodynamic design using CFD and control theory, PhD thesis of NUAA, P.R. China, 2000. [39] Z. Tang, J. Périaux, Uncertainty based robust optimization method for drag minimization problems in aerodynamics, Comput. Methods Appl. Mech. Eng. 217–220 (2012) 12–24. [40] T. Kipouros, Stochastic Optimisation in Computational Engineering Design, Adv. Intell. Soft Comput. 175 (2012) 475–490. [41] T. Ghisu, G.T. Parks, J.P. Jarrett, P.J. Clarkson, An integrated system for the aerodynamic design of compression systems: part II: application, ASME J. Turbomach. 133 (1) (2011), http://dx.doi.org/10.1115/1.4000535. [42] V. Lattarulo, T. Kipouros, G.T. Parks, Application of the multi-objective alliance algorithm to a benchmark aerodynamic optimization problem, in: IEEE Congress on Evolutionary Computation, IEEE, 2013, pp. 3182–3189.