Evolutionary Design of Engineering Constructions

The basic information required to utilize one of possible computation tools/algorithms (mainly the evolution strategy) to solve a wide class of real practical engineering optimization problems is presented and discussed in the present paper. The effectiveness of the considered method is demonstrated by the possibility of the use of different form of objective functions, various and numerous nonlinear constraints and different types of design variables (continuous, discrete, real, integer). The sensitivity of the algorithm to the choice of the evolution strategy parameters is also discussed herein. The generality of the evolution strategy is illustrated by the analysis of various examples dealing with: the design of helical springs, the buckling of cylindrical composite panels and the buckling of pressure vessels with domed heads.


I. Introduction
It is well-known that the necessity of optimal design can arise in a natural way in various scientific/research fields. The search for better solutions is perhaps as old as mankind itself. Therefore the formulation of an optimization problem is well-known and well-defined in mathematical terms. However, there is still an open problem dealing with the selection of the most flexible and convenient tool for searching for the optimum which is usually called as an optimization algorithm. The algorithms for solving the real engineering problems are very complicated and in order to obtain practical solutions a significant computational effort is needed. Therefore, complex practical problems cannot be found with the use of classical methods. Over the last thirty years, different heuristic methods were formulated. H.Wilf [1] explains the sense of those algorithms in the following way: "…methods that seem to work well in practice, for reasons nobody understands…".
Optimization algorithms (methods) can be classified in different groups. One of possible classifications is given below. It combines both classical and modern approaches, including different types of heuristics.
Combinatorial methods [2] -all possible solutions (combinations -each variable takes one of a finite set value) of the problem are analyzed.
Deterministic optimization problems [3] (such as e.g. hill climbing and others) are commonly well-defined in terms of rigorous mathematical analysis. They are based on a series of axioms dealing with the linearity, convexity, unimodality (one optimum), differentiability and connectivity of the problem. Usually, the objective is given in an analytical form, and the solution is a unique one. In addition, the design space is a small one.
Mixed algorithms (se for instance Refs [22][23][24][25][26][27][28][29][30] and especially chaotic neural networks [31][32][33][34][35][36]) combine the best features of the deterministic and stochastic analysis. The optimal solutions are independent on an initial information. They are effective in the sense of possibility of finding global unimodal or multi-modal optimal solutions for extremely hard, often poorly specified problems. The analytical description of the problem is not required in order to find a global optimum -they may work as a black box. On the other hand, the computational effort is reduced to the minimum.
In the late 60s several researchers focused their attention and invented a new type of optimization methods that are called now as evolutionary algorithms (EAs). Now, those techniques are also known under the name metaheuristics [37]. Therefore, the proposed herein method of the classification is not a unique, and, of course, is one of possible proposals. However, it demonstrates evidently the present tendency of the optimization strategy in the sense of the search for the convenient mixed algorithms having the highest effectiveness and robustness in the multimodal optimization problems with objective functions evaluated in a numerical way. Such a group of algorithms is also suitable for solving discrete optimization problems arising with the use of design codes constituting a discrete space of design variables. It is well-known that no optimization method is better to any other. A broader discussion of those problems is presented e.g. in Ref [38]. However, not trying to propose better methods is an erroneous conclusion.
Taking the above remarks into consideration we intend to build an optimization algorithm that can solve a broad class of engineering problems encountered in practice, possessing a large number of design parameters (integer, discrete, continuous and integer/discrete-continuous) and a large number of constraints. Genetic algorithms (GA) has been highly successful as one of evolutionary computation techniques in searching for a broad class of stacking sequence, size, topology optimization problems for composite structures (see Refs [39][40][41]). However, the application of that optimization technique is limited to a very few inequality constraints only. Now, our attention will be particularly focused on the formulation of the alternative to GA evolutionary computation optimization technique. It is based on the idea of the evolution strategy so that the proposed algorithm allows us to eliminate the problems associated with taking into account a large number of constraints.
In details, the main objective of the present considerations will be two-fold: 1. To present the classification of evolutionary algorithms and the actual terminology dealing with the application of evolutionary algorithms in the structural optimization problems in order to emphasize their differences and similarities and on the other hand to eliminate ambiguities and misunderstandings arising with the use of the improper term "evolutionary optimization" -see e.g. Ref [42]. In fact, the term evolutionary refers to one individual part of the optimization process and directly to an optimization algorithm only what is clearly underlined in Ref [43] 2. To illustrate the generality, effectiveness and accuracy of the described herein, completely new evolution strategy that is based on the appropriate analysis of constraints (feasible regions).

