A Python script for adaptive layout optimization of trusses

Numerical layout optimization employing an adaptive ‘member adding’ solution scheme provides a computationally efficient means of generating (near-)optimum trusses for problems involving single or multiple load cases. To encourage usage of the method, a Python script is presented, allowing medium to large-scale problems to be solved efficiently. As well as handling multiple load cases, the short (98 line) script presented can tackle truss optimization problems involving unequal limiting tensile and compressive stresses, joint costs, and non-convex polygonal domains, with or without holes. Various numerical examples are used to demonstrate the efficacy of the script presented.


Introduction
Truss layout optimization using the 'ground structure' approach provides a fully automated means of identifying (near-)optimum truss structures. Although first proposed in the 1960s by Dorn et al. (1964), it has not been widely used in practice for various reasons. For example, the solutions generated can appear impractical, at least when conventional manufacturing techniques are envisaged. However, in recent decades new manufacturing techniques have been developing apace, e.g. additive manufacturing, and means of rationalizing the solutions obtained via layout optimization have also been developed (e.g. He and Gilbert 2015). Also, in 2003, an adaptive 'member adding' scheme was proposed by Gilbert and Tyas (2003), which allowed solutions to be obtained more quickly, and also for extremely large problems (e.g. >1,000,000,000 members) to be solved. This means that layout optimization can be used to obtain highly accurate benchmark solutions that can be used to validate analytical solutions (e.g. Sokół 2011b; Sokół and Rozvany 2012), and also to discover new optimum structural forms (e.g. Darwich et al. 2010;Tyas et al. 2011;Fairclough et al. 2018). To increase the practical applicability of solutions, Pritchard et al. (2005) extended the adaptive 'member adding' scheme to solve problems involving multiple load cases whilst Tyas et al. (2006) incorporated structural stability considerations. Also, Smith et al. (2016) utilized layout optimization to design truss-like metallic components suitable for fabrication via additive manufacturing, with the techniques involved recently extended by He et al. (2018). In parallel, the layout optimization technique has been successfully transferred to other applications. For example, Bolbotowski et al. (2018) utilized layout optimization to design (near-)optimum grillages and Smith and Gilbert (2007) and Gilbert et al. (2014) proposed the 'discontinuity layout optimization' (DLO) method to identify critical patterns of discontinuities at collapse for in-plane and out-of-plane plasticity problems.
However, although interest in numerical layout optimization has grown in recent years, the power of the method still appears under-appreciated in the structural optimization research community. One possible reason for this is that accessible educational resources for use by new researchers have been lacking. Although a script was made available by Sokół (2011a), the Mathematica language employed is not commonly used in the community. Matlab scripts developed by Zegard and Paulino (2014) and Zegard and Paulino (2015) helped in this respect, although these contain relatively complex functions that could potentially reduce uptake. Also, most significantly, although proposed in 2003, the 'member adding' scheme has not yet been incorporated in any of the aforementioned publicly available scripts, potentially leading researchers to underestimate the potential computational efficiency of the layout optimization method. In contrast, the 99-line Matlab script for SIMP (Sigmund 2001) helped pave the way for worldwide research in the field of continuum topology optimization. As an indication of the dominance of SIMP over numerical layout optimization in terms of published educational source codes, only one of the 22 scripts listed by Wei et al. (2018) employ the latter.
In order to provide an accessible yet computationally efficient educational implementation of numerical layout optimization, this paper introduces a simple Python script. With 98 lines, it incorporates the efficient adaptive 'member adding' scheme for 2D problems subject to multiple load cases, with the potential for unequal limiting tensile and compressive stresses, joint costs, and non-convex polygonal domains (with or without holes). The paper is organized as follows: first, the mathematical layout optimization formulation is introduced; second, important code sections are explained; third, numerical examples are shown to demonstrate the efficacy of the script; finally, conclusions are drawn.

