The purpose of this paper is to improve the efficiency of multiobjective topology optimization using a genetic algorithm (GA) with barsystem representation. We proposed a new GA using an elite initial population obtained from a Solid Isotropic Material with Penalization (SIMP) using a weighted sum method. SIMP with a weighted sum method is one of the most established methods using sensitivity analysis. Although the implementation of the SIMP method is straightforward and computationally effective, it may be difficult to find a complete Paretooptimal set in a multiobjective optimization problem. In this study, to build a more convergent and diverse global Paretooptimal set and reduce the GA computational cost, some individuals, with similar topology to the local optimum solution obtained from the SIMP using the weighted sum method, were introduced for the initial population of the GA. The proposed method was applied to a structural topology optimization example and the results of the proposed method were compared with those of the traditional method using standard random initialization for the initial population of the GA.
1. Introduction
Topology optimization is one of the useful methods for finding the layout of an optimal structure in a given design domain for the loadbearing structures. Generally, because there are usually two or more loadings on a structure in a structural topology optimization problem, it becomes a multiobjective optimization problem
[1]
. In multiobjective optimization problems, solutions consist of a Pareto optimal set, which can only be improved by degrading at least one of its other objectives. The objective in solving multiobjective topology optimization problems is to find complete Paretooptima solutions
[2]
.
For the topology optimization, various methods have been developed
[2

8]
. One popular approach for multiobjective topology optimization is the solid isotropic material with the penalization (SIMP) method (Bendsøe,1989) using the weighted sum method
[4]
. The SIMP method redistributes the material throughout the specified design domain using an optimality criterion
[3]
. The SIMP method is extended to multiobjective optimization using the weighted sum method, because this method is relatively easy to implement and finds Paretooptimal solutions quickly. However, the SIMP method has multiple optima because most topology design problems have nonconvex objective functions
[4]
. Thus, it is difficult to find a single optimal solution. Moreover, the weighted sum method cannot find a complete Paretooptimal solution set if the objective space is nonconvex
[2]
.
To perform global searching and find the complete Paretooptimal set, genetic algorithms (GAs) have been increasingly used in the multiobjective topology optimization problems. In the GAs, individuals of the GA are evolved simultaneously by cooperation and competition. Thus, a GA can find Paretooptimal solutions in a single simulation. However, because a GA requires a large number of evaluations of objective functions, its computational cost becomes high. This computational cost problem is the most important factor in optimization problems using GAs
[8]
.
In this study, for an efficient multiobjective topology optimization method using a GA, the elite initial population method was proposed. The elite initial population which is similar to local optimum solutions obtained from the SIMP with the weighted sum method was used for the initial population of the GA. The SIMP with the weighted sum method can provide well converged and distributed individuals in the objective function space. A local optimum solution used as a good initial guess helps to improve the efficiency of the GA. In this study, to verify the convergence and diversity of Paretooptimal solutions of the proposed method, a standard initial population and an elite initial population were compared.
2. Elite initial population
 2.1 SIMP with the weighted sum method for a multiload case problem
In this study, local optimum solutions, which are obtained by SIMP with the weighted sum method, were used for an initial population of GA. In the SIMP with the weighted sum method, the problem of optimizing the structure topology for minimum compliance can be formulated in a discrete form as follows.
Comparison of the conceptual diagrams of a traditional GA method and the proposed method
where C is the compliance,
x_{e}
is the e element density,
k_{e}
is the
e
element stiffness matrix,
u_{i,e}
is the
e
element displacement for the i loadcase, V is the maximum volume, K is the global stiffness matrix, U
_{i}
is the
i
loadcase displacement vector and F
_{i}
is the i loadcase vector. And
p
denotes the penalty parameter,
n
is the total number of loadcases, and
m
is the total number of finite elements in the design domain.
In the SIMP using the weighted sum method, Paretooptimal solutions can be easily obtained by using a varied set of weights such that
and w
_{i}
> 0. The computational efficiency of the SIMP method is greater than a GAbased method, because it is a gradientbased method. However, if the objective function space is multimodal (nonconvex), the optimization problem has nonunique solutions when choosing different initial points and design parameters for the algorithm
[4]
. Moreover, although implementation of the weighted sum method is very easy, it cannot find complete Paretooptimal solutions
[2]
. Conversely, the GA topology optimization can find complete Paretooptimal solutions more effectively. However, it is necessary to reduce the computational cost
[8]
.
Therefore, in this study, local optimum solutions obtained from the SIMP with the weighted sum method, which is more computationally effective than simply using a GA, were employed for the initial population of the topology optimization to improve the efficiency when using a GA. Well converged and distributed local optimum solutions become good initial guesses for the GA, and lead to more rapid convergence to Paretooptimal solutions as shown in
Fig. 1
.
 2.2 Elite initial population