II. Evolutionary Computation
Now, in all branches of sciences an evolutionary computation plays a very special and significant role among various optimization algorithms. By many researchers the evolutionary computation area is divided in four specific domains: (i) genetic algorithms (GAs), (ii) evolutionary programming (EP), (iii) evolution strategies (ESs), (iv) genetic programming (GP).
In fact, now, hybrid of the above-mentioned four algorithms (metaheuristics) are employed in the analysis. Implementing evolutionary algorithms a series of operations should be taken into considerations. They are summarized in Table 1. However, it is necessary to point out that all operations that are conducted in four groups of evolutionary computations have their own meaning and are different. For instance, it is impossible to give the identical meaning of the reproduction or the crossover or the recombination operations. In genetic algorithms reproduction means the selection of individuals from the initial populations using only the fitness value. Thus, for instance having the size of the initial population equal to 10 and having randomly generated 10 different individuals, after reproduction 5 different individuals only may create a population having the same size as previously, i.e. equal to 10. Then, the crossover operation may be conducted. As it may be seen those operations are not connected with the form of chromosomes. In addition, it is worth to emphasize that according to the standard [43] recombination has a different meaning than crossover since crossover is limited to so-called sexual recombination only. In general, the use of different terms in the brief presentation of evolutionary algorithms should demonstrate different operations and the difference in specific groups of evolutionary computations.
The foundation of genetic algorithms (in 1950s) is connected with the name of the biologist A.S. Fraser [44]. J. H. Holland [45] (1962) with the students (K.A. De Jong, D. E. Goldberg) have had the greatest influence on the developments of GAs.
In 1966 L. J. Fogel et al. [46] presented a new idea called now as evolutionary programming. The lack of a crossover operation is a characteristic feature of this method.
Working on configurations of hinged panels in a aerodynamic tunnel I. Rechenberg [47] and H.P. Schwefel [48] invented a new optimization method known now as the evolution strategy.
Latin American Journal of Solids and Structures, 2018, 15(4), e87 3/21 The last major method of EC was founded by R. M. Friedberg [49] (called as genetic programming) and it is joined directly and particularly with the work on the length of computer programs. As it may be noticed GA, EP, ES and GP can be characterized by a variety of qualities in common, although they also differ in a sequence of operations and the termination of operations. Differences between all four of them seem to become more and more indistinct.

III. Optimization Problem Formulation
An arbitrary engineering problem dealing with the optimal design has three fundamental common features, such as: 1) selection one or more design variables, 2) choosing an objective function, 3) identification a set of equality or inequality constraints. In our problem design variables can be described as the vector x  having I components xi (i=1,2,…,I). It is also assumed that each component, dependent on the problem considered, can belong to a set of natural, discrete or real numbers. A discrete set of numbers corresponds directly to constraints existing in various codes where dimensions of real products are not continuous but discrete real numbers. All components of the vector x  must be positive as variables corresponding to real dimensions of constructions.
They are constrained by lower and upper limits -Eq. (1), i.e.: To simplify numerical computations they are normalized in the following manner -Eq. (2): The constraints are given in the inequality or equality form -Eq. (3): Now, one objective function is analyzed only and it is always formulated as searching for the maximum of the objective function F -Eq. (4), i.e.: In a classical manner all Min or MinMax optimization problems can be transformed to the above.

