An engineering method for complex structural optimization involving both size and topology design variables

This work presents an engineering method for optimizing structures made of bars, beams, plates, or a combination of those components. Corresponding problems involve both continuous (size) and discrete (topology) variables. Using a branched multipoint approximate function, which involves such mixed variables, a series of sequential approximate problems are constructed to make the primal problem explicit. To solve the approximate problems, genetic algorithm (GA) is utilized to optimize discrete variables, and when calculating individual fitness values in GA, a second‐level approximate problem only involving retained continuous variables is built to optimize continuous variables. The solution to the second‐level approximate problem can be easily obtained with dual methods. Structural analyses are only needed before improving the branched approximate functions in the iteration cycles. The method aims at optimal design of discrete structures consisting of bars, beams, plates, or other components. Numerical examples are given to illustrate its effectiveness, including frame topology optimization, layout optimization of stiffeners modeled with beams or shells, concurrent layout optimization of beam and shell components, and an application in a microsatellite structure. Optimization results show that the number of structural analyses is dramatically decreased when compared with pure GA while even comparable to pure sizing optimization.

as many truss elements as possible. The existence or nonexistence of each element is then determined by an optimization technique. In the preliminary research works, the cross-sectional areas of the truss elements were regarded as continuous variables, and mathematical programming methods were adopted to solve these problems. 3 When the cross-sectional areas reached an extremely small value or a value of zero, corresponding truss elements were removed from the ground structure, and thus an optimal layout for these truss elements could be obtained. Not until the year of 1968 when Sved and Ginos 4 found the existence of singular optimal solutions in truss topology optimization was it recognized that those methods would lose their effectiveness in finding the singular solutions because they were treating the cross-sectional areas as continuous variables. Cheng and Guo 5 and Guo et al 6 explained the essence of the singularity phenomenon and proposed a relaxation method to handle this problem with stress and local buckling constraints. However, as stated in the work of Liu et al, 7 that method was found difficult to handle frame topology optimization with frequency constraints due to the existence of disjoint feasible domains. Therefore, it can be seen that the method is limited to specific problems, and it is hard to make it generalized in common complex optimizations. Meanwhile, by introducing integer variables 0/1 to represent the absence/presence of each truss element, the original problem that has singularity could be transformed into a new one that contains discrete variables 0/1, which is actually the main form of the mathematical model for the truss topology optimization, and it can be tackled with some optimization tools that are suitable for discrete variables. As a relatively effective way in dealing with discrete variables and evolutionary algorithms, genetic algorithms (GAs), in particular, are widely adopted in truss topology optimal designs. One major shortcoming of such algorithms is the dramatically large number of required structural function evaluations, for it is hard to establish the approximate problem of the truss topology design, which contains discrete variables 0/1, and commonly the approximate concept can only be applied to continuous optimization problems rather than problems containing integer variables. Therefore, those algorithms are not acceptable for solving practical large-scale problems because of the large number of structural analyses. 8 Thus, it can be seen that, in terms of efficiency and availability, those optimization methods are mostly effective in dealing with small-scale truss topology designs or limited to specific problems but not applicable for large-scale and complex engineering problems that aim at simultaneous optimization of trusses, beams, and/or shell components in a structural system.
To avoid using discrete variables in truss topology optimization, the concept of continuum topology optimization was put forward. 9 This technology is to seek optimal material distribution in a design space based on load and boundary conditions, resulting in an innovative design. Compared with the discrete structural topology optimization, a more efficient structure is expected to be achieved as more design space is explored, and in the meantime, traditional 0/1 discrete variables are replaced with continuous variables such as material densities, resulting in the improvement of computational efficiency. The most popularly adopted method is the SIMP method, ie, the solid isotropic material with penalization method, 10 in which the material densities for all FEs are treated as continuous design variables. Besides the SIMP model, there are also some other common methods in continuum topology optimization, like homogenization method, 11 level set methods, 12,13 and evolutionary structural optimization (ESO)/bidirectional ESO (BESO) method. [14][15][16] Techniques to solve these optimization models include optimality criterion, mathematical programming, and heuristic algorithms like GA. In most cases, as the obtained result is always in a skeletal configuration, it can theoretically be regarded as the result or a reference for discrete structure topology optimization. However, in density-based methods, like homogenization method and the SIMP method, sometimes we cannot obtain a clear configuration as many elements with intermediate densities exist. By combining additive manufacturing technologies with the continuum topology optimization, 17 these elements can possibly be fabricated. However, the fabrication of the complex design is currently infeasible, and it is believed that large-scale additive manufacturing will be a possibility in the future. 18 Additionally, aiming to produce intermediate densities as few as possible, the penalty approach is adopted to remove many of those elements. With the penalty approach, the boundaries of the newly obtained result are always coarse to treat it as real structure, and therefore, postprocessing techniques are necessarily employed to make it have clear and smooth boundaries. When using level set-based topology optimization method, one can obtain clear boundaries without any postprocessing. However, it is generally believed that the applicability of topology optimization is limited to the design of single type of components or simple structures. 18 When using high-performance computing, 18 large-scale designing problems can be handled, but, to the best knowledge of the authors' knowledge, it has not yet been found that those methods are applicable for simultaneous optimization of trusses, beams, and/or shell components in a structural system, ie, simultaneous optimization of different types of structural components.
Besides the truss structures, the frame structure is another typical form of a discrete structure. Compared with truss structures, frame structures modeled with beam elements consider both axial and flexural behaviors, which makes the topology optimization problems more difficult and challenging. Since flexural behaviors are considered, frequency-constrained frame topology design problems become singular, similar to the problem encountered in stress or local buckling constrained truss topology optimizations. Even though the relaxed method proposed by Cheng and Guo 5 and Guo et al 6 is effective in dealing with truss topology optimizations subject to stress constraints, it could be difficult to address problems with disjoint feasible domains in frequency constrained problems. 7 Therefore, some special techniques were proposed to tackle frame topology optimizations under frequency constraints. In the work of Ni et al, 19 some particular forms of area/moment of inertia-density interpolation schemes were proposed to restore the connectedness of the feasible domain; in the work of Yamada and Kanno, 20 the frequency constraint was formulated as a positive semidefinite constraint of a certain symmetric matrix, and then the constraint was relaxed to make the feasible set connected. Besides truss and frame structures, practical complex structures may also contain shells, plates, or a combination of some or all of them. For topology (or layout) optimization of complex structures, which can be a problem involving both size and topology design variables, multiple types of structural components (such as beams and plates) may need to be considered simultaneously. Although the methods mentioned in the preceding paragraphs are applicable in the optimization design of single type of structural component, they would lose their effectiveness or efficiency in handling complex structures by simultaneously optimizing multiple types of structural components, not to mention large scale problems. In fact, due to the complexity of the problem and the poor efficiency of the solving process, only a few research works on topology optimization for frames were conducted, 21 and reports on more general structures were rarely reported.
Initially proposed to address truss topology optimizations, the authors presented a two-level multipoint approximate method in previous studies. 22,23 In the approach, a series of special approximate problems related to the primal optimization problem were established by using a special approximate function, which involved both discrete (corresponding to the existence of each truss element) and continuous (corresponding to the cross-sectional area) variables. To solve these approximate problems, a genetic algorithm was adopted to optimize the discrete variables, and when calculating individual fitness values in GA, the second-level approximate problems to be solved with dual methods were built by involving retained continuous variables. As a result, structural analyses were dramatically reduced and only needed before improving the proposed approximate functions in the iteration cycles. As no intermediate variables were introduced in the process, the optimization results could be directly used without further postprocessing. Furthermore, the basic idea of this method has been extended and developed for stacking sequence optimization of composite laminates, 24-28 actuator placement optimization design in adaptive trusses, 29 and frame topology optimization, 30 and numerical examples have shown the effectiveness and efficiency of this method.
The main objective of this work is to present a method for simultaneous topology and size optimization of complex structure that may contain bars, beams, shells, and their combinations, which actually is the extension of the previous work, 30 to form a general method for common engineering problems. The two-level multipoint approximate method initially proposed for truss topology optimization design is further developed, and a generalized form of the optimization method for tackling such optimizations with continuous (size) and discrete (topology) variables, or mixed variables, is presented. The core idea is that a type of branched explicit functions is constructed to approximate the primal functions that are implicit and involve mixed variables. A genetic algorithm and the dual method are introduced to optimize discrete variables and continuous variables, respectively, in the approximate problem. The primal problem is then solved after multiple iterations. With numerical examples, this method is verified in the integrated sizing and topology optimization of frames, layout optimization of stiffeners for structures modeled with beams, shells, or a combination of them. It can be found that the number of structural analyses required in the optimization process is even comparable to size-variable optimization. By applying this method in a microsatellite structure design, its applicability in practical engineering problems is then illustrated. The structure of this paper is organized as follows. In Section 2, the paper starts by giving the generalized form of problem formulations for structural optimizations involving mixed design variables. In Section 3, the process for the two-level multipoint approximate method is elaborated, followed by some typical examples presented in Section 4. Some concluding remarks are drawn in Section 5.

