Optimization of slender reinforced concrete columns subjected to biaxial bending using genetic algorithms

: In this paper, a computational tool was developed to optimize the design of slender reinforced concrete columns subjected to biaxial bending considering the material and geometric nonlinearities rigorously. The optimization process utilizes the technique of genetic algorithms to find the best cross- sectional dimension and the best distribution and amount of reinforcement, to minimize the cost of the column subject to certain constraints of strength, stability, feasibility and regulatory. The analysis applies to rectangular cross-section of columns


INTRODUCTION
In an earlier paper, Pires and Silva [1] developed a numerical procedure for the design of slender columns of reinforced concrete subjected to uniaxial bending. In the design of slender columns, displacements arising from secondorder effects gradually increase until the column find a deformed position to establish the equilibrium condition. However how to find the steel area of a column, whose internal loads are not known? How to find the second-order internal loads of a column whose steel area is unknown? To answer these questions, the authors developed an iterative process that establishes the smallest area of steel necessary to equilibrate the column, based on a previously chosen reinforcement distribution.
A natural progression for this procedure would be to develop it to design slender columns subjected to biaxial bending. However, there would be nothing new on this subject, as many researchers have studied the nonlinear behavior of slender reinforced concrete columns subjected to uniaxial and biaxial bending [2]- [12].
Therefore, the idea arose to optimize the design of slender columns considering the concrete area, the steel area, and the reinforcement distribution as variables, so that the column is equilibrated in the safest and most economical possible way.
These studies confirm the efficiency of GA for reinforced concrete structures. In particular, studies found on the optimization of reinforced concrete columns, isolated or as frame members, provide answers using, in general, the design solution by simplified methods. The aim of this paper is to expand the use of AG for reinforced concrete columns by introducing a more comprehensive procedure in which material and geometric nonlinearities are strictly considered. This procedure allows the analysis of reinforced concrete columns from the shortest to the slenderest with the same rigor.
In the present study, the results provided by GA are compared with the solutions obtained by the procedure called Total Research, which provides the exact optimal results for the design of reinforced concrete columns. It should be noted that the Total Search procedure, although providing the exact solution, leads to an extremely high computational cost, and for this reason this procedure is used only to validate the results obtained by GA, but it is not economically viable as a computational tool.
In summary, this paper presents the systematization of the optimal design of slender reinforced concrete columns. The optimization process uses the technique of genetic algorithms that allow to find the best dimension of the cross section and the best distribution of the reinforcement so that the cost of the column, subject to certain constraints of strength, stability, feasibility, and regulatory specifications, is minimized. Columns are treated as individuals that are coded according to their characteristics, and are assessed by a cost function that is penalized according to the imposed constraints. The optimal column will correspond, at the end of the procedure, to the individual who best meets the safety criteria and, at the same time, has the lowest cost.
The present analysis covers isolated reinforced concrete columns of rectangular cross section submitted to biaxial bending, considering any type of support, and characteristic compressive concrete strengths between 20 MPa and 50 MPa. Material and geometric nonlinearities are strictly considered.

Uniaxial constitutive models
The constitutive model for concrete in compression defined by the Brazilian standard NBR 6118 [29] is considered. The concrete stress-strain relationship is given by Equations 1 to 3: where σ c = compressive stress of concrete; ε c = compressive strain of concrete; and f cd = design compressive strength of concrete for concretes from 20 MPa to 50 MPa. The constitutive model for steel adopted by the Brazilian Standard NBR 6118 [29] applies to tension and compression, and the stress-strain relationship is given by Equations 4 and 5: where σ s = stress of steel; E s = elastic modulus of steel; ε s = strain of steel; f yd = design strength of steel, and ε yd = yield strain of steel.

