Modified Proportional Topology Optimization Algorithm for Multiple Optimization Problems

material interpolation scheme and introduce the Heaviside threshold function based on the PTO al-gorithm. To confirm the effectiveness and superiority of the presented algorithms, multiple optimization problems for the classical MBB beam are solved, and the original PTO algorithm is compared with the new algorithms. Numerical examples show that MPTOc, MPTOs and MPTOm enjoy distinct advantages over the PTO algorithm in the matter of convergence efficiency and the ability to obtain distinct topology structure without redundancy. Moreover, MPTOc has the fastest convergence speed among these algorithms and can acquire the smallest (best) compliance value. In addition, the new algorithms are also superior to PTO concerning suppressing gray-scale elements.


Introduction
Structural lightweight design is a significant and effective way to achieve energy conservation and environmental protection [1].There are a wide variety of approaches to achieve the lightweight design of structures, and the mainstream approaches widely studied by scholars mainly include the following three types: shape, sizing and topology optimization [2,3].Topology optimization is the optimal way to make full use of materials among the three methods.Generally, topology optimization is considered as an excellent approach to find the main load transfer path via scientific calculation within a prescribed design space and under specified constraints [4].
It is well known that the material interpolation scheme plays a key during the process of topology optimization.the scheme is to penalize the intermediate densities in the design variables, so that the densities of more design variables tend to approach 0 or 1.Nevertheless, whatever SIMP (solid isotropic material penalty) or RAMP (rational approximation of material properties) is used, topology optimization is often accompanied by some numerical instability phenomena, mainly including mesh dependence, checkerboard pattern, and too many gray-scale elements [5].In order to solve these problems, Araujo et al. [6] used a topology optimization algorithm derived from generalized finite-volume theory, which is implemented via applying mesh independent filter, to avoid mesh dependency.Cui et al. [7] presented the modified optimality criterion (OC) method related to the density filtering based on tanh-function to solve the SIMP model, so as to suppress the generation of intermediate density elements as well as improve the computational efficiency.
To get an ideal topology optimization result, choosing an appropriate numerical optimization algorithm is also an essential part of the optimization process.After decades of development in topology optimization, various numerical optimization algorithms have been proposed.Peng et al. [8] proposed the OC method for multi-constraint topology optimization.Gonçalves et al. [9] used a successive linear programming (SLP) approach to optimize the placement and shape of traditionally embedded piezoelectric actuators.The control variables about the iron density distribution on the cross section of the magnetic channel are optimized by using Method of Moving Asymptotes (MMA) [10].All of the above algorithms belong to gradient-based optimization methods, which are essential for calculating the sensitivities of objective function and constraint function during the process of solving corresponding optimal problems.In this way, the computation about sensitivity will produce computational inconvenience and complexity during the iterative process of topology optimization to a certain extent, especially when dealing with some complicated optimization problems with more rigorous and complex sensitivity calculations [11].For example, in the exterior acoustic problem, the sensitivity analysis becomes complicated as a result of the variation of the interface between acoustic medium and structure [12].The implementation of the sensitivity analysis may become infeasible owing to a great many stress constraints applied [13].Besides, the researchers can read the paper opened by Sigmund [14], in which the gradient versus non-gradient approaches are discussed in detail.
In the light of the drawbacks of the gradient-based optimization method, nowadays a growing number of scholars have gradually focused on non-gradient-based optimization method.Dinh et al. [15] employed a multi-phases optimization design approach applied to compliant mechanisms.In this approach, the topology optimization is combined with neural network algorithm, intelligent modelling as well as finite element method.Millán Páramo and Begambre Carrillo [16] used an improved simulated annealing algorithm to deal with structural optimization problems.Biyikli and To [11] suggested a new non-sensitivity optimization approach named proportional topology optimization (PTO) algorithm.And this algorithm is classified as a heuristic approach that distributes density variables proportionally to elements by prescribed rules.
For topology optimization problems, most of the topology optimization research work aims at solving the minimum compliance optimization problems under material volume constraint and the stress-constraints optimization problems [17].For instance, an effective algorithm combining density filter and SIMP method for addressing the stressconstrained topology optimization problem was presented by Le et al. [18].Mela and Koski [19] developed topology optimization for trusses with multi-load conditions, considering and tackling the optimization problems of minimizing compliance under material volume constraint and the minimum weight problem under stress constraints, respectively.Tovar and Khandelwal [20] introduced a control-based optimization algorithm utilized to resolve the minimum compliance problem.In addition, some scholars use the displacement, the frequency and temperature as constraints in topology optimization algorithms [21,22].According to the published literature, there are few reports on solving topol-ogy optimization problems under both stress and compliance constraints.
In this work, we propose modified proportional topology optimization (MPTO) algorithms by modifying the material interpolation scheme in the original PTO algorithm and adding the Heaviside threshold function.MPTOc and MPTOs in the offered algorithms are designed to tackle the minimum compliance problem and the stress-constraint problem, respectively, while MPTOm is employed to deal with the minimum volume fraction problem where compliance and stress are simultaneously constrained.Furthermore, the effectiveness of each new algorithm is validated.Compared with the original PTO algorithm, the new algorithms enjoy the great advantage in convergence efficiency and distinct topology structures without redundancy.

