Optimization of branching structures for free-form surfaces using force density method

ABSTRACT Branching structures are mechanically efficient in supporting large-span structures, such as free-form roofs. To support a roof with a specified geometry, we present a novel shape and topology optimization method to find the optimal branching structure in this paper. In the proposed method, the branching structure is modelled as a cable-net, while the reaction forces from the roof are taken as external loads. The force densities of the members are the design variables. The optimal branching structure can be obtained by minimizing one of the several proposed objective functions. The shape of the branching structure represented by the nodal coordinates is determined by solving the linear equilibrium equations. The topology is optimized by removing the members with small axial forces and incorporating the closely spaced nodes. The cross-sectional areas can be easily calculated, if the allowable stress is assigned. Hence, it is very convenient to simultaneously optimize the cross-section, shape, and topology of a branching structure. Numerical examples show that this method can be easily applied to a 2D problem. For a 3D problem, the constraints on the reaction forces should be relaxed. Considering the roof supports as variables is also an effective solution for 3D problems. Graphical Abstract

Optimization method for branching structures bracing free-form surfaces Shape optimization by force density method Topology optimization by removing the members with small internal force (The first innovation point)

Verification by numerical examples
The method is effective for 2D problem The method is effective for 3D problem after optimizing supports on surface (The third innovation point)

Discussions on objective functions
Sum of the strain energy Sum of the absolute value of the internal axial forces Sum of the absolute value of the force densities

Variations in member lengths
More uniform distribution of member lengths (The second innovation point)

