{'fraction_non_alphanumeric': 0.04818515739168214, 'fraction_numerical': 0.015028239661913977, 'mean_word_length': 4.282599624452327, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On Topology Optimization and Canonical Duality Method

Topology optimization for general materials is correctly formulated as a bi-level knapsack problem, which is considered to be NP-hard in global optimization and computer science. By using canonical duality theory (CDT) developed by the author, the linear knapsack problem can be solved analytically to obtain global optimal solution at each design iteration. Both uniqueness, existence, and NP-hardness are discussed. The novel CDT method for general topology optimization is refined and tested by both 2-D and 3-D benchmark problems. Numerical results show that without using filter and any other artificial technique, the CDT method can produce exactly 0-1 optimal density distribution with almost no checkerboard pattern. Its performance and novelty are compared with the popular SIMP and BESO approaches. Additionally, some mathematical and conceptual mistakes in literature are explicitly pointed out. A brief review on the canonical duality theory for solving a unified problem in multi-scale nonconvex/discrete systems is given in Appendix.


Introduction
Topology optimization is a mathematical tool that optimizes material layout within a prescribed design domain in order to obtain the best structural performance under a given set of loads and geometric/physical constraints. Due to its broad applications, this tool has been studied extensively by both engineers and mathematicians for more than 40 years (see the comparative reviews [55,56]). Generally speaking, a topology optimization problem involves both continuous state variable (such as the deformation field u(x)) and the mechanically sound structural topology. Second, the Newton-type algorithms can be used only for convex minimization. For nonconvex problems, numerical results obtained by these algorithms depend sensitively on initial data and numerical precisions adopted. It was discovered by Gao and his co-workers that the global optimal solutions are usually nonsmooth not only for coupled optimal design problems (see Chapter 7, [18]), but also for general nonconvex variational problems [29]. These nonsmooth solutions can't be captured by any Newton-type algorithms. By the fact that the SIMP is not a mathematically correct penalty method, this most popular engineering approach can never produce exact integer solution for any given real-world problem. The existence of gray scale elements and appearance of checkerboards patterns are the SIMP's two major intrinsic problems [12,57]. Although the commercial code by the BESO can produce integer solutions, it was discovered recently [26] that this popular method is actually a direct approach for solving a knapsacktype problem and it is not a polynomial-time algorithm. This the reason why the BESO is computationally expensive and can be used only for small sized problems.
Duality approaches for topology optimization have been studied via the traditional Lagrangian duality theory [4,5,39,40,41,61]. However, the Lagrange multiplier method can be used mainly for convex problems with equality constraints [44]. For nonconvex problems, the Lagrangian L(x, y) is usually not a saddle function. By the fact that min x max y L(x, y) ≥ max y min x L(x, y), the Lagrangian duality theory produces a so-called duality gap at each iteration. In order to reduce this duality gap, much effort has been made by mathematicians during the past 30 years [28,48]. For inequality constraints, both the Lagrange multiplier and the constraint must satisfy the KKT conditions. The associated complementarity problem is very difficult even for linearly constrained problems in continuous space [38]. Although the augmented Lagrange multiplier method can be used for solving both equality and inequality constrained problems, the constraints must be linear since any simple nonlinear constraint could lead to a nonconvex minimization problem [44]. Unfortunately, all these mathematical difficulties were not correctly addressed in the topology optimization literature (see [39,40,41]).
Canonical duality theory (CDT) is a precise methodological theory, which can be used not only for modeling complex systems within a unified framework [23], but also for solving a large class of challenging problems in nonconvex analysis and global optimization [18,25,27]. Application of this theory to general topology optimization was given recently [24,25]. It was discovered that the 0-1 integer programming in topology optimization for linear elasticity is actually equivalent to the well-known Knapsack problem, which can be solved analytically by the CDT [20]. The main goal of this paper is to present a detailed study on the canonical duality approach for solving general topology optimization problems with applications to 2-D and 3-D linear elastic structures. In the next section, the general topology optimization problem and its challenges are discussed. A mathematically correct topology optimization problem is formulated as a coupled bilevel knapsack problem. The conceptual mistakes in topology optimization and mathematical difficulties in SIMP method are discussed. Section 3 shows that the knapsack problem can be solved analytically to obtain global optimal solution at each iteration. A canonical dual algorithm for computing the globally optimal dual solution is explained in Section 4. Some fundamental issues on challenges in topology optimization and NP-hardness in computational complexity are addressed in Section 5. Numerical examples are shown in Section 6. Conclusion remarks are given in Section 7. A brief review of the canonical duality theory is provided in Appendix.

On Mathematical Models and Challenges
Let us consider an elastically deformable body that in an undeformed configuration occupies an open domain Ω ⊂ R d (d = 2, 3) with boundary Γ = ∂Ω. We assume that the body is subjected to a body force b (per unit mass) in the reference domain Ω and a given surface traction t(x) of dead-load type on the boundary Γ t ⊂ ∂Ω, while the body is fixed on the remaining boundary Γ u = ∂Ω \ Γ t . The total potential energy of this deformed body has the following standard form: where the displacement u : Ω → R d is a continuous field variable, the mass density ρ : Ω → {0, 1} is a discrete design variable; the stored energy density W (F) is an objective function (see Appendix) of the deformation gradient F = ∇u. It should be emphasized that the design variable ρ is in a discrete space subjected to a so-called knapsack condition Ω ρ(x)dΩ ≤ V c (V c > 0 is a given desired volume), classical variational method can't be used to obtain the criticality condition for ρ. Therefore, analytical methods for studying general topology optimization problems are fundamentally difficult.
By using finite element method, the domain Ω is divided into n elements {Ω e } and in each Ω e , the unknown fields can be numerically discretized as where N e is an interpolation matrix, u e is a nodal displacement vector, the binary design variable ρ e ∈ {0, 1} is used for determining whether the element Ω e is a void (ρ e = 0) or a solid (ρ e = 1). Let U a ⊂ R m be a kinetically admissible nodal displacement space, v e > 0 be the volume of the e-th element Ω e , and Thus, by substituting (2) into (1), the total potential energy functional can be numerically reformulated as a real-valued function Π h : U a × Z a → R: where C(ρ, u) = ρ T c(u), By the facts that ρ ∈ Z a is the main design variable, the displacement u depends on each given domain Ω, and ρ p = ρ ∀p ∈ R, ρ ∈ Z a , the topology optimization for general nonlinearly deformed structures should be formulated as a so-called bi-level mixed integer nonlinear programming [25]: s.t.ū(ρ) = arg min u∈Ua Π h (u, ρ).
In this formulation, Φ p (ρ, u) represents the upper-level cost function and the total potential energy Π h (u, ρ) represents the lower-level cost function. By the fact that for each givenū, the upper-level optimization is actually a typical knapsack problem, (P bk ) is essentially a coupled bilevel knapsack problem (BKP). Bilevel programming is known to be strongly NP-hard [35], and it has been proven that merely evaluating a solution for optimality is also a NPhard task [65]. Even in the simplest case of linear bilevel programs, where the lower level problem has a unique optimal solution for all the parameters, it is not likely to find a polynomial algorithm that is capable of solving the linear bilevel program to global optimality [58]. The proof for the non-existence of a polynomial time algorithm for linear bilevel problems can be found in [11].
For large deformation problems, the total potential energy Π h is usually a nonconvex function of u. Therefore, this bi-level optimization should be the most challenging problem so far in global optimization and computational mechanics.
In reality, the topology optimization is a design process, an alternative iteration method can be naturally used to solve the bilevel optimization problem, i.e.
(i) For a given ρ k−1 ∈ Z a , to solve the lower-level problem (8) first for (ii) Then, for the fixed u k ∈ U a , to solve the upper-level integer minimization problem (7) for ρ k ∈ Z a such that The canonical duality and finite element method for solving (9) has been studied extensively for computational mechanics and global optimization [14,27]. This paper will focus only on the upper-level problem (10).
Let w = c(u) ∈ R n and p = 1. The upper-level optimization (10) is equivalent to a standared linear knapsack problem ((P kp ) for short) [24]: This well-known problem in decision science makes a perfect sense in topology optimization, i.e. among all elements {Ω e }, one should keep only those who stored more deformation energy {w e } = {c e (u)}. Due to the integer constraint, even this linear 0-1 programming is listed as one of Karp's 21 NP-complete problems [43]. However, this challenging problem can be solved analytically by using the canonical duality theory.
For linear elastic structures without body force, the total potential energy is simply a quadratic function: where K(ρ) = {ρ e K e } ∈ R n×n is the overall stiffness matrix, obtained by assembling the sub-matrix ρ e K e for each element Ω e . For a given ρ ∈ Z a , the global optimal solution for the lower-level minimization problem (8) is Then the topology optimization problem for linear elastic structures can be written in the following form s.t.ū(ρ) = arg min This is a typical bilevel knapsack-quadratic optimization. By the fact that the lower-level problem has a unique solution governed by K(ρ)u = f , the single-level reduction for solving this problem leads to Remark 1 (On Minimum Compliance Problem and SIMP Method) Instead of (P le ) or (P sl ), the topology optimization problem in literature is usually formulated as the so-called minimum compliance problem [6,56]: where the linear cost 1 2 u T f is called the "mean compliance" in topology optimization [56]. If the state variable is replaced by u = K(ρ) −1 f , then (P c ) can be written as (P c sl ) : min ρ∈Za Clearly, this problem is equivalent to (P sl ) under the regularity condition, i.e. u = K(ρ) −1 f is well-defined on Z a . Instead, the given force is replaced by f = K(ρ)u such that (P c ) is commonly written in the minimum strain energy form (P s ) : min Clearly, this problem contradicts (P sl ) in the sense that the alternative iteration for solving (P s ) leads to an anti-Knapsack problem By the fact that w = c(u) = { 1 2 u T e K e u e } ≥ 0 ∈ R n ∀u ∈ U a , this antiknapsack problem has only trivial solution.
The compliance C is a well-defined concept in engineering mechanics, which is complementary to the stiffness K in the sense of C = K −1 . In continuum physics, the linear scalar-valued function u T f ∈ R is called the external (or input) energy, which is not an objective function (see Appendix). Since f is a given force, it can't be replaced by K(ρ)u. Although the cost function P c (ρ) = 1 2 f T C(ρ)f can be called as the mean compliance, it is not an objective function either. Thus, the problem (P c ) works only for linear elastic structures. Its complementary form can be called a maximum stiffness problem, which is equivalent to (P le ) in the sense that both problems produce the same results by the alternative iteration method. Therefore, it is a conceptual mistake to call the strain energy 1 2 u T K(ρ)u as the mean compliance and (P s ) as the compliance minimization. 1 Also, the compliance can be defined only for linear elasticity. For nonlinear elasticity or plasticity, even if the stiffness can be defined as the Hessian matrix K(F) = ∇ 2 W (F), the associated compliance C can't be well-defined since K(F) is usually not invertible due to the nonlinearity/nonconvexity of the strain energy W (F).
The problem (P s ) has been used extensively by many well-known methods in topology optimization, including the most popular SIMP [2,56,69]: Clearly, the integer constraint ρ ∈ {0, 1} n in (P s ) is artificially replaced by a box constrain ρ ∈ (0, 1] n in (P simp ) via the so-called power-law ρ p = {ρ p i }. Although p > 1 is called the penalization parameter, the SIMP is not a mathematically correct penalty method for solving either (P c ) or (P s ). In order to avoid the embarrassed anti-knapsack problem, the alternative iteration is not allowed by the SIMP method. Therefore, using u = C(ρ p )f the strain energy in (P simp ) is written as P s (ρ p ) = 1 2 f T C(ρ p )f . Since P s (ρ p ) is not coercive on its domain, unless some artificial techniques are adopted, the global minimum solution of (P simp ) can be achieved only on its boundary but can never be 0 due to the restriction ρ > 0. The so-called "magic number" p = 3 works only for certain materials. This is the reason why the SIMP method suffers from the fatal drawbacks of gray scale elements and checkerboard patterns. ♣ To understand this remark, let us consider a 2-D problem (see [60] and Example 2 in Section 5) with where a, b ≥ 0 are material constants. For a given f = {f i } ∈ R 2 , the so-called penalized mean compliance is Fig. 1 shows that P s (ρ p ) is not a coercive function for p = 1, 3, and its global min can be obtained only at the boundary ρ 1 + ρ 2 = 1. Fig. 2 shows clearly that P c (ρ) is strictly convex on its boundary ρ 1 + ρ 2 = 1. The feasible solutions for (P c ) can be achieved only on the corners ρ = (1, 0) and ρ = (0, 1), both of them are global maximizers. Therefore, the problem (P s ) is considered to be NP-hard by traditional theories. For the SIMP problem (P simp ) with p = 2, its global minimizer is ρ = (0.5, 0.5), which is not a solution to the topology optimization problem. Although for p = 3 the global min of P s (ρ p ) can be achieved at either ρ = (1, 0) or ρ = (0, 1) (see Fig. 2 (a)), these two solutions can't be obtained by any Newton-type method since they are not a critical points. Moreover, if the constant a = (2 − √ 2)/4, the global min of P s (ρ 3 ) is again ρ = (0.5, 0.5) (see Fig. 2(b)). This shows that p = 3 in the SIMP method is a "magic number" only for certain materials/structures. By the fact that P s (ρ p ) is nonconvex for any given p > 1, the SIMP problem (P simp ) is a typical "NP-hard" box constrained nonconvex minimization subjected to the knapsack condition, which can't be solved deterministically by all traditional methods in polynomial time.

Canonical Dual Solution to Knapsack Problem
The canonical duality theory for solving constrained quadratic 0-1 integer programming problems was first proposed by Gao in 2007 [20]. Applications have been given to general problems in operations research [21] and recently to topology optimization for general materials [24]. In this paper, we focus mainly on linear elastic structures.
Following the standard procedure in the canonical duality theory (see [20,30]), the canonical measure for 0-1 integer programming problem (P kp ) can be given as which is a convex and lower semi-continuous (l.s.c) function in R n+1 . By this function, the knapsack problem (P kp ) can be relaxed in the following unconstrained minimization form: By the convexity of Ψ(ξ), its conjugate function can be defined uniquely via the Fenchel transformation: where Thus, by using the Fenchel-Young equality the function Ψ(ρ) in (27) can be written as the Gao-Strang total complementary function [33], Based on this function, the canonical dual of Π u (ρ) can be defined by where in which, By the standard complementarydual principle in the canonical duality theory, we have the main result: Proof: By the convexity of Ψ(ξ), we have the following canonical duality relations [17,18]: where is the sub-differential of Ψ. Thus, in terms of ξ = Λ(ρ) and ς = {σ, τ }, the canonical duality relations (35) can be equivalently written as These are exactly the KKT conditions 2 for the inequality constraints ρ The equality (34) holds due to the canonical duality relations in (35). 2

By the complementarity condition
Then for any given ς = {σ, τ } ∈ S + a , the function Ξ(·, ς) : R n → R is convex and the canonical dual function of P u in (P kp ) can be well-defined by Thus, the canonical dual problem of (P kp ) can be proposed as the following: This is a concave maximization problem over a convex subset in continuous space, which can be solved via well-developed convex optimization methods to obtain global optimal solution. Thus, whence a canonical dual solutionς is obtained, the solutionρ to the primal problem can be defined in an analytical form ofς.
Theorem 2 (Analytical Solution) For any given u ∈ U a such that w = c(u) ∈ R n + , ifς = (σ,τ ) ∈ S + a is a solution to (P d u ), then is a unique global optimal solution to (P kp ) and Proof : It is easy to prove that for any given u ∈ U a , the canonical dual function Ifς is a KKT point of P d u (ς), then it must be a unique global maximizer of P d u (ς) on S + a . By Theorem 1 we know that ifς = {σ,τ } ∈ S + a is a KKT point of Π d u (ς), thenρ = ρ(ς) defined by (41) must be a KKT point of Π u (ρ) (see [20,24]). Since Ξ(ρ, ς) is a saddle function on R n × S + a , we have Sinceσ > 0, the complementarity condition in (36) leads tō Thus, we have Theorem 2 shows that although the canonical dual problem is a concave maximization in continuous space, it produces the analytical solution (41) to the well-known integer Knapsack problem (P kp )! This analytical solution was first proposed by Gao in 2007 for general quadratic integer programming problems (see Theorem 3, [20]). The indicator function of a convex set in (26) and its sub-differential were first introduced by J.J. Moreau in 1963 in his study on unilateral constrained problems in contact mechanics [50]. His pioneering work laid a foundation for modern analysis and the canonical duality theory.
In solid mechanics, the indicator of a plastic yield condition is also called a super-potential. Its sub-differential leads to a general constitutive law and a unified pan-penalty finite element method in plastic limit analysis [13]. In mathematical programming, the canonical duality leads to a unified framework for nonlinear constrained optimization problems in multi-scale systems [32,44].

Canonical Penalty-Duality Method and Algorithm
According to Theorem 2, the global optimal solution to (P kp ) can be obtained by solving its canonical dual problem (P d kp ). However, the rate of convergence could be very slow since P d u (σ, τ ) is nearly a linear function of σ ∈ R n + when σ is far from its origin. In order to overcome this problem, a so-called βperturbed canonical dual method was proposed by Gao and Ruan in integer programming [30]. This β-perturbation is actually based on a so-called canonical penalty-duality method, i.e. the integer constraint in (P kp ) is first written in the canonical form = ρ•ρ−ρ = 0, then is relaxed by the external penalty method, while the volume constraint is simply relaxed by the Lagrange multiplier τ ≥ 0, thus, the knapsack problem can be reformulated as the following canonical penalty-duality form where β > 0 is a penalty parameter. As we can see clearly that due to the nonlinearity of the canonical constraint (ρ) = 0, the canonical penalty function is nonconvex in ρ. This is the reason why the standard penalty method can't be used for solving general nonlinearly constrained problems. However, by using the canonical duality theory, this nonconvex min-max problem can be equivalently reformulated as a canonical dual problem: which is exactly the same canonical penalty-duality problem proposed by Gao for solving the knapsack problem in topology optimization [24].
Theorem 3 (Perturbed Solution to Knapsack Problem) For any given u ∈ U a , V c > 0, and β > 0 such that w = c(u) ∈ R n + , the problem (P d β )) has at most one solution ς β = (σ β , τ β ) ∈ S + a . Moreover, there exists a β c 0 such that for any given β ∈ [β c , +∞), the vector is a global optimal solution to (P kp ).
Proof. It is easy to show that for any given β > 0, is strictly concave on the open convex set S + a . Thus, (P d β ) has a unique solution if w e − τ v e = 0 ∀e = 1, . . . , n. Indeed, the criticality condition ∇P d β (ς) = 0 leads to the following canonical dual algebraic equations: It was proved by the author (see Section 3.4.3, [18] and the Appendix of this paper) that for any given β > 0 and θ e (τ ) = τ v e − w e = 0, e = 1, . . . , n, the canonical dual algebraic equation (47) has a unique positive real solution where and φ c e is the complex conjugate of φ e , i.e. φ e φ c e = 1. Also, the canonical dual algebraic equation (48) has a unique solution This shows that the perturbed canonical dual problem (P d β ) has a unique solution ς β in S + a if θ e (τ ) = τ v e − w e = 0, e = 1, . . . , n. Thus, the density distribution ρ β can be analytically obtained by substituting (49) and (50) into (45). By the fact that lim β→∞ Π d β (ς) = Π d u (ς), there must exists a β c > 0 such that By Theorem 2 we know that this perturbed solution must be a global minimum solution to the knapsack problem. A similar proof of this theorem was given in [30]. 2 By Theorems 2 and 3 we know that for a given desired volume V c > 0, the optimal density distribution can be analytically obtained in terms of its canonical dual solution in continuous space. By the fact that (P bk ) is a bilevel mixed integer nonlinear programming, numerical optimization depends sensitively on the initial volume V 0 . If µ c = V c /V 0 1, any given iteration method could lead to unreasonable numerical solutions. In order to resolve this problem, a volume reduction control parameter µ ∈ (µ c , 1) was introduced in [24] to produce a volume sequence V γ = µV γ−1 (γ = 1, . . . , γ c ) such that V γc = V c and for any given V γ ∈ [V c , V 0 ], the problem (P bl ) is replaced by The initial values for solving this γ-th problem are V γ−1 , u γ−1 , ρ γ−1 . Based on the above strategies, the canonical duality algorithm (i.e. CDT [24]) for solving the general topology optimization problem (P bk ) can be proposed in Algorithm 1.

