An SQP Algorithm for Structural Topology Optimization Based on Majorization–Minimization Method

: When applying the sequential quadratic programming (SQP) algorithm to topology optimization, using the quasi-Newton methods or calculating the Hessian matrix directly will result in a considerable amount of calculation, making it computationally infeasible when the number of optimization variables is large. To solve the above problems, this paper creatively proposes a method for calculating the approximate Hessian matrix for structural topology optimization with minimum compliance problems. Then, the second-order Taylor expansion transforms the original problem into a series of separable and easy-to-solve convex quadratic programming (QP) subproblems. Finally, the quadratic programming optimality criteria (QPOC) method and the QP solver of MATLAB are used to solve the subproblems. Compared with other sequential quadratic programming methods, the advantage of the proposed method is that the Hessian matrix is diagonally positive deﬁnite and its calculation is simple. Numerical experiments on an MBB beam and cantilever beam verify the feasibility and efﬁciency of the proposed method.


Introduction
Structural topology optimization aims to optimize the material layout in a given design space under a given set of loads, boundary conditions, and constraints to achieve the optimal performance of the structure. At present, structural topology optimization has been widely used in related technical fields, such as aerospace [1], 3D print [2], and reliability study of structures [3]. Since Bendsøe and Kikuchi conducted seminal work [4] on numerical topology optimization in 1988, more and more structural topology optimization methods have been proposed, such as density-based method [5], level set method (LSM) [6], moving morphable components (MMC) [7], and bidirectional evolutionary structural optimization (BESO) [8]. Among them, the density-based method has received wide attention due to its strong adaptability and simplicity of implementation. Accordingly, a series of mature solution algorithms have been developed. These algorithms can be divided into two main categories: one is the optimality criteria (OC) [9] algorithm and its extensions, such as the PDOC [10] method based on the idea of proportional differential control. The advantages of this type of method are the high computational efficiency, and the speed of the solving is less sensitive to the number of design variables. However, this type of method can be used for single-constraint optimization problems, while it is difficult to be extended to multi-constraint optimization problems. Another class of algorithms to solve the densitybased method is represented by the approximation algorithms. The main idea of this class of methods is to transform the original problem into a simpler subproblem, and then use the solution of the subproblem as the initial point for the next iteration until convergence. The main approximation algorithms are sequential linear programming (SLP) [11] SQP [12], the method of moving asymptote (MMA) [13], and its globally convergent counterpart (GCMMA) [14]. This type of method can solve multi-constraint problems, but its solving speed is slower than that of the previous methods.
The SQP algorithm has been one of the most successful general methods for solving nonlinear constrained optimization problems. The solving of nonlinear problems with SQP is divided into two main steps. The first step is calculating the Hessian matrix to obtain the subproblem, and the second step is solving this subproblem. Compared with other nonlinear programming problems, QP problems are more attractive. Etman et al. [15] transformed the sequential convex programming subproblem into a Newton Lagrange QP subproblem, and Rong et al. [16] directly used the Hessian matrix in the MMA subproblem to construct the quadratic programming subproblem. In practice, the Hessian matrix of most optimization problems is difficult to calculate. The SQP algorithms generally use quasi-Newton [17] methods to approximate the Hessian matrix to obtain the QP subproblem. In topology optimization, TopSQP [18] uses information from part of the second-order derivatives to ensure the positive definiteness of the Hessian matrix.
Albeit successful in solving nonlinear optimization problems, the SQP method is not as widely used as the MMA and OC methods in structural topology optimization. The main reason is that the calculation of the Hessian matrix may become a heavy burden for largescale topology optimization problems. On the one hand, the existence of the equilibrium equation makes it difficult to calculate the exact Hessian matrix, which may also not be positive definite; on the other hand, using the quasi-Newton method to approximate the Hessian matrix will lead to a massive amount of computation and storage. Facing these obstacles, using the problem-driven majorization-minimization (MM) [19] algorithm to compute the approximate Hessian matrix would be a good breakthrough point. This paper proposes a method to calculate the approximate Hessian matrix according to the MM algorithm, which takes advantage of the problem structure. Then, the original problem is expanded by second-order Taylor expansion to form a QP subproblem. Finally, to verify the correctness of the approximate Hessian matrix method, this paper uses the proposed QPOC method and the QP solver of MATLAB to solve the subproblem. As the numerical experiments show, the proposed method in this work has certain advantages in terms of convergence rate, cost time, and optimal solution compared with the OC and MMA methods.
The paper is organized as follows: Section 2 describes the MM algorithm and how to calculate the approximate Hessian matrix. The algorithms to solve the subproblem are presented in Section 3. Section 4 provides two benchmark problems to validate the effectiveness of the proposed method. Concluding remarks are provided in Section 5.

