Efficient Combination of Topology and Parameter Optimization

This paper presents a combination method of Particle Swarm Optimization (PSO) and topology optimization. With this method a better result can be achieved compared with the sequential application of these two optimization methods. It inherits the ability in finding global optimum from PSO and also suits for discretized design domain. Some special schemes are used in order to provide higher computation efficiency. This method has only been tested with a convex optimization problem. The application in case of a concave problem will be a future study.


Introduction
Parameter optimization and topology optimization are two powerful tools in mechanical design, which affect remarkably the shape and the performance of a structure.Traditionally, the usage of these two procedures is separated and sequential.The structure is often parameter-optimized firstly, and then topology-optimized.However, the sequential method has potential insufficiencies.The general size of a structure, which is usually regarded as a design domain in topology optimization, is always changed due to the parameter optimization.The change in geometric size is sometimes useless or even conflicting as a preparation to topology optimization.In this paper, an efficient method of combined optimization is presented.It executes the parameter optimization with continuous or discretized variables and the topology optimization simultaneously, and thus achieves a better result than the sequential optimization.The paper is organized as follows: The theoretical basis of topology and parameter optimizations is introduced in Section 2. For a better presentation and explanation, an optimized case with variable bearing hole position is shown in Section 3. The implementation of the combined optimization and the methods for increasing the efficiency are proposed in Section 4. Finally, conclusion and future work are provided.

Finite Element Method
The finite element method (FEM) is a numerical method that divides a domain into small connected sub-domains (elements), such that differential equations can be solved by solving a set of polynomial equations.With the development of computer science, more and more complicated mechanic problems in different fields are now solvable with FEM.One feature of FEM we need here is the definition of parameters for each element individually.This is important for the implementation of topology optimization, since the design variables should be set for every single element locally.

Topology Optimization
Topology optimization is a mathematical approach that optimizes the material layout within a given design space, for a given set of load and boundary conditions.It is widely used in the field of shape or structure design, especially for the lightweight design without losing much stiffness.Different kinds of topology optimization algorithms have already been carried out.Xie and Steven introduced the Evolutionary Structural Optimization (ESO) to eliminate the less effecting elements of a structure [1] [2].Querin et al. proposed Bi-directional Evolutionary Structural Optimization (BESO) as an improvement of ESO that can also add elements in a design domain to compensate the structure defect caused by a single directional elimination of elements [3].Bendsøe presented the Solid Isotropic Material with Penalization (SIMP) method to update the element-density-dependent material properties in order to reach the final optimum [4].The penalization factor is applied to reduce the amount of intermediate densities.However, the introduction of penalization usually has a negative effect on both computational time and final results.In order to reduce such influences, Dadalau [5] introduced SIMP with an adaptive penalization scheme.In this paper the SIMP-model is used.As shown in Equation (1), it regards every element density e ρ as a single optimization variable.All the variables are updated with a specific scheme which is derived by Lagrange multipliers [6].
Here k ρ denotes the density at k-th step.Variables η and ζ control the possible change in one iteration with the typical values of 0.5 and 0.2.This scheme adds material to the structure if where k B is the displacement at k-th step, which is calculated from the equilibrium function of the structure, Λ is a Lagrange multiplier and p is the penalization factor.

Parameter Optimization
Classical methods of parameter optimization, which are based on sensitivity calculation and the gradient descent method like the Newton method, are already well developed with high accuracy.However, the drawbacks are also obvious.It is expensive to calculate the needed sensitivities, if the variables are not explicit in the objective function.Moreover, in many discretized situations, the loss of continuity is usually problematic.Furthermore, the global optimum is difficult to reach, since it can be easily trapped by a local optimum.The Particle Swarm Optimization (PSO) method, known as a purpose optimization method, was originally proposed by Kennedy and Eberhart [7].It optimizes a problem by looking for the best candidate solution in a search domain with a set of particles.The general equations for moving particles can be written as follows, ( where , i k v is the search velocity of i-th particle for k-th step, , i k x is the position of i-th particle for k-th step, ω is the inertia of particle, i p is the best point i-th particle ever met, i g is the global best point among per- sonal bests, p φ and g φ are weight parameters.Particles are randomly generated at first and they find the next step with the influence of personal and global best points.After several iterations, all the particles could gradually converge to the optimal design point.Various researches about the behavior due to parameter selection in PSO have been presented, such as Shi and Eberhart [8], Van Den Bergh [9] and Trelea [10].As the choice of position for the next iteration is relatively random, this method can be easily applied in both the continuous and the discretized problem.x y

Problem Setting
where c is the compliance as the objective function, c x and c y are the coordinates of the hole center, V is the maximum accepted volume after optimization, min ρ is a value close to zero to avoid numerical problems due to the zero density and N is the element number.The density for each element e ρ is related to the Young's modulus of each element with penalization factor.
In the sequential method the best hole position is found by parameter optimization before executing the topology optimization.The final result is calculated by applying topology optimization on this hole position.The question is, whether the best position found by parameter optimization is still the best one after the topology optimization.Two enumeration tests were made as verification, which cover all the possible results in this optimi-Figure1.A simple structure with bearing hole to be optimized.) with the value of 63.63, which is significantly lower than the former one.This means, at least in some cases, the results from the sequential usage of parameter and topology optimization are far from the real optimum.For this reason, we proposed a new method to find the real optimum.