Mathematical model of topology optimization problems
Mathematical models of three kinds of optimization problems, i.e., the minimum compliance problem under volume constraint, the minimum volume fraction problem under stress constraint, and minimum volume fraction problem under compliance and stress constraints, are defined in the following sections.

Minimum compliance problem under volume constraint
In the light of the published literature [11,23], the optimization model of minimum compliance under volume constraint can be described: where C denotes the compliance (and also the objective function), U and F indicate the global displacement vector and external force vectors of the structure respectively, K represents the global stiffness matrix, N indicates the number of all elements, E(xi) is the Young's modulus whose value can be obtained from the Eq. ( 10) in this paper, ui is the elemental displacement vector, k0 denotes the elemental stiffness matrix, vi is elemental volume, V represents the total material volume in the design domain, xi is the elemental density (and also the design variable), and xmax and xmin are its upper and lower boundaries respectively.Typically, xmax =1 while xmin needs to be given a proper value to prevent the occurrence of stiffness singularities.

Minimum volume fraction problem under stress constraint
The stress problem is to minimize the material volume under stress constraint, the optimization problem can be mathematically expressed as: , if 0 01 where, besides the nomenclature provided for the minimum compliance problem, σi indicates the elemental stress measure, considered here to be the von Mises of the elements, and σlim is regarded as the upper limit of stress herein.Just to make it clear, the elemental von Mises in two dimensional examples can be formulated: where, σx and σy are the components of the stress tensor σ in the x and y directions, respectively, and τxy denotes the shear stress.
The stress tensor can be described: And attained by where D denotes the constitutive matrix shown in Eq. ( 6), B indicates the derivative matrix for shape function shown in Eq. ( 7), and u represents the element displacement vector shown in Eq. ( 8).
where v is the Poisson's ratio, and L is the length of the side of the discrete element (square element).More information on the theoretical derivation of the stress tensor and the finite element (FE) analysis can be found in the article opened by Biyikli and To [11].

Minimum volume fraction problem under compliance and stress constraints
As with the optimization model described above, single-constraint problems in topology optimization have been studied by many scholars, and few literatures have reported on multi-constraint problems (both stress and compliance are used as constraints in this paper).Consequently, a mathematical model for proportional topology optimization that minimizes volume fraction while satisfying compliance and stress constraints is proposed in this paper.And it can be constructed as follows: , if 0 01 where, besides the nomenclature given for Eqs. ( 1) and ( 2), Clim denotes the compliance limit, which can usually be set for different engineering problems.

The modified proportional topology optimization algorithm
The essential content and implementing process of the modified proportional topology optimization algorithm are briefly introduced from the following aspects.

