Computational Design of Lightweight Trusses

TL;DR

This paper presents a novel computational framework for designing lightweight trusses that optimizes geometry and topology, achieving smaller volumes and significantly faster results than previous methods.

Contribution

It introduces an alternating linear programming formulation for geometry optimization and new topological operations, including a subdivision scheme inspired by Michell's theory.

Findings

Achieves smaller truss volumes

Over two orders of magnitude faster than recent methods

Handles various functional specifications

Abstract

Trusses are load-carrying light-weight structures consisting of bars connected at joints ubiquitously applied in a variety of engineering scenarios. Designing optimal trusses that satisfy functional specifications with a minimal amount of material has interested both theoreticians and practitioners for more than a century. In this paper, we introduce two main ideas to improve upon the state of the art. First, we formulate an alternating linear programming problem for geometry optimization. Second, we introduce two sets of complementary topological operations, including a novel subdivision scheme for global topology refinement inspired by Michell's famed theoretical study. Based on these two ideas, we build an efficient computational framework for the design of lightweight trusses. \AD{We illustrate our framework with a variety of functional specifications and extensions. We show that our method achieves trusses with smaller volumes and is over two orders of magnitude faster compared with recent state-of-the-art approaches.

Figures

Tables

Table 1. Table 1. Comparison of the achieved volume and running time for different numerical optimization methods. We only compare the result of the geometry optimization part without performing topology operations or subdivision.

Fig.	GSM		ALP		SQP		SQP*		GD		GD*		[Jiang et al. 2017]
	Volume	Time(s)	Volume	Time(s)	Volume	Time(s)	Volume	Time(s)	Volume	Time(s)	Volume	Time(s)	Volume	Time(s)
16(1)	98.000	0.25	85.416	2.35	87.420	6.10	87.331	6.24	97.731	1.93	96.993	1.88	92.992	40.80
16(2)	23.868	0.38	22.318	4.01	22.318	11.47	22.318	9.17	22.662	2.53	22.638	2.33	22.790	42.50
22	6.516	0.23	6.183	2.99	6.482	1.68	6.399	2.16	6.656	1.13	6.516	1.07	6.396	29.56
23	3.125	0.23	3.041	6.52	3.569	2.23	3.099	0.77	3.101	1.39	3.100	1.15	3.124	33.13
24	3.872	0.23	3.555	5.76	4.064	3.68	3.609	1.15	4.197	1.38	3.872	1.17	3.783	31.62
25(a1)	3.375	0.48	3.238	1.50	3.536	5.89	3.238	0.72	3.371	2.03	3.356	1.22	3.343	33.26
27	19.510	0.22	19.081	0.79	21.478	2.14	19.141	0.50	19.338	1.86	19.358	1.37	19.268	48.11
28	13.891	0.36	11.885	0.54	12.000	6.57	11.901	17.40	12.767	3.70	13.217	3.48	13.653	57.30
29	15.429	0.35	14.627	1.11	17.282	20.36	14.627	17.37	15.135	3.43	15.173	2.79	16.722	70.10

Table 2. Table 2. Statistics of results presented in this paper: For each truss design, we report the number of bars, the total volume of material consumption, as well as the computational time for different stages. Please note that the results from optimal coarse trusses are already better than the results from the ground structure method without introducing additional joints and bars. Topological refinement through subdivision provides additional reduction in material consumption generally in less than one minute in our current implementation which still has a significant room for further improvement.

Fig.	Initial Truss		Optimal Coarse Truss			Truss Subdivision 1			Truss Subdivision 2			Truss Subdivision 3
	Bars	Volume GSM	Bars	Volume	Time(s)	Bars	Volume	Time(s)	Bars	Volume	Time(s)	Bars	Volume	Time(s)
22	800	6.516	14	5.830	3.76	26	5.732	0.31	50	5.709	1.41	98	5.703	2.97
23	404	3.125	22	2.967	5.86	46	2.913	5.85	118	2.909	8.86	358	2.907	15.16
24	124	3.872	16	3.365	0.64	36	3.248	0.88	100	3.203	5.46	324	3.188	17.28
25(a1)	835	3.375	19	3.210	5.02	45	3.183	1.43	133	3.170	8.79	453	3.163	29.56
25(b1)	835	3.335	25	3.114	5.42	67	3.041	1.81	211	3.021	7.77	739	3.016	17.84
26	466	4.340	25	4.294	4.92	67	4.289	3.22	211	4.288	9.68	739	4.287	35.63
27	2061	19.510	9	19.081	4.28	18	18.700	0.46	37	18.610	2.39	76	18.587	5.56
28	1118	13.891	24	11.742	5.87	60	11.724	7.73	180	11.718	19.09	612	11.713	57.41
29	3674	15.429	10	14.628	1.61	16	14.416	1.53	28	14.325	1.97	52	14.316	2.11
30	3674	30.850	24	29.049	8.70	58	28.507	1.96	174	28.285	14.75	598	28.217	64.01
31	3674	20.545	10	18.900	3.41	16	18.542	0.4	30	18.412	6.22	59	18.375	14.10

Table 3. Table 3. Comparisons with previous work. Compared with previous approaches, our framework consistently creates truss designs with smaller volumes with significantly shorter computational times.

Fig.	Method	Bars (initial)	Bars (final)	Time	Volume
1	D2013	258	96	1376s	408.807
	Ours	258	201	7.7s	333.395
17	D2013	31	19	n/a	34.977
	Ours	804	25	3.4s	34.593
18	G2003	>1 billion	n/a	>6h	4.4998
	Ours	105	2178	30s	4.4986
19	H2015	12,456,601	4244	4875s	4.3228
	Ours	105	2178	30s	4.3223
20	S2017	>7 billion	40	>1h	n/a
	Ours	1118	24	13.6s	11.742

Equations

$$\underset{a_i,s_i}{\mathrm{minimize}}$$

$${\mathbf{B}}^T\mathbf{s}=-\mathbf{f},$$

$$a_i+s_i\ge 0,$$

$$a_i-s_i\ge 0,$$

$$\underset{s_i}{\mathrm{minimize}}$$

$${\mathbf{B}}^T\mathbf{s}=-\mathbf{f}.$$