Considerations on stability of slender reinforced concrete columns subjected to biaxial bending
The design of slender reinforced concrete columns requires the study of the element's stability, which implies in an analysis of the second order effects.
A bar subjected to biaxial bending caused by an eccentric load F d , with an eccentricity e1, will have its axis deformed. In the case of slender bars, the second order eccentricity (e 2 ) is added to the transverse displacements. Axial loading and biaxial bending considering second order effects are given by Equations 6 to 9: In Equations 6 to 9, the numerical values used for the variables correspond to respective vector modules. When the axial force N d and the bending moment Md (components M zd and M yd ) are considered as ultimate values (N R , M RZ , M RY ), they produce strains and stresses corresponding to the ultimate limit state of the cross section. The relationships between the ultimate values (N R , M Rz , M Ry ) and the stresses in concrete and steel are given by Equations 10 a 12:  (12) It is necessary to know the deformed configuration of the column to solve Equations 10 to 12, which implies knowing the curvature of the column axis and the stresses and strains in each column cross section.

Assumptions and ultimate limit states
The assumptions admitted in this paper are: (1) the plane sections remain plane after the element is strained; (2) the strain in a generic fiber of the cross section is directly proportional to its distance from the neutral axis; (3) there is a perfect bonding between the reinforcement bars and the concrete that surrounds them; (4) the tensile strength of concrete is totally neglected; (5) the parabola-rectangle diagram is adopted to represent the stress-strain relationship of concrete; (6) a perfect elastic-plastic diagram is adopted to represent the stress-strain relationship of steel; (7) the concrete cross section is considered to be non-cracked when calculating the stiffness matrix; (8) the conventional value of 0.35% is taken as the strain limit for shortening concrete in partially compressed sections; concrete strain limits vary from 0.2% to 0.35% under non-uniform compression. In this case, the strain related to the fiber positioned at 3/7 of the total section height from the most compressed edge remains unchanged and equal to 0.2%; (9) the maximum elongation allowed in the tensile reinforcement is 1%.
The small displacements assumption is considered. Thus, the curvature (χ) of the column axis is given by where w = transverse displacement. The variable x is measured along the non-deformed axis, and the axial load remains constant regardless of column deformations.
The safety of reinforced concrete columns is verified for the following ultimate limit states: (1) the ultimate limit state of loss of equilibrium of the structure; (2) the ultimate limit state related to the load capacity of the structure considering the second order effects.

Calculation of displacements
A local coordinate system for the bar is adopted to define the displacements (Figure 1a). The bar is subjected to biaxial bending. Internal forces occur in the XZ and XY planes. The axis of the bar is displaced u 0 (x) in the x direction. In the XZ plane, the bar undergoes a displacement w(x) and the cross-section undergoes a θ z (x) rotation ( Figure 1b). In the XY plane, the bar undergoes a displacement v(x) and the cross-section undergoes a θ y (x) rotation ( Figure 1c). The transverse displacements w(x) and v(x) are positive in the direction of the local axis, and the bending rotations, θ z (x) and θ y (x), are positive in the clockwise direction. The displacement u(x, y, z) in a generic fiber of the cross section, according to Chen and Atsuta [30], is given by Equation 14: The longitudinal strain is given by Equation 15: The axial strain ε o , the χ y curvature in the XY plane, and χ z curvature in the XZ plane are given by Equations 16, 17 and 18, respectively: In the expression for the longitudinal strain (ε x ), the non-linear relationship is given by the axial strain (ε o ).
In the analysis of slender columns, material and geometric nonlinearities must be considered so that the displacements are calculated accurately. The numerical procedure developed for calculating displacements is based on the Finite Element Method using the Virtual Work Principle. The column axis is discretized into small elements of length L that are connected by their nodes. Each node has an axial displacement, a transverse displacement in the z direction, a transverse displacement in the y direction, a rotation in the XZ plane, and a rotation in the XY plane. The axis displacements of the element (uo, w, v) are given by Equations 19, 20 and 21, respectively: The interpolation functions ϕ i are given by Equations 22 to 27: The nodal nonlinear forces in each element ( ) are given by Equations 28 to 37: 2 4 3 5 5 9 6 10 6 10 0 0 The Gauss-Legendre quadrature rule is used to solve these equations. The nodal displacements, the six interpolation functions, and the ten equilibrium equations, combined, establish the system of Equation 38: where Fn (e) = load vector; K (e) = stiffness matrix; and U (e) = displacement vector of the e th finite element. Geometric nonlinearity is included in the last term of Equations 29 to 32 and Equations 34 to 37. The modified Newton-Raphson method with constant stiffness is used to solve the previously established nonlinear equations.
Material nonlinearity is considered when calculating the axial load N d given by Equation 10, and when calculating the bending moments M zd and M yd obtained by Equations 11 and 12, respectively. The procedure proposed by Campos [31] is used to calculate these internal loads.
The global equation system is made by assembling the respective vectors and matrices of the finite elements that make up the column, and is given by Equation 39: where Fn = global load vector; K = global stiffness matrix; and U = global displacement vector.

