Computational Design of Lightweight Trusses
Caigui Jiang, Chengcheng Tang, Hans-Peter Seidel, Renjie Chen, Peter Wonka
CComputational Design of Lightweight Trusses
CAIGUI JIANG,
MPI for Informatics, Germany and KAUST, Saudi Arabia
CHENGCHENG TANG,
Stanford University, USA
HANS-PETER SEIDEL,
MPI for Informatics, Germany
RENJIE CHEN,
MPI for Informatics, Germany
PETER WONKA,
KAUST, Saudi Arabia ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) Fig. 1. A bridge design problem from [Descamps and Coelho 2013]. The input is an initial structure of 258-bars with supporting points (red) and two sets offorces: (a) downward loads of magnitude 1 and (b) horizontal loads of magnitude 0.2 perpendicular to the bridge’s main direction. The result in [Descamps andCoelho 2013] has a total volume of 408.8 and took over 1000s to compute. (c) Our optimal geometry and topology viewed from different angles. The totalvolume is 333.4 and the running time is less than 10s. (d) Further structure refinement based on (c) to achieve a total volume of 331.9.
Trusses are load-carrying light-weight structures consisting of bars con-nected at joints ubiquitously applied in a variety of engineering scenarios.Designing optimal trusses that satisfy functional specifications with a mini-mal amount of material has interested both theoreticians and practitionersfor more than a century. In this paper, we introduce two main ideas to im-prove upon the state of the art. First, we formulate an alternating linearprogramming problem for geometry optimization. Second, we introduce twosets of complementary topological operations, including a novel subdivisionscheme for global topology refinement inspired by Michell’s famed theo-retical study. Based on these two ideas, we build an efficient computationalframework for the design of lightweight trusses. We illustrate our frameworkwith a variety of functional specifications and extensions. We show that ourmethod achieves trusses with smaller volumes and is over two orders ofmagnitude faster compared with recent state-of-the-art approaches.CCS Concepts: •
Computing methodologies → Shape modeling .Additional Key Words and Phrases: truss, topology optimization, geometry
Trusses are crucial and fundamental structures in multiple mod-ern engineering domains. They consist of bar elements that areconnected by pin joints. Because of their efficiency and lightweightnature, trusses see considerable amount of usage in industrial designand architectural construction, e.g., for support structures of build-ings, bridges, transmission towers, or even domes in playgrounds.Designing a lightweight truss typically starts with a functionalspecification, e.g., in the form of external forces that the structurehas to withstand. The design problem can then be formulated as an optimization problem to determine the geometry, topology, and thecross-sections of the truss. In other words, we have to find answersto the following questions: Where to put the intermediate joints?How to connect the joints with bars? What are the cross-sectionareas of the bars? These tasks are notoriously challenging becausethe optimization of geometry, topology, and cross-sections is inter-related, and there exists an infinite number of possible topologieswhich are difficult to classify and quantify. Even for a simple case,the optimal topology is not intuitive. As shown in Figure 2, to sup-port two pairs of opposing forces lying on two straight lines, thesimplest truss on the left with two bars, one in tension (blue) andanother in compression (red), may be intuitively considered as thelightest truss. However, a lighter design with more intermediatejoints and connections can be found as shown in the Figure 2 right.Another simple functional specification problem is called the threeforces problem (3FP) [Chan 1966; Sokół and Lewiński 2010]. 3FP isformulated as follows: find the lightest fully stressed truss transmit-ting three self-equilibrated co-planar forces. Although there are onlythree forces, the problem is still unsolved analytically for generalcases.In this paper, we mainly follow previous work and take functionalspecifications in the form of supporting points and applied forces asinput. Our goal is then to construct a lightweight truss with optimaljoint positions, topology, and cross-sections. As an example shownin Figure 3 left, two supporting points and one external force aregiven as inputs. Our computational method generates the topology a r X i v : . [ c s . G R ] J a n utomatically and optimizes the nodal positions and cross-sectionareas as shown in Figure 3 right. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) Fig. 2. Two truss designs for the same functional specification. Left: astraightforward design with two bars. Right: a more complex and less intu-itive design with less material usage. inputinputinputinputinputinputinputinputinputinputinputinputinputinputinputinputinput outputoutputoutputoutputoutputoutputoutputoutputoutputoutputoutputoutputoutputoutputoutputoutputoutput
Fig. 3. Left: a functional specification including two supporting points andone external force. Right: an optimal truss.
There are two major strategies to tackle this problem in previouswork. The first strategy is to start with a densely-connected structureand to subsequently identify which bars to remove (e.g., the groundstructure method [Dorn 1964] and its variations [Gilbert and Tyas2003; Sokół 2017]). The main limitation of this strategy is its sub-optimality due to its heavy dependence on the initialization as it isextremely unlikely that the optimal structure is a sub-structure ofthe initial one. The other strategy is to start with a sparse structureand to iteratively add new joints and bars. One of the most famousmethods in [Martinez et al. 2007] adds one joint at a time and canonly deal with one single load 2D problem.Even though topology optimization is such a longstanding andfundamental problem in structural engineering, one can identifylarge possible improvements to the current state of the art. First, thesearch space for truss topology is not properly explored by previousalgorithms. We could observe that it is very difficult to optimize thetopology in a single stage. Much better results can be achieved byproceeding in two stages of topology optimization: computing acoarse truss and truss subdvision. Our novel subdivision approachis inspired by Michell’s pioneering theoretical treatment of optimaltruss design in [Michell 1904]. Second, the geometry optimizationused in previous work is not efficient. To tackle this problem wedecompose geometry optimization into alternating linear program-ming formulations to reduce the running time.Our main contributions are as described in the following • We propose two categories of complementary topology opera-tions, local and global. While local operations have been usedin previous work, our global operations based on subdivisionare our original contribution. • We introduce a novel algorithm for geometry optimizationbased on alternating linear programming (ALP) that jointlyoptimizes joint positions and bar cross sections. • Based on these two technical contributions, we build a frame-work for lightweight truss design, a longstanding and impor-tant problem in structural engineering, architecture, graph-ics, and design. Compared with recent state-of-the-art ap-proaches, our method creates trusses with smaller volumes,can handle more complex functional specifications, and isover two orders of magnitude faster.
In recent years, combining geometric modeling together with real-istic engineering considerations, especially static equilibrium andmanufacturability, have attracted the interests of many researchersin the graphics community. Beyond applications in the virtual world[Smith et al. 2002], those previous works enable novel and func-tional designs manufacturable with 3D printing [Wang et al. 2013;Zhou et al. 2013], laser cutting [Martínez et al. 2015], masonry struc-ture [Block and Ochsendorf 2007; de Goes et al. 2013; Liu et al.2013; Panozzo et al. 2013; Tang et al. 2014; Vouga et al. 2012], fortoys [Bächer et al. 2014; Prévost et al. 2013], furniture [Umetani et al.2012; Yao et al. 2017], and architecture [Jiang et al. 2015; Pietroniet 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 designand optimize space structures where only a small set of cross-sectionareas are allowed. Therefore, Jiang et al. [Jiang et al. 2017] computea specialized form of truss, but we focus on the classical problem oftruss design without discrete restrictions on the cross sections. Themain practical difference is that our proposed method can generatea truss from scratch, whereas Jiang et al. relies on a reasonable trussbeing given as input. In our results we also demonstrate that ourproposed optimization technique ALP produces better results onour problem formulation than the geometric optimization techniqueused in [Jiang et al. 2017]. Kilian et al. [Kilian et al. 2017] provide aninteresting geometric understanding of "optimality" of surface-likelightweight structures. Compared with their work, we tackle a prob-lem for common and general trusses in both 2D and 3D, instead offocusing on load-carrying surfaces.The problem of designing a truss with a minimal volume ofmaterial that supports imposed external forces was first studiedin [Michell 1904]. In the milestone paper, Michell proved that anoptimal truss must follow orthogonal networks of lines of maxi-mal and minimal strains in a constant-magnitude strain field. Anoptimal truss is usually called a Michell truss. Following his work,research on the topic of optimal truss can be divided into two cate-gories: exact-analytical formulations and approximate-discretizedformulations.An exact-analytical formulation assumes that the truss is a con-tinuum structure connected by an infinite number of bars withinfinitesimally small cross-sections. In analytical formulations, thetheoretical optimal design is determined exactly through the simul-taneous solution of a system of equations expressing the condi-tions for optimality. The basic principles were establish in [Maxwell1870] and [Michell 1904] and a more general treatment was outlinedin [Hemp 1973] and [Prager and Rozvany 1977]. Recent works on de-riving exact solutions were presented in [Rozvany 1998], [Lewińskiand Rozvany 2007], [Lewiński and Rozvany 2008b], and [Lewiński nd Rozvany 2008a] for a series of benchmark problems. Basically,the analytical solutions are very hard to obtain and only availablefor some special boundary conditions. While they are less practicalin most of the generic scenarios, these solutions could be used asreferences to verify the performance of numerical methods.Discretized numerical formulations are more practical and effi-cient approaches for structural design tasks presented in the realworld. The most influential method is the ground structure method(GSM) which was first proposed in [Dorn 1964]. This method con-sists of generating a fixed grid of joints and adding bars in some orall of the possible connections between the joints as potential struc-tural or vanishing bars. The optimized structure for the imposedfunctional specification is found using the cross-section areas asdesign variables, and the whole problem is formulated as a linearprogramming problem. Its optimal topology is achieved by elimi-nating the zero-area cross sections. The ground structure methodhas been recently improved in [Gilbert and Tyas 2003] and [Sokół2017].Besides GSM, some other numerical methods are proposed re-cently, such as the method in [Martinez et al. 2007], carries outgeometry optimization in conjunction with a heuristic ‘joint adding’algorithm, generating an increasingly complex truss structure froma relatively simple initial layout. However, this algorithm can onlyadd one joint per time and only works for single load cases. Anefficient algorithm proposed by He and Gibbert [2015] combineslayout optimization with geometry optimization. Similar to GSM,its layout optimization starts from a densely connected truss andis formulated by a linear programming problem, and its geome-try optimization is formulated by a non-linear optimization as apost-processing step. Our framework has the following major components: • Functional specification (C1). The input to our frameworkis the functional specifications including the external forcesand supporting points together with a set of structural con-straints, e.g., design regions, geometric obstacles, and materialproperties. (See Section 4.1) • Initialization (C2). To obtain an initial truss, we create a gridof intermediate joints and densely connect them. This grid islocated inside the design region and its size is proportionalto the bounding box of the points in the input specification.(See Section 4.2) • Local topology operations (C3). We locally manipulate thetopology through some geometry operations such as remov-ing bars with vanishing cross-section areas and joints withoutany connection, merging close joints, etc. (See Section 4.3) • Global topology refinement using subdivision (C4). Usean optimized coarse truss as input, we further refine the trussthrough subdivision. (See Section 4.4 ) • Geometric optimization using ALP (C5). Given a fixedtopology, we propose an alternating linear programming al-gorithm (ALP) to reduce the total volume of the truss byadjusting the joint positions and cross-section areas of bars. This algorithm is an essential component and its details areintroduced in Section 5.
Pipeline Overview.
The individual components work together asfollows: We start from an input specification to compute an initialtruss. Then, we proceed in two phases. In the first phase, coarse trussoptimization, we interleave geometric optimization (C5) with localtopology operations (C3). In the second phase, structure refinementthrough subdivision, we interleave global topology refinement (C4)with geometric optimization (C5). See Figure 4 for an overview ofour framework. We discuss each individual component and theoverall framework in detail in the next sections.
We provide a framework for the computational design of light-weight trusses. In this section, we describe the input specification,the initialization, local topology operations, and global topologyrefinement.
The input to our framework is the functional specification includingthe external forces and supporting points together with a set ofstructural constraints, e.g., design regions, geometric obstacles, andmaterial properties. Throughout the paper, we visualize supportingjoints as red dots, joints with active forces as blue dots, and inter-mediate joints as yellow dots. We also visualize bars in tension inblue and bars in compression in red. In addition, the thickness ofthe bars is visualized in proportional to the computed cross-sectionareas. 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.
We build on previous work to compute an initial truss. There aretwo approaches to tackle this problem. One simply adds connectionsbetween provided joints in the functional specification (supportingjoints and joints with active forces). For example, as shown in Figure5 left, two bars connecting the joints with active forces (blue) andthe supporting joints (red) are set as an initial truss. In some cases,this initialization is too simple to construct an equilibrium forcesystem. Another method adds a grid of intermediate points overthe design region and densely connects them as shown in Figure 5right. The increasing method such as the work in [Martinez et al.2007] used the first initialization. The GSM usually uses the latterone with a large number of intermediate joints. ocal Topology Operations (C3)Geometry Optimization (C5) Subdivision (C4)Initialization Input specification
Geometry Optimization (C5)Local Topology Operations (C3) Geometry Optimization (C5)Subdivision (C4)Geometry Optimization (C5) Ouput Truss ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e ) Fig. 4. Framework overview: Based on an input functional specification (a), our system creates an initial truss (b). Then, we proceed in two phases, coarsetruss optimization (c) and structure refinement through subdivision (d). In the first phase, we interleave geometric optimization using ALP with local topologyoperations. In phase 2, we interleave subdivision with geometry operation using ALP. The output truss is shown in (e).Fig. 5. Two kinds of initializations. Left: connect force application pointsand supporting points. Right: a densely connected initial structure for GSM.
In our framework, we first add some intermediate joints andconnections. For instance, a size of n × n
2D grid points or n × n × n n is a user specifiedparameter. The default value of n is the number of joints specifiedin the functional specification. This is quite similar to the GSM, thedifference is that the number of new joints that we add is usuallymuch less. We use the following local topology operations: • Removing the bars with cross-section areas less than a smallthreshold ϵ . • Removing joints without any attached bars. • Merging joints that are closer to each other than a smallthreshold ϵ . • Removing intermediate joints with valence two, as shown inFigure 6(a). • Deleting the longest bar of a long narrow triangle as shownin Figure 6(b). • Adding a new joint for each pair of intersecting bars. Splitthis pair of bars into four new bars and connect them at thenew joint as shown in Figure 6(c). • Fixing non-boundary T-junctions by adding a bar and a newjoint connecting the new bar and the original truss as shownin Figure 6(d). The new joint is created at the point closest tothe extension of the existing bar creating the T-junction. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) Fig. 6. Four types of local topology operations: (a) Delete a joint of valencetwo. (b) Remove an ill-shaped narrow triangle. (c) Add an additional jointfor a pair of intersecting bars.(d) Fix non-boundary T-junction.
These operations change the local topology and update jointpositions.
The main idea of global topology refinement is to add joints andbars to the truss to be able to reduce its volume after geometryoptimization. While previous work, e.g., [Martinez et al. 2007], alsoproposes to add joints and bars, they add only one joint at a time bytesting a large number of candidate locations. This results in a veryexpensive algorithm. By contrast, we propose to add new joints andbars based in a systematic manner. Our algorithm is inspired by twoobservations. First, Michell’s theory [Michell 1904] concludes thatthe minimum-weight truss should follow two families of continuescurves which are orthogonal to each other, one in tension and onein compression. Second, an interesting aspect of truss design is thattrusses with more bars can often be lighter than trusses with fewerbars, as more degrees of freedom are provided to approximate ananalytical limit. Our algorithm refines the discrete equivalent ofsuch families of curves by subdivision in an efficient and coordinatedmanner. Most importantly, we insert multiple bars in one step.We first calculate a pair of tension-compression directions at eachjoint. As shown in Figure 7(a), for each joint, we separately averagethe bars connected with this joint according to their force signs(+ for compression (red) and - for tension (blue)) with their forcemagnitudes as weights. Figure 7(b) shows the calculated nearly-orthogonal directions on a truss. Then we calculate the new jointsfor the bars which are estimated to be split. Take the bar in Figure7(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 abar in tension, we calculate the two directions, v i and v j , which areorthogonal 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 ofthe curve as the new joint. For a bar in compression, we follow asimilar procedure. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) p i p i p i p i p i p i p i p i p i p i p i p i p i p i p i p i p i p j p j p j p j p j p j p j p j p j p j p j p j p j p j p j p j p j v i v i v i v i v i v i v i v i v i v i v i v i v i v i v i v i v i v j v j v j v j v j v j v j v j v j v j v j v j v j v j v j v j v j Fig. 7. Subdivision of a truss: (a) construction of compression-tensiondirections at each joint; (b) a compression-tension field for a truss; (c) thestrategy to calculated an edge mid-point.
The purpose of truss subdivision is to improve the orthogonalityof the bars in tension and in compression. Given an initial coarsetruss, 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 specifiedin 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 differentfrom the other two and connect the new joint with the oppositejoint. A quadrilateral is subdivided if it has two non-adjacent barsin compression and the other two non-adjacent bars in tension. Asshown in Figure 8 upper row, we add four new joints for its fourbars and one more joint at its face center initialized as the average ofthe previous four, and connect the face-center joint with each edge-middle joint. In the subdivided truss, we remove each bar where anew joint is added, and connect its two ends with the new joint asshown in Figure 8. Figure 9 show the results of different levels ofsubdivisions using the functional specification in Figure 4(a).
Fig. 8. Truss subdivision strategies for quadrilaterals (top) and triangles(bottom).
Although the above illustrations are for 2D cases, we can use thesame subdivision strategies for 3D trusses. We extract all trianglesand quads and test if they should be subdivided. In 3D, we use thesame conditions as in 2D (see Fig. 8). ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) Fig. 9. Optimal truss designs in different subdivision levels. The inputfunctional specification is given in Figure 4.
In this section, we introduce the alternating linear programming(ALP) step which optimizes joint positions and cross-section areasof bars for a given topology. ALP serves as the backbone of theproposed approach. Both the local and global topology operationsin the previous section are based on optimization of joint positionsand cross-section areas. Directly optimizing joint positions andaxial forces is a highly nonlinear problem. Therefore, we split theproblem into two linear problems. In Algorithm a, we solve forcedensities alone without changing the joint positions by the groundstructure method. In Algorithm b, we update joint positions andforce densities jointly based on results from Algorithm a.
Let us first recall the basic plastic formulation of the ground structuremethod [Zegard and Paulino 2014], which solves a continuous linearprogramming problem to minimize the total volume of materialunder the premise of force balance with feasible axial forces:minimize a i , s i | E | (cid:213) i = l i a i , (1)subject to B T s = − f , (1a) a i + s i ≥ , i = , . . . , | E | (1b) a i − s i ≥ , i = , . . . , | E | (1c)where each scalar a i is the cross-section area of the i -th bar. B T isthe nodal equilibrium matrix, built from the directional cosines ofthe bars. More details about this matrix are given in the additionalmaterials. s is a vector with the internal (axial) force for all bars,and | E | is the number of bars. f is a vector of the external forcefor all joints. The internal force s i should be within the range ofadmissible axial forces [− σ T a i , σ C a i ] . Here, we assume that themaximal compressive and tensile strains are the same, σ C = σ T = σ ,which is a constant value. We set σ = a i ≥ | s i | . Asthe length of each bar, l i , is positive, l i >
0, the objective functionrequires the cross-section area of each bar, a i , to be its smallestpermissible value, just enough to support the actual axial force ofthat bar. Then, we have a i = | s i | and the following formulation.minimize s i | E | (cid:213) i = l i | s i | , (2)subject to B T s = − f . (2a) he 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 formulationis transformed to:minimize w i | E | (cid:213) i = l i | w i | , (3)subject to C T w = − f . (3a)Here, in Equation 3a, the matrix C is a simpler expression than B because its elements are linear combinations of joint positions. Moredetails about the matrix C are given in the additional materials. In Algorithm a, the joint positions are assumed to be fixed and theaxial forces are the only variables. To further reduce the total volumeof material, we complement it with Algorithm b and calculate thedisplacements of joints to leverage more degrees of freedom. Weassume the initial values of force densities, w , are known by solvingan LP problem in Equation 1 and set the difference of joint positions, u , and the difference of force densities of bars, ∆ w , as variables. Bydirectly rewriting Equation 3, we haveminimize u i , ∆ w i | E | (cid:213) i = sgn ( w i + ∆ w i )( l i + ∆ l i ) ( w i + ∆ w i ) , (4)subject to ( C + ∆ C ) T ( w + ∆ w ) = − f . (4a)Here, we assume that the change of force densities, ∆ w , is small andthat the signs of force densities remain the same, sgn ( w i + ∆ w i ) = sgn ( w i ) . As Algorithm b is applied after Algorithm a, the values, l i , w i , sgn ( w i ) , and C , are all known.To simplify the problem which has a cubic objective functionand quadratic constraints, our goal is to approximate the aboveformulation with a linear programming problem and solve it in asequential manner in conjunction with Algorithm a. By expandingthe objective function, we have ( l i + ∆ l i ) ( w i + ∆ w i ) = ( l i w i + l i ∆ l i w i + ∆ l i w i + l i ∆ w i + l i ∆ l i ∆ w i + ∆ l i ∆ w i ) ≈ ( l i w i + l i ∆ l i w i + l i ∆ w i ) . Here, we remove the higher order terms, and use the factthat l i w i is constant. The objective function is approximated by (cid:205) | E | i = sgn ( w i )( l i ∆ l i w i + l i ∆ w i ) . As shown in Figure 10, ∆ l i ≈ d i ·( u i − u i ) , where d i is the unit direction vector of the i-th barconnecting the joints i i
2. The objective function is linear with u j and ∆ w i as variables, where j = , . . . , | V | and i = , . . . , | E | .The force equilibrium constraint, ( C + ∆ C ) T ( w + ∆ w ) = − f , isequivalent to C T ∆ w + ∆ C T w + ∆ C T ∆ w = because C T w = − f isensured by Algorithm a. Here we remove the small higher orderterm ∆ C T ∆ w . Then the force balance constraint is linearized as C T ∆ w + ∆ C T w = . The matrix C and the vector w are known fromAlgorithm a and the elements in matrix ∆ C are linear combinationsof nodal variations u j = ( u jx , u jy , u jz ) . Then the constraint is alsolinear with respect to u j and ∆ w i .Finally, the formulation for Algorithm b is written as i i i i i i i i i i i i i i i i i iiiiiiiiiiiiiiiii i i i i i i i i i i i i i i i i i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i ) ∆ l i = d i · ( u i − u i )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i u i Fig. 10. (a) The i -th bar connects the joints i and i ; u i and u i arejoint displacements. (b) The length change along the bar direction, ∆ l i ≈ d i · ( u i − u i ) . minimize u i , ∆ w i | E | (cid:213) i = sgn ( w i )( l i w i d i · ( u i − u i ) + l i ∆ w i ) , (5)subject to C T ∆ w + ∆ C T w = , (5a) − δ i ≤ ∆ w i ≤ δ i ; i = , . . . , | E | , (5b) − λ j ≤ u jx , u jy , u jz ≤ λ j ; j = , . . . , | V | . (5c)where λ j and δ i are the bounds of the variables u i and ∆ w i . Ineach iteration, we set small values for these bounds, e.g., δ i = . | w i | and λ i = . l , where ¯ l is the average length of all bars. The above two algorithms are formulated as two LP problems inEquation 1 and Equation 5. The inputs to Algorithm a are the jointpositions, p , and the functional specification such as the externalforces, LOAD , and the supporting points,
SUPP , and the outputs arethe force densities of bars, w . The algorithm a is written as [ w , V ] = ALGa ( p , LOAD , SUPP ) , where V is the total volume of materials.The inputs of Algorithm b are the initial force densities, w , the initialjoint positions, p , and the same functional specification. The outputsare the changing values of joint positions and force densities. Then,Algorithm b is written as [ u , ∆ w ] = ALGb ( p , w , LOAD , SUPP ) . Inthe whole algorithm, we organize them in an alternating way asshown in Algorithm 1. N max is the maximum iteration number and S max is the maximum line search step.In Figure 11, we show the effectiveness of ALP for truss optimiza-tion with three sets of load specifications involving torques. The generality of our formulation allows extensions of our approachfor a broad range of scenarios and applications, such as tacklingmultiple load specifications, respecting stability analysis, and incor-porating project-specific fabrication constraints.
For the input specification with multiple sets of external forces, theALP algorithm is adjusted accordingly. Static equilibrium is requiredfor each set of external loads with generally different internal forces.Thus, assuming trusses are required to withstand K sets of externalforces f , ... , f K , the formulation of Algorithm a is rewritten as: a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) Fig. 11. Three trusses optimized through ALP for different input specifications, each with eight supporting joints and eight external loads along circles. Fromleft to right, (a1), (b1), and (c1) are the initial trusses with a total volume of material consumption being 103.93, 105.34 and 109.58 respectively. (a2), (b2), and(c2) are the optimal trusses using ALP with a total volumes of 69.88, 68.53 and 69.21, respectively.
Algorithm 1
Alternating LP for truss geometry optimization procedure Alternating LP Initial joint positions p ; LOAD and
SUPP ; [ w , V ] = ALGa ( p , LOAD , SUPP ) ; Flag ← True ; N ← while Flag do [ u , ∆ w ] = ALGb ( p , w , LOAD , SUPP ) ; procedure Line search for j = S max do s ← − j ; ˆ p ← p + s u ; [ ˆ w , ˆ V ] = ALGa ( ˆ p , LOAD , SUPP ) ; if ˆ V < V then V ← ˆ V ; p ← ˆ p ; Break ; else if j == S max then Flag = False ; endprocedure N ← N + if N > N max then Flag = False ; endwhile minimize a i , s ki | E | (cid:213) i = l i a i , subject to B T s = − f , a i + s i ≥ , i = , . . . , | E | a i − s i ≥ , i = , . . . , | E | ... B T s K = − f K , a i + s Ki ≥ , i = , . . . , | E | a i − s Ki ≥ , i = , . . . , | E | where k = , . . . , K . Force equilibrium constraints similar to Equa-tion 1 are required for each set of external forces. It is also worthnoting that this set of equations is sufficient to guarantee force equi-librium in response to linear interpolation of the sets of specifiedexternal forces. As each bar needs to support the maximal axialforces from each set of the reaction forces, the cross-section of the i -th bar, a i = max {| s i | , ..., | s i K |} . We denote m i ∈ { , ..., k } as the set index for which the i -th bar attains its maximal axial force, | s m i i | = max {| s i | , ..., | s i K |} . As those indices, m i , could be easilyfound from the result of Algorithm a, corresponding cross sectionsfollow directly, a i = | s m i i | . Similarly, by defining w m i i = s m i i / l i , theformulation of Algorithm b for the multiple-load case can be writtenas:minimize u i , ∆ w i k | E | (cid:213) i = sgn ( w m i i )( l i w m i i d i · ( u i − u i ) + l i ∆ w m i i ) , subject to C T ∆ w + ∆ C T w = ,... C T ∆ w K + ∆ C T w K = , − δ i ≤ ∆ w ki ≤ δ i , − λ j ≤ u jx , u jy , u jz ≤ λ j , where i = , . . . , | E | , k = , . . . , K , and j = , . . . , | V | . Here, onlythe nodal variations, u i , and the change of force densities, ∆ w ki ,are variables, others are known from Algorithm a. Using the samealternating scheme in Section 5.3 by adjusting the Algorithm a andb accordingly, our method can tackle cases of multiple-load inputspecifications. See Figure 12 for an example. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) Fig. 12. A truss designed to support three sets of different external loads.Here, red joints are fixed and blue joints have external loads applied. Withthree sets of external loads shown in (a)-(c), our method creates an optimaltruss that supports all of them, as shown in (d).
Truss stability, in particular, the external stability, has been investi-gated extensively for its importance in keeping its rigid shape. Bycomparing the number of degrees of freedom against constraintcounts, a condition for external stability is given in [Kassimali 2011]as | E | + r ≥ d | V | . ith | E | being the number of bars and r being the count of reactioncomponents, the sum on the left-hand-side includes the number ofboth internal (e.g., fixed edge lengths) and external (e.g., supportedvertices) constraints. The right-hand-side counts the total numberof degrees of freedom, as a product of the number of joints, | V | ,and the number of degrees of freedom for each joint, d , which is2 for plane trusses and 3 for spatial trusses. A truss is externallystable if the condition is met and unstable otherwise. There aremultiple approaches to integrate stability analysis to our framework,e.g., through a post-processing step checking the aforementionedcondition and adding necessary auxiliary bars with a minimumcross-section to increase the number of | E | on the left side and thusachieving external stability. Figure 13 provides an example. Similarlyto stability analysis, buckling analysis through simulation packageslike ANSYS are applicable as a feedback procedure to improve thetruss structural properties of interest in an interactive manner. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) Fig. 13. (a) With | E | =26, r=4, | V | =15, the truss is externally stable. (b1) With | E | =27, r=4, | V | =16, the truss is is externally unstable. (b2) With one morebar added, | E | =28, r=4, | V | =16, the truss in (b1) becomes externally stable. Besides general engineering settings, our system is able to incorpo-rate project-specific fabrication constraints through minor adjust-ments or complementary post-processing procedures. For projectsrequiring a limited number of parts, bounding the number of sub-division steps controls the amount of bars and joints. Meanwhile,while our algorithm reduces the total volume of material consump-tion based on a continuous selection of cross-sections in the defaultsetting, for engineering practices which demand viable bars to bechosen from a set of predefined category or a discrete number oftypes, our method can be combined with other published algorithms.We show an example in limiting the types of cross sections by ap-plying [Jiang et al. 2017] directly as a post-processing procedure inFigure 14. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) Fig. 14. A truss optimized by our method in Figure 24, with a total volumeof 3.188, is post-processed by a discrete optimization algorithm in [Jianget al. 2017] to obtain trusses constructed by (a) three types of cross-sectionareas with a total volume of 4.674, and (b) five types of cross-section areaswith a total volume of 3.832.
In this section, we illustrate truss designs using our framework fordifferent types of input specifications and compare the results withstate-of-the-art methods on selected benchmark design problems.
We show the results of our method for different types of functionalspecifications. We present 2D truss designs with a parallel equi-librium force system in Figure 22, a concurrent equilibrium forcesystem in Figure 23, and a non-concurrent equilibrium force systemin Figure 24. We illustrate examples of designs for the same inputexternal forces and supporting joints but different design regions inFigure 25. For 3D trusses, we demonstrate the results for input of aparallel equilibrium force system in Figure 27, a concurrent equilib-rium force system in Figure 28, and a non-concurrent equilibriumforce system in Figure 29. In addition, we show examples based onreal functional requirements such as a 2D bike frame in Figure 26, a3D cantilever in Figure 30, and a 3D bridge in Figure 1.
Our framework is implemented in Matlab R2016 on a workstationwith an Intel Xeon X5550 2.67 GHz processor. We use Mosek [ApS2015] as the solver for linear programming. For the ALP algorithm,we set the maximum iteration number N max =500 and the maxi-mum line search step S max =10. For the number of grid points inthe initialization, n , we use the default value, the number of jointsspecified in the functional specification. We set the number of maxiteration of Phase 1, P max1 =
5. The number of subdivisions in Phase2, P max2 , controls the trade off between number of bars and trussweight and is set depending on the context. We set the thresholds ϵ and ϵ in Section 4.3 as ϵ = .
002 ¯ a , ϵ = .
01 ¯ d , where ¯ a is theaverage cross-section, and ¯ d is the average value of distances be-tween joints specified in the functional specification. In Table 2, foreach optimized truss, we report parameters of trusses such as thenumber of bars, the total volume of material, and the computationtime of coarse truss optimization and different levels of structurerefinement. To test the sensitivity of our framework to the initialization, weshow optimized trusses starting from initializations with differ-ent numbers of intermediate joints and bars (before subdivisionoperations are conducted). Figure 15 shows that adding differentintermediate joints and bars results in different trusses. However,these trusses have similar structure and they all function as robustdiscrete approximations of the optimal truss for further subdivision.
We compare the pro-posed 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 theimplementation of SQP and GDM are given in the appendix. Jianget al.’s framework contains multiple parts and we only compare tothe solver that is comparable to our proposed ALP algorithm (Theobjective function in Jiang et al. is different though). For variousinput specifications, we compare the quality (volume) of the out-put. Fig. 16 and Table 1 show the results. The results demonstratethat our proposed optimization converges to a better solution thancompeting approaches, especially GD struggles to find meaningful ig. 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 Table 1. Comparison of the achieved volume and running time for different numerical optimization methods. We only compare the result of the geometryoptimization part without performing topology operations or subdivision. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) Fig. 15. Comparing different initializations. From left to right we initializewith 4 ×
4, 6 ×
6, 8 ×
8, and 10 ×
10 intermediate joints and 124, 364, 732, and1228 intermediate bars. solutions. We noticed that competing algorithms can be improvedwhen using the ground structure method (GSM) for initialization.In the table, we use a * to indicate the version of SQP and GD thathas been initialized using GSM. In our comparison, our method gen-erates the best results for all input specifications. SQP* can matchour result in some cases. We can also observe that our method isfast and better scales to larger inputs than SQP. The reason for thepoor performance of Jiang et al. [Jiang et al. 2017] is that they usesoft constraints instead of hard constraints in the formulation.
We comparethe performance of our method with several previous methods usingthe functional specifications, solutions, and running times providedin their papers. We denote the methods in [Descamps and Coelho2013], [Gilbert and Tyas 2003], [He and Gilbert 2015], and [Sokół2017] as D2013, G2003, H2015 and S2017, respectively. In Figure 1and Figure 17, we compare our method with [Descamps and Coelho2013] for a 2D and a 3D bridge design problem. In Figure 18 and 19,we compare our method with [He and Gilbert 2015] and [Gilbertand 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 shownin Table 3. We can observe that our method is orders of magnitudefaster, even though we can achieve a lower volume than previouswork. Unfortunately, these methods do not have code or executablespublicly available, so we cannot test them on the same machine. However, it seems unlikely that this significant difference in runningtime can be overcome by slightly faster hardware.
For the cantileverdesign problems in Figure 18 and 19, the total volumes of analyticalsolutions are 4.498115 and 4.232168. The volumes of our discretedesigns for these two cases are 4.498635 and 4.3223, which are closerto the analytical solutions than the previous work, [Gilbert and Tyas2003] and [He and Gilbert 2015]. For the three force problem, wecompare our result with the analytical solution presented in [Sokółand Lewiński 2010]. In Figure 21, we show the computed discretetruss (130 bars, volume: 6.838) on the right which is visually similarto the analytical solution on the left (volume: 6.831).
Discussion and Limitations.
Compared with previous work, thegeometry optimization in terms of both the axial forces and jointpositions through two linear programming problems alternativelyapplied in ALP provides more degrees of freedom compared with theoriginal formulation of the ground structure method which solvesa single linear programming problem. Splitting a highly nonlinearprogramming problem into two linear ones also attains better effi-ciency compared with the original nonlinear formulation. Moreover,the two categories of topological operations, local and global, com-plementarily allow both flexible yet stable topology changes andmanipulation. Most importantly, the subdivision approach, whichhas been overlooked by previous work, is a natural choice for topol-ogy refinement from coarse to fine, which creates valid topologies atdifferent levels both efficiently and robustly. Despite the efficiencyand efficacy, our optimization framework cannot guarantee a globaloptimum. However, in simple special cases where the analyticaloptimum is known, we observe that the method almost reachesthe known global optimum. As our method is based on subdivisionof edge-networks on surfaces, the optimal trusses in 3D also needto constitute sheets of surfaces. For general 3D specifications, thesubdivision approach requires an initial surface-like structure, ora structure that consists of multiple sheets, which is a challengingproblem for future work. The proposed ALP algorithm will workfor general 3D structures, however. a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) ( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e ) ( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f ) ( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) ( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e ) ( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f )( f ) ( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д )( д ) Fig. 16. Volume optimization comparison with different methods. From left to right: (a) the ground structure method (GSM), (b) ALP (ours), (c) SQP withinitialization of cross-section of constant values, (d) SQP with initialization of cross-section from (a), (e) the gradient descent method (GDM) with initializationof cross-section of constant values, (f) uses the gradient descent method with initialization of cross-section from (a). (g) the method used in [Jiang et al. 2017].Table 1 shows our ALP method has better performance.
Fig. Initial Truss Optimal Coarse Truss Truss Subdivision 1 Truss Subdivision 2 Truss Subdivision 3Bars VolumeGSM 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.9723 404 3.125 22 2.967 5.86 46 2.913 5.85 118 2.909 8.86 358 2.907 15.1624 124 3.872 16 3.365 0.64 36 3.248 0.88 100 3.203 5.46 324 3.188 17.2825(a1) 835 3.375 19 3.210 5.02 45 3.183 1.43 133 3.170 8.79 453 3.163 29.5625(b1) 835 3.335 25 3.114 5.42 67 3.041 1.81 211 3.021 7.77 739 3.016 17.8426 466 4.340 25 4.294 4.92 67 4.289 3.22 211 4.288 9.68 739 4.287 35.6327 2061 19.510 9 19.081 4.28 18 18.700 0.46 37 18.610 2.39 76 18.587 5.5628 1118 13.891 24 11.742 5.87 60 11.724 7.73 180 11.718 19.09 612 11.713 57.4129 3674 15.429 10 14.628 1.61 16 14.416 1.53 28 14.325 1.97 52 14.316 2.1130 3674 30.850 24 29.049 8.70 58 28.507 1.96 174 28.285 14.75 598 28.217 64.0131 3674 20.545 10 18.900 3.41 16 18.542 0.4 30 18.412 6.22 59 18.375 14.10
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, aswell as the computational time for different stages. Please note that the results from optimal coarse trusses are already better than the results from theground structure method without introducing additional joints and bars. Topological refinement through subdivision provides additional reduction in materialconsumption generally in less than one minute in our current implementation which still has a significant room for further improvement.
Fig. Method Bars(initial) Bars(final) Time Volume1 D2013 258 96 1376s 408.807Ours 258 201
17 D2013 31 19 n/a 34.977Ours 804 25
18 G2003 >1 billion n/a >6h 4.4998Ours 105 2178
30s 4.4986
19 H2015 12,456,601 4244 4875s 4.3228Ours 105 2178
30s 4.3223
20 S2017 >7 billion 40 >1h n/aOurs 1118 24
Table 3. Comparisons with previous work. Compared with previous ap-proaches, our framework consistently creates truss designs with smallervolumes with significantly shorter computational times.
We present a method for the design of optimal trusses satisfyingfunctional specifications with minimized material consumption. Thecore components of the proposed approach include an alternating ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) ( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e ) Fig. 17. A comparison with the method in [Descamps and Coelho 2013]on a 2D bridge model. (a) The input functional specification and initialtruss for the method in [Descamps and Coelho 2013]. (b) Optimal truss in[Descamps and Coelho 2013] (number of bars: 19, volume: 34.977). (c) Ouroptimal coarse truss (number of bars: 25, volume: 34.593). Refined structuresafter 1 and 2 rounds of subdivisions are shown in (d) and (e). ig. 18. A comparison with the method in [Gilbert and Tyas 2003] on abenchmark design problem (Hemp cantilever). Compared with the bestresult presented in [Gilbert and Tyas 2003], obtained with 6h50m of com-putation, shown in (a), which used an initial truss of 116,288,875 bars andobtain a total volume of 4.499827, we obtain a result, shown in (b), with 2178bars in 30s, which achieves a total volume of 4.498635. For this problem, ananalytical solution with an infinite amount of bars exists, which has a totalvolume of 4.498115.Fig. 19. A comparison with the method in [He and Gilbert 2015] on abenchmark design problem (Hemp cantilever). (a) The optimal truss of [Heand Gilbert 2015] (number of bars: 4244, running time: 4875s, volume: 4.3228).(b) Our optimized truss (number of bars: 2178, running time: 30s, volume:4.3223). The analytical solution of optimal volume is 4.3217. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) Fig. 20. A comparison with the method in [Sokół 2017] on a 3D model. (a)is the optimal structure in [Sokół 2017]. To get this topology, they used a 3Dgrid of 50X50X50 joints and 7,318,049,198 bars, the running time is close to 2hours. (b) and (c) are results of our method. To get the initial structure in (b),we use 1118 bars and running time is less than 10s. To get the optimizationresult in (c) our running time is less than 20s. linear programming formulation for geometry optimization and twosets of topological operations. The subdivision scheme inspired byMichell’s theoretical studies utilized in the global topology refine-ment step plays a crucial role for the efficiency and efficacy of theproposed approach. The performance of our framework is validatedby comparisons with multiple previous studies in different scenar-ios, which indicate that our method creates trusses with smallervolumes and is over two orders of magnitude faster in terms ofcomputational speed. For future work, it would be exciting to studydynamic structures created by trusses with movable parts that have ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) Fig. 21. A comparison with an analytic solution in [Sokół and Lewiński2010]. The volume of the analytic solution with an infinite amount of barsis 6.831 and the volume of our optimized discrete truss (130 bars) is 6.838. multiple configurations, e.g., robotic and mechatronic systems withtrusses. Moreover, a better theoretical understanding of optimaltrusses in 3D would be inspiring for the entire field.
REFERENCES
MOSEK ApS. 2015.
The MOSEK optimization toolbox for MATLAB manual. Version 7.1(Revision 28). http://docs.mosek.com/7.1/toolbox/index.htmlMoritz Bächer, Emily Whiting, Bernd Bickel, and Olga Sorkine-Hornung. 2014. Spin-it:optimizing moment of inertia for spinnable objects.
ACM Transactions on Graphics(TOG)
33, 4 (2014), 96.Phillipe Block and John Ochsendorf. 2007. Thrust Network Analysis: A New Methodol-ogy For Three-Dimensional Equilibrium.
J. Int. Assoc. Shell and Spatial Structures
48, 3 (2007), 167–173.HSY Chan. 1966. Minimum weight cantilever frames with specified reactions. Universityof Oxford, Department of Engineering Science, Eng. Laboratory, Parks Road.Fernando de Goes, Pierre Alliez, Houman Owhadi, and Mathieu Desbrun. 2013. On theEquilibrium of Simplicial Masonry Structures.
ACM Trans. Graph.
32 (2013).Benoît Descamps and Rajan Filomeno Coelho. 2013. A lower-bound formulation forthe geometry and topology optimization of truss structures under multiple loading.
Structural and multidisciplinary optimization
48, 1 (2013), 49–58.William S Dorn. 1964. Automatic design of optimal structures.
Journal de mecanique
Engineering Computations
20, 8 (2003), 1044–1064.L He and M Gilbert. 2015. Rationalization of trusses generated via layout optimization.
Structural and Multidisciplinary Optimization
52, 4 (2015), 677–694.WS Hemp. 1973. Michell’s structural continua.
Optimum Structures (1973), 70–101.Caigui Jiang, Chengcheng Tang, Hans-Peter Seidel, and Peter Wonka. 2017. Designand Volume Optimization of Space Structures.
ACM Transactions on Graphics
36, 4(2017), 1–14.Caigui Jiang, Chengcheng Tang, Marko Tomičí, Johannes Wallner, and HelmutPottmann. 2015. Interactive Modeling of Architectural Freeform Structures: Combin-ing Geometry with Fabrication and Statics. In
Advances in Architectural Geometry2014 . Springer, 95–108.Aslam Kassimali. 2011.
Structural analysis . Cengage Learning.Martin Kilian, Davide Pellis, Johannes Wallner, and Helmut Pottmann. 2017. Material-minimizing forms and structures.
ACM Trans. Graphics
36, 6 (2017), article 173.https://doi.org/10.1145/3130800.3130827 Proc. SIGGRAPH Asia.T Lewiński and GIN Rozvany. 2007. Exact analytical solutions for some popularbenchmark problems in topology optimization II: three-sided polygonal supports.
Structural and Multidisciplinary Optimization
33, 4-5 (2007), 337–349.T Lewiński and GIN Rozvany. 2008a. Analytical benchmarks for topological optimiza-tion IV: square-shaped line support.
Structural and Multidisciplinary Optimization
36, 2 (2008), 143–158.T Lewiński and GIN Rozvany. 2008b. Exact analytical solutions for some popularbenchmark problems in topology optimization III: L-shaped domains.
Structuraland Multidisciplinary Optimization
35, 2 (2008), 165–174.Yang Liu, Hao Pan, John Snyder, Wenping Wang, and Baining Guo. 2013. ComputingSelf-Supporting Surfaces By Regular Triangulation.
ACM Trans. Graph.
32 (2013).Jonàs Martínez, Jérémie Dumas, Sylvain Lefebvre, and Li-Yi Wei. 2015. Structureand appearance optimization for controllable shape design.
ACM Transactions onGraphics (TOG)
34, 6 (2015), 229.P Martinez, P Marti, and OM Querin. 2007. Growth method for size, topology, and geom-etry optimization of truss structures.
Structural and Multidisciplinary Optimization
33, 1 (2007), 13–26.J Clerk Maxwell. 1870. I.-On Reciprocal Figures, Frames, and Diagrams of Forces.
Earthand Environmental Science Transactions of the Royal Society of Edinburgh
26, 1 (1870), a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) Fig. 22. A truss design with the input being a 2D parallel force system. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) Fig. 23. A truss design with the input being a 2D concurrent force system. Note that the concave quadrilateral in (a) is not subdivided. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) Fig. 24. A truss design with the input being a 2D non-concurrent force system. Note that the orthogonality of the two families of trusses, in compression andtension, is improved with subdivision. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) Fig. 25. Truss designs for different design regions. Top: truss designs for a functional specification that the design region is the entire plane. Bottom: trussdesigns for the same functional specification except that the design region is the upper half plane, with the lower half plane being an obstacle to be avoided.
The London, Edinburgh, and Dublin Philosophical Magazine and Journalof Science
ACM Trans. Graph.
32 (2013).Nico Pietroni, Davide Tonelli, Enrico Puppo, Maurizio Froli, Roberto Scopigno, andPaolo Cignoni. 2015. Statics aware grid shells. In
Computer Graphics Forum , Vol. 34.Wiley Online Library, 627–641.William Prager and George IN Rozvany. 1977. Optimization of structural geometry. In
Dynamical systems . Elsevier, 265–293.Romain Prévost, Emily Whiting, Sylvain Lefebvre, and Olga Sorkine-Hornung. 2013.Make it stand: balancing shapes for 3D fabrication.
ACM Transactions on Graphics(TOG)
32, 4 (2013), 81.GIN Rozvany. 1998. Exact analytical solutions for some popular benchmark problemsin topology optimization.
Structural and Multidisciplinary Optimization
15, 1 (1998), 42–48.Jeffrey Smith, Jessica Hodgins, Irving Oppenheim, and Andrew Witkin. 2002. Creatingmodels of truss structures with optimization.
ACM Trans. Graph.
21 (2002).Tomasz Sokół. 2017. On the Numerical Approximation of Michell Trusses and theImproved Ground Structure Method. In
World Congress of Structural and Multidisci-plinary Optimisation . Springer, 1411–1417.Tomasz Sokół and Tomasz Lewiński. 2010. On the solution of the three forces problemand its application in optimal designing of a class of symmetric plane frameworksof least weight.
Structural and Multidisciplinary Optimization
42, 6 (2010), 835–853.Chengcheng Tang, Xiang Sun, Alexandra Gomes, Johannes Wallner, and HelmutPottmann. 2014. Form-finding with Polyhedral Meshes Made Simple.
ACM Tran.Graph.
33 (2014).Nobuyuki Umetani, Takeo Igarashi, and Niloy J. Mitra. 2012. Guided Explorationof Physically Valid Shapes for Furniture Design.
ACM Transactions on Graphics(Proceedings of SIGGRAPH 2012)
31, 4 (2012). a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) ( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e ) Fig. 26. A bike frame design. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) ( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e ) Fig. 27. A truss design with the input being a 3D parallel force system. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) Fig. 28. A truss design with the input being a 3D concurrent force system. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) Fig. 29. A truss design with the input being a 3D non-concurrent force system.
Etienne Vouga, Mathias Höbinger, Johannes Wallner, and Helmut Pottmann. 2012.Design of self-supporting surfaces.
ACM Trans. Graph.
31, 4 (2012).Weiming Wang, Tuanfeng Y Wang, Zhouwang Yang, Ligang Liu, Xin Tong, WeihuaTong, Jiansong Deng, Falai Chen, and Xiuping Liu. 2013. Cost-effective printing of3D objects with skin-frame structures.
ACM Trans. Graph.
32, 6 (2013).Jiaxian Yao, Danny M. Kaufman, Yotam Gingold, and Maneesh Agrawala. 2017. Interac-tive Design and Stability Analysis of Decorative Joinery for Furniture.
ACM Trans.Graph.
36, 2, Article 20 (March 2017), 16 pages. https://doi.org/10.1145/3054740Tomás Zegard and Glaucio H Paulino. 2014. GRAND-Ground structure based topologyoptimization for arbitrary 2D domains using MATLAB.
Structural and Multidisci-plinary Optimization
50 (2014).Qingnan Zhou, Julian Panetta, and Denis Zorin. 2013. Worst-case structural analysis.
ACM Trans. Graph.
32, 4, Article 137 (July 2013), 12 pages.
APPENDIX8.1 Equilibrium Force Systems
The input functional specification is essentially an equilibriumforce system. For a system with forces F , F ,..., F n acting at points P , P ,..., P n , there are two constraints of equilibrium: the balance offorces, (cid:205) ni = F i = (cid:174)
0, and the balance of torques, (cid:205) ni = F i × P i = (cid:174)
0. Inmechanical engineering, the equilibrium force systems are usuallyclassified into three types: parallel force systems, concurrent forcesystems, and non-concurrent force systems as shown in Figure 32. a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) ( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e ) Fig. 30. A 3D cantilever design. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) ( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d )( d ) ( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e )( e ) Fig. 31. A 3D example.
First, lets recall the formulation in Eq.3. In this formulation, jointpositions are considered as known values. The objective functionand constraints in the above formulation are both linear. If the jointpositions u j , j = , ..., | V | are considered as design variables, the for-mulation is a nonlinear programming with nonlinear constraints. Tosolve such a problem by SQP, we reformulate the problem similarlyas in [Descamps and Coelho 2013]. We introduce slack variables, w + i and w − i , which represent the magnitudes of compressive andtensile forces respectively: w + i = max ( , w i ) , w − i = max ( , − w i ) .Then, we have w i = w + i − w − i , and | w i | = w + i + w − i . Here we denotethe vectors consisted of all w + i and w + i as w + and w − . We assumethe i th bar is connected by joints i i
2, then l i = ( u i − u i ) .The elements in matrix C are linear combinations of joint positions.The whole problem is reformulated asminimize u j , w + i , w − i | E | (cid:213) i = ( u i − u i ) ( w + i + w − i ) , (6)subject to C T ( w + − w − ) = − f , (6a) w + i ≥ , (6b) w − i ≥ , (6c)The algorithm in [Descamps and Coelho 2013] used the aboveformulation by SQP method. We implement the algorithm using the f mincon routine in Matlab. To use the GDM, we first transfer the constrained optimizationproblem in Eq. 6 into an unconstrained formulation by adding anextra penalty term into the objective function as f = λ E volume + λ E static , (7) where E volume = | E | (cid:213) i = ( u i − u i ) ( w + i + w − i ) , and E static = ( C T ( w + − w − ) + f ) . This objective function is defined similarly in [Jiang et al. 2017]. Theconstraint of w + i ≥ w − i ≥ w + i and w − i with ( ϱ + i ) and ( ϱ + i ) , where ϱ + i and ϱ − i are new variables withoutany constraints. Usually, the weight of penalty term λ should bemuch larger than λ . The above object function is differentiable. Wedo the comparison using standard gradient decent method. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) Fig. 32. Three equilibrium force systems. (a) a parallel force system. (b) aconcurrent force system. (c) a non-concurrent force system. DDITIONAL MATERIALS8.4 Special Cases
The theorem of Maxwell(1872) states that (cid:205) f i l i = C , where C isthe same constant for any statically admissible truss layout consid-ering a given set of external forces. As a consequence, if it is foundpossible to design a structure all of whose members are in tension,or alternatively compression, then the optimum design has beenachieved, because (cid:205) | f i | l i = | C | . The above statement assumes theinput force system is precisely given. For example, a self-equilibriumforce system with no supporting points, such as the 3-force problemin 2D and 4-force problem in 3D. However, it is not correct whensome supporting points are given without prescribing their reactionforces. For example, as shown in Figure 33, we can find a structureof two bars with compression on the left, but it is not an optimalstructure because reaction forces at the boundary could be changedto achieve a better results as show on the right. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) Fig. 33. The inputs are two supporting points and one external load. (a) atruss of two bars in compression. (b) a better design as reaction forces atsupporting points are changed to obtain a less material structure.
According to the theorem of Maxwell, we can easily find somespecial cases. For example, given a concurrent force systems withforces F , F ,..., F n through a common point O and their correspond-ing action points P , P ,..., P n , we have F i = λ i −−→ OP i . If all the λ i havethe same sign, either positive or negative, then the truss connect all P i and O is an optimum design. Note that the optimum for thesespecial boundary conditions is usually not unique. Figure 34 showsmultiple optimum designs for the same boundary condition where3 forces are self-equilibrium in 2D plane. From a design standpoint,the non-uniqueness provide a flexibility of topology and geometrydesign for these special boundary conditions. ( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a )( a ) ( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b )( b ) ( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c )( c ) Fig. 34. Three different truss designs for the same special boundary condi-tion. These truss use the same amount of material. There are infinite numberof such optimal trusses for this special boundary condition.
To systematically investigate the boundary condition under whicha design with all bars in compression or tension can be found isinteresting but out of the scope of this paper. Our numerical compu-tational method can be use to detect such kind of boundary condi-tion in the initial truss generation stage. Under such kind of specialboundary conditions, further subdivision is not necessary.
The nodal equilibrium matrix, B , transforms the magnitudes ofinternal forces, s = ( s , s , . . . , s | E | ), at each bar, into forces alongthe x , y and z axes at each joint. If a truss is in equilibrium, then eachjoint must be in equilibrium. The equilibrium equation is written as B T s = − f . (8)Let’s take a look at the equilibrium equations of a particular i -thjoint with a valence of n i , connected to nodes p i , j , j = , . . . , n i through bars b i , j , j = , . . . , n i respectively. Then the equilibriumequations at the joint i are n i (cid:213) j = x p i , j − x i l b i , j s b i , j = − f ix , n i (cid:213) j = y p i , j − y i l b i , j s b i , j = − f iy , n i (cid:213) j = z p i , j − z i l b i , j s b i , j = − f iz , where l b i , j = (cid:113) ( x p i , j − x i ) + ( y p i , j − y i ) + ( z p i , j − z i ) is the lengthof 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 linearcombinations of the directional cosines of the bars.If we use force densities w i = s i / l i as variables, the above equa-tions are transformed to simpler forms: n i (cid:213) j = ( x p i , j − x i ) w b i , j = − f ix , n i (cid:213) j = ( y p i , j − y i ) w b i , j = − f iy , n i (cid:213) j = ( z p i , j − z i ) w b i , j = − f iz . The above equilibrium equations could be written in a similarmatrix form as C T w = − f . (9)All the elements in the matrix C are linear combinations of the jointcoordinates. Similarly, the elements in the matrix ∆ C in Equation5a are linear combinations of nodal variations.in Equation5a are linear combinations of nodal variations.