$$\underset{w_i}{\mathrm{minimize}}$$

$${\mathbf{C}}^T\mathbf{w}=-\mathbf{f}.$$

$$\underset{u_i,\Delta w_i}{\mathrm{minimize}}$$

$$(\mathbf{C}+\Delta \mathbf{C}{)}^T(\mathbf{w}+\Delta \mathbf{w})=-\mathbf{f}.$$

$$\underset{u_i,\Delta w_i}{\mathrm{minimize}}$$

$${\mathbf{C}}^T\Delta \mathbf{w}+\Delta {\mathbf{C}}^T\mathbf{w}=0,$$

$$-{\delta }_i\le \Delta w_i\le {\delta }_i;i=1,\dots ,|E|,$$

$$-{\lambda }_j\le u_{jx},u_{jy},u_{jz}\le {\lambda }_j;j=1,\dots ,|V|.$$

$$\underset{a_i,s_i^k}{\mathrm{minimize}}$$

$${\mathbf{B}}^T{\mathbf{s}}^1=-{\mathbf{f}}^1,$$

$$a_i+s_i^1\ge 0,$$

$$a_i-s_i^1\ge 0,$$

$$...$$

$${\mathbf{B}}^T{\mathbf{s}}^K=-{\mathbf{f}}^K,$$

$$a_i+s_i^K\ge 0,$$

$$a_i-s_i^K\ge 0,$$

$$\underset{u_i,\Delta {w_i}^k}{\mathrm{minimize}}$$

$${\mathbf{C}}^T\Delta {\mathbf{w}}^1+\Delta {\mathbf{C}}^T{\mathbf{w}}^1=0,$$

$$...$$

$${\mathbf{C}}^T\Delta {\mathbf{w}}^K+\Delta {\mathbf{C}}^T{\mathbf{w}}^K=0,$$

$$-{\delta }_i\le \Delta w_i^k\le {\delta }_i,$$

$$-{\lambda }_j\le u_{jx},u_{jy},u_{jz}\le {\lambda }_j,$$

$$|E|+r\ge d|V|.$$

$$\underset{{\mathbf{u}}_j,{\mathbf{w}}_i^{+},{\mathbf{w}}_i^{-}}{\mathrm{minimize}}$$

$${\mathbf{C}}^T({\mathbf{w}}^{+}-{\mathbf{w}}^{-})=-\mathbf{f},$$

$$w_i^{+}\ge 0,$$

$$w_i^{-}\ge 0,$$

$$f={\lambda }_1{E}_{volume}+{\lambda }_2{E}_{static},$$

$$E_{volume}=\sum _{i=1}^{|E|}({\mathbf{u}}_{i1}-{\mathbf{u}}_{i2}{)}^2(w_i^{+}+w_i^{-}),$$

$$E_{static}=({\mathbf{C}}^T({\mathbf{w}}^{+}-{\mathbf{w}}^{-})+\mathbf{f}{)}^2.$$

$${\mathbf{B}}^T\mathbf{s}=-\mathbf{f}.$$

$$\sum _{j=1}^{n_i}\frac{x_{p_{i,j}}-x_i}{l_{b_{i,j}}}{s}_{b_{i,j}}$$

$$\sum _{j=1}^{n_i}\frac{y_{p_{i,j}}-y_i}{l_{b_{i,j}}}{s}_{b_{i,j}}$$

$$\sum _{j=1}^{n_i}\frac{z_{p_{i,j}}-z_i}{l_{b_{i,j}}}{s}_{b_{i,j}}$$

$$\sum _{j=1}^{n_i}(x_{p_{i,j}}-x_i)w_{b_{i,j}}$$

$$\sum _{j=1}^{n_i}(y_{p_{i,j}}-y_i)w_{b_{i,j}}$$

$$\sum _{j=1}^{n_i}(z_{p_{i,j}}-z_i)w_{b_{i,j}}$$

$${\mathbf{C}}^T\mathbf{w}=-\mathbf{f}.$$

Full text

truss, topology optimization, geometry

††price: 15.00††ccs: Computing methodologies Shape modeling

1. Introduction

Trusses are crucial and fundamental structures in multiple modern engineering domains. They consist of bar elements that are connected by pin joints. Because of their efficiency and lightweight nature, trusses see considerable amount of usage in industrial design and architectural construction, e.g., for support structures of buildings, bridges, transmission towers, or even domes in playgrounds.

Designing a lightweight truss typically starts with a functional specification, e.g., in the form of external forces that the structure has to withstand. The design problem can then be formulated as an optimization problem to determine the geometry, topology, and the cross-sections of the truss. In other words, we have to find answers to the following questions: Where to put the intermediate joints? How to connect the joints with bars? What are the cross-section areas of the bars? These tasks are notoriously challenging because the optimization of geometry, topology, and cross-sections is interrelated, and there exists an infinite number of possible topologies which are difficult to classify and quantify. Even for a simple case, the optimal topology is not intuitive. As shown in Figure 2, to support two pairs of opposing forces lying on two straight lines, the simplest truss on the left with two bars, one in tension (blue) and another in compression (red), may be intuitively considered as the lightest truss. However, a lighter design with more intermediate joints and connections can be found as shown in the Figure 2 right. Another simple functional specification problem is called the three forces problem (3FP) [Chan, 1966; Sokół and Lewiński, 2010]. 3FP is formulated as follows: find the lightest fully stressed truss transmitting three self-equilibrated co-planar forces. Although there are only three forces, the problem is still unsolved analytically for general cases.

In this paper, we mainly follow previous work and take functional specifications in the form of supporting points and applied forces as input. Our goal is then to construct a lightweight truss with optimal joint positions, topology, and cross-sections. As an example shown in Figure 3 left, two supporting points and one external force are given as inputs. Our computational method generates the topology automatically and optimizes the nodal positions and cross-section areas as shown in Figure 3 right.

There are two major strategies to tackle this problem in previous work. The first strategy is to start with a densely-connected structure and to subsequently identify which bars to remove (e.g., the ground structure method [Dorn, 1964] and its variations [Gilbert and Tyas, 2003; Sokół, 2017]). The main limitation of this strategy is its sub-optimality due to its heavy dependence on the initialization as it is extremely unlikely that the optimal structure is a sub-structure of the initial one. The other strategy is to start with a sparse structure and to iteratively add new joints and bars. One of the most famous methods in [Martinez et al., 2007] adds one joint at a time and can only deal with one single load 2D problem.