Procedure for the optimal design of the column
The aim of optimization is to select the lowest cost column that meets the criteria established by the Brazilian Standard NBR 6118 [29]. Optimization of reinforced concrete column cross sections starts from a set of randomly created cross sections, and is carried out over generations. During the evolutionary process, individuals (cross sections) are coded by a set of bits (binary coding), and represented by chromosomes. Each value assigned to variables b, h, dicam13, dicam24, ncam13 and ncam24 is related to an integer and corresponds to a binary number according to Table 1.  Table 2 shows an example of coding the cross section of reinforced concrete columns, and Figure 3 illustrates the resulting chromosome.   After coding, each individual is evaluated during the evolutionary process by means of a fitness function, Fa (x). This function is used to rank the best solutions, indicating the chances of survival and permanence of good features over generations. The fitness function, Fa (x), is defined by Equation 40: where F(x) = objective function; and P(x) = penalty function.
In order to work with the term fitness of the individual, the cost minimization problem is transformed into a fitness maximization problem. The objective function F(x) is defined as the inverse of the cost function Cost(x) in Equation 41: The cost function Cost(x) is given by Equation 42: where C c = cost of concrete per unit volume; C s = cost of steel per unit volume; C f = cost of the formwork per unit area; l = length of the column; b and h = dimensions of the column cross section; and A s tot = total steel area of the cross section of the column. The values for C c , C s and C f were obtained on the website of the Foundation for the Development of Education [32]. The penalty function P(x) represents the different constraints of the problem, and is given by Equation 43: where α = coefficient used to control the degree of penalty; t = current generation; r = number of inequality constraints; and m = number of equality constraints.
In this paper, there are no equality constraints h j (x), and there are seven-inequality constraints g i (x), namely: the maximum strains of concrete and steel, the stability of the column, the minimum spacing of the reinforcement, the maximum spacing of the reinforcement, the minimum reinforcement rate, and the maximum reinforcement rate.
The magnitude of the coefficient α, in Equation 43, defines the degree of penalty that will be imposed on the individual. The degree of penalty must be sufficient to make the best feasible solution of the generation have the highest fitness value after the penalty. In this paper, the procedure proposed by Nanakorn and Meesomklin [33] is used to adjust the coefficient α.
In the optimization process, whenever an individual violates a constraint, whether strength, stability, or regulatory, that individual reduces its chances of being selected for the next generation. The individuals selected for the next generation are recombined or mutated and will form a new population that, in turn, will be used as input for the next iteration of the optimization algorithm. This procedure is repeated until a solution is found that satisfies the convergence criteria. The developed algorithm was implemented in a computer program called GENETIC ALGORITHM.
A second computer program called TOTAL SEARCH was developed to validate the GENETIC ALGORITHM computer program. The role of TOTAL SEARCH is to obtain the optimal design solution for a given column, using all possible combinations of variables, which in this paper is a total of 1,048,576 individuals (Figure 2b). In the TOTAL SEARCH program, each individual is checked once again for the same constraints used in the GENETIC ALGORITHM program. Among all the individuals, the one that has the lowest cost and satisfies all the constraints is chosen.
The TOTAL SEARCH program provides the exact answer to the problem, but its processing time is significantly longer than the time used by the GENETIC ALGORITHM program. Therefore, the results obtained by the TOTAL SEARCH program are only used as reference values to confirm the efficiency of the GENETIC ALGORITHM program.
A third computer program called COLUMN PROCESSING was developed to check the constraints of the column regarding its strength, stability, and compliance with normative specifications.
The development and details of the TOTAL SEARCH, GENETIC ALGORITHM and COLUMN PROCESSING computer programs are found in Pires [34].