Remark 2 (Volume Reduction and Computational Complexity)
Theoretically speaking, for any given sequence {V γ } we should have In reality, different sequence {V γ } may produce totally different structural topology. This is an intrinsic difficulty for all bi-level optimal design problems. By the facts that there are only two loops in the CDT algorithm, i.e. the γloop and the k-loop, and the canonical dual solution is analytically given in the k-loop, the main computing is the m × m matrix inversion in the γ-loop.
The complexity for the Gauss-Jordan elimination is O(m 3 ). Therefore, the CDT is a polynomial-time algorithm.
The optimization problem of BESO as formulated in [36] is posed in the form of minimization of mean compliance, i.e. the problem (P c ). Since the alternative iteration is adopted by BESO, which leads to an anti-Knapsack problem Algorithm 1 Canonical Dual Algorithm for Topology Optimization (CDT) 1: Input parameters: µ, β and error allowances: ω 1 , ω 2 > 0. 2: Initiate ρ 0 = {1} ∈ R n , V 0 , τ 0 > 0, and let γ = 1. 3: Solve the lower-level optimization problem (9): (5), V γ = max{V c , µV γ−1 }. Let k = 1. 5: Compute the canonical dual solutions σ k = {σ k e } and τ k by , then let k = k + 1, go to Step 5; Otherwise, continue. 7: Compute the upper-level solution ρ γ by by (19) , therefore, the BESO should theoretically produce only trivial solution at each volume evolution. However, instead of solving the anti-Knapsack problem (19), a comparison method is used to keep those elements which store more strain energy. So, the BESO is actually a direct method for solving the knapsack problem (P kp ). This is the reason why the numerical results obtained by BESO are similar to that by CDT (see Section 6). But, the direct method is not a polynomial-time algorithm. Due to the combinatorial complexity, this popular method is computationally expensive and can be used only for small sized problems. This is the very reason why the knapsack problem has been considered as NP-complete for all existing direct approaches. ♣