PROBLEM FORMULATION
The concept of the ground structure is adopted to solve complex structure optimization involving continuous/size and discrete/topology variables. The optimization starts from an initial design with the topology of ground structure, which includes all possible components. The final optimal design only keeps some of the necessary components in the ground structure, whereas the unnecessary components are removed. To keep the meshes of the FE model unchanged during the optimization process, the unnecessary components are modeled with thin elements, which are assigned with very small values of cross-sectional dimensions. These thin elements in the FE model are called unnecessary elements that have almost no load-bearing capabilities, whereas the other normal elements modeling the necessary components are called necessary elements. Although the final design consists of both necessary and unnecessary elements in the FE model, it is expected that the physical behavior of the design is controlled solely by the necessary elements, and this means that the mechanical properties of a structure consisting of only necessary elements are almost the same with those of a structure containing both necessary and unnecessary elements, indicating that the existence or absence of the unnecessary elements will not cause much influence on the structural behaviors. This assumption is always accepted in the topology optimization of discrete structures of trusses or continuum structures. However, for frequency-constrained topology optimization problems, local modes associated with the unnecessary elements may occasionally appear. 31 One way to solve this problem is to keep the local modes from being observed during the optimization process, 32 and a similar approach was used in the works of Tcherniak 33 and Takezawa et al 34 by letting the mass value of the unnecessary elements equal zero. Likewise, in the present work, the material density values of the unnecessary elements are assigned with an extremely small value, whereas the densities for the retained elements keep their original value.
Based on the idea described earlier, structure optimization involving continuous/size and discrete/topology variables can be formulated as follows: where X is the continuous/size variable vector, x k is the element in that vector, and m is the number of such continuous variables; x L k and x U k are the lower and upper bounds on the continuous/size variable x k , respectively, if the related structural component is retained; k is a small value to represent the property of the structural component when it is not necessary; is the discrete/topology design variable vector and i denotes the presence (1) or absence (0) for the associated component; n is the total number of structural components whose presences are to be decided; f (X, ) and g j (X, ) are the objective function and the jth constraint function, respectively, and J 0 is the number of considered constraints; X I k and X k (k = 1, … , m) are two domains that are defined according to the values of the continuous/size variables, and they are formulated as follows: When i = 1, the ith structural component is necessary and its related continuous/size variables x k are changed within given bounds, and the values of these continuous variables belong to X I k . Likewise, if i = 0, the ith structural component is to be deleted and the associated continuous/size variables are assigned with small values, and these values belong to X k . Additionally, in cases where some of the structural components are constantly present in the structure and only their cross-sectional dimensions are to be optimally determined as continuous variables, for such components, their corresponding size variables always fall into the domain of X I k . To illustrate the formulation established earlier, let us consider the optimization design of stiffener layout for plate/shell structures as an example (as shown in Figure 1). The stiffened cylinder made of metallic material is fixed at the left end as the boundary condition. Eight stiffeners with solid-rod cross-sections are uniformly distributed along the longitudinal direction, which is considered as the initial design. For this problem, the existence of each stiffener can be optimally determined, and the radii or areas for the retained stiffeners need to be optimized simultaneously. The objective is to minimize the structural mass, and the constraints could be the requirement that the fundamental frequency should be more than a given value.
Based on the problem definition in Equation (1), the design variables for stiffener layout designs can be specified as X is the continuous size variable vector, consisting of cross-sectional dimension variables x i , and x i is the l i th cross-sectional dimension variable for the ith stiffener in the initial structure. L i is the number of cross-sectional dimensions for the ith structural component group (there are more than one component in a group when design variable FIGURE 1 A stiffened cylinder link is introduced). is the discrete design variable vector and i denotes the presence (1) or absence (0) for the ith stiffener. n is the total number of stiffeners to be decided in the initial structure, and it is also the number of discrete variables. m is the number of continuous variables and m = ∑ n i=1 L i . According to the descriptions earlier, for the design problem in Figure 1, there are eight continuous size variables, ie, X = {x 1 , x 2 , … , x 8 } T , indicating the cross-sectional dimensions (ie, radii or areas) of the eight stiffeners (where L i = 1, i = 1, … , 8), and eight discrete variables, ie, = { 1 , 2 , … , 8 } T , which take the value of 1 or 0 to represent the presence or absence of the eight stiffeners. In a certain case, besides the layout of the stiffeners, the shell thickness of the cylinder is also treated as a continuous design variable, but the cylinder shell is supposed not to be removed. In this case, the continuous size variables should be updated as X = {x 1 , x 2 , … , x 8 , x 9 } T , where x 9 is the cylinder thickness variable and it belongs to the domain X I k . Meanwhile, the discrete variables should still be defined as { 1 , 2 , … , 8 } T because the cylinder cannot be removed.