Material interpolation scheme
In terms of topology optimization, choosing an appropriate material interpolation model is the necessary guarantee for the precise solution of the optimization model.To tackle the topology optimization problem for continuum structure with respect to 0-1 discrete variable, the SIMP method have been extensively applied [24,25].As a separate note, the modified SIMP approach is applied in PTO method [11].We use a new material interpolation scheme derived from the improved SIMP approach, offered by Nie [26] for the gradient algorithm of topology optimization, which is written as: where Emin indicates the Young's modulus, which is set an appropriate value to void elements for the sake of warding off the singular stiffness matrix, E0 represents the initial Young's modulus, and p (p≥1) is called penalty factor.In some literature [27,28], p = 3, but the value of p is not invariable, which depends on the specific physical problems.To be clear, although the presented algorithm in this paper belongs to the non-gradient algorithm, that is, it does not need to calculate the derivative of the design variable.The material interpolation model can also be applied to the gradient algorithm.The first partial derivative of E concerning xi reads: ) where there is a constant 1/e p in Eq. (11).The existence of 1/e p can significantly decrease the number of elements with derivative value of zero, which correspondingly decreases the number of elements with elastic modulus equal to zero, that is, compared with the traditional SIMP model, the possibility of optimization is increased, and the optimization result may be better [26].

Density filter
On account of the characteristic of non-sensitivity in MPTO algorithms, it is very important to take some measures in the process of optimization to preclude the numerical instability phenomena (i.e., mesh dependence and checkboard).Referring to the literatures [11,29], the density filtering technique is applied to MPTO algorithms, the core of which is to replace the initial central elemental density with the weighted average of all elemental densities in the filtering radius.It is given: where x1i indicates the density of element i after filtering calculation, Ni denotes the neighbourhood assembly of elements made up of the elements j within the radius rmin of element i, measured from their respective geometry center [18].ωi,j represents the weight coefficient, which can be acquired by: , max(0, ( , )) where Δ(i, j) is the distance between elements i and j, which is measured from center to center of the elements.

Heaviside threshold function
In order to obtain the 0/1 solution more easily and efficiently, the Heaviside threshold function is fully utilized to improve the optimization results in topology optimization.For example, Xu et al. [30] proposed a volume-preserving density filter in the light of the Heaviside threshold function, which makes the value of more elemental densities tend to 0 or 1. Li and Khandelwal [31] obtained binary discrete topologies by using the continuation scheme of material penalty coefficient and filter parameters based on the scheme of volume preserving Heaviside filter.Fu et al. [32] discussed in detail the advantages of applying the Heaviside smoothing function over the Heaviside step function in the topology optimization algorithm.To improve the optimization results for the concerning three problems in this paper, the Heaviside threshold function is used, which is written as: where x2i indicates the physical density calculated by the Heaviside filter, β represents a scaling parameter that controls the rate of approximation, η is a threshold value (typically 0.5).In the optimization process, this scheme plays a significant role in decreasing the intermediate densities: For the Heaviside filter, β increases gradually during the iterations.It can be seen from Eq. ( 15) that the density variable can be completely equal to 0 or 1 when β increases to infinity.

Procedure for modified proportional topology optimization algorithm
The execution procedure for MPTOc is presented below: Step 1. Define the corresponding parameters of the algorithm, such as the upper bound xmax and lower bound xmin for the density variable, proportion exponent (q), history coefficient (α), filter radius rmin, penalty coefficient p, and Young's modulus E0 and Emin.
Step 2. Set up vectors and matrices prepared for FE analysis and density filtering, and calculate the elemental weight factor based on Eq. ( 12).
Step 3. Perform the FE analysis, compliance and stress calculation with Eqs. ( 1), ( 3) and (10).Judge whether the main loop is terminated according to the termination condition; if it is satisfied, jump out of the loop and output the optimization results, or else execute step 4.
Step 4. Initialize the remaining material amounts (VRM) based on the target material amounts (VTM) by the volume constraint.Calculate the proportion of structural compliance Cpro by: where vj is elemental volume and Ci indicates the elemental compliance.
Step 5. Enter the internal loop.The VRM is assigned to elements proportionally according to their objective function values, which can be expressed: where x pro i denotes the elemental density after performing the proportion calculation.
After that, filter the elements based on density filter with Eqs. ( 12) and ( 13), and the value of density variable is limited to an interval between the upper and lower boundaries.Then, calculate the filtered actual material amounts (VAM), and update the VRM obtained by subtracting the actual material amounts from the target material amounts.When the termination condition (the value of VRM is less than a prescribed tolerance) of the internal loop is met, that is to say, the VRM is almost allocated, output the density variable x pro i and carry out step 6.If no, execute step 5 again.
Step 6. Exit the internal loop and update the density variable, which can be obtained by: ( 1) where, x new i denotes the new density of element i after executing the internal loop in the algorithm, x pre i indicates the density of element i from the previous iteration, and x pro i indicates the optimized elemental density in the current iteration, α represents the history coefficient.
For a more comprehensive implementation of the MPTOc algorithm, readers can refer to the original PTO algorithm in the literature [11].As a note, the main difference between NPTO and PTO is that NPTO adopts an interpolation scheme different from PTO, and the Heaviside threshold function is added on this basis.
The general procedure for MPTOs follows the MPTOc.The main differences between MPTOs and MPTOc algorithms are as follows: MPTOs adds a parameter (i.e., the stress limit σlim) in the step 1; in the step 3, replace Eq. ( 1) with Eq. (2); in the step 4 and 5, the value of TM is calculated by: , lim 1 , lim 1 0.001 , for max( ) > 0.001 , for max( ) where  , denotes the von Mises of element i and N indicates the total number of elements.And MPTOs algorithm distributes the RM proportionally based on the proportion of structural stress σpro instead of the proportion of structural compliance Cpro.The distribution equation read: The implementation process of the MPTOm algorithm mainly follows the MPTOs algorithm, except for the following points: add the compliance limit Clim as an input parameter in the step 1; in the step 3, replace Eq. (2) with Eq. (9).