Symmetry and NP-Hardness for Knapsack Problem
There is a hidden condition in the proof of Theorem 3, which is actually the limitation of the canonical duality theory.
the canonical dual feasible set S + a = ∅ and the knapsack problem (P kp ) has a unique solution. Otherwise, if θ e (τ c ) = 0 for at least one e ∈ {1, . . . , n}, then S + a = ∅ and (P kp ) has at least two solutions.
Proof. From the proof of Theorem 3 we know that for any given β > 0, if thre exists a τ c > 0 such that the conditions in (55) hold, the canonical dual algebraic equation (47) has a unique σ e > 0 for every e = 1, . . . , n, which can be obtained analytically by (49) see Theorem 3.4.4 in [18]). Therefore, the canonical dual problem (P d kp ) has a unique solution in S + a . Correspondingly, the primal problem (P kp ) has a unique solution defined by either (41) or (45).
If θ e (τ c ) = 0 for at least one e = 1, . . . , n, the equation 2β −1 σ 3 e + σ 2 e = 0 has two solutions σ e = 0. In this case, the KKT points of (P d β ) are located on the boundary of the open set S + a . Therefore, S + a = ∅ and (P kp ) has at least two solutions. 2 Remark 3 (NP-Hard Conjecture, Symmetry, and Linear Perturbation) Theorem 4 shows that under the condition (55), the well-known knapsack problem is not NP-hard and can be solved analytically by the canonical duality theory. It is discovered recently [25] that the critical value τ c in (55) can be deterministically given by Otherwise, as long as θ e (τ c ) = 0 for at least one e = 1, . . . , n, the solution to (P kp ) can't be written in the analytical form (45), and the knapsack problem could be really NP-hard. Actually, it is a conjecture first proposed by the author in 2007 [20], i.e.

