6 Optimization of truss structures

In this chapter, we focus on the optimization of discrete truss structures. A general procedure to solve such problems involves the following steps:

1. Definition of the truss with all nodes, elements, material properties, constraints and initial design choice

2. Solution of unknown displacements for each node in the truss given the current design

3. Evaluation of the objective function and its gradients w.r.t. the design variables

4. Approximation of the problem with MMA at the current design point

5. Solution of the approximated dual problem to find the next design point

6. Repetition of steps 2,3,4 and 5 until a maximum of iterations is reached or until there is no improvement of the objective function any more

6.1 Size Optimization

One possible objective in truss optimization is the search for the cross sectional areas \(\mathbf{a}\) resulting in the stiffest truss. The corresponding objective function could be simply the displacement at a specific node or a global displacement measure \(\mathbf{u} (\mathbf{a}) \cdot \mathbf{u} (\mathbf{a})\). However, we choose to minimize the compliance \(C: \mathcal{R}^{M} \rightarrow \mathcal{R}\) defined as \[C(\mathbf{a}) = \mathbf{f} \cdot \mathbf{u}(\mathbf{a})\] of the truss instead. The minimal compliance results also in the stiffest truss, but this formulation is preferred over the pure displacement formulation, because \(C(\mathbf{a})\) is convex (Svanberg 1985). The convexity is proven by showing that the Hessian \(\nabla^2 C(\mathbf{a})\) is positive semi-definite, a property that is inherited from the positive semi-definiteness of \(\mathbf{K}(\mathbf{a})\) (see Chapter 5.2.1 in (Christensen and Klarbring 2008) for details).

We employ box constraints \(\mathcal{A}\) on the upper and lower limit of the design variables. The upper limit \(\mathbf{a}^+\) could be a manufacturing constraint that prevents very thick elements. The lower limit \(\mathbf{a}^-\) should be larger than \(0\) to prevent deletion of truss elements. If the cross sectional area of all elements connected to a node would become \(0\), i.e. all trusses are deleted, the displacement of that node would be undetermined and the stiffness matrix would become singular. In this chapter, we prevent this by demanding \(a_j^- > 0\).

A typical inequality constraint is a volume constraint for the total truss volume used in the assembly. If we would not have a volume constraint, the optimization would predict a trivial solution with all cross sections at their maximum area. However, we are interested in solutions that put a limited material resource to its best use and thus limit the total volume to be \[\mathbf{a} \cdot \mathbf{l} \le V_0\] where \(V_0\) describes the maximum volume used in the design In order to find a solution at all, we must obviously choose \(V_0 \ge \mathbf{a}^- \cdot \mathbf{l}\), otherwise it would be impossible to find any \(\mathbf{a}\) fulfilling the constraint.

In conclusion, the size optimization problem for a stiff truss with limited volume may be noted as follows:

\[\begin{aligned} \min_{\mathbf{a}} \quad & C(\mathbf{a}) = \mathbf{f} \cdot \mathbf{u}(\mathbf{a})\\ \textrm{s.t.} \quad & \mathbf{a} \cdot \mathbf{l} - V_0 \le 0 \\ & \mathbf{a} \in \mathcal{A}\\ \end{aligned} \label{eq:size_optimization}\]