The first-level approximate problem
To solve the problem in Equation (1), which is always implicit and nonlinear, a first-level approximate problem is firstly constructed to transform it into an explicit problem. In the pth stage, the first-level approximate problem is where f (p) (X, ) and g ( ) (X, ) are the approximate objective and constraint functions, respectively, created for the pth stage; J 1 is the number of active constraints related to the primal problem in Equation (1); andx U k(p) andx L k(p) are the move limits of x k at the pth stage.
The multipoint approximation method has proved to be one of the most effective approximation techniques for complicated structure optimization problems, 35 and a branched multipoint approximation (BMA) function is employed here to approximate the primal problem. For a unified formulation, the approximate functions f (p) (X, ) and g ( ) (X, ) in Equation (4) are substituted with the BMA function w (p) (X, ), as shown in the following: wherew and X t or t is the tth known point, h t (X, ) is the weighted function, and H is the number of points to be counted with the upper bound H max . When the number of known points is larger than H max , only the last H max points obtained are taken into account (in this work, H max = 5). It can be seen that the construction of the approximate function only requires the partial derivatives to the continuous variables, and discrete variables only act as the role of branching, deciding x k to belong to which domain, ie, X I k or X k . The form of the BMA function is similar to the Hermite interpolation function, and in the known points, the value of the approximate function as well as their derivatives are consistent with the values in the primal function. The approximation accuracy of BMA functions have been discussed in the works of Huang et al. 36,37 From Equation (7), the BMA function w (p) (X, ) can actually be regarded as the weighted representation of approximate functionsw t (X, )(t = 1, … , H). In fact, the weighted function h t (X, ) satisfies the following conditions: 0<h t (X, )<1 and ∑ H t=1 h t (X, ) = 1.w t (X, ) is similar to the Taylor expansion with respect to the continuous intermediate vari- } T , and meanwhile, with two branches of expressions, the approximation function is still available at the separated design point, which has extremely small values.
In Equation (9), the exponents r o,t and r m,t (t = 1, … ,H) are the adaptive parameters used to control the nonlinearity of the approximation. When t = 1, … ,H-1, r o,t and r m,t can be obtained by using the least-squares parameter estimation, as shown in the following: Since the first-level approximate problem involves both continuous and discrete variables, it can never be solved directly by using general mathematical programming methods. If GA is directly used and the continuous and discrete variables are encoded simultaneously, the scale of design variables and the computational costs will become tremendous. Thus, a layered optimization strategy is introduced. Discrete variables are optimized with GA in the external layer and continuous variables are optimized in the internal layer by using the dual method, which could significantly reduce the gene code length in GA while improving the optimization efficiency as well as accuracy.