Single-load case problem
The standard layout optimization process (Dorn et al. 1964;Gilbert and Tyas 2003;Pritchard et al. 2005) involves a number of steps, as shown in Fig. 1. Firstly the design domain, load and support conditions are specified, Fig. 1a; secondly, the design domain is discretized using nodes, Fig. 1 Steps in layout optimization: a specify design domain, loads and supports; b discretize domain using nodes; c interconnect nodes with potential truss members; d use optimization to identify the optimal structural layout (a) (b) (c) (d) Fig. 1b; thirdly, these nodes are interconnected with potential members to create a 'ground structure', Fig. 1c; finally, an optimum layout is identified by solving the optimization problem below (Fig. 1d): where, V is the structural volume, a = [a 1 , a 2 , ..., a m ] T is a vector containing member cross-sectional areas, with m denoting the number of members. l = [l 1 , l 2 , ..., l m ] T is a vector of member lengths. B is a 2n × m equilibrium matrix comprising direction cosines, with n denoting the number of nodes and q is vector containing the internal member forces; a simple example of the equilibrium constraint is presented in Fig. 2. f is a vector containing the external forces; σ + and σ − are limiting tensile and compressive stresses respectively.

Multiple load case problem
Problem (1) can easily be extended to obtain solutions for multiple load case plastic design problems. In this case the equilibrium constraint (1b), stress limits (1c) and (1d) must be satisfied in each load case, e.g. the former gives: Bq k = f k , for k = 1, 2, ..., p, where, k is the load case identifier and p represents the total number of load cases.

Joint costs
It was found in previous studies (e.g. Gilbert and Tyas 2003;Pritchard et al. 2005) that layout optimization will often generate structures that are complex in form, containing large numbers of short members. A simple means of addressing this is to include a notional joint cost (Parkes 1978), which has the effect of penalizing short members and hence simplifying the solution in many cases:l = [l 1 + s, l 2 +s, ..., l m +s] T in the objective function (1a), wherel is the augmented member length vector with predefined joint cost (or length) s.

Full problem
Considering multiple load cases and joint costs, the full problem can now be written as: which is, like (1), a linear programming (LP) problem that can be solved efficiently using modern optimization solvers. However, when using the ground structure approach the number of potential numbers grows rapidly with the number of nodes employed, with up to n(n−1) 2 members required. This can lead to extremely large optimization problems that are computationally expensive to solve. To address this, a 'member adding' scheme, which is based on the 'column generation' approach (Dantzig and Wolfe 1960) was proposed by Gilbert and Tyas (2003). This allows problem (3) to be decomposed into a number of smaller subproblems, with information from the dual problem used to guide the process.

Adaptive 'member adding' scheme
Using the adaptive 'member adding' scheme (Gilbert and Tyas 2003;Pritchard et al. 2005), initially only a reduced set of members are used to solve problem (3), and remaining potential members in the 'ground structure' are gradually added to the set, until the following constraint is satisfied for all potential members (readers who are interested in derivation of Eq. 4 can refer to Appendix A): where i is the summed maximum virtual strain of member i in all load cases; u k i is a vector containing virtual displacement of nodes connected by member i, under load case k. Note that, when a primal-dual interior point method is used to solve problem (3), u k i is obtained automatically with no extra cost. Although only a reduced set of members are used in Eq. 3, a full virtual displacement field u k of all nodes can still be generated. If any violation of Eq. 4 is found in the virtual displacement field, the most violated potential members can be added to the reduced set of members. Then (3) can be solved again, this time generating an improved virtual displacement field in which (4) is no longer violated for the newly added members. This adaptive 'member adding' process continues until no violation is detected; according to duality theory, the solution obtained will then be equivalent to that obtained by solving the full problem from the outset. The entire process is shown in Fig. 3.

Python implementation
Python is an open-source interpreted high-level programming language that is becoming increasingly popular for solving scientific and engineering problems. Unlike Matlab, where built-in toolboxes are automatically loaded when started, in Python such tools are imported using the keyword at the beginning of the script. These extra lines allow dependencies on external tools to be clearly seen.
In this paper, a number of freely available tools are used, including mathematical tools such as scipy and numpy, the convex optimization tool cvxpy (Diamond and Boyd 2016;Agrawal et al. 2018), and the geometry tool shapely.
In this section the various functions used to construct and solve the layout optimization problem (3) are first Fig. 3 Flowchart of the adaptive 'member adding' process introduced; then the main workflow is outlined. As Python is, like Matlab, an interpreted language, loop structures are generally inefficient. Therefore, to improve execution speed, vectorized operations are utilized in place of loop structures wherever possible.

