Abstract
Both the entire weight and joint displacements of grid structures are minimized at the same time in this study. Four multiobjective optimization algorithms, NSGAII, SPEAII, PESAII, and AbYSS are employed to perform computational procedures related to optimization processes. The design constraints related to serviceability and ultimate strength of grid structure are implemented from Load and Resistance Factor DesignAmerican Institute of Steel Constructions (LRFDAISC Ver.13). Hence, while the computational performances of these four optimization algorithms are compared using different combinations of optimizerrelated parameters, the various strengths of grid members are also evaluated. For this purpose, multiobjective optimization algorithms (MOAs) employed are applied to the design optimization of three application examples and achieved to generate various optimal designations using different combinations of optimizerrelated parameters. According to assessment of these optimal designations considering various quality indicators, IGD, HV, and spread, AbYSSS shows a better performance comparatively to the other three proposed MOAs, NSGAII, SPEAII, and PESAII.
1. Introduction
The grillage systems utilized in different structures like bridge or ship decks, building floors and space buildings, and so forth, contain traverse and longitudinal beams, which are made of available steel profiles with different crosssections. The optimal selection of steel crosssections from a discrete set of practically available steel profiles provides a big contribution to constructing cost of a grid structure. Therefore, either weight of grid structure or deflection of its joints is minimized according to certain design limitations prescribed by any code of practice, such as LRFD. During the design optimization of grillage systems, designer is frequently faced with a problem related to making a decision about determination of the most appropriate one between these two conflicting and commensurable objective functions. Although a displacementrelated constraint is imposed as a (max span/300) according to the provisions of LRFDAISC specification, the safety margin on displacement constraint is large when taking into account the grid structures with higher sensitivity against displacement, such as ship decks and floors of industrial buildings which bears special machines required an horizontally balanced position for a regular work. This task has been easily overcome in a way of introducing the concept of multiobjective optimization to the design applications of grid systems.
Preliminary multiobjective optimization techniques, which their fundamentals were constituted on mathematical programming principles, were dated back to the 1950s. However, mathematical programming techniques adjust the decision variables of continuous type according to gradient information computed by use of objective functions. Furthermore, they fail when search space is concave and discontinuous. In order to deal with these tasks, alternative optimization procedures, which mimic various natural events, for example, evolutionary systems, immune systems, social behaviors of ants, insects, and animals, have been developed. The most preferable ones are the evolutionarybased algorithms. These biological evolutionary models utilize characteristic properties of nature, for example, heredity, selection, and so forth, to create populations with higher qualities. Particularly, genetic algorithms (GAs) are the most flexible multiobjective evolutionary tools due to allowing an implementation of various operators for its evolutionary computation. Therefore, it has been hybridized with different local search techniques.
First study on the multiobjective evolutionary optimization, called Vector Evaluated Genetic Algorithm (VEGA) was performed by use of GAโs principle [1]. VEGA employs a number of subpopulations to search the solution space considering a modified selection mechanism. Following the emergence of VEGA, two featured approaches, linear aggregating and lexicographic ordering of objective functions were developed [2, 3]. They transformed all objective functions into a single objective function by way of optimizing each objective function without decreasing their solution qualities.
As an alternative to early attempts mentioned above, a paretobased evolutionary approach was developed to increase the population diversity [4]. Some of the preliminary paretobased multiobjective optimization approaches are Nondominated Sorting Genetic Algorithm (NSGA) by Srinivas and Deb [5], a Niched Pareto Genetic Algorithm (NPGA) by Horn et al. [6], a Multiobjective Genetic Algorithm (MOGA) by Fonseca and Fleming [7], and a Multiobjective Evolutionary Algorithm (MOEA) by Tanaka and Tanino [8]. In order to improve their optimal results, these algorithms have been developed by either enhancing their current optimization strategies, such as Nondominated Sorting Genetic Algorithm II (NSGA II), Improved Strength Pareto Evolutionary Algorithm II (SPEA II), Improved Pareto EnvelopeBased Selection Algorithm (RegionBased Selection) II (PESA II), or adapting a competitive Search Technique, such as Adapting Scatter Search AbYSS.
Although it is clear that there are a number of evolutionary methods known in the field of evolutionary multiobjective optimization, in this paper, an exhaustive literature review on this field is omitted. Instead, representative works related to four multiobjective optimization techniques employed are summarized. Furthermore, the recent multiobjective approaches utilized in the field of structural engineering design are also reviewed.
In this regard, paper is organized as firstly introducing the first steps in MOAs including applications in the field of structural engineering after a brief introduction to the multiobjective optimization problem and concepts. The computational procedures of proposed MOAs, NSGAII, SPEA II, PESA II, and AbYSS are presented in Section 3. The design requirements prescribed by LRFDAISC Ver.13 and optimal design procedure are given in Section 4 prior to the introduction of search methodology located in Section 5. Following the discussion of results along presented in Section 6, final remarks are summarized in Section 7.
2. Background
2.1. Multiobjective Optimization: Problem and Concepts
A general multiobjective optimization problem consisted of objective functions and constraints, defined by decision variables, is represented as follows:
bounded by an upper and lower value, and are used to represent a decision variable set defined in decision variable space (DS) and computation of both objective functions and constrains and in a solution space (SS). In order to explore the optimal solutions (designations) located in feasible region (FR) which contains the unpenalized solutions of SS, the obtained solutions are penalized when they do not satisfy the constraint conditions. Then, the penalized values denoted by are included into their related objective functions .
At each run of an evolutionary optimization algorithm, a random solutions set is obtained. Some of them are Nondominated solutions (none is better for all objectives) and referred to as โpareto solutionโ defined in a concept named as domination [9]. Thus, the pareto solutions are used to form โpareto frontโ which determines bounds of Nondominated solutions.
2.2. First Steps in MOAs
After introduction of the evolutionary mechanism to the multiobjective optimization problem, the development of new multiobjective approaches has been accelerated [4]. Primary one of these approaches is Nondominated Sorting Genetic Algorithm (NSGA) proposed by Srinivas and Deb [5]. NSGA determines the Nondominated solutions according to ranks of their reproductive potentials.
Although paretoranking procedure gives a guarantee for transmission of elite individuals to next generations, an excessive repetition of ranking procedure causes to lose promising genetic material due to the genetically distortion of migrated data. In order to diminish this negative effect, Niched Pareto Genetic Algorithm (NPGA) that determined the Nondominated solutions by tournament selection method was suggested [6]. Afterwards, Fonseca and Fleming [7] suggested a penalization process that generated the promising pareto solutions considering their crowding densities.
Using an external archive to store Nondominated solutions leads to an increase in the capabilities of multiobjective optimization tools, such as in an evolutionary search proposed by Zitzler and Thiele, called Strength Pareto Evolutionary Algorithm (SPEA) [10]. The members of external archive are chosen according to their closeness to pareto front. However, enlargement of the external archive makes the convergence speed of its evolutionary search to be poor due to a decrease in its selection pressure. In order to deal with this difficulty, Knowles and Corne [11] suggested a grid system, called Pareto Archived Evolutionary Strategy (PAES) to compute the optimizationrelated procedures. Hence, distributing the entire population to this grid system by an adaptive mapping process, which of each node was used to represent an individual, makes it easier to maintain the diversity in pareto set.
2.3. An Overview of MOAs Applied in Field of Structural Engineering
Although MOAs mentioned above are successful for a solution space represented by design variable of continuous type, they fail to explore optimal designations in nonconcave and discontinuous solution space of structural design problems represented by design variables of discrete type. Hence, the preliminary studies in the field of structural engineering with multiple objectives were developed using weighting [12], goal programming [13, 14], and modified game theory methods [15]. Sunar and Kahraman [16] compared the computational performances of these algorithms considering optimal designations of a space truss with 25 bars and a satellite system and showed that modified game theory and goal programming were superior to weighting approach. Although it was reported that the weighting approach failed to explore Nondominated solutions on nonconvex parts of pareto front [17โ19], it has been improved due to its easy adaptable mechanism. One of these attempts was based on a systematic alteration of objective function weights [20]. This adaptive approach was applied to design optimization of a truss problem with 3 bars and achieved to obtain a welldistributed pareto set in a nonconcave solution space.
Realworld engineering structures are represented by a large number of discrete design variables. Hence, the evolutionary search is easily misguided due to an increase in the computational effort. Therefore, new optimization algorithms inspired by some biological events, such as particle swarm, differential evolution [21], artificial immune system [22], ant colony [23], and some other evolutionary algorithms, such as microgenetic algorithm [24] have been utilized to solve optimization problems with multiple objectives.
Janga and Nagesh proposed an evolutionary technique, called elitistmutated particle swarm optimization and applied to three test problems: design optimizations of twobar truss, Ibeam, and welded beam [25].
Cooperate and coevolutionary strategies were introduced to the evolutionary search mechanism to increase diversity within the set of Nondominated solutions stored [26]. For this purpose, chromosomes divided into different species are recombined and evolved to create a pareto set. However, the number of chromosomes chosen may exceed their predetermined number. In order to deal with this difficulty, the exceeded number of chromosomes is reduced according to the qualities of their crowding distances and redistributed using an elitism strategy. The fundamentals of NSGA [5], Niched Pareto Genetic Algorithm (NPGA) [27, 28], and Controlled Elitist Nondominated Genetic Algorithm (CNSGA) [5] were constituted on this approach. The proposed evolutionary approach, which of main evolutionary operations has similarities to the optimization algorithms named NSGA, NPGA, and CNSGA, were utilized for topology optimization of 2D heat transfer structures. It was demonstrated that using a decreased size of species increased the computational performance of the proposed optimization procedure.
The other promising approaches were developed as hybridizing these methods with each other. One of substantial attempts is the integration of neural network and fuzzy systems to provide a control mechanism for an evolutionary mechanism [29]. While a neural network is utilized to predict an individual with higher quality, operator parameters of evolutionary algorithm are updated according to the rules of fuzzy logic. This hybrid system was applied for design optimization of several composite beams with three layers, piezoelectric bimorph beam, a truss structure, and airplane wing and displayed to obtain more converged optimal designations comparing to a pure and independent usage of each algorithms without any hybrid implementation.
3. Introduction of MOAs, NSGAII, SPEA II, PESA II, and AbYSS for Design Optimization of Grillage Systems Utilizing Multiple Objectives
Although it is displayed that the evolutionary algorithms have been successfully utilized as optimization tools, some computational bottlenecks cause to obtain a poordistributed pareto front. In order to increase their ability in generation of promising designations for a diverse pareto front, the conventional MOAs have been improved and/or developed. The promising versions are NSGA II, SPEA II, PESA II, and AbYSS. In this study, these improved or developed MOAs, are utilized for the optimal designations of grillage systems noting that the design optimization of grillage systems was carried out using a single objective by up to now [30โ32] and a displacementbased matrix analysis approach [33]. In order to execute their optimization procedures, JMETAL coded in a programming language Java is employed [34]. Both JMETAL and the proposed four MOAs are introduced in the following subsections.
A Brief Introduction of JMETAL
JMETAL coded in a platform of objectoriented Java is used as an optimization tool to solve the multiobjective optimization problems [34]. It contains a number of classes which represent the building blocks of various multiobjective algorithms. However, their basic evolutionaryrelated elements are the same. Therefore, the architecture of JMETAL is constituted on a simple but an interdependently higher framework. JMETAL is an open source project. The computational procedures of various MOAs including the ones employed in this study are extensively documented in [34]. Therefore, the description of its base classes is not presented. Instead, the computational order of the preliminary classes that contain the fundamental parameters is briefly depicted by the pseudocodes. Whereas the class names are presented by use of a mark โ โ, their related parameters are defined by use of italic characters.
In general, JMETAL contains six packages, named โbase, experiments, metaheuristics, problems, quality indicator and utilโ. While the package of โexperimentsโ includes configurations of various metaheuristic procedures, such as AbYSS and NSGAII, their default configurations are located in the โsettingsโ package. Design problems with various complexities are comprised in the โproblemsโ package. The run of JMETAL firstly begins by execution of a purposed metaheuristic class, such as AbYSS.Java, located in the โexperimentsโ package. This class invokes the experiment class which is responsible for activation of proposed metaheuristic algorithm. The purposed metaheuristic class invokes a class of design problem to be solved. Basic computation processes, named โevaluate () and evaluateconstaint (),โ are carried out in this class, named design problem. Due to the use of binary coding scheme for the evolutionary computations, generated individuals can exceed the predefined upper bounds of design variables. In order to deal with inappropriate binary strings, the purposed metaheuristic class named AbYSS.Java, and its related subclasses are extended by including a new class named โrepairing mechanismโ which is responsible to correct the decoded values of binary strings according to the limits of design variables.
3.1. Nondominated Sorting Genetic Algorithm II, NSGA II
NSGA was firstly designed by Srinivas and Deb [5]. The evolutionary computation of NSGA is managed by classified โpopulation sizeโ individuals. Classification process begins firstly ranking the population in order to determine Nondominated individuals. Then, fitness values of individuals are shared considering their niching measures. The individuals with shared fitness values are selected by use of a selection method, called โstochastic universal samplingโ and regenerated using mutation operator with a probability mutation_probablity and combination operator with a probability crossover_probablity until completing a predetermined evolution number max_evolution. The enhanced version of NSGA, NSGA II (see pseudocode in Algorithm 1) was developed to increase diversity among individuals [35, 36]. Although NSGA II utilizes the evolutionary operators for generation of individuals as in the computational procedure of NSGA, the generated individuals stored in an archive with population_size individuals are used to obtain a pareto front considering crowding distances of individuals. Furthermore, a hypercube, which is formed by use of the Nondominated individuals, is utilized to compute a hypervolume.

