MILP‑based discrete sizing and topology optimization of truss structures: new formulation and benchmarking

Discrete sizing and topology optimization of truss structures subject to stress and displacement constraints has been formulated as a Mixed-Integer Linear Programming (MILP) problem. The computation time to solve a MILP problem to global optimality via a branch-and-cut solver highly depends on the problem size, the choice of design variables, and the quality of optimization constraint formulations. This paper presents a new formulation for discrete sizing and topology optimization of truss structures, which is benchmarked against two well-known existing formulations. Benchmarking is carried out through case studies to evaluate the influence of the number of structural members, candidate cross sections, load cases, and design constraints (e.g., stress and displacement limits) on computational performance. Results show that one of the existing formulations performs significantly worse than all other formulations. In most cases, the new formulation proposed in this work performs best to obtain near-optimal solutions and verify global optimality in the shortest computation time.


Previous work
Optimization has been successfully employed to improve the performance of a structure subject to a given set of loads and boundary conditions. Computing power as well as the availability of optimization solvers have increased significantly over the last decades, which have enabled structural optimization to become increasingly adopted in practice. Maximization of material usage efficiency through minimization of the structure weight is a common objective function. In this context, structural optimization may also contribute to reducing environmental impacts (Lagaros 2018).
Optimal layouts can be obtained through optimization of the structural geometry, topology, and member crosssection sizing (Rozvany et al. 1995). Typical methods to optimize truss systems are based on the ground structure approach introduced by Dorn et al. (1964). A ground structure consists of a network of members of which a subset will be selected for the optimal design. Typically, structural optimization is carried out using continuous design variables for the member cross-section areas (Gilbert and Tyas 2003;Ohsaki 2017;Rozvany et al. 1995). The adoption of continuous variables typically simplifies the problem and allows the use of efficient numerical techniques (e.g., linear and non-linear programming). Digital fabrication techniques such as additive manufacturing enable the custom design of structures that are made of non-standard components (Pellens et al. 2019;Smith et al. 2016). However, conventional building construction is typically restricted to standard sets of discrete cross sections (e.g., I-beams, hollow sections, etc.), thus, making the design variables discrete, which has led to the development of discrete sizing and topology optimization methods (Huang and Arora 1997;Toakley 1968). Stolpe (2016) provides an extensive review of structural Meta-heuristic and stochastic methods including evolutionary and genetic algorithms (Deb and Gulati 2001;Kaveh and Kalatjari 2002;Rajan 1995;Rajeev and Krishnamoorthy 1992), simulated annealing (Bennage and Dhingra 1995;Kripka 2004;Shea and Smith 2006), ant-colony optimization (Camp et al. 2005), and particle swarm optimization (Li et al. 2009) have been applied. Typical disadvantages of meta-heuristic algorithms are a large number of structural analyses necessary to achieve convergence and that solution optimality cannot be guaranteed (Stolpe 2016). A different approach is to formulate discrete sizing and topology optimization as a Mixed-Integer Linear Programming (MILP) problem, which can be solved to global optimality through combinatorial optimization techniques such as branch-and-bound methods (Bertsimas and Tsitsiklis 1997;Nemhauser and Wolsey 1999;Stolpe and Svanberg 2003). Mathematical formulations for discrete sizing and topology optimization of trusses based on MILP have been presented by Ghattas and Grossmann (1991) and Rasmussen and Stolpe (2008) among others. Mela (2014) extended the formulation of Rasmussen and Stolpe (2008) by adding decision variables and logic constraints to avoid mechanisms. Shahabsafa et al. (2018) presented a formulation for cross-section size optimization that combines MILP with neighborhood search to reduce computation time. MILPbased methods guarantee solution global optimality, which is a significant benefit as it allows for rigorous benchmarking between solutions that are subject to different boundary conditions. Recent work shows that MILP formulations have been employed more frequently for the optimal design of reticular structures. This is because MILP-based truss optimization enables ease of integration of "real-world" constraints into the problem formulation. For example, member buckling based on code regulations (Mela 2014), constraints on the number of structural nodes and members (Fairclough and Gilbert 2020), construction-related constraints such as avoiding member overlapping and intersection (Fairclough and Gilbert 2020;Ohsaki and Katoh 2005), and joint capacity constraints computed based on code regulations (Van Mellaert et al. 2016). Discrete sizing and topology optimization has been extended to design truss and frame structures from a limited stock of reclaimed elements (Brütting et al. 2020(Brütting et al. , 2019. In these studies, the objective is the minimization of environmental impact quantified through Life-Cycle Assessment. Generally, previous work has shown that discrete structural optimization based on MILP is relevant to research and practice, and thus, strengthens the need for new formulations with improved computational efficiency. MILP problems have been successfully solved using branch-and-bound algorithms whose computational performance is significantly influenced by the problem formulation, i.e., by the chosen set of design variables and optimization constraints. Quoting Bertsimas and Tsitsiklis (1997, p. 461), "In linear programming, a good formulation is one that has a small number […] of variables and constraints […] because the computational complexity of the problem grows polynomially […]. In addition, given the availability of several efficient algorithms for linear programming, the choice of a formulation, although important, does not critically affect our ability to solve the problem. The situation in integer programming is drastically different. Extensive computational experience suggests that the choice of a formulation is crucial." In MILP-based formulations, the computation time to obtain the global optimum does not necessarily increase with the increase in the number of variables and constraints (Bertsimas and Tsitsiklis 1997).

New contributions of this work
The computation time to solve a MILP problem to global optimality depends on the problem size and problem formulation (Bertsimas and Tsitsiklis 1997;Nemhauser and Wolsey 1999), i.e., on the number of variables, the number of constraints, and the quality of the constraint formulations. Different parameters, such as the number of structural members, number of candidate cross sections, number of load cases, as well as the formulation of stress and displacement constraints significantly influence the computation time to obtain globally optimal solutions. No comprehensive study exists in the literature that evaluates the influence of such parameters on the computational performance of different MILP formulations. New contributions offered by this work are 1. Analysis of two well-known discrete sizing and topology optimization formulations (MILP-based). 2. The introduction of a new MILP formulation for sizing and topology optimization of truss structures. 3. A comparison of the problem size for all presented MILP formulations. 4. A computational performance benchmark by applying the MILP formulations to case studies. 5. Recommendations to choose a specific formulation depending on design constraints.
The benchmark provides a qualitative and quantitative evaluation of the influence of formulation specificities on computational performance. Results show that one of the existing formulations performs significantly worse compared to the other formulations. In most cases, the new formulation proposed in this work performs best to obtain near-optimal solutions and verify global optimality in the shortest computation time.

Outline
This paper is organized as follows. Section 2 starts with a review of existing MILP formulations for discrete truss sizing and topology optimization subject to ultimate and serviceability limit states (stress constraints, nodal displacement limits). Sections 2.2 and 2.3 give the problem statement of two well-known formulations, (F RS ) by Rasmussen and Stolpe (2008) and (F GG ) by Ghattas and Grossmann (1991). Drawing from (F RS ) and (F GG ), a new MILP formulation is proposed denoted as (F BSF ).
The new formulation comes into two variants depending on the set of employed state variables: ( F BSF 1 ) uses member forces and elongations (Sect. 2.4.1); ( F BSF 2 ) employs member elongations only (Sect. 2.4.2). Section 2.5 compares the problem size, i.e., the number of design variables and constraints, for all MILP formulations. A brief overview of the branch-and-bound procedure to solve MILP problems to global optimality as well as of the solver adopted in this work are given in Sect. 2.6. In Sect. 3, the MILP formulations are applied to the design of planar (2D) and spatial (3D) structural configurations and benchmarks related to computational performance are provided. Section 4 discusses results and gives recommendations for choosing a specific formulation depending on design constraints. Section 5 concludes this paper.

General formulation
This section presents a basic formulation for discrete crosssection sizing and topology optimization of truss structures to minimize the structure volume. It is assumed, for simplicity that all structural members have the same material, and thus, the objective function is equivalent to the minimization of the structure weight. The formulation considers a truss ground structure with a total of m members. Each truss member is denoted with index i. The catalog (or set) of discrete candidate cross sections contains n elements. Reference to a distinct cross section within the set is given by index j. Assignment of a cross section j at a member position i is expressed through a binary assignment variable t ij ∈ {0,1}: Cross-section assignment is combined with the evaluation of structural constraints following the simultaneous 1 if cross-section j is assigned to member position i 0 if cross-section j is not assigned to member position i analysis and design (SAND) approach (Arora and Wang 2005;Haftka 1985), whereby assignment variables and state variables (e.g., member forces and nodal displacements) are treated as optimization variables. SAND is different from nested analysis and design (NAND) because it avoids the need to perform explicit structural analysis (i.e., matrix inversion) and design sensitivity analysis at each optimization iteration (Arora and Wang 2005). The SAND approach is a key to the formulation of discrete sizing and topology optimization as a MILP problem and to solving it to global optimality (Rasmussen and Stolpe 2008). Structural constraints include the equilibrium of forces at nodes, geometric compatibility of member elongations with nodal displacements as well as stress and displacement limits. All methods given in this work are formulated with the assumption of linear-elastic material behavior, small strains, and small displacements, i.e., no second-order effects such as geometric and material nonlinearity are considered. Since the formulations by Rasmussen and Stolpe (2008) and Ghattas and Grossmann (1991) do not consider member buckling constraints, for ease of comparison, the formulations given in this work do the same. Member buckling is a linear constraint in MILP formulations (Mela 2014) and, hence, could be added to the formulations if required. Similarly, no constraint is added to avoid the emergence of mechanisms through topology optimization. For simplicity, the formulations are presented for one load case. Extension to multiple load cases is straightforward by adding state variables and structural analysis constraints for each load case.
The basic problem formulation for minimum-volume discrete sizing and topology optimization is denoted with (F 0 ): The design variables of (F 0 ) are the binary assignment variables t ij collated in the vector t ∈ {0,1} mn . The column vector t is assembled as t = [[t 11 , …, t 1j ,…, t 1n ], …, [t i1 , t i2 ,…, t in ],…, [t m1 , t m2 ,…, t mn ]] T . The continuous state variables are member internal forces p ∈ ℝ m and nodal displacements u ∈ ℝ d , where d denotes the number of free degrees of freedom (i.e., nodal displacements). The optimization objective (Eq. (2)) is the minimization of the structure volume obtained from the product of member length l i and assigned cross-section area a j . Equation (3) ensures that at most one cross-section j is assigned to member position i. The structure topology changes if no cross section is assigned to a member position. Static equilibrium of forces at nodes is ensured by Eq. (4), where B ∈ ℝ d×m is the reduced equilibrium matrix and f ∈ ℝ d is the static external load. B concatenates the direction-cosine vectors b i ∈ ℝ d of all m members. Equation (5) combines geometric compatibility with constitutive laws by relating nodal displacements with member forces through element elongations. In Eq. (5), e i is the Young's modulus for member i. Stress constraints in Eq. (6) limit the member forces to the admissible member stress σ i min in compression (negative number) and σ i max in tension (positive number). The displacements u are bounded by lower and upper limits u min and u max (Eq. (7)) as, for example, required by serviceability limits.
In (F 0 ), all equations, except the compatibility constraint (Eq. (5)), are linear functions of the optimization variables t, p, and u. The compatibility constraint contains a product of binary assignment variables t ij and continuous displacement variables u, which is a bi-linear term and, therefore, (F 0 ) is not in the standard MILP problem form. To obtain a standard MILP form, the bi-linear constraints must be reformulated into linear constraints through "big-M" techniques and the introduction of auxiliary variables (Glover 1984(Glover , 1975). The following sections present different ways to reformulate (F 0 ) into an equivalent MILP in standard form. The first two are the well-known formulations (F RS ) and (F GG ) by Rasmussen and Stolpe (2008) and Ghattas and Grossmann (1991), respectively. Although (F GG ) was presented 17 years before (F RS ), for methodological reasons, they are here discussed in anti-chronological order. Thereafter, the new formulation offered by this work is introduced.

Formulation ( F RS ) by Rasmussen and Stolpe
To linearize (F 0 ), Rasmussen and Stolpe (2008) extend the member force vector p ∈ ℝ m to p' ∈ ℝ mn which holds one continuous variable p ′ ij for each combination of member i and cross section j, i.e., for each assignment variable t ij as assembled in the vector t. The vector p' is referred to as the extended member force vector. Accordingly, the compatibility constraint (5) is adapted to Note that Eq. (8) still contains a bi-linear product of t ij and u. As shown by Rasmussen and Stolpe (2008), the following condition must hold: This condition can be formulated through the following inequality constraints that do not contain bi-linear products: Equations (10) and (11) are referred to as big-M formulations, where c min ij and c max ij are the big-M constants. From Eq. (10), it follows that p ′ ij is set to zero when t ij = 0, and that p ′ ij is bounded between c min ij and c max ij when t ij = 1. At the same time, when t ij = 1, the constraint in Eq. (11) becomes b T i u e i a j l i = p � ij as required by Eq. (9). A sufficiently large number should be set for the big-M constants to avoid truncating the solution domain and excluding feasible solutions. Rasmussen and Stolpe (2008) show that when the displacements u are bounded (Eq. (7)), the constants c min ij and c max ij , which are forces in this context, can be computed from constitutive laws based on the maximum admissible nodal displacements: The right-side terms in Eqs. (12) and (13) are linear programming problems with the nodal displacements u as design variables. These linear programming problems can be solved analytically due to the simple bounds on the admissible displacement limits (u min ≤ u ≤ u max ) (Rasmussen and Stolpe 2008): In Eqs. (14) and (15), b ir denotes the rth entry in the column vector b i of the equilibrium matrix B ∈ ℝ d×m . u min r and u max r denote the respective entries in the nodal displacement bounds u min and u max . The notation "|" expresses the condition to consider only those entries where b ir is larger or smaller than zero, respectively. Note that Eqs. (14) and (15) are precomputed for each assignment t ij before optimization (Eqs. (10) and (11)). The so-obtained big-M constants c min ij and c max ij are then employed in the constraints of the global optimization problem. The MILP formulation (F RS ) as given by Rasmussen and Stolpe (2008) is: Since the force variables are extended from p ∈ ℝ m to p' ∈ ℝ mn , the equilibrium matrix B ∈ ℝ d×m must be equivalently extended to B' ∈ ℝ d×mn by horizontally concatenating each direction-cosine vector b i for n times (Rasmussen and Stolpe 2008). Note that Rasmussen and Stolpe (2008) do not explicitly include the big-M Eq. (10) into (F RS ) since the admissible stress constraint (Eq. (20)) serves the same purpose of setting p ij to zero when t ij = 0.

Formulation ( F GG ) by Ghattas and Grossmann
An alternative reformulation of (F 0 ) is given by Ghattas and Grossmann (1991). The compatibility constraints in Eq. (5) between member elongations and forces are separated into two equality constraints: (1) geometric compatibility of nodal displacements u with member elongation i , and (2) constitutive relation between member elongation i and internal force p i : Equation (23) contains a bi-linear product between t ij and i . Ghattas and Grossmann (1991) replace each member elongation i with auxiliary continuous variables v ij for each combination of member i and cross section j, i.e., for each assignment variable t ij . The vector v is denoted as the extended member elongation vector. Accordingly, the constitutive relations are stated as follows: In addition, the following big-M condition must hold: Equation (25)  The complete MILP formulation (F GG ) as presented by Ghattas and Grossmann (1991) is: The vector v ∈ ℝ mn collates each member elongation variable v ij in the same order as the assignment variables are assembled in t (Sect. 2.2). Note that (F GG ) includes member stresses σ ∈ ℝ m as continuous design variables (one variable per member i) in Eq. (33) which is employed to relate the elongation v ij and stress i variables through Hooke's law. The member stress variables σ are bounded by admissible stresses min and max in compression and tension. The vectors min and max collate the admissible member stresses also employed in the big-M Eqs. (26) and (27).
Different to (F RS ), in (F GG ) cross-section assignment Eq. (29) is an equality constraint, which ensures that exactly one element is assigned to member position i. This avoids that the elongations v ij of a member i are set to zero through Eq. (34) when the assignment variables t ij are zero. If this was allowed, from geometric compatibility (Eq. (31)), it would follow that b T i u = 0, which would set the displacements of the nodes connected to member i also to zero or to a displacement that does not result in member elongation. However, this would violate compatibility conditions for other members that connect to the same nodes.
Hence, since assignment Eq. (29) is an equality constraint, a zero-area element with a j = 0 must be added to the set of candidate cross sections to enable topology optimization (Ghattas and Grossmann 1991). When a zero-area element is assigned to member i, which in practice means member i is removed, Eq. (32) automatically sets the member force p i to zero to satisfy the equilibrium constraint. However, after member i is removed, the related elongation v ij remains a design variable that must fulfill geometric compatibility (Eq. (31)) and the Big-M constraints (Eq. (34)). Therefore, for a zero-area element, the big-M constants must be set sufficiently large not to exclude feasible solutions. Although not explicitly mentioned by Ghattas and Grossmann (1991), the following conditions are considered in this work to enable topology optimization: When a j ≠ 0, the Big-M constants min ij and max ij are computed from the admissible stresses in compression ( min i ) and tension ( max i ), respectively. When a j = 0, i.e., the member is removed, the elongation variable v ij should be unbounded and, therefore, min ij and max ij are set to − ∞ and + ∞, respectively. For the MILP solver employed in this work, a numerical value of >10 30 is treated as infinity. That being said, even when unbounded, v ij will not take an infinite value due to the geometric compatibility constraint that relates it to the nodal displacements, which themselves are bounded. In addition, only the elongations v ij of members with a non-zero area (a j ≠ 0) are considered in Eq. (33); hence, the constraint is not active for removed members.

New formulation
Drawing from the formulations described in Sects. 2.2 and 2.3, a new MILP formulation is proposed. The new formulation comes in two variants depending on the set of employed state variables: ( F BSF 1 ) uses member forces and elongations (Sect. 2.4.1) and ( F BSF 2 ) member elongations only (Sect. 2.4.2).

Formulation ( F BSF 1 ) with member force and elongation variables
This section introduces the first step of a new discrete structural optimization formulation ( F BSF 1 ) that draws from (F RS ) and (F GG ) but differs in the choice of design variables and constraints. The new formulation ( F BSF 1 ) adopts a similar approach for the formulation of geometric compatibility constraints. However, the redundant stress variables adopted in (F GG ) are excluded in ( F BSF 1 ). As mentioned by Bollapragada et al. (2001), to reduce (F GG ) problem size, member stress variables can be substituted with elongations as combines the big-M formulations from (F RS ) and (F GG ) by bounding the member elongation variables through stress constraints, as done in (F GG ), as well as through displacement constraints adapting (F RS ) formulation. ( F BSF 1 ) includes assignment variables t ∈ {0,1} mn as well as continuous state variables for member forces p ∈ ℝ m , member elongations v ∈ ℝ mn , and nodal displacements u ∈ ℝ d .
Similar to (F GG ), in ( F BSF 1 ) the extended member elongations v ij are treated as state variables. This way, the relation of member elongations with forces in Eq. (5) can be separated into geometric compatibility of nodal displacements u with member elongation v ij in Eq. (42), and constitutive relations between member elongation v ij and internal force p i in Eq. (43). The big-M constraints in Eq. (44) ensure that the member elongation v ij = 0 when t ij = 0, and otherwise, it is bounded by Δ min ij and Δ max ij when t ij = 1. The bounds Δ min ij and Δ max ij are computed from both material ( ) and geometric ( ) conditions: The lower bound Δ min ij for elongation v ij is set to the largest among the two values min ij and min ij . The upper bound Δ max ij is set to the smallest among of the two values max ij and max ij . Note that these bounds are exclusively based on the displacement limits u min and u max , and thus, they are independent of cross-section assignment (no index j). The bounds in Eqs. (51) and (52) are related to member elongations rather than to internal forces as it is done in (F RS ). The linear programming problems in Eqs. (51) and (52) can be solved analytically before optimization, as done for Eqs. (12) and (13). Which of the bounds or is tighter, usually depends on material strength and displacement limits, i.e., whether the problem is governed by stress or stiffness.

Formulation ( F BSF 2 ) with member elongation variables
In ( F BSF 2 ), the number of state variables is reduced further from ( F BSF 1 ) by expressing member forces as resulting from elongations: Since the force variables p ∈ ℝ m are replaced with member elongation variables v ∈ ℝ mn , the equilibrium matrix B ∈ ℝ d×m must be extended to B V ∈ ℝ d×mn by horizontally concatenating each direction-cosine vector b i for n times and multiplying each entry with the respective member stiffness e i a j /l i . Note that for zero-area members, the entry becomes zero (a j = 0), and hence, a member that is removed does not contribute to nodal equilibrium, as required. The complete formulation ( F BFS 2 ), using only member elongations and nodal displacements as continuous state variables, is: Geometric compatibility and big-M constraints (Eqs. (57) and (58)) are identical to those in formulation ( F BSF 1 ). Stress constraints that limit the member forces are not explicitly added in ( F BFS 2 ), but they are included via the big-M constraints (cf. Eqs. (47) and (48)).

Comparison of MILP problem sizes
A MILP problem has the following standard form: The column vector of design variables x = [x bin , x con ] T concatenates the binary variables x bin ∈ {0, 1} n var,bin (here representing the assignment variables t) and all continuous state variables x con ∈ ℝ n var,con (e.g., member forces p, elongations v, and displacements u). The size of x is n var = n var,bin + n var,con . The objective function value is the inner product of the cost vector c ∈ ℝ n var and design variables x. For the weight optimization formulations considered in this work, all costs associated with continuous variables are zero. All formulation constraints are collated in the constraint matrix A ∈ ℝ n con ×n var , where n con is the total number of constraints. The right-hand side vector is denoted as b ∈ ℝ n con . For notation simplicity, the inequality sign ' ≤ ' in Eq. (61) represents both equality and inequality constraints.
As an example, Fig. 1 shows schematically the assembly of the structural optimization constraints for ( F BSF 1 ) into the linear constraint matrix A and right-hand side vector b. The black lines within the gray blocks indicate the contribution of each constraint and, thus, the non-zero entries. The linear constraint matrix typically is highly sparse. Equilibrium and compatibility constraints are included in A, and thus, B is enclosed in A.
The MILP formulations described in Sects. 2.2 to 2.4 differ in the set of design variables and constraints. Table 1 gives information about the optimization problem size for each formulation, in particular the number of binary design variables n var,bin , the number of continuous design variables n var,con , the number of linear constraints n con , and the number of non-zeros n nz entries in the linear constraint matrix A. The number of non-zero entries in A gives an indication of the amount of data that needs to be processed by the MILP solver. The number of variables n var and constraints n con depends on the number of members m in the structure, the number of candidate cross sections n, and the number of free nodal displacements d. Formulation (F RS ) contains m·n binary assignment variables in the vector t, whereas formulations (F GG ), ( F BSF 1 ), and ( F BSF 2 ) contain m·(n + 1) binary assignment variables. The notation (n + 1) refers to adding a zero-area element to the catalog of n cross sections to allow element removal for topology optimization.
The number of members m and candidate cross sections n are the two most significant parameters. Further, the number of constraints n con and non-zeros n nz depend on the number of free nodal displacements d and the number of non-zeros n nz,B in the equilibrium matrix B. As shown in Fig. 1, B is enclosed in A and d is the number of equilibrium constraints since it is the number of free degrees of freedom (i.e., free nodal displacements). Typically, d and n nz,B depend on how many members are aligned with the coordinate axis (direction cosines) and how many members are connected at nodes (node valence). Referring to Table 1, the total number of variables is 2·m·n + d for all four formulations, plus an additional term 4·m, 3·m, and 2·m for (F GG ), ( F BSF 1 ), and ( F BSF 2 ), respectively. This indicates that the difference in the number of variables is not significant since the dominant term 2·m·n is identical for all formulations. Instead, the number of constraints increases much faster with m an n in the case of formulation (F RS ) compared to the other formulations. The dominant term for the number of constraints in (F RS ) is 4·m·n, which is double compared to the other formulations. Formulation (F RS ) has the highest number of non-zeros n nz since the m·n term is multiplied by the highest scalar and the n nz,B term scales with n while in (F GG ) and ( F BSF 1 ) the n nz,B term is constant. A visual aid of the problem size in relation to the number of members and candidate cross sections is given in Fig. 2. The plot has been obtained based on a study of ground structures with members connecting only neighboring nodes and expressing d and n nz,B as functions of m.

MILP solver
MILP problems can be solved to global optimality using combinatorial optimization such as branch-and-bound algorithms (Nemhauser and Wolsey 1999; Bertsimas and Tsitsiklis 1997). A branch-and-bound algorithm explores a search  44) and (45)). 1 and 0 denote vectors of ones and zeros, respectively tree while solving continuous linear programming (LP) relaxations of the original (mixed-) integer programming problem (Bertsimas and Tsitsiklis 1997). Since all variables are treated as continuous, solving the LP relaxation gives a lower bound f LB of the MILP objective function value. At the same time, when a solution of an LP relaxation is integer feasible, i.e., all binary variables take a value of 0 or 1, it is denoted as the upper bound f UB . Once the first upper bound f UB is obtained, the so-called optimality-or "MIP-gap" is computed as follows: The MIP-gap provides a quality measure for the integerfeasible solutions obtained during the branching process. At any step of the process, the expected worst-case difference between the global optimum and the best integer-feasible solution obtained (f UB ) is known. The availability of information on solution quality is a clear advantage compared to other optimization methods such as meta-heuristic algorithms. When the gap is zero (f UB = f LB ), the solution global optimality is verified. In this work, all MILP problems are solved through Gurobi (Gurobi Optimization, LLC 2019). Gurobi combines a branch-and-bound solver with cutting plane methods into a so-called "branch-and-cut" algorithm (Bertsimas and Tsitsiklis 1997; Gurobi Optimization, LLC 2019). Adding cutting planes, i.e., linear inequality constraints, is automatically carried out by the solver in order to reduce the size of the feasible domain (i.e., the size of the search space) of the problem. Cutting planes decrease the feasible domain size of the relaxed LP problem but do not change the feasible domain of the integer problem. This way, the bounds obtained from LP relaxations can be tightened which, generally, reduces the size of the search tree. In this work, Gurobi 9.1.2 has been used with default settings except for IntFeasTol (integrality) and FeasibilityTol (solution feasibility) which have been reduced from 10 -5 to 10 -6 and from 10 -6 to 10 -7 , respectively. This change has been made to ensure solution feasibility for all case studies. The optimization problems have been modeled via the

Benchmark studies
The four formulations (F RS ), (F GG ), ( F BSF 1 ) and ( F BSF 2 ) presented in Sect. 2 differ in the set of variables and constraints, which influences the convergence speed of the branch-andbound solver. In this section, the computational performance of the four formulations is benchmarked through case studies.
To investigate the influence of the Big-M constraints Eqs. (47) to (52), formulation variants are also considered. In particular, the elongation bounds based on displacement constraints (Eqs. (51) and (52)) are omitted for ( F BSF 1 ) and ( F BSF 2 ) and they are added to (F GG ). These three variants are denoted with (F GG )*, ( F BSF 1 )*, and ( F BSF 2 )*. Note that the displacement-based big-M constants min i and max i affect only zero-area elements (a j = 0). For all other cross-section sizes, the admissible stress/strain relation (cf. min ij and max ij ) governs (Eqs. (47) and (48)) as displacement bounds are set to a relatively larger value in the following case studies. For ( F BSF 1 )* and ( F BSF 2 )*, this means that the elongation variables v ij are unbounded (± ∞) when a j = 0, whereas they are bounded by min

Case study 1: L-truss
The first case study, which is taken from Rasmussen and Stolpe (2008), is an L-truss ground structure consisting of m = 54 members. Dimensions, support locations, and loading are indicated in Fig. 3. The structural elements material is assumed to be aluminum with Young's modulus of 70,000 MPa and yield strength of 170 MPa. The set of candidate cross-section areas is a = {5, 10} × 10 -3 m 2 . The minimum and maximum admissible displacements of all unsupported nodes are set to ± 2.0 m in both the X-and Y-direction. Two cases are investigated: Sizing and topology optimization are carried out simultaneously. Although the optimized configurations are in equilibrium with the applied loads, they might not be kinematically stable (Sect. 2.1).  Figure 4 shows the optimal topology and sizing for both cases, (a) and (b). In Fig. 4, the cross-section area is indicated by line thickness and numeric labels. Table 2 provides optimization metrics on problem size and branch-and-bound process for the four formulations. The solution produced by all formulations has a total volume of 0.0466 m 3 in case (a) and 0.0572 m 3 in case (b), which are identical to the values reported by Rasmussen and Stolpe (2008). Generally, including two load cases in (b) results in a larger problem size than in (a) ( Table 2). The number of explored nodes in the search tree corresponds to the number of LP relaxations that are solved during the branch-and-bound process. The fewest nodes have been explored by ( F BSF 1 ). In case (b), the number of simplex iterations that have been performed to solve the continuous LP relaxations is much lower for (F GG ), ( F BSF 1 ), and ( F BSF 2 ) and their variants (*) compared to (F RS ). The computation time to obtain the integer-feasible solution (upper bound), which is denoted as "Optimum obtained after" in Table 2, is less than 1 s for all formulations in case (a). The remaining computation time is required to prove global optimality of the upper-bound solution by obtaining increasing lower bounds until the MIP-gap = 0%. Figures 5 and 6 show the convergence plot for all formulations in cases (a) and (b), respectively. For all formulations, the objective value of the continuous root relaxation (the very first lower bound solution) is 0.02779 m 3 and 0.0287 m 3 in cases (a) and (b), respectively, and it is obtained in less than 10 ms. Generally, the optimal solution (upper bound, continuous lines) is obtained relatively early in the branch-an-bound process. The lower bounds (dashed lines) increase from the root relaxation to the global optimum. Global optimality is verified when lower and upper bounds intersect. For (F RS ), after approximately 1/4th of the total time, the lower bounds stagnate and increase very slowly. (F RS ) requires approximately 3 to 12 times longer computation time to prove global optimality  than the fastest solution, which is obtained using ( F BSF 2 ) for case (a) and ( F BSF 2 )* for case (b). In both cases, adding displacement-based elongations in (F GG )* (dark green) increases the computation time to prove global optimality compared to (F GG ) (light green). Instead, in both cases, the upper bound is obtained in the shortest computation time by (F GG )*. The removal of displacement-based elongation bounds (big-M constants), helps improve the performance of ( F BSF 2 )* only in case (b). Instead, removing displacement-based elongation bounds increases the computation time to prove global optimality using ( F BSF 1 )* in case (b).

Case study 2: 3D-cantilever
The second case study, which is taken from Rasmussen and Stolpe (2008), is a 3D-cantilever truss consisting of m = 40 members. Dimensions, support locations and loading are given in Fig. 7  • Case (b): the set of 20 candidate cross-section areas is a b = {0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10} × 10 -3 m 2 Sizing and topology optimization are carried out simultaneously. Figure 8 shows the optimal topology and sizing for both cases (a) and (b). Table 3 provides optimization metrics on the problem size and branch-and-bound process for all four formulations. The solution produced by all four formulations has a volume of 0.656 m 3 in case (a) and 0.546 m 3 in case (b). The solution for (a) is identical to the one reported by Rasmussen and Stolpe (2008). Case (b) has not been investigated in previous work. As for Case Study 1 (Sect. 3.1), formulation variants (F GG )*, ( F BSF 1 )*, and ( F BSF 2 )* are tested in addition to the original formulations. As before, due to the large displacement bounds (± 2.0 m), the variation only affects the bounds for the elongation variables related to the zero-area elements.
Generally, including significantly more candidate cross sections in case (b) results in a larger problem size than in case (a), and hence, requires longer computation time (Table 3). Figures 9 and 10 show the convergence plots for all formulations and cases. In cases (a), (F RS ) requires the longest computation time to verify global optimality. In case (b), (F RS ) performs significantly worse compared to all other formulations. (F GG )* verifies global optimality in the shortest computation time in both cases followed closely by ( F BSF 2 )*.
• Case (e): Topology and sizing optimization. Displacement limit set to ± 2 in. • Case (f): Sizing optimization. Displacement limit set to ± 2 in.
These cases are selected to study the effect of topology optimization and displacement bounds on the computation time. As suggested by Bollapragada et al. (2001) and Van Mellaert (2017), discrete truss optimization problems may become significantly harder to solve in the presence of tight displacement constraints. Formulation variants (F GG )*, ( F BSF 1 )*, and ( F BSF 2 )* (cf. Sect. 3.1) are only considered for cases (a), (b), and (c) where topology optimization is included. In the other cases (b), (d), and (f), no zero-area elements can be assigned and due to stress limits (i.e., strains), the elongation variables v ij are always constrained by stressbased bounds rather than by displacement-based bounds (cf. Eqs. (47) and (48)). Figure 12 shows the optimal topology and sizing for all cases (a)-(f). For cases (a) and (b), all formulations produce identical, globally optimal designs. In case (b), the solution has a weight of 1856.7 lb which is identical is different than that obtained by the other formulations, as shown in Fig. 12(c) and (d).
Generally, including topology optimization only marginally increases the problem size: only a zero-area element is added to the other 42 candidate cross sections (cf. cases (a) and (b), Table 4). Figure 13 , 14, 15, 16, 17, and 18 show the convergence plots for cases (a) -(f), respectively. In cases (a) -(d), ( F BSF 1 ) (*) is the fastest to verify global optimality. The lower bounds obtained through ( F BSF 1 )* increase significantly faster than through the other formulations, thus, allowing the verification of global optimality in a shorter computation time. The lower bounds of (F RS ) practically stagnate in (c)  displacement bounds are tightened, none of the formulations can verify global optimality within the set time limit of 5 h. The lower bounds increase steeply early in the branchand-bound process and stagnate after 5000 s. The remaining MIP-gap (distance between upper and lower bounds) at termination is the smallest for ( F BSF 1 )* and ( F BSF 1 ) in cases (e) and (f), respectively. For cases (a) and (b) with large displacement bounds, topology optimization has practically no effect on computation time. Allowing topology optimization in case (c) produces solutions faster than in case (d) in which the topology is fixed. The problem size does not change when the displacement limits u min and u max are tighter. Only the elongation bounds and big-M constants change.
Hence, it is clear that the computational performance is not only affected by the MILP problem size but also by numerical specificities. Generally, the formulations that employ member elongations as design (state) variables (F GG , F BSF 1 , and F BSF 2 ) are able to obtain better feasible solutions (upper bounds) in less time than (F RS ). This might be due to the direct correlation of member elongations and nodal displacements through geometric compatibility. The upper-bound solutions of 4962.1 and 5490.7 lb that have been obtained within a few seconds through (F GG ), ( F BSF 1 ), and ( F BSF 2 ) in cases (e) and (f) have the same quality as the best solutions reported in previous work in which solutions have been obtained through (meta-)heuristic algorithms (Camp et al. 2005;Juang et al. 2003;Kripka 2004;Rajan 1995;Schmid 1993;Sonmez 2011;Zhang et al. 2014). This shows that the MILP formulations are effective to produce high-quality solutions in a short computation time, even though global optimality might not be verified.

Case study 4: 25-bar tower
The fourth case study is a "25-bar tower" that has been often studied in previous work (Rajeev and Krishnamoorthy 1992;Zhang et al. 2014). Dimensions, support locations, and loading are indicated in Fig. 19. The tower members are divided into eight groups as indicated in Table 5. The labels "S-E" refer to the start and end node of each member. Node numbering is shown in Fig. 19. The structural elements' material is assumed to have Young's modulus of 10 7 psi and yield strength of 40 × 10 3 psi. Weight minimization is carried out by factoring the structure volume (Eq. (2)) with a material density of 0.1 lb/in 3 . The set of 30 candidate cross-section areas comprises a = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.8, 3.0, 3.2, 3.4} in 2 (Rajeev and Krishnamoorthy 1992). For all members within a group, the same cross section is assigned. Two cases are investigated: • Case (a): the X-and Y-displacements of Nodes 1 and 2 are limited to ± 1.0 in, the Z-displacements to ± 50.0 in • Case (b): the X-and Y-displacements of Nodes 1 and 2 are limited to ± 0.35 in, the Z-displacements to ± 50.0 in For all other nodes, the X-, Y-, and Z-displacements are limited to ± 50.0 in. The structure topology is kept constant and only cross-section sizing is performed. These cases are considered to study the effect of reducing the displacement bounds on computational performance for a more complex configuration (3D) compared to previous cases. Formulation variants (F GG )*, ( F BSF 1 )*, and ( F BSF 2 )* (cf. Sect. 3.1) are not considered for this case study because the topology is invariant and the elongation variables v ij are constrained by   (47) and (48)). Table 6 gives the optimization problem size for all four formulations. The problem size is identical for both cases (a) and (b) as only the displacement bounds are changed. Table 7 gives optimization metrics for the branch-and-bound process and Table 8 indicates the optimal cross sections assigned to each member group. Figures 20 and 21 show the convergence plots of the four formulations for cases (a) and (b), respectively. As in previous sections, formulations (F GG ), ( F BSF 1 ), and ( F BSF 2 ) obtain faster increasing lower bounds. In case (a), (F RS ) is significantly slower compared to the other formulations to verify solution global optimality. In case (b), the global optimum with a weight of 484.3 lb. has been obtained by all formulations. ( F BSF 2 ) obtains the optimal solution and verifies global optimality in the shortest computation time. Remarkably, ( F BSF 2 ) reaches the optimal upper bound after 219 s whereas it takes the other formulations 1 -5 h. In case (b), global optimality cannot be verified through (F RS ) within the set time limit of 5 h (remaining MIP-gap = 58%). The globally optimal design in (b) is identical to the best solution obtained through meta-heuristics in the literature (Zhang et al. 2014). To the authors' knowledge, it is the first time that solution global optimality for this 25-bar tower truss design is verified through a deterministic method (MILP).