Equilibrium matrix B
For a member i, interconnecting nodes I and I I , its contribution to equilibrium matrix B can be written as: where, X i = x I I i − x I i and Y i = y I I i − y I i are the projected length of the member of length l i in the x and y axis directions respectively, and where [x I i , y I i ] and [x I I i , y I I i ] are the co-ordinates of nodes I and I I respectively. If a connected node is supported, then the coefficients of the corresponding row are set to zero (i.e. supported degrees of freedom are removed). Thus, taking into account boundary conditions, the above formulation can be expressed as: where d 0 , d 1 , d 2 , d 3 are 0-1 coefficients describing the degree of freedom in the x and y directions for the connected nodes The equilibrium matrix B is then assembled using the B i matrices for all members. B is a sparse matrix, which can be expressed using a row-column-value format comprising only non-zero entries, constructed by concatenating the coefficients calculated in Eq. 6 as well as their corresponding row (degree of freedom at nodes) and column (member) numbers.

Solving the LP problem
In the script, problem (3) is solved in function solveLP.
In Python, the use of object oriented programming (OOP) allows problem (3) to be written in a format that is very similar to its natural expression, improving readability. For example, using the optimization tool cvxpy, the stress constraints in Eq. 3b are expressed as q[k] >= -sc * a and q[k] <= st * a.

Check dual violation
When the primal problem (3) is solved via the primal-dual interior point method, the virtual displacements in u for each load case are obtained from the equilibrium constraints defined in eqn. Therefore, violation of constraint (4) Fig. 4 Cantilever example: problem definition can be calculated for every member in the full 'ground structure' (except for those already present in the reduced set, where this is unnecessary as there will be no violation). All potential members are then sorted by their constraint violation numbers, and the most violating ones are added to the primal problem, with the number of new members calculated using: m P and m V are, respectively, the numbers of remaining members in the potential member list (i.e. members not yet included in Eq. 3), and those violating (4). α and β are control parameters, both set to 0.05 in this paper. Using the adaptive 'member adding' scheme, the size of problem (3) gradually increases, with Eq. 7 used to control the growth rate. When number of violations is large, only a small proportion are added into (3); on the other hand, when the number of violations is small, most (or all) are added in a given iteration. This balances problem size with the number of iterations required in order to achieve an efficient iterative process. Various alternative heuristics can be used to determine the priority and number of members to be added at each iteration, e.g. the two-stage high-pass filter in Gilbert and Tyas (2003) and the active member set method described in Sokół (2014). It is important to note that, although a heuristic algorithm is used here to improve computational efficiency, the final solution, obtained when (4) is satisfied for all members in the ground structure, is equivalent to the solution obtained by solving the full problem, due to the strong duality nature of LP problems (Vanderbei 2001).

Visualization
The solutions are visualized in function plotTrusses. Members with areas greater than a threshold number (e.g. 1 × 10 −3 of the maximum member area) are plotted in red or blue, depending on the sign of their internal forces. When the signs of internal forces vary in different load cases, grey is used instead.

Main workflow
The main trussopt function contains several steps to define the design domain, load, and support conditions for any given layout optimization problem, to create the 'ground structure', and to solve the LP problem (3c) using the adaptive 'member adding' scheme.
-Firstly, a polygonal design domain is specified using the Polygon function imported from shapely, and then its convexity checked. (c) (d) -Secondly, a uniformly distributed grid of nodes is created using the meshgrid function, with points lying inside the design domain added to node list Nd; then support and load conditions are added. -Thirdly, the 'ground structure' is created by generating all valid potential members, excluding overlapping lines (checked using the greatest common divider, function gcd, when joint costs are not specified), and lines crossing the polygonal boundary (using function contains from shapely). -Finally, the initial reduced set of members is created by including only short member connections. The adaptive 'member adding' process then starts, in which potential members are added to this reduced set until no violation of Eq. 4 is found in any member in the 'ground structure'.

Numerical examples
A number of examples serve to demonstrate the efficacy of the script. Firstly, a simple cantilever problem (Fig. 4) can be solved by running the following command in a terminal:   Alternatively, the script can be imported as a module in Python using: And then function trussopt can be called to perform a layout optimization with the following input arguments: which solves the simple cantilever problem shown in Fig. 4 with width = 20, height = 10, limiting stress in tension = 1, in compression = 1 and joint cost = 0. In this section, all line numbers refer to the script in Appendix C, and all quoted CPU times were obtained using a laptop PC equipped with an Intel I7-7700HQ CPU and running 64-bit Windows 10.

