- Email: [email protected]
Topology optimization for minimum compliance using a control strategy
Engineering Structures 48 (2013) 674–682
Contents lists available at SciVerse ScienceDirect
Engineering Structures journal homepage: www.elsevier.com/locate/engstruct
Topology optimization for minimum compliance using a control strategy Andrés Tovar a,⇑, Kapil Khandelwal b a Department of Mechanical Engineering, Purdue School of Engineering & Technology, Indiana University–Purdue University Indianapolis, 723 West Michigan Street, SL 260 N, Indianapolis, IN 46202, USA b Department of Civil Engineering and Geological Sciences, University of Notre Dame, 156 Fitzpatrick Hall of Engineering, Notre Dame, IN 46556, USA
a r t i c l e
i n f o
Article history: Received 14 June 2012 Revised 8 December 2012 Accepted 10 December 2012 Available online 11 January 2013 Keywords: Structural optimization Finite element analysis Control Compliance design
a b s t r a c t This paper introduces a control-based optimization algorithm to solve topology optimization problems for structures of minimum compliance. In this approach, the iterative solution process is expressed as a multivariable control system. The elements comprising the structure are numerically incorporated with sensors, controllers, and actuators. The sensors determine a response signal as a function of the problem’s sensitivity coefficients. The controllers minimize the error between this response and a corresponding setpoint obtained from the problem’s optimality conditions. The actuators modify the design variables according to the control signal while satisfying all constraints. A proportional–integral–derivative controller is shown to be computationally efficient. Numerical issues involving local minima, mesh dependency, checkerboard patterns, and intermediate densities are tackled using continuation, filtering, and penalization methods. The performance of this algorithm is demonstrated for unpenalized and penalized topology optimization problems. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Topology optimization may be described as a binary programming problem in which the objective is to find a two-phase binary structure composed of solid and void elements. In a discretized design domain, the design variables correspond to the presence and absence of material at each discrete location. Even though (stochastic) combinatorial approaches are sometimes used in topology optimization [1], they remain restricted to small-size test problems due to the associated computational cost [2]. In order to avoid the use of such expensive combinatorial methods, the binary problem is relaxed to allow elements partially filled with material. The relaxed problem formulation requires the use of interpolation functions and penalization strategies to avoid the intermediate material condition in the final design [3]. However, the relaxed condition allows the implementation of gradient-based optimization algorithms. Some of the most traditional gradient-based algorithms in topology optimization include Optimality Criterion (OC) methods [4], convex linearization (CONLIN) [5], and the method of moving asymptotes (MMA) [6], among other (dual) sequential approximate optimization (SAO) algorithms [7]. Heuristic approaches that make use of gradient information include the evolutionary structural optimization (ESO) method [8]. Difficulties of such ⇑ Corresponding author. Tel./fax: +1 317 278 7090. E-mail addresses: [email protected] (A. Tovar), [email protected] (K. Khandelwal). 0141-0296/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engstruct.2012.12.008
method arise from the greater number of iterations as compared to traditional gradient-based approaches and convergence to an optimal solution [9]. Problems involving one single constraint can be efficiently solved using OC methods [10]. These methods have demonstrated their relevance in engineering problems involving stress, displacement, natural frequency, and stability constraints [11]. This investigation presents an alternative that is numerically more efficient than OC methods and easier to implement than SAO algorithms. This approach makes use of distributed control rules and is referred to as control-based optimization (CBO). The premise of CBO is that a distributed feedback control system efficiently reduces the error between a current state and an optimum condition. This paper incorporates this premise using distributed proportional–integral–derivative (PID) controllers. The control gains are tuned using variations of the popular Ziegler–Nichols rules [12]. Minimum compliance problems are used to demonstrate the efficiency of the CBO method. The results are compared with available algorithms including OC, MMA, and other SAO approaches [7]. Section 2 of this paper describes the CBO algorithm, including the derivation of the expressions for the structure’s response, setpoint, and PID control gains. Section 3 describes the implementation of CBO in minimum compliance problems along with the implementation to mitigate numerical issues such as local minima, mesh dependency, checkerboard patterns, and intermediate densities. Section 4 presents numerical experiments, and Section 5 contains the final remarks.
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
675
2. Control-based optimization 2.1. KKT conditions and control signals The goal of the control-based optimization (CBO) method is to drive every element of the discretized design domain towards optimality satisfying Karush–Kuhn–Tucker (KKT) conditions. To this end, the structural optimization problem is transformed into a distributed control system composed of local controllers. Each controller has the simple task of minimizing the deviation of a local response with respect to a given setpoint. To illustrate this process, consider the topology optimization problem described by
min
06xi 61;i¼1;...;n
subject to
f ðxÞ ð1Þ
Fig. 1. CBO algorithm described as a feedback control system operating in every discrete element comprising the design domain.
hðxÞ ¼ 0
where x = [x1, x2, . . . , xn]T is the vector of design variables with side constraints 0 and 1, f(x) is the objective function, and h(x) an equality constraint. The KKT conditions of (1) are explicitly stated as
ri f ðx Þ þ l ri hðx Þ þ kUi kLi ¼ 0;
ð2Þ
hðx Þ ¼ 0 and 0 6 xi 6 1 for i ¼ 1; . . . ; n; kLi xi ¼ 0 and kUi 1 xi ¼ 0; for i ¼ 1; . . . ; n;
ð4Þ
kLi P 0 and kUi P 0 for i ¼ 1; . . . ; n;
ð5Þ
ð3Þ
and referred to as optimality (2), feasibility (3), complementary (4) and non-negativity (5) conditions. The Lagrange multipliers associated with the equality and side constraints are represented by l ; kUi and kLi , respectively. For elements with intermediate density, 0 < xi < 1, the complementary conditions (4) yield kUi ¼ kLi ¼ 0. For a feasible design (3), the KKT conditions are satisfied when
ri f ðx Þ þ l ri hðx Þ ¼ 0:
ð6Þ
For void elements, xi ¼ 0, the complementary conditions (4) yield kUi ¼ 0 while the non-negativity condition (5) states kLi P 0. For a feasible design, the KKT conditions are satisfied when
ri f ðx Þ þ l ri hðx Þ ¼ kLi P 0: Finally, following a similar analysis for solid elements, sible design satisfies KKT conditions when
ri f ðx Þ þ l ri hðx Þ ¼ kui 6 0:
¼ 1, a fea-
ð8Þ
ð9Þ
and the setpoint is defined by the product of the gradient of the constraint and its corresponding Lagrange multiplier,
ski ¼ lk ri hðxk Þ;
ð10Þ
then the error signal can be defined as
eki ¼ ski þ yki :
ð11Þ
From (6) follows that the KKT conditions for an element with intermediate density are satisfied when ei ¼ 0. In the same way, conditions for void (7) and solid elements (8) are satisfied when ei 6 0 and ei P 0, respectively. In the proposed CBO method, a local controller is implemented in each element in order to vary the design variable xi according to error signal ei. The iterative scheme is described as
xikþ1 ¼ xki þ Dxki ;
2.2. Feasibility The setpoint defined by (10) is a function of the Lagrange multiplier lk that is not known a priori. A feasible design xk+1 can be defined only once this multiplier is known. The functional constraint h in (1) is set to be satisfied at every iteration k, from which h(xk) = 0 or h(xk+1) = 0. At a given iteration k only lk is unknown. Therefore, the value of the Lagrange multiplier lk can be determined as the root of h(lk) = 0. The root-finding loop is implemented in the CBO algorithm as depicted in Fig. 1. The root-finding loop uses of the bisection method [13]. Since the Lagrange multiplier lk is associated to an equality constraint, there is no restriction in sign, i.e., lk 2 (1, 1). Using the following change of variable:
ð7Þ xi
If the system’s response at the discrete location i and discrete time k is defined by the gradient of the objective function,
yki ¼ ri f ðxk Þ;
obtained from the system’s response and the setpoint ski . The Lagrange multiplier lk that defines the setpoint ski in (10) is iteratively updated in order to satisfy the functional constraint h(xk+1) = 0. The initial value of the setpoint may be any constant, e.g., l0 = 1. The expected result xkþ1 is a KKT design. The root finding algorithm and i the control strategies are described in the following sections.
ð12Þ
where Dxki is a function of the error signal eki . The CBO algorithm is described by Fig. 1, in which the design xki is analyzed using the FE method to obtain the system’s response yki . The error signal eki is
lk ¼
1 xk
xk
ð13Þ
;
the root-finding problem is reduced to determining the value of xk 2 [1, 0) [ (0, 1]. The root-finding algorithm locates the value of xk, such that h(xk) = 0. The value of lk defines the setpoint skI at every iteration according to (10). 2.3. Control actions The use of a distributed controller iteratively decreases the error (11) by adjusting the material distribution. A proportional– integral–derivative (PID) controller is the most commonly used approach in feedback engineering systems [14]. The PID controller involves three independent terms, namely: proportional, integral and derivative. The proportional term calculates the change in the design variable that is proportional to the current error value. The contribution of the integral term is proportional to both the error signal and its duration in time. The derivative term is proportional to the slope of the current error signal. In a discrete dynamic system with a PID controller, the output takes the form
xkþ1 ¼ K pro eki þ K int i
k1 X
esi þ K der eki ek1 i
ð14Þ
s¼0
where Kpro, Kint, and Kder are nonnegative scalars referred to as the proportional, integral, and derivative gains, respectively. The expression (14) is referred to as a ‘‘position’’ PID algorithm [14]. The relative change in the system is determined from a ‘‘velocity’’ algorithm. To derive it, consider the output (14) defined at k,
676
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
xki ¼ K pro eik1 þ K int
k2 X esi þ K der eik1 eik2 :
ð15Þ
s¼0
Subtracting (14) from (15) yields
Dxki ¼ K 1 eki þ K 2 ek1 þ K 3 eik2 ; i
ð16Þ
where
K 1 ¼ K pro þ K der ;
K 2 ¼ K int 2K der K pro ;
¼ K der :
and K 3 ð17Þ
The velocity algorithm (16) is computationally more efficient than the position algorithm (14), as the integration is iteratively carried out with no additional memory use. For topology optimization, the iterative scheme after incorporating the side constraints is
xikþ1 ¼ min max 0; xki þ Dxki ; 1 :
ð18Þ
If (18) is a contraction mapping, then its corresponding fixed point satisfies the problem’s KKT conditions [15]. If a general PID control is considered as in (16), then the algorithm makes use of gradient information from two previous iterations, i.e., ek1 and ek2 . The i i straightforward implementation of (18) and (16) may not result in a good controller. The implementation of an effective controller requires the consideration of windup and tuning [16].
2.5. Tuning The control gains in (16) are selected to meet specific process performance. Their selection is a compromise solution between fast system response and steady-state stability. Control tuning involves the selection of the optimum values of control gains in order to achieve the maximum rate of convergence [17]. The best known tuning methods are those developed by Ziegler and Nichols (Z–N) [12]. The recommended settings for PI and PID controllers are summarized in Table 1. In this table, Ku and Tu correspond to the ultimate gain and period, respectively. The ultimate gain Ku is determined by setting Kint and Kder to zero, and increasing Kpro until the system starts oscillating. This critical value of Kpro corresponds to Ku and the oscillation period is Tu. In the proposed CBO algorithm, the oscillation is observed when the error signal eki changes in sign at some discrete location i in consecutive iterations. The ultimate period corresponds to a unit value, Tu = 1. For determining Ku using Z–N method, Kint = Kder = 0 and from (14) the proportional position controller has the form
xkþ1 ¼ min max 0; xki K pro eki ; 1 : i
ð20Þ
In (20), the search direction corresponds to eki while xki K pro corresponds to the step size. The optimum step size is given by
K pro ¼ argmin ff : h ¼ 0 for k ¼ 0g:
ð21Þ
The line search proposed in (21), which is required only for the initial design x0, is not required to be exact. Numerical experience demonstrates that the ultimate gain Ku is a fraction of K pro of the order of 10–20%, e.g., K u ¼ 0:15K pro .
2.4. Windup 2.6. Termination Windup is a phenomenon caused by the interaction of saturations and error accumulation. In topology optimization, saturations occur whenever a design variable reaches one of its limits, 0 or 1. When this happens the design variable may remain at its limit even if the magnitude of the error is significant. In a PID controller such as the one used in (16), this error will continue to be integrated over time, thereby becoming very large; in other words, the controller ‘‘winds up’’. One approach to dealing with windup is the use of back-calculation and tracking [16]. The objective of this approach is to compare the output signal from the controller Dxki with the actual response of the structure xikþ1 . Using this principle, xki is used as a tracking signal by the PID controller yielding outstanding performance. The anti-windup algorithm is depicted in Fig. 2. The anti-windup algorithm in velocity form is expressed as
xikþ1 ¼ min max 0; xki þ xki Dxki ; 1 ;
Based on the change in design variable, a suitable termination criterion is
kxkþ1 xk k1 6 ex ;
ð22Þ
where the tolerance value ex is a small positive number, e.g., ex = 102. Once the termination criterion (22) is satisfied, the dynamic system is considered to be in a steady state and no further changes in the design variables are required. The termination criterion for the root-finding loop described in Section 2.2 is
jhðxkþ1 Þj 6 eh ;
ð23Þ
where the tolerance value eh is a small positive number, e.g., eh = ex 104. 2.7. Algorithm
ð19Þ The CBO algorithm can be summarized as:
where Dxki is defined by (16).
Step 1. Initial design: Definition of the initial design x0 and control parameters following Table 1. The initial value of the Lagrange multiplier is set to l0 = 1 or x0 = 0.5 in (13). Step 2. System’s response: Evaluation of the system’s response yki (9). Step 3. New feasible design xk+1:
Table 1 Recommended PI and PID control gains by Ziegler and Nichols. Controller
Fig. 2. Anti-windup CBO algorithm: Tracking of the design variable is used to mitigate windup.
Gain Kpro
Kint
Kder
PI
0.45Ku
1:2
0.0
PID
0.60Ku
K pro Tu K pro 2:0 T u
1 8 K pro T u
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
Step 3.1 Evaluate setpoint ski (10) and error signal eki (11). Step 3.2 Obtain candidate design xkþ1 using the PID controller i (16). Step 3.3 If the functional constraint h(xk) = 0 is not satisfied within a reasonable tolerance (23), update the Lagrange multiplier lk using xk in (13) within the interval [1, 0) [ (0, 1]. Go to Step 3.1. Step 4. Convergence check: If the termination criterion of the design variable (22) is satisfied, then x⁄ = xk+1; otherwise, k k + 1 and go to Step 2. 3. Implementation 3.1. Problem definition The implementation of the CBO method is illustrated using a minimum compliance problem described by
min
06xi 61;i¼1;...;n
subject to
f ðxÞ ¼ F T U ^ ¼ m1 m hðxÞ ¼ mðxÞ m0 0
n X ^ ¼ 0; xi q0 v i m
ð24Þ
i¼1
where the compliance f(x) is the scalar product of the structure’s force and displacement vectors, F and U, respectively. The FE equilibrium equation is KU = F, where K is the structure’s stiffness matrix. The constraint is defined in terms of the mass of the structure ^ In (24) q0 and m0 represent the m(x) limited by the mass fraction m. density and the mass of the solid structure, respectively. Hereafter, sub-index 0 will refer to solid material. The volume of each element is defined by vi. The n design variables xi correspond to the relative density of each finite element within the design domain. This is xi = qi/q0, where qi is the element density and q0 is the (solid) material density. Using the modified solid isotropic material with penalization (SIMP) interpolation function [18], the element stiffness matrix is expressed as
h
i K i ¼ Ei K 1i ¼ Emin þ xpi E0 Emin K 1i ; i i
ð25Þ
which is a constant value for all i = 1, . . . , n. Using (10), the setpoint is defined as
si ¼ l
1 q v i: m0 0
ð30Þ
3.3. Regularization schemes Topology optimization problems typically portray numerical complications such as local minima, mesh dependency, checkerboard patterns, and intermediate densities [19]. These complications require regularization schemes including continuation, filtering, and penalization [18]. In the context of the SIMP interpolation function (25), the compliance problem (24) is convex for a penalization power p = 1 and non-convex for p > 1. In order to penalize intermediate densities and mitigate the premature convergence to one of the multiple local minima when solving the non-convex problem, this investigation makes use of the continuation method. The continuation strategy incorporated into the CBO algorithm, previously presented by [7], makes pk = 1 for 20 iterations followed by pk = min{3.0, 1.02pk1}. This methodology is not proven to converge to the global optimum, but it regularizes the algorithm and allows the comparison of different optimization strategies. The sensitivity mesh-independent filtering incorporated into the CBO algorithm, previously suggested by [18], is applied on the system’s response as
P ^i ¼ y
j2N i wj xj yj
xi
P
j2Ni wj
;
Ni ¼ fj : dði; jÞ 6 Rg;
3.2. Sensitivity analysis
4.1. The MBB problem
f ðxÞ ¼ U T KU ¼
U Ti K i U i ;
ð26Þ
i¼1
where Ui and Ki are the element displacement vector and the element stiffness matrix, respectively. For finite element force control, the sensitivity of f can be obtained using direct differentiation and is given by
ri f ðxÞ ¼ U T ðri KÞU
ð27Þ
ð32Þ
where d( ) denotes distance. The weighting function wj is given by the linearly decaying (cone-shape) function
wj ¼ maxf0; R dði; jÞg;
n X
ð31Þ
where Ni is the neighborhood of the ith element and wj is a weighting function. The neighborhood is defined by the elements that have centers within a given filter radius R of the center of the ith element. This is
where K1i is the element stiffness matrix for unit Young’s modulus Ei = 1, E0 is the Young’s modulus of the solid material, and Emin is i a minimum allowed value, e.g., Emin ¼ E0 109 . The penalization i power p > 1 is used to make elements with intermediate densities inefficient as it artificially increases their compliance per unit mass.
The objective function f in (24) is described by
677
ð33Þ
which shares similarities with the one suggested by [20]. 4. Numerical experiments
The Messerschmidt–Bölkow–Blohm (MBB) beam problem (Fig. 3) is considered to demonstrate the performance of the CBO method. The linear isotropic material model in this problem assumes Young’s modulus E0 = 1.0 and Poisson’s ratio 0.3. The units are omitted but a consistent unit system is assumed. The con^ ¼ 0:5. The filter rastraint is defined by a mass fraction limit of m dius R is fixed at 8% of the base length L. Due to the symmetry of the problem, only half is analyzed. A mesh of 25 75 finite elements is used in the numerical analysis. Four-node bilinear quadrilateral elements are used in the finite element analysis. This
Using (9) and (25), the system’s response is expressed in a closed form as
yi ¼ pxp1 E0 Emin U Ti K 1i U i : i i
ð28Þ
The sensitivity of the constraint in (24) is
ri hðxÞ ¼
1 q v i; m0 0
ð29Þ
Fig. 3. MBB problem 6L L. Due to symmetry conditions, only half of 3L 3 is analyzed.
678
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
Table 2 Number of iterations of the MBB problem within a prescribed tolerance ex for p = 1 using different optimization algorithms. Results from OC, R, T2:R, E, T2:E, and MMA are presented in [7]. Tolerance ex
Algorithm
CBO: PI (Ziegler–Nichols) CBO: PID (Ziegler–Nichols) CBO: PID (Manual tuning) Optimality criterion (OC) Reciprocal approximation (R) Taylor series approximation to the reciprocal approximation (T2:R) Exponential approximation (E) Taylor series approximation to the exponential approximation (T2:E) Method of moving asymptotes (MMA)
102
104
106
108
16 14 11 12 26 25
64 33 22 44 59 58
114 50 39 77 92 90
164 67 48 111 143 105
18 17
33 35
50 51
67 62
34
68
99
153
(a)
(a)
(b) Fig. 5. Convergence of the MBB problem using the continuation method: (a) value of the penalization power p and (b) evolution of the objective function.
setting allows for the comparison of the performance of the CBO method and other topology optimization algorithms [7].
(b) Fig. 4. Convergence of the MBB problem using the CBO algorithm with p = 1: (a) evolution of the objective function and (b) change in the design as a function of the iteration number.
4.1.1. Unpenalized optimization with p = 1 The first problem consists of minimizing compliance using a penalization power p = 1. This is equivalent to using thickness as the design variable in two-dimensional problems. The initial design corresponds to a homogeneous density (or thickness) distribution. Table 2 shows the performance of several topology optimization algorithms including the proposed control-based optimization (CBO) method, optimality criterion (OC), reciprocal approximation (R), Taylor series approximation to the reciprocal approximation (T2:R), exponential approximation (E), Taylor series approximation to the exponential approximation (T2:E), and method of moving asymptotes (MMA).
Table 3 Evolution of the MBB topology as a function of the iteration number k with p = 1. Topology
Iteration k
Penalization p
Objective f
0
1.0000
253.39
11
1.0000
165.52
48
1.0000
165.49
679
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682 Table 4 Evolution of the topology as a function of the iteration number k and the penalization power p using PID (N–Z) control. Topology
Iteration k
Penalization p
Objective f
0
1.0000
253.39
20
1.0000
165.49
60
2.2080
224.99
70
2.6916
220.67
80
3.0000
204.52
100
3.0000
204.26
Table 5 Mesh independent results.
The performance is measured in terms of the number of iterations for a given prescribed tolerance ex. These results are a rather fair reflection of theses algorithms’ behavior since, with exception of MMA, they all converge almost monotonically [7]. For the CBO method, the ultimate gain corresponds to Ku = 9.5. Using Table 1, the control gains for PI control are Kpro = 4.275, Kint = 5.13, and Kder = 0.0. For PID control, the gains are Kpro = 7.5, Kint = 11.4, and Kder = 0.7125. A third CBO alternative corresponds to manual tuning with values Kpro = 9.0, Kint = 15.0, and Kder = 2.0. These manually tuned values depict the best performance for this particular problem. Clearly, for this MBB test problem the CBO method outperforms other available algorithms. The monotonic convergence of the CBO strategies is observed in Fig. 4. Table 3 shows the evolution of the design with the corresponding objective values. The differences between the solution for ex = 102 and ex = 108 are not perceptible to the naked eye. With p = 1 the solution corresponds to a gray (thickness) distribution of use in the conceptual design of the structure. Fig. 4a and b shows that both PI and PID controllers have generally monotonic convergence; however, the convergence rate is higher for the PID controller.
4.1.2. Topology optimization with continuation The continuation method along with PI and PID control strategies allows the CBO algorithm to obtain 0–1 solutions. This is demonstrated through the MBB problem together with penalization strategy described in Section 3.3. Fig. 5a shows the evolution of the penalization coefficient p, and Fig. 5b shows the corresponding convergence history of the objective function. Due to the penalization change, the value of the objective function increases before the algorithm converges. The evolution of the objective function monotonically decreases when the value of p remains constant. The initial value of the compliance is f(x0) = 253.4 for p = 1 and f(x0) = 1013.5 for p = 3. Using p = 3, there is an 80% improvement in the objective function since f(x100)/f(x0) = 0.2015. The PID controller demonstrates a higher convergence rate than the PI controller. P control is not used in this work due to poor convergence and residual steady-state error [16]. The integral contribution in the PI control is proportional to the error and its duration over time. The PI controller improves the system’s response and eliminates the residual steady-state error, but it creates
680
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
Table 6 Two-dimensional examples: wheel, cantilever, and L-shape problems.
an ‘‘overshoot’’. The resulting system’s response shown in Fig. 5b corresponds to the sum of all element responses. The derivative contribution in the PID controller slows the rate of change and reduces the magnitude of the PI overshoot. This action improves the system’s response and the rate of convergence. The evolution of the topology at discrete time intervals is also presented in Table 4. It can be observed that using the above continuation method within the CBO algorithm leads to a 0–1 design. 4.1.3. Effect of the filter size and mesh-refinement The regulation schemes (i.e., continuation and filter) used with the CBO algorithm mitigate mesh dependency phenomena. Mesh independency can be seen from the results presented in Table 5, where two filter sizes are used. Table 5 shows the topologies obtained after 100 iterations. It can be observed that topology changes with the filter size as expected. The topologies are insensitive to the mesh refinement for a given filter size. Therefore, the resulting topologies in this example are not affected by the successive mesh refinement; with appropriate filtering scheme, the CBO algorithm can produce mesh-independent results. 4.2. Two-dimensional examples The application of the CBO algorithm to 2D compliance minimization problem is shown through three examples referred to as the wheel, cantilever, and L-shape problems (Table 6). All problems assume a material model with Young’s modulus E0 = 1.0 and Poisson’s ratio of 0.3. The loads are of unit magnitude. The constraint is defined by a mass fraction limit of
^ ¼ 0:50. The solution to these problems is obtained using the m continuation method with a filter radius R. The mesh size employed for each case is also shown in Table 6. The application of the CBO method makes use of the tuning approach proposed by Ziegler–Nichols and described in Section 2.5. The ultimate gain Ku is evaluated for the initial design. The presented results are obtained after 20 and 100 iterations, respectively. Monotonic convergence is obtained for every problem. Of interest in every case are the convergence criterion at iteration 20 and the final compliance reduction at iteration 100. In every problem, the change in the density distribution jDx20j for k = 20 and p = 1 is less than 102. The compliance of the final design at k = 100 and p = 3 depicts an improvement that greater than 80% or f(x100)/f(x0) < 0.2. Similar test problems have been considered using ant colony optimization [21]. Even though the final topologies are alike, the CBO method allows including a considerably larger number of design variables at a lower computational cost.
4.3. Three-dimensional examples The application of the CBO algorithm to 3D problems is demonstrated with problems related to the ones previously shown. The examples in this section include the 3D versions of the MBB, wheel, cantilever, and L-shape problems (Table 7). The underlined triangle represents a three-dimensional planar joint. The linear isotropic material model assumes Young’s modulus E0 = 1.0 and Poisson’s ratio of 0.3. Eight-node hexahedral ‘‘brick’’ elements are used in the finite element model. The loads are of unit magnitude. The
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
681
Table 7 Three-dimensional examples: MBB, wheel, cantilever, and L-shape problems.
mesh size employed for each case is also shown in Table 7. The ^ ¼ 0:25. constraint is defined by a mass fraction limit of m The solution to these problems is obtained after 100 iterations using a filter radius R and the continuation method. Monotonic convergence is obtained for each problem. Of interest in every case is the final compliance reduction. In every problem, the ration between the initial the final compliance at k = 100 and p = 3 is
f(x100)/f(x0) < 0.04, which is an improvement greater than 96% in the structure’s performance. 5. Discussion and conclusions This paper introduces an optimization algorithm based on control theory. This algorithm, referred to as control-based optimization
682
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
(CBO), has been employed to solve topology optimization problems for structures of minimum compliance in 2D and 3D. The CBO algorithm makes use of a proportional–integral–derivate (PID) control action. The purpose of such a controller is to iteratively reduce the error eki between the system’s response yki and a setpoint ski . The problem’s KKT conditions are used to define yki and ski . The design variables xki are updated satisfying all constraints at every iteration. The arising algorithm is simple and numerically efficient. When convergence is achieved, the result is an optimum structure satisfying KKT conditions. In order to mitigate numerical complications such as mesh dependency, checkerboard patterns, intermediate densities, and multiple local minima, the CBO algorithm incorporates a sensitivity filter and the continuation method. While the application of optimization theory to control is well established [22], the use of control theory in optimization is emerging. In the last few years, PID control has been used within heuristic [23] and stochastic design optimization algorithms [24], as well as in optimization-based models of biological structures [25]. To the best of our knowledge, our work is among the first contributions of control theory to gradient-based topology optimization. The CBO method proposed in this paper shows similarities with the Optimality Criterion (OC) method for solving compliance problems [11]. Both methods make use of a closed form expression for the sensitivities to update the design variables. While the OC method uses the ratio between the sensitivity coefficients of the objective function and the functional constraint, the CBO method drives the Lagrangian sensitivity to zero using a control algorithm. A PID controller incorporated into the CBO algorithm shows improved convergence for the MBB test problem. The control gains in the CBO method are determined using Ziegler–Nichols tuning for PI and PID controllers. It is shown that if the control gains are properly tuned, the solution monotonically converges. For the test problem presented, the CBO method finds the optimum design in fewer iterations than other optimization strategies. Manual tuning may present the best performance but this procedure is problem dependent. Ziegler–Nichols tuning has demonstrated to be suitable for all test problems presented in this work. The PID controller has a better convergence rate than the PI alternative. Research on dynamic, adaptive, and self-tuning control [26–28] may be applicable to address tuning challenges in the CBO method. The CBO algorithm is simpler to implement than most sequential approximate optimization methods (SAO), e.g., reciprocal and exponential approximations. While SAO methods can be applied to a larger variety of problems (e.g., problems involving multiple constraints), they also involve iterative stages of a problem’s approximation and solution. The CBO method proposed for compliance problem does not require such approximations making its implementation simpler. In terms of the convergence rate for the test problem considered, the performance of the CBO method is shown to be comparable or superior to the performance of SAO algorithms. In the CBO algorithm, no constraints or conditions are applied among consecutive search directions, except that every direction makes the design feasible. However, one of the limitations of the proposed approach is that only one functional constraint is handled. Further work is in progress to incorporate multiple design constraints into the CBO framework.
Acknowledgements Honda R & D Americas, the Indiana Space Grant Consortium (INSGC), and the National Science Foundation (NSF CMMI1055314) supported this research. Any opinions, findings, conclusions, and recommendations expressed in this paper are those of the writers and do not necessarily reflect the views of the sponsors. References [1] Kicinger R, Arciszewski T, De Jong K. Evolutionary computation and structural design: a survey of the state-of-the-art. Comput Struct 2005;83(23-24): 1943–78. [2] Sigmund O. On the usefulness of non-gradient approaches in topology optimization. Struct Multidisc Optim 2011;43(Compendex):589–96. [3] Bendsøe MP, Sigmund O. Material interpolation schemes in topology optimization. Arch Appl Mech 1999;69(9–10):635–54. [4] Venkayya VB. Optimality criteria: a basis for multidisciplinary design optimization. Comput Mech 1989;5(1):1–21. [5] Fleury C. CONLIN: an efficient dual optimizer based on convex approximation concepts. Struct Optim 1989;1(2):81–9. [6] Svanberg K. The method of moving asymptotes – a new method for structural optimization. Int J Numer Meth Eng 1987;24(2):359–73. [7] Groenwold AA, Etman LFP. A quadratic approximation for structural topology optimization. Int J Numer Meth Eng 2010;82(4):505–24. [8] Huang X, Xie YM. Topology optimization of nonlinear structures under displacement loading. Eng Struct 2008;30(7):2057–68. [9] Rozvany GIN. A critical review of established methods of structural topology optimization. Struct Multidisc Optim 2009;37(3):217–37. [10] Patel NM et al. Comparative study of topology optimization techniques. AIAA J 2008;46(8):1963–75. [11] Bendsøe MP, Sigmund O. Topology optimization: theory. Methods and applications. Springer; 2003. [12] Ziegler JG, Nichols NB. Optimum settings for automatic controllers. Trans ASME 1942;64:759–68. [13] Mathews JH. Numerical methods: for mathematics, science, and engineering. 2nd ed. NJ: Englewood Cliffs; 1992. p. 646. [14] Shearer JL, Kulakowski BT, Gardner JF. Dynamic modeling and control of engineering systems. 2nd ed. United States of America: Prentice Hall; 1996. [15] Edwards CH. Advanced calculus of several variables. Mineola (NY): Dover; 1994. [16] Åström KJ, Hägglund T. PID controllers. International Society for Measurement and Control; 1995. [17] Penninger CL et al. Convergence analysis of hybrid cellular automata for topology optimization. Struct Multidisc Optim 2010;40(1–6):271–82. [18] Sigmund O. Morphology-based black and white filters for topology optimization. Struct Multidisc Optim 2007;33(4–5):401–24. [19] Sigmund O, Petersson J. Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Struct Optim 1998;16(1):68–75. [20] Bruns TE, Tortorelli DA. An element removal and reintroduction strategy for the topology optimization of structures and compliant mechanisms. Int J Numer Meth Eng 2003;57(10):1413–30. [21] Kaveh A et al. Structural topology optimization using ant colony methodology. Eng Struct 2008;30(9):2559–65. [22] Sage AP. Optimum systems control. Prentice Hall; 1977. p. 372. [23] Tovar A et al. Topology optimization using a hybrid cellular automaton method with local control rules. J Mech Des 2006;128(6):1205–16. [24] Cui Z et al. PID-controlled particle swarm optimization. J Multiple-Valued Logic Soft Comput 2010;16(6):585–609. [25] Andreaus U, Colloca M, Iacoviello D. An optimal control procedure for bone adaptation under mechanical stimulus. Control Eng Practice 2012;20(6): 575–83. [26] Astrom KJ, Hagglund T. Revisiting the Ziegler–Nichols step response method for PID control. J Process Control 2004;14(6):635–50. [27] Kao C-C, Chuang C-W, Fung R-F. The self-tuning PID control in a slider-crank mechanism system by applying particle swarm optimization approach. Mechatronics 2006;16(8):513–22. [28] Romero JA, Sanchis R, Balaguer P. PI and PID auto-tuning procedure based on simplified single parameter optimization. J Process Control 2011;21(6): 840–51.
Contents lists available at SciVerse ScienceDirect
Engineering Structures journal homepage: www.elsevier.com/locate/engstruct
Topology optimization for minimum compliance using a control strategy Andrés Tovar a,⇑, Kapil Khandelwal b a Department of Mechanical Engineering, Purdue School of Engineering & Technology, Indiana University–Purdue University Indianapolis, 723 West Michigan Street, SL 260 N, Indianapolis, IN 46202, USA b Department of Civil Engineering and Geological Sciences, University of Notre Dame, 156 Fitzpatrick Hall of Engineering, Notre Dame, IN 46556, USA
a r t i c l e
i n f o
Article history: Received 14 June 2012 Revised 8 December 2012 Accepted 10 December 2012 Available online 11 January 2013 Keywords: Structural optimization Finite element analysis Control Compliance design
a b s t r a c t This paper introduces a control-based optimization algorithm to solve topology optimization problems for structures of minimum compliance. In this approach, the iterative solution process is expressed as a multivariable control system. The elements comprising the structure are numerically incorporated with sensors, controllers, and actuators. The sensors determine a response signal as a function of the problem’s sensitivity coefficients. The controllers minimize the error between this response and a corresponding setpoint obtained from the problem’s optimality conditions. The actuators modify the design variables according to the control signal while satisfying all constraints. A proportional–integral–derivative controller is shown to be computationally efficient. Numerical issues involving local minima, mesh dependency, checkerboard patterns, and intermediate densities are tackled using continuation, filtering, and penalization methods. The performance of this algorithm is demonstrated for unpenalized and penalized topology optimization problems. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Topology optimization may be described as a binary programming problem in which the objective is to find a two-phase binary structure composed of solid and void elements. In a discretized design domain, the design variables correspond to the presence and absence of material at each discrete location. Even though (stochastic) combinatorial approaches are sometimes used in topology optimization [1], they remain restricted to small-size test problems due to the associated computational cost [2]. In order to avoid the use of such expensive combinatorial methods, the binary problem is relaxed to allow elements partially filled with material. The relaxed problem formulation requires the use of interpolation functions and penalization strategies to avoid the intermediate material condition in the final design [3]. However, the relaxed condition allows the implementation of gradient-based optimization algorithms. Some of the most traditional gradient-based algorithms in topology optimization include Optimality Criterion (OC) methods [4], convex linearization (CONLIN) [5], and the method of moving asymptotes (MMA) [6], among other (dual) sequential approximate optimization (SAO) algorithms [7]. Heuristic approaches that make use of gradient information include the evolutionary structural optimization (ESO) method [8]. Difficulties of such ⇑ Corresponding author. Tel./fax: +1 317 278 7090. E-mail addresses: [email protected] (A. Tovar), [email protected] (K. Khandelwal). 0141-0296/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engstruct.2012.12.008
method arise from the greater number of iterations as compared to traditional gradient-based approaches and convergence to an optimal solution [9]. Problems involving one single constraint can be efficiently solved using OC methods [10]. These methods have demonstrated their relevance in engineering problems involving stress, displacement, natural frequency, and stability constraints [11]. This investigation presents an alternative that is numerically more efficient than OC methods and easier to implement than SAO algorithms. This approach makes use of distributed control rules and is referred to as control-based optimization (CBO). The premise of CBO is that a distributed feedback control system efficiently reduces the error between a current state and an optimum condition. This paper incorporates this premise using distributed proportional–integral–derivative (PID) controllers. The control gains are tuned using variations of the popular Ziegler–Nichols rules [12]. Minimum compliance problems are used to demonstrate the efficiency of the CBO method. The results are compared with available algorithms including OC, MMA, and other SAO approaches [7]. Section 2 of this paper describes the CBO algorithm, including the derivation of the expressions for the structure’s response, setpoint, and PID control gains. Section 3 describes the implementation of CBO in minimum compliance problems along with the implementation to mitigate numerical issues such as local minima, mesh dependency, checkerboard patterns, and intermediate densities. Section 4 presents numerical experiments, and Section 5 contains the final remarks.
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
675
2. Control-based optimization 2.1. KKT conditions and control signals The goal of the control-based optimization (CBO) method is to drive every element of the discretized design domain towards optimality satisfying Karush–Kuhn–Tucker (KKT) conditions. To this end, the structural optimization problem is transformed into a distributed control system composed of local controllers. Each controller has the simple task of minimizing the deviation of a local response with respect to a given setpoint. To illustrate this process, consider the topology optimization problem described by
min
06xi 61;i¼1;...;n
subject to
f ðxÞ ð1Þ
Fig. 1. CBO algorithm described as a feedback control system operating in every discrete element comprising the design domain.
hðxÞ ¼ 0
where x = [x1, x2, . . . , xn]T is the vector of design variables with side constraints 0 and 1, f(x) is the objective function, and h(x) an equality constraint. The KKT conditions of (1) are explicitly stated as
ri f ðx Þ þ l ri hðx Þ þ kUi kLi ¼ 0;
ð2Þ
hðx Þ ¼ 0 and 0 6 xi 6 1 for i ¼ 1; . . . ; n; kLi xi ¼ 0 and kUi 1 xi ¼ 0; for i ¼ 1; . . . ; n;
ð4Þ
kLi P 0 and kUi P 0 for i ¼ 1; . . . ; n;
ð5Þ
ð3Þ
and referred to as optimality (2), feasibility (3), complementary (4) and non-negativity (5) conditions. The Lagrange multipliers associated with the equality and side constraints are represented by l ; kUi and kLi , respectively. For elements with intermediate density, 0 < xi < 1, the complementary conditions (4) yield kUi ¼ kLi ¼ 0. For a feasible design (3), the KKT conditions are satisfied when
ri f ðx Þ þ l ri hðx Þ ¼ 0:
ð6Þ
For void elements, xi ¼ 0, the complementary conditions (4) yield kUi ¼ 0 while the non-negativity condition (5) states kLi P 0. For a feasible design, the KKT conditions are satisfied when
ri f ðx Þ þ l ri hðx Þ ¼ kLi P 0: Finally, following a similar analysis for solid elements, sible design satisfies KKT conditions when
ri f ðx Þ þ l ri hðx Þ ¼ kui 6 0:
¼ 1, a fea-
ð8Þ
ð9Þ
and the setpoint is defined by the product of the gradient of the constraint and its corresponding Lagrange multiplier,
ski ¼ lk ri hðxk Þ;
ð10Þ
then the error signal can be defined as
eki ¼ ski þ yki :
ð11Þ
From (6) follows that the KKT conditions for an element with intermediate density are satisfied when ei ¼ 0. In the same way, conditions for void (7) and solid elements (8) are satisfied when ei 6 0 and ei P 0, respectively. In the proposed CBO method, a local controller is implemented in each element in order to vary the design variable xi according to error signal ei. The iterative scheme is described as
xikþ1 ¼ xki þ Dxki ;
2.2. Feasibility The setpoint defined by (10) is a function of the Lagrange multiplier lk that is not known a priori. A feasible design xk+1 can be defined only once this multiplier is known. The functional constraint h in (1) is set to be satisfied at every iteration k, from which h(xk) = 0 or h(xk+1) = 0. At a given iteration k only lk is unknown. Therefore, the value of the Lagrange multiplier lk can be determined as the root of h(lk) = 0. The root-finding loop is implemented in the CBO algorithm as depicted in Fig. 1. The root-finding loop uses of the bisection method [13]. Since the Lagrange multiplier lk is associated to an equality constraint, there is no restriction in sign, i.e., lk 2 (1, 1). Using the following change of variable:
ð7Þ xi
If the system’s response at the discrete location i and discrete time k is defined by the gradient of the objective function,
yki ¼ ri f ðxk Þ;
obtained from the system’s response and the setpoint ski . The Lagrange multiplier lk that defines the setpoint ski in (10) is iteratively updated in order to satisfy the functional constraint h(xk+1) = 0. The initial value of the setpoint may be any constant, e.g., l0 = 1. The expected result xkþ1 is a KKT design. The root finding algorithm and i the control strategies are described in the following sections.
ð12Þ
where Dxki is a function of the error signal eki . The CBO algorithm is described by Fig. 1, in which the design xki is analyzed using the FE method to obtain the system’s response yki . The error signal eki is
lk ¼
1 xk
xk
ð13Þ
;
the root-finding problem is reduced to determining the value of xk 2 [1, 0) [ (0, 1]. The root-finding algorithm locates the value of xk, such that h(xk) = 0. The value of lk defines the setpoint skI at every iteration according to (10). 2.3. Control actions The use of a distributed controller iteratively decreases the error (11) by adjusting the material distribution. A proportional– integral–derivative (PID) controller is the most commonly used approach in feedback engineering systems [14]. The PID controller involves three independent terms, namely: proportional, integral and derivative. The proportional term calculates the change in the design variable that is proportional to the current error value. The contribution of the integral term is proportional to both the error signal and its duration in time. The derivative term is proportional to the slope of the current error signal. In a discrete dynamic system with a PID controller, the output takes the form
xkþ1 ¼ K pro eki þ K int i
k1 X
esi þ K der eki ek1 i
ð14Þ
s¼0
where Kpro, Kint, and Kder are nonnegative scalars referred to as the proportional, integral, and derivative gains, respectively. The expression (14) is referred to as a ‘‘position’’ PID algorithm [14]. The relative change in the system is determined from a ‘‘velocity’’ algorithm. To derive it, consider the output (14) defined at k,
676
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
xki ¼ K pro eik1 þ K int
k2 X esi þ K der eik1 eik2 :
ð15Þ
s¼0
Subtracting (14) from (15) yields
Dxki ¼ K 1 eki þ K 2 ek1 þ K 3 eik2 ; i
ð16Þ
where
K 1 ¼ K pro þ K der ;
K 2 ¼ K int 2K der K pro ;
¼ K der :
and K 3 ð17Þ
The velocity algorithm (16) is computationally more efficient than the position algorithm (14), as the integration is iteratively carried out with no additional memory use. For topology optimization, the iterative scheme after incorporating the side constraints is
xikþ1 ¼ min max 0; xki þ Dxki ; 1 :
ð18Þ
If (18) is a contraction mapping, then its corresponding fixed point satisfies the problem’s KKT conditions [15]. If a general PID control is considered as in (16), then the algorithm makes use of gradient information from two previous iterations, i.e., ek1 and ek2 . The i i straightforward implementation of (18) and (16) may not result in a good controller. The implementation of an effective controller requires the consideration of windup and tuning [16].
2.5. Tuning The control gains in (16) are selected to meet specific process performance. Their selection is a compromise solution between fast system response and steady-state stability. Control tuning involves the selection of the optimum values of control gains in order to achieve the maximum rate of convergence [17]. The best known tuning methods are those developed by Ziegler and Nichols (Z–N) [12]. The recommended settings for PI and PID controllers are summarized in Table 1. In this table, Ku and Tu correspond to the ultimate gain and period, respectively. The ultimate gain Ku is determined by setting Kint and Kder to zero, and increasing Kpro until the system starts oscillating. This critical value of Kpro corresponds to Ku and the oscillation period is Tu. In the proposed CBO algorithm, the oscillation is observed when the error signal eki changes in sign at some discrete location i in consecutive iterations. The ultimate period corresponds to a unit value, Tu = 1. For determining Ku using Z–N method, Kint = Kder = 0 and from (14) the proportional position controller has the form
xkþ1 ¼ min max 0; xki K pro eki ; 1 : i
ð20Þ
In (20), the search direction corresponds to eki while xki K pro corresponds to the step size. The optimum step size is given by
K pro ¼ argmin ff : h ¼ 0 for k ¼ 0g:
ð21Þ
The line search proposed in (21), which is required only for the initial design x0, is not required to be exact. Numerical experience demonstrates that the ultimate gain Ku is a fraction of K pro of the order of 10–20%, e.g., K u ¼ 0:15K pro .
2.4. Windup 2.6. Termination Windup is a phenomenon caused by the interaction of saturations and error accumulation. In topology optimization, saturations occur whenever a design variable reaches one of its limits, 0 or 1. When this happens the design variable may remain at its limit even if the magnitude of the error is significant. In a PID controller such as the one used in (16), this error will continue to be integrated over time, thereby becoming very large; in other words, the controller ‘‘winds up’’. One approach to dealing with windup is the use of back-calculation and tracking [16]. The objective of this approach is to compare the output signal from the controller Dxki with the actual response of the structure xikþ1 . Using this principle, xki is used as a tracking signal by the PID controller yielding outstanding performance. The anti-windup algorithm is depicted in Fig. 2. The anti-windup algorithm in velocity form is expressed as
xikþ1 ¼ min max 0; xki þ xki Dxki ; 1 ;
Based on the change in design variable, a suitable termination criterion is
kxkþ1 xk k1 6 ex ;
ð22Þ
where the tolerance value ex is a small positive number, e.g., ex = 102. Once the termination criterion (22) is satisfied, the dynamic system is considered to be in a steady state and no further changes in the design variables are required. The termination criterion for the root-finding loop described in Section 2.2 is
jhðxkþ1 Þj 6 eh ;
ð23Þ
where the tolerance value eh is a small positive number, e.g., eh = ex 104. 2.7. Algorithm
ð19Þ The CBO algorithm can be summarized as:
where Dxki is defined by (16).
Step 1. Initial design: Definition of the initial design x0 and control parameters following Table 1. The initial value of the Lagrange multiplier is set to l0 = 1 or x0 = 0.5 in (13). Step 2. System’s response: Evaluation of the system’s response yki (9). Step 3. New feasible design xk+1:
Table 1 Recommended PI and PID control gains by Ziegler and Nichols. Controller
Fig. 2. Anti-windup CBO algorithm: Tracking of the design variable is used to mitigate windup.
Gain Kpro
Kint
Kder
PI
0.45Ku
1:2
0.0
PID
0.60Ku
K pro Tu K pro 2:0 T u
1 8 K pro T u
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
Step 3.1 Evaluate setpoint ski (10) and error signal eki (11). Step 3.2 Obtain candidate design xkþ1 using the PID controller i (16). Step 3.3 If the functional constraint h(xk) = 0 is not satisfied within a reasonable tolerance (23), update the Lagrange multiplier lk using xk in (13) within the interval [1, 0) [ (0, 1]. Go to Step 3.1. Step 4. Convergence check: If the termination criterion of the design variable (22) is satisfied, then x⁄ = xk+1; otherwise, k k + 1 and go to Step 2. 3. Implementation 3.1. Problem definition The implementation of the CBO method is illustrated using a minimum compliance problem described by
min
06xi 61;i¼1;...;n
subject to
f ðxÞ ¼ F T U ^ ¼ m1 m hðxÞ ¼ mðxÞ m0 0
n X ^ ¼ 0; xi q0 v i m
ð24Þ
i¼1
where the compliance f(x) is the scalar product of the structure’s force and displacement vectors, F and U, respectively. The FE equilibrium equation is KU = F, where K is the structure’s stiffness matrix. The constraint is defined in terms of the mass of the structure ^ In (24) q0 and m0 represent the m(x) limited by the mass fraction m. density and the mass of the solid structure, respectively. Hereafter, sub-index 0 will refer to solid material. The volume of each element is defined by vi. The n design variables xi correspond to the relative density of each finite element within the design domain. This is xi = qi/q0, where qi is the element density and q0 is the (solid) material density. Using the modified solid isotropic material with penalization (SIMP) interpolation function [18], the element stiffness matrix is expressed as
h
i K i ¼ Ei K 1i ¼ Emin þ xpi E0 Emin K 1i ; i i
ð25Þ
which is a constant value for all i = 1, . . . , n. Using (10), the setpoint is defined as
si ¼ l
1 q v i: m0 0
ð30Þ
3.3. Regularization schemes Topology optimization problems typically portray numerical complications such as local minima, mesh dependency, checkerboard patterns, and intermediate densities [19]. These complications require regularization schemes including continuation, filtering, and penalization [18]. In the context of the SIMP interpolation function (25), the compliance problem (24) is convex for a penalization power p = 1 and non-convex for p > 1. In order to penalize intermediate densities and mitigate the premature convergence to one of the multiple local minima when solving the non-convex problem, this investigation makes use of the continuation method. The continuation strategy incorporated into the CBO algorithm, previously presented by [7], makes pk = 1 for 20 iterations followed by pk = min{3.0, 1.02pk1}. This methodology is not proven to converge to the global optimum, but it regularizes the algorithm and allows the comparison of different optimization strategies. The sensitivity mesh-independent filtering incorporated into the CBO algorithm, previously suggested by [18], is applied on the system’s response as
P ^i ¼ y
j2N i wj xj yj
xi
P
j2Ni wj
;
Ni ¼ fj : dði; jÞ 6 Rg;
3.2. Sensitivity analysis
4.1. The MBB problem
f ðxÞ ¼ U T KU ¼
U Ti K i U i ;
ð26Þ
i¼1
where Ui and Ki are the element displacement vector and the element stiffness matrix, respectively. For finite element force control, the sensitivity of f can be obtained using direct differentiation and is given by
ri f ðxÞ ¼ U T ðri KÞU
ð27Þ
ð32Þ
where d( ) denotes distance. The weighting function wj is given by the linearly decaying (cone-shape) function
wj ¼ maxf0; R dði; jÞg;
n X
ð31Þ
where Ni is the neighborhood of the ith element and wj is a weighting function. The neighborhood is defined by the elements that have centers within a given filter radius R of the center of the ith element. This is
where K1i is the element stiffness matrix for unit Young’s modulus Ei = 1, E0 is the Young’s modulus of the solid material, and Emin is i a minimum allowed value, e.g., Emin ¼ E0 109 . The penalization i power p > 1 is used to make elements with intermediate densities inefficient as it artificially increases their compliance per unit mass.
The objective function f in (24) is described by
677
ð33Þ
which shares similarities with the one suggested by [20]. 4. Numerical experiments
The Messerschmidt–Bölkow–Blohm (MBB) beam problem (Fig. 3) is considered to demonstrate the performance of the CBO method. The linear isotropic material model in this problem assumes Young’s modulus E0 = 1.0 and Poisson’s ratio 0.3. The units are omitted but a consistent unit system is assumed. The con^ ¼ 0:5. The filter rastraint is defined by a mass fraction limit of m dius R is fixed at 8% of the base length L. Due to the symmetry of the problem, only half is analyzed. A mesh of 25 75 finite elements is used in the numerical analysis. Four-node bilinear quadrilateral elements are used in the finite element analysis. This
Using (9) and (25), the system’s response is expressed in a closed form as
yi ¼ pxp1 E0 Emin U Ti K 1i U i : i i
ð28Þ
The sensitivity of the constraint in (24) is
ri hðxÞ ¼
1 q v i; m0 0
ð29Þ
Fig. 3. MBB problem 6L L. Due to symmetry conditions, only half of 3L 3 is analyzed.
678
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
Table 2 Number of iterations of the MBB problem within a prescribed tolerance ex for p = 1 using different optimization algorithms. Results from OC, R, T2:R, E, T2:E, and MMA are presented in [7]. Tolerance ex
Algorithm
CBO: PI (Ziegler–Nichols) CBO: PID (Ziegler–Nichols) CBO: PID (Manual tuning) Optimality criterion (OC) Reciprocal approximation (R) Taylor series approximation to the reciprocal approximation (T2:R) Exponential approximation (E) Taylor series approximation to the exponential approximation (T2:E) Method of moving asymptotes (MMA)
102
104
106
108
16 14 11 12 26 25
64 33 22 44 59 58
114 50 39 77 92 90
164 67 48 111 143 105
18 17
33 35
50 51
67 62
34
68
99
153
(a)
(a)
(b) Fig. 5. Convergence of the MBB problem using the continuation method: (a) value of the penalization power p and (b) evolution of the objective function.
setting allows for the comparison of the performance of the CBO method and other topology optimization algorithms [7].
(b) Fig. 4. Convergence of the MBB problem using the CBO algorithm with p = 1: (a) evolution of the objective function and (b) change in the design as a function of the iteration number.
4.1.1. Unpenalized optimization with p = 1 The first problem consists of minimizing compliance using a penalization power p = 1. This is equivalent to using thickness as the design variable in two-dimensional problems. The initial design corresponds to a homogeneous density (or thickness) distribution. Table 2 shows the performance of several topology optimization algorithms including the proposed control-based optimization (CBO) method, optimality criterion (OC), reciprocal approximation (R), Taylor series approximation to the reciprocal approximation (T2:R), exponential approximation (E), Taylor series approximation to the exponential approximation (T2:E), and method of moving asymptotes (MMA).
Table 3 Evolution of the MBB topology as a function of the iteration number k with p = 1. Topology
Iteration k
Penalization p
Objective f
0
1.0000
253.39
11
1.0000
165.52
48
1.0000
165.49
679
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682 Table 4 Evolution of the topology as a function of the iteration number k and the penalization power p using PID (N–Z) control. Topology
Iteration k
Penalization p
Objective f
0
1.0000
253.39
20
1.0000
165.49
60
2.2080
224.99
70
2.6916
220.67
80
3.0000
204.52
100
3.0000
204.26
Table 5 Mesh independent results.
The performance is measured in terms of the number of iterations for a given prescribed tolerance ex. These results are a rather fair reflection of theses algorithms’ behavior since, with exception of MMA, they all converge almost monotonically [7]. For the CBO method, the ultimate gain corresponds to Ku = 9.5. Using Table 1, the control gains for PI control are Kpro = 4.275, Kint = 5.13, and Kder = 0.0. For PID control, the gains are Kpro = 7.5, Kint = 11.4, and Kder = 0.7125. A third CBO alternative corresponds to manual tuning with values Kpro = 9.0, Kint = 15.0, and Kder = 2.0. These manually tuned values depict the best performance for this particular problem. Clearly, for this MBB test problem the CBO method outperforms other available algorithms. The monotonic convergence of the CBO strategies is observed in Fig. 4. Table 3 shows the evolution of the design with the corresponding objective values. The differences between the solution for ex = 102 and ex = 108 are not perceptible to the naked eye. With p = 1 the solution corresponds to a gray (thickness) distribution of use in the conceptual design of the structure. Fig. 4a and b shows that both PI and PID controllers have generally monotonic convergence; however, the convergence rate is higher for the PID controller.
4.1.2. Topology optimization with continuation The continuation method along with PI and PID control strategies allows the CBO algorithm to obtain 0–1 solutions. This is demonstrated through the MBB problem together with penalization strategy described in Section 3.3. Fig. 5a shows the evolution of the penalization coefficient p, and Fig. 5b shows the corresponding convergence history of the objective function. Due to the penalization change, the value of the objective function increases before the algorithm converges. The evolution of the objective function monotonically decreases when the value of p remains constant. The initial value of the compliance is f(x0) = 253.4 for p = 1 and f(x0) = 1013.5 for p = 3. Using p = 3, there is an 80% improvement in the objective function since f(x100)/f(x0) = 0.2015. The PID controller demonstrates a higher convergence rate than the PI controller. P control is not used in this work due to poor convergence and residual steady-state error [16]. The integral contribution in the PI control is proportional to the error and its duration over time. The PI controller improves the system’s response and eliminates the residual steady-state error, but it creates
680
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
Table 6 Two-dimensional examples: wheel, cantilever, and L-shape problems.
an ‘‘overshoot’’. The resulting system’s response shown in Fig. 5b corresponds to the sum of all element responses. The derivative contribution in the PID controller slows the rate of change and reduces the magnitude of the PI overshoot. This action improves the system’s response and the rate of convergence. The evolution of the topology at discrete time intervals is also presented in Table 4. It can be observed that using the above continuation method within the CBO algorithm leads to a 0–1 design. 4.1.3. Effect of the filter size and mesh-refinement The regulation schemes (i.e., continuation and filter) used with the CBO algorithm mitigate mesh dependency phenomena. Mesh independency can be seen from the results presented in Table 5, where two filter sizes are used. Table 5 shows the topologies obtained after 100 iterations. It can be observed that topology changes with the filter size as expected. The topologies are insensitive to the mesh refinement for a given filter size. Therefore, the resulting topologies in this example are not affected by the successive mesh refinement; with appropriate filtering scheme, the CBO algorithm can produce mesh-independent results. 4.2. Two-dimensional examples The application of the CBO algorithm to 2D compliance minimization problem is shown through three examples referred to as the wheel, cantilever, and L-shape problems (Table 6). All problems assume a material model with Young’s modulus E0 = 1.0 and Poisson’s ratio of 0.3. The loads are of unit magnitude. The constraint is defined by a mass fraction limit of
^ ¼ 0:50. The solution to these problems is obtained using the m continuation method with a filter radius R. The mesh size employed for each case is also shown in Table 6. The application of the CBO method makes use of the tuning approach proposed by Ziegler–Nichols and described in Section 2.5. The ultimate gain Ku is evaluated for the initial design. The presented results are obtained after 20 and 100 iterations, respectively. Monotonic convergence is obtained for every problem. Of interest in every case are the convergence criterion at iteration 20 and the final compliance reduction at iteration 100. In every problem, the change in the density distribution jDx20j for k = 20 and p = 1 is less than 102. The compliance of the final design at k = 100 and p = 3 depicts an improvement that greater than 80% or f(x100)/f(x0) < 0.2. Similar test problems have been considered using ant colony optimization [21]. Even though the final topologies are alike, the CBO method allows including a considerably larger number of design variables at a lower computational cost.
4.3. Three-dimensional examples The application of the CBO algorithm to 3D problems is demonstrated with problems related to the ones previously shown. The examples in this section include the 3D versions of the MBB, wheel, cantilever, and L-shape problems (Table 7). The underlined triangle represents a three-dimensional planar joint. The linear isotropic material model assumes Young’s modulus E0 = 1.0 and Poisson’s ratio of 0.3. Eight-node hexahedral ‘‘brick’’ elements are used in the finite element model. The loads are of unit magnitude. The
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
681
Table 7 Three-dimensional examples: MBB, wheel, cantilever, and L-shape problems.
mesh size employed for each case is also shown in Table 7. The ^ ¼ 0:25. constraint is defined by a mass fraction limit of m The solution to these problems is obtained after 100 iterations using a filter radius R and the continuation method. Monotonic convergence is obtained for each problem. Of interest in every case is the final compliance reduction. In every problem, the ration between the initial the final compliance at k = 100 and p = 3 is
f(x100)/f(x0) < 0.04, which is an improvement greater than 96% in the structure’s performance. 5. Discussion and conclusions This paper introduces an optimization algorithm based on control theory. This algorithm, referred to as control-based optimization
682
A. Tovar, K. Khandelwal / Engineering Structures 48 (2013) 674–682
(CBO), has been employed to solve topology optimization problems for structures of minimum compliance in 2D and 3D. The CBO algorithm makes use of a proportional–integral–derivate (PID) control action. The purpose of such a controller is to iteratively reduce the error eki between the system’s response yki and a setpoint ski . The problem’s KKT conditions are used to define yki and ski . The design variables xki are updated satisfying all constraints at every iteration. The arising algorithm is simple and numerically efficient. When convergence is achieved, the result is an optimum structure satisfying KKT conditions. In order to mitigate numerical complications such as mesh dependency, checkerboard patterns, intermediate densities, and multiple local minima, the CBO algorithm incorporates a sensitivity filter and the continuation method. While the application of optimization theory to control is well established [22], the use of control theory in optimization is emerging. In the last few years, PID control has been used within heuristic [23] and stochastic design optimization algorithms [24], as well as in optimization-based models of biological structures [25]. To the best of our knowledge, our work is among the first contributions of control theory to gradient-based topology optimization. The CBO method proposed in this paper shows similarities with the Optimality Criterion (OC) method for solving compliance problems [11]. Both methods make use of a closed form expression for the sensitivities to update the design variables. While the OC method uses the ratio between the sensitivity coefficients of the objective function and the functional constraint, the CBO method drives the Lagrangian sensitivity to zero using a control algorithm. A PID controller incorporated into the CBO algorithm shows improved convergence for the MBB test problem. The control gains in the CBO method are determined using Ziegler–Nichols tuning for PI and PID controllers. It is shown that if the control gains are properly tuned, the solution monotonically converges. For the test problem presented, the CBO method finds the optimum design in fewer iterations than other optimization strategies. Manual tuning may present the best performance but this procedure is problem dependent. Ziegler–Nichols tuning has demonstrated to be suitable for all test problems presented in this work. The PID controller has a better convergence rate than the PI alternative. Research on dynamic, adaptive, and self-tuning control [26–28] may be applicable to address tuning challenges in the CBO method. The CBO algorithm is simpler to implement than most sequential approximate optimization methods (SAO), e.g., reciprocal and exponential approximations. While SAO methods can be applied to a larger variety of problems (e.g., problems involving multiple constraints), they also involve iterative stages of a problem’s approximation and solution. The CBO method proposed for compliance problem does not require such approximations making its implementation simpler. In terms of the convergence rate for the test problem considered, the performance of the CBO method is shown to be comparable or superior to the performance of SAO algorithms. In the CBO algorithm, no constraints or conditions are applied among consecutive search directions, except that every direction makes the design feasible. However, one of the limitations of the proposed approach is that only one functional constraint is handled. Further work is in progress to incorporate multiple design constraints into the CBO framework.
Acknowledgements Honda R & D Americas, the Indiana Space Grant Consortium (INSGC), and the National Science Foundation (NSF CMMI1055314) supported this research. Any opinions, findings, conclusions, and recommendations expressed in this paper are those of the writers and do not necessarily reflect the views of the sponsors. References [1] Kicinger R, Arciszewski T, De Jong K. Evolutionary computation and structural design: a survey of the state-of-the-art. Comput Struct 2005;83(23-24): 1943–78. [2] Sigmund O. On the usefulness of non-gradient approaches in topology optimization. Struct Multidisc Optim 2011;43(Compendex):589–96. [3] Bendsøe MP, Sigmund O. Material interpolation schemes in topology optimization. Arch Appl Mech 1999;69(9–10):635–54. [4] Venkayya VB. Optimality criteria: a basis for multidisciplinary design optimization. Comput Mech 1989;5(1):1–21. [5] Fleury C. CONLIN: an efficient dual optimizer based on convex approximation concepts. Struct Optim 1989;1(2):81–9. [6] Svanberg K. The method of moving asymptotes – a new method for structural optimization. Int J Numer Meth Eng 1987;24(2):359–73. [7] Groenwold AA, Etman LFP. A quadratic approximation for structural topology optimization. Int J Numer Meth Eng 2010;82(4):505–24. [8] Huang X, Xie YM. Topology optimization of nonlinear structures under displacement loading. Eng Struct 2008;30(7):2057–68. [9] Rozvany GIN. A critical review of established methods of structural topology optimization. Struct Multidisc Optim 2009;37(3):217–37. [10] Patel NM et al. Comparative study of topology optimization techniques. AIAA J 2008;46(8):1963–75. [11] Bendsøe MP, Sigmund O. Topology optimization: theory. Methods and applications. Springer; 2003. [12] Ziegler JG, Nichols NB. Optimum settings for automatic controllers. Trans ASME 1942;64:759–68. [13] Mathews JH. Numerical methods: for mathematics, science, and engineering. 2nd ed. NJ: Englewood Cliffs; 1992. p. 646. [14] Shearer JL, Kulakowski BT, Gardner JF. Dynamic modeling and control of engineering systems. 2nd ed. United States of America: Prentice Hall; 1996. [15] Edwards CH. Advanced calculus of several variables. Mineola (NY): Dover; 1994. [16] Åström KJ, Hägglund T. PID controllers. International Society for Measurement and Control; 1995. [17] Penninger CL et al. Convergence analysis of hybrid cellular automata for topology optimization. Struct Multidisc Optim 2010;40(1–6):271–82. [18] Sigmund O. Morphology-based black and white filters for topology optimization. Struct Multidisc Optim 2007;33(4–5):401–24. [19] Sigmund O, Petersson J. Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Struct Optim 1998;16(1):68–75. [20] Bruns TE, Tortorelli DA. An element removal and reintroduction strategy for the topology optimization of structures and compliant mechanisms. Int J Numer Meth Eng 2003;57(10):1413–30. [21] Kaveh A et al. Structural topology optimization using ant colony methodology. Eng Struct 2008;30(9):2559–65. [22] Sage AP. Optimum systems control. Prentice Hall; 1977. p. 372. [23] Tovar A et al. Topology optimization using a hybrid cellular automaton method with local control rules. J Mech Des 2006;128(6):1205–16. [24] Cui Z et al. PID-controlled particle swarm optimization. J Multiple-Valued Logic Soft Comput 2010;16(6):585–609. [25] Andreaus U, Colloca M, Iacoviello D. An optimal control procedure for bone adaptation under mechanical stimulus. Control Eng Practice 2012;20(6): 575–83. [26] Astrom KJ, Hagglund T. Revisiting the Ziegler–Nichols step response method for PID control. J Process Control 2004;14(6):635–50. [27] Kao C-C, Chuang C-W, Fung R-F. The self-tuning PID control in a slider-crank mechanism system by applying particle swarm optimization approach. Mechatronics 2006;16(8):513–22. [28] Romero JA, Sanchis R, Balaguer P. PI and PID auto-tuning procedure based on simplified single parameter optimization. J Process Control 2011;21(6): 840–51.