MILP formulation performance
Different MILP formulations (F RS ), (F GG ), ( F BSF 1 ), and ( F BSF 2 ) for discrete sizing and topology optimization of *Premature termination after the computation time limit of 5 h is exceeded 3. The benchmark studies have shown that, generally, a slightly larger number of variables, constraints, and non-zeros entries in the constraint matrix may not necessarily result in a significantly longer computation time. This is the case for formulations (F GG ), ( F BSF 1 ), and ( F BSF 2 ). Differently, formulation (F RS ) involves a significantly larger number of variables, constraints, and non-zero entries in the constraint matrix (Sect. 2.5) which requires exploring more nodes of the search tree and thus solving more LP relaxations. In addition, the LP relaxations are harder to solve because they require more Simplex iterations due to the higher number of variables and constraints. For these reasons, (F RS ) requires a significantly longer computation time to reach the optimal solution (upper bound) and verify global optimality compared to the other three formulations.
It is more challenging to assess how the mathematical structure of a formulation affects computational performance. In the attempt to reach a holistic understanding of 3-10 6-7 4-9 5-8 3-8 4-7 6-9 5-10 3-7 4-8 5-9 6-10   the performance of the MILP formulations, a synthesis of the case studies discussed in (Sects. 3.1 -3.4) is offered in this section by focusing on three benchmark metrics m k∈{1,2,3} : m 1 total computation time to verify global optimality, m 2 the number of nodes explored in the search tree, and m 3 the computation time required to obtain upperbound solutions that are global optima (although global optimality has not been verified). The formulations are evaluated by ranking their performance for all cases and through performance profiles (Dolan and Moré 2002). Performance profiles have been employed to benchmark the performance of optimization solvers for topology optimization problems by Rojas-Labanda and Stolpe (2015) among others.
In this work, performance profiles are employed to evaluate the performance of different MILP formulations to solve volume minimization through discrete sizing and topology optimization. Since the problem is identical and all case studies are carried out using the same optimization solver, the focus of the assessment is on the performance of the mathematical formulation of the problem.  The performance profiles are built for each metric of interest m k c,f and for each formulation f ∈ F = {1, 2, ...7} , including the (F)* variants, to solve problem cases c ∈ C = {1, 2..n c } . The ratio of performance r k c,f compares the performance of a formulation against the best performance among all formulations for a specific case: When r k c,f = 1 , formulation f performs best for case study c with respect to metric m k , otherwise r k c,f quantifies formulation f performance ratio with respect to the best performing formulation (Dolan and Moré 2002). The performance of a formulation f with respect to metric m k is defined as the probability that the performance ratio r k c,f is at most a factor τ to the best performance ratio: In other words, k f ( ) is the probability that, in the worst case, formulation f performance is lower by a factor τ compared to the best performing formulation (Dolan and Moré 2002). n c denotes the total number of problem cases. Table 9 gives the formulation ranking with respect to m 1 for all n c = 10 cases where global solution optimality has been proven (i.e., all cases except (e) and (f) in Sect. 3.3). The figures indicate in how many cases a formulation f has performed at a certain rank. Figure 22 gives the corresponding performance profile. Generally, it is desirable to observe a significant performance increase for small variations of τ. When a performance profile grows fast as τ increases, it indicates that the corresponding formulation performs close to the best one over the entire test set. (F GG ) has 90% probability to perform with a ratio of at most 1.8 times the computation time of the fastest formulation. ( F BSF 1 )* has the highest probability to perform with a ratio of 1.5 (70%) and 1.0 (30%). Adding the displacement-based elongation bounds (big-M constants) brings no significant improvement in (F GG )*. Removing the displacement-based elongation bounds is (64) k f ( ) = 1 n c size{c ∈ C ∶ r k c,f ≤ τ} beneficial to ( F BSF 1 )* and ( F BSF 2 )* as their performance grows faster for small values of τ ≤ 3. Generally, formulations (F GG ), ( F BSF 1 ), and ( F BSF 2 ) perform well for strength-governed problems. As mentioned in (Bollapragada et al. 2001;Van Mellaert 2017) and shown in Sects. 3.3 and 3.4, displacement bounds can strongly influence computational performance. Formulations (F GG ), ( F BSF 1 ), and ( F BSF 2 ) require significantly reduced computation time compared to (F RS ) when moderate displacement bounds are set. However, when displacement constraints are tight, resulting in low utilization of element capacity (i.e., stiffness-governed problems), it becomes difficult to verify global optimality, regardless of the formulation (cases (e) and (f), Sect. 3.3). In those cases where global optimality has not been verified within the time limit, employing ( F BSF 1 ) has resulted in the smallest MIP-gaps, i.e., it provides the tightest lower bounds. Table 10 gives the formulation ranking with respect to m 2 , the number of explored search tree nodes. It is clear that ( F BSF 1 ) requires by far the least number of continuous LP relaxations to be solved. This indicates the efficiency of

