Highly Efficient Feasible Direction Method (HEFDiM) for Structural Topology Optimization

01/07/2020 ∙ by Zhi Zeng, et al. ∙ Xidian University 0

Feasible Direction Method (FDM) is a concise yet rigorous mathematical method for structural topology optimization, which can be easily applied to different types of problems with less modification. In addition, the FDM always converges to a near optimum rapidly. However, the problem of inefficiency stays unsolved. In this work, we advance the state-of-the-art by proposing a highly efficient feasible direction method (HEFDiM), which substantially improves the efficiency of the FDM with negligible loss of accuracy. The proposed method can benefit us in at least four aspects: 1) Analytical gradient projection; 2) Fewer heuristics and clear physical meaning; 3) Negligible memory and time-cost for updating; 4) Directly applied to different problems without extra efforts. In particular, we address problems including 1) Efficient determination of effective constraints for gradient projection; 2) Acceleration of null-space projection calculation; 3) Avoidance of time costing 1-D searching; 4) Elimination of split-stepping and zig-zag problems. Benchmark problems, including the MBB, the force inverter mechanism, and the 3D cantilever beam are used to validate the effectiveness of the method. Specifically, the results show that the updating speed of the HEFDiM is approximately 10 times higher (even faster for larger-scale optimization problems) with lower objective value than that of the classical efficient 88-line MATLAB code (Andreassen et, al. 2011). The HEFDiM is implemented in MATLAB which is open-sourced for educational usage.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 5

page 7

page 8

page 10

Code Repositories

HEFDiM-preview-version

Highly Efficient Feasible Direction Method (HEFDiM) for Structural Topology Optimization (Preview Version)


view repo
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Since the introduction of topology optimization to structural design, it has been successfully applied to many different types of structural design problems through various optimization schemes including density methods Bendsoe1989a , boundary variation methods (like the level-set approach Allaire2002 ; Wang2003 ; Zhu2015 ), intelligent optimization methods (like the genetic evolutionary method Yang1999 ), etc. A comprehensive review can be found from Ref.  Sigmund2013 ; Deaton2014 ; bendsoe2009topology ; zhang2018topology .

The Solid Isotropic Material with Penalization for intermediate densities method (SIMP)Bendsoe1989a ; bendsoe2009topology

is considered the most effective material interpolation scheme thus has been widely implemented in industrial applications. The scheme is formulated as follows:

(1)

where is the Young s modulus,

is the vector of design variables which are constrained to

, is the stiffness of the material, is the penalization factor. Using this scheme, different approaches, including the Optimality Criteria (OC) method Yin2001 ; Ananiev2005 , the Method of Moving Asymptote (MMA) Svanberg1987

, the Sequential Linear Programming (SLP) method 

Sigmund1997 ; Fleury1989 , and the Feasible Direction Method (FDM) Vanderplaats1984 ; Vanderplaats1973 ; Chang2014 ; Tavakoli2012 , etc., have been proposed to solve the structural topology optimization problem and each approach has its own desirable characteristics.

The OC method is one of the most efficient and widely used method. It tends to convert the optimization problem to an equation solving problem using the KT condition, which is attacked by constructing a contraction mapping whose fix-point is the solution. While, the criteria (as well as the contraction mapping) for different problems may be derived case by case (for example, the contraction mapping for compliance minimization problems howell2001compliant is different from that for compliant mechanism design problems Sigmund2001 ). Especially, for more complicated cost functions, the heuristic updating scheme is difficult to design patnaik1995merits .

The MMA converts the implicit constraints to explicit ones with theoretical guarantee, and transforms the original problem into a series of localized strictly convex approximating subproblems which are solved by a dual method Svanberg1987 . Similarly, the SLP method involves sequentially solving a serials of approximated linear subproblems using the linearized objective and constraint functions Sigmund1997 ; Fleury1989 . These methods are also widely used in topology optimization, while they are usually less efficient comparing with the OC method.