Introduction
Branching structures are very useful for supporting architectural roofs for a large space, especially when combined with free-form surfaces. A branching structure can provide a large space near the ground for human activities using a small number of columns, which is preferable in view of architectural design. One famous tree-like supporting structure is the terminal building at Stuttgart International Airport, designed by German architect/engineer Frei Otto and built in 1960 (Cui, Zhou, and Ohsaki 2016). Recently, such kind of supporting structures has been widely applied in many public buildings. Nature is an inspiration source of architectural design. Obviously, the branching structures in architecture come from branches of trees, while many different approaches have also been undertaken by architects and engineers (Rian and Sassone 2014). The branches and nodes are the essential components for branching structures.
In the early 1990s, physical experiments by thread models have been used to explore structural forms for planar or spatial branching structures. One of the experimental techniques (Buelow 2007) uses dry strings connected with beads. The beads are held in place simply by friction against the strings. The strings can be fixed in a frame or tensioned with suspended weights. The beads allow repositioning of the nodes, and the lengths of the threads can also be adjusted. Obviously, in this model, different positions of the nodes result in a different tree-like structural form, although locations of the boundary nodes and topology of the structure are fixed. In the thread and bead model, the whole structure is in a tensile state.
With the development of computer technology, many researchers put forward many numerical methods, such as dynamic relaxation method (Barnes 1977;Barnes, 1999), genetic algorithms (Buelow 2007), formfinding methods with finite elements (Haug et al. 2009), sensitivity method (Cui, Zhou, and Ohsaki 2016;Cui and Jiang 2014), etc., to deal with disadvantages of the physical modeling techniques. Barnes (Barnes 1977;Barnes, 1999) introduced the numerical procedures of form finding method based on the method of dynamic relaxation with kinetic damping for long-span cable nets and grid shells, etc. Haug et al. (2009) summarized the status of design and analysis of industrial membrane structures in architecture, automotive, and space industries by several types of computational software. One of them is the Lightweight Structures Analysis, based on the Ph.D. work by E.J. Haug, 1969Haug, -1971 at UC Berkeley, which is an effective finite element code for shape finding and static/ dynamic analyses of the flexible membranes, beams, and cable structures. Wu et al. (Wu, Zhang, and Cao 2011) proposed the inverse-hanging recursive method, in which the primary idea is to find the load-bearing centers of the branch roofs at each level. The loadbearing centers in the adjacent level are then connected by lines to generate the geometric form of branching structures progressively. However, this method considered only the resultant forces of the loads on the specified roof areas, and it did not consider the interaction between the structural shape and the loads. Hunt et al. (Hunt, Haase, and Sobek 2009) proposed a design tool for two-dimensional tree structures, where the inner nodes are allowed to move only in the vertical direction so as to stabilize the analysis process. The topology is specified in advance and cannot change in the optimization process. Cui et al. (Cui, Zhou, and Ohsaki 2016) presented an optimal design method for a branching structure, generated by fractal geometry, to support a free-form shell. Only the shape of the branching structure was optimized by minimizing the strain energy, ie, the topology of the branching structure cannot change in the optimization process. Cui and Jiang (Cui and Jiang 2014) minimized the strain energy of the structure using the sensitivity coefficients, and presented a method for simultaneously optimizing the shape, topology, and crosssectional areas of framed structures, including branching structures. However, they did not show the effectiveness of finding the optimal shape of the branching structures supporting free-form surfaces.
In this paper, we propose a method for simultaneous optimization of cross-section, shape, and topology of the branching structures supporting free-form surfaces using the force density method (FDM), which is primarily applied to form-finding of cable-nets and tensegrity structures. The nonlinear equilibrium equations with respect to the unknown nodal coordinates are converted into a set of linear equations by introducing the concept of force density (Gruendig and Singer 2000;Ohsaki 2006, 2015), which is defined as the ratio of member force to its length. The force density method was published for the first time in the article in 1971 by Linkwitz and Schek (Linkwitz and Schek 1971). Schek (Schek 1974) expanded the method in 1974 for equidistant square cable nets. The methods nowadays used for form finding of tensile structures are mostly based on his method. Descamps et al. (Descamps et al. 2011) developed a constrained FDM to enforce geometric restrictions by a local approach based on static equilibrium equations. The geometry constraints were imposed on the positions of loaded nodes. Miki and Kawaguchi (Miki and Kawaguchi 2010) extended the FDM by minimizing various objective functions for finding the forms of complex tension structures composed of cables, membranes, and compressive members. Malerba et al. (Malerba, Patelli, and Quagliaroli 2012) proposed an extended FDM by setting conditions in terms of nodal reactions to deal with form-finding problems of cable-nets. Lee et al. (Lee, Lee, and Kang 2018) presented a numerical method for form-finding of tensegrity structures with multiple states of self-stress by using the force density method combined with a genetic algorithm. The design variable can be uniquely defined in the case of multiple states of selfstress using only the constraint of member types. Xu et al. (Xu, Wang, and Luo 2018) proposed an optimization approach for tensegrity structures using mixed integer nonlinear programming. Zhao et al. (Zhao et al. 2020) presented a strong coupled form-finding and optimization algorithm for reticulated shell structures. Wang et al. (Wang, Xu, and Luo 2021) introduced a general computational framework for the form-finding of tensegrity structures. These methods can handle tensegrity structures that have a force density matrix with a rank deficiency greater than the required minimum value.
The optimal topology of a truss is often obtained by removing the members with small cross-sectional areas from the initial ground structure that has many redundant members (Hagishita and Ohsaki 2009), while the locations of the nodes are fixed. The optimal shape and topology can be simultaneously obtained by considering nodal locations as additional variables. Recently, Ohsaki and Hayashi (Ohsaki and Hayashi 2017) explored the merits of FDM in shape and topology optimization of trusses; the objective function, namely, the compliance, and the constraint function, namely, the structural volume, are expressed explicitly by the force densities based on the fact that the optimal truss is statically determinate with the same absolute value of stress in existing members. Shen and Ohsaki (Shen and Ohsaki 2020) extended this method to simultaneous shape and topology optimization of planar frames. In their method, the nodal locations are expressed as functions of force densities of an auxiliary truss or a cable-net, and the cross-sectional areas of members of the primary frame are determined by solving a nonlinear programming problem. The members with small cross-sectional areas are removed to obtain an optimal topology. In other words, the topology of the structure remains unchanged in the optimization process, and the inefficient components are deleted after optimization.
From the inspiration that the thread and bead model is in a tensile state, we use the FDM in this paper to find the optimal shape of the branching structure. The design variables are only the force densities of the members. The initial branching structure is not a structure like a fractal tree, but a grid connected by bars. Concentrated loads are applied to represent the reaction forces from the free-form roof. For a 2D problem, the locations of loaded nodes on the roof are fixed in the optimization process, and the members with small internal forces are removed directly at each step in the process of topology optimization. The structural shape gradually changes in the direction of increasing the minimum internal force and improving the force transmission efficiency of components. This process is similar to the idea of evolutionary structural optimization (Huang and Xie 2010), fully stressed design, and the optimality criteria approach.
The method for 2D problems is then extended to 3D problems. To support a free-form surface, the branching structure should be in equilibrium with the reaction forces in z-direction, which are transferred from the roof surface. We first find the optimal locations of supports on the roof surface, which are kept unchanged while applying the reaction forces to the branching structure to find its optimal shape. However, the 3D problem is more complex than the 2D problem, because the vectors of reaction at different surface supports do not always intersect with each other. To resolve this difficulty, locations of the surface supports and the branching structure are optimized at the same time. Numerical examples demonstrate that this turns out to be easier to find an optimal branching structure.
The paper is organized as follows: In Section 2, we introduce the method of shape and topology optimization by FDM and discuss the properties of the method through examples of planar trusses. In Section 3, we solve the problem of optimizing the branching structure supporting a free-form surface. In Section 4, numerical examples are utilized to validate the effectiveness of the method. In Section 5, the conclusions are given.
Throughout the paper, the ith component of a vector a is written as a i , and the (i, j) component of a matrix B is denoted by B ij . All vectors are column vectors.