Conjecture of NP-Hardness:
A global optimization problem is NPhard if its canonical dual problem has no solution in S + a .
It is also an open problem left in [21,30]. The reason for NP-hard problems and possible solutions were discussed recently in [23].
Geometrically speaking, the reason for multiple solutions of (P kp ) is due to certain symmetry on the modeling, boundary condition and the external load. By the fact that nothing is perfect in this real world, a perfect symmetry is not allowed for any real-world problem. Mathematically speaking, if a problem has multiple solutions, this problem is not well-posed [22]. In order to solve such NP-hard problems, the key idea is to break the symmetry. A linear perturbation method has been proposed by the author and his co-workers with successful applications in hard cases of trust region method [9], nonconvex constrained optimization [51], and integer programming [67]. ♣ The symmetry plays a fundamental rule in mathematical modeling. But, it is also a main reason that leads to chaos in nonlinear dynamics, post-buckling in large deformation mechanics, NP-hard problems in complex systems [27,45], and the well-known paradox of Buridan's ass 3 .

Example 1 (Buridan's Ass)
A donkey facing two identical hay piles starves to death because reason provides no grounds for choosing to eat one rather than the other. Mathematically, this is a knapsack problem: Due to the symmetries: v 1 = v 2 = 1, and w 1 = w 2 = w, the solution to (56) is τ c = w. Therefore, θ e (τ c ) = 0 ∀e = 1, 2 and by Theorem 4 this problem has multiple (two) solutions, which is NP-Hard to this donkey.

Applications to Benchmark Problems and Novelty
The proposed CDT algorithm for topology optimization has been implemented in Matlab. The CDT code for 2-D topology optimization is based on the popular 88-line SIMP code (TOP88) proposed by Andreassen et al [2]. By the facts that the density distribution is solved analytically at each iteration and no density filter is needed, the CDT code has only 66-lines. The CDT code for 3-D topology optimization is based on the TOP3D code proposed by Liu and Tovar [42]. The CDT code has been performed for various numerical examples to test its performance. For the purpose of illustration, the applied load and geometry data are chosen as dimensionless. Young's modulus and Poisson's ratio of the material are taken as E = 1 and ν = 0.3, respectively. The stiffness matrix of the structure in CDT algorithm is given by Clearly, we have K(1) = n e=1 EK e and K(0) = n e=1 E min K e . The reason for choosing E min = 0 is to avoid singularity in computation. To compare with other approaches, the parameters penal = 3, rmin = 1.5, and ft=1.0 are used in the SIMP 88-line code, BESO code, and the TOP3D code. The error allowances are set to be ω 1 = 2e − 16 for CDT algorithm and ω 2 = 10 −2 for all methods (SIMP is usally failed to converge if ω is too small). The initial value for τ used in CDT is τ 0 = 1. We take the design domain V 0 = 1, the initial design variable ρ 0 = {1} n for both CDT and BESO algorithms. All computations are performed by a HP laptop computer with Processor Intel Core I7-4810, CPU @ 2.80GHz and memory 2.80 GB.

MBB Beam Problem
The first example is the well-known benchmark Messerschmitt-Bölkow-Blohm (MBB) beam problem in topology optimization (see Fig. 4). The design domain is L × 2h = 180 × 60. Performance of the CDT method is first tasted for different mesh resolutions. Results in Fig. 5 show that for any given mesh resolutions, the CDT method produces precise integer solutions without using filter. Clearly, the finer the resolution, the smaller the compliance with better result (almost no checkerboard). To compare with the SIMP and BESO methods, we use the same mesh resolution of 180×60 but with different volume fractions. The volume reduction rate is fixed to be µ = 0.975 for both CDT and BESO. Computational results are reported in Figure 6, which show clearly that for any given volume fraction µ c , the CDT method produces better results within the significantly short times.

2-D Cantilever Beam Problem
The second example is the 2-D classical long cantilever beam (see Fig. 7). For Fig. 7. The design domain for a long cantilever beam with external load this benchmark problem, we first let the volume fraction µ c = 0.5. Computational results obtained by the CDT and by SIMP and BESO are summarized in Fig. 8. Clearly, the precise solid-void solution produced by the CDT method is much better than that by other methods. The second test for this problem is the comparison of the three methods with different volume fractions. From Fig. 9 one can see clearly that the CDT produces mechanically sound structures with the shortest computing time. A large range of checkerboards is observed in the results by SIMP, a small range of checkerboards is observed in the results by BESO, while almost no such pattern for the proposed CDT method. For µ c = 0.4, the SIMP code does not converge and the result reported in Fig. 9(a) is the output at the 629th iteration. This test also revealed an important truth, i.e. as the volume fraction µ c is decreasing, the strain energies produced by all three methods are increasing instead of decreasing. This shows clearly that the minimum strain energy problem is incorrect for topology optimization.
By the fact that the BESO is a direct method for solving the correct knapsack problem, although it can produce the results similar to those by the CDT, it is not a polynomial-time algorithm. Fig. 10 verified this truth. For the given mesh resolution 180 × 60, and V c = 0.5, µ = 0.975, both methods produce similar target function (i.e. strain energy C(ρ, u) = 1 2 u T K(ρ)u) (see Fig. 10(a,b)), but the BESO's computing time is exponentially blowing up as the increase in the mesh numbers (see (see Fig. 10(c)). Results in Fig. 10(a,b) show that both BESO and CDT are maximizing the strain energy, i.e. they are solving the knapsack problem (P kp ), while the SIMP method is minimizing the strain energy, i.e. it is solving (P s ).
The reduction rate µ in Algorithm 1 plays an important role for both convergence and final result. We compare in Fig 11 the computational results at different evolutionary rate from µ = 0.925 to µ = 0.985. As we can see that big µ produces more delicate structure but requires more computational time. Results in Fig. 11 verify a truth in iteration method for bilevel optimization, i.e. the optimal solution depends strongly on the parameter µ.
By the fact that the alternative iteration is adopted in CDT method, the solution also depends sensitively on mesh resolution. Fig 12 shows the CDT solutions to the cantilever beam with different mesh resolutions from 20 × 8 to 200 × 80. Clearly, different mesh resolution leads to different topology. This is completely reasonable since different mesh size leads to totally different

3-D Cantilever Beam
A CDT 3-D code (CDT3D) is developed based on the 169 line 3d topology optimization code (TOP3D) by Liu and Tovar [42]. Here we only provide a simple application to the cantilever beam. Since the BESO is not a polynomialtime algorithm, we only compare the CDT with the SIMP. The parameters in TOP3D are chosen to be penal=3, rmin=1.5. The volume reduction rate used in CDT3D is µ = 0.95. For a given volume fraction µ c = 0.3, Fig. 13(b) shows the optimized beams by TOP3D. Figures 14 shows the solution by CDT3D. Fig 15 shows the CDT3D solutions to the cantilever beam with µ = 0.9 and µ = 0.94, which verified again a truth in iteration method for solving a coupled optimization problem, i.e. the optimal solution depends strongly on the parameter µ.
In order to have a closed look at inside of the 3-D beam, we increase the beam  thickness from nez = 4 to nez = 10 and decrease the volume fraction from µ = 0.5 to µ c = 0.1. Fig. 16 and Fig. 16 show results obtained by the TOP3D and CDT (with µ = 0.93) methods. We can see clearly that TOP3D produces many gray elements. While the optimal structure by the CDT is elegant and

Concluding remarks and open problems
Topology optimization has been studied via the canonical duality theory. Theoretical results presented in this paper show that this methodological theory is indeed powerful not only for solving the most challenging topology optimization problems, but also for correctly understanding and modeling multi-scale problems in complex systems. The numerical results verified that the CPD method can produce mechanically sound optimal topology, also it is much more faster than the popular SIMP and BESO methods. Specific conclusions are given below.
(1) The mathematical problem for general topology optimization should be formulated as a bi-level knapsack programming (P bk ). This model works for both linearly and nonlinearly deformed elasto-plastic structures. (2) The CDT is a polynomial-time algorithm, which can solve general topology optimization problem to obtain global optimal solution at each volume iteration. (3) The β-perturbation for solving integer programming problems is actually the canonical penalty-duality method proposed in [13]. Theorem 3 is an application of the pure complementary energy principle in nonlinear elasticity [16]. This principle plays an important role not only in nonconvex analysis and computational mechanics, but also in topology optimization, especially for large deformed structures. (4) Unless a magic method is proposed, the volume reduction method is necessary for solving general knapsack-type problems if µ c = V c /V 0 1. But the global optimal solution depends sensitively on the reduction rate µ ∈ [µ c , 1). (5) The so-called compliance minimization problem (P c ) is actually a singlelevel reduction for solving (P bk ). Alternative iteration for the strain energy minimization (P s ) leads to an anti-knapsack problem, which has only trivial solution. (6) The problem (P simp ) is a box-constrained nonconvex minimization subjected to a knapsack condition. The SIMP is not a mathematically correct penalty method for solving either (P c ) or (P s ). Even if the magic number p = 3 works for certain materials, this method can't produce correct integer solutions for any given p > 1. (7) The BESO method is actually a combination of the volume reduction and a direct method for solving the knapsack problem (P kp ). For small-scale problems, the BESO method can produce reasonable results better than by SIMP. But it is not a polynomial-time algorithm.
By the fact that the general topology optimization problem (P bk ) is a coupled, mixed integer nonlinear programming, although for each given lowerlevel solution u ∈ U a the upper-level problem can be solved analytically by the canonical duality theory, the final result produced by the CDT may not be the global optimal solution to (P bk ) since the alternative iterations for ρ and u are adopted in Algorithm 1. This is the main reason why the numerical solution produced by the CDT depends on mesh resolution and volume reduction rate µ. The main open problems include the optimal parameter µ in order to ensure the fast convergence rate with the optimal results, the influence of the mesh resolution to the optimal topology, the existence and uniqueness of a global optimization solution to (P bk ) for a given design domain V c .
This paper is just a simple application of the canonical duality theory for linear elastic topology optimization. A 66-line Matlab code for topology optimization will be available soon to public. Also a refined algorithm with a 44-line Matlab code has been developed. The canonical duality theory is particularly useful for studying nonconvex, nonsmooth, nonconservative large deformed dynamical systems [19]. Therefore, the future works include the CDT method for solving general topology optimization problems of large deformed elastoplastic structures subjected to dynamical loads.

Appendix: A Brief Review on Canonical Duality Theory
The canonical duality theory proposed in [18] comprises mainly i) a canonical dual transformation, ii) a complementary-dual principle, and iii) a triality theory.
The canonical dual transformation is a versatile methodological method which can be used to model complex systems within a unified framework and to formulate perfect dual problems without a duality gap. The complementary-dual principle presents a unified analytic solution form for general problems in continuous and discrete systems. The triality theory reveals an intrinsic duality pattern in multi-scale systems, which can be used to identify both global and local extrema, and to develop deterministic algorithms for effectively solving a wide class of nonconvex/nonsmooth/discrete optimization/variational problems.
The canonical duality theory was developed from Gao and Strang's original work [33] for solving the following nonconvex/nonsmooth variational problem where F (u) is a generalized external energy, which is a linear functional of u on its effective domain U a ⊂ U; the internal energy density W (F) is the so-called super-potential in the sense of J.J. Moreau [50], which is a nonconvex/nonsmooth function of the deformation gradient F = ∇u. Numerical discretization of (60) leads to a global optimization problem in finite dimensional space X ⊂ R n . Canonical dual finite element method for solving the nonconvex/nonsmooth mechanics problem (60) was first proposed in 1996 [14].