3.2. Improved Strength Pareto Evolutionary Algorithm, SPEA II
The primary version of SPEA approach was firstly suggested by Zitzler and Thiele [10]. Its evolutionary processes are managed by two populations. While one of these populations, called regular population is utilized to generate offspring, other population, called an archive with archive_size individuals is employed to preserve the evolutionary information of pareto front. In the beginning of evolutionary process, the archive is empty and filled by promising individuals. The exceed number of individuals are reduced giving a higher chance to the individuals of archive. For this purpose, a clustering technique which is based on an assignment of strength value to the individuals and an assessment of these individuals according to their strength values is utilized to discard the related individuals. However, a decrease in variation among individuals of population causes to increase the randomness in the evolutionary search and hence decrease the accuracy degree of density estimation, which is utilized by clustering process [37]. This leads to vanishing of the promising solutions located on the pareto front. In order to deal with this negativity, SPEA II is developed. SPEA II estimates density of strength values using th nearest neighbor method (see Algorithm 2). Furthermore, a new selection mechanism, called โenvironmental selectionโ, is used to update the archive and to preserve promising solutions during truncation of archive.

3.3. Improved Pareto EnvelopeBased Selection Algorithm (RegionBased Selection), PESA II
Pareto EnvelopeBased Selection Algorithm (PESA) was firstly presented by Corne et al. [38]. Having similarities to basic features of SPEA II, PESA uses two populations for its evolutionary processes. However, in the estimation of strengthvalue density, PESA uses a measure called โsequence factor.โ An extended version of PESA, PESA II, utilizes hyperboxes which are obtained by dividing the entire search space into small ones (see Algorithm 3). Thus, the number of individuals contained in hyperboxes is used to determine a sequence factor. Following the determination of subregion numbers defined by bisections, the archive is created in the module named โAdaptiveGridArchive.โ A population called โsolutionโ is initialized and evolved by the application of evolutionary operators until being completed a fixed evolution number. Although SPEA II was shown to be computationally faster than NSGA II and SPEA II, some complicating tasks, such as possibility of existing dominated individuals in any hyperboxes, keeping a fixed number of subregions throughout the evolutionary search, and so forth, must be overcome in order to obtain satisfactory optimal designations [39, 40].