Numerical examples
In this part, the algorithms for MPTOc, MPTOs and MPTOm concerning the MBB beam are validated by numerical examples.Fig. 1 shows the dimensions, boundary conditions and external load of full design domain (Fig. 1,  a) and half design domain (Fig. 1, b) for MBB beam.In view of the symmetry of MBB beam structure in this optimization problem, we adopt the MBB beam model with half design domain in the numerical example.The aspect ratio in the design domain is set to 120:40, that is to say, the design domain is divided into 120×40 meshes, and the thickness of the domain is given 1.The main parameters used in this paper are the same as the literature published by Biyikli and To [11], for the convenience of calculation and comparison, which are provided as follows: the external load F = 1, Poisson's ratio v = 0.3, the Young's modulus Emin = 1.0×10 -9 , E0 = 1, and the number of loaded elements Ne = 3.The other parameters, the filter radius rmin, penalty factor p, control parameters (α and q), will be discussed in detail in section 4.1 and 4.2.The corresponding parameter units are as follows: the units are Newton for load, mm for length, MPa for stress, and Nmm for compliance.As a note, nelx, nely, and rmin are in units of element (e.g., nelx indicates the number of elements in the x direction), and the thickness of elements and element edge length are considered to be unity.The material volume fraction (V) represents the proportion of the final optimized solid material volume to the full solid material volume.In the topology optimization algorithm, the parameters of penalty factor p and filter radius rmin with different values will have a strong effect on the optimization results.Because the material interpolation scheme in MPTO algorithms is different from that in PTO algorithm, and the Heaviside threshold function is introduced in this algorithm, it is necessary to discuss the parameters p and rmin with different value so as to attain the optimal topology optimization results, that is, a clear topology, smaller compliance values and material volume fractions and excellent convergence performance can be obtained.It should be noted that the algorithm in this paper is an improvement and extension based on the PTO algorithm, and the computer programs partially inherits from the MATLAB code by Biyikli and To [11].Therefore, when discussing the parameters p and rmin, the values of q and α are assigned the values with the optimal performance in the PTOc and PTOs algorithms respectively, that is, in the MPTOc algorithm, α = 0.5, q = 1; In the MPTOs and MPTOm algorithms, α =0, q = 2.
After numerical examples verification, taking into account the objective function value and topological configuration, we can conclude that when p = 3, rmin = 3, the optimization result of MPTOc algorithm performs the best, when p = 5, rmin = 2.5, the MPTOm algorithm can obtain the optimal optimization results, and when p = 5, rmin = 2.5, MPTOm algorithm can attain the optimal optimization results.