About developed programs
The development of the GENETIC ALGORITHM program is derived from the LGADOS program developed by Coley [35]. In present paper, all necessary adjustments were inserted for adapting the LGADOS program to the proposed problem.
The analysis of the constraints in the GENETIC ALGORITHM program is done by two approaches. The first approach applies the Death Penalty Method, in which individuals (cross sections of the column) that do not meet one or more constraints are eliminated in the next generation. The second approach considers the inclusion of individuals (cross sections of the column) that have violated one or more constraints. This second approach applies the Penalty Method according to the procedure proposed by Nanakorn and Meesomklin [33].
The options adopted for the development of the GENETIC ALGORITHM program are presented below: selection -roulette; crossing -1 point (randomly selected point); elitism -yes (one individual); number of generations -100; stopping criterion -number of generations; population size (NPOP) -60; crossover probability (P c ) -0.6; mutation probability (P m ) -0.02. For the Death Penalty Method, linear scaling is adopted, with a scale constant (C) equal to 2. For the Penalty Method, bilinear scaling is adopted, with constants C, φ and Z equal to 2, 1 and 5, respectively.
The COLUMN PROCESSING program calculates the ultimate load for a reinforced concrete column subjected to biaxial bending, considering the material and geometric nonlinearities rigorously. This program also checks the constraints imposed on the objective function F(x).

Analysis of the COLUMN PROCESSING computer program
The COLUMN PROCESSING program, before being used as a computational tool to support the optimization programs GENETIC ALGORITHM and TOTAL SEARCH, had its efficiency verified through a comparative analysis. Two experimental studies were used: Kim and Yang [36], and Claeson and Gylltoft [37]. The following variables were considered in this comparative analysis: c is the concrete cover of the column; L is the length of the column; fc is the compressive strength of concrete; f y is the yield strength of steel; P u,a is the ultimate load calculated by COLUMN PROCESSING program; P u,t is the ultimate load of the experimental study.
Kim and Yang [36] performed a series of tests on reinforced concrete columns to verify the effects of concrete compressive strength, slenderness, and reinforcement ratio on ultimate load and on the relationship between axial force and bending moment. The authors used different values for concrete compressive strength f c (25.5 MPa, 63.5 MPa, and 86.2 MPa), reinforcement rate ρ (1.98% and 3.95%), and slenderness ratio λ (10, 60, and 100). In this comparative analysis, two columns (1 and 2) were considered, and only concrete compressive strength equal to 25.5 MPa was used. The rate of change (Δ) of the values of the ultimate loads calculated by the COLUMN PROCESSING program in relation to the experimental ultimate loads obtained by Kim and Yang [36] ranged from 0.89 to 1.09. Table 3 presents the characteristics and results of those analyzed columns. Claeson and Gylltoft [37] carried out an experimental study to analyze the behavior of reinforced concrete columns. The authors analyzed the effects of slenderness and concrete compressive strength on the ultimate load of the columns. The rate of change (Δ) of the values of the ultimate loads calculated by the COLUMN PROCESSING program in relation to the experimental ultimate loads obtained by Claeson and Gylltoft [37] ranged from 0.95 and 1.07. Table 4 presents the characteristics and results of those analyzed columns. The comparative study showed that there is good agreement between the experimental results and the values obtained by the COLUMN PROCESSING program.