The FDM, proposed by Zoutendijk Vajda1962 , was firstly introduced to structural topology optimization by Vanderplaats Vanderplaats1973 . By utilizing the gradient at the current point, the method can provide a feasible search direction via solving a linear programing problem at each iteration. The FDM is a concise yet rigorous mathematical method for structural topology optimization, which can be easily applied to different types of problems with less modification. In addition, FDM always converges to a near optimum rapidly. While, the efficiency of state-of-the-art cannot be guaranteed since a two-stage iteration is included and the simplex method used in the second stage has a polynomial complexity on average. As a promising special case of the FDM, the gradient projection method (GPM), proposed by RosenRosen1960 ; Rosen1961

, projects the gradient onto the tangent hyperplane in the constraint space. The original GPM was found to be inadequate for structural optimization since how to find out the projection is still a bottleneck. To alleviate this problem, a modified method called Direct Gradient Projection Method (DGPM)is proposed by Chang 

Chang2014 . By transforming the design variable, only one effective constraint exists. This method simplifies the optimization problem but at the cost of changing the constraints to nonlinear ones. Since a series of correction steps are required, tradeoffs between efficiency and reliability still have to be made.

As discussed above, although different attempts have been made to employ the FDM and its derivatives to structural topology optimization, the problem of inefficiency remains unsolved. After extensive investigations, the essential issues that hinder the efficiency are obtained and summarized as:

  • Effective constraints for gradient projection are not easy to be determined

  • Null-space projection is time-costing

  • Obtaining the optimal step via 1-D searching along the feasible direction is inefficient

  • Numerical problems, including split-stepping , zig-zag and numerical imprecision, reduce the convergence speed.

Aiming to address the above listed problems, the current study provides a highly efficient feasible direction method (HEFDiM), which substantially improve the efficiency of the FDM with negligible loss of accuracy. The scope of the current study is given as follows. The problem statement as well as the paradigm for attacking these problems are presented in Section 2. The HEFDiM is presented in Section 3. Section 4 provides three benchmark problems to validate the effectiveness of the proposed method. Concluding remarks are provided in Section 5.

2 Problem statements

One typical large class of topology optimization problems have the structure that an energy-like objective function is minimized over a feasible region defined by a second order elliptic PDE together with bilateral bound and a single equality constraint. In this paper, we focus on solving this class of optimal design problems which have a wide range of applications in engineering design tavakoli2010nonmonotone . The optimization problem can be mathematically expressed as:

(2)

where is the vector of design variables. is the number of elements used to discretize the design domain. is the objective function, and are the global displacements and the force vector respectively. is the global stiffness matrix. and are the material volume and design domain volume, respectively. Finally, is the prescribed volume fraction.

In this section, the FDM is briefly reviewed first. Then, essential problems that hinder its efficiency are investigated. Finally, the paradigm for attacking these problems are provided.

2.1 Feasible Direction Method (FDM)

Let be the design variable, and be a move direction of and no normalization is required here. The FDM seeks to update the design variable in an iterative way as follows:

(3)

subjecting to

(4)

in which . The superscript denotes the iteration number (Eq. (4) also holds for ). The constraints for form a convex polyhedron in . For the inequality constraint on , we say that it is effective if either or holds, and let denote the set of indices of all such effective constraints.

The move direction has to be USABLE, which means

(5)

Otherwise, the objective function would not decrease. The constraints on are dependent. That is, for , we have

(6)

and

(7)

The constraints for form a polyhedral convex cone (the feasible region) in . For the inequality constraint on , we say that it is effective if both (or ) and holds. The move direction is said to be FEASIBLE if . Note that the superscript is dropped from here in order to keep the conciseness of the content.

A good choice of the feasible direction usually minimizes the left part of Eq. (5), which can be obtained by solving a linear programming problem using the simplex method Sharma2014DeterminationOF . However, this approach is usually time-costing especially when large-scale problems are encountered. An promising way to solve this problem is to set as the projection of on . We denote this projection as .