Fig. 23
Performance profile 2 f , number of explored search tree nodes ( F BSF 1 ) to provide fast increasing lower bounds. In general, solving LP relaxations of different search tree nodes can be carried out independently on different computing cores (Gendron and Crainic 1994), which might further improve the performance of ( F BSF 1 ) when the number of CPU cores is increased. Figure 23 gives the corresponding performance profile. In this case, only the variants (F GG )* and ( F BSF 2 )* perform slightly better than the corresponding original formulations. In all cases, it has been observed that (near-)optimal integer-feasible solutions (upper bounds) are obtained very early in the branch-and-bound process and that most of the total computation time is spent to verify global optimality through increasing lower bounds. Table 11 gives the formulation ranking with respect to m 3 and Fig. 24 the corresponding performance profile. In most cases, (F RS ) has the worst performance. In this case, the addition of displacement-based elongation bounds (big-M constants) improves (F GG )* performance significantly, which is the fastest to obtain high-quality upper-bound solutions in 60% of all cases (τ = 1). Adding the displacement-based elongation bounds halves the performance ratio reducing it to 2.0 from 4.3 for a probability of 90%. On the contrary, removing the displacement-based elongation bounds is beneficial to the performance of ( F BSF 2 )*.
In general, this analysis confirms that the MILP formulations investigated in this work are effective to produce high-quality solutions in a short computation time. The study of the variants (F)* has shown that when the difference in the number of variables and constraints is small, the mathematical structure of the problem has a significant influence on the formulation performance.
Over the last few years, the performance of MILP solvers has improved significantly due to advances in method implementation and increased computing power. This is evident when comparing the computation times reported in this work with those given in the literature (cf. case studies 1 and 2 in Sects. 3.1 and 3.2 with those reported by Rasmussen and Stolpe (2008)). To the authors' knowledge, it is the first time that solution global optimality has been verified for the 25-bar tower truss design (Sect. 3.4).
From the cases studied in this paper, it can be inferred that (F GG ), ( F BSF 1 ), and ( F BSF 2 ) will perform better for larger-scale problems compared to (F RS ). ( F BSF 1 ) (and its variant) has the highest probability to verify global optimality in the shortest computation time. Nevertheless, solving MILP to global optimality with branch-and-bound methods is at worst-case exponentially complex, which might limit the application of the presented formulations to small-to medium-scale problems (< 200 structural members and < 50 candidate cross-section areas).