Optimization Model
To produce an almost void-and-solid design, the density-based method uses the solid isotropic material with penalization (SIMP) [20] material interpolation model to penalize elements with intermediate densities. The relationship between Young's modulus and density is defined as where E 0 is the stiffness of the material, x i is the design variable, E min is a very small positive number to avoid singularity of the stiffness matrix, and p is a penalization factor introduced to ensure black-and-white solutions. The mathematical formulation of the minimum compliance problem can be expressed as where C is compliance, F and U are the force and global displacement vectors, respectively, N is the number of elements used to discretize the design domain, K is the global stiffness matrix, u i is the element displacement vector, k 0 is the element stiffness matrix for an element with unit Young's modulus, V 0 and V(x) are the design domain volume and material volume, respectively, and f is the prescribed volume fraction. The derivative of the objective function is found as Deriving the derivatives of the equilibrium equation with respect to the design variables and manipulating it mathematically, one can write Substituting Equation (4) into Equation (3) and then generalizing the sensitivity of the compliance in relation to the "N" design variables yields In this way, the sensitivity of the objective function can be rewritten as The exact Hessian matrix of the objective function can be expressed as Here, the matrices A and B are both positive definite symmetric matrices; A is a diagonal matrix whose each term on the diagonal can be expressed as

Majorization-Minimization Algorithm
The MM algorithm is an algorithmic framework commonly used in machine learning and signal processing, and many algorithms can be interpreted by the MM algorithm framework, such as the EM [21] algorithm. The main idea of the MM algorithm is to find a surrogate function at the current iteration point and then use the extreme point of the surrogate function for the next iteration point, which is similar to the MMA method. The function g(x, x k ) is called the surrogate function of the original function f (x) if it satisfies the following two properties: where the x k is the current iteration point. The procedure of the MM is shown pictorially in Figure 1. Generally speaking, there are two principles for finding a surrogate function. One is that the surrogate function should approximate the original function as closely as possible to achieve faster convergence. The other is that the surrogate function should be as simple as possible to reduce the computational cost. For more methods of constructing surrogate functions, the reader is referred to [22]. The MM algorithm is locally convergent and converges to the stationary point, while its good results in practice have led to its wide application. function ( , ) k g x x is called the surrogate function of the original function () f x isfies the following two properties: where the k x is the current iteration point. The procedure of the MM is shown pi in Figure 1. Generally speaking, there are two principles for finding a surrogate f One is that the surrogate function should approximate the original function as cl possible to achieve faster convergence. The other is that the surrogate function sh as simple as possible to reduce the computational cost. For more methods of cons surrogate functions, the reader is referred to [22]. The MM algorithm is locally con and converges to the stationary point, while its good results in practice have l wide application. Figure 1. MM algorithm procedure.

