Adaptive solution of truss layout optimization problems with global stability constraints

Truss layout optimization problems with global stability constraints are nonlinear and nonconvex and hence very challenging to solve, particularly when problems become large. In this paper, a relaxation of the nonlinear problem is modelled as a (linear) semidefinite programming problem for which we describe an efficient primal-dual interior point method capable of solving problems of a scale that would be prohibitively expensive to solve using standard methods. The proposed method exploits the sparse structure and low-rank property of the stiffness matrices involved, greatly reducing the computational effort required to process the associated linear systems. Moreover, an adaptive ‘member adding’ technique is employed which involves solving a sequence of much smaller problems, with the process ultimately converging on the solution for the original problem. Finally, a warm-start strategy is used when successive problems display sufficient similarity, leading to fewer interior point iterations being required. We perform several numerical experiments to show the efficiency of the method and discuss the status of the solutions obtained.

The aforementioned contributions which include Euler buckling constraints are intended to avoid slender bars in compression being included in the solution while including nominal forces are designed to ensure nodes connecting bars in compression are adequately braced. Formulations which directly include global stability constraints are concerned with ensuring the stability of the whole structure. This is of interest because, even if a structure is well-braced, it may fail as a result of insufficient overall elastic stiffness. Note that global stability problem formulations implicitly address the nodal instability problem, though do not take into account local instabilities of the sort dealt with by the Euler formula (e.g. Levy et al. 2004).
The focus of this paper is on problems with global stability constraints. These problems are in general formulated as nonlinear and nonconvex semidefinite programs (Ben-Tal et al. 2000;Levy et al. 2004;Stingl 2006;Evgrafov 2005). Such problems are computationally challenging and for large-scale structures are usually considered numerically intractable, though software capable of solving small problems is available (Fiala et al. 2013;Kočvara and Stingl 2003). For this reason, some studies formulate the nonlinear semidefinite programming problem as an equivalent nonlinear programming problem (Tugilimana et al. 2018). Nevertheless, in all approaches, the size of problems that have been solved thus far in the literature has either been small or otherwise limited, e.g. by only specifying minimum connectivity between the nodes in the design space. This is in stark contrast to the plastic design formulation, solvable via linear programming, when full nodal connectivity problems can be solved even for high nodal densities.
In this article, we propose a relaxation of the nonlinear and nonconvex semidefinite programming formulation, for which we develop an efficient optimization algorithm based on interior point methods. The method is coupled with other novel techniques to make it capable of solving large-scale problems and with full nodal connectivity. The relaxed problem is still a semidefinite program but ignores the kinematic compatibility constraints present in standard elastic formulations. We observe huge computational gains by solving the relaxed formulation and provide lower bounds for the associated general nonlinear problems. Moreover, we report an estimation of the violation of the removed kinematics compatibility equations by solving an associated least-squares problem. For some small-scale benchmark problems, we additionally make comparisons between solutions of the relaxed and original nonlinear problems. In general, the error due to ignoring the kinematic compatibility equations are observed to be very small for reasonable values of the stability load factor, suggesting that solutions to the relaxed problems are acceptable. However, as we increase the value of the stability load factor beyond practically realistic values, we observe a high degree of violation in the compatibility equations and significant differences in the optimal designs. For some of the examples presented in this article, we also calculate the violation of the elastic stress constraints for the solution of the relaxed semidefinite program. The numerical results once again confirm the usefulness of the solutions obtained, when stability load factors of practical interest are used. Next, we describe the techniques that contribute to the efficiency of the proposed solution algorithm.
Firstly, we employ the adaptive 'member adding' approach, previously used to solve plastic truss layout optimization problems (Gilbert and Tyas 2003;Sokół and Rozvany 2013;Weldeyesus and Gondzio 2018) via linear programming, though now solving the problems of interest here via semidefinite programming. It is a procedure in which we approximate the large-scale original problem by a sequence of smaller subproblems, the solutions of which ultimately converge to that of the original problem. In the context of truss optimization, this is done by first solving the problem with minimum nodal connectivity, followed by generating and adding more members/bars based on degree of violation of constraints in the dual problem. The procedure continues until the solution to the original problems is obtained. This procedure greatly reduces the memory required to solve a given problem and, even using a standard desktop computer, we have managed to solve problems that otherwise would require hundreds of GB memory. Detailed statistics are presented in Section 6; see in particular the large-scale bridge example problem described in Section 6.3.1.
Secondly, similar to Ben-Tal and Nemirovski (1997), we explicitly utilize the structures of the problems, i.e. the high degree of sparsity and low-rank property of the element stiffness matrices (Bendsøe et al. 1994;Ben-Tal 1993b;Achtziger et al. 1992;Bendsøe and Sigmund 2003), to address the computational bottle-neck associated with using the interior point method to solve semidefinite programming problems. This determines the coefficient matrix of the linear system originating in the algorithm. Roughly speaking, instead of performing O(mn 3 +m 2 n 2 ) arithmetic operations (Fujisawa et al. 2000), by using standard and straightforward expressions to determine the matrices involved in the linear systems, we perform O(m 2 n) arithmetic operations, where m is the number of bars and n is number of nodal degrees of freedom. Note that the sparsity of the element matrices is also effectively used in performing matrix inner products in the adaptive member adding procedure, as described in Remark 6.
Finally, as when solving the plastic truss layout optimization using the interior point method (Weldeyesus and Gondzio 2018), we apply a warm-start strategy to solve some of the subsequent problems, determining an initial point that reduces the number of interior point iterations and overall improves the convergence characteristics of the optimization process. The technique relies on an observation that the number of newly added bars decreases towards the end of the adaptive member adding procedure and therefore the degree of similarity between successive subproblems increases at this stage.
The paper is organized as follows. In Section 2, we present the general nonlinear and nonconvex semidefinite programming model of the truss layout optimization problem with global stability constraints, its relaxation, and the least-squares problem used to estimate violation of the kinematic compatibility constraints. We describe the general framework of the primal-dual interior point method and exploitation of the structure of the matrices in Section 3 and the adaptive member adding procedure in Section 4. The warm-start strategy and related mathematical analysis are presented in Section 5 and numerical experiments are described in Section 6. Finally, conclusions and possible future research directions are presented in Section 7.