Effectiveness of the MPTO algorithms based on control parameters (q and α)
In this section, the effectiveness of the MPTO algorithms considering the control parameters (proportion exponent q and history coefficient α) is validated.Based on the values of rmin and p obtained in Section 4.1 (i.e., rmin = 3, p = 3 in MPTOc algorithm; rmin = 3, p = 5 in MPTOs algorithm and rmin = 2.5, p = 5 in MPTOm algorithm), α values from 0 to 0.9 at 0.1 intervals and q values from 0.5 to 2.0 at 0.25 intervals.
Tables 1, 3, and 5 show the topology structures for the MBB beam gained by each algorithm with above parameters.Table 2 lists the compliances (C) and iteration numbers (Iter) of MPTOc algorithm, Table 4 represents the iteration numbers and material volume fraction (V) for MPTOs algorithm, and Table 6 displays the iteration numbers and material volume fraction for MPTOm algorithm.It is worth noting that the results marked with the symbol "/" are considered as invalid topology structures for non-convergent simulations, which run over 1000 iterations.Taking the topology structures and optimized data into account in Tables 1-6, the optimal topology optimization results are obtained when q = 1, α =0.5 in MPTOc algorithm and q = 2, α = 0 in MPTOs algorithm and MPTOm algorithm.To some extent, the optimization results fully verify the values of q and α in Section 4.1.

Comparison of MPTO to PTO
In this section, we compare the MPTO to PTO concerning the MBB beam exhibited in Fig. 1, b.The topology optimization problem is solved based on setting the following control parameters: α from 0 to 0.9 in increments 0.1, q from 0.5 to 2 in increments 0.25.
Table 7 demonstrates the final optimal results (contain iteration numbers, material volume fraction, compliance, max stress and topology structure) gained by each algorithm in accordance with their respective optimal parameter settings.Additionally, the stress distributions for PTOs, MTPOs and MPTOm are also presented in Table 7.It should be pointed out that the results of algorithms PTOc and PTOs are from the paper opened by Biyikli and To [11] while the data of algorithms MPTOc, MPTOs and MPTOm are derived from the new algorithms suggested in this article.And the values of the constraints are parenthesized.
As can be seen from Table 7, comparing MPTOc with PTOc, under the constraint of material volume fraction, the compliance of MPTOc (C = 261.27) is smaller (better) than that of PTOc (C = 266.61),and the iteration numbers of MPTOc (Iter = 36) is much less than that of PTOc (Iter = = 170), and the topology structure of MPTOc performs better than that of PTOc (there exists obvious gray-scale elements and redundancy in the topology structure of PTOc).Besides, the maximum stress value (1.08) of the structure obtained by the two algorithms is the same.
It should be clarified that MPTOm is obtained by adding compliance constraint and extending on the basis of MPTOs.Therefore, MPTOm and MPTOs are compared with PTOs in Table 7.It is observed that the three algorithms attain the same material volume fraction (0.31), and it is undeniable that the compliance of PTOs (C = 294.92) is slightly smaller (better) than that of MPTOs (C = 299.42).However, MPTOm can obtain the smallest (best) compliance (291.96) and the maximum stress (1.0), it should be noted here that although the compliance and stress values (Clim = 294.92,σlim =1.08) are taken as the constraints of MPTOm, the algorithm converges when the algorithm is not completely equal to the constraints value during the execution process.In addition, MPTOs and MPTOm have the same iteration numbers (Iter = 63), which are far smaller than that of PTOs (Iter = 206).Observing the stress nephograms of PTOs, MPTOs and MPTOm, the stress distributions in most areas of them are similar, and they are all within the allowable stress.Whereas there is no redundancy and gray-scale elements in the optimized structure for MPTOs and MPTOm, which is greatly superior to the optimized structure for PTOs.
In summary, based on the PTO algorithm, because the interpolation function is modified and the Heaviside threshold function is added in the MPTO algorithms, the convergence speed of MPTO is much faster than PTO, and the topology structure of MPTO also has prominent advantages over that of PTO with redundancy and gray-scale elements.Thus, the MPTO algorithms can obtain better optimization results than PTO.
In order to better check the convergence performance of MPTO algorithms, their performance is discussed in detail as the control parameters (proportion exponent q and history coefficient α) take different values.
Table 1 Topology structures for the MBB beam gained by MPTOc considering the parameters: rmin = 3, p = 3, q = 0.5−2, α = 0−0.9q \ Algorithm \ α 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Table 2 Compliances and iteration numbers for the MBB beam gained by MPTOc considering the parameters: rmin = 3, p = 3, q = 0.5−2, α = 0−0.9Figs. 2 and 3 show the topology optimization results (contain iteration numbers and topology structures with the best performance of each algorithm in terms of comprehensive indicators) attained by PTO (including PTOc and PTOs) and MPTO (including MPTOc, MPTOs and MPTOm).Examining Fig. 2, a-g, it is observed that the iteration numbers of algorithms PTOc and MPTOc will change with different values of the parameter α when parameter q takes a fixed value.Obviously, when both the algorithms MPTOc and PTOc converge, the iteration numbers of MPTOc are much lower than that of PTOc.In other words, the convergence speed of MPTOc is faster than that of PTOc.Furthermore, the topology structures of PTOc are inferior to that of MTPOc due to the redundancy attached to PTOc.
Fig. 3, a-g displayed that the iteration numbers of algorithms PTOs, MPTOs and MPTOm are also subjected to α for a fixed value of q.On the whole, the average number of iterations of the MPTOs is the least among the three algorithms, that is to say, the MPTOs has the fastest convergence speed, followed by MPTOm, and PTOs is the slowest.Optimized results for the MBB beam acquired by PTOc and MPTOc: a -q = 0.5, b -q = 0.75, c -q = 1, d -q = 1.25, e -q = 1.5, f -q = 1.75, g -q = 2 a b c d e f g Fig. 3 Optimized results for the MBB beam gained by PTOs, MPTOs and MPTOm: a -q = 0.5, b -q = 0.75, c -q = 1, d -q = 1.25, e -q = 1.5, f -q = 1.75, g -q = 2 In addition, it can be seen from the change curve that except when q = 0.5, the iteration numbers for PTOc increase as the increase of α while MPTOs and MPTOm change relatively gently.In terms of topology structures of these algorithms, when q = 0.5−1, there are more or less redundancy in the topology structures of the three algorithms; when q = = 1.25−2, the topology structures of MPTOs and MPTOm have notable advantages over those of PTOs which have some redundancy in most cases.