Quadratic Programming Subproblem
The method [22] of combining the MM algorithm with the SQP algorithm i an approximate Hessian matrix H , where − H B A . Since matrices A and positive definite, the most straightforward approach is to take = H B like the method, which can ensure that the matrix ( ) − − H B A is positive definite. Howe trix B needs to calculate the inverse of the global stiffness matrix, which is very sive for large-scale topology optimization problems and does not meet the requi of low computational cost. Since the matrix A is diagonally positive definite, the exist a curvature parameter ( 0) qq , such that q − A B A . Then, the surrogate of Equation (2) can be expressed as The solution of Equation (10) is given by In this paper, we do not directly calculate the value of q , but use numerical expe to verify the effect of different q on the convergence result. The QP subproblem

Quadratic Programming Subproblem
The method [22] of combining the MM algorithm with the SQP algorithm is to find an approximate Hessian matrix H, where H B − A. Since matrices A and B are positive definite, the most straightforward approach is to take H = B like the TopSQP method, which can ensure that the matrix H − (B − A) is positive definite. However, matrix B needs to calculate the inverse of the global stiffness matrix, which is very expensive for large-scale topology optimization problems and does not meet the requirements of low computational cost. Since the matrix A is diagonally positive definite, there must exist a curvature parameter q(q > 0), such that qA B − A. Then, the surrogate function of Equation (2) can be expressed as The solution of Equation (10) is given by Appl. Sci. 2022, 12, 6304 5 of 13 In this paper, we do not directly calculate the value of q, but use numerical experiments to verify the effect of different q on the convergence result. The QP subproblem of Equation (2) is as follows: where where move is the moving limit, which is usually set in topology optimization to ensure the stability of the solutions. The surrogate function has several properties desired by the literature [22]: convex separable and the existence of a closed-form minimizer. The design of efficient and neat solution algorithms to verify the feasibility of the surrogate function is presented in the next section.

QPOC Method
Since the subproblem in Equation (12) is convex and separable, the algorithm for solving it should take full advantage of this property. First, before solving Equation (12), we construct an easy to solve quadratic programming problem, which is where λ is a scalar parameter. Because A is a diagonal positive definite matrix, Equation (14) can be decomposed into N problems; the optimal solution of Φ(x; λ) is Constrained by upper and lower bounds, the solution of Equation (14) is given by Then, we solve the univariate nonlinear equation, According to Equations (15) and (16), x i (λ) and r(λ) are monotonically decreasing continuous functions. Note that r(0) > 0, r(∞) < 0; hence, Equation (17) must have a unique solution. The literature [23] has a detailed description of the solution of single-dimensional functions and gives several efficient solutions, such as the dichotomous method, secant method, and brent method. To better combine the 88 lines of code (OC) [24], the dichotomous method was selected in this paper. Because the KT conditions of Equation (14) and Equation (17) constitute the KT conditions of Equation (12), once λ * is located such that Equation (17) holds, x(λ * ) is a global minimizer since it meets the KT conditions and convexity of Equation (12). The efficiency of the proposed method is mainly reflected in two aspects. One is the low computation and storage cost of the Hessian matrix. The other is to transform the N-dimensional QP subproblem into the root-finding problem of a one-dimensional function, which significantly accelerates the solving speed.
In fact, the update formula of the OC method is also explicitly expressed as Equation (18) also uses the dichotomous method to solve the Lagrange multipliers, which means that reproducing the above algorithm requires only initializing q and modifying the update formula for the design variables in 88 lines of code. In summary, the proposed algorithm is more like the application of the OC method in quadratic programming; therefore, it is called the quadratic programming optimality criteria (QPOC) method.

IPC Method
The QPOC method can only be used to solve single linearly constrained quadratic programs subject to lower and upper bounds. For topology optimization, the constraints are not only the volume constraint; fixed value constraints and perimeter constraints [25] can also be added. At present, the mathematical community has systematically studied the solution algorithm of multi-constraint quadratic programming, mainly including the interior point method [26], trust-region method [27], and gradient projection method [28]. To verify the feasibility of the QP subproblem for multiple constraints, we used the interiorpoint convex (IPC) algorithm that comes with the built-in function quadprog of MATLAB.