The problem formulation with stability constraints
In this section, we describe the formulation for the layout optimization of trusses with global stability constraints problem. We use the 'ground structure' approach (Dorn et al. 1964) to formulate all problems. This is done by distributing a finite set of nodes, say d, across the design space and connecting these nodes by all possible potential bars, including overlapping ones. Hence, we have m = d(d − 1)/2 bars, where clearly m d. We denote the cross-sectional areas of the bars by a i , i = 1, ..., m. Let f ∈ R n , = 1, . . . , n L be a set of external forces applied to the structure where n (≈ Nd, N is the dimension of the design space, i.e. 2 or 3) is the number of the non-fixed degrees of freedom. Then, the associated (nodal) displacements u ∈ R n , = 1, . . . , n L satisfy the elastic stiffness equation where the stiffness matrix K(a) is computed as and the element stiffness matrices K i 's are given by with γ i ∈ R n being the vector of direction cosines for the ith bar. Introducing the axial forces q in member i which are given by allows us to rewrite (1) as where B = (γ 1 , · · · , γ m ) ∈ R n×m . Next, we define the geometric stiffness matrix G(q ) as given by where in which the vectors δ i , η i are determined so that γ i , δ i , η i are mutually orthogonal (see Kočvara (2002) for details). The vectors δ i and η i are not necessarily unique. These are chosen in Kočvara (2002) as the orthogonal basis of the null space of γ T i and we follow similar approach in our implementation. Now, the multiple load case minimum weight (or, strictly speaking, minimum volume) truss layout optimization problem with global stability constraints can be formulated as minimize a,q ,u l T a subject to where l ∈ R m is a vector of lengths of the bars, and σ + > 0 and σ − > 0 are the material yield stresses in tension and compression, respectively. Note that we can find a similar formulation to problem (8) in Tugilimana et al. (2018) with additional constraints enforcing local (Euler) buckling. The parameter τ can be interpreted as a stability load factor and must be set toτ ≥ 1, ∀ to indicate that the resulting optimal structure is stable for the loads f , ∈ {1, . . . , n L }.
It is worth mentioning that problem formulation (8) does not address local buckling, as addressed, e.g. by the Euler buckling equation. Therefore, we can expect that the optimal design obtained by solving problem (8) could potentially include slender bars (Levy et al. 2004).
Due to the inclusion of the nonlinear kinematic compatibility equations (4), problem (8) is a nonlinear and nonconvex semidefinite programming problem. Such problems are in general very difficult to solve. In Fiala et al. (2013) and Stingl (2006), methods for treating nonlinear semidefinite programming problems are described in which a variant of formulation (8) is solved. These formulations are very attractive, but the challenge of solving large-scale problems still remains. In Tugilimana et al. (2018), problem (8) is transformed from a semidefinite programming problem to a standard nonlinear programming problem.
In this paper, we relax problem (8) by ignoring the kinematic compatibility constraint (4), and solve a semidefinite programming problem of very large dimension because we allow full nodal connectivity. Namely, we consider minimize a,q l T a subject to which is a (linear) semidefinite program and can be interpreted as the plastic design formulation with global stability constraints. In our numerical experiments described in Section 6, we additionally report the maximum violation of the kinematic compatibility constraints for the optimal designs obtained by solving the relaxed problem (9). This violation is estimated by solving the least-squares problem where a * and q * are the solutions of the relaxed problem (9). Note that the special case τ = 0, ∀ , or in other words, excluding the matrix inequality constraints, reduces problem (9) to the plastic layout optimization problem, which is a linear program that can be solved efficiently by an interior point method (Weldeyesus and Gondzio 2018).