Implementation of Combination Method
The general idea of the combination method is based on parameter optimization with the PSO method.The topology optimization is used as a way for calculating compliance, which is also the objective function of PSO.The algorithm is shown in the flowchart in Figure 3.The method is implemented with a Matlab program based on the finite element example code in [11] and the topology optimization code in [12].However, due to the low  efficiency of Matlab with its default large-scale matrix operations, the matrix storage scheme and solving method is modified according to [13] to accelerate the whole process.
In the first step, the bearing hole is set randomly on several positions.On those positions, the fast topology optimization is executed in Step 2. "Fast" means less iterations are needed to converge compared with normal topology optimization by some special schemes, which are mentioned in the following sections.The new positions are generated in Step 4 through Equation ( 3), until all the particles (hole positions) converge to the same position (Step 3).Finally, a slower topology optimization is applied on the global best position to achieve the optimum structure (Step 5).

Passive Element Method
The task of moving the hole position during the optimization process must be simple to handle.Normally, the empty area is not included in mesh.When the hole position is changed, meshing and assembling of the stiffness matrix need to be repeated, which decreases the efficiency greatly.This problem can be solved by the passive element method, with which the hole in the design domain can be defined as "zero" density element ( min ρ in the program) instead of a description of the geometric information.The mechanism of the passive element is to set the density in the hole area from the last iteration to one and set those in a new hole area to zero, when the position is changed during the searching procedure.The nodes in the lower part of the hole can also be easily selected to add the load.With the passive element method, only the material density parameter is applied to evaluate the global stiffness matrix, the remeshing process is no longer needed.

Unified Penalization Factor
The penalization factor p plays an important role in topology optimization.By 1 p = , which means the density and the Young's modulus of an element are exactly linearly related, the optimality of the topology optimization is guaranteed.However, the optimized structure will be filled with large numbers of intermediate densities, which are not realistic to be manufactured.With a large penalization factor, such as 3 p = , the structure will be more realistic and closer to a 0-1 structure.Nevertheless, high penalization factors also have some disadvantages.Table 1 shows that the compliance by 3 p = is worse, and it converges more slowly than that by 1 p = .Due to the desired efficiency, the topology optimization during the PSO should be fast, thus the penalization factor within the fast topology optimization (Step 2 in Figure 3) is fixed to one.

Simplified Convergence Condition
The convergence condition also has a great effect on the convergence speed.Normally, the topology optimization is terminated when the largest change of the element density is smaller than a preset value [6] [12].Under this condition, a very fine result can be achieved.Nevertheless, the change of density can fluctuate on some positions of the domain, which makes the process difficult or sometimes impossible to converge, although the change of compliance is already very small.The simplified convergence condition terminates the optimization loop, if the change of compliance is smaller than a preset value.According to Table 1, it decreases the convergence loops greatly without losing much accuracy.Through the simplified convergence condition and unity penalization, the fast topology optimization will converge in around 10 loops instead of almost 100 loops each time.Along with the flexibility of the program, the strict condition with a high penalization value can still be used to get the accurate final structure (Step 5 in Figure 3).

Test Results
The optimized hole position found by the combination method is converged on (93, 40) with the compliance value of 63.64, which is just one element width beside the best position of the reference (Figure 2(b)).Actually, due to the property of the PSO method, the converged result is not always the best but very close to it.Figure 4 describes the two structures after two separate steps in a sequential method.The final structure resulting from the combined method is shown in Figure 5.

Conclusion
In this paper, an algorithm is implemented to combine the topology optimization and PSO-based parameter optimization.Different schemes are applied to increase the computation efficiency.This algorithm was tested by an example case with a movable bearing hole.It finds the acceptable optimum result, whose compliance after the topology optimization is significantly lower than the sequential optimization.As the calculation of compliance for each position in one PSO iteration is absolutely independent, this algorithm can be easily extended to parallel computation, with which the computation speed can be further increased.It is suitable for the structure design with a convex design domain.How to extend the algorithm to the case with a concave boundary will be studied in the future.

Figure 1
Figure1shows an example case to be optimized.The length of this structure is 130, the height is 100.It has a bearing hole with a radius of 15 and the wall thickness of the hole is 10.The hole can be moved in a 51 × 51 rectangle.The structure is meshed with 2D bilinear quadrilateral elements (Q4) with an element size of 1. Thus the design domain is also discretized into 51 × 51 nodes.Along the lower half of the hole, a distributed force of 3 is added in −Y direction.At the right bottom corner, two unit forces in +X and +Y direction are applied.The left boundary is fixed in both X and Y direction.The task is to reduce the volume by 50% and to optimize the hole position, such that the compliance is minimized.The volume of the hole is also included as a part of the deleted volume.Mathematically the optimization problem can be expressed as follows:( ) In the first test, the compliance of the structure is calculated with the hole position on each node of the domain without topology optimization.The value distribution is shown in Figure 2(a).The second test shown in Figure 2(b) is with topology optimization with penalization factor 3 p = .As can be seen from Figure 2(a), the best hole position is located on (86, 73) before the topology optimization, which means this position would be found by the sequential method and used as the hole position for further topology optimization.And from Figure 2(b) we can easily see that the final compliance on (86, 73) is 71.51.But the optimized position from Figure 2(b) is (94, 41

Figure 2 .
Figure 2. Results of the objective function raster.This is a mapping of compliance value within the design domain (zone within dash line in Figure 1).(a) stands for the compliance mapping before topology optimization while (b) stands for the one after topology optimization.

Figure 4 .
Figure 4. Structures before and after topology optimization in sequential method.The left figure shows the actual structure after getting the best hole position by parameter optimization.The right figure is the final structure acquired by sequential method.

Figure 5 .
Figure 5. Final structure from the combination optimization.

Table 1 .
Compliance and convergence speed test.