Marco Montemurro

Angela Vincenti

Professor, Paolo Vannucci
email: paolo.vannucci@uvsq.fr

The Automatic Dynamic Penalisation (ADP) Constraint-Handling Technique for Engineering Optimisation Problems

research documents, whether they are published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

During the last four decades, an increasing interest in problem-solving strategies based on the principles of evolution and heredity has emerged. Among these systems we can include Evolution Algorithms (EAs) [START_REF] Rechenberg | Evolutionsstrategie: Optimierung technischer Systeme nach Prinzipien der biologischen Evolution[END_REF][START_REF] Schwefel | Numerical optimization for computer models[END_REF][START_REF] Michalewicz | Genetic algorithms + data structures = evolution programs[END_REF]], Fogel's Evolutionary Programming [START_REF] Fogel | Artificial intelligence through simulated evolution[END_REF], Glover's Scatter Search techniques [START_REF] Glover | Heuristics for integer programming using surrogate constraints[END_REF]. Besides these techniques, the most popular and well-known evolution-based strategies are Genetic Algorithms (GAs) [START_REF] Holland | Adaptation in natural and artificial systems[END_REF][START_REF] Goldberg | Genetic algorithms[END_REF]. A useful and common term often used for all the evolution-based systems cited beforehand is Evolution Programs (EPs), see [START_REF] Michalewicz | Genetic algorithms + data structures = evolution programs[END_REF].

The reason that underlies the success of such approaches is that traditional non-linear programming methods are local search techniques which tend to converge towards local optima. This is a serious drawback in many practical structural optimisation problems where the design space could be non-convex or even disjoint. In addition, such techniques are very popular due to their ability to naturally handle discrete design variables.

However, the main drawback of EPs concerns the treatment of constraints. GAs, as well as other EPs, do not take into account constraints directly. This means that standard genetic operators, such as mutation, crossover and so on, are "blind" with respect to constraints. Consequently, despite the parents satisfy some constraints, they can generate offspring which violate them. Technically, this means that GAs (and, more generally, EPs) perform an unconstrained search. Therefore, the major concern for researchers in this field consists in developing more and more accurate strategies in order to incorporate constraints directly into the objective and/or the fitness function.

Several authors put an effort in developing appropriate and effective strategies, in the framework of EPs, in order to deal with constrained optimisation problems. A certain number of surveys on constraint-handling techniques is available in the specialised literature, see for example [START_REF] Michalewicz | Genetic algorithms + data structures = evolution programs[END_REF][START_REF] Michalewicz | Evolutionary Algorithms in Engineering Applications[END_REF][START_REF] Michalewicz | Evolutionary algorithms for constrained parameter optimization problems[END_REF][START_REF] Coello Coello | Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art[END_REF].

The most common approach to handle constraints (particularly, inequality constraints) in the EP community consists in using penalties. Penalty functions were originally proposed by Courant in the 1940s [START_REF] Courant | Variational methods for the solution of problems of equilibrium and vibrations[END_REF] and later generalised by Carroll [START_REF] Carroll | The created response surface technique for optimizing nonlinear restrained systems[END_REF] and Fiacco and McCormick [15]. The idea that underlies such approaches consists in transforming the constrained optimisation problem into an unconstrained one by adding some given values to the objective function, based on the amount of constraint violation, for each considered point within the search space.

Penalty-based method are well suited for a wide class of optimisation problems, but it is difficult to find a proper penalty factor when balancing between feasible and infeasible solutions [START_REF] Coello Coello | Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art[END_REF]. To avoid such an issue, Runarsson and Yao [START_REF] Runarsson | Stochastic ranking for constrained evolutionary optimization[END_REF] developed a stochastic ranking technique wherein the objective values were used to rank the solutions within the infeasible region of the search space. A probability parameter was used to determine the fitting of two individuals belonging to the infeasible space, when compared with each other. This method was proved to be effective in solving a wide range of constrained optimisation problems, but the quality of the final solution depends upon the choice of the probability parameter. Based on the value of the objective function of both infeasible and feasible solutions, Farmani and Wright [START_REF] Farmani | Self-adaptive fitness formulation for constrained optimization[END_REF] proposed a two-stage self-adaptive fitness formulation method without the need of tuning any parameter. Lin and Wu [START_REF] Lin | Self-organizing adaptive penalty strategy in constrained genetic search[END_REF] proposed a self-organising adaptive penalty strategy having the aim of providing a parameter-free method with self-adaptive selection pressure and, simultaneously, an automatic normalisation of the constraints. Huang et al. [START_REF] Huang | An effective co-evolutionary differential evolution for constrained optimization[END_REF] proposed a co-evolutionary differential evolution for constrained optimisation through the use of two populations evolving simultaneously in order to update penalty factors and decision solutions interactively.

Repair methods are also widely employed as constraint-handling techniques. The original focus of such methods consists in transforming an infeasible solution into a feasible one through specially dedicated operator.

On one hand, repair methods ensure the feasibility of the solutions, but on the other hand they introduce a strong bias in the search process [START_REF] Coello Coello | Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art[END_REF]. To overcome such an issue, Orvosh and Davis [START_REF] Orvosh | Shall we repair? genetic algorithms, combinatorial optimization and feasibility constraints[END_REF][START_REF] Orvosh | Using a genetic algorithm to optimize problems with feasibility constraints[END_REF] developed the so-called "5% rule" for combinatorial optimisation, where 5% of the repaired chromosomes were used to substitute the infeasible original ones. Chootinan and Chen [START_REF] Chootinan | Constraint handling in genetic algorithms using a gradient-based repair method[END_REF] proposed a gradient-based repair operator for continuous optimisation, where the infeasible solution was repaired with a certain probability. Takahama and Sakai [START_REF] Takahama | Constrained optimization by the a constrained differential evolution with gradient-based mutation and feasible elites[END_REF] used the gradient-based mutation to change the infeasible solutions with a given probability.

Conceiving special selection operators represents a sound alternative to handle constraints. Deb [START_REF] Deb | An efficient constraint handling method for genetic algorithms[END_REF] proposed a pair-wise comparison in tournament selection without any penalty parameter. In particular, in the context of such a selection strategy, feasible solutions are always better than infeasible ones, the comparison between infeasible solutions depends on the amount of constraints violation, and comparison between feasible solutions depends on the corresponding objectives values. Wang et al. [START_REF] Wang | Constrained optimization based on hybrid evolutionary algorithm and adaptive constraint handling technique[END_REF] proposed an adaptive constrainthandling technique based on the number of feasible solution restrained in the population. More recently, Wang and Li [START_REF] Wang | An effective differential evolution with level comparison for constrained engineering design[END_REF] proposed an effective differential evolution with a level comparison, where the single level is dynamically updated to exploit the information coming from infeasible solutions at the beginning of the search process and forcing the feasibility at later stages.

Recently, multi-objective optimisation techniques have been employed to deal with constrained optimisation problems. The main idea of such methods is to treat the constraints as objective functions, thus the constraints and the original objective are optimised simultaneously. Surry and Radcliffe [START_REF] Surry | The COMOGA method: constrained optimisation by multi-objective genetic algorithms[END_REF] considered the constrained optimisation problem firstly as a constrained satisfaction problem (neglecting at this stage the objective function) and secondly as an unconstrained optimisation problem, by ignoring the constraints.

Venter and Haftka [START_REF] Venter | Constrained particle swarm optimization using a biobjective formulation[END_REF] treated all the constraints as an equivalent single-objective function and solved the optimisation as a bi-objective optimisation problem by means of a specialised bi-objective particle-swarm optimisation algorithm.