Since is a special case of , the constraints given in Eq. (6-7) are also applicable to . We denote xp as the set of indices of all effective constraints on , and suppose that there are elements in xp. Each effective constraint corresponds to a hyper-plane in , whose normal vector can be written as

(8)

for Eq. (6), where and

(9)

for Eq. (7), respectively. By concatenating all these normal vectors, we obtain a matrix as follows:

(10)

We denote the bases of the orthogonal complement of the subspace spanned by as . Then, the following equation holds:

(11)

is also called the feasible steepest descendent direction.

At this stage, a natural question arises that how could we determine , which can be formulated as follows.

(12)

subjecting to

(13)

where is an arbitrary subset of . Actually, is the smallest subset that satisfying the constraint (13). Otherwise, one can always release a redundant constraint and make larger.

2.2 Problem Statements

Although the FDM introduced in the previous section has a concise form and rigorous in mathematics, it is seldomly used in large scale practical applications for the sake of its low efficiency. The essential issues that hinder the efficiency can be attributed to:

  1. Directly solving Eq. (12) is an NP hard problem

    The number of candidate is , which increases rapidly with respect to . i.e. when the dimension of the design variable is large, the effective constraints on the feasible steepest decedent direction is hard to determine.

  2. The calculation of and Eq. (11) is time and memory-costing

    The null space calculation involves matrix decomposition which is extremely inefficient when the scale of the problem is quite large. In addition, the matrix calculation would also be time and memory-costing when large scale problems are encountered.

  3. Obtaining the optimal step is inefficient

    Suppose the above problems are all addressed, we still have efficiency problems regarding the calculation of , which usually involves 1-D searching. Since we have to calculate the objective function for each iteration, this process is inefficient.

  4. Numerical problems, including split-stepping , zig-zag and numerical imprecision, also reduce the convergence speed

For the first problem, existing methods for solving Eq. (12) are equivalent to solving a linear programming problem using the simplex method Sharma2014DeterminationOF . Thus, it has a polynomial complexity on average. Instead of solving an optimization problem, we directly expand iteratively. As for the second problem, we propose an analytical solution to the null-space projection, since the general form of constraints used in topology optimization is special in mathematics. The details will be demonstrated in the next section.

3 Highly Efficient Feasible Direction Method (HEFDiM)

In this section, the detailed procedure of the HEFDiM is provided. Four essential issues that hinder the efficiency are addressed.

3.1 Effective Constraints Determination

A two-stage iterative strategy is proposed to efficiently determine the effective constraints for gradient projection. We reduce the searching space from to . This is reasonable because an effective constraint on the feasible steepest descendent direction is also effective on . Mathematically, this can be expressed as

(14)

First, a smaller subset inside of that containing is generated. We set as the initial set and expand iteratively according to Algorithm until the expanded set satisfies the constraint (13).

  Initially,, let and
  repeat
     
  • Let , where is the component of

  • set

  • Increase by one

  until  satisfy the constraint (13)
Algorithm 1

The iteration is indispensable because even the component of is not conflict with , it may still be an effective constraint con the projection. This situation is illustrated in Fig. 1.

Figure 1: Why should we use an iterative process to obtain

As shown in Fig. 1, , , and are three components of the design variable . , , and are hyper-planes represented by constraints given in Eq. (4). In particular, corresponds to the equality constraint. is the feasible region. is a move direction (being the gradient or not) which is in confliction with only. The projection of on is (the small white dot). Since is not effective on , one may project onto the intersection of and , resulting , which conflicts with . Therefore, both and are effective on .

Second, there may exists a few redundant constraints in . We can get rid of them using Algorithm 2 as follows:

  For each index in
  repeat
     If satisfies the constraint (13), remove from
  until endFinally  one gets
Algorithm 2

For most problems in practice, the procedure given in Algorithm 1 can be accomplished within a few steps. In addition, an observation of less than 1% redundant constraints existing after running Algorithm 1. Thus, the procedure given in Algorithm 2 can be omitted to improve the efficiency without loss of accuracy.

3.2 Efficient Null-space Projection