Numerical Examples
In this paper, two benchmark examples of an MBB beam and cantilever beam are used to demonstrate the performance of the proposed algorithm. The Young's modulus E 0 = 1, E min = 0, and penalty factor p = 3; to avoid the singular value phenomenon, x min = 0.001. Regarding the filters, the QPOC, OC, and MMA methods use sensitivity filters, while the IPC method uses density filters that do not use chain rule, yielding better results. As for the convergence criterion, QPOC uses the same convergence criterion as the OC method. For MMA and IPC, due to the large change of design variables, it is necessary to set new convergence criteria mean(sum(abs(x k+1 − x k ))) < 4 × 10 −4 and iterations < 200. To measure whether an optimized design has converged to a discrete solution, a measurement index [29] is introduced.
If all elements' densities are equal to 1 or 0, the optimized design is fully discrete, and M nd is 0%. A smaller value of M nd means that there are fewer intermediate densities in the structure, and the diagram of the optimized structure will be clearer. In addition, all the experiments were executed on a personal computer with an Intel CORE i7 11700 processor, Windows 10 (64 bit), and MATLAB R2021a.

MBB Beam
The MBB beam is a typical benchmark example used in topology optimization. In accordance with [24], the design domain, the boundary conditions, and the external load for the MBB beam are shown in Figure 2. The design domain was discretized with 120 × 40 square elements, and the volume fraction f was set to 0.5. The optimized designs obtained using the OC, IPC, and QPOC methods are illustrated in Figure 3.

MBB Beam
The MBB beam is a typical benchmark example used in topology optimization. In accordance with [24], the design domain, the boundary conditions, and the external load for the MBB beam are shown in Figure 2. The design domain was discretized with 120 × 40 square elements, and the volume fraction f was set to 0.5. The optimized designs obtained using the OC, IPC, and QPOC methods are illustrated in Figure 3.  An excellent solution algorithm is mainly reflected in the optimal solution, convergence rate, and calculation cost. The advantages and disadvantages of the method proposed in this paper was compared with the OC method from these aspects.

Optimal Solution and Optimized Design
It can be seen in Figure 3 that, by selecting appropriate curvature parameters and filter radius, the QPOC method and IPC method could obtain the optimized structures similar to the OC method, which verifies the feasibility of the proposed method. The optimized structures obtained using the QPOC method and OC method had the phenomenon of fuzzy boundary, while the optimized structures (Figure 3h,i) obtained using the IPC method were clearer, and the structure contour was more distinct.
The curvature parameter and the filter radius have a crucial influence on the structure, and it is necessary to explore the optimal combination of the two parameters. Because QPOC is a derivative algorithm of the OC method, its filter radius should be chosen consistently with the OC method, which is also verified by the fact that Figure 3b-d exhibits The MBB beam is a typical benchmark example used in topology optimization. In accordance with [24], the design domain, the boundary conditions, and the external load for the MBB beam are shown in Figure 2. The design domain was discretized with 120 × 40 square elements, and the volume fraction f was set to 0.5. The optimized designs obtained using the OC, IPC, and QPOC methods are illustrated in Figure 3. An excellent solution algorithm is mainly reflected in the optimal solution, convergence rate, and calculation cost. The advantages and disadvantages of the method proposed in this paper was compared with the OC method from these aspects.

Optimal Solution and Optimized Design
It can be seen in Figure 3 that, by selecting appropriate curvature parameters and filter radius, the QPOC method and IPC method could obtain the optimized structures similar to the OC method, which verifies the feasibility of the proposed method. The optimized structures obtained using the QPOC method and OC method had the phenomenon of fuzzy boundary, while the optimized structures (Figure 3h,i) obtained using the IPC method were clearer, and the structure contour was more distinct.
The curvature parameter and the filter radius have a crucial influence on the structure, and it is necessary to explore the optimal combination of the two parameters. Because QPOC is a derivative algorithm of the OC method, its filter radius should be chosen consistently with the OC method, which is also verified by the fact that Figure 3b-d exhibits An excellent solution algorithm is mainly reflected in the optimal solution, convergence rate, and calculation cost. The advantages and disadvantages of the method proposed in this paper was compared with the OC method from these aspects.