In this paper we propose a new penalty-based approach, i.e. the Automatic Dynamic Penalisation (ADP) technique, for constrained optimisation problems. Such a method belongs to the category of exterior penaltybased strategies. Our aim consists in providing a simple and effective parameter-free constraint-handling technique: indeed, with the ADP strategy we do not need the intervention of the user to tune the penalty coefficients values for any given optimisation problem. The ADP strategy is implemented within the GA BIANCA which is a multi-population GA for engineering optimisation problems which allows the parallel evolution of individuals belonging to different species, thanks to the introduction of some special genetic operators for crossover and mutation among species. For more details on the GA BIANCA the reader is addressed to [START_REF] Vincenti | BIANCA: a genetic algorithm to solve hard combinatorial optimisation problems in engineering[END_REF][START_REF] Montemurro | A two-level procedure for the global optimum design of composite modular structures -application to the design of an aircraft wing. part 1: theoretical formulation[END_REF][START_REF] Montemurro | A two-level procedure for the global optimum design of composite modular structures -application to the design of an aircraft wing. part 2: numerical aspects and examples[END_REF].

The paper is organised as follows: firstly the ADP technique is described in Sec. 2, then in Sec. 3 the ADP strategy is applied to three well-known benchmark problems in order to show its robustness and effectiveness and the obtained results are then compared to those available in literature. Finally, in Sec 4 two examples of real-word engineering applications are considered (i.e. the optimal design of damping properties of hybrid elastomer/composite laminates and the maximisation of the first buckling load of composite laminates with given elastic symmetries) in order to show the effectiveness of the proposed approach when dealing with more complex constrained optimisation problems.

The Automatic Dynamic Penalisation (ADP) strategy

Let us consider a constrained optimisation problem stated as follows:

min x Φ (x) , subject to :        g i (x) ≤ 0 i = 1, ..., r , h j (x) = 0 j = 1, ..., m , x L ≤ x ≤ x U , ( 1
)
where vectors and matrix terms are marked in bold typeface. In this formulation x is the n-dimensional vector of design variables, while x L and x U are the n-dimensional vectors representing the lower and upper bounds of the design variables, i.e. the whole definition domain or search space. The full set of equality and inequality constraints along with the box constraints represents the feasible domain or design space Ω . Design variables can be of different type: continuous, regular discrete, scattered (i.e. discrete variables without a discretisation step) or "grouped", these last being a sort of "abstract" variables representing a group of different variables, such as, for example, in the case of the constitutive material of a structure, when the material is chosen within a database: once a particular material is associated to a part, the whole set of the properties of the material are determined, i.e. elastic moduli, mass density and so on.

The goal of the optimisation consists in minimising the objective function Φ (x) subject to a given number of constraints: g i (x) (i = 1, ..., r) are the r functions of inequality constraints, while h j (x) (j = 1, ..., m) are the m functions of equality constraints. For any solution x in the feasible domain Ω , all equality constraints h j (x) = 0, are active at all points of Ω , and an inequality constraint that satisfies g i (x) = 0, is said an active constraint.

An optimisation problem can be characterised either by the type of constraints in the problem formulation or by the linearity or non-linearity of the objective and constraint functions. A problem where at least some of the objective and constraint functions are non-linear is called non-linear programming problem (NLPP).

In the framework of the penalty-based approaches the constrained NLPP of Eq.( 1) is transformed into an unconstrained one, by defining a new modified objective function as follows:

min x Φ p (x) ,
where :

Φ p (x) =                        Φ (x) if g i (x) ≤ 0 and h j (x) = 0 , Φ (x) + r ∑ i=1 c i G i (x) if g i (x) > 0 and h j (x) = 0 , Φ (x) + m ∑ j=1 q j H j (x) if g i (x) ≤ 0 and h j (x) ̸ = 0 , Φ (x) + r ∑ i=1 c i G i (x) + m ∑ j=1 q j H j (x) if g i (x) > 0 and h j (x) ̸ = 0 , i = 1, ..., r , j = 1, ..., m , (2)
where Φ p is the penalised (or expanded) objective function, while c i and q j are the penalty coefficients for inequality and equality constraints, respectively. The quantities G i (x) and H j (x) are defined as:

G i (x) = max{0, g i (x)} i = 1, ..., r , H j (x) = max{0, |h j (x)| -ϵ} j = 1, ..., m . (3)
It can be noticed that in Eq.( 2) and (3) equality constraints were transformed (as normally done in the literature) into inequality ones having the form:

|h j (x)| -ϵ ≤ 0 , (4)
which is numerically "acceptable" until the allowed tolerance ϵ assumes a sufficient small value that does not affect the quality of the solution of the problem. The advantage is that at this point all constraints are of the same nature (inequalities) and they can be treated by the same technique of constraint-handling.

As specified in [START_REF] Coello Coello | Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art[END_REF], the penalty should be kept as low as possible, just above the limit below which infeasible solutions are optimal (this is called, the minimum penalty rule [START_REF] Davis | Genetic Algorithms and Simulated Annealing[END_REF]). This is due to the fact that if the penalty is too high or too low, then the GA might have some difficulties in finding an appropriate feasible optimal solution [START_REF] Michalewicz | Genetic algorithms + data structures = evolution programs[END_REF][START_REF] Davis | Genetic Algorithms and Simulated Annealing[END_REF].

In fact, on one hand, if the penalty is too high and the optimum point is placed on the boundary of the feasible region, the GA is pushed inside the feasible region very quickly, and it is not able to move back towards the boundary between the feasible and infeasible regions, i.e. a large penalty discourages the exploration of the infeasible region since the beginning of the search process. On the other hand, if the penalty is too low, a lot of the search time is spent in exploring the infeasible region because the penalty term is often negligible with respect to the objective function and the algorithm might converge to an optimum outside the feasible domain. These issues are very important in GAs, because in several problems the optimal solution lies close to or on the boundary between the feasible and the infeasible regions.

The ADP strategy is an original method, firstly presented in [START_REF] Vincenti | BIANCA: a genetic algorithm to solve hard combinatorial optimisation problems in engineering[END_REF], that we have developed and generalised in this work for automatically choosing and updating the penalty coefficients. The basic idea is that some infeasible individuals can be anyway important to drive the exploration towards interesting zones of the feasible domain, namely when the optimum point lies on its boundary, i.e. on an active constraint. For this reason, in the context of the ADP strategy, infeasible points are not automatically excluded from the population and are used to dynamically update the penalty coefficients in an automatic way, i.e. without the intervention of the user. This is especially important at the early stages of the search in order to widely explore the whole search space.

Concerning the penalty coefficients c i and q j of Eq.( 2), in classical penalty-based methods, the user have to properly set their values in order to ensure that the search of solutions is forced within the feasible domain.

Nevertheless, the choice of these coefficients is very difficult and it is common practice to estimate their values by trial and error. Moreover, it could be useful to adjust penalty pressure along the generations by tuning these coefficients, but this is directly linked to a deep knowledge of the nature of the optimisation problem at hand.

The main concept that underlies the ADP strategy is that it is possible to exploit the information restrained in the whole population (also the infeasible part of it), at the current generation, to better guide the search through the whole search space. Generally, at the first generation the population is randomly generated. The individuals are more or less uniformly distributed over both feasible and infeasible regions of the definition domain and the corresponding values of objective and constraints functions can be used to estimate an appropriate level of penalisation, i.e. the values of the penalty coefficients c i and q j .

At each current generation, inside the population we can separate feasible and infeasible individuals, see Fig. 1, and we can also classify each one of the two groups in terms either of the values of the objective function or of the amount of constraint violation. Thus, the best individual of each group is the potential candidate to be solution of the optimisation problem on the feasible and infeasible side of the domain, respectively. Then, we choose the following definition of the penalty coefficients:

c i (t) = |Φ F best -Φ N F best | (G i ) N F best i = 1, ..., r , q j (t) = |Φ F best -Φ F best | (H j ) N F best j = 1, ..., m .
(5)

In Eq.( 5) the coefficients c i and q j are evaluated at the current generation t, while the superscripts F and N F stand for feasible and non-feasible, respectively. Φ F best and Φ N F best are the values of the objective function for the best individuals within the feasible and the infeasible side of the domain, respectively, whilst (G i ) N F best and (H i ) N F best represent the violated inequality and equality constraints for the best infeasible solution. The main reason that underlies the definition of Eq.( 5) is that some performing infeasible solutions (in terms of objective function values, i.e. infeasible minima) are retained within the population and act as "attraction points", improving in this way the search properties of the GA. In other words, the presence of such infeasible points within the population improves the exploration of the whole search domain, particularly around the boundaries between feasible and infeasible sides. In particular, it is worth noting that substituting Eq.( 5) into Eq.( 2), the value of the objective function for the best infeasible individual is forced to be equal to that of the best feasible individual. Such a situation is very convenient, mostly for what concerns the case of non-linear, non-convex optimisation problems wherein the global unconstrained minimum is located within the infeasible region and, potentially, the global constrained minimum lies on the boundary. In such conditions, the best individual of the infeasible side owns the same fitness as the best feasible individual, and thus the same probability to be selected in order to take part into the genetic operations of crossover and mutation. As a consequence, its genetic patrimony is likely to be inherited by the following generation.

Consequently, the GA can handle the information coming from the infeasible region in order to drive the search towards more convenient sub-domains (namely towards the boundaries between feasible and infeasible sides).

It appears clearly that the estimation of the penalty factors, according to the Eq.( 5), can be repeated at each generation, thus tuning the appropriate penalty pressure on the current population. The main advantages of such an approach are substantially two:

• this procedure is automatic and problem-independent because the GA can automatically calculate the values of the penalty coefficients without the intervention of the user by simply exploiting the values of the objective and constraint functions in the current population;

• the method is dynamic since the evaluation of the penalty level is updated at each generation, and this allows the values of the penalty coefficients to be the most suitable to the current distribution of feasible and infeasible individuals in the population, the expected effect being eventually to extinguish the infeasible group in the population or to limit infeasible individuals to regions close to the boundary between feasible and infeasible domains. 4. the penalty coefficient c k (or q k ) is then evaluated according to Eq.( 5).

• Case 2 : no feasible individuals (with respect to the k th constraint function) within the current population 5. the penalty coefficient c k (or q k ) is then evaluated according to Eq.( 5);

It is worth noting that, in the context of the ADP strategy, each constraint is treated separately and independently from each other.

To understand the way whereby the ADP strategy acts on the individuals within the population, and also to show its effectiveness, let us consider the following optimisation problem:

min x Φ (x 1 , x 2 ) = -e ka √ x1+x2 sin (ax 1 ) cos (2bx 2 ) , subject to :        g (x 1 , x 2 ) = e cx 2 1 -1 -x 2 ≤ 0 , 0 ≤ x 1 ≤ 4π , 0 ≤ x 2 ≤ 2π , ( 6
)
where

Benchmark Problems

To have an idea of the effectiveness of the ADP strategy, we consider three benchmark problems belonging to the engineering world, which were extensively studied in the literature.

Such benchmark problems are: the welded beam problem, originally proposed by Rao [START_REF] Rao | Engineering Optimization[END_REF], the pressure vessel problem, firstly studied by Kannan and Kramer [START_REF] Kannan | An augmented lagrange multiplier based method for mixed integer discrete continuous optimization and its applications to mechanical design[END_REF], and the tension-compression spring problem taken from Arora [START_REF] Arora | Introduction to Optimum Design[END_REF] and Belugundu [START_REF] Belegundu | A study of mathematical programming methods for structural optimization[END_REF].

Concerning the genetic parameters employed in all the simulations, the population is composed of N ind individuals evolving along a fixed maximum number of generation N gen . For each considered benchmark we perform 30 runs of our GA, and, for each run, 80000 fitness evaluations are carried out. This implies that such a number of fitness evaluations can be obtained with various combinations of number of individuals and maximum number of generations, i.e. the parameters N ind and N gen must satisfy the following relationship:

N gen × N ind = 80000 . ( 7
)
In addition, the crossover and mutation probabilities are p cross = 0.85 and p mut = 1/N ind , respectively. The selection is performed by the roulette-wheel operator, the single-individual elitism is active and, of course, the ADP method was used for handling constraints.

For each test-case, we compare our results with the ones reported in the literature and obtained by other researchers by the use of several EA-based methods which employ different constraint-handling techniques.

The welded beam problem

The welded beam problem was originally studied by Rao [START_REF] Rao | Engineering Optimization[END_REF]. In such a problem a welded beam is designed for minimum cost subject to constraints on shear stress (τ ), bending stress in the beam (σ), buckling load on the bar (P c ), end deflection of the beam (δ), and side constraints. There are four design variables as shown in Fig. 4 :

x 1 = h, x 2 = l, x 3 = t and x 4 = b.
The problem can be stated as follows:

min x Φ (x) = 1.10471x 2 1 x 2 + 0.04811x 3 x 4 (14.0 + x 2 ) , subject to :                                                        g 1 (x) = τ (x) -τ max ≤ 0 , g 2 (x) = σ (x) -σ max ≤ 0 , g 3 (x) = x 1 -x 4 ≤ 0 , g 4 (x) = 0.10471x 2 1 + 0.04811x 3 x 4 (14.0 + x 2 ) -5.0 ≤ 0 , g 5 (x) = 0.125 -x 1 ≤ 0 , g 6 (x) = δ (x) -δ max ≤ 0 , g 7 (x) = P -P c (x) ≤ 0 , 0.1 ≤ x 1 ≤ 2.0 , 0.1 ≤ x 2 ≤ 10.0 , 0.1 ≤ x 3 ≤ 10.0 , 0.1 ≤ x 4 ≤ 2.0 , (8)
where:

τ (x) = √ τ 2 1 + 2τ 1 τ 2 x 2 2R + τ 2 2 , τ 1 = P √ 2x 1 x 2 , τ 2 = M R J , M = P (L + x 2 2 ) , R = √ ( x 2 2 ) 2 + ( x 1 + x 3 2 ) 2 , J = 2 { √ 2x 1 x 2 [ x 2 2 12 + ( x 1 + x 3 2 ) 2 ]} , σ (x) = 6 P L x 4 x 2 3 , δ (x) = 4P L 3 Ex 4 x 3 3 , P c (x) = 4.013E √ x 2 3 x 6 4 /36 L 2 ( 1 - x 3 2L √ E 4G ) . ( 9
)
The quantities P , L, δ max , τ max , σ max , E and G are constant parameters and their numerical values, in the appropriate units [START_REF] Coello Coello | Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art[END_REF], are: P = 6000, L = 14, δ max = 0.25, τ max = 13600, σ max = 30000, E = 30 × 10 6 and G = 12 × 10 6 .

Figure 4: Rough sketch of the welded beam considered within the optimisation problem (8), taken from [START_REF] Coello Coello | Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art[END_REF] The approaches employed to deal with such a problem are: geometric programming (Ragsdell and Phillips [START_REF] Ragsdell | Optimal design of a class of welded structures using geometric programming[END_REF]), standard GA with static penalty function (Deb [11]), an improved GA with a feasibility-based tournament selection scheme inspired by the multi-objective optimisation techniques (Coello and Montes [START_REF] Coello | Constraint-handling in genetic algorithms through the use of dominance-based tournament selection[END_REF]), a co-evolutionary particle-swarm strategy for constrained optimisation problems (He and Wang [START_REF] He | An effective co-evolutionary particle swarm optimization for constrained engineering design problems[END_REF]), a hybrid particle-swarm optimisation with feasibility rules (He and Wang [START_REF] He | A hybrid particle swarm optimization with a feasibility-based rule for constrained optimization[END_REF]) and a hybrid GA with flexible allowance technique (Zhao et al. [START_REF] Zhao | An effective hybrid genetic algorithm with flexible allowance technique for constrained engineering design optimization[END_REF]).