Limitations and future work
This work offers a benchmark of a new MILP-based formulation for sizing and topology optimization of truss structures against two well-known formulations. The benchmark

Fig. 24
Performance profile 3 f , computation time to reach upperbound solutions that are global optima is based on four problem cases (although with sub-cases to evaluate the influence of displacement constraints) of small and medium sizes. Future studies could deepen the generalization of the conclusions reached in this work by extending the benchmark to other problems. The purpose of the benchmark given in this work is to investigate the specificities of discrete sizing and topology optimization subject to stress and displacement limits and without the influence of additional constraints such as global stability (including mechanism avoidance), member overlapping, and buckling. In addition, such constraints are not included in the benchmark formulations. Other formulations that include constraints on global stability, member overlapping, or member buckling have been developed (Fairclough and Gilbert 2020;Mela 2014;Ohsaki and Katoh 2005). Future work could extend the formulation given in this paper by implementing such constraints and investigating their influence on computational performance and solution quality.
Solving MILPs through branch-and-bound algorithms can be parallelized to multiple computing cores (Bertsimas and Tsitsiklis 1997;Gendron and Crainic 1994;Gurobi Optimization, LLC 2019). Although computation time might not reduce significantly for large-scale problems owing to the exponential complexity of MILP (Nemhauser and Wolsey 1999), future work could look into evaluating the potential improvement margin through parallelization that depends on the formulation type. The examples analyzed in this work have shown that verifying global optimality is harder than obtaining optimal or near-optimal upper-bound solutions, especially when displacements are limited by tight bounds. Future work could (1) study the correlation between the MILP mathematical model and the effect of tight displacement bounds, (2) quantify further the formulation performances through additional studies, and (3) develop formulations that enable tighter lower bounds to reduce the MIP-gap faster. The latter might be achieved by adding other sets of redundant design and state variables as well as redundant constraints. Another approach that might lead to improvement could be that proposed by (Shahabsafa et al. 2018) who developed an incremental assignment formulation for discrete cross sections, which has shown to reduce computation time. It has been observed that certain settings of the Gurobi solver (Gurobi Optimization, LLC 2019) employed in this work may affect computational performance. In this work, most of the solver parameters have been kept at the default value. Future work could carry out parameter tuning to improve computational performance and benchmark results with other software (e.g., CPLEX (IBM 2020)).

