Improved Topological Optimization Method Based on Particle Swarm Optimization Algorithm

To solve the grid-dependent phenomena, such as singular structure, checkerboard structure, and gray element phenomenon, that often appear in the topology optimization of the SIMP interpolation model of the variable density method, the particle swarm optimization (PSO) algorithm was introduced to the SIMP model. In addition, a new inertia weight was designed to calculate the fitness value of structural flexibility for the redivision population of the elements after sensitivity filtering. The MBB beam was designed by different topology optimization methods. The results show that the improved algorithm can effectively avoid iteration falling into a locally optimal solution, and the grid-dependent and gray element phenomenon of the optimized structure are suppressed to a certain extent. The optimized structure has good stability, the number of gray elements is reduced by 6%~ 17%, and the boundaries of the structure are clearer.


I. INTRODUCTION
Topology optimization is a mathematical method to optimize material distribution in the specified areas according to given load conditions, constraints, and performance indexes to achieve ''good use of things''. It is an important approach to achieve lightweight and structural integration design of additive manufacturing products. The variable density method is the most common topology optimization method at present, among which the Solid Isotropic Penalty Microstructure Model with Penalization (SIMP) is widely used because of its simple concept and easy implementation. However, in the SIMP, grid-dependence and the checkerboard phenomenon may occur due to different discretization degrees of the structure, or a large number of gray elements may appear due to the average effect, which brings great difficulties to manufacturing.
According to Wen [1], grid-dependent phenomena, such as checkerboard and other singular structures, are caused by different degrees of mesh division, and they usually occur together. Methods that can suppress singular structures can also restrain the checkerboard phenomenon.
To suppress grid-dependence, Acar [2] proposed a microstructural design method using a discrete adjoint sen-The associate editor coordinating the review of this manuscript and approving it for publication was Kuo-Ching Ying . sitivity analysis scheme. The design optimization of the microstructure was performed with the accompanying sensitivity analysis, and then the design optimization based on the gradient was performed. The output of the sensitivity analysis was used to enhance the material properties and reduce the dependence of modeling on the finite element mesh. In the study of Li [3], parallel computing and grid adaptive technologies were properly combined, and a parallel distributed open-source framework for topology optimization of the fullsize 3D structure was proposed. The optimization algorithm reduced the dependence on initial guessing grid resolution to a certain extent. Qi [4] proposed an improved SIMP interpolation model to solve the problem of grid-dependence under a large-scale grid, which improved the stiffness of the structure, especially when the grid size was small. Guo [5], [6] introduced a new convolution factor into the sensitivity filtering method to weaken the boundary diffusion problem and improved the checkerboard phenomenon in topology optimization. Rui [7] proposed a layered double punishment method and combined it with the layered mesh subdivision strategy to improve the grid-dependence and improved the convergence rate.
In the improvement of the gray element problem, Fuchs et al. [8] introduced inverted variables and SRV (the sum of the recip rocal variables) constraints into the SIMP method to reduce the generation of gray elements in the VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ optimization process. Based on Fuchs [8], SIMP and SRV methods were combined by Xiang [9], which satisfied the shortcoming of SRV requiring a high initial value of design variables, and further proposed a method combining SIMP and SSV. According to the generation principle of grayscale elements, Zhi [10] designed a dual SIMP method, taking the optimization result with sensitivity filtering as the intermediate value, and then using the SIMP method without sensitivity filtering for the second optimization, which effectively suppressed the generation of gray elements and obtained a clear structure. However, the number of iterations increases obviously, which increases the amount of calculation. To suppress gray elements in view of the low computational efficiency caused by the Heaviside projection function frequently called density matrix, Xiang et al. [11] designed a new method, which could significantly improve the phenomenon of grayscale elements. Mao [12] proposed a new optimization strategy by combining the topology optimization and shape optimization of the structure, reducing the proportion of gray elements in the topology optimization.
In the improvement of the topology optimization algorithm, many solutions have been proposed by scholars and some results have been achieved, but few studies have been done on the combination of intelligent algorithm, especially particle swarm optimization algorithm. Most scholars apply swarm intelligence algorithms in structural layout and site selection, but seldom in structural optimization, or the improvement effect is relatively simple, even at the cost of other advantages. Therefore it is still a hot research topic to explore a better topology optimization algorithm to promote its application in production. The main contributions of this study are as follows: (1) The PSO algorithm was introduced into the traditional variable density method SIMP interpolation model by MAT-LAB programming.
(2) Linear inertial weights are introduced into the PSO to balance the capabilities of global and local search in iterative calculations.
(3) A new formula of inertia weight is proposed to solve the problem that the linear decreasing inertia weight has insufficient punishment for element density in the iterative process.
(4) The improved algorithm is effective in reducing griddependence and gray element phenomenon in structural topology optimization.