Concerning the GA BIANCA, in order to satisfy the condition on the maximum number of fitness evaluations of Eq.( 7), we consider a population of N ind = 250 individuals evolving along N gen = 320 generations.

The best solution found by BIANCA as well as the best solutions obtained by the aforementioned approaches are listed in Table 1, while the statistical results for each considered strategy are detailed in Table 2. Fig. 5 shows the variation of the best solution along the generations. It can be noticed that the global minimum is found after 260 generations. Ragsdell and et al. [START_REF] Zhao | An effective hybrid genetic algorithm with flexible allowance technique for constrained engineering design optimization[END_REF] Wang [START_REF] He | A hybrid particle swarm optimization with a feasibility-based rule for constrained optimization[END_REF] Wang [START_REF] He | An effective co-evolutionary particle swarm optimization for constrained engineering design problems[END_REF] Montes [START_REF] Coello | Constraint-handling in genetic algorithms through the use of dominance-based tournament selection[END_REF] Phillips [

The pressure vessel problem

A cylindrical vessel is capped at both ends by hemispherical heads as shown in Fig. 6. The goal of this problems consists in minimising the total cost of the structure, including the cost of the material, forming and welding. The design variables are: the thickness of the shell x 1 = T s , the thickness of the head x 2 = T h , the inner radius x 3 = R and the length of the cylindrical section of the vessel (not including the head) x 4 = L.

Moreover, T s and T h are real discrete design variables discretised with a precision ∆x 1 = ∆x 2 = 0.0625 (which corresponds to the available thickness of rolled steel plates), while R and L are continuous.

Using the same notation adopted by Kannan and Kramer [START_REF] Kannan | An augmented lagrange multiplier based method for mixed integer discrete continuous optimization and its applications to mechanical design[END_REF], the problem can be stated as follows: The approaches used in the literature to solve this problem are: an augmented Lagrangian multiplier approach (Kannan and Kramer [START_REF] Kannan | An augmented lagrange multiplier based method for mixed integer discrete continuous optimization and its applications to mechanical design[END_REF]), a genetic adaptive search (Deb [START_REF] Deb | GeneAS: a robust optimal design technique for mechanical component design[END_REF]), and, again, the aforementioned approaches of Coello and Montes [START_REF] Coello | Constraint-handling in genetic algorithms through the use of dominance-based tournament selection[END_REF], He and Wang [START_REF] He | An effective co-evolutionary particle swarm optimization for constrained engineering design problems[END_REF][START_REF] He | A hybrid particle swarm optimization with a feasibility-based rule for constrained optimization[END_REF] and Zhao et al. [START_REF] Zhao | An effective hybrid genetic algorithm with flexible allowance technique for constrained engineering design optimization[END_REF].

min x Φ (x) = 0.6224x 1 x 3 x 4 + 1.7781x 2 x 2 3 + 3.1661x 4 x 2 1 + 19.84x 3 x 2 1 , subject to :                                        g 1 (x) = -x 1 + 0.0193x 3 ≤ 0 , g 2 (x) = -x 2 + 0.00954x 3 ≤ 0 , g 3 (x) = -πx 2 3 x 4 -
In order to satisfy the condition on the maximum number of fitness evaluations of Eq.( 7), the size of the population N ind and the maximum number of generations N gen in BIANCA are chosen equal to 400 and 200, respectively.

The best solution found by BIANCA as well as the best solutions obtained by the aforementioned approaches are listed in Table 3, while the statistical results for each considered strategy are detailed in Table 4. Fig. 7 shows the variation of the best solution along the generations. It can be noticed that the global minimum is found after only 20 generations.

The tension-compression spring problem

The design problem of a tension-compression spring was firstly studied by Arora [START_REF] Arora | Introduction to Optimum Design[END_REF] and Belegundu [START_REF] Belegundu | A study of mathematical programming methods for structural optimization[END_REF].

The main goal is to minimise the weight of the tension-compression spring (as shown in Fig. 8) subject to constraints on the minimum deflection, the shear stress, the surge frequency and the outside diameter. The design variables are: the wire diameter x 1 = d, the mean coil diameter x 2 = D and the number of active coils

x 3 = N coil .
Adopting the same notation of Arora [START_REF] Arora | Introduction to Optimum Design[END_REF], the problem can be stated as follows:

min x Φ (x) = (x 3 + 2) x 2 x 2 1 , subject to :                                            g 1 (x) = 1 - x 3 2 x 3 71785.0x 4 1 ≤ 0 , g 2 (x) = 4x 2 2 -x 1 x 2 12566.0(x 3 1 x 2 -x 4 1 ) + 1 5108.0x 2 1 -1 ≤ 0 , g 3 (x) = 1 - 140.45x 1 x 2 2 x 3 ≤ 0 , g 4 (x) = x 1 + x 2 1.5 -1 ≤ 0 , 0.05 ≤ x 1 ≤ 2.0 , 0.25 ≤ x 2 ≤ 1.3 , 2.0 ≤ x 3 ≤ 15.0 , ( 11
)
Figure 8: Rough sketch of tension-compression spring considered within the optimisation problem [START_REF] Deb | Optimal design of a welded beam via genetic algorithms[END_REF], taken from [START_REF] He | An effective co-evolutionary particle swarm optimization for constrained engineering design problems[END_REF] We compare our results with the ones carried out through the following strategies taken from the literature:

the numerical optimisation technique proposed by Belegundu [START_REF] Belegundu | A study of mathematical programming methods for structural optimization[END_REF], a numerical optimisation technique called "constraint correction at constant cost" (Arora [START_REF] Arora | Introduction to Optimum Design[END_REF]) and, again, the aforementioned approaches of Coello and

Montes [START_REF] Coello | Constraint-handling in genetic algorithms through the use of dominance-based tournament selection[END_REF], He and Wang [START_REF] He | An effective co-evolutionary particle swarm optimization for constrained engineering design problems[END_REF][START_REF] He | A hybrid particle swarm optimization with a feasibility-based rule for constrained optimization[END_REF] and Zhao et al. [START_REF] Zhao | An effective hybrid genetic algorithm with flexible allowance technique for constrained engineering design optimization[END_REF].

Concerning the GA BIANCA, in order to satisfy the condition on the maximum number of fitness evaluations of Eq.( 7), we consider a population of N ind = 320 individuals evolving along N gen = 250 generations.

The best solution found by BIANCA as well as the best solutions obtained by the aforementioned approaches are listed in Table 5, while the statistical results for each considered strategy are detailed in Table 6. Fig. 9 shows the variation of the best solution along the generations. It can be noticed that the global minimum is found after only 50 generations.

Discussion of results

Concerning the effectiveness and robustness of the GA BIANCA, we can see that for all the three considered benchmark problems, the quality of the results found using BIANCA is, practically, of the same order of that obtained via the hybrid strategies of He and Wang [START_REF] He | A hybrid particle swarm optimization with a feasibility-based rule for constrained optimization[END_REF] and Zhao et al. [START_REF] Zhao | An effective hybrid genetic algorithm with flexible allowance technique for constrained engineering design optimization[END_REF]. Indeed, the relative errors (evaluated with respect to the solutions found by Zhao et al. [START_REF] Zhao | An effective hybrid genetic algorithm with flexible allowance technique for constrained engineering design optimization[END_REF]) are 0.0034% for the welded beam design problem, 0.0004% for the pressure vessel design problem and 0.0047% for the tension-compression string design problem. In addition, one can notice that the average searching quality and the standard deviation of the results found by BIANCA in 30 independent runs are of the same order as the other methods (with the exception of the hybrid strategies).

Concerning the pressure vessel design problem, it can be noticed that the strategy of Kannan and Kramer [START_REF] Kannan | An augmented lagrange multiplier based method for mixed integer discrete continuous optimization and its applications to mechanical design[END_REF] produces a solution with a significantly lower value of L. This solution is, however, not feasible since the first constraint is slightly violated. The results produced by the other methods (including the ones found using BIANCA) indicate that it is more reasonable to variate the other design variables, allowing larger values of L because this leads to find feasible designs having a lower cost.

Concerning the tension-compression spring design problem, it can be noticed that the solution found by

Arora [START_REF] Arora | Introduction to Optimum Design[END_REF] has a lower value of the number of active coils N coil . Indeed, such a solution is infeasible since the first constraint is slightly violated. We can see that the results produced by the other approaches (including the ones found using BIANCA) indicate that allowing greater values of the number of active coils N coil leads the algorithm to find feasible designs showing a lower weight.

The hybrid GA with flexible allowance technique developed by Zhao et al. [START_REF] Zhao | An effective hybrid genetic algorithm with flexible allowance technique for constrained engineering design optimization[END_REF] is a hybrid GA with Levenberg-Marquardt mutation operator which creates new feasible individuals (offspring) by considering the auxiliary information coming from the evaluation of the constraints gradient. Such a strategy allows to obtain a numerically "exact" minimum for all the considered benchmarks. The authors assert that such solutions were found after only 20000 evaluations of the objective function (which corresponds to 1/4 of the number of evaluations carried out by the other considered methods).

Nevertheless, for such benchmark problems the objective and constraint functions are available in a closed algebraic form. Generally, this is not the case when dealing with optimisation of complex engineering systems (which require, for example, finite element calculations for the constraints and/or the objective function evaluation). In such cases, the results obtained via a hybrid GA, as the one proposed in [START_REF] Zhao | An effective hybrid genetic algorithm with flexible allowance technique for constrained engineering design optimization[END_REF], could be affected by the way wherein the derivatives of the constraints and objective function are evaluated (which represent a key-point for the Levenberg-Marquardt mutation operator). Moreover a more accurate analysis of the effectiveness and performances of such techniques (in terms of time spent to find a solution) could be made on the basis of the effective computing time rather than on the basis on the number of fitness evaluations, because the effective calculation cost includes also the evaluation of derivatives for the gradient-based operations of constraint repair.

Finally, we can assert that the main advantage of our approach relies on the fact that it remains a "purely" genetic approach: on one side, we do not need an estimation of auxiliary quantities, such as the derivatives of the constraint or objective functions, while on the other side, we practically have the same quality as the hybrid strategies in finding optimal solutions.

Two applications of the ADP strategy

In this section we discuss two examples of real-world engineering constrained optimisation problems to which the ADP strategy has been successfully applied. Such applications are: (a) the optimal design of damping properties of hybrid elastomer/composite laminates (see [START_REF] Montemurro | Design of damping properties of hybrid laminates through a global optimisation strategy[END_REF] for more details) and (b) the maximisation of the first buckling load of composite laminates with given elastic symmetries (see [START_REF] Vincenti | BIANCA: a genetic algorithm to solve hard combinatorial optimisation problems in engineering[END_REF] for more details). Both applications correspond to highly non-linear and non-convex constrained optimisation problems: in the first one, only inequality constraints are considered, whilst in the second problem also a difficult equality constraint is taken into account.

Further examples of real-world engineering applications whereto our strategy was successfully applied include the least-weight design of a stiffened wing-box section for an aircraft wing [START_REF] Montemurro | A two-level procedure for the global optimum design of composite modular structures -application to the design of an aircraft wing. part 1: theoretical formulation[END_REF][START_REF] Montemurro | A two-level procedure for the global optimum design of composite modular structures -application to the design of an aircraft wing. part 2: numerical aspects and examples[END_REF] and the inverse problem concerning the identification of electromechanical properties of composite laminated plates equipped with piezoelectric patches [START_REF] Montemurro | Identification of electromechanical properties of piezoelectric structures through evolutionary optimisation techniques[END_REF].

Optimal design of damping properties of hybrid elastomer/composite laminates

Problem description

In this first application we want to find a solution for the problem of designing the damping properties of hybrid elastomer/composite laminates and we applied our strategy to a rectangular hybrid plate, whose dimensions are depicted in Fig. 10.

Figure 10: Geometry of the hybrid plate.

Concerning the typical dimensions of the plate, the thickness of each layer is constrained to remain sufficiently small compared to both the width and the length of the plate, in order to keep valid the assumptions of the thin plate model. Moreover, we assume that the fibre-reinforced plies have linear elastic orthotropic behaviour. A glass-epoxy lamina is employed for the generic elastic ply, whose material properties (issued from [START_REF] Berthelot | Damping analysis of composite material and structures[END_REF]) are listed in Table 7.

Young The material used for the viscoelastic layers is a rubber-like material having linear isotropic behaviour.

In addition, the properties of that material are considered dependent upon the loading frequency f and are taken from [START_REF] Le Maoût | Simultaneous geometrical and material optimal design of hybrid elastomer/composite sandwich plates[END_REF]. The variation of the Young's modulus with the frequency is expressed as:

E(f ) = E s + E d log ( f f ) , ( 12
)
where E s = 0.0041 GPa is the steady-state value of the Young's modulus, E d = 0.0322 GPa is the amplitude of the part that depends upon the frequency, while f = 1 Hz is a reference value for the frequency. The Poisson's ratio and the material loss factor are kept constants and equal to ν = 0.3 and η v = 0.3, respectively, whilst the density is ρ = 968.1 kg m - 3 . More details about this problem can be found in [START_REF] Montemurro | Design of damping properties of hybrid laminates through a global optimisation strategy[END_REF].

Finite element model of the hybrid plate and loading conditions

In order to correctly predict the different mechanisms of dissipating the vibrational energy that underlie the viscoelastically damped structures (see [START_REF] Nakra | Vibration control in machines and structures using viscoelastic damping[END_REF]), 3D solid elements with 8 nodes and 3 degrees of freedom (DOFs) per node were employed to model either the elastic or the rubber layers. The FE model is built in the ANSYS environment.

The dynamic response of the structure is evaluated through a classical free-vibration analysis. Only the first N = 5 non-rigid modes are calculated considering free-displacement boundary conditions on the edges of the plate. It is worth noting that, since the material properties of the viscoelastic layers depend upon the frequency, the calculation of the eigenfrequencies (as well as that of the modal loss factors) is iterative for each eigenfrequency and was carried out using the well-known Iterative Modal Strain Energy (IMSE) method, see [START_REF] Montemurro | Design of damping properties of hybrid laminates through a global optimisation strategy[END_REF].

Once the convergence on the i th undamped natural frequency is reached, the corresponding modal loss factor η i is evaluated as:

η i = η v (f i ) W v (f i ) W tot (f i ) , ( 13
)
where η v (f i ) is the material loss factor at the current frequency, while W v (f i ) and W tot (f i ) are the strain energy of the viscoleastic layers and the total strain energy of the structure for the i th mode, respectively.

Along with the increase of damping capability, the structure must withstand the application of static loads, i.e. the structure has to exhibit good properties in terms of stiffness [START_REF] Montemurro | Design of damping properties of hybrid laminates through a global optimisation strategy[END_REF]. In particular we consider three different static loading conditions imposing a prescribed unitary displacement, δ = 1 mm, in some well-chose regions of the plate and we measure the in-plane and out-of-plane stiffness in terms of the reaction forces R x , R y and R z .

Mathematical statement of the damping properties design problem

The problem of designing the damping properties of a hybrid plate is stated as a constrained optimisation problem. The goal of our strategy consists in maximising the first N modal loss factors of the structure, without degrading the stiffness properties of the plate and without increasing too much its weight.

The problem is stated in the most general case, thus the design variables are:

• the total number of layers (both elastic and viscoelastic), n;

• the position and the number of the viscoelastic layers within the stack, which are directly linked to the variable ID k , (k = 1, ..., n), that identifies the nature of the k th ply, i.e. ID k = 1 if the k th ply is viscoelastic, ID k = 0 otherwise;

• the thickness of each layer, t k (k = 1, ..., n);

• the fiber orientation of the elastic plies, δ k (k = 1, ..., n).

The optimisation problem can now be formulated. The maximisation of the first N modal factors can be expressed as the minimisation of the following objective function:

Φ = - N ∑ i=1 η i ( 14
)
that represents the opposite of the sum of the first N modal loss factors. Moreover, the constraints on the maximum decrease of the stiffness properties and on the maximum increase of the mass of the plate have to be considered. Therefore, the constrained minimisation problem can be stated as a classical non-linear programming problem (NLPP) as follows:

min Φ (n, ID k , t k , δ k ) (with k = 1, ..., n) , subject to :                                g 1 (n, ID k , t k , δ k ) = R x ref -R x (n, ID k , t k , δ k ) R x ref -ϵ x ≤ 0 , g 2 (n, ID k , t k , δ k ) = R y ref -R y (n, ID k , t k , δ k ) R y ref -ϵ y ≤ 0 , g 3 (n, ID k , t k , δ k ) = R z ref -R z (n, ID k , t k , δ k ) R z ref -ϵ z ≤ 0 , g 4 (n, ID k , t k ) = M (n, ID k , t k ) -M ref M ref -ϵ M ≤ 0 . (15)
In Eq. ( 15) R x , R y and R z are the reactions of the plate, which represent a measure of the stiffness of the structure, while M is the mass of the plate. The superscript ref stands for reference value. The reference values of the reactions and mass are calculated, before the optimisation process, on a reference undamped structure, i.e. a plate without elastomeric layers. The quantities ϵ x , ϵ y , ϵ z and ϵ M are the user-defined tolerances on each constraint.

It can be noticed that the NLPP of Eq. ( 15) is highly non-linear and non-convex in the space of design variables (in fact, for a given eigenfrequency the strain energy, and hence the loss factor of the structure, depends upon circular functions of the plies orientation). In addition, we can notice that the objective as well as the constraint functions are not available in closed algebraic form, being the outcome of the FE calculations.

It is worth noting that, since the number of layers n is included among the optimisation variables, the total number of design variables of the whole optimisation process can change for each possible point-solution in the design space, or, in other words, the procedure determines by itself the optimal number of design variables.

Indeed, in the objective and constraint functions, we can enumerate three variables for each ply (the nature of the ply, ID k , the thickness t k and the orientation δ k ) besides the variable number n of layers. Therefore, in the most general case, the overall number of design variables for the problem ( 15) is 3n + 1.

The reference values for the reactions and the mass of the plate as well as those of the user-defined tolerances on the constraints of the problem (15) are issued from [START_REF] Montemurro | Design of damping properties of hybrid laminates through a global optimisation strategy[END_REF] and are evaluated considering a reference plate composed of 6 glass-epoxy laminae having a thickness equal to 0.3 mm and with the following stacking sequence: [90/45/0] S .

Numerical results

Here we show only the results concerning the most general case for which no simplifying assumptions on the position of the viscoelastic layers within the laminate stack are made.

The design variables, their nature and bounds are detailed in This set of values is problem-independent and demonstrated its validity over the large range of optimisation problems solved by BIANCA [START_REF] Vincenti | BIANCA: a genetic algorithm to solve hard combinatorial optimisation problems in engineering[END_REF][START_REF] Montemurro | A two-level procedure for the global optimum design of composite modular structures -application to the design of an aircraft wing. part 1: theoretical formulation[END_REF][START_REF] Montemurro | A two-level procedure for the global optimum design of composite modular structures -application to the design of an aircraft wing. part 2: numerical aspects and examples[END_REF][START_REF] Montemurro | Design of damping properties of hybrid laminates through a global optimisation strategy[END_REF][START_REF] Montemurro | Identification of electromechanical properties of piezoelectric structures through evolutionary optimisation techniques[END_REF] The best solution found by BIANCA is shown in Table 9. The optimal number of plies is 6. Moreover, for this case, the optimal configuration of the hybrid plate shows 2 viscoelastic plies on the top of the structure.

Comparing the plots of Fig. 11 (a) and 11 (b), one can notice something that happens systematically: convergence towards the best value of the number of modules (here, the number of plies) and of the objective function are independent. They never occur at the same time, and the optimisation of the number of modules happens always before that of the objective function. In other words, the strategy used in BIANCA for evolving simultaneously species and individuals, normally lets attain first the best species, and then continues to evolve individuals within the best species towards the best individual. In fact, as shown in Fig. 11, the global constrained minimum has been found after 62 generations, whilst the optimal number of plies is found after only 8 generations.

Fig. 12 shows the evolution of the species restrained in the whole population from the initial until the final generation. We can notice that in the first generation the number of individuals belonging to all the different species is equiprobable, whilst some species extinguish during the generations and cannot be found within the final population: only solutions with 5 to 8 plies can be found within the population in the last generation. in fact more and more new individuals are generated by the GA during the generations and this results in fitter individuals that gradually push the search process towards more convenient feasible regions of the whole search space. This is a very important fact: the ADP method seems to be a powerful strategy able to explore in an extensive way the whole search space, either from the feasible or the infeasible side, even if the initial population is made of infeasible individuals.

Maximisation of the first buckling load of composite laminates with given elastic symmetries

We consider the case of a rectangular simply-supported uncoupled orthotropic laminated plate, whose side lengths are a and b, subject to compression in-plane loads having resultant forces per unit length N x and N y , respectively, as shown in Fig. 13.

Figure 13: Rectangular laminated plate subject to in-plane forces per unit length N x and N y .

For such a case, buckling shapes are sinusoidal and, being m and n the number of half-waves along x and y directions, respectively, the critical value of the load multiplier λ crit (m, n) for the buckling load can be expressed as follows:

λ crit (m, n) = π 2 D xx ( m a ) 4 + 2 (D xy + 2D ss ) ( m a ) 2 ( n b ) 2 + D yy ( n b ) 4 N x ( m a ) 2 + N y ( n b ) 2 . ( 16
)
In Eq.( 16), D ij (i, j = x, y, s) are the Cartesian components of the bending stiffness tensor expressed in the global reference system of the plate R = {O; x, y, z}. It can be noticed that expression ( 16) can be applied only to orthotropic uncoupled plates, and therefore bending orthotropy as well as elastic uncoupling represent necessary conditions that the laminate must satisfy. In addition we want that the laminated plate satisfies the membrane orthotropy requirement too. Moreover, as tacitly assumed in Eq.( 16), the principal axes of orthotropy (for both membrane and bending behaviour) have to be coincident with the plate axes.

Along with the previous aspects, we look for a highly stiff orthotropic plate. Thus, constraints have to be imposed on the extension Young moduli along the orthotropy axes, namely E A x and E A y . Generally speaking, the problem of maximising the buckling load of a simply-supported laminated plate subject to the aforementioned constraints can be stated as:

max δ ( min m,n λ crit (δ) ) , subject to :            g 1 (δ) = E A x -E A x ≤ 0 , g 2 (δ) = E A y -E A y ≤ 0 , h 1 (δ) = I (P (δ)) = 0 . ( 17
)
In Eq.( 17), δ ∈ R N is the vector of layers orientations (the design variables of the present problem), while N is the number of plies. The objective function λ crit (δ) is highly non-linear and non-convex in the space of the layer orientations (this is due to the fact that the components of the bending stiffness tensor D ij depend upon circular functions of the orientations). g 1 (δ) and g 2 (δ) are the constraints imposed on the values of the extension Young moduli along the orthotropic axes of the plate, whilst E A x = 60000 MPa and E A y = 30000 MPa represent the target values of such quantities along x and y axis, respectively. The equality constraint h 1 (δ) = I (P (δ)) is actually the square of a dimensionless tensor distance (see [START_REF] Vannucci | Designing the elastic properties of laminates as an optimisation problem: a unified approach based on polar tensor invariants[END_REF][START_REF] Julien | Conception optimale de l'anisotropie dans les structures stratifiées à rigidité variable par la mèthode polaire-génétiques[END_REF] for more details). The quadratic form I (P (δ)) of Eq.( 17) is a non-dimensional, positive semi-definite function of the so-called polar parameters of the laminate [START_REF] Vannucci | Designing the elastic properties of laminates as an optimisation problem: a unified approach based on polar tensor invariants[END_REF]. Such a function depends on all the mechanical and geometrical properties of the laminate, i.e. stacking sequence, ply orientations, material and thickness of the plies. In addition, I (P (δ)) is non-convex in the space of layer orientations, since also the polar parameters of the laminate depend upon circular functions of the orientations. The quadratic form I (P (δ)) describes the requirements on the elastic symmetries that the laminate must satisfy (for the present case they are elastic uncoupling, bending as well as in-plane orthotropy with the main orthotropy axes aligned with the axes of the plate).

It can be noticed that the equality constraint on the quadratic form I (P (δ)) was transformed into an inequality one by defining a small tolerance ε I = 10 -4 for what concerns the precision of the solution in terms of elastic symmetries, i.e. I (P (δ)) -ε I ≤ 0.

Here, as a numerical example, we consider a square laminated plate with side length a = b = 100 mm, subject to pressure loads N x = N y = 1 N/mm. The plate is made of highly anisotropic unidirectional carbonepoxy layers T300/5280 [START_REF] Tsai | Introduction to composite materials[END_REF], whose material properties are shown in Table 10. We fixed the number of layers N = 16, whilst the optimisation variables are the plies orientation angles, which can assume all values between -90 • and 90 • , with a discretisation step of 1 • .

Concerning the genetic parameters that regulate the iterations of the GA BIANCA, the population size is set to N ind = 500 and the maximum number of generations is assumed equal to N gen = 120. The crossover and mutation probabilities are p cross = 0.85 and p mut = 1/N ind , respectively. Selection is performed by the roulette-wheel operator, the single-individual elitism is active and the ADP method is used for handling constraints. 10: Technical moduli for unidirectional plies of carbon-epoxy T300/5208, from [START_REF] Tsai | Introduction to composite materials[END_REF] Table 11 shows the best stacking sequence found using BIANCA for problem [START_REF] Glover | Heuristics for integer programming using surrogate constraints[END_REF] along with the values of the buckling load multiplier λ crit , the extension Young's moduli E A x and E A y as well as that of the function I (P (δ)). As conclusive remark, it can be noticed that the ADP strategy leads us to obtain very good and competitive results also when we deal with highly non-linear non-convex optimisation problem subject to non-convex equality constraints, as the one considered in this Section.

Conclusions

Traditional penalty-based approaches are very often used for handling constrained non-linear optimisation problems, due to their relative simplicity and robustness. Nevertheless, the tuning process on penalty param- eters is often delicate, as the optimal parameters themselves are problem dependent.

We have presented in this paper a new Automatic Dynamic Penalisation (ADP) constraint-handling technique, i.e. a new parameter-free method which is implemented within the GA BIANCA. Such an approach seems to be very effective when applied to the considered benchmark problems, showing very good GA search capabilities compared to the best solutions provided in literature. In addition, our approach was successfully applied to more complex real-world engineering problems showing very competitive performances.

In particular, for what concerns the engineering applications considered in this paper, i.e. the optimal design of damping properties of hybrid elastomer/composite laminates and the maximisation of the first buckling load of composite laminates with given elastic symmetries, the use of an evolutionary strategy along with the fact that the problem is stated in the most general case, leads us to find some non-conventional configurations which show better performances when compared to classical solutions available in the literature.

In addition, the ADP strategy showed to be very effective when dealing with both inequality and equality constraints, also when such constraints are expressed by non-linear, non-convex functions.

The proposed numerical tool seems to be very flexible and applicable to various engineering problems, especially when the results are given by complex and expensive models and a high number of analyses is necessary to reach a suitable optimum.

Figure 1 :

 1 Figure 1: Feasible and infeasible regions of the definition domain

1 .

 1 the individuals of the whole population are firstly classified with respect to the violation or nonviolation of the k th constraint function; 2. the individuals are then sorted into two different groups: the individuals having smaller values of the k th violated constraint are grouped in a sub-space of "virtually feasible" individuals (with respect to the k th constraint function), while the rest are grouped in the sub-space of infeasible individuals. The number of individuals grouped in the "virtually feasible" region corresponds to 10% of the population size; 3. the remaining 90% individuals in the population, which are considered "effectively infeasible" with respect to the k th constraint function are ranked in terms of their objective function values: the objective function of the best individual of such a sub-space is Φ N F best , while the k th violated constraint function is (G k ) N F best (or (H k ) N F best in the case of equality constraints); 4. the individuals within the "virtually feasible" region (with respect to the k th constraint function) are ranked with respect to their objective function values: the objective function of the best individual of such a group is Φ F best ;

  a = 1, b = 0.6, c = 0.012 and k = 0.2 are constant parameters. A 3D plot of the objective function Φ (x 1 , x 2 ) and of the constraint g (x 1 , x 2 ) is given in Fig. 2. It can be noticed that such a function is highly non-linear and non-convex with several local minima into the feasible region, while the global unconstrained minimum is placed on the infeasible side. However, the global constrained minimum lies on the boundary between the two regions and very close to the global unconstrained (infeasible) minimum.

Figure 2 :

 2 Figure 2: a) 3D plot and b) contour plot of the objective and constraint functions on the definition domain