GA to optimize discrete variables
Genetic algorithm imitates the evolutionary process of living organisms and starts with a gene coding operator that encodes an individual, ie, a design point in the search space, as a gene string. 7 The discrete variable described in Equation (1) { 1 , 2 , … , n } T takes 1 or 0 to represent the existence or nonexistence of each component in the initial design. Accordingly, one string, ie, str = 1 2 · · · n , is used as gene code in GA. After the random generation of the initial population, the GA procedure is conducted, following the fitness calculation for each member in the population. We assume that the number of individuals in the population is N, and the maximum generation number is given as MaxG.
Based on the adaptive scheme, 25 a constrained minimization problem using an exterior penalty function is established for the fitness calculation of the Sth individual where In addition, ⟨ f (p) (X, )⟩ represents the average value of the objective function over the current population, g (p) is the violation of the jth constraint averaged over the current population, f (p) (X S , S ) and g (p) (X S , S ) are the objective value and constraint values with respect to the optimal continuous values X S under a given configuration S , and all of these values will be obtained by solving the second-level approximate problem, which will be detailed in the following subsection. is a small fraction of the maximum constraint to be subtracted from the objective to discriminate between multiple designs with the same objective value, all of which satisfy all the constraints. When the fitness calculation is completed, a simulated roulette wheel selection method is used as the reproduction operator, followed by the crossover and mutation operations. A one-point crossover method is adopted, and a simple mutation method is used as the mutation operator.
The GA is executed sequentially by using the operations earlier. It does not stop until the number of evolution generations reaches the maximum generation number, ie, MaxG, which is preset.

The second-level approximate problem to optimize continuous variables
For each member in GA, the discrete variables in the first-level approximate problem are fixed as s . The continuous variables associated with the structural components whose corresponding discrete variables are zero will not change and they are fixed with small values. Therefore, only such continuous variables that are related to the necessary components are left in the first-level approximate problem. Considering that the number of retained continuous variables is usually more than the number of active constraints, it is reasonable to use the dual method to solve this problem. However, the first-level approximate problem is still in a complex form, and we cannot easily get a clear relationship between design variables and dual variables. Thus, a second-level approximate problem that can be solved by the dual method is established to approach the optimum of the first-level approximate problem. For a more common use, the objective and constraint functions in the first-level approximate problem are expanded into the linear Taylor series with respect to the continuous variables and their reciprocal variables, respectively, to construct the second-level approximate problem in a so-called variable separable form, only with which form the dual method can be effectively used. Under a given discrete variable vector s , this approximate problem in the qth step is formulated as follows: whereX is the continuous variable vector related to the necessary components, and I is the number of those retained continuous variables;̃( q) (X , s ) is the objective function at the qth step;g (q) (X , s ) is the constraint function, and J 2 is the number of the retained constraints;x U k(q) andx L k(q) are the move limits ofx k at the qth step; and x U k(q) and x L k(q) are the upper and lower bounds onx k . Then, the established second-level approximate problem earlier can be easily solved with dual methods.
The dual problem for the approximate problem Equation (15) can be stated as follows: where The relationship between the continuous design variables and the dual variables can be established as follows: where

Flow chart and convergence criteria
The flow chart of the optimization method described earlier is shown schematically in Figure 2. It should be noted that the adopted multipoint approximation function is gradually established well with the known design points that are obtained from the iteration process, and the number of required structural analyses is equal to the iteration number. Therefore, it is quite different from the common response surface function, which needs to be built by conducting structural analyses among a series of testing points before starting the optimization. Moreover, as the dramatic function evaluations required in GA are achieved by using the established approximation function instead of direct structural analysis, the computational costs could be reduced significantly. In addition, some methods such as hard-kill ESO/BESO methods would cause some elements to be deleted and never make them reversed, resulting in an extremely nonoptimal design, 38