II. SIMP METHOD A. SIMP MODEL
The SIMP interpolation model is an interpolation model based on the variable density method proposed by the homogenization method, as shown in (1). By introducing penalty factor P, the intermediate density values between 0 and 1 converge to 0 or 1 so that the topology optimization model of continuous variables can better approximate 0-1 discrete variables [13].
where E (x i ) represents the elastic modulus of the material after interpolation calculation. E 0 and E min represent the elastic modulus of the solid part and pore part respectively, and the pore part is generally 1 of the solid part. x i represents the cell density, which is between 0 and 1. At the beginning of the calculation, the initial volume fraction is defined as the density of each element. Under different initial volume fractions, different penalty factors are used to punish the cell density, as shown in Fig. 1.
According to Fig. 1, as the penalty factor P increases in the SIMP interpolation model, the variation trend of the elastic modulus of the element calculated by interpolation is closer to the step function, which means that an increasing number of elements will be in the range of truncation at the beginning of the iteration (the elastic modulus is lower than the value of the empty element). This also increases the chance that iterations will fall into local optimality or cause numerical instability. In practical examples, the value of P is generally between 3 and 7, which also leads to insufficient punishment of intermediate density units, resulting in a large number of gray elements. More elements should be allowed to participate in the iterative process.

B. FILTRATION
The grid-dependency problem is mainly solved by the filtering method in the existing topology optimization algorithms. Its essence is to use a certain calculation method to carry out a weighted average of a physical value of all elements within a circle and replace the value of the central element with the average value. The circle is centered on the central element and has a filter radius of R min . A schematic diagram of the filtration method is shown in Fig. 2. When R min is smaller than the element size, only the central element is selected. The change rate of the filtered objective function is equal to the original one, so the filtering is invalid, showing the checkerboard phenomenon. In practical examples, 1∼3 times the element size is usually taken as R min . Filtering greatly reduces the probability that the density of the center goes to 0 or 1, and effectively inhibits the checkerboard phenomenon caused by the ''0-1'' distribution. However, due to the obvious randomness of the sensitivity value or density value of adjacent elements, the weighted average value cannot be guaranteed to develop in the same direction all the time and will even further aggravate the generation of gray elements. The suppression effect on gray elements is reflected by the proportion of gray elements in all elements.

C. PROJECTION
The projection method is a common method to solve the gray element phenomenon. A method of nonlinear mapping of filtered density values using step functions was proposed by Guest in 2004 [14]. To avoid the optimization falling into the locally optimal solution, the projection function commonly used is shown in (2).
where ρ e represents the density of unit e after projection, and β stands for projection slope. The larger β is, the closer the projection function is to the step function. θ β represents the projection point. When θ β = 0.5, the curve of (2) under different β is shown in Fig. 3. It is a linear function when β = 1. With the increase in β, the function curve is closer to the step.
In the projection method, the density values of grayscale elements are divided into two intervals and concentrated to 0 or 1, which will neither make a large number of elements empty, resulting in an insufficient structural stiffness, nor make a large number of elements solid, resulting in insignificant topology optimization effect. However, the obvious disadvantage is that the effect of the intermediate density cells to 0 or 1 concentration is not obvious when β is not large enough, and further operations need to be carried out  by other punishment methods. When β is too large, the step property is enhanced, and too many elements exit the iteration at the beginning of the iteration, which makes the calculation fall into the locally optimal solution and result in a singular structure. In this paper, the inhibition of grid-dependence is represented by a specific form of topology.