Effect of using the adaptive 'member adding' scheme
Using the adaptive 'member adding' scheme, a large-scale layout optimization problem is solved as a series of much smaller LP problems solved in succession; e.g. Fig. 5 shows steps in the solution of the cantilever design problem shown in Fig. 4.
When solving medium or large-scale problems, the adaptive 'member adding' scheme reduces memory consumption and can result in significant speed improvements, as indicated in Table 1. Whilst the resulting structural volumes are identical (to five significant figures), the latter becomes increasingly efficient with increasing problem size, e.g. 8.2 times faster when 60 × 30 nodal divisions are used. Note that 'member adding' can be switched off by increasing the maximum initial member connection distance in line 85, for example changing this to: for the purposes of this paper.

Alternative LP solvers
By default, the LP problem (3c) is solved using the opensource solver ECOS (Domahidi et al. 2013) distributed with cvxpy. Alternatively, a range of other, potentially more efficient, LP solvers such as GUROBI and MOSEK can be used with cvxpy. To use MOSEK, for example, line 29 is replaced with the following (note that this requires MOSEK has already been added to Python environment; if not, guidance can be found in Appendix B): The MOSEK parameter "MSK IPAR INTPNT BASIS" is optional; however this disables the unnecessary basis identification step to improve speed. Using MOSEK, the CPU cost required to solve problems is significantly reduced, making the layout optimization process very efficient. Table 2 shows sample solutions for the cantilever

Joint cost
To approximately account for the costs associated with joints, a joint cost can be directly specified when calling the main function trussopt. For example, a joint cost of s = 1 unit length can be added using the following command; the simpler solution then generated is shown on Fig. 6.

Unequal stress limits
Unequal stress limits can also be directly specified. For example, a compressive stress limit can be set to σ − = 0.1σ 0 with the following input arguments: Alternatively, a different tensile stress limit can be set if desired. The corresponding results are shown in Fig. 7.

Change load and support conditions
Load and support conditions can be changed by modifying lines 73 and 74. For example, to generate the 'half wheel' example shown in Fig. 8, the pin and roller supports can be applied by replacing line 73 with: which locates the two corner nodes at the base of the domain and then sets the appropriate degrees of freedom. Similarly, the point load can be applied by replacing the original line 74 with:

Multiple load cases
Additional load cases can be included by appending load vector f after line 74. For the two load case cantilever example problem shown in Fig. 9, load case Q is added using the following lines: And the point load P is moved to the bottom-right corner by replacing the original line 74 with: Initially all load cases are concatenated in vector f, which is then reformatted to a p × 2n array in line 83.

Non-convex domain
Various polygonal domains can be specified using shapely. For example, to add a hole to the design domain, as shown in Fig. 10, the following lines can be added after line 65: This subtracts a rectangular region from the original design domain in Fig. 4. Note that, when a non-convex domain is used, lines that intersect domain boundary should be excluded from the 'ground structure'. Line 80 checks such intersections using the contains function from shapely. Since this can become computationally expensive when the number of potential connections is large, it is skipped when a convex domain is used.

Conclusions
A simple Python implementation of truss layout optimization using an adaptive 'member adding' scheme has been presented for educational use. Various examples are used to demonstrate the efficacy of the script for problems involving single and multiple load cases, unequal limiting stresses in tension and compression, joint costs, and non-convex polygonal domains, with or without holes. The solutions obtained are globally optimal for the prescribed nodal discretization and can be used to benchmark other methods or to provide inspiration for structural designers.
The 98-line Python script described is listed for reference in Appendix B and is also provided in the form of downloadable source code; see Supplementary material section.

B.2 Installation with Python 2
It is possible to run the script with Python 2 (2.6 and later) by making minor modifications to the script, as detailed below.
First, remove the gcd method from math in the first line: And instead import it from fractions instead, by adding the following new line: Also, line 47 is replaced with: to take into account the fact that function ceil returns a floating point rather than integer value in Python 2.

B.3 Alternative solvers
To use an alternative solver, as described in Section 4.2, this also needs be installed. For example, MOSEK can be installed using the following command: In addition, a licence file is required, which can be freely obtained for academic use from the MOSEK website (https://www.mosek.com).

B.4 Python IDE
For beginners, the Spyder IDE (Integrated Development Environment) bundled with the Anaconda package provides an easy-to-use means of editing and running the script.