Shape and topology optimization problem
Shape optimization by the FDM is introduced first to obtain the optimal shape of branching structures as shown in Figure 1. It is easy to optimize the crosssection after obtaining the optimal shape. However, the cross-section optimization is not our concern in this paper, and all the existing members have the same cross-sectional area. For topology optimization, the members with small internal forces will be removed at every step of shape optimization.

Shape optimization by force density method
In this subsection, the FDM proposed in Ref. (Ohsaki and Hayashi 2017) for shape and topology optimization of trusses is presented for completeness of the paper. Consider a pin-jointed structure, which consists of m members and n nodes. The axial force of the ith member is F i , its length is l i , and the force density q i is defined as The force density vector q is given as Let x∈< n , y∈< n and z∈< n denote the vectors of x-, y-, and z-coordinates, respectively, of all nodes. The loaded nodes and the supports in the branching structure are classified as fixed nodes, and the other nodes are free nodes, indicated by the subscripts "fix" and "free", respectively. Moreover, the nodes are labeled such that the free nodes precede the fixed nodes, and we have Let C∈< m×n denote the connectivity matrix, or incidence matrix of the structure, which defines connectivity between the members and the nodes. If member i is connected to nodes j and k (j < k), then the components C ij of C is defined as follows (Zhang and Ohsaki 2015): while the other components in the ith row of C are 0. C can be divided into two parts as follows: The force density matrix Q for expressing the equilibrium in x-, y-, or z-direction is formulated as follows, which can be divided into four parts as Equation (6): Let p x ∈< n , p y ∈< n and p z ∈< n denote the nodal load vectors including the reaction forces in x-, y-, and z-directions, respectively. The equilibrium equation in x-direction, for instance, can be written in terms of nodal coordinates as where p free x 2 R n free ,p fix x 2 R n fix are the nodal loads applied at the free and fixed nodes, respectively. Note that the locations of loaded nodes are fixed in the following optimization problems; therefore, they are included in the fixed nodes and p x free is a zero vector. The following equations are obtained from Equation (7): where R x is the reaction force vector. If node k is the surface support as specified in Figure 1, a constraint is given so that the reaction force R xk is equal to the specified load P k from the roof surface. If the force densities of all members and the locations of fixed nodes are assigned, then the locations of free nodes are obtained from the set of linear equations in Equation (8). Therefore, x free is considered as a function of q. For the equilibrium in y-and z-directions, we have similar equations. Let R ∈< L denote the vector consisting of the reaction components corresponding to the specified load vector P, ie, the relation R = P should be satisfied. Here, L = 3 c for 3D problem, and L = 2 c for 2D problem, where c is the number of surface supports. Form finding based on the force density will offer the mechanical information to the topology optimization in Section 2.2 and the calculation of the objective function in Section 2.3. The form-finding process of the branching structure is summarized as follows: (1) Formulate the connectivity matrix C by Equation (4).
(2) Obtain the force density matrix Q for the given force density vector q, and also Q free , Q fix , Q link by Equation (6). (3) Compute the coordinates of free nodes x free , y free, and z free by Equation (8), where p x free ¼ 0, p y free ¼ 0, and p z free ¼ 0 are satisfied. (4) Obtain the reaction forces R x , R y , and R z by Equation (9). (5) Find the force density vector q by solving an optimization problem under constraints of R = P. More details can be found in Section 2.3.
After obtaining the force densities resulting in the specified reaction forces, we can easily determine shape of the structure, namely, the coordinates of all nodes, by solving Equation (8). Furthermore, the axial forces of the members can be obtained according to the definition of the force density in Equation (1).