III. MBB OPTIMIZATION
As shown in Fig. 4, the ratio of length to width of the MBB beam is 4 : 1. Fixed constraints are adopted in the lower left and right corners, and a downward unit load is applied at the midpoint of the upper boundary.
In this example, the element density is taken as the design variable, the minimum flexibility (i.e., the maximum stiffness) is taken as the optimization objective, and the geometric volume is taken as the constraint condition, which was set as 0.4.
The constraint conditions are described as follows, where nelx and nely are the number of elements in the horizontal and vertical directions of the beam, respectively.
Topology optimization using the SIMP interpolation model requires R min , P, and mesh density. In this paper, the suppression effect of the improved algorithm on grid-dependence and gray elements is mainly considered, while the selection of parameter collocation is secondary, so the coupling effect between parameters is not discussed. TABLE 1 ∼ TABLE 3 show the influence of R min , P, and mesh density on the topological structure. TABLE 1 shows the evolution of the topology optimization structure when P = 3, the mesh density nelx * nely = 80 * 20, and the filter radius R min changes between 1 and 3. When R min = 1, a large number of checkerboard phenomena appear in the structure, and when R min is too small or too large, singular structures appear. When R min ≥ 1.1, the checkerboard phenomena disappeared, but more ''0'' or ''1'' elements participated in the weighted average and evolved into gray elements, increasing the proportion of gray elements in the optimized structure.
Next, the effect of the penalty factor on topology is considered when R min is 1.1. TABLE 2 shows the evolution of the topology optimization structure when R min = 1.1, the mesh density nelx * nely = 80 * 20, and P changes between 1 and 3.5. When P is small, a large number of gray elements appear in the structure. With the increase in P, the punishment intensity of the model on the intermediate density element gradually increases, resulting in the number of gray elements decreasing. When P > 3, too  many elements withdraw from iteration prematurely, leading to failure to jump out of the locally optimal solution, and some singular structures appear in the structure, showing a certain grid-dependence.Then, the effect of the mesh partition density on the topological structure is considered when P = 3 and R min = 1. 1.  TABLE 3 shows the evolution trend of the topology optimization structure with the increasing density of meshing when P = 3 and R min = 1.1. With the increase in mesh density, the model structure is more detailed, and the number of gray elements increases, but the proportion decreases under the same parameters. With the increase of grid density, the structure pattern changes greatly, and the emergence of fine structure also gradually increases, reflecting the grid dependence.
According to TABLE 1 ∼ TABLE 3, the filtering method and projection method can improve the numerical instability of topology optimization, but they have certain limitations. For example, the filtering method eliminates griddependence, but also increases the number of gray elements. The numerical instability in the SIMP topology optimization method can be improved by adjusting parameters, but it is always accompanied by increasing gray elements.

IV. IMPROVED SIMP
According to the above analysis, the main task of this section is to find an appropriate method based on sensitivity filtering so that most elements can participate in the iterative calculation. Based on this idea, the specific gravity function w(x i ) used in the filtering method was taken into account, with different weights set according to the distance between element i and central element e. Here, the weight can be set to make the weight related to the number of iterations. Most elements participate in the calculation at the beginning of the iteration, while only some units remain and participate in the calculation at the end of the iteration, which is in line with the topology optimization theory.

A. IMPROVEMENT BASED ON PSO
PSO was introduced into the SIMP to improve the algorithm, and the algorithm flow is shown in Fig. 5. To realize particle updating, the element compliance values calculated by the SIMP method are compared with those calculated by the PSO algorithm.
In PSO, the number of elements is regarded as the population size. The range of the particle's position X is the range of cell density, 0 ≤ X ≤ 1. The velocity of the particle has a range of −X max ≤ V ≤ X max . The learning factor c 1 = c 2 = 2. During the iterative calculation: (1) Initialize the position and velocity of the particle; (2) Sets the current particle position value to the individual optimal value and the current column position value to the global optimal value; (3) The flexibility of the current particle is calculated and compared with that calculated by the SIMP method; (4) The position and velocity are updated until the termination condition is satisfied; (5) Individual and global extremes are obtained.

B. IMPROVED INERTIA WEIGHT
In PSO, inertia weight reflects the ability of particles to maintain or inherit the previous state of motion. To balance the ability of global and local search, the linear decreasing  inertia weight was first proposed by Shi et al. and applied to the particle swarm optimization algorithm [15]. The inertia weight of linear decline is shown in (3), denoted as LIN.
where, ω(t) represents the weight coefficient related to the number of iterations, ω 1 represents the initial inertia weight, ω 2 represents the inertia weight at the maximum number of iterations, t represents the current number of iterations, and T max represents the maximum number of iterations. In general, the algorithm performance is better when ω 1 is large and ω 2 is small. In this way, as the iteration progresses, the inertia weight decreases slowly. Which not only maintains the strong global search capability at the beginning of iteration, making it more likely to traverse the solution space, and avoid falling into the locally optimal solution, but also retains strong local search ability at the end of iteration, which enables it to lock the optimal solution more quickly and accurately.Therefore, a new weight formula is proposed by adding exponential α to the terms related to iteration times according to the evolution trend of the exponential function, as shown in (4). where, the weight formulas for exponents 2 and 3 are denoted EXP-2 and EXP-3, respectively.  Set ω 1 = 0.9, ω 2 = 0.4, and T max = 300. The variation rules of inertial weight values of the LIN, EXP-2, and EXP-3 with the number of iterations are shown in Fig. 6.
At the beginning of the iteration, the weights of EXP-2 and EXP-3 are larger and decrease slowly, which is suitable for global search to avoid falling into the locally optimal solution. while at the end of the iteration, the weight value decreases rapidly, which is suitable for local search to obtain an accurate solution.
When EXP-2 is used for topology optimization of structures with different mesh densities, the resulting topology is shown in Fig. 7. It can be seen that when the index α = 2, there is a great difference in topology structure, and the inhibition ability of grid-dependence is insufficient, so the index should be further increased.