Even though topology optimization is such a longstanding and fundamental problem in structural engineering, one can identify large possible improvements to the current state of the art. First, the search space for truss topology is not properly explored by previous algorithms. We could observe that it is very difficult to optimize the topology in a single stage. Much better results can be achieved by proceeding in two stages of topology optimization: computing a coarse truss and truss subdvision. Our novel subdivision approach is inspired by Michell’s pioneering theoretical treatment of optimal truss design in [Michell, 1904]. Second, the geometry optimization used in previous work is not efficient. To tackle this problem we decompose geometry optimization into alternating linear programming formulations to reduce the running time.

Our main contributions are as described in the following

• •

We propose two categories of complementary topology operations, local and global. While local operations have been used in previous work, our global operations based on subdivision are our original contribution.

• •

We introduce a novel algorithm for geometry optimization based on alternating linear programming (ALP) that jointly optimizes joint positions and bar cross sections.

• •

Based on these two technical contributions, we build a framework for lightweight truss design, a longstanding and important problem in structural engineering, architecture, graphics, and design. Compared with recent state-of-the-art approaches, our method creates trusses with smaller volumes, can handle more complex functional specifications, and is over two orders of magnitude faster.

2. Previous Work

In recent years, combining geometric modeling together with realistic engineering considerations, especially static equilibrium and manufacturability, have attracted the interests of many researchers in the graphics community. Beyond applications in the virtual world [Smith et al., 2002], those previous works enable novel and functional designs manufacturable with 3D printing  [Wang et al., 2013; Zhou et al., 2013], laser cutting  [Martínez et al., 2015], masonry structure  [Block and Ochsendorf, 2007; Vouga et al., 2012; Liu et al., 2013; Panozzo et al., 2013; de Goes et al., 2013; Tang et al., 2014], for toys  [Prévost et al., 2013; Bächer et al., 2014], furniture  [Umetani et al., 2012; Yao et al., 2017], and architecture  [Jiang et al., 2015; Pietroni et al., 2015]. The most relevant works to ours are  [Jiang et al., 2017] and  [Kilian et al., 2017]. Jiang et al. propose a framework to design and optimize space structures where only a small set of cross-section areas are allowed. Therefore, Jiang et al. [Jiang et al., 2017] compute a specialized form of truss, but we focus on the classical problem of truss design without discrete restrictions on the cross sections. The main practical difference is that our proposed method can generate a truss from scratch, whereas Jiang et al. relies on a reasonable truss being given as input. In our results we also demonstrate that our proposed optimization technique ALP produces better results on our problem formulation than the geometric optimization technique used in [Jiang et al., 2017]. Kilian et al. [Kilian et al., 2017] provide an interesting geometric understanding of ”optimality” of surface-like lightweight structures. Compared with their work, we tackle a problem for common and general trusses in both 2D and 3D, instead of focusing on load-carrying surfaces.

The problem of designing a truss with a minimal volume of material that supports imposed external forces was first studied in  [Michell, 1904]. In the milestone paper, Michell proved that an optimal truss must follow orthogonal networks of lines of maximal and minimal strains in a constant-magnitude strain field. An optimal truss is usually called a Michell truss. Following his work, research on the topic of optimal truss can be divided into two categories: exact-analytical formulations and approximate-discretized formulations.

An exact-analytical formulation assumes that the truss is a continuum structure connected by an infinite number of bars with infinitesimally small cross-sections. In analytical formulations, the theoretical optimal design is determined exactly through the simultaneous solution of a system of equations expressing the conditions for optimality. The basic principles were establish in [Maxwell, 1870] and [Michell, 1904] and a more general treatment was outlined in [Hemp, 1973] and  [Prager and Rozvany, 1977]. Recent works on deriving exact solutions were presented in [Rozvany, 1998], [Lewiński and Rozvany, 2007], [Lewiński and Rozvany, 2008b], and [Lewiński and Rozvany, 2008a] for a series of benchmark problems. Basically, the analytical solutions are very hard to obtain and only available for some special boundary conditions. While they are less practical in most of the generic scenarios, these solutions could be used as references to verify the performance of numerical methods.

Discretized numerical formulations are more practical and efficient approaches for structural design tasks presented in the real world. The most influential method is the ground structure method (GSM) which was first proposed in [Dorn, 1964]. This method consists of generating a fixed grid of joints and adding bars in some or all of the possible connections between the joints as potential structural or vanishing bars. The optimized structure for the imposed functional specification is found using the cross-section areas as design variables, and the whole problem is formulated as a linear programming problem. Its optimal topology is achieved by eliminating the zero-area cross sections. The ground structure method has been recently improved in [Gilbert and Tyas, 2003] and [Sokół, 2017].

Besides GSM, some other numerical methods are proposed recently, such as the method in [Martinez et al., 2007], carries out geometry optimization in conjunction with a heuristic ‘joint adding’ algorithm, generating an increasingly complex truss structure from a relatively simple initial layout. However, this algorithm can only add one joint per time and only works for single load cases. An efficient algorithm proposed by He and Gibbert  [2015] combines layout optimization with geometry optimization. Similar to GSM, its layout optimization starts from a densely connected truss and is formulated by a linear programming problem, and its geometry optimization is formulated by a non-linear optimization as a post-processing step.

3. Overview

Our framework has the following major components:

• •

Functional specification (C1). The input to our framework is the functional specifications including the external forces and supporting points together with a set of structural constraints, e.g., design regions, geometric obstacles, and material properties. (See Section 4.1)

• •

Initialization (C2). To obtain an initial truss, we create a grid of intermediate joints and densely connect them. This grid is located inside the design region and its size is proportional to the bounding box of the points in the input specification. (See Section 4.2)

• •

Local topology operations (C3). We locally manipulate the topology through some geometry operations such as removing bars with vanishing cross-section areas and joints without any connection, merging close joints, etc. (See Section 4.3)

• •

Global topology refinement using subdivision (C4). Use an optimized coarse truss as input, we further refine the truss through subdivision. (See Section 4.4 )