Figure 1
shows illustrations of the traditional and proposed methods. Generally, in the traditional method for GAs, a standard random initializing procedure is used for the initial population. Its initial population is completely random as shown in
Fig. 1
(a). Conversely, the initial population of the proposed method already has well converged and distributed solutions, because local optimum solutions were used as initial populations as shown in
Fig. 1
(b). By comparing
Figs. 1
(a) and
(b)
, the proposed method is expected to converge more rapidly than the traditional method.
The flowchart of the proposed method is presented in
Fig. 2
and summarized here:
1) Generate
M
local optimum solutions using the SIMP with the weighted sum method.
2) Create
M
×
N
individuals that have topologies similar to each
M
local optimum solution of Step 1 using topology optimization using the GA.
3) Randomly select L individuals obtained from Step 2 for the initial population of topology optimization using the GA.
In topology optimization using the GA in Step 2 and 4, the barsystem representation method
[8]
was used for the topology representation. It is well known that the choice of a topology representation method is critical to the process. In the barsystem representation, the positions of the vertices and thickness of the bars are used as design parameters in the topology optimization using a GA. Mapping on the
Flowchart of the proposed method and comparison of the proposed and traditional methods.
FE mesh is simple in the barsystem representation. If the center of gravity of that element lies inside bars, the density of the element becomes 1. Thus, the density of an element is only 0 or 1.
Figure 3
shows the topology example for a barsystem representation.
Figure 3
(a) and
(b)
shows an example of the genotype (artificial chromosome) and phenotype (the decoded parameter), respectively.
However, it is difficult to use local optimum solutions directly for the initial population of a GA, because properly decoded GA genotype must be found, corresponding to local optimum solutions with barsystem representation. Thus, topology optimization using a GA was used in the Step 2 to create individuals corresponding to local optimum solutions. We called these individuals the “elite initial individuals”. However, it is difficult to find the elite initial individual with exactly the same phenotype as the local optimum solution. Moreover, there are many genotypes corresponding to one phenotype. Thus,
N
elite initial individuals that have a similar phenotype but different genotypes for each local optimum solution were found. We present the procedure to create the elite initial population as follows.
First, in the FEM element of the local optimal solutions of Step 1, if an element has a density below 0.05 it will be modified to 0.5. The densities of other elements remain as they are. However, if a volume constraint is introduced, it is not necessary to modify the density of any elements in this step. The density values 0.05 and 0.5 were determined empirically.
Second, to find the proper genotypes corresponding to the local optimal solutions, the density distribution of the modified local optimal solutions, and the phenotype transformed from the genotype using the barsystem representation, were compared in Step 2 topology optimization using a GA with formula (6) as the objective function.
Genotype and phenotype in the barsystem representation
where
m
is the total number of finite elements in the design domain,
x_{s,i}
is the ith element density of the density modified local optimal solution, and
x_{b,i}
is the
i
element density of the individual in the GA. By using formula (6), the most correct genotype, which corresponds to the local optimal solution, can be obtained.
Formula (6) is used as the objective function of topology optimization using the GA of Step 2. However, the computational cost must be investigated. Because formula (6) is a simple sum of the products of a matrix, its computational cost will be low. We present an example of formula (6) as shown in
Fig. 4
. Because
m
(total number of elements) is 9, the sum of the product of the reference individual (local optimal solution) and density matrix of the individual of the bar system is 3.5, the result of formula (6) will be 5.5.
3. Multiobjective optimization problem method using GA
 3.1 Multiobjective optimization problem
Generally, a multiobjective optimization problem can be formulated in a discrete form as follows:
Calculation example for the proposed objective function
where f
_{m}
is the objective function, g
_{j}
is the inequality constraint, h
_{k}
is the equality constraint and
x_{i}
is the design variable.
 3.2 Elitist nondominated sorting GA with archiving