Optimal Solution and Optimized Design
It can be seen in Figure 3 that, by selecting appropriate curvature parameters and filter radius, the QPOC method and IPC method could obtain the optimized structures similar to the OC method, which verifies the feasibility of the proposed method. The optimized structures obtained using the QPOC method and OC method had the phenomenon of fuzzy boundary, while the optimized structures (Figure 3h,i) obtained using the IPC method were clearer, and the structure contour was more distinct.
The curvature parameter and the filter radius have a crucial influence on the structure, and it is necessary to explore the optimal combination of the two parameters. Because QPOC is a derivative algorithm of the OC method, its filter radius should be chosen consistently with the OC method, which is also verified by the fact that Figure 3b-d exhibits a similar structure with the same filter radius as the OC method. For the IPC method, reducing the filter radius can obtain a clearer structure as shown in Figure 3g,h, but too small a filter radius will lead to fine branches in the optimized structure and affect the fabrication, as shown in Figure 3i. It can be seen from Figure 3e-g that when the curvature parameter was increased, the filter radius could be appropriately reduced to obtain a clearer structure without fine branches.
To further analyze the influence of the QPOC and IPC methods and parameters on the optimized structure, the compliance and discrete index M nd were calculated, and the results are shown in Table 1. As can be seen from Table 1, the values of M nd for the QPOC and OC methods were comparable, while the M nd for the IPC method was much smaller than that of the other two methods. When q = 2 and r = 2, the M nd of the IPC method was the smallest, 8.5%, approximately 72% smaller than that of the OC method, which indicates that the IPC method could significantly reduce the intermediate densities and reduce the boundary-blurring effect. In addition, the compliance of QPOC and OC methods was approximately equal; the compliance obtained by the IPC method was smaller than that of the previous two methods. It is worth noting that the compliance of the IPC (q = 2, r = 2.4) method was about 12% smaller than that of the OC method. Because the optimized structures of both had similar geometries, the reason for this phenomenon was that the M nd of the IPC method was smaller, which means that the optimized structure had fewer intermediate densities, while SIMP penalized the stiffness of elements with intermediate density, resulting in a decrease in the integral structural stiffness.
In summary, the optimized structure and compliance of the QPOC method were consistent with those of the OC method, and the selection of the filter radius should be the same as that of the OC method. The optimized structure of the IPC was clearer, with a distinct boundary and smaller compliance. The selection of filter radius was negatively correlated with curvature parameters. Figure 4a shows the convergence curves of the compliance under different methods and parameters. For the QPOC method, when the curvature parameter decreased, the corresponding purple (q = 3), red (q = 2), and black (q = 1) convergence curves were below the previous curve, and the same phenomenon was also observed for the IPC method. This indicates that the curvature parameter was negatively related to the convergence rate. However, this does not mean that the curvature parameter can be arbitrarily small; otherwise, it may lead to the surrogate function being smaller than the original function resulting in non-convergence of the algorithm. Figure 4b illustrates that the curvature parameter and the average iteration step of the design variables were negatively correlated, and a larger average iteration step could enhance its capability to find the extreme points and improve the convergence rate. It is worth noting that Figure 4b also indirectly illustrates that the selection of the filter radius in the IPC method was positively correlated with the average iteration step of the design variables.

Convergence Rate
in the IPC method was positively correlated with the average iteration step of the design variables.
Overall, the IPC ( 1 q = ) method had the fastest convergence rate, and the OC and QPOC ( 1 q = ) converged at about the same rate. The compliance of all the methods rapidly decreased in the first few steps, and then the rate of decrease slowed down, with the compliance of all methods reaching convergence in fewer than 60 steps.