• •

Geometric optimization using ALP (C5). Given a fixed topology, we propose an alternating linear programming algorithm (ALP) to reduce the total volume of the truss by adjusting the joint positions and cross-section areas of bars. This algorithm is an essential component and its details are introduced in Section 5.

Pipeline Overview

The individual components work together as follows: We start from an input specification to compute an initial truss. Then, we proceed in two phases. In the first phase, coarse truss optimization, we interleave geometric optimization (C5) with local topology operations (C3). In the second phase, structure refinement through subdivision, we interleave global topology refinement (C4) with geometric optimization (C5). See Figure 4 for an overview of our framework. We discuss each individual component and the overall framework in detail in the next sections.

4. Design Framework

We provide a framework for the computational design of lightweight trusses. In this section, we describe the input specification, the initialization, local topology operations, and global topology refinement.

4.1. Functional Specification

The input to our framework is the functional specification including the external forces and supporting points together with a set of structural constraints, e.g., design regions, geometric obstacles, and material properties. Throughout the paper, we visualize supporting joints as red dots, joints with active forces as blue dots, and intermediate joints as yellow dots. We also visualize bars in tension in blue and bars in compression in red. In addition, the thickness of the bars is visualized in proportional to the computed cross-section areas. Note that when the external forces are in self-equilibrium, the input specification may have no supporting points. For example, the three forces problem (3FP) in 2D.

4.2. Truss Initialization

We build on previous work to compute an initial truss. There are two approaches to tackle this problem. One simply adds connections between provided joints in the functional specification (supporting joints and joints with active forces). For example, as shown in Figure 5 left, two bars connecting the joints with active forces (blue) and the supporting joints (red) are set as an initial truss. In some cases, this initialization is too simple to construct an equilibrium force system. Another method adds a grid of intermediate points over the design region and densely connects them as shown in Figure 5 right. The increasing method such as the work in  [Martinez et al., 2007] used the first initialization. The GSM usually uses the latter one with a large number of intermediate joints.

In our framework, we first add some intermediate joints and connections. For instance, a size of $n\times n$ 2D grid points or $n\times n\times n$ 3D grid points and their dense connections, where $n$ is a user specified parameter. The default value of $n$ is the number of joints specified in the functional specification. This is quite similar to the GSM, the difference is that the number of new joints that we add is usually much less.

4.3. Local Topology Operations

We use the following local topology operations:

• •

Removing the bars with cross-section areas less than a small threshold $\epsilon_{1}$.

• •

Removing joints without any attached bars.

• •

Merging joints that are closer to each other than a small threshold $\epsilon_{2}$.

• •

Removing intermediate joints with valence two, as shown in Figure 6(a).

• •

Deleting the longest bar of a long narrow triangle as shown in Figure 6(b).

• •

Adding a new joint for each pair of intersecting bars. Split this pair of bars into four new bars and connect them at the new joint as shown in Figure 6(c).

• •

Fixing non-boundary T-junctions by adding a bar and a new joint connecting the new bar and the original truss as shown in Figure 6(d). The new joint is created at the point closest to the extension of the existing bar creating the T-junction.

These operations change the local topology and update joint positions.

4.4. Global Topology Operation — Subdivision

The main idea of global topology refinement is to add joints and bars to the truss to be able to reduce its volume after geometry optimization. While previous work, e.g., [Martinez et al., 2007], also proposes to add joints and bars, they add only one joint at a time by testing a large number of candidate locations. This results in a very expensive algorithm. By contrast, we propose to add new joints and bars based in a systematic manner. Our algorithm is inspired by two observations. First, Michell’s theory [Michell, 1904] concludes that the minimum-weight truss should follow two families of continues curves which are orthogonal to each other, one in tension and one in compression. Second, an interesting aspect of truss design is that trusses with more bars can often be lighter than trusses with fewer bars, as more degrees of freedom are provided to approximate an analytical limit. Our algorithm refines the discrete equivalent of such families of curves by subdivision in an efficient and coordinated manner. Most importantly, we insert multiple bars in one step.

We first calculate a pair of tension-compression directions at each joint. As shown in Figure 7(a), for each joint, we separately average the bars connected with this joint according to their force signs (+ for compression (red) and - for tension (blue)) with their force magnitudes as weights. Figure 7(b) shows the calculated nearly-orthogonal directions on a truss. Then we calculate the new joints for the bars which are estimated to be split. Take the bar in Figure 7(c) for example, we know the coordinates of its two ends, $p_{i}$ and $p_{j}$, and the tension-compression directions at its two ends. For a bar in tension, we calculate the two directions, $v_{i}$ and $v_{j}$, which are orthogonal to the compression directions (red) at its two ends. Using $(p_{i},p_{j},v_{i},v_{j})$, we calculate a Bézier curve and set the mid-point of the curve as the new joint. For a bar in compression, we follow a similar procedure.

The purpose of truss subdivision is to improve the orthogonality of the bars in tension and in compression. Given an initial coarse truss, we know its geometry, topology, and the axial force of each bar. Consider the truss as a graph, we extract triangles and quadrilaterals. The triangles are usually formed by bars connected to joints specified in the functional specification. For a triangle as shown in Figure 8, lower row, we add a new joint for the bar whose force sign is different from the other two and connect the new joint with the opposite joint. A quadrilateral is subdivided if it has two non-adjacent bars in compression and the other two non-adjacent bars in tension. As shown in Figure 8 upper row, we add four new joints for its four bars and one more joint at its face center initialized as the average of the previous four, and connect the face-center joint with each edge-middle joint. In the subdivided truss, we remove each bar where a new joint is added, and connect its two ends with the new joint as shown in Figure 8. Figure 9 show the results of different levels of subdivisions using the functional specification in Figure 4(a).

Although the above illustrations are for 2D cases, we can use the same subdivision strategies for 3D trusses. We extract all triangles and quads and test if they should be subdivided. In 3D, we use the same conditions as in 2D (see Fig. 8).

5. Alternating LP