Topology optimization of branching structures
The topology optimization has two kinds of operations: member removal and node incorporation, which are implemented after shape optimization at each step during the optimization process.

Member removal
After obtaining the nodal locations using the FDM explained in Section 2.1, the members with small absolute values of axial forces will be deleted from the branching structure except the members between the two loaded nodes (surface supports). When there are many members with zero internal forces, these members are eliminated firstly since they have little effect on the structural performance. Let F min denote the minimum absolute value; ie, where l' is the number of members in the branching structure except those between the surface supports. A reference value F c for the member elimination is set as to determine which members will be deleted, where γ is a parameter controlling the rate of convergence. Namely, a member is removed if its absolute axial force is not larger than the reference value F c , ie, the condition for member removal is given as The value of γ has some influence on the optimization results, because it determines the number of members to eliminate at each step. γis effective when there are few zero internal force members. A large value will make the evolution direction uncontrollable and cannot obtain a proper structural form. For instance, if it is close to 1.0, the internal forces of members will not become smaller any more in the final stage and maybe only one member is deleted in one step. Hence, γ should be a small value to make F c appropriately small. The convergence condition for optimization is assigned as where λ is a small positive number between 0.01 and 0.10. If λ is too large, the force transmission path of the structure may be destroyed, resulting in an unreasonable structural form. This threshold is difficult to determine because it depends on the structural type. The internal forces of branching structure are generally non-uniform; so λ should be small enough.

Node incorporation
When distance among some nodes is smaller than a specified small value, they will be incorporated into one node. The new coordinates are the geometry center of these nodes. Overlapping members caused by node incorporation will be incorporated as one member. Unless there are a lot of redundant members whose internal forces are very small, the elimination of members will not lead to removal of too many members. However, the node incorporation, which is a simple geometry operation, may delete many members as a result. Hence, the node incorporation will have a great influence on the optimization result when the initial structure has dense members.
With the reduction in the number of members, the structure will approach a statically determinate structure. If the stress constraints are specified, the crosssectional areas of the members will be easily determined by a two-level optimization approach, where the cross-sectional areas are modified in the lowerlevel optimization problem. However, the crosssection optimization is not our concern in this paper.

Objective functions
To let the branching structure efficiently support the surface, we are to find a branching structure with external forces P k being equal to the reaction forces R k as shown in Equation (14). The topology operation in Section 2.2 will remove inefficient members and maximize the minimum internal force of the components. In addition, the purpose of this paper is to optimize the structure for architectural design; therefore, we must pay attention to the aesthetic effect of the structural form in addition to structural performance. Different objective functions and the corresponding optimal structural forms are investigated for better structural performance while maintaining aesthetic effect.
In the following, we use a group of specified loads P, for simplicity, representing the reaction forces from the upper roof in the planar problems. The first optimization problem is given as Equation (14).
Opt2 : minimize Opt3 : minimize Opt4 : minimize where t c is the total number of the loaded nodes, and q L i and q U i are the lower and upper bounds of q i . To discuss their effect on the shape of the branching structure, we have three more objective functions as in Equations (15), (16) and (17). K is the set of the loaded nodes, and E i is the strain energy of the ith member.
Opt1 is just to confirm whether the force density method can find a solution that satisfies the force balance between the surface and the branching structure. If the problem does not result in a zero objective function, then the solution is infeasible because the reaction forces are not balanced. For more sophisticated optimal design of branching structures, we will consider other metrics in Opt2 -Opt5.The problem Opt2 is to minimize the sum of the absolute values of the internal forces, which corresponds to minimization of the number of members with non-zero internal forces. The objective function of Opt3 is minimizing the total strain energy. For a pin-jointed structure that has only axial forces, the strain energy of a member is proportional to the squares of the internal forces and to the member length. The objective function of Opt4 is the sum of absolute values of force densities, which is the internal force divided by the member length.
The bounds for q i are defined as follows: Here, q max is the maximum force density when the surface is supported by only one member connected directly to the base.P z k is the kth specified external force in positive z-direction for 3D problem and in positive y-direction for 2D problem. H is the height of the structure. k c denotes the number of layers of the branching structure. When an optimal force density q is obtained, the coordinates of the free nodes can be calculated by Equation (8) to obtain the optimal shape. The flowchart of the proposed method is illustrated in Figure 2.