IV. Description of the Evolution Strategy
According to Rechenberg's [50] conclusions the field of evolution strategies can be characterized as the evolution of evolution. Usually, the evolution strategies utilize three types of operations: 1) recombination -Eqs (5), (6), 2) mutation -Eq. (7) and 3) selection.
During the recombination operation (see Fig.1 R denotes a function and vel, migr are called as strategy parameters. All new children (particles in Fig.1) should satisfy all equality or inequality constraints (3). In the present formulation the restoration method is used -see Fig.2. Other methods are described e.g. in Ref [51]. In general, the restoration method is based on returning back new produced children to the boundary by a simple operation that can be written in the following way:

Figure 2. The restoration method
In the evolutionary computation methods there is a problem with handling of any type of constraints. Some of the authors (e.g. Ref [52]) proposed that in the optimal analysis infeasible solutions are not allowed.
Latin American Journal of Solids and Structures, 2018, 15(4), e87 5/21 The optimization process is terminated as for the assumed objective function F the following condition is satisfied: for all members of the population. The broader description of all assumptions and procedures introduced in the proposed version of the evolution strategy can be found e.g. in Ref [53].

V. Applications of the Proposed Algorithm -Numerical Examples
In this section we intend to show how the proposed algorithm works in some engineering problems. They are limited to the analysis of isotropic and 2D composite structures and mainly devoted to the cost optimization problems understood in the sense of the volume optimization (minimization) problems. In general, we intend to demonstrate how the algorithm works as the various types of inequality constraints exists. They number is even higher than the number of design variables but such a problems are typically encountered in engineering practice. It should be mentioned also that the field of possible applications is (and can be) much higher.

A. Design of Helical Springs
In various engineering applications helical compression springs constitute an important element. It is wellknown that they not only can exert external compressive force but also provides the structural flexibility and can store or absorb energy. In the machine design such elements should satisfy a lot of various mechanical requirements. For helical springs the optimization problem can be formulated in the following manner: to minimize the total volume (weight) of the spring. The volume of the spring is the sum of the volume of active and inactive coils.
where the symbols N, Q denotes the number of active coils, and the number of inactive coils, respectively. Q is a constant and it is assumed to be equal to 2. D is the mean coil diameter, whereas d is the wire diameter (see Fig.  3). In order to use the consistent notation three variables (N, D, d) in Eq. (8) will be treated as the design variables and are denoted by xi (i=1,2,3), respectively. There are seven constraints for this problem -Eqs (9)-(16): maximum allowable stresses: allowable length: lower limit on wire diameter: maximum allowable outside diameter: allowable spring index: pocket length: clash allowance: where Kw is the Wahl correction factor, Fmax=4. 45   The bounds for the analyzed problem have been chosen in the following way:   The above optimization problem have been introduced and analyzed by Sandgren [54]. As it may be noticed the present results obtained with the use of the introduced evolution strategy is identical to that demonstrated in Ref [54]. -see Table 2. However, it should be emphasized that the analysis deals with continuous, real numbers (design variables). In fact, in engineering applications the number of active coils is an integer number, whereas the wire diameters d are described by a discrete set of real numbers corresponding to those determined by design codes. The discrete set of allowable wire diameters is shown below (in inches). Comparing the above values and the constraint (11) one may notice that 11 values of the wire diameters can be considered in the optimization problem. Fig. 4 demonstrates the convergence of the optimization algorithm.
Beginning from the 40 th generation there are individuals characterizing the optimal solution presented in Table 2. However, all individuals have the identical optimal values at 125 th generation and the accepted error between them (7) is equal to 10 -7 -see Fig. 4 d.
The solution obtained with the use of the presented method for the new set of design variables (N -an integer number, D -a real continuous number, d -a real discrete number) is demonstrated in Table 3 and compared with the data available in the literature. Similarly as previously the optimization algorithm allows us to obtain very good results. However, the accepted error is reached later than previously (Fig.4) after 230 th generations assuming the identical, as for continuous design variables, values of parameters characterizing the evolution strategy.  Although, as it is proved above, the proposed evolution strategy allows us to obtain very good, even excellent numerical results it is necessary to mention about the influence and the role of strategy parameters, i.e. the total number of individuals in the population (pop), the velocity of the strategy (vel) and the number of migrations (migr). They can affect significantly the convergence as well as the efficiency of the optimization algorithm since they have a great influence on the so-called oscillations of solutions -visible also in Fig.5. When vel is very high it hits the boundary on nearly every iteration and this ineffective trajectory searches the same points repeatedly. As vel is low the particle explores the optimum very slowly -see  2  2  33  23  11 23  12 13  13  12 23  22 13  11 22  12   2  2  3  2  2  2  11  11  66  12  12  66  13  12  11  12  66  22  22  66   2  3  23  12  66  22 Assuming that the term s in Eq. (18) is identically equal to zero the relation (18) what means that the second term (membrane) in Eq. (3), denoted by the symbol s, plays a significant role. It is seen that the critical load is a function of geometrical ratios t/R and L/R. For isotropic shells to find the optimal buckling force let us assume that the derivative of Pxcr with respect to L is equal to zero and the corresponding minimum buckling load -Eq. (20) is: For composites the relation between the bending and membrane states depends not only on the values of the geometrical ratios t/R and L/R, but also on the mechanical properties and laminate configurations. Therefore, the identical procedure to that mentioned for isotropic structures cannot be conducted. It is seen in Min P (21) subjected to six constraints: Using such a representation of the design variables it is possible to write the membrane part in Eq. (18) as the polynomial function of two variables A A x x 3 1 , . The contour plots demonstrate that the location of the maximum is dependent on the geometrical ratios so that it is impossible to start the optimization procedure from the derivation of the optimal shell length.
where Z1, Z2 and Z3 are constants in the optimization problem. Since the constants are positive the maxima always occurs at the boundaries ( ). It is worth to emphasize that the contour plot solution gives directly the number of plies having orientations 0 0 , ±45 0 , 90 0 but not their location in the laminate. The location of those plies can be found by the analysis of the values D D x x 3 1 , . However, having an information about the positions of the optimal stacking sequences for bending states it is possible to find independently the positions of maximal buckling loads (the optimal stacking sequence) for membrane states only, and then to join the results. With the use of the decoding procedure is possible to find a series of the corresponding to values

C. Buckling of Cylindrical Shells with Domed Heads
In engineering constructions composite cylindrical shells are closed by domed heads having various forms as: flat plates, ellipsoidal, torispherical, spherical, paraboloidal or even conical. We intend to analyse the influence heads on buckling loads and forms. The presented analysis is a parametric study of various geometrical effects and it is treated as an introduction to optimal design in the form similar to that studied by Muc [53]. In the present optimal design we tend to minimize the geometrical ratio L/D where the value of the buckling external pressure (pbuckl) does not exceed the allowable value of the pressure (pallow) -Eq. (24), i.e.: buckl allow p p  (24) In particular the considerations deal with composite pressure vessels having the domed heads in the form of: a) shallow torisphere (R/D=1, r/D=0.1) with a cylindrical part - Figure 7a, b) deep paraboloid (f/D=1.5, i.e. z=6r 2 /D) with a cylindrical part - Figure 7b. It is assumed that the shells are made of carbon/epoxy resin and have quasi-isotropic material properties. The numerical investigations are carried out for axi-symmetric structures taking into account geometrical pre-buckling nonlinearities. The results of buckling investigations for spherical and ellipsoidal domes were studied by Muc [54,55].
For identical thicknesses of domed heads and cylindrical parts the buckling loads decrease with the increasing length of the cylindrical part - Figure 8. If the thickness of the head is lower than in the cylindrical part one can observe the region as buckling pressures arenot sensitive to the change of the geometrical L/D ratio. In the latter case the loss of stability in the head preceeds the buckling of the whole structure. The effect of the geometry is much more manifestant for deep paraboloidal shells than for shallow torispheres. The results of buckling investigations for spherical and ellipsoidal domes with and without cylindrical portions and the possibility of optimal design were studied in details by Muc [57,58]. For a infinite laminated plates subjected to axial/biaxial compression/tension or in-plane shear and having a circular/ elliptical hole the optimum design of laminate stacking sequence is presented by Sharma et al. [59]. The authors of the paper adapted the formulation and the solution given by Lekhnitski [60], Savin [61] for stresses around holes in an infinite angle-ply balanced symmetric laminated plates under in-plane loading (Fig.1) Tij represent the classical terms of the compliance matrix for the plane stress problem.
where α is the angular coordinate measured along the ellipse from the axis 0-x, (σy=σ∞ =N0/t at the direction y - Fig.11), and s1, s2 are the complex roots of the equation (25). For orthotropic or cross-ply structures the roots can be expressed by the relations (27): where: 11 and Aij denote the in plane stiffness of the laminate.
In an infinite orthotropic plate the analytical solution for the stress component σy at the edge of a hole and α=0 0 is given by Eq. (28): However, the above relations cannot be used in the optimization analysis since for finite width of the plate the stress distributions/concentrations around the holes are different than those compared with the results for infinite plates. It was demonstrated by Muc, Ulatowska [62] and Xu Xiwu et al. [63]. Tan [64] introduced the approximate formulas for the stress concentration factors for the orthotropic laminated plates valid in the range 0≤b/a≤1 -Eq. (29) where Kt∞ is the stress concentration factor for the infinite plate that can be derived from the following formula and Aij denotes the effective laminate in-plane stiffnesses with 1 and 2 corresponding to parallel and perpendicular to the loading directions (Fig.11), i.e. to the y and x directions, respectively. The optimization problem is formulated in the following way: to maximize the tensile load -Eq. (31), i.e.: subjected to the constraints corresponding to each of the individual plies (i) (i=1, 2,…, N -total number of plies) in the orthotropic plate -Eqs (32):  2   1  1  1  ,  ,  ,  1  1  1 1  1  ,  ,  ,  2   xx  yy  ss  t  c  t c  x  y  xy  xx yy  t  c The constraints represent the First-Ply-Failure (FPF) in the form proposed by Tsai-Wu. The criterion is formulated in the local coordinate system (1-2) defined previously. Xt, Xc, Yt, Yc are allowable stresses for tension (t), compression (c) in the directions parallel (X) and perpendicular (Y) to fibres, and S denotes the allowable inplane shear stresses. For each of the plies the transformation rule between the local (1-2) and the global coordinate system is defined by Eq. (33) below. The transformation is valid for the (ith) ply, however the symbol (i) is omitted.
Similarly as in the previous section for the assumed discrete fibre orientations 0 0 , ±45 0 , 90 0 in the laminate let us use the design variables in the form of two integers    It was found that the increase of the curvature resulted in the increase of the maximal failure strength, i.e. the highest failure strength was observed for the horizontal ellipses and the lowest for the vertical elliptical holes. The obtained optimal fibre orientations presented by Sharma et al [59] are different than those presented above since the mentioned authors analysed infinite plates only.
For discrete fibre orientations the distributions of the strength were analysed in the triangular space of design variables. The results are plotted in Fig. 12 for circular hole. The maximum of the strength occurred for fibres oriented at 0 (8 plies) and 90 (8 plies). The plies oriented at 45 were not present in the optimal stacking sequence. The sequence of the plies oriented at 0 and 90 is not important in the analysis since the values of the terms Aij do not vary with the change of the stacking sequences.
For laminated multilayered plates with a cut-out more information about the stress concentration problems can be found in Refs [65][66][67].