Analysis of the GENETIC ALGORITHM program
Three columns named P1, P2 and P3 were used to evaluate the effectiveness of the Genetic Algorithm program. Table 5 shows the geometric characteristics of these columns and the characteristics of the materials used. The structural analysis was performed using FEM. Column P1 is subject to uniaxial bending. Columns P2 and P3 are subject to biaxial bending, with column P3 subject to much larger first-order moments (M z and M y ) than those applied to column P2. Each column was composed of 10 bars and 11 nodes. Nodes 1 and 11 are the end nodes. Node 1 was constrained in the x, y and z directions, and node 11 was constrained in the y and z directions. The axial force applied in the x direction at nodes 1 and 11 was equal to 1019 kN at column P1, 1230 kN at column P2, and 2000 kN at column P3. The bending moment M z applied to nodes 1 and 11 of column P1 was equal to 2038 kN.cm. The bending moments M z and M y applied to nodes 1 and 11 of column P2 were equal to 1291.5 and 701.1, respectively. The bending moments M z and M y applied at nodes 1 and 11 of column P3 were equal to 60000 kN.cm.
Initially, the TOTAL SEARCH program is used to calculate the exact global minimum cost (C global ) for each column analyzed.
Next, the GENETIC ALGORITHM program is used according to the following steps: each column is calculated 10 times by applying the Death Penalty Method to the constraint analysis, and the column is again calculated 10 times by applying the Death Penalty Method. The minimum cost values (C min ) are obtained. Tables 6, 7  The results obtained for column P1 show that the relative error (E rel ) ranged from 0 to 3%, with seven values converging to the optimal value calculated by the TOTAL SEARCH program, when the Penalty Method was adopted for the constraints analysis. The relative error (E rel ) ranged between 0 and 2%, and six values converged to the optimum value, when the Death Penalty Method was used.
The results obtained for column P2 show that the relative error (E rel ) ranged from 0 to 4%, with two values converging to the optimal value calculated by the TOTAL SEARCH program, when the Penalty Method was adopted for the constraints analysis. The relative error (E rel ) ranged between 0 and 4%, and two values converged to the optimum value, when the Death Penalty Method was used.
The results obtained for column P3 show that the relative error (E rel ) ranged from 0 to 5%, with one value converging to the optimal value calculated by the TOTAL SEARCH program, when the Penalty Method was adopted for the constraints analysis. The relative error (E rel ) ranged between 0 and 4%, and one value converged to the optimum value, when the Death Penalty Method was used.
In general, it was found that the results obtained by the GENETIC ALGORITHM program converged to the optimal solution or to values close to the optimal solution, taking the values provided by the TOTAL SEARCH program as reference. This finding confirms the efficiency of the GENETIC ALGORITHM program, since the TOTAL SEARCH program provides the exact value of the optimization process. It should also be noted that the processing time spent by the GENETIC ALGORITHM program is significantly lesser than the processing time spent by the TOTAL SEARCH program. The respective processing times spent on the analysis of columns P1, P2 and P3 are shown in Tables 6, 7 and 8. Finally, it was found that no significant differences were found when comparing the results obtained by the two methods used in the GENETIC ALGORITHM program for the constraints analysis, the Penalty Method and the Death Penalty Method. Both methods showed consistent and satisfactory results.