Remark 1
The relaxed problem (9) belongs to the class of linear semidefinite programming problems. Hence, any of its solutions are also globally optimal solutions. Moreover, this provides a (strict) lower bound to nonconvex problem (8) for τ > 0.

Remark 2
The relaxed problem (9) can be solved very efficiently by extending the adaptive 'member adding' scheme which has been used previously to solve large-scale plastic layout optimization of trusses formulated as linear programs (Gilbert and Tyas 2003;Sokół and Rozvany 2013;Weldeyesus and Gondzio 2018).
Remark 3 For some of the small-scale examples in Section 6, we additionally solve the nonlinear semidefinite program (8) and compare the solutions obtained with those obtained using the linear SDP relaxation (9). To solve the nonlinear semidefinite program (8), we have implemented a standard interior point method that uses outer and inner loops, where the inner loop is used to solve first-order optimality conditions for a fixed value of a barrier parameter. Then, we allow a very small reduction in the barrier parameter in the outer loop. This of course means that many iterations are required and convergence is slow. Nevertheless, this has proved effective for the small-scale problems considered herein. Note that the member adding procedure described in Section 4 is not used when solving the nonlinear semidefinite program (8) as the assumptions needed to apply the techniques of Section 4 in general only hold for the linear SDP relaxation (9).
Remark 4 The least-squares problem (10) always has an objective value of 0 for a single-load case problem when τ = 0 in (9). This is because for τ = 0, the problem (9) precisely reduces to the so-called least-weight (or minimum volume) plastic design problem which has indeed been shown to be equivalent to the elastic minimum compliance problem (Hemp 1973;Achtziger et al. 1992;Ben-Tal and Bendsøe 1993;Bendsøe and Sigmund 2003;Achtziger 1996;Stolpe and Svanberg 2004).

The primal-dual interior point framework
We adopt the Mehrotra-type primal-dual predictor-corrector interior point method (Fujisawa et al. 2000) for semidefinite programming.
Introducing the slack variables s − , s + ∈ R m + , and S ∈ S n + , we rewrite (8) as and write down the following dual: where λ ∈ R n denotes the virtual nodal displacement, x + , x + ∈ R m , = 1, . . . , n L , x a ∈ R m , X ∈ S n + , and the Remark 5 The primal problem formulation (11) is traditionally referred to as the dual problem, and the dual problem formulation (12) is traditionally referred to as the primal problem in literature on semidefinite programming (Wolkowicz et al. 2000).
where (borrowing MATLAB notation)B= blkdiag(B,..., B), and D kl are diagonal matrices. The vector (ξ 1 , ξ 2 , ξ 3 ) T is the resulting right-hand side. For a complete description of the interior point method, we refer the reader to Fujisawa et al. (2000). The rest of this section is dedicated to the computational difficulties associated with the interior point method for semidefinite programming and the techniques we use to resolve these. There are several computational challenges associated with the linear system (14). Firstly, all block matrices A kl , k, l = 1, ..., n L are dense and require a large amount of memory to store them (see Fig. 1a for the sparsity structure of the coefficient matrix for a two-load case problem). Secondly, the straightforward computation of the coefficient matrix requires O(mn 3 + m 2 n 2 ) operations (Fujisawa et al. 2000).
The second challenge can be easily resolved by exploiting the low-rank property and sparsity of the data matrices K i , G i , i = 1, ..., m (Ben-Tal and Nemirovski 1997; Bendsøe et al. 1994;Ben-Tal 1993b;Achtziger et al. 1992;Bendsøe and Sigmund 2003). From (3) and (7), it can be seen that the rank of the element stiffness matrices K i , i = 1, ..., m is always 1 and the rank of the element geometric stiffness matrices G i , i = 1, ..., m is 1 for two-dimensional problems and 2 for three-dimensional problems. The direction cosine vectors γ i , δ i , and η i are all very sparse with at most 4 or 6 nonzero entries for twoand three-dimensional problems, respectively. Therefore, we utilize this property to compute the coefficient matrix  (14) for a two-load case problem. a Without adaptive member adding (number of non-zeroes = 73,893). b With adaptive member adding, and in the final SDP iteration (number of non-zeroes = 24,126) efficiently. For example, consider the block matrix A 11 (single-load case for notation simplicity). We have which can be computed in O(n) arithmetic operations and hence the computation of the coefficient matrix can be brought down to O(m 2 n).
In the next section, we discuss the novel approach employed to deal with the first challenge, i.e. the large memory requirements.