FIGURE 2
Schematic overview of the optimization strategy and with a refined mesh, an optimal design is still possibly obtained. 39,40 Unlike those methods, the method proposed in this work makes it possible that the removed elements in the optimization process can be regained. That is because the sensitivities on those elements might be very large when they are represented by a small value, and the large values of sensitivities could be expressed in the first-level approximate problem, which is built based on the actual response values and sensitivities so that the optimization process would make the removed elements re-emerge. Within the optimization calculation procedure, the convergence criteria for both the first-level and second-level approximation problems are applied as follows.
1. The first-level approximation problem convergence evaluation.
There are three conditions involved. or where DELTA1, DELTAX, and DELTAC are convergence control parameters, and LPMAX1 is the maximum iteration number of the first-level approximation problem. To evaluate the convergence, conditions a, b, and c should be satisfied simultaneously.

The second-level approximation problem convergence evaluation.
Two conditions are involved.
where DELTA and NMAX1 are the convergence precision and maximum iteration number of the second-level approximation problem. When either of the two conditions is satisfied, the second-level approximation problem is converged. In this work, DELTA1 = 0.001, DELTAX = 0.001, DELTAC = 0.01, LPMAX1 = 50, DELTA = 0.001, and NMAX1 = 10.

NUMERICAL EXAMPLES
Numerical examples are given in this section to illustrate the effectiveness and efficiency of the proposed method in handling problems with discrete/topology design variables. The first example is to optimize the topology of a frame structure, and the computational cost is compared with pure GA so that the efficiency of the proposed method is highlighted. In the second and third examples, the stiffener layouts for a plate or a shell structure are optimally designed, and the stiffeners in the second example are modeled with beam elements, whereas they are modeled with shell elements in the third example. Example 4 is to concurrently optimize the layouts for both beam and shell components in a single structure, and Example 5 is an engineering application of a microsatellite, where the layouts of stiffeners, which are modeled with beam elements, and the thicknesses of plates that are not considered to be removed from the structure system are simultaneously optimized. Both Examples 4 and 5 are to show the effectiveness of the method in dealing with multiple types of structural components. Meanwhile, with the same design parameters, optimization results from pure sizing optimization are also given in Examples 1 to 4. The optimization method used for pure sizing optimization is the two-level multipoint approximation approach, 35 where the primal problem is approximated and made explicit with the first branch of Equation (9). This approach has been applied in benchmark studies with high accuracy and efficiency, and also been successfully used for practical and complex engineering problems. 41

Example 1: topology optimization for a frame structure
A five-beam frame structure is shown in Figure 3, with a mass value of 454 kg for each lumped mass point. All beams have solid round sections and the material density is 2770 kg/m 3 . Young's modulus is 6.89×10 10 Pa, with a Poisson's ratio of 0.33. This structure is expected to be designed by optimizing its topology and corresponding cross-sectional sizes. Under the constraint that the fundamental frequency should be more than 14 Hz, the structural mass was minimized with an improved version of genetic algorithm based on the work of Cheng and Liu, 43 and the results are given in Table 1. By using MSC.Patran/Nastran 2005, 44 this structure is analyzed with the obtained results from the aforementioned work, 43 and it is found that the fundamental frequency is 15.76 Hz. In order to make a suitable comparison, in this study, the lower bound on the frequency constraint is replaced with 15.76 Hz instead of 14 Hz. The initial designs for all radii are 0.5 m, and the lower and upper bounds are 0.2 m and 2.5 m, respectively. With the proposed method, the population size and the maximum generation number are both 50 in GA, and the optimization result is given in Table 1. By considering structural symmetry, the other case of the optimum structure can be obtained by deleting beam 5 while retaining beam 4 with the optimum sizes unchanged, as demonstrated in the work of Cheng and Liu. 43 Thus, the obtained topology configuration is rational in this work. As the optimum radii are quite different from those in the aforementioned work, 43 the structural mass is quite smaller, whereas the fundamental frequency keeps unchanged. With the use of GA in the work of Cheng and Liu, 43 a population size of 100 and a maximum generation number of 200 were used to obtain the optimal solution. As structural analysis was needed for each member in the population at each generation, the number of structural analysis shown in Table 1 is calculated via multiplying the population size by the generation number. With the same initial design and lower and upper bounds on size variables, a design case of pure sizing optimization (Case A) is conducted, and its optimization results is also given in Table 1. It can be seen that the radius of beam 3 reaches its lower bound, ie, 0.2 m. By decreasing all lower bounds to 0.01 m, another case (Case B) of pure sizing optimization is carried out, and the design results are also listed in Table 1. The optimization results in Case B does not change too much compared with the results in Case A, and the structural masses obtained in pure sizing optimizations are larger than those achieved in topology and sizing optimizations. Therefore, those solutions in pure sizing optimization can be seen as local optima. The optimum topology design cannot be obtained from pure sizing optimization because frequency-constrained frame topology optimization is strongly singular in the sense that its feasible domain is disconnected and the global optimal solution is located at the tip of a separated low-dimensional subdomain. 19  In this work, by introducing discrete variables in the problem formulation and using GA to optimize discrete variable so that each possible subdomain is searched, the possibility to achieve the global optimal solution is realized, as seen in Table 1. On the other hand, the number of required structural analyses in GA is around 2000, whereas it is only 14 in the present method, which is comparable to pure sizing optimization, as shown in Table 1. The iteration history for the structural mass is given in Figure 4. Here, an iteration refers to a solution of the first-level approximation problem, which is obtained after many times of GA execution, and after each iteration, a structural analysis and sensitivity analysis would be implemented to check the convergence of the first-level approximation problem, as described in Figure 2. A complete independent run consists of several iterations of such execution process. The optimization result shown in Table 1 and Figure 4 is a complete independent optimization run, in which 14 iterations are involved. After comparisons as stated previously, it can be seen that the proposed method can approach reasonable results with remarkable efficiency.