In this section, we introduce the alternating linear programming (ALP) step which optimizes joint positions and cross-section areas of bars for a given topology. ALP serves as the backbone of the proposed approach. Both the local and global topology operations in the previous section are based on optimization of joint positions and cross-section areas. Directly optimizing joint positions and axial forces is a highly nonlinear problem. Therefore, we split the problem into two linear problems. In Algorithm a, we solve force densities alone without changing the joint positions by the ground structure method. In Algorithm b, we update joint positions and force densities jointly based on results from Algorithm a.

5.1. Algorithm a: The Ground Structure Method

Let us first recall the basic plastic formulation of the ground structure method  [Zegard and Paulino, 2014], which solves a continuous linear programming problem to minimize the total volume of material under the premise of force balance with feasible axial forces:

[TABLE]

where each scalar $a_{i}$ is the cross-section area of the $i$-th bar. $\mathbf{B}^{T}$ is the nodal equilibrium matrix, built from the directional cosines of the bars. More details about this matrix are given in the additional materials. ${\mathbf{s}}$ is a vector with the internal (axial) force for all bars, and $|E|$ is the number of bars. ${\mathbf{f}}$ is a vector of the external force for all joints. The internal force $s_{i}$ should be within the range of admissible axial forces $[-\sigma_{T}a_{i},\sigma_{C}a_{i}]$. Here, we assume that the maximal compressive and tensile strains are the same, $\sigma_{C}$=$\sigma_{T}$=$\sigma$, which is a constant value. We set $\sigma=1$ in the formulation. The inequality constraints, 1b and 1c are equivalent to $a_{i}\geq|s_{i}|$. As the length of each bar, $l_{i}$, is positive, $l_{i}>0$, the objective function requires the cross-section area of each bar, $a_{i}$, to be its smallest permissible value, just enough to support the actual axial force of that bar. Then, we have $a_{i}=|s_{i}|$ and the following formulation.

[TABLE]

The above formulation is equivalent when we use force densities $w_{i}=s_{i}/l_{i}$ as variables instead of axis forces $s_{i}$. The new formulation is transformed to:

[TABLE]

Here, in Equation 3a, the matrix $\mathbf{C}$ is a simpler expression than $\mathbf{B}$ because its elements are linear combinations of joint positions. More details about the matrix $\mathbf{C}$ are given in the additional materials.

5.2. Algorithm b: Relocation of Joints

In Algorithm a, the joint positions are assumed to be fixed and the axial forces are the only variables. To further reduce the total volume of material, we complement it with Algorithm b and calculate the displacements of joints to leverage more degrees of freedom. We assume the initial values of force densities, ${\mathbf{w}}$, are known by solving an LP problem in Equation 1 and set the difference of joint positions, ${\mathbf{u}}$, and the difference of force densities of bars, $\Delta{\mathbf{w}}$, as variables. By directly rewriting Equation 3, we have

[TABLE]

Here, we assume that the change of force densities, $\Delta{\mathbf{w}}$, is small and that the signs of force densities remain the same, $\mathbf{sgn}(w_{i}+\Delta w_{i})=\mathbf{sgn}(w_{i})$. As Algorithm b is applied after Algorithm a, the values, $l_{i}$, $w_{i}$, $\mathbf{sgn}(w_{i})$, and $\mathbf{C}$, are all known.

To simplify the problem which has a cubic objective function and quadratic constraints, our goal is to approximate the above formulation with a linear programming problem and solve it in a sequential manner in conjunction with Algorithm a. By expanding the objective function, we have $(l_{i}+\Delta l_{i})^{2}(w_{i}+\Delta w_{i})=(l_{i}^{2}w_{i}+2l_{i}\Delta l_{i}w_{i}+\Delta l_{i}^{2}w_{i}+l_{i}^{2}\Delta w_{i}+2l_{i}\Delta l_{i}\Delta w_{i}+\Delta l_{i}^{2}\Delta w_{i})\approx(l_{i}^{2}w_{i}+2l_{i}\Delta l_{i}w_{i}+l_{i}^{2}\Delta w_{i})$. Here, we remove the higher order terms, and use the fact that $l_{i}^{2}w_{i}$ is constant. The objective function is approximated by $\sum_{i=1}^{|E|}\mathbf{sgn}(w_{i})(2l_{i}\Delta l_{i}w_{i}+l_{i}^{2}\Delta w_{i})$. As shown in Figure 10, $\Delta l_{i}\approx{\mathbf{d}}_{i}\cdot({\mathbf{u}}_{i2}-{\mathbf{u}}_{i1})$, where ${\mathbf{d}}_{i}$ is the unit direction vector of the i-th bar connecting the joints $i1$ and $i2$. The objective function is linear with ${\mathbf{u}}_{j}$ and $\Delta w_{i}$ as variables, where $j=1,\ldots,|V|$ and $i=1,\ldots,|E|$.

The force equilibrium constraint, $(\mathbf{C}+\Delta\mathbf{C})^{T}({\mathbf{w}}+\Delta{\mathbf{w}})=-{\mathbf{f}}$, is equivalent to $\mathbf{C}^{T}\Delta{\mathbf{w}}+\Delta\mathbf{C}^{T}{\mathbf{w}}+\Delta\mathbf{C}^{T}\Delta{\mathbf{w}}=\mathbf{0}$ because $\mathbf{C}^{T}{\mathbf{w}}=-{\mathbf{f}}$ is ensured by Algorithm a. Here we remove the small higher order term $\Delta\mathbf{C}^{T}\Delta{\mathbf{w}}$. Then the force balance constraint is linearized as $\mathbf{C}^{T}\Delta{\mathbf{w}}+\Delta\mathbf{C}^{T}{\mathbf{w}}=\mathbf{0}$. The matrix $\mathbf{C}$ and the vector ${\mathbf{w}}$ are known from Algorithm a and the elements in matrix $\Delta\mathbf{C}$ are linear combinations of nodal variations ${\mathbf{u}}_{j}=(u_{jx},u_{jy},u_{jz})$. Then the constraint is also linear with respect to ${\mathbf{u}}_{j}$ and $\Delta w_{i}$.

Finally, the formulation for Algorithm b is written as

[TABLE]

where $\lambda_{j}$ and $\delta_{i}$ are the bounds of the variables $u_{i}$ and $\Delta w_{i}$. In each iteration, we set small values for these bounds, e.g., $\delta_{i}=0.1|w_{i}|$ and $\lambda_{i}=0.1\bar{l}$, where $\bar{l}$ is the average length of all bars.