Adaptive 'member adding'
In problem formulation (9) we consider fully connected ground structures. Hence, for a N-dimensional problem comprising d nodes, the coefficient matrix (KKT) of the reduced Newton system (14) has dimension where m = d(d − 1)/2 and n ≈ Nd. Moreover, all of the larger blocks of the coefficient matrix, each with dimension m × m, are full matrices. This indicates that it would be computationally prohibitive to store or factorize the matrix (see also the large-scale bridge example problem described in Section 6.3.1).
In order to overcome this we extend the adaptive 'member adding' approach initially proposed for linear plastic truss layout optimization problems by Gilbert and Tyas (2003) and later used in other studies, for example by Sokół and Rozvany (2013) and Weldeyesus and Gondzio (2018). It is a strategy whereby the original large problem is solved by successively solving a number of smaller subproblems, i.e. problem instances with fewer bars, saȳ m. In practice,m m, hence each of the m × m dense blocks in the KKT matrix in (14) reduces to smallerm ×m blocks. Moreover, as the problem size increases (i.e. when higher nodal densities are used), the fraction of the bars m/m used in the SDPs corresponding to the final member adding iterations decreases (see Remark 10 in Section 6.2). An overview of the member adding strategy, which can only be applied to the relaxed linear SDP (9), is provided below.
First, we start with a structure constituting minimum connectivity, for example with the structure shown in Fig. 2a for two-dimensional problems and Fig. 2b for threedimensional problems, and let m 0 be the number of bars in the the initial structure. We denote by K 0 ⊂ {1, . . . , m} the set of indices of the bars for which the primal problem (11) and its dual (12) are currently solved. Next, we compute the dual violations using only variables λ and X in (12) which are described below.
For any member i to be dual feasible (see (12)), we need Now, since x + ≥ 0 and x − ≥ 0, from (17a), we have and from (17b), we have Combining (18) and (19), we get for a member i to be dual feasible. Any member that violates (20) is said to be dual infeasible. Now, solving the problem for members with indices in K 0 , we use (20) to generate the set K as where with λ * and X * being optimal values, and β > 0 some prescribed tolerance. Then, we identify the bars with indices in K, filter them, and then finally add then to form the next problem. The purpose of the filtering is to limit the number of bars to be added in order to prevent fast growth of the size of the problem (for details of heuristic filtering approaches, see Weldeyesus and Gondzio 2018). Note that use of a different filtering approach may affect the size of each problem to be solved as part of the adaptive member adding process, and the number of iterations required, but not the final solution. For the numerical experiments in Section 6, we use the member bar length approach described as filtering strategy AP3 in Weldeyesus and Gondzio (2018). The member adding procedure terminates when K = ∅.
Remark 6 In our implementation, the sparsity of the data matrices K j and G j is exploited to determine the set K in (21) while performing the operations K j • X * and G j • X * . Hence, this step becomes inexpensive. The CPU times reported for the numerical experiments in Section 6 include this procedure.
Remark 7 In Fig. 1, we present the sparsity and size of the coefficient matrix of the reduced Newton system (14) for a small problem. Figure 1a shows the situation when the problem is solved for all potential bars, and Fig. 1b shows the situation when we apply the adaptive member adding strategy. The sparsity structures may look similar but the size is reduced. Moreover, this reduction in size becomes even more significant for larger problems (see the largescale bridge example problem described in Section 6.3.1).

Warm-start strategy
After performing several member adding iterations, the subsequent subproblems start to become more and more similar. Therefore, we use a warm-start strategy and determine an initial point that can reduce the number of interior point iterations required to obtain a solution. This has been used for the basic truss layout optimization problem formulated as a linear program in Weldeyesus and Gondzio (2018) and is now applied to the semidefinite programming formulations presented in this paper. The discussion in this section and the mathematical analysis closely follow Section 6 of Weldeyesus and Gondzio (2018).
As described in Section 4, we generate the set K in (21) at every member adding iteration. If K = ∅, then we form the new problem in which the variables are where the vectors with the super-bar, all in R k , k = |K|, represent the new variables corresponding to the newly added bars. We assume that the old solution (the left-hand side in (22)) was feasible in the previous instance of the adaptive member adding scheme.