Figure 3 :

 3 Figure 3: Distribution of the individuals over the search space along the generations for problem (6): a) initial generation, b) after 20 generations, c) after 50 generations, d) final generation

Figure 5 :

 5 Figure 5: Best values of the objective function along the generations for the optimisation problem (8)

4 10 )Figure 6 :

 4106 Figure 6: Rough sketch of the pressure vessel considered within the optimisation problem[START_REF] Davis | Genetic Algorithms and Simulated Annealing[END_REF], taken from[START_REF] Coello Coello | Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art[END_REF]

Figure 7 :

 7 Figure 7: Best values of the objective function along the generations for the optimisation problem[START_REF] Davis | Genetic Algorithms and Simulated Annealing[END_REF]

Figure 9 :

 9 Figure 9: Best values of the objective function along the generations for the optimisation problem[START_REF] Deb | Optimal design of a welded beam via genetic algorithms[END_REF]

  Concerning the genetic parameters, the population size is set to N ind = 60, while the maximum number of generations is N gen = 100. We ran the GA BIANCA with the following values of genetic operators probabilities: crossover and mutation probabilities are p cross = 0.85 and p mut = 1/N ind , respectively, while the shift operator and chromosomes number mutation probabilities are p shif t = 0.5 and (p mut ) chrom = (n max -n min )/N ind , where n max and n min are the upper and lower bounds on the number of layers, i.e. the maximum and the minimum number of chromosomes for the generic individual. Selection is performed by the roulette-wheel operator, the single-individual elitism is active and the ADP method is used for handling constraints.