E. Vibration Frequency of the Cantilever Beam-Mass System
Let us consider the cantilever beam system shown in Fig.13. The beam cross section is rectangular. The objective function is to find the optimal cross-sectional dimensions a and b that minimize the weight of the beam (the beam length is constant) while keeping the fundamental vibration frequency larger than 8 rad/sec. Thus, the fundamental optimization problem can be represented by Eqs (34) Figure 13. The cantilever beam-mass system The analytical solution of that problems exists and the optimum can be found at the boundary:   The use of the proposed optimization algorithm allows us to obtain the identical results as above very quickly, after 50 generations with the error ε equals to 10 -7 . For different number of generations numerical results are plotted in Fig.14 -dimensionless design variables. As it may be seen they present an excellent convergence and accuracy.

F. Optimal Proportion of the I-shaped Cross Section
A lot of modern sheet isotropic or composite constructions require stiffeners that are normally of I-, Z-or Cshaped cross sections. An I-shaped cross section, shown in Fig. 15, has been selected in our example. Since a large number of these stiffeners are employed, it is important to find the optimal proportion of its dimensions x1, x2 and x3 being design variables in our problem. A stiffener is treated as an axially loaded column, as shown in Fig. 15. The goal is to have a minimal volume of the the column (understood in the sense of the area A since the column length L is constant): while meeting the following yield stress and buckling conditions: yield stress: overall buckling load:  (41)   3  2  3  3  3  1  1  1  1  1  2 1  1 2  3  1  1  3  1  1 2  1  3  1  6  2  12  6 12 x y flange buckling load: web buckling load:  Figure 15. The axially loaded column having the I-shaped cross-section The problem considered (Eq. (39)) is much more complicated than the previous one since now we have 3 design variables, 3 bounds on design variables (Eq. (45)) and 4 inequality constraints -Eqs (40)- (44). It cannot be solved in an analytical way. Using the proposed optimization algorithm we can find the optimal dimensions of the I-shaped cross-section. Their values are presented in Fig.16 and Table 4 as a function of the column length L. For short columns the optimum occurs at the bounds of the cross-section thickness and width (compare Table 2 with Eq. (45)) and the growth of the optimal area is observed with the increase of the column length (Fig.16). Further growth of the column length L results in the reduction of the optimal cross section area A (comparing with the previous values obtained for the lower values of the column length) and in addition all three variables describing the I shaped cross section vary. For all cases considered herein the optimum is observed at the boundary characterized by Eq. (41) -the overall buckling.