The effect of slenderness on the optimal design of reinforced concrete columns
The effect of slenderness ratio (λ) on the optimal design of reinforced concrete columns was also analyzed. Twelve columns were adopted, based on the characteristics of P2 column. These columns were called P2.i, with i ranging from 1 to 12. Column P2 was selected as a reference because it was subjected to biaxial bending and not very high first-order moments, which allowed the length of the column to be varied until very high values of λ were reached.
Values between 300 cm and 1400 cm were considered for the column length (l), and the smallest dimension (h) of the column cross section was fixed at 25 centimeters in order to induce high values for λ.
The slenderness ratio of columns P2.1 to P2.12 were calculated according to the following criteria of code NBR 6118/2014: columns must have a slenderness ratio less than or equal to 200 (λ ≤ 200); second order local internal loads on isolated members can be neglected when slenderness ratio is less than the limit value (λ 1 ). Slenderness ratio (λ) and limit slenderness ratio (λ 1 ) are obtained by Equations 46 and 47, respectively: where l = column length; h = size of the cross section parallel to the plane of the moment acting on the column; e 1 = smallest value of the first-order eccentricity at the cross section under consideration; and α b = coefficient that depends on the bending moment distribution in the column.  The optimal design results shown in Table 9 were obtained by the TOTAL SEARCH program. This computer program was used in this step because it provides the exact solution in the search space. Figure 4 elaborated from the values in Table 9 shows the relationship between critical slenderness ratio and column cost. In this analysis, the critical slenderness ratio always occurred in the XZ plane, that is, λ z . The slenderness ratio λ z ranged from 41.57 to 193.99, and costs ranged from R$ 546.40 to R$ 9843.42. These relationships show that the more slender a reinforced concrete column is, the higher its cost.
The classification of the columns according to slenderness used in this work is shown in Table 10, and the results found for slenderness ratios and reinforcement rates for columns P2.1 to P2.12 are shown in Table 11.   By initially analyzing the complete set of columns from P2.1 to P2.12, it can be seen that as the slenderness ratio increases, the concrete area (A c ) and the steel area (A s ) increase. However, while A c increases by at most 50% (from 1,000 cm2 to 1,500 cm2), A s increases about 20 times (from 4.71 cm2 to 99.65 cm2). This behavior shows that the increased bending due to second order effects caused by the increased slenderness is now absorbed predominantly by the reinforcement.
It is also possible to analyze columns P2.1 to P2.12 according to the slenderness classification shown in Table 10. The columns were divided into three groups: moderately slender columns, with λ 1 <λ ≤ 90 (P2.1 to P2.4); slender columns, with 90 <λ ≤ 140 (P2.5 to P2.8), and very slender columns, with 140 <λ ≤ 200 (P2.9 to P2.12). In the moderately slender columns group, the concrete area varied around 37% (from 1,000 to 1,375 cm 2 ), while the reinforcement rate (ρ) varied around 54% (from 0.46% to 0.71%). In the slender columns group, the concrete area varied around 9% (from 1,375 to 1,500 cm 2 ) while the reinforcement rate varied more than 400% (from 0.42% to 2.40%). In the very slender columns group, the concrete area varied 10% -very close to the variation in the previous group -and the reinforcement rate varied around 130% (from 3.43% to 7.97%). These results show that in the optimization of moderately slender columns there was an equilibrium between concrete and steel. On the other hand, in the search for the optimal solution for slender and very slender columns, the concrete area oscillated very little, while the steel area played a predominant role. Figure 5 elaborated from the values in Table 11 shows the relationship between the reinforcement rate and the slenderness ratio for columns P2.1 to P2.12.

CONCLUSIONS
In this paper, a procedure was developed to choose, within a range of possible solutions, the column that best meets the requirements of safety, economy and regulation. From this procedure, a computational tool was developed to optimize the cross section of slender reinforced concrete columns subjected to biaxial bending, strictly considering material and geometric nonlinearities, i.e., simplified methods are not used.
The computational tool composed of the GENETIC ALGORITHM and COLUMN PROCESSING programs provided answers that always converged to an optimal solution or to a solution located in the vicinity of the optimal solution. It is relevant to note that, of the more than one million possible solutions that make up the search space of the problem, the Genetic Algorithm technique allowed a satisfactory answer to be found by exploring only two thousand solutions.
It is also possible to conclude, analyzing each program separately, that: -the COLUMN PROCESSING program, before being used to support the optimization programs, had its efficiency proven through a comparative study with experimental results obtained in the current literature; -the efficiency obtained by the GENETIC ALGORITHM program showed that the chosen optimization technique is not limited to a series of random crossings. Rather, the technique uses evolutionary crossings that always result in an efficient solution to the problem; -the Genetic Algorithm technique and the associated program developed in this study resulted in a robust procedure that works with discrete variables, and proved to be quite adaptable to the proposed problem; -no significant differences were found when comparing the results obtained by the two methods used in the GENETIC ALGORITHM program for constraint analysis, the Penalty Method and the Death Penalty Method. Both methods showed consistent and satisfactory results. Finally, the GENETIC ALGORITHM program combined with the COLUMN PROCESSING program proved to be an efficient computational tool, providing safe and economical solutions that guarantee the best cost-benefit ratio for the design.
Regarding the effect of the slenderness coefficient (λ) on the optimal design of reinforced concrete columns, the analysis showed that the more slender a reinforced concrete column is, the higher its cost will be. The analysis also showed that to establish the equilibrium of the structural element, the increase in bending due to second-order effects caused by the increase in slenderness starts to be absorbed predominantly by the reinforcement.