Optimization of 3D structures
When supporting a free-form surface, the branching structure bears the reaction forces from the roof surface. We first find the optimal locations of supports of the roof surface under the negative z-directional uniform vertical load. Then the reaction forces are applied in the opposite direction to the branching structure to make it in a tensile state and find a reasonable shape by solving one of the optimization problems in Equations (14)-(17) using the FDM. However, the spatial (3D) problems are more complex than the planar (2D) problems. It will be challenging to find an optimal  solution in which the loads at the surface support nodes of the branching structure are equal to the specified reaction forces from the surface. That is to say, for the 3D problem, it is difficult to find an optimal solution which can simultaneously satisfy the constraints on the reaction forces in x-, y-and z-directions, when the locations of the surface supports are fixed.
The reason for this difficulty can be explained as illustrated in Figure 3. For a 2D problem, the reaction force vectors at two adjacent surface supports will always intersect at a point. However, for a 3D problem, the reaction force vectors of two adjacent surface supports are generally not co-planar, and there could be no intersection point. Namely, the intersection is a special case, and in general the reaction forces are not co-planar. Therefore, it is difficult to find an optimal force density q to realize the balance of the structure, especially when the locations of surface support nodes are fixed.
To resolve this difficulty, we investigate two types of optimization problems: one is to allow the surface supports to move on the surface by minimizing the strain energy, and the structural analysis includes the surface and the branching structure; the other is the alternate optimization of the locations of the surface supports and the shape and topology of the branching structure by the FDM under updated reactions from the surface in each optimization step. Examples show that the latter approach is preferable in view of convergence property.

Optimization of locations of surface supports
The locations of the surface supports are usually assigned following architectural requirements and preferences of architects. If the roof surface is supported by only vertical columns, it may be difficult to change their locations due to the requirements of architectural planning. However, for the branching structure, it will be possible and mechanically meaningful to adjust locations of the surface supports.
The shape of a continuous shell structure supported by the branching structure is described by using B-spline surface in this paper. The x-, y-, z-coordinates of the surface are derived from the following equations: In Equation (19), the parameters s and t are the orders of the B-spline basis functions for u and v (u, v∈ [0, Cui, Zhou, and Ohsaki 2016]), and m'+1 and n'+1 are the numbers of knots. α i;j ; β i;j ; γ i;j denote the (i, j) combination coefficients, namely, the control points, of the B-spline basis functions for the (x, y, z) coordinates. It is notable that the surface might not pass through the control points, making it difficult for the designer to intuitively describe the initial surface shape. Hence, we utilize the locations of key points Q i,j (x, y, z) (i = 0, 1, . . ., m'; j = 0, 1, . . ., n'), which are located on the surface shown in Figure 4, as design variables in this study. z-axis is defined in the right-hand system as in Figure 4. When the orders s and t of the surface are specified, the knot vectors for the parameters u and v will be determined by the chord length method. Then, the B-spline functions will be determined by the coordinates of the key points uniquely; ie, the combination coefficients α i;j ; β i;j ; γ i;j À � of the B-spline basis functions can be calculated according to the global surface interpolation method introduced in Chapter 9 of Ref. (Piegl. and Tiller 1996). Hence, a set of key points on the surface uniquely determines the shape of B-spline surface.
The surface is discretized into shell elements. The roof surface is first supported by the pin supports, and static responses are computed using the standard approach of finite element analysis. The objective function C, which is the total strain energy, is defined as Equation (20).
where F is the nodal force vector, which is computed from the uniform vertical loads on the surface, and U is the nodal displacement vector of the surface. The optimization problem for locations of surface support nodes is given as Opt5 : where X fix is the coordinate vector of the support nodes on the B-spline surface, which is provided by the branching structure. X fix i is the coordinate vector of ith node in X fix , whose coordinates in (u, v) space is (u i , v i ). The (x, y) coordinates of the supports (namely, the fixed nodes on the top of the branching structure) for the input of the optimization are given in accordance with the preference by the designer and converted into (u, v) coordinates in the optimization process by numerical approximation method. The optimal positions of the fixed nodes are found in (u, v)-space and converted into (x, y, z) coordinates to prepare the data in next optimization stage.