VI. Concluding Remarks
The present review of existing works in the open literature demonstrates evidently that the number of results rapidly increases without the proper comparison of them. Therefore, it is practically impossible to draw more general conclusions even in the specific area of interest. Although, the number of references in each of the papers also grows it is very difficult to track and compare the validity of the published results in view of various assumptions made in different works. Such a situation may lead to misunderstanding and pitfalls for researchers.
We do not intend to discuss herein the problem of the choice of the design variables (see Ref[68].). However, the investigations are mainly focused on the assumptions that may affect significantly on the results and then on the comparison of results.
The detailed analysis demonstrated evidently that the proposed algorithm works well for different optimization problems, various number of design variables and equality or inequality constraints. More information about those problems can be also found in Ref [53].
For the analyzed optimization problems the effectiveness of the use of the evolution strategy is demonstrated since after few iterations particles (individuals) in the population cluster in one (sometimes several) region of the search (design) space. These clusters indicate the presence of the optima where idividuals' relatively good performances have caused them to attract their neighbors, who in moving toward the optimal regions improved their own performances, attracting their neighbors and so on. The efficiency of the algorithm can be evaluated by the mean position of a clustered population and the relative distance from an individual particle position to the mean value. As it is demonstrated the use of the proposed algorithm allows us to find the position of global optima with the required accuracy for different engineering problems having a lot of nonlinear constraints in the equality or inequality form. The generality of engineering problems (understood in the sense of the form of the objective function as well as of the form of constraints) that can be solved with the use of the proposed methods seems to be their fundamental advantage. According to the authors' knowledge and experience the algorithm works much better than genetic algorithms. In addition the method can be applicable not only for isotropic but also for multilayered composite structures. The effectiveness of the algorithm is strongly dependent on the proper choice of the evolution strategy parameters.