Remark 4 (Objectivity and Conceptual Mistakes)
Objectivity is a central concept in our daily life, related to reality and truth. In philosophy, it means the state or quality of being true even outside a subject's individual biases, interpretations, feelings, and imaginings 4 . In science, the objectivity is often attributed to the property of scientific measurement, as the accuracy of a measurement can be tested independent from the individual scientist who first reports it 5 . In continuum mechanics, a real-valued function The objectivity plays a fundamental rule in mathematical modeling, which is also called the principle of material frame indifference [52].
However, this important concept has been misused in optimization literature. As a result, the general problem in the so-called nonlinear programming is formulated as where f (x) is called the "objective function", 6 which is allowed to be any arbitrarily given function, even a linear function. Clearly, this mathematical problem is too abstract. Although it enables one to model a very wide range of problems, it comes at a price: many global optimization problems are considered to be NP-hard. Without detailed information on these arbitrarily given functions, it is impossible to have a powerful theory for solving the artificially given constrained problem (62). ♣ In linguistics, a complete and grammatically correct sentence should be composed by at least three words: subject, object, and a predicate. As a language of science, the mathematics should follow this rule. Based on the excellent works by Oden-Reddy, Strang 7 , and Tonti [52,53,62,64], a general mathematical problem was proposed in [23]: where D is a linear operator (∇ is a special case), which maps the configuration space X into a so-called intermediate space W [52,53,64]; W : W → R∪{+∞} is an objective function such that its sub-differential leads to the constitutive duality w * ∈ ∂ − W (w) ⊂ W * , which is governed by the constitutive law and all possible physical constraints; while F : X → R ∪ {−∞} is a so-called subjective function such that its over-differential leads to the action-reaction duality x * ∈ ∂ + F (x) ⊂ X * including all possible boundary conditions and geometrical constraints [17]. This subjective function must be linear on its effective domain X a = {x ∈ X | F (x) > −∞}, i.e. F (x) = x, f ∀x ∈ X a , where f ∈ X * is an given action, the feasible set X a includes all possible geometrical constraints. The predicate in (P g ) is the operator "−" and the difference Π(x) is called the target function in general problems. The object and subject are in balance only at the optimal states. By the fact that X and W can be in different dimensional spaces, this unified problem can be used to model multi-scale systems [18,34].
The unified form (P g ) covers general constrained nonconvex/nonsmooth/discrete variational and optimization problems in multi-scale complex systems [27,34] as well as all equilibrium problems in linear elasto-plasticity and mathematical physics [53,62].