3.4. Adaptive Scatter Search for Multiobjective Optimization, AbYSS
AbYSS can be categorized as an evolutionarysearchbased optimization algorithm derived by use of principle features of NSGA II, PESA II, and SEPA II [41]. A pseudocode for ABYSS is presented in Algorithm 4. After initializing a population with population_size individuals, the population is firstly regenerated using gridbased search technique named โdiversification generation.โ Then, the regenerated population is maintained by discarding the mutated individuals with poor qualities according to a dominancebased comparison test. This elimination process named โimprovementโ is inspired from NSGA II approach. Then, a search process named โreference updateโ is invoked to firstly construct a subpopulation RefSet2 with Ref_Set_size individuals from a subpopulation RefSet1 with Ref_Set2_size individuals and then, update the subpopulation RefSet2. The subpopulation RefSet2 is built by the individuals of RefSet1 with minimum euclidian distance. After generation of reference sets, the Nondominated individuals extracted from the reference sets are stored in an external population called archive with archive_size individuals using density estimation of individuals. In this regard, the existence of individuals in a densest region is identified considering their niching measures. Its features of utilizing the niching measures and density estimation for the evolutionary search are inspired from adaptive grid method by PAES II and selection strategy by SPEA II, respectively.

4. Optimum Design of Grillage Systems
The proposed four MOAs are presented in the previous section. These MOAs are evaluated for the optimal designations of grillage systems in this study. For this purpose, the objective functions, their related constraints, and structural analysis procedures are coded in Java. In this regard, the design requirements used in constraints and design optimization procedure are introduced in Sections 4.1 and 4.2.
4.1. Design Requirements of Grillage Systems according to LRFDAISC Ver.13
Grillage systems comprise a number of lateral braced beams. If the beams loaded in plane of lateral system have no sufficient lateral stiffness, then they are buckled out of plane of loading. This case is called lateraltorsional buckling. The lateraltorsional buckling strength varies depending on the unbraced length and compactness of beam plays an important role in the load carrying capacity of beam. If a compact beam determined according to its web and flange dimensions has a sufficient unbraced length, then nominal flexural strength is calculated in an elastic domain, otherwise an inelastic one. In the inelastic case, a short and unbraced beam length causes to yield its outer fibers before attaining elastic buckling load. The formulation of nominal flexural strength M_{n} that is managed by limit states of yielding, lateral torsional, and flange local buckling is presented in the following part as defined in AISCLRFD Ver.13 (see Figure 1). For simplicity, two distinct figures used to depict the lateraltorsional and flange local buckling depicted in AISCLRFD Ver.13 is coarsely combined, but equation numbers corresponding to formulations for limit states are presented in a separate parenthesis.
The limit states of yielding of beam crosssection are written as
In elastic and inelastic domains, two unbraced lengths and are used to determine the compactness of sections manage the strength of lateraltorsional buckling. Nominalflexural strength is computed as follows.
For inelastictorsional buckling limited by ,where,
For elastictorsional buckling occurred in a segment limited by (see Figure 1), where,
In (6), moment modification factor is utilized to consider the effect of lower torsional buckling arisen from a nonuniform distribution of moment. Therefore, associating with moment diagrams, it is computed as follows.
The local buckling of flanges for noncompact section is governed by two parameters, , . Depending on these parameters, nominalflexural strength is computed as follows.
For inelasticflange buckling occurred in a segment limited by ,
For elasticflange buckling occurred in a segment limited by ,
In addition, nominalshear strength is computed as follows (see the limit states of shear strength in Figure 2): where, In (11), webplate buckling coefficient is equal to 5 for unstiffened webs with .
4.2. Design Procedure for Optimal Grillage Systems
In this work, a design problem of grillage system is represented by two objectives, entire weight of grillage system, and joint deflection, and expressed as follows: Subject to
The term is total weight of all grid members and computed using and which are unit weight to be selected from Wsections list of LRFDAISC Ver.13 and length of a grid member. While is termed as a joint displacement corresponding to related degree of freedom which is denoted , total numbers of joint and grid member are indicated by and . is taken as (max span/300). In constraint inequalities, while displacements of joints are constrained by an upper limit , bendingmoment strength of grid members is limited by allowable nominalmoment strength . Shear strength of grid members is limited by allowable nominalshear strength (see in (11)). In (14) and (15), and are resistance factors for moment and shear and taken as 0.9. Furthermore, it must be noted that indicated in (14) manages total three strength requirements: yielding (2), flangelocal buckling (9) and (10), and lateraltorsional buckling (3) and (5).
5. Search Methodology
Due to the fact that application problems are chosen from realworld engineering design problems with design variables of discrete type, a reliable and consistent search strategy must be established for an evaluation of proposed MOAsโ search capabilities. After generating optimal designations, the computational performances of MOAs are assessed according to the closeness of these designations to a pareto front known beforehand thereby using a number of quality measuring metrics. Furthermore, accuracy of assessment of multiobjective optimization algorithms must be confirmed by outcomes of statistical tests. The other important difficulty is how to adopt a common methodology for a conventional mathematical model, which is performed its computational procedures by usage of design variables of continuous type. Therefore, a reasonable approach is to obtain a pareto front by running current discrete optimization model in bigger and repeated generation numbers. In this regard, independent 10 runs of proposed MOAs are executed by use of both increased and decreased size of generation and population. Then, optimal designations obtained are utilized in computation of qualitymeasuring metrics, such as hyper volume and generational distance. Due to the nature of stochastic algorithms, a statistical test for analysis of these qualitymeasuring metrics computed must be performed with a certain level of confidence. Moreover, in order to decrease the effect of parameter values of evolutionary operators on MOAsโ performance evaluation, various parameter combinations are also considered. Details about qualitymeasuring metrics and statistical tests employed for MOAsโ performance assessment are presented in Sections 5.1 and 5.2.
5.1. QualityMeasuring Metrics
Differentiation in MOAs architecture prevents to lay down the different aspects of MOAsโ performance. Therefore, quality indicators have a big impact on accurately evaluation of MOAs performance. In this study, three quality indicators, hyper volume ratio (Hv), inverted generational distance (IGD), and spread are employed.
Hyper Volume Ratio
Hyper volume (HV) is an indicator which defines a volume covered by Nondominated set of solutions included in a region bounded by a pareto front (see (17)). For this purpose, a reference point chosen among Nondominated solutions with worst objective is utilized in computation of hypercube of each Nondominated solutions,
Hyper volume ratio (HVR) is an indicator which shows a ratio of current hyper volume to the true hyper volume computed by use of true pareto front and nondominated solutions true pareto front (see (18)),
Higher values of HVR indicate a large coverage of Nondominated solutions in a solution space.
Inverted Generational Distance
Inverted generational distance (IGD) estimates the far of Nondominated solutions included in current pareto front generated by the proposed MOA, from those included in true pareto front (see (19)),
where is number of Nondominated solutions found by proposed MOA and is Euclidian distance between each of these and nearest member of true pareto front. A lower value of IGD indicates an increase in approximation of current pareto front obtained to the true pareto front in terms of convergence.
Spread
This quality metric abbreviated as () is used to measure an expanding spread exhibited by Nondominated solutions obtained and computed as
where is Euclidian distance between consecutive Nondominated solutions, is mean of these distances, and and are distances to extreme solutions of current pareto front. A lower value points out a better distribution among Nondominated solutions. In other words, it is implied that Nondominated solutions are located in different positions.
5.2. Statistical Tests
The quality indicators mentioned above are utilized in comparison of MOAsโ distribution qualities. Therefore, after computing means and standard deviations of quality indicators obtained by in the end of independent 10 executions, consistency of these results are checked through performing a statistical analysis in a certain level of confidence. If a probability value resulted from a statistical testing procedure satisfies a userdefined significance level, then it is said that distribution of current MOAs approximation set is acceptable.
Computational procedures of the statistical analysis are performed in MATLAB [42]. Firstly, lillie test is carried out to check quality indicator values for whether to be exhibited a normal distribution (if the null distribution is completely specified, then KolmogorovSmirnov test is more appropriate). Then, existence of a variance homogeneity is controlled through Leveneโ test. If homogeneity in the variance exists, Welch test is performed, otherwise Anova test. In order to compare the statistical outputs acquired from different algorithms, a post hoc testing is performed through โmulticompareโ function coded in MATLAB.
6. Discussion of Results
It is known that the computational performances of MOAs vary depending on the values of their interacting parameters. Therefore, in order to provide an accurate evaluation for their computational performances, some combination sets of parameter values are chosen. Due to a variety on the operator parameters of MOAs, various parameter combinations and the assigned values for their evolutionary operators are summarized in Table 1. The combination numbers (C. No.) are denoted by superscripts attached to related parameter values. In order to provide unbiased competition for the performance evaluation of MOAs employed, the parameter values of evolutionary operators are kept for each MOA. In this regard, the parameter combinations are sorted into two main groups denoted by numbers (1โ4) and (5โ8). Furthermore, each of main groups contains both upper and lower value sets of mutation and crossover distribution indexes in order to provide an intensive mutation or crossover effect for evolutionary search. For example, in first case of AbYSS (shown by italic characters in Table 1), Crossover and Mutation Probability values are 1.00 and 1.00; the values of Crossover and Mutation Distribution Indexes are 30 and 30; the values of Reference set IโII and archive size are 20, 20, and 40. Hence, reproducibility of related MOAs by use of these parameters is ensured. The detailed descriptions of these parameters are found in Section 3.
The penaltyrelated parameters used by penalty functions of weight and displacement are taken as and 0.005, and 2, and and 1 for and (see (17) and (18)), respectively.
Using these different parameter combinations, an optimal design of grillage systems is carried out by four MOAs, NSGA II, SPEA II, PESA II, and AbYSS according to optimum design procedure mentioned previously. The yield stress, elastic modulus, and shear modulus of steel material used to construct the grillage system are taken as 50โksi (345โMPa), 29,000โksi (200โMPa), and 14,500โksi (100โMPa). Crosssectional properties of grid members are chosen from a discrete set with 274 Wsections. Sequence number of each crosssections included in this profile list database is the same as given in LRFDAISC Ver.13 and contains all sectional properties (such as area, inertia moments in all directions). Design variables are represented by binary strings. Thus, a binary length of with possible gene combinations will be adequate to represent 274 ready sections. Profile list database, finite element attributes of design examples and a structural analysis formulation for the grillage system are coded in Java in order to appropriately compile with the optimization tool named JMETAL. In this regard, a userdefined input data coded in Java Language and related parts are presented for the design optimization of example 1 by utilizing AbYSS algorithmbased optimizer (see Algorithm 5). Details about class names and their related packages in Algorithm 4 are given in Section 3.