The objective function \(C(\mathbf{a})\) is convex, the inequality constraint \(\mathbf{a} \cdot \mathbf{l} - V_0\) is convex, and the set \(\mathcal{A}\) is compact and convex. Therefore, finding a KKT point is equivalent to finding the global unique solution of the problem. We can denote the Lagrangian \[\mathcal{L}(\mathbf{a}, \mu) = C(\mathbf{a}) + \mu \left( \mathbf{a} \cdot \mathbf{l} - V_0 \right)\] and need to compute the gradient \(\frac{\partial \mathcal{L}}{\partial a_j}\) to solve the optimization problem. It is \[\frac{\partial \mathcal{L} (\mathbf{a}, \mu)}{\partial a_j} = \frac{\partial C (\mathbf{a})}{\partial a_j} + \mu l_j = \mathbf{f} \cdot \frac{\partial \mathbf{u} (\mathbf{a})}{\partial a_j} + \mu l_j, \label{eq:lagrange_truss_problem}\] where the gradient \(\frac{\partial \mathbf{u} (\mathbf{a})}{\partial a_j}\) remains to be determined. It can be obtained by differentiating Equation [eq:global_stiffness] w.r.t \(\mathbf{a}\) \[\begin{aligned} \frac{\partial}{\partial a_j} \mathbf{f} &= \frac{\partial}{\partial a_j} \left( \mathbf{K} (\mathbf{a}) \cdot \mathbf{u} (\mathbf{a}) \right) \\ 0 &= \frac{\partial \mathbf{K} (\mathbf{a})}{\partial a_j} \cdot \mathbf{u} (\mathbf{a}) + \mathbf{K} (\mathbf{a}) \cdot \frac{\partial \mathbf{u} (\mathbf{a})}{\partial a_j} \end{aligned}\] and rearranging the terms to yield \[\frac{\partial \mathbf{u} (\mathbf{a})}{\partial a_j} = - \mathbf{K}^{-1}(\mathbf{a}) \cdot \frac{\partial \mathbf{K} (\mathbf{a})}{\partial a_j} \cdot \mathbf{u} (\mathbf{a}).\] Substituting this expression into the previous Equation [eq:lagrange_truss_problem] gives us the sensitivity of the compliance with respect to the design variables \[\begin{aligned} \frac{\partial C (\mathbf{a})}{\partial a_j} &= - \mathbf{f} \cdot \mathbf{K}^{-1}(\mathbf{a}) \cdot \frac{\partial \mathbf{K} (\mathbf{a})}{\partial a_j} \cdot \mathbf{u} (\mathbf{a}) \\ &= - \mathbf{u} (\mathbf{a}) \cdot \frac{\partial \mathbf{K} (\mathbf{a})}{\partial a_j} \cdot \mathbf{u} (\mathbf{a}) \\ &= - \mathbf{u}_j (\mathbf{a}) \cdot \frac{\partial \mathbf{k}_j(a_j)}{\partial a_j} \cdot \mathbf{u}_j (\mathbf{a}) \\ &= - \mathbf{u}_j (\mathbf{a}) \cdot \mathbf{k}^0_j \cdot \mathbf{u}_j (\mathbf{a}) \\ &= - 2 w_j (\mathbf{a}) \label{eq:compliance_sensitivity} \end{aligned}\] with \[\mathbf{k}_j^0 = \frac{E}{l_j} \begin{pmatrix} \cos{\phi}^2 & \cos{\phi}\sin{\phi} & -\cos{\phi}^2 & -\cos{\phi}\sin{\phi} \\ \cos{\phi}\sin{\phi} & \sin{\phi}^2 & -\cos{\phi}\sin{\phi} & -\sin{\phi}^2 \\ -\cos{\phi}^2 & \cos{\phi}\sin{\phi} & \cos{\phi}^2 &\cos{\phi}\sin{\phi} \\ -\cos{\phi}\sin{\phi} & -\sin{\phi}^2 & \cos{\phi}\sin{\phi} & \sin{\phi}^2 \\ \end{pmatrix}\] and the strain energy density in each element \[w_j = \frac{1}{2} \mathbf{u}_j (\mathbf{a}) \cdot \mathbf{k}^0_j \cdot \mathbf{u}_j (\mathbf{a}). \label{eq:element_strain_energy}\] The strain energy density is always greater or equal to zero due to the positive semi-definiteness of \(\mathbf{k}^0_j\). What we have just done is called a sensitivity analysis - we obtained an expression of how much a change of variable \(a_j\) contributes to the change of the objective function \(\frac{\partial C (\mathbf{a})}{\partial a_j}\). It is very valuable to have an analytical expression of the sensitivity, as solving it numerically via finite differences causes a high computational cost and loss in precision. Finally, the gradient of the Lagrangian becomes \[\frac{\partial \mathcal{L} (\mathbf{a}, \mu)}{\partial a_j} = - 2 w_j (\mathbf{a}) + \mu l_j \label{eq:lagrange_sensitivity}\]

If we could solve this expression for \(\mathbf{a}^*(\mu)\) at the stationary point of the Lagrangian, we could use Lagrangian duality to compute the solution easily. Unfortunately, the solution is not easy, as \(w_j\) depends on the global design \(\mathbf{a}\) and not just on the design variable in its element \(a_j\). For this reason, we will use MMA to solve the problem.