Computing a warm-start point
We set the initial point for the variables without the superbar in the right-hand side of (22) to the solution of the previous problem instance obtained with a loose tolerance. Following Gondzio (1998), the interior point algorithm should not be initialized at a point too close to the boundary of the feasible region. For the newly added variables, i.e. those with bars in (22), we use a specialized initialization procedure given below. We first setx + andx − as where μ 0 is the value of the barrier parameter corresponding to the saved solution of the previous problem instance. Its value is computed as where m 0 = |K 0 |. Then, we set the new dual slack variable as Finally, new primal variables are defined bȳ Now, similar to Weldeyesus and Gondzio (2018), we estimate the bounds on the primal and dual infeasibilities, and the violation in complementarity slackness conditions induced by the new variables.
Now, we determine the bound for the last primal infeasibility ξ p 4, (13f).
This is because where the last inequality above holds sincex a,i ≥ μ 1 2 0 by definition. Moreover, ||γ iγ T i || ∞ ≤ 2N . This is because, for example, when N = 2 the non-zeros entries of the direction cosine γ i are γ i = (− (l x ) i l i , − (l y ) i l i , (l x ) i l i , (l y ) i l i ) which implies ||γ iγ T i || ∞ ≤ 4 . Therefore, the expressions in (27) and (28) demonstrate that primal infeasibility is at worst proportional to μ 1 2 0 , and hence insignificant.
Using the definition ofx a in (25) and the fact that r − |r| ≤ 2|r|, r ∈ R, we have where σ max = max{σ − , σ + }, and ε − and ε + are given as (21). Hence, where (K X ) j =K j •X , ∀j ∈ K. This expression reveals that there can be a considerable violation of the first dual constraint in (12) and so we apply the warm-starting routine presented by Gondzio (1998), Gondzio and González-Brevis (2015) to address this.

Centrality
We compute complementarity products for all newly added variables to evaluate the centrality of the new point. Note that the last centrality in condition (3) is automatically satisfied. Moreover, the pairs (ā,x a ) are μ 0 centred from (26). Since Next, finding upper bound on (x a ) j Then, using (33) and (34) σ Similarly, The bounds in (35) and (36) imply that the pairs (x − ,s − ) and (x + ,s + ), ∈ {1, ..., n L } have no prominent outliers from the μ 0 -centrality. This is due to the shift-like terms involving (ε − j + ε + j ) in the right-hand side being multiplied by μ 1 2 0 , thus reducing the induced violation of the μ 0centrality.

Numerical results
The interior point method has been implemented in MATLAB (R2016a). All numerical experiments have been performed using a PC equipped with an Intel Core™ i5-4590T CPU running at 2.00 GHz with 16 GB RAM. For the case of the cold-start runs, the initial points a 0 , s +0 , s −0 , x +0 , x −0 , x + a are set to unity, S , X to the identity matrix I , and q 0 , λ 0 to zero. The interior point algorithm terminates when wheref = (f 1 , ..., f n L ) T , matrix terms are vectorized, and the primal and dual residuals ξ p and ξ d are given by (3). Note that, for feasible primal and dual points, the duality gap can easily be written as by performing some elementary algebraic operations. The primal and dual relative feasibility tolerances are set to p = d = 10 −6 .
For the optimality tolerance, we use loose tolerances in the first few member adding iterations since the first few subproblems should not have to be solved to optimality, and then tighter tolerances in the final iterations, i.e. opt = [10 −2 , 10 −2 , 10 −3 , 10 −4 ] and then always 10 −5 . The reported CPU times correspond to the entire solution process, including the member adding computations.
In the original problems, we consider all the potential bars, including overlapping bars. At the start of the solution process, we begin with the structure shown in Fig. 2a for two-dimensional problems and Fig. 2b for threedimensional problems, and use β = 0.001 to generate the set of members in K given in (21) that are dual infeasible. If the warm-start strategy is used, then it is activated at the fourth member adding iteration or else before this if (m k − m k−1 )/m k ≤ 0.12, where m k is the number of bars used in the kth member adding iteration. In applying this strategy, we use the solutions obtained with tolerance ε opt = 0.1 in the preceding problem instance to determine the initial point of the subsequent problem.
For all examples, λ min denotes the smallest positive eigenvalue of the generalized eigenvalue problem Moreover, we set τ = 0 in (8) and (9) for problems without stability constraints and τ ≥ 1 for problems with stability constraints. Its specific values are mentioned in the examples below. Additionally, when reporting the solution of the SDP relaxation (9), we also provide an estimate of the violation of the compatibility equations by solving the least-squares problem (10). Finally, we use Young's modulus E = 210 GPa, and equal tensile and compressive strengths of 350 MPa. In the plots of the optimal designs, bars in tension are shown in red and bars in compression are shown in blue (except for sake of clarity in the case of Fig. 13). In all cases, the bars shown are those with cross-sectional area ≥ 0.001a max and the dark dots are the active nodes connecting these bars.