In the multiload case problem, the elitist nondominated sorting genetic algorithm (NSGAIIa)
[9]
which is one of the best known MultiObjective Evolutionary Algorithms (MOEAs), was used to minimize the compliance of each loadcase.
The flow of NSGAII is similar to the flow of a common GA
[10]
as shown in
Fig. 5
. However, to improve diversity and convergence of solutions in this algorithm, the crowding distance metric and the crowded tournament selection are introduced. The crowding distance metric is the relative average distance to other individuals with the same rank in the objective space. The crowded tournament selection is used to compare two individuals and returns the winner of the tournament according to the rank and crowding distance between them. Moreover, an elitist strategy was applied to avoid the problem of losing Paretooptimal solutions. If the number of Paretosolutions exceeds the population size, some Paretooptimal solutions can be eliminated. Therefore, an archiving method was used for the solutions preserving method. For the crossover method of a GA, the parent centric crossover (PCX)
[11]
and the binarylike crossover
[10]
methods were applied simultaneously, and for the mutation method of a GA, the polynomial mutation
[12]
was applied.
Flowchart of the genetic algorithm
4. Verification of the proposed method
To verify the efficiency of the proposed method, a simply supported beam structure compliance minimization problem was used. In this section, the traditional and two proposed methods were compared to perform a comparative study as follows:
Method1: Standard random initialization method (Traditional method)
Method2: 20×
N
elite initial population (Proposed method with
M
=20)
Method3: 40×
N
elite initial population (Proposed method with
M
=40)
where
N
is the population size of the topology optimization using a GA in Step 2 and
M
is the number of local optimum solutions obtained from the SIMP with the weighted sum method. For the proposed method, two methods with different parameter values for
M
(the number of local optimum solutions) were used to investigate the influence of parameter
M
. Parameter
N
of Step 2 was set to 50 empirically. The
N
=50 was adequate for finding the various genotypes corresponding to one phenotype. Twenty independent runs for each method were performed.
To verify the efficiency of the proposed method in the multiobjective topology optimization problem, the convergence and diversity were compared using the hyper volume ratio (HVR) and Paretooptimal frontier. By using HVR, it was possible to provide a qualitative measure of the convergence and diversity of the Paretooptimal solutions. Moreover, the computational cost of each method was compared to verify the efficiency of the proposed method.
Design domain and bar system for the simply supported beam
 4.1 Compliance minimization problem for the simply supported beam
A simply supported beam structure as shown in
Fig. 6
. was adopted as the verification problem for the proposed method.
Figure 6
shows the design domain and bar system of the structure. A design domain of 2H×H discretized into a 40×20 mesh was used, and the loadings F1=1 and F2=1 were assumed. The basic parameters assumed were Young’s modulus E=1.0, Poisson’s ratio
ρ
=0.3, density
t
=1.0, and the thickness of the structure t=1.0. In this simulation, a graph with four chain subgraphs was used, each of which consisted of seven vertices and six edges as shown in
Fig. 6
. The positions of the starting vertices were fixed at the support point and the position of the end vertices were fixed at the loading point, and all the other vertices were constrained to be within the design domain.
For the GA parameters of Step 2, a population size of 50, crossover probability of 0.9, and mutation probability of 0.2 were used. Total number of generations was 100. The optimization problem of topology optimization using a GA was to minimize formula (6).
On the other hand, for the Step 3 GA, a population size of 100, crossover probability of 0.9, and mutation probability of 0.2 were used. The optimization problem of topology optimization using a GA was to minimize compliances for the given two load cases of the structure with a 0.4 volume fraction.
 4.2 HVR
Hyper volume is generally used for the qualitative measurement of convergence and diversity in a multiobjective optimization problem
[2]
. Hyper volume is constructed with Paretooptimal solutions and a reference point in the objective space. The reference point was set as a vector constructed using worst objective function values. The metric HVR is usually introduced to calculate a normalized HV value. This HVR is the ratio of the HV of Q and of P
^{*}
as follows.
where HV(Q) is the HV calculated for the Paretooptimal solutions and HV(P
^{*}
) is the maximum (ideal) HV value. When Q=P
^{*}
, HVR=1. If Paretooptimal solution convergence and diversity are improved, HVR is close to 1. For the calculation of HVR in this study, the results of Method 1 with 600 generations, Method 2 with 300 generations, and Method 3 with 300 generations were used as the ideal HV(P
^{*}
), because generally the ideal HV(P
^{*}
) of a formula (8) is not obvious.
5. Result of verification
 5.1 Computational cost