Successive optimization of support locations and shape and topology of branching structure
The initial locations of the pin supports of the free-form surface are uniformly distributed. The initial branching structure (namely, the cuboid grid) is also included in the initial structure. Next, the surface support nodes, namely, the top nodes of the branching structure, are taken as design variables. The strain energy is minimized to find the optimal support locations by solving Opt5. The reaction forces in z-direction at the supports are converted into upward nodal forces (reactions in x-, y-direction are converted simultaneously) and applied at the surface support nodes of the branching structure to optimize the shape and topology of the branching structure by FDM as described in Section 2. This process is repeated until the convergence condition in Equation (13) is satisfied. The flowchart is shown in Figure 5.
As indicated in Figure 5, Opt5 can change the locations of the supports, which cannot be dealt with by the force density method itself. From the point of view of the optimization cycle, Opt5 will improve the position of fixed nodes in next optimization cycle after the shape and topology optimization by force density method. The whole optimization process is a step-bystep iterative approximation process. Finally, a better solution will be found, not the optimal structure.
The applied nodal loads, namely, reaction forces from the surface, are updated at every optimization step of the branching structure. These reaction forces may change due to the stiffness change of the branching structure, when the shape and topology are updated.
The alternate process of analysis and optimization can reflect the influence of such change in the optimization process. The reaction forces at the surface support nodes will not vary too much and remain positive.

Numerical examples
The effects of different objective functions are first investigated in Example 1. Secondly, a branching structure with two supports is used in Example 2 to validate the effectiveness of the proposed method in finding a reasonable structural form with multiple supports. Finally, a spatial branching structure is optimized in Example 3. In these examples, γ is 1.06, and λ is 0.1 if they are not specified. The self-weight is not considered for the branching structure. We solve the optimization problem using the nonlinear programming library fmincon provided in the Optimization Toolbox of MATLAB R2018a using the interior-point method. Structural analysis, with the truss elements for the branching structure and the shell elements for the roof surface, is carried out using OpenSees (2020).