Calculation Cost
As can be seen from Table 1, the single-step computation time of OC and QPOC was the same, about one-fifth of the single-step computation time of the IPC method. The reason for the basic agreement between the two computational costs is that both update formulas can be expressed explicitly and both use the dichotomous method to compute Lagrange multipliers. The QPOC method transforms the QP subproblem into a root-finding problem for a single-dimensional function, which is much simpler than solving a system of nonlinear equations using the Newton-Raphson method in the IPC method. This is the reason why the computational cost of QPOC is much cheaper than that of the IPC method.
The comparison between different methods of compliance and calculation time and discrete indices can be found in Figure 5, where the data in the figure are all normalized. In general, the OC and QPOC methods had an advantage over the IPC method in terms of computation time, and the QPOC could be adjusted by the curvature parameter ( 2 q = ) such that its computation time was less than that of the OC method. IPC method could obtain better optimized structures, and it could be used to solve multi-constrained problems. Overall, the IPC (q = 1) method had the fastest convergence rate, and the OC and QPOC (q = 1) converged at about the same rate. The compliance of all the methods rapidly decreased in the first few steps, and then the rate of decrease slowed down, with the compliance of all methods reaching convergence in fewer than 60 steps.

Calculation Cost
As can be seen from Table 1, the single-step computation time of OC and QPOC was the same, about one-fifth of the single-step computation time of the IPC method. The reason for the basic agreement between the two computational costs is that both update formulas can be expressed explicitly and both use the dichotomous method to compute Lagrange multipliers. The QPOC method transforms the QP subproblem into a root-finding problem for a single-dimensional function, which is much simpler than solving a system of nonlinear equations using the Newton-Raphson method in the IPC method. This is the reason why the computational cost of QPOC is much cheaper than that of the IPC method.
The comparison between different methods of compliance and calculation time and discrete indices can be found in Figure 5, where the data in the figure are all normalized. In general, the OC and QPOC methods had an advantage over the IPC method in terms of computation time, and the QPOC could be adjusted by the curvature parameter (q = 2) such that its computation time was less than that of the OC method. IPC method could obtain better optimized structures, and it could be used to solve multi-constrained problems.

Cantilever Beam
To demonstrate the effectiveness of the surrogate function under multiple con straints, a cantilever beam problem is provided in this section. The design domain of the cantilever beam is illustrated in Figure 6. It was fully constrained at the left end and sub jected to a downward unit load at the right center. To set multiple constraints, a fixed

Cantilever Beam
To demonstrate the effectiveness of the surrogate function under multiple constraints, a cantilever beam problem is provided in this section. The design domain of the cantilever beam is illustrated in Figure 6. It was fully constrained at the left end and subjected to a downward unit load at the right center. To set multiple constraints, a fixed value constraint was added, and the virtual densities corresponding to the first four columns of finite elements in the design area were set to 1, expressed as follows: (20) Figure 5. Comparison of optimization indices obtained using different methods.

Cantilever Beam
To demonstrate the effectiveness of the surrogate function under multiple constraints, a cantilever beam problem is provided in this section. The design domain of the cantilever beam is illustrated in Figure 6. It was fully constrained at the left end and subjected to a downward unit load at the right center. To set multiple constraints, a fixed value constraint was added, and the virtual densities corresponding to the first four columns of finite elements in the design area were set to 1, expressed as follows: In terms of the other parameters, 1, 2, 0.5 The design domain was discretized with 60 × 30, 80 × 40, and 120 × 60 square elements. Both the IPC and MMA methods were applied to this problem, and the optimization results are shown in Figure 7.   In terms of the other parameters, move = 1, q = 2, f = 0.5. The design domain was discretized with 60 × 30, 80 × 40, and 120 × 60 square elements. Both the IPC and MMA methods were applied to this problem, and the optimization results are shown in Figure 7.