Benchmark example problems
The objective of the examples in this section is to provide an insight into the solution obtained when a linear SDP relaxation (9) is solved for benchmark problems reported in the literature. This is done by comparing solutions with those obtained using the nonlinear SDP (8), which includes compatibility constraints. We display the corresponding optimal designs in Figs. 3 and 5 for τ = 1 and τ = 10 that are likely to be of interest to engineering practitioners. However, we also show numerical results in Tables 1 and 2 for larger values of τ , to show the growth of the violation of (4) in extreme cases.
For the examples described in this section, we also calculate the violation of the elastic stress constraints for the solution of the relaxed SDP (9). This is in addition to the least-squares approach for estimating violation of the kinematic compatibility equation, and is done by computing the stress values as σ * ,i = E l i γ T i u * , where u * = K −1 (a * )f , with a * being the solution of the relaxed SDP (9).

L-shaped truss example
We solve the benchmark L-shaped truss example problem shown in Fig. 3a, comprising 132 bars. It has dimensions 1 m × 3 m × 4 m (including the null region of dimensions 1 m × 2 m × 3 m), and each of the applied nodal loads is 350 kN, applied simultaneously. The optimal designs are given in Fig. 3b-f and resemble the solutions presented in Levy et al. (2004) and Tugilimana et al. (2018), who solved various problem formulations incorporating global stability constraints, and those presented in Tyas et al. (2006) and Descamps and Coelho (2014), who solved problems incorporating destabilizing nodal forces.
When the problem is solved without stability constraints, the solution shown in Fig. 3b is obtained which comprises two parallel planar trusses. In that case, the optimal design has volume 0.0620 m 3 and λ min = 1.2e−05 < 1; hence, it is not stable. Next, solving the relaxed problem (9) with stability constraints for τ = 1, we obtain the solution shown in Fig. 3c, where a connection between the two parallel planes is now established. The volume of the structure is 0.062217 m 3 and λ min = 1. It is useful to establish an estimation of the violation of the kinematic compatibility equation, obtained from (10); this is found to have a value of 4.9624e−06. Moreover, we can compare the solution shown in Fig. 3c to that of the optimal design shown in Fig. 3d, obtained when solving the standard nonlinear SDP formulation (8). In this case, the design has a marginally higher volume of 0.062241 m 3 . In order to evaluate the result obtained by solving the relaxed SDP (9) and the nonlinear problem (8) in more detail, we also solve the problems for a larger value of the loading factor τ . Thus for τ = 10, the solution to the relaxed SDP gives the design shown in Fig. 3e which has a volume 0.064333 m 3 and the violation of the kinematic compatibility equation is equal to 5.3177e−04 . This is larger than in the case above when τ = 1. For τ = 10, the nonlinear SDP gives the design shown in Fig. 3f which has a somewhat larger volume, of 0.064639 m 3 . Table 1 shows the behaviour for even higher values of τ , i.e. τ = 20, 30, ..., 90. This indicates that when the value of τ is increased, the magnitude of the violation of the kinematic compatibility constraints increases. Nevertheless, the results seem to agree for small values of τ and especially for the required minimum value τ = 1, so that the structure remains stable when the load is applied.
Next, we calculate violation of the elastic stress constraints in the relaxed SDP (9) solutions following the procedure described above. We found maximum violation of the elastic stress constraints to be 0.35% when τ = 1 and 3.6% when τ = 10. Moreover, the resulting load factors (the smallest positive eigenvalues of (38)) were found to be 0.9983 and 9.8310, respectively, as shown in Fig. 4a and b. The elastic stress distributions are plotted in Fig. 4, where the values of the stress are within the specified limits in the grey bars, but violated slightly in the blue (in compression) and red (in tension) bars.   Table 2 Tower example: comparison of volumes obtained obtained by solving the nonlinear SDP (8) and the SDP relaxation (9), and violation of the compatibility constraints (4) estimated by solving the least-squares problem (10) Remark 8 The solution to the problem without stability constraints presented in Fig. 3a constitutes not only two independent planar trusses but also unstable nodes connecting bars that are in compression. The unstable nodes are stabilized in Fig. 3c-f with bracing bars.