Example 1: Planar structure with a single support
A 4 × 4 m rectangular grid with 1 m spacing distance is shown in Figure 6. The central point at the bottom is fixed in x-and y-directions. The concentrated y-directional loads P i = 10 N are applied at the five nodes in the top layer. The nodal load in x-direction is equal to 0. The initial force densities for all the members are 1.0 N/m, and the bounds for the force densities are specified as 0:01 � q i � 5q max for the branching structure, and À    5q max � q i � 5q max for the members in the top layer. Here, q max is determined by Equation (18). Optimal shapes for the problems Opt1-4 are shown below.
Opt1: Even if P x = 0 and P y = 10 N are specified simultaneously, a branching structure as shown in Figure 7 can be obtained. This optimization problem directly finds a solution satisfying the constraints by minimizing the sum of the errors. There are many members in this result, and the force densities and the lengths of the members are not uniform.
Opt2: Using the sum of absolute values of the member forces as the objective function, we can obtain the optimal shape as shown in Figure 8, for which the member lengths are not uniform. The number of members is smaller than the optimal solution for Opt1, because fewer members directly lead to a smaller objective function value. This objective function does not consider the length of members, and the structure often has long and short members, which is considered to be an unreasonable form for practical applications.
Opt3: All of the members have the same crosssectional area A = 0.02 m 2 and elastic modulus E = 2.1 × 10 4 N/m 2 . Note that this elastic modulus is not related to a practical material. The optimization result is shown in Figure 9. The objective function, which is the sum of strain energy in members, has a similar role as the objective function of Opt2, where the member forces have a significant influence. On the other hand, minimizing the objective function reduced the total lengths of members, and consequently, leads to a structure with fewer members. By observing the optimization results in Figures 8 and 9, we can see that the optimal structure of Opt3 does not have the short members in the middle, which exist in the optimal structure of Opt2.
Opt4: By minimizing the sum of absolute values of the force densities, the optimal solution has a nearly uniform force density distribution compared to the previous three cases. We notice that the maximum value of the force densities is smaller than those in other cases. The reason is that there are three members connected to the fixed node at the bottom of the structure. If the force density vector q in the objective function does not include the top horizontal members, the solution will not be unique because the force densities in the top members are the important for the whole structure's balance and affected by those in the branching structure. Moreover, the member lengths are neither too long nor too short, and the loaded nodes tend to connect the fixed node at the bottom directly as shown in Figure 10. The short members appear in Figure 8 for Opt2, but disappear in the optimization result of Opt4 in Figure 10. However, it has longer members than those in Figure 9. Hence, it will be used in the spatial problem.
The parameters γ and λ have great impact on the optimization results. The reference value γ for member removal should be small to reduce the number of removed members at each step. The elimination process sometimes removes only one member, and this will result in an asymmetric structure. The horizontal member on the top of the branching structure is important to ensure existence of the optimal solution. If the horizontal member is removed, it will be difficult to find a solution which satisfies the reaction constraints. Only the solution satisfying the constraint Py = 10 N can be found without satisfying the x-directional reaction constraint Px = 0.
Example 2: Planar (2D) structure with multiple supports and curved roof  The loaded curved line is generated by the parabola y = 0.25(x-5) 2 + 3.5, which has a span of 10 m. The concentrated load P k = 10 N is applied at each node in y-direction. There are two pinned supports at the bottom, whose distance is 6.0 m. The initial structure is shown in Figure 11. The initial force densities for all the members are 1.0 N/m, the lower and upper bounds of the force density are the same as Example 1. The Opt1 is just to verify to find a solution that meets the balance, so we did not apply Opt1 to example 2. The optimal shapes for Op2-Opt4 are shown in Figure 12. The shapes of the optimal structures Opt2 and Opt3 are similar. The force densities in Opt4 are smaller and more uniform than those in Opt2 and Op3. This demonstrates that the method effectively finds the optimal branching structure with multiple supports by solving Opt4.