5.3. Alternating Scheme

The above two algorithms are formulated as two LP problems in Equation 1 and Equation 5. The inputs to Algorithm a are the joint positions, ${\mathbf{p}}$, and the functional specification such as the external forces, $\mathbf{LOAD}$, and the supporting points, $\mathbf{SUPP}$, and the outputs are the force densities of bars, ${\mathbf{w}}$. The algorithm a is written as $[{\mathbf{w}},V]=\mathbf{ALGa}({\mathbf{p}},\mathbf{LOAD},\mathbf{SUPP})$, where $V$ is the total volume of materials. The inputs of Algorithm b are the initial force densities, ${\mathbf{w}}$, the initial joint positions, ${\mathbf{p}}$, and the same functional specification. The outputs are the changing values of joint positions and force densities. Then, Algorithm b is written as $[{\mathbf{u}},\Delta{\mathbf{w}}]=\mathbf{ALGb}({\mathbf{p}},{\mathbf{w}},\mathbf{LOAD},\mathbf{SUPP})$. In the whole algorithm, we organize them in an alternating way as shown in Algorithm 1. $N_{\mathrm{max}}$ is the maximum iteration number and $S_{\mathrm{max}}$ is the maximum line search step.

In Figure 11, we show the effectiveness of ALP for truss optimization with three sets of load specifications involving torques.

6. Extensions

The generality of our formulation allows extensions of our approach for a broad range of scenarios and applications, such as tackling multiple load specifications, respecting stability analysis, and incorporating project-specific fabrication constraints.

6.1. Multiple Load Specifications

For the input specification with multiple sets of external forces, the ALP algorithm is adjusted accordingly. Static equilibrium is required for each set of external loads with generally different internal forces. Thus, assuming trusses are required to withstand $K$ sets of external forces ${{\mathbf{f}}}^{1}$, … , ${{\mathbf{f}}}^{K}$, the formulation of Algorithm a is rewritten as:

[TABLE]

where $k=1,\ldots,K$. Force equilibrium constraints similar to Equation 1 are required for each set of external forces. It is also worth noting that this set of equations is sufficient to guarantee force equilibrium in response to linear interpolation of the sets of specified external forces. As each bar needs to support the maximal axial forces from each set of the reaction forces, the cross-section of the $i$-th bar, $a_{i}=max\{|{s_{i}}^{1}|,...,|{s_{i}}^{K}|\}$. We denote $m_{i}\in\{1,...,k\}$ as the set index for which the $i$-th bar attains its maximal axial force, $|s_{i}^{m_{i}}|=max\{|{s_{i}}^{1}|,...,|{s_{i}}^{K}|\}$. As those indices, $m_{i}$, could be easily found from the result of Algorithm a, corresponding cross sections follow directly, $a_{i}=|s_{i}^{m_{i}}|$. Similarly, by defining $w_{i}^{m_{i}}=s_{i}^{m_{i}}/l_{i}$, the formulation of Algorithm b for the multiple-load case can be written as:

[TABLE]

where $i=1,\ldots,|E|$, $k=1,\ldots,K$, and $j=1,\ldots,|V|$. Here, only the nodal variations, $u_{i}$, and the change of force densities, $\Delta w_{i}^{k}$, are variables, others are known from Algorithm a. Using the same alternating scheme in Section 5.3 by adjusting the Algorithm a and b accordingly, our method can tackle cases of multiple-load input specifications. See Figure 12 for an example.

6.2. Stability and Buckling

Truss stability, in particular, the external stability, has been investigated extensively for its importance in keeping its rigid shape. By comparing the number of degrees of freedom against constraint counts, a condition for external stability is given in  [Kassimali, 2011] as

[TABLE]

With $|E|$ being the number of bars and $r$ being the count of reaction components, the sum on the left-hand-side includes the number of both internal (e.g., fixed edge lengths) and external (e.g., supported vertices) constraints. The right-hand-side counts the total number of degrees of freedom, as a product of the number of joints, $|V|$, and the number of degrees of freedom for each joint, $d$, which is 2 for plane trusses and 3 for spatial trusses. A truss is externally stable if the condition is met and unstable otherwise. There are multiple approaches to integrate stability analysis to our framework, e.g., through a post-processing step checking the aforementioned condition and adding necessary auxiliary bars with a minimum cross-section to increase the number of $|E|$ on the left side and thus achieving external stability. Figure 13 provides an example. Similarly to stability analysis, buckling analysis through simulation packages like ANSYS are applicable as a feedback procedure to improve the truss structural properties of interest in an interactive manner.

6.3. Fabrication Constraints

Besides general engineering settings, our system is able to incorporate project-specific fabrication constraints through minor adjustments or complementary post-processing procedures. For projects requiring a limited number of parts, bounding the number of subdivision steps controls the amount of bars and joints. Meanwhile, while our algorithm reduces the total volume of material consumption based on a continuous selection of cross-sections in the default setting, for engineering practices which demand viable bars to be chosen from a set of predefined category or a discrete number of types, our method can be combined with other published algorithms. We show an example in limiting the types of cross sections by applying  [Jiang et al., 2017] directly as a post-processing procedure in Figure 14.

7. Results

In this section, we illustrate truss designs using our framework for different types of input specifications and compare the results with state-of-the-art methods on selected benchmark design problems.

7.1. Example Designs

We show the results of our method for different types of functional specifications. We present 2D truss designs with a parallel equilibrium force system in Figure 22, a concurrent equilibrium force system in Figure 23, and a non-concurrent equilibrium force system in Figure 24. We illustrate examples of designs for the same input external forces and supporting joints but different design regions in Figure 25. For 3D trusses, we demonstrate the results for input of a parallel equilibrium force system in Figure 27, a concurrent equilibrium force system in Figure 28, and a non-concurrent equilibrium force system in Figure 29. In addition, we show examples based on real functional requirements such as a 2D bike frame in Figure 26, a 3D cantilever in Figure 30, and a 3D bridge in Figure 1.

7.2. Quantitative Evaluation