Figure 11 :

 11 Figure 11: (a) Best values of the objective function and (b) number of layers along generations for problem (15).

Figure 12 :

 12 Figure 12: Evolution of species: number of individuals belonging to the same species in initial and final populations for problem (15).

Fig. 14 (

 14 Fig. 14 (a) shows the variation of the best solution along the generations: the global constrained minimum has been found after 55 generations. Fig. 14 (b) shows the first component of the homogenised stiffness tensors of the laminate, i.e. A * , B * and D * [44, 23]. We can see that both laminates are uncoupled (the curve representing B * 11 is reduced to a small black point in the center of the plot), orthotropic in extension and bending (the curves of A * 11 and D * 11 have two orthogonal axes of symmetry). Moreover, the main orthotropy axes of A * and D * are aligned with the x axis of the plate, in fact they are oriented at 0 • . Similar considerations can be done for the other components of these tensors, not shown in Fig. 14 (b) for the sake of brevity.

Figure 14 :

 14 Figure 14: a) best values of the objective function during iterations and b) first component of the homogenised stiffness tensors of the laminate, in MPa.

Table 1 :

 1 Comparison between the best solutions found with different penalty-based approaches for the optimisation problem (8) (N.A. stands for "Not Available").

	36]

Table 2 :

 2

Statistical results found with different penalty-based approaches for the optimisation problem (8) (N.A. stands for "Not Available").