Once the effective constraints on are determined, the feasible steepest decedent direction can be easily calculated using Eq. (11). During the practice, we found that the calculation of is time-costing especially when the scale of the problem is quite large, since matrix decompositions are involved. In addition, the calculation of Eq. (11) is also time and memory-costing when large scale problems are encountered. In order to tackle these problems, an analytical expression of the projection is proposed, which can greatly reduce the calculation. Details are given as follows.

Note that given proper shuffling of the indices of the components of , the matrix can always be written in the following form.

(15)

then, can be expressed as follows.

(16)

where is the matrix consisting the bases of the orthogonal complement of the subspace spanned by . This can be easily verified by noting that and .

In order to calculate , we proposed an analytical expression as follows:

(17)

where

(18)

in which is the number of columns in . can be further decomposed into a low rank matrix and a spare one, which is

(19)

where

(20)

Hence, we have

(21)

The projection can be calculated by

(22)

where

(23)

Suppose there are elements in , then

(24)

Note that the proposed method gives an explicit analytical expression for the null space projection and no matrix multiplication exists. Therefore, for large scale problems, the proposed method is efficient both in memory and in calculation.

3.3 Updating the Design Variables

After we have obtained the gradient projection , the next step is to find out an optimal searching step , where is the maximum step, which is defined as:

(25)

in which

(26)

The above equation shows that a new updating iteration begins whenever a new constraint for an element becomes effective. Thus, the maximum step is small enough so that the objective function is monotonous in . Therefore, the optimum step is simply .

Once is determined, the design variable can be updated using Eq. (3). The updating procedure is repeated iteratively until the change of

is below a prescribed tolerance. It should be noted that during this process, numerical problems including split-stepping , zig-zag and numerical imprecision, may reduce the convergence speed. This is because the algorithm begins a new loop whenever a new constraint is encountered. Thus, the optimized structure would be binarized pixel by pixel. A binarization acceleration mechanism is proposed to alleviate this problem, in which we let gray pixels near either

or be exactly or . Other pixels are simply shifted to maintain the volume fraction. The motivation is similar to that of the ’constraint thickness’ in Ref. Vanderplaats1973 .

4 Numerical Examples

Three benchmark examples including an MBB beam in 2D (Minimum compliance design problem), a compliant force inverter mechanism in 2D (Compliant mechanism design problem) and a cantilever beam in 3D are provided in this section to demonstrate the performance of the HEFDiM. The OC method is also applied to the MBB beam and the 3D cantilever beam for the purpose of comparison. All the experiments are executed on a personal computer with an Intel CORE i9 9900X processor, 128 GB memory, Windows 10 (64-bit), and MATLAB R2019b.

4.1 MBB Beam

Figure 2: Design domain of the MBB beam

The benchmark problem of finding the optimal material distribution, associated with the MBB beam, in terms of minimum compliance, with a constraint on the total amount of material, is presented in this section as an example. In accordance with Ref. Andreassen2011 , the design domain, the boundary conditions, and the external load for the MBB beam are shown in Figure 2.

The design domain is discretized with , and square elements, respectively. The volume fraction is set to 0.5. The Young s modules and the penalization factor is . Both the HEFDiM and the OC method are implemented in MATLAB and applied to this problem (the OC method is implemented using Andreassen s 88 lines of code in MATLAB Andreassen2011 ). As suggested by the authors, density and sensitivity filters with a radius of = 2.4, 6 and 12, respectively, are used in the OC method to eliminate the numerical difficulties. Similarly, in our work, a Gaussian filter with a radius of = 1.1, 2.0 and 4.0 for sensitivity and = 0.55, 1.0 and 2.0 for density is utilized in the HEFDiM. It should be noted that the filter sizes for both methods can not be further reduced or local minima would emerge Sigmund2007 . The normalized compliance as the function of iterations for the two method are plotted in Figure 3.

Figure 3: Convergence curves for the MBB beam compliance minimization problem using two different methods. Snapshots for the HEFDiM are above the lines, and snapshots for the OC method are below the lines.