Our framework is implemented in Matlab R2016 on a workstation with an Intel Xeon X5550 2.67 GHz processor. We use Mosek [ApS, 2015] as the solver for linear programming. For the ALP algorithm, we set the maximum iteration number $N_{\mathrm{max}}$=500 and the maximum line search step $S_{\mathrm{max}}$=10. For the number of grid points in the initialization, $n$, we use the default value, the number of joints specified in the functional specification. We set the number of max iteration of Phase 1, $P_{\mathrm{max1}}=5$. The number of subdivisions in Phase 2, $P_{\mathrm{max2}}$, controls the trade off between number of bars and truss weight and is set depending on the context. We set the thresholds $\epsilon_{1}$ and $\epsilon_{2}$ in Section 4.3 as $\epsilon_{1}=0.002\bar{a}$, $\epsilon_{2}=0.01\bar{d}$, where $\bar{a}$ is the average cross-section, and $\bar{d}$ is the average value of distances between joints specified in the functional specification. In Table 2, for each optimized truss, we report parameters of trusses such as the number of bars, the total volume of material, and the computation time of coarse truss optimization and different levels of structure refinement.

7.3. Evaluation of the Initialization

To test the sensitivity of our framework to the initialization, we show optimized trusses starting from initializations with different numbers of intermediate joints and bars (before subdivision operations are conducted). Figure 15 shows that adding different intermediate joints and bars results in different trusses. However, these trusses have similar structure and they all function as robust discrete approximations of the optimal truss for further subdivision.

7.4. Comparisons

7.4.1. Comparing ALP to alternative solvers

We compare the proposed ALP algorithm with three alternative optimization methods: sequential quadratic programming (SQP), gradient descent (GD), and the method of Jiang et al. [Jiang et al., 2017]. The details of the implementation of SQP and GDM are given in the appendix. Jiang et al.’s framework contains multiple parts and we only compare to the solver that is comparable to our proposed ALP algorithm (The objective function in Jiang et al. is different though). For various input specifications, we compare the quality (volume) of the output. Fig. 16 and Table 1 show the results. The results demonstrate that our proposed optimization converges to a better solution than competing approaches, especially GD struggles to find meaningful solutions. We noticed that competing algorithms can be improved when using the ground structure method (GSM) for initialization. In the table, we use a * to indicate the version of SQP and GD that has been initialized using GSM. In our comparison, our method generates the best results for all input specifications. SQP* can match our result in some cases. We can also observe that our method is fast and better scales to larger inputs than SQP. The reason for the poor performance of Jiang et al. [Jiang et al., 2017] is that they use soft constraints instead of hard constraints in the formulation.

7.4.2. Comparisons with Previous Numerical Methods

We compare the performance of our method with several previous methods using the functional specifications, solutions, and running times provided in their papers. We denote the methods in [Descamps and Coelho, 2013], [Gilbert and Tyas, 2003], [He and Gilbert, 2015], and [Sokół, 2017] as D2013, G2003, H2015 and S2017, respectively. In Figure 1 and Figure 17, we compare our method with [Descamps and Coelho, 2013] for a 2D and a 3D bridge design problem. In Figure 18 and 19, we compare our method with [He and Gilbert, 2015] and [Gilbert and Tyas, 2003] on a benchmark problem—Hemp cantilever design. For 3D truss optimization, we compare our method with [Sokół, 2017] on a simple 3D model in Figure 20. The comparison is shown in Table 3. We can observe that our method is orders of magnitude faster, even though we can achieve a lower volume than previous work. Unfortunately, these methods do not have code or executables publicly available, so we cannot test them on the same machine. However, it seems unlikely that this significant difference in running time can be overcome by slightly faster hardware.

7.4.3. Comparisons with Analytical Solutions

For the cantilever design problems in Figure 18 and 19, the total volumes of analytical solutions are 4.498115 and 4.232168. The volumes of our discrete designs for these two cases are 4.498635 and 4.3223, which are closer to the analytical solutions than the previous work, [Gilbert and Tyas, 2003] and [He and Gilbert, 2015]. For the three force problem, we compare our result with the analytical solution presented in  [Sokół and Lewiński, 2010]. In Figure 21, we show the computed discrete truss (130 bars, volume: 6.838) on the right which is visually similar to the analytical solution on the left (volume: 6.831).

Discussion and Limitations.

Compared with previous work, the geometry optimization in terms of both the axial forces and joint positions through two linear programming problems alternatively applied in ALP provides more degrees of freedom compared with the original formulation of the ground structure method which solves a single linear programming problem. Splitting a highly nonlinear programming problem into two linear ones also attains better efficiency compared with the original nonlinear formulation. Moreover, the two categories of topological operations, local and global, complementarily allow both flexible yet stable topology changes and manipulation. Most importantly, the subdivision approach, which has been overlooked by previous work, is a natural choice for topology refinement from coarse to fine, which creates valid topologies at different levels both efficiently and robustly. Despite the efficiency and efficacy, our optimization framework cannot guarantee a global optimum. However, in simple special cases where the analytical optimum is known, we observe that the method almost reaches the known global optimum. As our method is based on subdivision of edge-networks on surfaces, the optimal trusses in 3D also need to constitute sheets of surfaces. For general 3D specifications, the subdivision approach requires an initial surface-like structure, or a structure that consists of multiple sheets, which is a challenging problem for future work. The proposed ALP algorithm will work for general 3D structures, however.

8. Conclusions and Future Work

We present a method for the design of optimal trusses satisfying functional specifications with minimized material consumption. The core components of the proposed approach include an alternating linear programming formulation for geometry optimization and two sets of topological operations. The subdivision scheme inspired by Michell’s theoretical studies utilized in the global topology refinement step plays a crucial role for the efficiency and efficacy of the proposed approach. The performance of our framework is validated by comparisons with multiple previous studies in different scenarios , which indicate that our method creates trusses with smaller volumes and is over two orders of magnitude faster in terms of computational speed. For future work, it would be exciting to study dynamic structures created by trusses with movable parts that have multiple configurations, e.g., robotic and mechatronic systems with trusses. Moreover, a better theoretical understanding of optimal trusses in 3D would be inspiring for the entire field.

Appendix

8.1. Equilibrium Force Systems