Table 3 :

 3 Comparison

	Design variables BIANCA	Zhao	He and	He and	Coello and	Deb [12]	Kannan and
			et al. [49]	Wang [20]	Wang [19]	Montes [6]		Kramer [24]
	x1	0.812500	0.812500	0.812500	0.812500	0.812500	0.937500	1.125000
	x2	0.437500	0.437500	0.437500	0.437500	0.437500	0.500000	0.625000
	x3	42.096800	42.098456	42.098456	42.091266	42.097398	48.329000	58.291000
	x4	176.658000	176.636596	176.636596	176.746500	176.654050	112.679000	43.690000
	Constraints
	g1(x)	-0.000032	N.A.	N.A.	-0.000139	-0.000020	-0.004750	0.000016
	g2(x)	-0.035897	N.A.	N.A.	-0.035949	-0.035891	-0.038941	-0.068904
	g3(x)	-5.631534	N.A.	N.A.	-116.382700	-27.886075	-3652.876838	-21.220104
	g4(x)	-63.342000	N.A.	N.A.	-63.253500	-63.345953	-127.321000	-196.310000
	Objective
	Φ(x)	6059.9384	6059.7143	6059.7143	6061.0777	6059.9463	6410.3811	7198.0428
		Method		Best	Mean	Worst	Standard deviation
		BIANCA		6059.9384 6182.0022 6447.3251		122.3256
		Zhao et al. [49]	6059.7143 6059.7143 6059.7143		2.8 × 10 -12
		He and Wang [20]	6059.7143 6099.9323 6288.6770		86.2000
		He and Wang [19]	6061.0777 6147.1332 6363.8041		86.4545
		Coello and Montes [6]	6059.9463 6177.2533 6469.3220		130.9297
		Deb [12]		6410.3811	N.A.	N.A.		N.A.
		Kannan and Kramer [24] 7198.0428	N.A.	N.A.		N.A.