Tower example with downward vertical load
We solve the tower example problem shown in Fig. 5a, comprising 1953 bars in the fully connected ground structure. This is motivated by the similar problems solved by Stingl (2006) and Tyas et al. (2006). In this example, the tower is for sake of simplicity assumed to have dimensions of 1 m × 1 m × 3 m, is fixed at its base, and is subjected to a downwards vertical load of 350 kN at the centre of its upper surface.
The optimal design is shown in Fig. 5b for the problem without stability constraints, which turns out to comprise six vertical inline bars with no bracing elements. Its volume is 0.00300 m 3 and λ min = 2.3277e−04; clearly, this structure is not stable. Now, setting τ = 1 and solving the relaxed problem with stability constraints (9), we obtain the solution shown in Fig. 5c with no intermediate unstable nodes, with  (9). The values of the stress are within the specified limits in the grey bars, but violated in the blue (in compression) and red (in tension) bars by at most 0.35% when τ = 1 and 3.6% when τ = 10 (actual load factors for relaxed SDP structures a 0.9983, b 9.8310) bracing bars connecting the loaded node. The stable design has a volume of 0.003010 m 3 . Estimating the violation of the kinematic compatibility equation, we solve problem (10) and get 3.6617e−06. The nonlinear formulation (8) produces the solution shown in Fig. 5d and has a volume of 0.003020 m 3 .
For a higher value of τ = 10, the solution to the relaxed SDP (9) returns the design shown in Fig. 5e, with a volume 0.003102m 3 and compatibility constraint violation of 5.0965e−05, and the nonlinear SDP (8) returns the design presented in Fig. 5f with volume 0.003200 m 3 . Note that the results are in broad agreement with the results obtained by Stingl (2006) and Tyas et al. (2006), where the tower problem is respectively solved using a nonlinear semidefinite formulation with global stability constraints and with the introduction of destabilizing nodal forces. In Table 2, results are presented for higher values of τ =  Finally, we calculate violation of the elastic stress constraints in the relaxed SDP (9) solutions. We found maximum violation of the elastic stress constraints to be 0.33% when τ = 1 and 0.6% when τ = 10. Moreover, the resulting load factors (the smallest positive eigenvalues of (38)) were found to be 0.999965 and 9.9994, respectively, as shown in Fig. 6a and b. The elastic stress distributions are plotted in Fig. 6, where the values of the stress are within the limit in the grey bars, but violated slightly in the blue (in compression) bars.
Remark 9 As mentioned in Section 2, models (8) and (9) do not address local buckling. This is shown in Fig. 5c and f where long bars in compression are used as bracing or as means of stabilizing otherwise unstable nodes.

Tower example with upwards vertical load
The purpose of this example is simply to demonstrate that if the bars in the optimal design are all in tension, then the solution obtained with or without the global stability constraints are identical. To show this, we solve the problem in Example 6.1.2 but with the direction of the load reversed, as shown in Fig. 7a. The optimal design is shown in Fig. 7b for the problem without stability constraints and once again comprises six vertical inline bars, all in tension and with no bracing elements. Its volume is 0.00300 m 3 and λ min = 426.5575 > 1. This shows that the design is already stable and setting τ = 1 and re-solving the problem with stability constraints (9) will change neither its volume (which is is 0.00300m 3 ) nor its geometry, as can be seen in Fig. 7c.

Adaptive 'member adding' problems
Here, we report on the efficacy of the adaptive member adding strategy described in Section 4 for the relaxed linear SDP (9). This is achieved by solving problems both with and without the strategy, verifying that the same solution is obtained, and reporting on comparative computational efficiency.

Bridge example (small-scale)
We solve the bridge-like example problem shown in Fig. 8a, comprising 3240 bars in the fully connected ground structure. The design domain has dimensions 8 m × 2 m × 2 m and has fixed pin supports at each of the four corner nodes. Vertical loads of magnitude 350 kN are applied to all nodes along each the two long edges at the base of the domain.
The solution obtained when stability constraints are not included is shown in Fig. 8b, comprising two parallel planar trusses. In this case, the optimal structure has a volume of 0.0540 m 3 and λ min = 3.8613e−08 (i.e. clearly not stable). Next, we solve the problem with the stability constraint (9) for τ = 1. Numerical results are shown in Table 3. Figure 8c shows the optimal design when solving the entire original problem and Fig. 8d shows the structure obtained when member adding is used. The optimal designs are clearly identical and have the same volume, equal to 0.05414 m 3 (see row 1 of Table 3). Moreover, the CPU times reported in the table illustrate the efficiency of the member adding scheme. In general, these efficiencies are much more pronounced for larger problems. Figure 9 illustrates the evolution of the solution when the adaptive member adding strategy is used, showing the potential bars and the corresponding optimal design for each member adding iteration. The violation of the kinematic compatibility constraint is found to equal 5.8336e−06 at the end of the process.  9 Bridge example (small-scale): potential bars and optimal designs illustrating the evolution of the optimal designs with respect to the member adding iterations reported in Table 3. Green bars represent newly added bars

Large-scale problems and warm-start strategy
We now solve large-scale problems in which we additionally demonstrate the numerical benefit of using the warm-start strategy described in Section 5.