Cantilever Beam
To demonstrate the effectiveness of the surrogate function under multiple constraints, a cantilever beam problem is provided in this section. The design domain of the cantilever beam is illustrated in Figure 6. It was fully constrained at the left end and subjected to a downward unit load at the right center. To set multiple constraints, a fixed value constraint was added, and the virtual densities corresponding to the first four columns of finite elements in the design area were set to 1, expressed as follows: In terms of the other parameters, 1, 2, 0.5 The design domain was discretized with 60 × 30, 80 × 40, and 120 × 60 square elements. Both the IPC and MMA methods were applied to this problem, and the optimization results are shown in Figure 7.    As shown in Figure 7, some checkerboard and porous elements appeared around the structure optimized using the MMA method, and the same phenomenon occurred in the literature [30]. Compared with the MMA method, the optimized structure of the IPC method was clearer, and the boundary contours were more distinct, with no checkerboard phenomenon. These benefits suggest that the optimized structure of the IPC method was superior to that of the MMA method in terms of manufacturing. Table 2 shows the optimization data of both methods for different sizes of design domains. The IPC method performed better than the MMA method. In the design domain, from small to large, the compliance of the IPC method was about 18.8%, 18.9%, and 13% smaller than that of MMA, respectively, while the number of iterations was 76%, 76%, and 42% smaller, respectively, and the overall calculation time was about 60%, 74%, and 4.2% smaller, respectively. It is worth noting that, in Figure 8a, the compliance curve of the MMA method showed a peak value, while the volume ratio convergence curve of MMA in Figure 8b showed a corresponding valley value. The decrease in volume ratio was the reason for the sudden increase of compliance, which also indicates that the MMA method was not as capable of maintaining constraints as the IPC method. superior to that of the MMA method in terms of manufacturing. Table 2 shows the optimization data of both methods for different sizes of design domains. The IPC method performed better than the MMA method. In the design domain, from small to large, the compliance of the IPC method was about 18.8%, 18.9%, and 13% smaller than that of MMA, respectively, while the number of iterations was 76%, 76%, and 42% smaller, respectively, and the overall calculation time was about 60%, 74%, and 4.2% smaller, respectively. It is worth noting that, in Figure 8a, the compliance curve of the MMA method showed a peak value, while the volume ratio convergence curve of MMA in Figure 8b showed a corresponding valley value. The decrease in volume ratio was the reason for the sudden increase of compliance, which also indicates that the MMA method was not as capable of maintaining constraints as the IPC method.

Conclusions and Discussion
This work introduced the MM algorithm to topology optimization, providing a new perspective for the application of the SQP algorithm in topology optimization. Numerical experiments on an MBB beam and cantilever beam verified the feasibility and efficiency of the proposed method. The contribution of our method can be concluded as follows: 1. We proposed a novel way to calculate the Hessian matrix for structural topology optimization with minimum compliance problems. The Hessian matrix could be given almost simultaneously when solving the first derivative, and the calculation cost was almost negligible.

Conclusions and Discussion
This work introduced the MM algorithm to topology optimization, providing a new perspective for the application of the SQP algorithm in topology optimization. Numerical experiments on an MBB beam and cantilever beam verified the feasibility and efficiency of the proposed method. The contribution of our method can be concluded as follows:

1.
We proposed a novel way to calculate the Hessian matrix for structural topology optimization with minimum compliance problems. The Hessian matrix could be given almost simultaneously when solving the first derivative, and the calculation cost was almost negligible.

2.
The QP subproblem is formed by the second-order Taylor expansion of the original problem, which is convex separable and easy to be imported into the QP solver. The convergence rate can be controlled by the curvature parameter.

3.
To solve the QP subproblem, we proposed the QPOC method. The derivation of the update formula of this method is supported by rigorous mathematical theory, while the update formula of the OC method is heuristic.
Compared to using general-purpose solvers, a specialized algorithm devised according to the structure of the problem might be a better option.