between the best solutions found with different penalty-based approaches for the optimisation problem (10) (N.A. stands for "Not Available").

Table 4 :

 4 Statistical results found with different penalty-based approaches for the optimisation problem (10) (N.A. stands for "Not Available").

Table 5 :

 5 Comparison

	Design variables BIANCA	Zhao	He and	He and	Coello and	Arora [1]	Belegundu [2]
			et al. [49]	Wang [20]	Wang [19]	Montes [6]
	x1	0.051613	0.051689	0.051689	0.051728	0.051989	0.053396	0.050000
	x2	0.354839	0.356717	0.356717	0.357644	0.363965	0.399180	0.315900
	x3	11.404700	11.288966	11.288966	11.244543	10.890522	9.185400	14.250000
	Constraints
	g1(x)	-0.000256	N.A.	N.A.	-0.000845	-0.000013	0.000019	-0.000014
	g2(x)	-0.000112	N.A.	N.A.	-0.000013	-0.000021	-0.000018	-0.003782
	g3(x)	-4.048164	N.A.	N.A.	-4.051300	-4.061338	-4.123832	-3.938302
	g4(x)	-0.729032	N.A.	N.A.	-0.727090	-0.722698	-0.698283	-0.756067
	Objective
	Φ(x)	0.012671	0.012665	0.012665	0.012675	0.012681	0.012730	0.012833
		Method		Best	Mean	Worst	Standard deviation
		BIANCA		0.012671 0.012681 0.012913	5.123200 × 10 -5
		Zhao et al. [49]	0.012665 0.012665 0.012665	3.200000 × 10 -7
		He and Wang [20]	0.012665 0.012707 0.012719	1.608500 × 10 -5
		He and Wang [19]	0.012675 0.012730 0.012924	5.198500 × 10 -5
		Coello and Montes [6] 0.012681 0.012742 0.012973	5.900000 × 10 -5
		Arora [1]		0.012730	N.A.	N.A.		N.A.
		Belegundu [2]		0.012833	N.A.	N.A.		N.A.

between the best solutions found with different penalty-based approaches for the optimisation problem (11) (N.A. stands for "Not Available").

Table 6 :

 6 Statistical results found with different penalty-based approaches for the optimisation problem (11) (N.A. stands for "Not Available").

Table 7 :

 7 Material properties for the glass-epoxy lamina

	's modulus E 1 [GPa] 29.9
	Young's modulus E 2 [GPa]	7.5
	Young's modulus E 3 [GPa]	7.5
	Shear modulus G 12 [GPa]	2.25
	Shear modulus G 23 [GPa]	2.25
	Shear modulus G 13 [GPa]	2.25
	Poisson's ratio ν 12	0.24
	Poisson's ratio ν 23	0.24
	Poisson's ratio ν 13	0.24
	Density ρ [kg m -3 ]	1500

Table 8 .

 8 The total number of design variables varies along with the number of plies: it can vary between 13 and 25. Even the number of DOFs of the FE model varies between 39975 and 71955.

	Design variable	Type	Lower bound Upper Bound Discretisation Step
	n	integer	4	8	-
	ID k	integer	0	1	-
	t k [mm]	discrete	0.1	1.0	0.01
	δ k [deg]	discrete	-75	90	15

Table 8 :

 8 Design variables for the optimisation problem[START_REF] Fiacco | Extensions of sumt for nonlinear programming: Equality constraints and extrapolation[END_REF].

Table 9 :

 9

Best solution found using BIANCA for the optimisation problem (15), case 3 (V denotes the position of the viscoelastic ply).

Table 11 :

 11 Best solution found using BIANCA for problem[START_REF] Glover | Heuristics for integer programming using surrogate constraints[END_REF].

Acknowledgements

The first author is grateful to the National Research Fund (FNR) in Luxembourg, for supporting this research work through Aides à la Formation Recherche Grant (PHD-09-139).