Table 1
shows the average number of evaluations of the objective function per single GA trial and the average evaluation time of the objective function per single
Comparison of the computational time
Average number of objective function evaluations per single GA trial and average evaluation time for the objective function (40×20 meshed beam)
Average number of objective function evaluations per single GA trial and average evaluation time for the objective function (40×20 meshed beam)
Computational time
evaluation. In
Table 1
, the average number of evaluations of compliance was calculated, because the SIMP with weighted sum method of Step 1 was used until convergence and might have obtained different solutions if choosing different initial value. Because the GA in Step 2 was implemented up to 50 generations with 100 populations, the number of evaluations of formula (6) at Step 2 can be calculated as 50 (total number of generations) × 100 (population) × M (the number of local optimum solutions in Step1). Therefore when
M
=20 and
M
=40, the number of evaluations of formula (6) at Step 2 evaluation. In Table 1, the average number of evaluations of compliance was calculated, because the SIMP with weighted sum method of Step 1 was used until convergence and might have obtained different solutions if choosing different initial value. Because the GA in Step 2 was implemented up to 50 generations with 100 populations, the number of evaluations of formula (6) at Step 2 can be calculated as 50 (total number of generations) × 100 (population) ×
M
(the number of local optimum solutions in Step1). Therefore when
M
=20 and
M
=40, the number of evaluations of formula (6) at Step 2 were 100,000 and 200,000, respectively. In Step 3, the GA wasimplemented up to 300 generations with 100 populations. Thus the number of evaluations of compliance 1 and 2 at Step 3 was 30,000 when the total number of generation was 300.
In the case of the 40×20 meshed model, formula (6) and compliance 1, 2 average calculation times were 0.00004 and 0.07 s, respectively, and all the CPU times were measured on a desktop computer with a 2.8 GHz Intel Core i7 processor clock speed.
Table 2
and
Fig. 7
show the total computation time foreach method (methods 13). As shown in
Table 2
and
Fig. 7
, the additional computation times for the proposed method (the cost of Steps 1 and 2) for
M
=20 and
M
=40 were not large. When
M
=20 and
M
=40, the additional computation times of Step 1 of the proposed method were 4.62 and 9.59 s, and the additional computation times of Step 2 of the proposed method were 4 and 8 s. This additional computation time for the proposed method was negligible compared with that for the GA in Step 2, because the GA of Step 3 usually has a higher computational cost than a gradientbased algorithm such as SIMP, which was used in Step 1. Moreover, because formula (6) is the simple sum of the products of a matrix, the
Comparison of average, maximum and minimum HVR with 300 generations
Results for the Hyper Volume Ratio
Results for the Hyper Volume Ratio
Step 2 GA computational cost was very low, as shown in
Fig. 7
. Therefore, the computational cost of the proposed method was almost equal to that of the traditional method.
 5.2 Convergence and diversity of Paretooptimal solutions