Example 2: layout optimization for stiffeners modeled with beams
The second example is a plate that is fixed at the corners and loaded at the center by a point load F = 200 N normal to its plane, as shown in Figure 5. Both the length and width of this plate are 0.1 m, and its thickness is 1×10 −3 m.  Young's modulus for the material is 2×10 11 Pa, and Poisson's ratio is 0.3, with a material density of 7800 kg/m 3 . This problem is to find the best stiffener layout for this base plate, and the stiffeners have a rectangular cross-section, as shown in Figure 6, with the same material of the plate.
For this design problem, the objective is to minimize the structural mass, and the design constraint requires that the maximum displacement should be less than 5×10 −4 m, ie, 0.5 mm. Three different initial designs, as shown in Figure 7A, C, and E, are considered, and all the stiffeners are modeled with beam elements for structural analysis. The number of stiffeners contained in each case is summarized in Table 2. By considering symmetry, only a quarter of the plate is designed, and the rest parts are obtained by using variable linking. As each stiffener corresponds to a discrete variable, which is used to decide their existence, there are 5, 20, and 72 discrete variables for each case, respectively, under the symmetry consideration. Meanwhile, the cross-sectional dimensions for the stiffeners, ie, the height (H) and the width (W), are also optimized as continuous variables. The number of both discrete and continuous variables for each case is also listed in Table 2. The initial values for H and W are both 4 mm, and their lower and upper bounds are 1 mm and 4 mm, respectively.
With the proposed method, the optimization results for the three cases are shown in Figure 7B, D, and F as well as Table 3. For Case 1, it can be observed that this plate appears to require only two stiffeners along two diagonal lines crossing the plate. In Case 2, besides the stiffeners along the diagonal lines, another four stiffeners are required around the center. Similar case happens in Case 3, where the required stiffeners are mainly located around the corners and the center where the force is applied. As the plate is fixed at four corners, which is similar to a clamped beam loaded by a concentrated force in the middle, the resistant bending moments in the center and corners are quite large. Therefore, the plate is mainly strengthened with stiffeners around the corners and the center. Meanwhile, as for the middle parts between the center and corners, the resistant bending moments are relatively small so the stiffeners becomes unnecessary at those parts. As a result, the stiffeners are disconnected in Case 3. For the three cases, with the increase of the amount of stiffeners, their initial structural mass is getting larger, whereas the obtained structural mass is getting smaller, which is due to the fact that the design space gets larger when the number of stiffeners increases. As for the computational costs, it can be found that the number of structural analysis for all three cases is less than 30, exhibiting high efficiency, and the iteration histories for the design objective are plotted in Figure 8. Corresponding pure sizing optimization results are also given in Table 3. For each case, the structural mass obtained in simultaneous layout and sizing optimization of stiffeners is smaller than that achieved from pure sizing optimization, and the number of required structural analysis is comparable to sizing optimization.
The optimization results are then compared with those from the references. In the work of Fatemi and Trompette, 45 a similar design problem was considered. Even though the geometry parameters in that literature are a bit different, the comparison between them is also possible. In the aforementioned work, 45 only two stiffeners were considered and the coordinates of the ends for both stiffeners were treated as design variables. By using GA as well as a local optimizer, the best layout for the two stiffeners were obtained, as shown in Figure 9A. From Figure 7B, it can be seen that the proposed method can achieve the same configuration. This problem was also studied by using continuum topology optimization method in the work of Zhang et al, 46 and the corresponding result is presented in Figure 9B. The material distribution is also along the two diagonal lines, and the obtained structure is strengthened in the center by having more materials. Similar results can also be found in Figure 7D and F. Additionally, for the result in Figure 7F, the stiffeners along the diagonal lines to connect the elements in the center with those around the corners are removed, and similar case can also be seen in Figure 9B, where corresponding regions require very few materials. After those comparisons previously, it can be said that the obtained results with the proposed method are quite reasonable.