Definition 1 (Properly and Well-Posed Problems) A problem is called
properly posed if for any given non-trivial input it has at least one non-trivial solution. It is called well-posed if the solution is unique.
Clearly, this definition is more general than Hadamard's well-posed problems in dynamic systems since the continuity condition is not included. Physically speaking, any real-world problems should be well-posed since all natural phenomena exist uniquely. But practically, it is difficult to model a real-world problem precisely. Therefore, properly posed problems are allowed for the canonical duality theory. This definition is important for understanding the triality theory and NP-hard problems.
It is well-known in finite deformation theory that the stored energy W (F) is an objective function if and only if there exists an objective strain tensor C = F T F and a function U (C) such that W (F) = U (C) [10]. In nonlinear elasticity, the function U (C) is usually convex (say the St Venant-Kirchhoff material) such that the duality relation C * ∈ ∂U (C) is invertible and the complementary energy U * (C * ) can be uniquely defined via the Legendre-Fenchel transformation U * (C * ) = sup C {C : C * − U (C)}. These basic truths in continuum physics laid a foundation for the canonical duality theory.
The key idea of the canonical transformation is to introduce a nonlinear operator ξ = Λ(x) : X → E and a convex, lower semi-continuous function Ψ(ξ) : E → R ∪ {+∞} such that the nonconvex function W (Dx) can be written in the canonical form and the following canonical duality relations hold where ξ; ς represents the bilinear form on E and its canonical dual E * . Thus, using the Fenchel-Young equality Ψ(ξ) = ξ; ς − Ψ * (ς), the nonconvex/nonsmooth minimization problem (63) can be equivalently written in the following min-max form where Ξ(x, ς) is the so-called total complementary function, first introduced by Gao and Strang in 1989 [33]. Particularly, if ε = Λ(x) is a homogeneous quadratic operator, then Λ(x); ς = 1 2 x T G(ς)x = G ap (x, ς) is the so-called complementary gap function, where the Hessian matrix G(ς) = ∇ 2 x G ap (x, ς) depends on the quadratic operator Λ(x). This gap function recovers the duality gap produced by the traditional Lagrange multiplier method and modern Fenchel-Moreau duality theory. By using Ξ, the canonical dual function Π d (ς) can be obtained as [21] where sta {f (x)| x ∈ X} stands for finding a stationary value of f (x), ∀x ∈ X, and G + represents a generalized inverse of G. Based on this canonical dual function, a general solution to the nonconvex problem (63) can be obtained [17]. This theorem was first proposed in 1997 for a large deformed beam theory [15]. In nonlinear elasticity, the canonical dual function Π d (ς) is the so-called pure complementary energy, first proposed in 1999 [16]. This principle solved a 50-year old open problem in finite deformation theory and is known as the Gao principle in literature [47]. Clearly, ς ∈ S + a if and only if G ap (x, ς) ≥ 0 ∀x ∈ X a . The following triality theory shows that this gap function provides a global optimality criterion in nonconvex mechanics and global optimization.
Ifς ∈ S − a = {ς ∈ E * | G(ς) ≺ 0}, then we have either This triality theory shows that the Gao-Strang gap function can be used to identify both global minimizer and the biggest local extrema. To understand the canonical duality theory, let us consider a simple nonconvex optimization in R n : min Π(x) = 1 2 β( 1 2 where α, λ > 0 are given parameters. The criticality condition ∇Π(x) = 0 leads to a nonlinear algebraic equation system β( 1 2 x 2 − λ)x = f . Clearly, this n-dimensional nonlinear algebraic equation may have multiple solutions. Due to the nonconvexity, traditional convex optimization theories can't be used to identify global minimizer. However, by the canonical dual transformation, this problem can be solved completely to obtain all extremum solutions.
This simple example reveals a fundament truth in global optimization, i.e., for a given nonconvex/discrete minimization problem in X ⊂ R n , its canonical dual is a concave maximization over a convex set S + a ⊂ R m , which can be solved easily to obtain a unique solution if S + a has no-empty interior. Very often n m. Otherwise, the primal problem could be really NP-hard. Therefore, the canonical duality theory can be used to identify if a global optimization problem is not NP-hard. Generally speaking, if a non-symmetrically distributed external force f = 0, the canonical dual feasible set S + a has noempty interior. It was proved in [21] that for quadratic 0-1 integer programming problem, if the source term f is bigger enough, the solution is simply [21]). The canonical duality theory has been used successfully for solving a large class of challenging problems in multidisciplinary fields of engineering and sciences. A comprehensive review and applications are given in the recent book [27].
Mathematics and mechanics have been complementary partners since the Newton times. Many fundamental ideas, concepts, and mathematical methods extensively used in calculus of variations and optimization are originally developed from mechanics. The canonical duality theory is a typical example developed from Gao and Strang's work on finite deformation theory, where the subjective function F (u) must be linear so that the input force f = ∂F (u) is independent of u (i.e. a dead-load). Dually, the objective function W ( ) must be nonlinear such that there exists an objective measure ξ = Λ(u) and a convex function Ψ(ξ), the canonical transformation W (Du) = Ψ(Λ(u)) holds for most real-world systems. This is the reason why the canonical duality theory can be used naturally to solve general challenging problems in multidisciplinary fields. However, the objectivity has been misused in engineering optimization. As a consequence, the subjective function F (u) = u T f has been extensively used as the objective function in topology optimization, the minimum compliance problem (P c ) has been mistakenly written as the minimum strain energy problem (P s ), and the stored energy W (Du) = 1 2 u T K(ρ)u has been confusedly called the mean compliance. Also, the objectivity has been misused in mathematical optimization. It turns out the canonical duality theory has been mistakenly challenged by M.D. Voisei and C. Zȃlinescu in sev-eral publications (cf. [27]). By oppositely choosing linear functions for W ( ) and nonlinear functions for F (u), they produced a list of "count-examples" and concluded: "a correction of this theory is impossible without falling into trivial". The conceptual mistakes in these challenges as well as in topology optimization revealed at least two important issues: 1) there exists a big gap between optimization and mechanics; 2) incorrectly using the well-defined concepts can lead to not only ridiculous arguments, but also wrong mathematical model. As V.I. Arnold indicated [3]: "In the middle of the twentieth century it was attempted to divide physics and mathematics. The consequences turned out to be catastrophic". Interested readers are recommended to read the recent papers [22,31] for further discussion.