Several conclusions can be drawn from this example. First of all, both methods almost converge in less than 30 iterations, which shows that both approaches are able to efficiently locate a good design from a uniform grey starting guess.

Second, for the OC method, although a steep drop is observed at the initial several steps, the convergence speed slows down dramatically after that and becomes even slower after 20 iterations. As for the HEFDiM, the convergence speed is a bit slower at the first several steps while surpasses the OC and converges to a smaller normalized compliance value. The reasons can be attribute to the fact that a heuristic contraction mapping is used in the OC method. While this may cause convergence problem when the variance of the distribution of the gradient with respect to different locations in the design domain is small.

Finally, a larger number of grey elements exist in the layout of the OC method than that of the HEFDiM, which shows that the HEFDiM is more likely to converge to a binarized results.

The optimized designs of the MBB beam and corresponding compliance obtained using different methods with different refinement are illustrated in Figure 4. The results demonstrate that, given proper filter settings, different refinements do not lead to different topologies.

Figure 4: Optimized design of the MBB beam and corresponding compliance obtained with the variant of the 88-line code (bottom) and the HEFDiM (top). A mesh with elements (left), elements (middle), and elements (right) has been used

Attentions should be paid that the choice of is not restricted to the gradient. Any usable move direction can be adopted in the proposed framework. For example, if we use an adaptively clipped gradient instead of the raw gradient in the HEFDiM for the MBB experiment, as shown in Fig. 5, dramatic improvement for the convergence speed and the reduction of the number of iterations for convergence can be obtained. The convergence curves for the HEFDiM using the raw gradient and that for the OC method are also plotted in Fig.  5 for comparisons. The way of choosing a better candidate is beyond the scope of this paper and will be tackled in future works.

Figure 5: Convergence curves for the MBB beam compliance minimization problem using three different method ( elements). The solid line and the dashed line are identical to the corresponding lines in Fig 4. The line with scattered points together with the frames show the updating process using the adaptively clipped gradient. Detailed settings are given as and .

Like the OC method, the time costs for each iteration depend on both the FEA process and an updating process. The updating speed for both methods in each iteration is plotted in Fig. 6. The results show that the updating speed of the HEFDiM is 10 times faster than that of the OC method and even faster when the scale of the problem is getting larger. Note also that the updating time is negligible comparing with that of the FEA process.

Figure 6: Computation time of the updating process for each iteration of the two methods

4.2 Compliant Force Inverter Mechanism

Figure 7: Design domain and boundary conditions of the compliant force inverter mechanism

To demonstrate the effectiveness of the HEFDiM in dealing with non-self adjoint problems, a benchmark example, associated with the compliant force converter mechanism, is provided.

Fig. 7 depicts the design domain of a compliant force inverter mechanism with single input-output behavior. The inverter outputs a displacement in an opposite direction to the actuating force. The fixed boundaries and the external force are denoted in the figure.

The design domain is discretized with square elements. The volume fraction is set to 0.2. An artificial spring with stiffness is attached to the output port and a dummy load is applied at the output port, which is expected to produce a horizontal displacement to the left. The goal is to maximize the geometric advantage ( ) of the mechanisms. The objective function used in this paper is formulated as li2014topology :

(27)

where is the global stiffness matrix of the structure. and represent the displacement fields of structures when only or the dummy load is applied, respectively. The HEFDiM is used to solve this problem. Similarly, a Gaussian filter with a radius of for density and for sensitivity is utilized in the HEFDiM.

Figure 8: Convergence curves for the complaint force inverter mechanism design problem with the set as the objective function

The as the function of iterations for the force inverter mechanism is plotted in Fig. 8. The iteration history and the associated layouts are also plotted in the figure. We can find from Fig. 8 that the problem almost converges in less than 100 iterations, which shows that the HEFDiM is able to quickly locate a good design from a uniform grey starting guess. By changing the artificial stiffness , the optimized design are illustrated in Fig. 9. The optimized designs of the inverter mechanism with different refinement are also illustrated. The results demonstrate that different refinements of the solution does not lead to a different topology.