Example 3: layout optimization for stiffeners modeled with shells
A simplified mold used for composite material component forming is designed in this example. This mold is a stiffened shell structure, where the stiffeners are modeled with shell elements, as shown in Figure 10. It is fixed at the corners and loaded with a distributed load of 1×10 5 N/m 2 on the top of the shell surface, normal to its plane. The length and width of the structure are 1.05 m and 0.5 m, respectively, and the height is 0.08 m. Young's modulus for the mold is 7×10 10 Pa, and Poisson's ratio is 0.33, with a material density of 2700 kg/m 3 . This design is to find the best stiffener layout for this shell structure. An initial design for the stiffener layout is shown in Figure 10, and the stiffeners are uniformly placed along the length and the width directions. The structural mass is to be minimized, and the design requirement is that the maximum displacement should be less than 5×10 −3 m. In the initial design, there are 58 stiffeners. By considering symmetry, only a quarter of the plate is designed, as shown in Figure 10B, and the remaining three-quarters are then obtained by using variable linking. Since a stiffener corresponds to a discrete variable to decide the existence, there are 17 discrete variables with the symmetry consideration. In the meantime, the shell thicknesses for each stiffener are also considered as continuous variables. Therefore, there are also 17 continuous size variables involved. The initial shell thickness for each stiffener is 5 mm, and its lower and upper bounds are given as 3 mm and 5 mm, respectively.
The optimization result is given in Figure 11 as well as in Table 4. After optimization, only 16 stiffeners are retained, whereas the considered constraint is satisfied. From the thicknesses that are given in Table 4, it can be found that the thickness values around the corners and the center are relatively larger than those of the other parts, which means the regions near the corners and the center should be strengthened. This is very similar to the results in the aforementioned example, ie, Example 2, where the optimal layout requires stiffeners around the corners and the center. Thus, it has shown that this result is also reasonable. The number of required structural analyses is only 11, and the iteration history for the objective is plotted in Figure 12. Pure sizing optimization results are also listed in Table 4. All size variables get to their lower bounds, but the structural mass is still larger than that from integrated layout and thickness optimization for stiffeners.

Example 4: concurrent layout optimization of beam and shell components
The structure considered in Case 2 in Example 2 is used in this example for concurrent layout optimization of beam and shell components in a single structure. Therefore, besides the beam stiffeners, the base plate is also considered to be designed. This plate is divided into 16 regions, as shown in Figure 13, and each region is possibly to be removed from the structure. Similarly, only a quarter of the plate is designed, and by considering symmetry, the rest parts are obtained by using variable linking. The structure geometry and the loading conditions, as shown in Figure 14, are the same with that in Example 2, and the objective as well as the design requirement is also the same.
For this problem, in addition to the 20 discrete variables that correspond to the existence of each stiffener, another four discrete variables, which are used to decide the existence of each region of the plate, are also involved. As for the continuous variables, there are 40 size variables for the beam cross-sectional dimensions (ie, H and W), and four thickness variables for the plate. The initial values and the lower and upper bounds on the beam cross-sectional dimensions are the same with those in Example 2. The initial value for the plate thickness is 1 mm, and the lower and upper bounds are given as 1 mm and 4 mm, respectively.
With the proposed method, the optimization results are shown in Figure 15 and Table 5. It is found that all regions of the plate are removed, and 14 stiffeners are retained in the optimization result. Compared with the result in Example 2, with the deletion of the base plate, Stiffeners number 14 instead of Stiffeners number 18 are required to form a closed-up structure for load transmission. Even though the number of required stiffeners is 14, ie, larger than 12 in Example 2, the obtained structural mass is 0.02470 kg, ie, quite smaller than 0.09630 kg, mainly caused by the removal of the base plate. The number of structural analyses is 26, and the objective iteration history is shown in Figure 16. Corresponding sizing optimization results are also shown in Table 5. Even though the maximum displacements in both cases are around the design boundary, the structural mass obtained in pure sizing optimization is much larger than that from concurrent layout optimization of beam and shell components.

Example 5: an application in a microsatellite optimization of plate thicknesses and the layout for stiffeners modeled with beams
This example deals with an engineering application, which is to optimally design a microsatellite. The FE model for the main body of this satellite is shown in Figure 17. Considering the connecting interface between the whole satellite structure and the launch vehicle, the bottom of the main body is fixed as the boundary condition, as shown in Figure 17. The satellite is divided into seven sections by eight clapboards, and these plates are made of aluminum with beam stiffeners. All these stiffeners have solid rectangular cross-sections, and the initial layout of these stiffeners on the plates is shown in Figure 18. The design problem is to optimize the plate thickness for the eight clapboards and the stiffeners layout on these plates. The initial structures are designed with empirical experiences, and the total mass of the whole satellite is 18.622 kg. The mass for the design domain, which consists of eight plates and their beam stiffeners, is around 1.183 kg, and mass values for each plate and its stiffeners are in the range of 0.1∼0.33 kg. As Plate_1 has the largest number of stiffeners in the initial design, the mass value for this plate and its stiffeners is the biggest among all the eight plates, which is 0.3293 kg. There are many holes in Plate_7 and the number of beam stiffeners is relatively small, so Plate_7 has the smallest mass value, which is 0.1078 kg. For this application problem, the structural mass is to be minimized, and the design requirements are that both the first and the second order of transverse bending vibration should be more than 25 Hz. Besides, the first order of the longitudinal vibration should be more than 50 Hz. The initial values for these vibration modes are 26.17 Hz, 26.32 Hz, and 58.44 Hz,