Conclusions
By modifying the interpolation function in the original PTO algorithm and adding the Heaviside threshold function, three algorithms are proposed in this paper, which are named MPTOc, MPTOs and MPTOm respectively.And the new algorithms are displayed by solving the topology optimization problems concerning MBB beam (minimum compliance problem under volume constraint, minimum volume fraction problem under stress constraint, and minimum volume fraction problem under compliance and stress constraints) and compared with PTO algorithm.Analysing and discussing the effectiveness and superiority of each algorithm based on control parameters (q and α), we can draw some conclusions as follows: These new algorithms (MPTOc, MPTOs and MPTOm) are able to solve their respective topology optimization problems and obtain effective topology structures for MBB beam.
MPTOc, MPTOs and MPTOm have obvious advantages over algorithms (PTOc and PTOs) concerning convergence efficiency (the convergence speed of MPTOc is almost 5 times that of PTOc, and the convergence speed of MTPOs and MPTOm is about 3 times that of PTOs) and acquiring clear topology structure without redundancy.In particular, MPTOc converges the fastest among these algorithms.On the other hand, MPTOs and MPTOm have similar performance in many aspects, such as iteration numbers, material volume fraction and topology structure.
The combined action of the Heaviside threshold function and the modified interpolation function in the new algorithms can effectively enhance the ability of the algorithms to attain better objective function values (e.g., the compliance value of MPTOc has been improved by approximately 2% compared to the compliance value of PTOc), better topology structure (more clear topology structure without redundancy) and noticeably improve the convergence efficiency of the algorithms.
The topology optimization results obtained by MPTOc, MPTOs and MPTOm have certain dependence on control parameters (q and α).To obtain better convergence performance and better topology structure, the parameter settings are recommended as following: q = 1, α = 0.5 in MPTOc and q = 2, α = 0 in MPTOs and MPTOm.

Fig. 1
Fig. 1 Dimensions, boundary conditions and external load for the MBB beam: a -full design domain, b -half design domain 4.1.Validation of the MPTO algorithms considering different values of filtering radius and penalty factor

Table 7
The final optimal results gained from each algorithm according to their respective optimal parameters