The problem in Equation [eq:size_optimization] is convex, but it is not separable and explicit, because we cannot find an explicit expression for \(\mathbf{a}^*(\mu)\) at the stationary point of the Lagrangian. Therefore we linearize the problem using MMA in \(1/(a_j-L_j^k)\). The linearization of the objective function becomes \[\begin{aligned} \tilde{C}(\mathbf{a}) &= C(\mathbf{a}^k) + \sum_j \frac{\partial C}{\partial a_j}(\mathbf{a}^k) (a^k_j-L^k_j) - \sum_j \frac{\partial C}{\partial a_j}(\mathbf{a}^k) \frac{(a^k_j-L^k_j)^2}{a_j-L^k_j}\\ &=C(\mathbf{a}^k) - \sum_j 2w_j (\mathbf{a}^k) (a^k_j-L^k_j) + \sum_j 2w_j (\mathbf{a}^k) \frac{(a^k_j-L^k_j)^2}{a_j-L^k_j}. \end{aligned}\] Note that this is a linearization using lower asymptotes only, because \[\frac{\partial C }{\partial a_j}(\mathbf{a}^k) = - 2 w_j (\mathbf{a}^k) < 0 \quad \forall \mathbf{a}^k\] and thus all terms in Equations [eq:mma_start] and [eq:mma_end] involving \(p_j^k\) vanish. The approximated optimization problem becomes \[\begin{aligned} \min_{\mathbf{a}} \quad & \tilde{C} (\mathbf{a})\\ \textrm{s.t.} \quad & \mathbf{a} \cdot \mathbf{l} - V_0 \le 0 \\ & \tilde{\mathbf{a}}^{-,k} \le \mathbf{a} \le \mathbf{a}^+\\ \end{aligned}\] with lower move limits \(\tilde{\mathbf{a}}^{-,k} = \max(\mathbf{a}^-, 0.9 \mathbf{L}^k + 0.1 \mathbf{a}^k)\) preventing division by zero. The constraint is already linear and thus we do not approximate it.

The Lagrangian of the approximated problem is \[\mathcal{L}(\mathbf{a}, \mu) = \tilde{C}(\mathbf{a}) + \mu \left( \mathbf{a} \cdot \mathbf{l} - V_0 \right)\] is separable with \[\begin{split} \mathcal{L}(\mathbf{a}, \mu) &= \underbrace{C(\mathbf{a}^k) - \sum_j 2w_j (\mathbf{a}^k) (a^k_j-L^k_j) - \mu V_0}_{\mathcal{L}^0} \\ &+ \sum_j \underbrace{\left(2 w_j (\mathbf{a}^k) \frac{(a^k_j-L^k_j)^2}{a_j-L^k_j} + \mu a_j l_j \right)}_{\mathcal{L}_j (a_j)}. \end{split}\] We want to find an explicit expression for \(\mathbf{a}^*(\mu)\). Therefore we compute the stationary point of the Lagrangian \[\frac{\partial \mathcal{L}_j}{\partial a_j} = -2 w_j (\mathbf{a}^k) \frac{(a^k_j-L^k_j)^2}{(\hat{a}_j-L^k_j)^2} + \mu l_j = 0\] and rearrange the equation to find \[\hat{a}_j(\mu) = L_j^k + \sqrt{\frac{2 w_j (\mathbf{a}^k) (L^k_j-a^k_j)^2}{\mu l_j}}.\] We simply clamp this result with box constraints to obtain the explicit relation \[\mathbf{a}^* (\mu) = \textrm{clamp}\left(\hat{\mathbf{a}}(\mu), \tilde{\mathbf{a}}^{-,k}, \mathbf{a}^+\right).\] Finally we can plug this result in the Lagrangian to get the dual function \[\underline{\mathcal{L}}(\mu) = \mathcal{L}(\mathbf{a}^* (\mu), \mu)\] and solve \[\max_{\mu} \underline{\mathcal{L}}(\mu)\] for \(\mu^*>0\). This solution procedure is simple - the dual Lagrangian depends only on one variable and the problem is concave. We can use a gradient decent method to solve it quickly and utilize this analytical expression for the gradient \[\frac{\partial \underline{\mathcal{L}}}{\partial \mu}(\mu) = \mathbf{a}^* (\mu) \cdot \mathbf{l} - V_0.\]

6.2 Topology optimization

In topology optimization of trusses, we want to identify those bars which should stay in our design and those trusses which may be removed. As seen in the last examples of the previous chapter, this is closely related to size optimization - the difference is that we seek a binary design (bar or no bar) instead of continuous cross section variables. Hence, we change the box constraint to a discrete constraint with a bar being either at maximum cross sectional area or minimum cross sectional area, but nothing in between.

The topology optimization problem of minimizing compliance for a given volume constraint has similar structure as the size optimization problem. The only change is that \(\mathbf{a}\) is now either the maximum or minimum cross section. Thus, we formally have \[\begin{aligned} \min_{\mathbf{a}} \quad & C(\mathbf{a}) = \mathbf{f} \cdot \mathbf{u}(\mathbf{a})\\ \textrm{s.t.} \quad & \mathbf{a} \cdot \mathbf{l} - V_0 \le 0 \\ & a_j \in \{a_j^-, a_j^+\}\\ \end{aligned}\]