V. OPTIMIZATION WITH IMPROVED ALGORITHM
In this section, the classical SIMP method and two improved algorithms, LIN and EXP-3, are used to recalculate the above examples when P = 3, R min = 1.1, and nelx * nely = 80 * 20. In the improved PSO algorithm, the initial inertia weight ω 1 = 0.9, the ending inertia weight ω 2 = 0.4, the maximum iteration number T max = 50, and the learning factor c 1 = c 2 = 2. The fitness curves are shown in Fig. 8. According to the results, when P = 3 and R min = 1.1, the number of iterations of the LIN method is approximately 7, and EXP-3 39. The EXP-3 method can effectively avoid iteration falling into a locally optimal solution.
The simulation results of the three topology optimization methods are shown in TABLE 4. According to the results, the grid-dependence caused by the parameter setting is improved. The flexibility of the structure is reduced, that is, the stiffness is increased, and the number of gray elements is also reduced after adopting PSO. When the LIN method is used, the grid-dependence has not been eliminated, and some singular structures still exist. When the EXP-3 method is used, the grid-dependence isfurther reduced, no singular structures appear, the number of gray elements is further reduced, and the flexibility is also slightly reduced; that is, the structural stiffness is improved.
Compared with the LIN method and the original method, the number of gray elements is reduced by 5.84% and 17.02%, respectively, when the EXP-3 method is used, which effectively improves the grid-dependence in the calculation and greatly reduces the generation of gray elements.
The structure dependence phenomenon caused by different mesh densities is studied below, and the results are shown in TABLE 5, where nelx * nely = 80 * 20 is shown in TABLE 4.
The comparison of gray unit proportion and structural flexibility is shown in Fig. 9.
According to the simulation results, the main structure obtained by using the improved inertia weight formula, i.e., EXP-3, is more obvious, and the tiny structures that is not conducive to manufacturing is avoided. With the increase in the number of elements, the topology structure obtained remains roughly unchanged and becomes more sta- ble, indicating that the improved algorithm can weaken the grid-dependent phenomenon caused by different mesh densities to a certain extent. In addition, when the improved algorithm is used for calculation, the number of gray elements is reduced by 6%∼17%, and the boundary of the topology structure is clearer. The flexibility of the structure is basically unchanged or slightly reduced; that is, the stiffness of the structure remains stable or increases slightly, and the mechanical properties have no tendency to decrease.
According to the above simulation results, a better result can be obtained when the index in the weight formula is 3, which meets the needs of structural optimization under the current constraints without further adjustment of the value of the index.

VI. CONCLUSION AND FUTURE DIRECTIONS
In this paper, PSO was introduced into the classical topology optimization algorithm by analyzing the shortcomings of the filtering method and projection method to solve the problem of grid-dependence and gray element instability when the SIMP interpolation model is used in topology optimization. On this basis, a new inertia weight was proposed, which makes more elements participate in the calculation at the beginning of the iteration to avoid falling into the locally optimal solution, and the elements decrease rapidly at the end of the iteration to speed up the iteration and find the exact solution. The main advantages of the improved algorithm are as follows: (1) Topology optimization using the improved algorithm weakens the generation of singular structures, and the stiffness of the topology structure increases slightly.
(2) The mesh adaptability of the improved algorithm is better. With the increase in the number of cells, the structure is more specific, but the topological structure is stable and the main structure is clear.
(3) The number of gray elements is reduced by 6%∼17%, and the topology structure is clearer. However, the improved algorithm still has some limitations as follows: (1) The improved algorithm has some defects in the adaptive selection of parameters, such as penalty factor P, filter radius R min , inertia weight ω and learning factor c.
(2) The improved algorithm is not suitable for topology optimization of 3D structures.
(3) At present, optimization problems in fluid, electromagnetic and other fields cannot be solved only by algorithms.
In future research, the adaptive particle swarm optimization algorithm should be considered and applied to the improved algorithm to adapt to the actual requirements of different optimization conditions. With the change in the evolutionary state, particle swarm weight values and learning factors are redefined to balance the global and local search capabilities of particle swarm optimization [16], [17]. Topological optimization of 3D structures can be extended on this basis. For optimization problems of other physical fields, cosimulation can be used to solve them, such as the cosimulation of MATLAB and COMSOL [18].
JIE GUAN is currently pursuing the master's degree in mechanical engineering with the College of Power Engineering, Naval University of Engineering, Wuhan, China. His research interests include additive manufacturing and topology optimization.
WENQUN ZHANG received the Ph.D. degree from the College of Weapons Engineering, Naval University of Engineering, Wuhan, China, in 2004. He is currently an Associate Professor with the Naval University of Engineering. His research interests include the optimized design of ship machinery, additive manufacturing, and topology optimization.