Design examples of realworld engineering structures are presented in an order of increasing size of their elements and joints. In order to display the change in the strength of grid members, both penalized and unpenalized optimal results corresponding to its maximum joint displacement or weight of grid structure are presented considering the joint and member numbers. Furthermore, convergence history obtained in the end of a complete evolutionary search is also presented. In order to provide an easier visualization for both weight and displacement values, the weight and displacement values are displayed for each of 100 equal segments obtained by a division of the maximum evolution number. Thus, the optimal designation located in one of these segments is easily determined.
6.1. Design Example 1: A Grid System with 4 and 3 Bays
This simple grillage system depicted in Figure 3 has four longitudinal and three lateral bays. Grid members of grillage system in and directions are linked into two separate groups: grid members (1โ16) as design variable 1, grid members (17โ32) as design variable 2, grid members (33โ44) as design variable 3, grid members (45โ56) as design variable 4, grid members (57โ68) as design variable 5. Lengths of spans are ft (3.2004โm), โft (3.5052โm), โft (3.6576โm), โft (3.6576โm), โft (4.0386โm), โft (3.81โm), and โft (3.048โm). Magnitudes of loads are taken as 17.9847โkipf (80.00โkN) for and , 20.2328โkipf (90.00โkN) for and , 17.9847โkipf (80.00โkN) for and , 15.7366โkipf (70.00โkN) for , and 16.8606โkipf (75.00โkN) for . The value of allowable displacement is constrained as (13.25/300 = 0.0441โft; 13.4416โmm).
The pareto fronts of a set of random Nondominated solutions obtained by both use of an increased population size and evolution number and AbYSS, NSGAII, PESAII, and SPEAII are depicted in Figures 4(a) and 4(b1)โ4(b4). The means of quality indicator values, standard deviations values corresponding to best means, and assessment of significance levels according to statistical tests are presented in Tables 2 and 3. Considering Tables 2 and 3, an increase in the evolution number and population size leads to an improvement in the values of quality indicator, in other words, an increase in the convergence degree of optimal designations. According to assessment of statistical analysis results corresponding to the lower evolution number and population size, it is clear that there is not a statistical confidence among proposed MOAs. The best values of quality indicators are obtained by AbYSS. SPEAII and PESAII exhibits better computational performance compared to NSGAII. It is also observed that an increase in the evolution number and population size forces MOAs to use higher values of operator parameter for improvement of their optimal designation qualities. This result is confirmed by examining the values of quality indicators with bold and italic characters. These higher quality indicator values are obtained by Casesโโ1โ4, which indicate a usage of lower evolution number and population size and Casesโโ5โ8, which indicate a usage of higher evolution number and population size. Considering the best spread values, the pareto fronts obtained by four MOAs employed are compared with both each other and the pareto front obtained by use of higher population size and evolution number (Figures 4(b1)โ4(b4)). The convergence history of the optimizer AbYSS (Caseโโ3) with the best spread value 0.8067 is presented in Figure 5 with two axes. A designation with (weight = 5446.1415โlb (2470.3282โkg) and displacement = 0.0406โft (12.3748โmm)) is obtained in a generation of no. = 4638 located in a segment of no. 92 with an interval .
The strength values of grid members are displayed considering penalized designation corresponding to a maximum displacement (weight = 18128.2500โlb (8222.835โkg) and displacement = 0.3977โft (121.2189โmm); see Figures 6(a1)โ6(a5)) and unpenalized designations corresponding to a maximum displacement (weight = 17095.5000โlb (7754.3883โkg) and displacement = 0.0068โft (2.0721โmm)) and a maximum weight (weight = 74672.5000โlb (33870.8762โkg) and displacement = 0.0002โft (0.0609โmm); see Figures 6(b1)โ6(b5) and 6(c1)โ6(c5)). Whereas the penalized designation corresponding to the maximum displacement contains a Wsection set (W30ร292, W6ร12, W27ร235, W6ร12, and W8ร15), the unpenalized designations corresponding to the maximum displacement and weight are represented by Wsection sets (W24ร68, W30ร124, W27ร129, and W30ร148) and (W14ร370, W14ร455, W36ร529, W36ร652, and W36ร361), respectively. The most critical strength values of penalized designation satisfied none of the constraints is obtained by both torsional and flange bucklingrelated constraints (see Figures 6(a1)โ6(a5)). Examining the strength values depicted in Figures 6(b1)โ6(b5) and 6(c1)โ6(c5), limit state values of unpenalized designations corresponding to the maximum weight are higher than one corresponding to the maximum displacement due to a usage of the bigger crosssectional properties for design variables. Furthermore, a common point of these strengths obtained by both panelized and unpenalized solutions corresponding to the maximum displacement is that the yielding limit state has a big impact on the flange buckling, torsional and yielding strengths (see Figures 6(a1), 6(a4), and 6(a5) and 6(b1), 6(b4), and 6(b5)).
6.2. Design Example 2: A Grid System with Five Bays
A grillage system with 160 grid members and 152 nodes is shown in Figure 7. This grillage system is almost 10times larger than first design example. Grid members are linked in 8 separate groups, resulting in total four design variables in direction and total four design variables in direction. According to this linkage system, members (1โ20), (21โ40), (41โ60), (61โ80), (81โ100), (101โ120), (121โ140), and (141โ160) are linked to represent design variables respectively. Magnitude of loads โ is equal to 49.46โkipf (219.998โkN). Lengths of spans are โft (2.200โm), โft (2.505โm), โft (2.398โm). The value of allowable displacement is constrained as (8.22/300 = 0.0274โft; 8.3515โmm).
The pareto fronts of a set of random Nondominated solutions obtained by both use of an increased population size and evolution number and AbYSS, NSGAII, PESAII, and SPEAII are depicted in Figures 8(a) and 8(b1)โ8(b4). The means of quality indicator values, standard deviations values corresponding to best means and assessment of significance levels according to statistical tests are reported in Tables 4 and 5. Examining the values of quality indicators tabulated in Tables 4 and 5, it is seen that the lower evolution number and population size decrease the quality degree of optimal designations leading to obtain an inconsistent statistical confidence among quality indicators. However, the quality degree of optimal designations is elevated associated by an increase in the evolution number and population size. According to the quality indicator values, AbYSS (Caseโโ1) succeeded to obtain a lower spread 0.8356 and a higher hypervolume value 0.9419 compared to NSGAII, PESAII, and SPEAII. It is also observed that a higher optimality degree is obtained by use of either lower parameter values of evolutionary operators in conjunction with a decreased evolution number and population size or higher ones in conjunction with an increased evolution number a population size. This claim is approved by examining the quality indicator values with bold and italic characters obtained by Casesโโ1โ4, which indicate a usage of lower evolution number and population size and Casesโโ5โ8, which indicate a usage of higher evolution number and population size (see Tables 4 and 5). Considering the best spread values, the pareto fronts obtained by four MOAs employed are compared with both each other and the pareto front obtained by use of higher population size and evolution number (Figure 8(b1)โ8(b4)). The convergence history of the optimizer AbYSS (Caseโโ1) with the best spread value 0.8067 is presented in Figure 9 with two axes. A designation with (weight = 13584.2465โlb (6161.7105โkg) and displacement = 0.0798โft (24.3230โmm)) is obtained in a generation of no. = 4257 located in a segment of no. 85 with an interval .
The strength values of grid members are displayed considering penalized designation corresponding to a maximum displacement (weight = 9888.9999โlb (4485.5749โkg) and displacement = 0.1633โft (49.7738โmm); see Figure 10(a1)โ10(a5)) and unpenalized designations corresponding to a maximum displacement (weight = 32718.4000โlb (14840.8165โkg) and displacement = 0.0125โft (3.81โmm)) and a maximum weight (weight = 102588.7500โlb (46533.4742โkg) and displacement = 0.0025โft (0.762โmm); see Figures 10(b1)โ10(b5) and 10(c1)โ10(c5)). Whereas the penalized designation corresponding to the maximum displacement contains a Wsection set (W14ร38, W12ร45, W12ร53, W18ร55, W12ร22, W6ร8.5, W6ร20, and W18ร71), the unpenalized designations corresponding to the maximum displacement and weight are represented by Wsection sets (W12ร87, W24ร131, W33ร152, W10ร112, W24ร131, W12ร170, W21ร132, W14ร132) and (W40ร149, W12ร210, W36ร441, W36ร529, W36ร395, W36ร330, W36ร361, and W27ร194), respectively. An increase in the displacement values causes convergence in the flange bucking, torsional, and yieldingrelated strength values to their related yielding limit values. Hence, the distribution of these strengths on the grid members becomes more close to their limit state values (see Figure 10(a1), 10(a4) and 10(a5) and 10(b1), 10(b4), and 10(b5)). An increase in the weight of grid structure leads to an elevation in the limit state values (see Figure 10(c1)โ10(c5)).
6.3. Design Example 3: A Grid System with Five and Six Bays
This grillage system with 236 members and 219 joint points has the highest complexity among design examples due to both lower number of support points constructed to carry a large grillage area and higher number of members and joints (see Figure 11). Grid members (1โ20), (21โ40), (41โ60), (61โ80), (81โ100), (101โ120), and (121โ140) are linked to represent design variables in direction while design variables , and 11 in direction are assigned to lattice beams which are denoted by (, and ), (, and ), (, and ), and (, and ). Magnitude of loads (โ and โ), (โ and โ), (โ and โ), and (โ) are taken as 24.730โkipf (109.999โkN), 35.970โkipf (159.995โkN), 38.210โkipf (169.959โkN) and 42.710โkipf (189.974โkN). Lengths of spans are equal to โft (2.30โm), โft (2.399โm), โft (2.499โm), โft (2.599โm), and โft (2.798โm). The value of allowable displacement is constrained as (9.18/300 = 0.0306โft; 9.3268โmm)
The pareto fronts of a set of random Nondominated solutions obtained by both use of an increased population size and evolution number and AbYSS, NSGAII, PESAII, and SPEAII are depicted in Figures 12(a) and 12(b1)โ12(b4). The quality indicator relatedquantities including statistical test results are reported in Tables 6 and 7. According to tabulated valued in Tables 6 and 7, a decrease in the evolution number and population size causes the quality indicator values to be poor and hence statistical confidence to be inconsistent. Considering the indicator values in Table 7, it is obvious that AbYSS (Caseโโ1) shows a better computational performance by obtaining a lower spread 0.8168 and a higher hypervolume value 0.9657 compared to NSGAII, PESAII, and SPEAII. Also, the quality indicators obtained by PESAII and SPEAII are better than NSGAII. According to the convergence history of AbYSS (Caseโโ1) with the best spread value 0.8168 presented in Figure 13, a designation with (weight = 22354.1046โlb (10139.6512โkg) and displacement = 0.2331โft (71.0488โmm)) is obtained in a generation of no. = 4461 located in a segment of no. 89 with a interval .
The strength values of grid members are displayed considering penalized designation corresponding to a maximum displacement (weight = 36375.3700โlb (16499.5902โkg) and displacement = 0.7651โft (233.2024โmm); see Figure 14(a1)โ14(a5)) and unpenalized designations corresponding to a maximum displacement (weight = 44028.4600โlb (19970.9735โkg) and displacement = 0.02142โft (6.528โmm)) and a maximum weight (weight = 176334.6800โlb (79984.0654โkg) and displacement = 0.0066โft (2.0116โmm); see Figures 14(b1)โ14(b5) and 14(c1)โ14(c5)). Whereas the penalized designation corresponding to the maximum displacement contains a Wsection set (W16ร40, W12ร45, W30ร391, W14ร30, W12ร30, W12ร19, W8ร15, W12ร96, W18ร175, W10ร39, and W30ร116), the unpenalized designations corresponding to the maximum displacement and weight are represented by Wsection sets (W18ร106, W18ร106, W18ร106, W12ร170, W18ร106, W18ร106, W18ร106, W18ร97, W21ร122, W21ร122, and W18ร97), and (W40ร503, W12ร305, W14ร500, W14ร730, W14ร550, W40ร278, W16ร100, W36ร361, W14ร665, W40ร593, and W14ร311), respectively. Examining Figure 14(a1)โ14(a5), it is seen that torsional and yielding strengths play a bigger and more important role in the bearing capacity of grid structure compared to the other strengths. A similar result as in the application of design examples 1 and 2 is obtained in this application example: a decrease in the displacement values causes the flange bucking, torsional, and yieldingrelated strengths to obtain highly near to their related yielding limit state values. Using the steel profiles with larger crosssectional properties leads to an increase in both weight of grid structure and limit state values.
7. Conclusion
In this study, the designs of grillage systems are optimized using the optimization tools named NSGAII, SPEAII, PESAII, and AbYSSS according to the design provisions of LRFDAISC Ver.13 specification. Hence, the computational performances of the MOAs employed are not only evaluated but the change in the strength of gird members in conjunction with joint displacements is also displayed. In order to assess the computational performance of MOAs, the various quality indicators, IGD, HV, and spread are computed considering the application of three design examples with an order of increasing complexity degree and evaluated their statistical confidences according to outcomes of statistical tests. According to the evaluation of quality indicators and optimal designations, the following results are observed.(i)According to the values of quality indicators, AbYSS shows better computational performance compared to SPEAII, PESAII, and NSGAII. Moreover, NSGAII obtains the worst values of quality indicators than SPEAII and PESAII, whose computational performances are almost equal.(ii)An efficient exploration of search space is only provided by usage of an increased evolution number and population size. Otherwise, the outcomes from statistical tests point out an inconsistent relation between quality indicators obtained by use of a decreased evolution number and population size.(iii)An increase in the parameter values of evolutionary operators in conjunction with the evolution number and population size leads to an improvement in the convergence degree of optimal designations.(iv)In conjunction with increasing the evolution number and populations size, the degrees of quality indicators are improved by use of (higher crossover probability and crossover distribution indexโhigher mutation probability and mutation distribution indexโhigher reference sets and archive size) for the evolutionary operators of AbYSS, and (lower crossover probability and crossover distribution indexโlower mutation probability and mutation distribution indexโhigher archive size) or (higher crossover probability and crossover distribution indexโlower mutation probability and mutation distribution indexโhigher archive size) for the evolutionary operators of SPEAII and PESAII. (v)A decrease in the joint displacement correspondingly causes the values of flange buckling, torsional, and yielding strengths to obtain near to their related limit state values. In this regard, the bearing capacity of grid structure has to be kept in the upper levels against any drastic decrease in joint displacements. Therefore, using a multiobjective optimization algorithm becomes a reasonable approach for the optimal design of grillage systems.
In the next work, the effect of using the different evolutionary parameter values on the quality degree of optimal designations will be evaluated. Furthermore, the proposed MOAs will be hybridized with both each other and the other evolutionary optimization algorithms in order to both reduce the evolutionary computational cost spent to obtain Nondominated solutions and improve the curves of pareto fronts.
Nomenclature for AISC_LRFD Ver.13 Specification
Flange width  
:  Lateraltorsional buckling modification factor for nonuniform moment diagrams 
:  Elasticity modules 
:  Critical stress 
:  Specified minimum yield stress 
:  Clear distance between flanges 
:  Distance between flange centroids 
:  Outofplane moment of Inertia 
:  Torsional constant 
:  Coefficient for slender unstiffened elements 
:  Webplate buckling coefficient 
:  Limiting laterally unbraced length for the limit state of yielding 
:  Distance between braces 
:  Limiting laterally unbraced length for the limit state of inelastic lateraltorsional buckling 
:  Nominal flexural strength 
:  Plastic bending moment 
:  Absolute value of monet at quarter point of unbraced segment 
:  Absolute value of monet at centerline point of unbraced segment 
:  Absolute value of monet at threequarter point of unbraced segment 
:  Radius of gyration about axes 
:  Effective radius of gyration used in the determination of for lateraltorsional buckling limit state 
:  Elastic section modules about principal axes 
:  Thickness of flange 
:  Web thickness 
:  Plastic section modules about principal axes 
:  Limiting slenderness parameter for compact element 
:  Limiting slenderness parameter for compact flange 
:  Limiting slenderness parameter for compact web 
:  Slenderness parameter. 
Acknowledgments
The author would like to thank Professor Dr. Antonio J. Nebro, Prof. Dr. Juan J. Durillo, and their team members for releasing the source code of JMETAL to execute computational procedures of multiobjective optimization methods. He thanks reviewers for their valuable comments and recommendations about revision of this paper. He also thanks Kadirli Municipality of Osmaniye for their support to this paper.