Unfortunately, the binary problem is notoriously hard to solve, because we cannot compute gradients on the solution and testing all solutions is computationally inaccessible. Hence we try to formulate a continuous relation between stiffness and a normalized variable \(\rho_j = a_j/a_j^+\) that still results in a binary result. One such formulation is called Solid Isotropic Material with Penalization (SIMP) and we denote it here as \[E(\rho_j)= \rho_j^p E\] with a penalization parameter \(p \ge 1\). The effect of penalization is shown in Figure 1.1 for typical values of \(p\). We may observe for \(p>1\) that the stiffness per invested material is unattractive for intermediate density values. For example, choosing \(\rho_j^+=0.5\) with \(p=2\) adds half an element volume to the total volume, but contributes only a quarter of the stiffness compared to a fully filled element. An optimization that tries to minimize compliance for a given volume will therefore rather favor elements that provide the full stiffness benefit or add only a minimal contribution to the total volume.

Hence, the topology optimization problem with SIMP becomes \[\begin{aligned} \min_{\pmb{\rho}} \quad & C(\pmb{\rho}) = \mathbf{f} \cdot \mathbf{u}(\pmb{\rho})\\ \textrm{s.t.} \quad & \pmb{\rho} \cdot \mathbf{V} - V_0 \le 0 \\ & \rho_j \in (0, 1] \end{aligned}\] where \(V_j=a_j l_j\) describes the element volumes and \(\rho_j >0\) prevents singularities. After introducing SIMP, the overall problem structure is (besides the different variable names) equivalent to the size optimization problem. The only change in comparison to size optimization is the computation of the sensitivity \[\begin{aligned} \frac{\partial C (\pmb{\rho})}{\partial \rho_j} &= - \mathbf{u}_j (\pmb{\rho}) \cdot \frac{\partial \mathbf{k}_j(\rho_j)}{\partial \rho_j} \cdot \mathbf{u}_j (\pmb{\rho}) \\ &= - p \rho_j^{p-1} \mathbf{u}_j (\pmb{\rho}) \cdot \mathbf{k}^0_j \cdot \mathbf{u}_j (\pmb{\rho}) \\ &= - 2 p \rho_j^{p-1} w_j (\pmb{\rho}), \end{aligned}\] which now accounts for the penalty factor. Note that for \(p=1\), this is identical to Equation [eq:compliance_sensitivity]. The entire remaining approximation procedure and solution procedure is identical to the size optimization discussed in the previous section.

6.3 Shape optimization

So far, we have performed size optimization of trusses by adjusting the cross sectional areas and we have performed topology optimization of trusses by allowing bars to disappear. However, the node positions remained fixed in all these cases. In this section, we want to optimize the position of nodes \(\mathbf{x}\) for a fixed topology \(\mathcal{E}\) and fixed cross sectional areas \(\mathbf{a}\).

The truss shape optimization for minimum compliance with a volume constraint has a structure similar to the previous problems. However, the design variable changes to the nodal positions \(\mathbf{x}\) or a subset of those nodal positions (if not all nodes are allowed to change position). While the cross sectional areas of the constraint are fixed, the element lengths are a function of design variables now. The full shape optimization problem may be denoted as follows: \[\begin{aligned} \min_{\mathbf{x}} \quad & C(\mathbf{x}) = \mathbf{f} \cdot \mathbf{u}(\mathbf{x})\\ \textrm{s.t.} \quad & g(\mathbf{x}) = \mathbf{a} \cdot \mathbf{l}(\mathbf{x}) - V_0 \le 0 \\ & \mathbf{x} \in \mathcal{X}\\ \end{aligned} \label{eq:shape_optimization}\]