The input functional specification is essentially an equilibrium force system. For a system with forces $F_{1}$,$F_{2}$,…,$F_{n}$ acting at points $P_{1}$,$P_{2}$,…,$P_{n}$, there are two constraints of equilibrium: the balance of forces, $\sum_{i=1}^{n}F_{i}=\vec{0}$, and the balance of torques, $\sum_{i=1}^{n}{F_{i}\times P_{i}}=\vec{0}$. In mechanical engineering, the equilibrium force systems are usually classified into three types: parallel force systems, concurrent force systems, and non-concurrent force systems as shown in Figure 32.

8.2. SQP

First, lets recall the formulation in Eq.3. In this formulation, joint positions are considered as known values. The objective function and constraints in the above formulation are both linear. If the joint positions ${\mathbf{u}}_{j},j=1,...,|V|$ are considered as design variables, the formulation is a nonlinear programming with nonlinear constraints. To solve such a problem by SQP, we reformulate the problem similarly as in  [Descamps and Coelho, 2013]. We introduce slack variables, $w_{i}^{+}$ and $w_{i}^{-}$, which represent the magnitudes of compressive and tensile forces respectively: $w_{i}^{+}=\max(0,w_{i})$, $w_{i}^{-}=\max(0,-w_{i})$. Then, we have $w_{i}=w_{i}^{+}-w_{i}^{-}$, and $|w_{i}|=w_{i}^{+}+w_{i}^{-}$. Here we denote the vectors consisted of all $w_{i}^{+}$ and $w_{i}^{+}$ as ${\mathbf{w}}^{+}$ and ${\mathbf{w}}^{-}$. We assume the $i$th bar is connected by joints $i1$ and $i2$, then $l^{2}_{i}=({\mathbf{u}}_{i1}-{\mathbf{u}}_{i2})^{2}$. The elements in matrix $\mathbf{C}$ are linear combinations of joint positions. The whole problem is reformulated as

[TABLE]

The algorithm in  [Descamps and Coelho, 2013] used the above formulation by SQP method. We implement the algorithm using the $fmincon$ routine in Matlab.

8.3. The Gradient Decent Method

To use the GDM, we first transfer the constrained optimization problem in Eq. 6 into an unconstrained formulation by adding an extra penalty term into the objective function as

[TABLE]

where

[TABLE]

and

[TABLE]

This objective function is defined similarly in  [Jiang et al., 2017]. The constraint of $w_{i}^{+}\geq 0$ and $w_{i}^{-}\geq 0$ are avoided by replacing $w_{i}^{+}$ and $w_{i}^{-}$ with $(\varrho_{i}^{+})^{2}$ and $(\varrho_{i}^{+})^{2}$, where $\varrho_{i}^{+}$ and $\varrho_{i}^{-}$ are new variables without any constraints. Usually, the weight of penalty term $\lambda_{2}$ should be much larger than $\lambda_{1}$. The above object function is differentiable. We do the comparison using standard gradient decent method.

Additional Materials

8.4. Special Cases

The theorem of Maxwell(1872) states that $\sum{f_{i}l_{i}}=C$, where $C$ is the same constant for any statically admissible truss layout considering a given set of external forces. As a consequence, if it is found possible to design a structure all of whose members are in tension, or alternatively compression, then the optimum design has been achieved, because $\sum{\left|f_{i}\right|l_{i}}=\left|C\right|$. The above statement assumes the input force system is precisely given. For example, a self-equilibrium force system with no supporting points, such as the 3-force problem in 2D and 4-force problem in 3D. However, it is not correct when some supporting points are given without prescribing their reaction forces. For example, as shown in Figure 33, we can find a structure of two bars with compression on the left, but it is not an optimal structure because reaction forces at the boundary could be changed to achieve a better results as show on the right.

According to the theorem of Maxwell, we can easily find some special cases. For example, given a concurrent force systems with forces $F_{1}$,$F_{2}$,…,$F_{n}$ through a common point $O$ and their corresponding action points $P_{1}$,$P_{2}$,…,$P_{n}$, we have $F_{i}=\lambda_{i}\overrightarrow{OP_{i}}$. If all the $\lambda_{i}$ have the same sign, either positive or negative, then the truss connect all $P_{i}$ and $O$ is an optimum design. Note that the optimum for these special boundary conditions is usually not unique. Figure 34 shows multiple optimum designs for the same boundary condition where 3 forces are self-equilibrium in 2D plane. From a design standpoint, the non-uniqueness provide a flexibility of topology and geometry design for these special boundary conditions.

To systematically investigate the boundary condition under which a design with all bars in compression or tension can be found is interesting but out of the scope of this paper. Our numerical computational method can be use to detect such kind of boundary condition in the initial truss generation stage. Under such kind of special boundary conditions, further subdivision is not necessary.

8.5. Nodal Equilibrium Matrix

The nodal equilibrium matrix, $\mathbf{B}$, transforms the magnitudes of internal forces, ${\mathbf{s}}=(s_{1},s_{2},\dots,s_{|E|}$), at each bar, into forces along the $x$, $y$ and $z$ axes at each joint. If a truss is in equilibrium, then each joint must be in equilibrium. The equilibrium equation is written as

[TABLE]

Let’s take a look at the equilibrium equations of a particular $i$-th joint with a valence of $n_{i}$, connected to nodes $p_{i,j},j=1,\dots,n_{i}$ through bars $b_{i,j},j=1,\dots,n_{i}$ respectively. Then the equilibrium equations at the joint $i$ are

[TABLE]

where $l_{b_{i,j}}=\sqrt{(x_{p_{i,j}}-x_{i})^{2}+(y_{p_{i,j}}-y_{i})^{2}+(z_{p_{i,j}}-z_{i})^{2}}$ is the length of bar $b_{i,j}$. By writing all the equilibrium equations in a matrix form, we get the Equation 8, where the element in the matrix $B$ are linear combinations of the directional cosines of the bars.

If we use force densities $w_{i}=s_{i}/l_{i}$ as variables, the above equations are transformed to simpler forms:

[TABLE]

The above equilibrium equations could be written in a similar matrix form as

[TABLE]

All the elements in the matrix $\mathbf{C}$ are linear combinations of the joint coordinates. Similarly, the elements in the matrix $\Delta\mathbf{C}$ in Equation 5a are linear combinations of nodal variations.