In this study, to verify the convergence and diversity of Paretooptimal solutions, the HVR was used.
Fig. 8
shows the comparison of the HVR of the Paretooptimal solutions of Step 3. In
Fig. 8
, HVR was calculated from the result of the topology optimization of the GA with 300 evolved generations. The HVR results obtained from the traditional method with 400, 500, and 600 evolved generations are also shown in
Table 3
.
In
Fig. 8
, it is observed that the average (AVG) and the maximum (MAX) of HVR when
M
=20 and
M
=40 are much closer to 1 than the result for the traditional method, even when 600 generations have evolved. This shows that the Paretooptimal solutions obtained from the proposed method has better convergence and diversity than the traditional method. Moreover, the variance of HVR in
Table 3
shows that both variances of the proposed method (
M
=20,
M
=40) were smaller than for the traditional method. This means that the convergence ability of the proposed method for global solutions was more stable than for the traditional method. These results show that the proposed method has better convergence ability and more diverse Paretooptimal solutions with almost the same computational cost as the traditional method, implying that the efficiency of the proposed method was greater than for the traditional method.
Next, we investigated the influence of the parameter
M
(the number of the local optimum solution in Step 1), which was introduced for the elite initial population construction. To investigate the influence of the parameter
M
, the results of the proposed method with
M
=20 and
M
=40 were compared.
By comparing the average and maximum of HVR from
Table 3
, we observed that the average and maximums value
An example of evolution direction of genotypes that have the same phenotype
Results for the Pareto frontier
for
M
=40 were higher than the results for
M
=20, implying that if the number of local optimum solutions for
M
increases, varied and more converged Paretooptimal solutions can be found. Since parameter
M
was the number of local optimal solutions for Step 1, it related to the diversity of the phenotypes of the elite initial population. Thus if a high value of
M
is introduced, more varied phenotypes can be generated. Therefore, the parameter
M
affects the diversity of the elite initial population.
On the other hand, if parameter
M
increases, some very similar topologies may be obtained in Step 1, because some phenotypes overlap in Step 1. However, even when similar phenotypes are used in Step 2, different genotypes can be obtained in Step 2. For example, as shown in
Fig. 9
, two different genotypes correspond to the same phenotype. Thus, more varied genotypes can be introduced in spite of the overlapping of phenotypes in Step 1. Moreover, various genotypes in the GA lead to more converged and diverse Paretooptimal solutions in Step 3, because different genotypes that have the same phenotype evolve in different directions as shown in
Fig. 9
. However, a high
M
parameter leads to increased computational cost. Thus a balance needs to be found between the results for Paretooptimal solutions and the computational cost.
The results of the traditional and proposed methods (
M
=20,
M
=40) are compared in
Fig. 9
(a). The traditional method with 300, 400 and 500 generations, and the proposed method (
M
=40) with 300 generations, were compared in
Fig. 9
(b).
Figure 10
shows the result of the Pareto frontier. In
Figure 10
(a), if the traditional and proposed methods have the same computational cost for Step 3, the proposed method with
M
=40 has more convergent results than the traditional and proposed methods with
M
=20. In
Figure 10
(b), we observed that when the maximum number of evolution generations for the traditional method was increased from 300 to 600, its Pareto frontier asymptotically approaches the Pareto frontier of the proposed method with
M
=40. This result can be also described quantitatively by the average and maximum HVR as shown in
Table 3
. These results imply that the proposed method had better convergence ability, because the additional computational cost was negligible.
Figure 11
shows the results for the topology obtained by the proposed method with
M
=40 and 300 generations. In
Fig. 11
, we found that structurally proper topologies were obtained for the Paretooptimal solutions by the proposed method.
Although the proposed method had some additional
Topology results of the proposed method with 300 generations
calculation cost for generating the elite initial population in Steps 1 and 2, we established that the overall computational efficiency of topology optimization using a GA can be greatly improved when compared to the traditional method.
6. Conclusions
In this study, for efficient multiobjective topology optimization using a GA with barsystem representation, local optimum solutions obtained from the SIMP with weighted sum method were used as an elite initial population for the GA. To verify the efficiency of the proposed method, the traditional method with standard random initialization and the proposed method were compared using the simply supported beam example, and the HVR was calculated to investigate the convergence and diversity of Paretooptimal solutions. As a result, the proposed method performed better than the traditional method for the average, maximum and variance of HVR measured by Paretooptimal solutions and the Paretofrontier, with almost the same computational cost as the traditional method.
Tada Y.
,
Seguchi Y.
,
Soh T.
1986
“Shape Determination Problem of Structures Considering Multiple Loading Conditions : Externsion of Inverse Variational Shape Determination Method”
Transactions of the Japan Society of Mechanical Engineers. A
[in Japanese]
52
(476)
233 
239
DOI : 10.1299/kikaia.52.233
Deb K.
2001
Multiobjective Optimization using Evolutionary Algorithms.
John Wiley and Son Ltd
New York
Bendsøe M. P.
,
Kikuchi N.
1998
“Generating optimal topologies in structural design using a homogenization method”
Computer Methods in Applied Mechanics and Engineering
71
(2)
197 
224
Bendsøe M. P.
,
Sigmund O.
2003
Topology Optimization, Theory, Methods and Applications
Springer
Berlin
Huang X.
,
Xie Y.M.
2010
Evolutionary topology optimization of continuum structures: methods and applications
John Wiley & Sons
Chichester
Yamada T.
,
Nishiwaki S.
,
Izui K.
,
Yoshimura M.
,
Takezawa A.
2011
“A Structural Optimization Method Incorporating Level Set Boundary Expressions Based on the Concept of the Phase Field Method.”
Journal of Environment and Engineering
6
(3)
567 
578
DOI : 10.1299/jee.6.567
Otomori M.
,
Yamada T.
,
Izui K.
,
Nishiwaki S.
2011
“Topology Optimization Based on the Level Set Method Using Mathematical Programming”
Transactions of Japan Society of Mechanical Engineers Series C
[in Japanese]
77
(783)
4001 
0414
DOI : 10.1299/kikaic.77.4001
Wang S. Y.
,
Tai K.
2005
“Barsystem representation method for topology optimization using the genetic algorithms”
Engineering Computations
22
(2)
206 
231
DOI : 10.1108/02644400510585493
Deb K.
,
Agrawal S.
,
Pratap A.
,
Meyarivan T.
2000
“A fast and elitist multiobjective genetic algorithm for multiobjective optimization: NSGAII”
in: Proceedings of the Parallel Problem Solving from Nature VI Conference
Paris
849 
858
Goldberg D.
1989
Genetic Algorithms in Search, Optimization, and Machine Learning
Addison Wesley Longman
Deb K.
,
Anand A.
,
Joshi D.
2002
“A computationally efficient evolutionary algorithm for realparameter optimization”
Evolutionary Computation
10
(4)
371 
395
DOI : 10.1162/106365602760972767
Deb K.
,
Goyal M.
1996
“A combined genetic adaptive search (GeneAS) for engineering design”
Computer Science and Informatics
26
(4)
30 
45