Both, the target function \(C(\mathbf{x})\) and the constraint \(g(\mathbf{x})\), are in general non-linear non-separable functions. Hence, we apply MMA to both functions to get sequential convex explicitly separable approximations. We can denote the Lagrangian of the approximated functions as \[\begin{aligned} \mathcal{L}(\mathbf{x}, \mu) &= \tilde{C}(\mathbf{x}) + \mu \tilde{g}(\mathbf{x}) \label{eq:shape_lagrangian} \end{aligned}\] which is a separable function and we can use it to solve the problem with the dual method. We need to distinguish cases for each separable term \(\mathcal{L}_i(x_i, \mu)\) as the MMA approximations are different depending on the signs of \(\frac{\partial C}{\partial x_i} (\mathbf{x}^k)\) and \(\frac{\partial g}{\partial x_i} (\mathbf{x}^k)\). Hence, we distinguish four cases for the gradient computation: \[\frac{\partial \mathcal{L}_i (\mathbf{x}, \mu)}{\partial x_i} = \begin{cases} \left(\frac{U_i^k-x_i^k}{U_i^k-x_i}\right)^2 \left(\frac{\partial C}{\partial x_i} + \mu \frac{\partial g}{\partial x_i} \right) &\textrm{if} \quad \frac{\partial C}{\partial x_i} > 0, \frac{\partial g}{\partial x_i} > 0 \\ \left(\frac{U_i^k-x_i^k}{U_i^k-x_i}\right)^2 \frac{\partial C}{\partial x_i} + \mu \left(\frac{x_i^k-L_i^k}{x_i-L_i^k}\right)^2 \frac{\partial g}{\partial x_i} &\textrm{if} \quad \frac{\partial C}{\partial x_i} > 0, \frac{\partial g}{\partial x_i} <0\\ \left(\frac{x_i^k-L_i^k}{x_i-L_i^k}\right)^2 \frac{\partial C}{\partial x_i} + \mu \left(\frac{U_i^k-x_i^k}{U_i^k-x_i}\right)^2\frac{\partial g}{\partial x_i} &\textrm{if} \quad \frac{\partial C}{\partial x_i} < 0, \frac{\partial g}{\partial x_i} > 0\\ \left(\frac{x_i^k-L_i^k}{x_i-L_i^k}\right)^2 \left(\frac{\partial C}{\partial x_i} + \mu \frac{\partial g}{\partial x_i} \right) &\textrm{if} \quad \frac{\partial C}{\partial x_i}< 0, \frac{\partial g}{\partial x_i}< 0 \end{cases}\] where all gradients are evaluated at \(\mathbf{x}^k\). To find the potential minimum \(\hat{x}_i(\mu)\), we investigate for each case individually, under which conditions it vanishes. For the first case, the second term is always positive considering the KKT condition \(\mu>0\). Therefore we find that in this case \(x\rightarrow -\infty\) would be a minimum. Using the box constraint this yields \(\hat{x}_i(\mu) = x^-_i\). Analogously, the last case yields \(\hat{x}_i(\mu) = x^+_i\). The other cases are solved simply by re-arranging the equations and finally we get \[\hat{x}_i (\mu) = \begin{cases} x^-_i &\textrm{if} \quad \frac{\partial C}{\partial x_i} > 0, \frac{\partial g}{\partial x_i} > 0 \\ \frac{U_i^k\omega + L_i^k}{1+\omega} \quad \text{with} \quad \omega = \sqrt{-\mu\frac{(x_i^k-L_i^k)^2\frac{\partial g}{\partial x_i}}{(U_i^k-x_i^k)^2\frac{\partial C}{\partial x_i}}} &\textrm{if} \quad \frac{\partial C}{\partial x_i} > 0, \frac{\partial g}{\partial x_i} <0\\ \frac{U_i^k-L_i^k\omega}{1+\omega} \quad \text{with} \quad \omega = \sqrt{-\mu\frac{(U_i^k-x_i^k)^2\frac{\partial g}{\partial x_i}}{(x_i^k-L_i^k)^2\frac{\partial C}{\partial x_i}}} &\textrm{if} \quad \frac{\partial C}{\partial x_i} < 0, \frac{\partial g}{\partial x_i} > 0\\ x^+_i &\textrm{if} \quad \frac{\partial C}{\partial x_i}< 0, \frac{\partial g}{\partial x_i} < 0. \end{cases}\] With move limits \(\tilde{\mathbf{x}}^{-,k}\) and \(\tilde{\mathbf{x}}^{+,k}\) according to Equation [eq:mma_move_limits], we can find the box-constrained minimum of the approximated Lagrangian w.r.t to \(\mathbf{x}\) as \[\mathbf{x}^*(\mu) = \textrm{clamp}\left(\hat{\mathbf{x}}(\mu), \tilde{\mathbf{x}}^{-,k}, \tilde{\mathbf{x}}^{+,k} \right).\] Finally we can plug this result in the Lagrangian to get the dual fucntion \[\underline{\mathcal{L}}(\mu) = \mathcal{L}(\mathbf{x}^* (\mu), \mu)\] and solve \[\max_{\mu} \underline{\mathcal{L}}(\mu) \label{eq:shape_dual_solution}\] for \(\mu^*>0\).