Figure 9: Optimum topologies for the inverter mechanism (a) = 0.01, mesh: (b) = 100, mesh: (c) = 0.01, mesh: (d) = 100, mesh: .

It should also be noted that the used algorithm is identical to those used in the MBB example. The only difference lies in the objective function, which means that the HEFDiM can be applied to different type of problems without modifications.

4.3 3D Cantilever Beam

Figure 10: Design domain of a 3D cantilever beam

A 3D cantilever beam problem is provided in this section. The design domain of the cantilever beam is illustrated in Fig. 10. The prismatic beam is fully constrained in one end and a unit distributed vertical load is applied downwards on the lower free edge of the other end. The dimension of the computational domain is . Both the OC method and the HEFDiM are used to solve this problem. The settings for the OC are identical to those used in Ref.liu2014efficient . The generalized compliance as a function of iteration for the two method are plotted in Fig. 11.

Figure 11: Convergence curves for the 3D cantilever beam using the OC and the HEFDiM. Snapshots for the HEFDiM are below the lines, and snapshots for the OC method are above the lines

We can find from the results that both methods converge rapidly at the first several steps. While, the OC method slows down quickly and converges to a relatively larger magnitude of generalized compliance, comparing to the results obtained by the HEFDiM.

Figure 12: Optimum topologies for the 3D cantilever beam with a mesh of . (a) the HEFDiM, (b) the OC method.

The optimized designs at different iterations for both methods are also illustrated in Fig. 11 for comparison. The figure shows that the two methods converge along different paths and to different final topologies. The final designs of the 3D cantilever beam using both methods are illustrated in Fig. 12. The topology shown in Fig. 12(a) is quite different from and more concise comparing with Fig. 12(b). The results show that the HEFDiM is more likely to locate an optimum result.

4.4 Additional Discussions

In the previous section, the HEFDiM has been validated to a high degree of efficiency and accuracy. During practical applications, one may find out that issues including checkerboard patterns, mesh-dependencies, and local minima may influence the results of topology optimization. However, these problems emerges mainly due to the finite-element discretization and can be eliminated using proper density or sensitivity filters. This has been investigated thoroughly by Sigmund Sigmund2007 thus are not covered in the current study.

The type of constraint is not restricted to the exact form indicated in this paper. Actually, the HEFDiM can be applied to any problems whose constraints can be transformed to those presented in this paper. The HEFDiM can be directly applied to different problems without extra efforts. Additionally, the HEFDiM provides the most adjacent feasible direction for an arbitrary move direction (not restrict to the gradient ) in the constraint space.

5 Conclusions

The feasible direction method has been tendentiously considered to be a rudimentary optimization scheme and inefficient when employed in structural topology optimization. But actually, as long as the efficiency problems were solved, FDM is a promising method since it has a rigorous mathematical base. Solving the efficiency problems is exactly the main contribution of the current study.

In this work, we proposed a highly efficient feasible direction method (HEFDiM) for structural topology optimization with its efficiency being substantially improved without loss of accuracy. Meanwhile, fewer heuristics and less hyper parameters are included, thus it is easy to be applied to different problems without extra efforts. Benchmark problems, including the MBB, the force inverter mechanism and the 3D cantilever beam, validate the effectiveness of the method.

The method is implemented in MATLAB and open sourced for educational usage. It is recommended that the readers try our code for a better understanding of the proposed method.

6 Replication of Results

The method proposed in this paper is implemented in MATLAB and open sourced on GitHub for educational usage https://github.com/zengzhi2015/HEFDiM-preview-version. For commercial usage, please contact the authors.

Acknowledgements.
The authors gratefully acknowledge the financial support from the National Natural Science Foundation of China under Grant No. 51805397 and No. 61805185 and the Open Fund of State Key Laboratory of Robotics and System (HIT) under Grant No. SKLRS-2019-KF-07. Instructive advices from Professor Guimin Chen are also acknowledged.

Conflict of interest

The authors declare that they have no conflict of interest.

References