Example 3: Spatial (3D) structure
In the 3D example as shown in Figure 13, the orders s and t of the surface are specified as s = t = 3. The key points on the roof are composed of 6 curved lines with the same shape; only the z-coordinates of the key points are different. Table 1 lists the coordinates of the key points of the first curve line. The second curve line is 0.2 m higher than the first. The third curve line is 0.1 m higher than the first. The fourth and sixth curve line is the same as the first. The fifth curve line is 0.1 m lower than the first. The spacing of key points on the top surface are 1.0 m in x-direction and 1.2 m in y-direction. The spacing of the grid under the surface is 2.0 m in both x-and y-directions. The coordinates (m) of the four supports at bottom are (2, 2, 0), (8, 2, 0), (2, 4, 0) and (8, 4, 0). The surface shape is shown in Figure 13(a).
A spatial grid is taken as the ground structure. Cuboids are added under the surface. The surface supports shown in Figure 13(b) are connected to the   top nodes in the cuboid grid, which is determined by the distance between the supports and the top nodes of the cuboid grid. The distance is equal to twice of the shortest distance between one surface support node and the nodes in the grid. The surface supports are also connected to each other by the members generated by the Delaunay triangles.
The surface bears uniform vertical loads. Some key works need to be done as follows: Firstly, the (u, v)coordinates corresponding to the (x, y)-coordinates of the supports are calculated by the distance approximation solved by the golden section method. Secondly, the mesh in (u, v)-space are generated by the Delaunay triangulation, the (x, y, z)-coordinates of the nodes in the mesh are calculated by Equation (19), and the nodal loads are obtained by covering the area at each node. Thirdly, the objective function is calculated by Equation (21).
The members of the branching structure are steel pipes with diameter D = 100 mm and thickness t -= 10 mm. All cross-sectional areas of the members are 28.274 cm 2 . The second moment of area is I y = I z = 289.912 cm 4 , the polar moment of area is I p = 579.624 cm 4 , Poisson's ratio is 0.3, and the mass density is 7.9 × 10 3 kg/m 3 . For the shell roof, its thickness is 0.1 m, the elastic modulus is 3.0 × 10 10 N/m 2 , Poisson's ratio is 0.2. The initial force density is 1000.0 N/m for all members. The bounds for the force densities are the same as Example 1. The vertical uniform load is equal to 1500 N/m 2 , which is the value per unit horizontal area in the negative z-direction. We find a feasible solution for the 3D structure satisfying R = P (for all x-, y-and z-directions) following the flowchart in Figure 5. The convergence condition is set as F min = 0.1F max .
To ensure convergence of the optimization process, the constraint tolerance value ε is set according to the average reaction force at the surface support nodes as where F T is the total load of the surface. A is the horizontal projected area of the surface. n top is the number of supports on the surface. The final result is shown in Figure 14. This example demonstrates that the optimization for finding a branching structure is a highly nonlinear problem. We used beam element or truss element to simulate the top member. The results showed that the bending moment in the top members using the beam elements will be transferred to the branching structure and cannot be ignored. When truss elements are used, instead, in the optimization process, the bending moment in the final confirmation analysis can be ignored, as shown in Table 2. In this table, N is the absolute value of the axial force; M is the maximum bending moment at the two ends of one member in the same direction with respect to its local coordinate system; M ymax and M zmax are the maximum bending moments for the two ends in its local y-and z-directions, respectively; e is the eccentricity of the axial force. This result shows that a mechanically efficient branching structure supporting the free-form surface with mainly axial forces has been generated through optimization.

Conclusion
For the purpose of supporting a free-form surface elegantly, a novel method for form-finding of the branching structure has been proposed in this study. In this method, Table 2. Larger eccentricity e caused by bending moment for the top members. In this table, N is the absolute value of the axial force; M is the maximum bending moment at the two ends of each member in local y-or z-direction. the shape and topology of the branching structure are simultaneously optimized by using force densities as design variables. If we impose stress constraints on the members, their cross-sectional areas can be determined from the force densities. Subsequently, the proposed method is a powerful tool for design of a branching structure with optimized cross-section, shape, and topology. At every step of the shape and topology optimization, the reaction forces in z-direction transmitted from the upper surface are converted into upward loads and applied as the external loads to the branching structure. The balanced shape of the branching structure is determined using the force density method (FDM), originally developed for tensile structures. The novelty of the paper is summarized as follows. First, the topology optimization is realized by directly removing the members with small absolute values of axial forces and merging the closely spaced nodes to a single node. This approach is useful when dealing with complex topology optimization problems, in which the initial structure has many nodes and members. Second, three kinds of objective functions are discussed by planar examples. The sum of the absolute value of force densities of all members is the most effective. Third, for the 3D problem, it might be difficult to obtain a shape of the branching structure that can be in equilibrium subjected to the specified reaction forces from the roof surface. To resolve this difficulty, we proposed a strategy of optimizing the locations of surface supports, in the shape and topology optimization of the branching structure. This method has been demonstrated to be an effective way in the 3D example.
There are two points of concern. First, the termination tolerance was relaxed. If we use a fixed value as the termination tolerance, it could lead to difficulty in convergence, especially when the objective function value is very large. To solve this convergence problem, the tolerance value of the constraint was relaxed according to the average nodal loads applied at the surface supports. Second, the top (auxiliary) members of the branching structure are beneficial to find an optimal solution in the 3D examples. The reason is that they can adjust the directions of reaction forces at the surface supports of the surface and make it easier to find the optimal shape satisfying the constraints on the reaction forces. The top members can be deleted after obtaining the optimal branching structure for practical purposes, because the axial forces of the top members will be undertaken by the in-plane forces of the free-form shell.