Conclusion
This work has introduced a new MILP-based formulation for sizing and topology optimization of truss structures. A benchmark study with two well-known existing formulations has shown that the computation time to obtain global optima via a branch-and-bound solver highly depends on the chosen set of variables and constraints. The benchmark given in this paper, which is the first of this kind, has helped to gain a deeper insight on how to improve MILP formulations. Through appropriate formulation variations, it is possible to significantly improve computational performance. One of the existing formulations (F RS ) performs significantly worse compared to all other formulations. Instead, formulations that employ member elongations as state variables (F GG , F BSF 1 , and F BSF 2 ) are able to obtain near-optimal upper bounds and global optima much faster. The formulation (F GG ), previously developed by Ghattas and Grossmann (1991), has been improved by adding displacement-based elongation bounds (big-M constants), which enables obtaining high-quality upper-bound solutions in a very short computation time. The new formulation ( F BSF 1 ) is, generally, faster to verify global optimality and also requires exploring significantly fewer nodes in the search tree.
Funding Open access funding provided by EPFL Lausanne. This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

Conflict of interest
The authors declare that there is no conflict of interest.

Replication of results
All formulation and case study details required for replication of results have been reported in the article. In addition, MILP model files (.mps format) for all case studies are provided as supplementary data. These files can be executed with Gurobi and other commercial MILP solvers. The employed Gurobi-C# implementation will be made available by the authors upon request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.