Bridge example (large-scale)
We consider again the bridge problem, though now with 90,100 bars, as shown in Fig. 12a; the loading conditions and dimensions are as described in Section 6.2.1. It is worth mentioning that if we attempted to solve the original problem without member adding, then we would need approx. 240 GB of memory to store the coefficient matrix in (14).
However, by applying the adaptive member adding technique we not only reduce the CPU time but also significantly reduce peak memory requirements. Numerical results are presented in Table 4. It is evident that the warm-start strategy reduces CPU time by approx. 30% when τ = 1 and by 53% when τ = 10 , which is achieved by cutting down the number of interior point iterations (see Figs. 10 and 11). Additionally, the number of member adding iterations is reduced by 1 when the warm-start strategy is used for the problem with τ = 10. In both cases, the warmstart is used starting in the fourth member adding iteration. The optimal designs are shown in Fig. 12, where Fig. 12b shows the solution without stability constraints, which has a volume of 0.05122 m 3 and λ min = 8.0376e−08 (i.e. is c Optimal designs with stability constraints. Note that in the plots, the member sizes are scaled for sake of visual clarity clearly not stable). Figure 12c and d show stabilized designs obtained, respectively, when τ = 1 and τ = 10. When τ = 1 the stable design has a volume 0.05147 m 3 with violation of the kinematic compatibility constraint equal to 5.2354e−06. When τ = 10 the stable design has a volume 0.05376 m 3 with a slightly larger violation of the kinematic compatibility constraint, equal to 5.3589e−04.
Remark 10 The small-scale bridge problem considered in Section 6.2.1 and the large-scale bridge problem considered in Section 6.3.1 demonstrate that when problem size increases, the fraction of the bars used in the final SDPs decreases. For the small-scale bridge problem with 3200 bars, 18.5% of the bars were needed to find the solution. However, for the large-scale bridge problem with 90,100 bars, only 5.6% of the bars were needed when τ = 1, and 8.1% of the bars when τ = 10. This indicates that the adaptive member adding procedure is likely to bring more and more benefit as problem size increases.

Stadium-roof application with multiple load cases
We solve the stadium roof design problem shown in Fig. 13a. The roof is subject to three load cases: LC1 = f 1 , LC2 = f 1 + f 2 , and LC3 = f 1 + f 3 , where the loads f 1 = 0.27 kN/m 2 , f 2 = 2.7 kN/m 2 , and f 3 = 0.75 kN/m 2 are uniformly distributed. Note that the loads and dimensions have been simplified in the interests of clarity. The roof spans 40 m in the y-direction, 80 m in the x-direction and 4.2 m in the the z-direction. Detailed dimensions are given in the caption of Fig. 13a. The layout optimization problem has 36, 856 potential members. We first solve the problem without stability constraints, obtaining the design shown in Fig. 13b, comprising parallel disconnected planar trusses with a volume 2.2992 m 3 and with minimum positive eigenvalues for the three load cases LC1, LC2, and LC3 of 7.0387e−04, 5.6185e−05, and 1.8554e−04, respectively. This indicates that the structure is not stable for all load cases. Note that since this is a multiple load case problem, we expect large violations of the kinematic compatibility constraint, even for problems without the stability constraints, as mentioned in Remark 4. In this case, the violation was 0.0011 even when τ = 0. Next, we set τ = 10, = 1, 2, 3, and solve problem (9) with stability constraints to obtain the design shown in Fig. 13c, where the parallel planar trusses are now connected. In this case, the volume of the structure is 2.3279 m 3 , only slightly higher than before. Moreover, the violation of the kinematic compatibility equation (10) by the stable design was found to be 0.3190, which is large compared with the single-load case examples considered previously. Computational details are reported in Table 5 and Fig. 14. The CPU time without the warm-start strategy was 3194 s, and this is improved by 26% when the warm-start strategy is used, for which the CPU time was 2365 s. Once again, this is achieved by cutting down the number of interior point iterations during warm-starting. In both cases, a total of 6 member adding iterations were required to obtain the solution and the final SDPs had approximately 7% of the entire 36,856 bars. c Problem sizes, with/without warm-start strategies. The warm-start was used from the 3rd member adding iteration We have solved the truss layout optimization problem with global stability constraints via linear semidefinite programming by relaxing the nonlinear kinematic compatibility constraint. A primal-dual interior point method has been used, tailored to solving these problems efficiently. The implementation utilizes the sparse structure and low-rank property of the element stiffness matrices to reduce computational complexity when determining the linear systems arising in the algorithm. Moreover, we have extended the range of application of the adaptive member adding and warm-start techniques previously applied to truss layout optimization problems formulated as linear programs, so these can now be applied to problems modelled as semidefinite programs. By doing so, we have been able to find solutions to large-scale problems that could not have been solved using previously available methods and standard desktop computers.
We have demonstrated the validity of the solutions obtained for the relaxed problem by comparing them with solutions obtained for the original nonlinear problem for values of the stability load factor that are of interest to engineering practitioners.
Finally, direct methods were used to solve the linear systems arising from the interior point algorithm. The computational effort might be further reduced by use of iterative methods.

Replication of results
The input and output data that are used for all of the examples described in Section 6 are explicitly provided there. The same material has been used in all examples and the material properties are reported directly before Section 6.1. Note that these values and the applied loads require appropriate scaling if one wishes to use standard SDP solvers.