FIGURE 12
Objective iteration history for Example 3 respectively, and it can be seen that this initial design is a feasible one. Since each stiffener corresponds to a discrete variable and there are totally 153 stiffeners for the eight plates, there are 153 discrete variables involved in the problem to decide the existence of each stiffener. As the cross-sectional dimensions (ie, height, H, and width, W) of the stiffeners are also considered, there are 306 (153*2) continuous variables associated with the stiffener cross-sectional dimensions.
(A) (B) Moreover, the thicknesses of the eight plates are also to be optimized as continuous variables. According to the problem formulations established in Section 2, eight continuous variables corresponding to plate thicknesses are also included. With the proposed method, the optimization results for the stiffener layout are shown in Figure 18. Together with their initial values and the lower and upper bounds, the optimization results on the continuous variables are given in the Supporting Information. After the optimization, the structural mass of the design domain is 0.808 kg, ie, a reduction of 0.375 kg with respect to the empirical design. If represented by percentage, this mass reduction is then 31.69%. The first and the second order of transverse bending vibrations are 28.64 Hz and 28.95 Hz, respectively, and the first order of the longitudinal vibration is 51.52 Hz, which all satisfy the design requirements.
Therefore, from the results earlier, it has been verified that this method is effective in dealing with problems where multiple types of design variables are involved, and this method has the capability to handle practical engineering problems so that the results can be provided to the engineers as a design reference.

Summary of numerical examples
From the results of the numerical applications, it can be found that the proposed method is effective in addressing problems for complex structures with discrete/topology design variables, including the frame topology optimization, optimal layout design for beam/shell stiffeners, and concurrent optimization of layouts for both beam and shell components. Pure sizing optimization results are also given in Examples 1 to 4, and the obtained structural masses are always larger than those achieved in design cases where topology (or layout) variables are also considered. For the examples of stiffener optimization, the cross-sectional sizes of the obtained stiffeners are not tapered, so these edges pose regions of stress concentration for practical applications. This is not the focus of this paper, and it will be considered in our future work. In sizing optimizations, some design variables get to their lower bounds. If the lower bounds are decreased, the size variables cannot yet reach zero or a very small value, as shown in Example 1. That is caused by the characteristics of singular optimum and disjoint feasible domains, which are always encountered in topology design problems. On the other hand, from the numerical examples previously, it can be observed that the number of structural analysis is even comparable to pure sizing optimization. Thence, the efficiency of this method is noticeable, especially when compared with pure GA.
In addition, for topology optimization designs, it might end up with a structure that contains too many thin members. Here, the "thin members" refers to the structural members whose relevant continuous variables, like cross-sectional dimensions, stay at the given lower bounds or so. This might happen, especially for large-scale problem with many variables, but it should be noted that this phenomenon is quite different from the topology optimization designs where the material densities or cross-sectional areas are treated as continuous variables. In those topology optimization problems, the result may contain many thin members, whose material densities for instance stay at intermediate values, which makes it difficult to decide whether to delete them from the structure or not. Under the formulation that is presented in this work, in order to prevent that result happening, which might contain many thin members, it is possible to add a constraint regarding the number of nonzeros, and similar constraints have been considered in the actuator optimal placement in adaptive trusses 29 to limit the number of actuators. Further investigations on this issue will be the focus in future work.

CONCLUSIONS
An optimization method that is suitable for tackling problems with mixed variables is presented in this paper. A generalized form of problem formulation for such optimizations with continuous/size and discrete/topology variables, or mixed variables, is firstly given and then an optimization process to solve this kind of problem is elaborated. The primal problem, which is implicit and involves discrete and continuous variables, is made explicit by using a type of branched multipoint functions, so the approximate problem is constructed, which also involves mixed variables. To solve this approximate problem, a genetic algorithm and the dual method are introduced to optimize discrete variables and continuous variables, respectively. The primal problem is then solved after multiple iterations. By applying this method in structural optimization problems, like frame topology optimization and layout optimization for stiffeners modeled with beam or shell elements, it can be found that the number of structure analysis required in the optimization process is even comparable to pure sizing optimization, and compared with a pure genetic algorithm, the computational costs can be dramatically saved. When using this method to solve the concurrent layout optimization of beam and shell components in a structure and to address a microsatellite design problem, this method has shown its capability in dealing with multiple types of design variables and handling practical engineering designs. Thence, the effectiveness and efficiency of the presented approach is verified with those numerical examples, and it will be expected to be used for other complex structural designs that involve discrete/topology variables.