A nodal‐based evolutionary optimization algorithm for frame structures

This work proposes a nodal‐based evolutionary design optimization algorithm to design frame structures whose edges are the Delaunay triangulation of homogeneously distributed nodes in the design domain. The remaining nodes can freely sway in the design domain except for the loading nodes and boundary nodes. As a result, it can extend the space of admissible solutions to this optimization problem and reduce the number of design variables. Then, the sensitivity of the objective function, namely, the sum of compliance and its volume, is input into the method of moving asymptotes to update the nodal coordinates and member thickness. The most inefficient node is deleted in each iteration based on the average nodal strain energy until its number reaches a prescribed limit. 2D numerical examples for the Michell arch, L‐shaped bracket, and Messerschmidt–Bölkow–Blohm beam show that the proposed algorithm can get the optimal structure within a few iterations. Compared with initial configurations in 3D, the optimal crane arm, transmission tower, and aquatic dome have less strain energy and fewer materials, showing a great potential of the proposed method to design actual space frames.


INTRODUCTION
Discrete structures such as trusses and frames are prevalent in engineering because of their excellent material saving capability and maximum load capacity utilization. Nowadays, structural engineers are eager to reduce weight without sacrificing its excellent mechanical properties and satisfying various stress, vibration, and buckling constraints. Structural topology optimization can lightweight structures, minimize strain energy, and tailor the effective material properties by designing the optimal material layout in a prescribed domain (Cheng & Olhoff, 1981;Zhou & Rozvany, 1991 by designing the void regions distributed in the design domain (Bendsoe & Kikuchi, 1988). The method of solid isotropic material with penalization (SIMP) uses the penalized element density as design variables to optimize structures with maximal stiffness (Bendsøe, 1989;Rozvany et al., 1992). Xie and Steven proposed the evolutionary structural optimization (ESO) method to generate holes in a uniform isotropic continuum by removing invalid elements based on some optimality criteria (Xie & Steven, 1993). M. Y. Wang et al. proposed the level set method to express structures as the zero-level contour of a level set function (M. Y. Wang et al., 2003). These methods have made stunning achievements in designing continuum structures. Unfortunately, they are not suitable for optimizing truss and frame structures featured with discretized members (Bendsoe & Sigmund, 2013). The investigations on the optimization of discrete structures are significantly earlier than the continuum structure. As early as 1904, an Australian mechanical engineer, Michell, tracked down the lightest or stiffest frames by barring unnecessary elements from the discrete elements (Michell, 1904). Unfortunately, discrete structural topology optimization lags continuum topology optimization because of various engineering constraints. Some optimum configurations attained through optimization are hard to be used in practice (Beghini et al., 2014). For instance, the American Institute of Steel Construction specifies the size and shape of steel members (American Institute of Steel Construction, 2005). The length and weight limitations of the steel bars in construction are essential in the optimal design (Guest et al., 2004). Also, it may be necessary to reserve pipes or holes in the steel bar for water pipes and electrical lines (Almeida et al., 2009). In addition, the architectural design requires aesthetic consideration, which puts forward higher requirements for discrete structural optimization (Gersborg & Andreasen, 2011).
Thus, optimization problems in civil engineering are often complex and non-linear, so these suitable optimization tools are required to find high-quality and viable solutions (Fister et al., 2013). The vibration of high-rising buildings under various accelerograms was well-controlled in multi-objective optimization, which employs a modified neural dynamic model of Adeli and Park (Gutierrez Soto & Adeli, 2017). They also used neural dynamics models to design space truss structures, successfully making a tradeoff between safety and cost . They adopted the data-parallel computing to use the neural dynamics model in a large-scale steel structural optimization with 20,096 members (Park & Adeli, 1997). They also applied the genetic algorithm (GA) to large realistic discrete structures (Adeli & Cheng, 1994). GAs were also employed to optimize skyscrapers and reduce costs by up to 15% (Aldwaik & Adeli, 2014). Kim and Adeli adopted the floating-point parameter instead of the binary parameter in the GA, which significantly reduced the computational cost of the cost optimization of composite floors (Kim & Adeli, 2001). A two-phase GA improved performance by imposing heuristic restrictions in designing large-scale steel structures while retaining the original aesthetics (Kociecki & Adeli, 2014). Inspired by their outstanding work, optimization techniques have been widely used in structural design. Mathakari et al. developed a GA to perform a reliability-based structural optimization for transmission towers (Mathakari et al., 2007), getting the Pareto-optimal set in one trial. Guo and Li systematically studied the effects of four discrete structural optimization algorithms on the design of power transmission towers (Guo & Li, 2011). They added topology optimization strategies (mirror node and layering) to meet particular constraints, including structural quarter symmetry, valid topology, and shape. París et al. (París et al., 2012) proposed the size/shape optimization framework of a transmission tower with multiple loading conditions based on the modern heuristic algorithm. A similar strategy was used to design transmission towers with pre-designed components (de Souza et al., 2016). The above cases illustrate that a discrete member optimization method using 1D elements is essential to practical application value. Beyond architectural design, Xia et al. used the discrete topology optimization methods for the strut-and-tie model generation, which can help clarify the transmission mechanism of force (Xia et al., 2021). Lu et al. combined multi-material topology optimization and cross-section size optimization to design a lighter and safer bus frame based on static stiffness and dynamic frequency stiffness constraints (Lu et al., 2019). Reintjes and Lorenz used discrete topology optimization to design additively manufactured lattice structures (Reintjes & Lorenz, 2021).
Since nature-inspired algorithms could adapt to almost any condition, they can be used in truss and frame optimization (Chalotra et al., 2015). Pioneering work in this field includes the university course scheduling problem solved by a selective search method based on particle swarm optimization (Hossain et al., 2019) and the traveling salesman problem using the discrete spider monkey optimization algorithm (Akhand et al., 2020). The spiral dynamics algorithm, used to simulate the spiral phenomenon (Tamura & Yasuda, 2011), was applied to solve a series of challenging engineering optimization problems (Siddique & Adeli, 2014a). The advantages of this algorithm lie in its ease of implementation, less control parameter requirement, and better reinforcement strategies. The water drop algorithm effectively solves combinatorial optimization problems by simulating the behavior of water droplets moving through the soil (Siddique & Adeli, 2014b). The simulated annealing algorithm is based on an analogy with the physical annealing of materials that avoids local minima problems (Delahaye et al., 2019). The simulated annealing algorithm is well-suited for solving optimization problems of objective functions based on complex simulation processes (Siddique & Adeli, 2016).
The most commonly used discrete structural design method is the ground structure method (Dorn, 1964). In this method, the number of members reflects the utmost connection possibilities of specified nodes in the design domain. Then, it calculates strain energy stored in each member to thicken vital members and remove unimportant ones accordingly (Hagishita & Ohsaki, 2009). Such a network of potential members named ground structure (or structural universe) should contain as many members as possible to cover all admissible optimal structures. For instance, it must have n(n − 1)/2 members for an n-node frame with non-overlapping members. The nodal number used in the ground structure method should be large enough to avoid missing the workable connection. It is possible to improve optimization efficiency by designing the nodal truss structures (Dobbs & Felton, 1969). However, due to the non-linearity of the combined size and geometric optimization problem, the initial nodal position and its potential connectivity still significantly affected eventual results (Bendsoe & Sigmund, 2013). Martinez et al. proposed a growth method to alleviate this shortcoming. The growth method used ordered nodes and elements to design discrete structures. It justified the analytical solution to the well-known Michell structures with the minimum weight (Martinez et al., 2007). However, this heuristic optimization algorithm may not always converge to a near-optimal solution, though it was relatively fast in some scenarios (He & Gilbert, 2015). Besides, this method became very inefficient when inserting a new joint for a structure with many members.
Researchers have used nodal positions as design variables to optimize discrete structures. Nan et al. used this concept to reduce the dependence of the ground method on the initial nodal position (Nan et al., 2020). Though this method produces selectable optimal solutions through Pareto optimal solutions, it still largely relies on the initial structure. Fenton et al. applied grammatical evolution, representing a variable nodal number and their locations on a continuum (Fenton et al., 2015). However, their nongradient-based method was inefficient. H. Wang et al. proposed a node-shifting scheme for grid-shell shape optimization (H. Wang et al., 2021). Some discrete topology optimization algorithms changed node connections by deleting nodes or elements. Xu et al. presented a practical method to finding optimal trusses by categorizing nodes and members and deleting the meaningless ones (Xu et al., 2003).
Haobin et al. proposed a non-gradient optimization method to delete nodes based on the number of connected members of trusses (Haobin et al., 2010). This idea of gradually deleting elements to get an optimal structure is like the ESO method, an efficient and easy-to-use continuum topology optimization method first used to optimize continuum structure. Chen et al. developed a parallel evolutionary continuum structural optimization algorithm to update nodal position and used a fixed-shaped circular cavity in the low-stress area to delete nodes (Chen et al., 2002). This excellent work enables structural topology to change naturally by gradually moving the node of the continuum element on the cavity boundaries (Chen et al., 2005). Huang and Xie gradually removed low-stress materials to shape the optimal mesh-independent structures based on nodal sensitivity other than elemental sensitivity (Huang & Xie, 2007). A similar concept was used in strutand-tie models to design the optimal truss concrete structures (Lanes et al., 2019). By combining the concurrent ESO algorithm and size optimization, Tomšič and Duhovnik obtained an excellent truss structure (Tomšič & Duhovnik, 2014). The ESO method also could quickly delete inefficient truss elements to get an optimized truss structure (Manickarajah et al., 2000). Combined with the ground structure method, this method could delete members with the lowest strain energy (Hu & Cheng, 2014;Kosaka et al., 2016).
Inspired by the ESO method, this work develops a nodalbased evolutionary frame optimization (NEFO) method to design frame structures efficiently. Besides the member thickness, this method uses nodal coordinates as design variables. The sensitivity of the objective function, namely, the sum of compliance and the frame volume, is analyzed using the adjoint variable method. After updating design variables using the method of moving asymptotes (MMA), the most inefficient node is removed based on the average nodal strain energy. In each iteration, the nodes only can shift within an acceptable margin to avoid numerical instability. Subsequently, the remaining nodes are re-meshed with the Delaunay triangulation algorithm to form the new frame members. Compared with the other techniques with fixed nodes, the proposed method can search for the optimal nodal position in the whole design domain. Thus, it extensively extends the admissible solution space and does not rely on the initial layouts.
Besides the Introduction section, the optimization problem and its sensitivity analysis are given in Section 2. Section 3 presents the element generation and nodal removing scheme. Section 4 describes the numerical implementation. Numerical examples in 2D and 3D are demonstrated in Sections 5 and 6, respectively. Finally, Section 7 concludes the article.

OPTIMIZATION DESIGN
Topology design optimization problems for discrete structures include truss structures and frame structures. This work focuses on frame structural optimization. Frame optimization attracts more attention than truss as they can withstand bending moments. In fact, in most cases, the bending of structural elements is inevitable (Xia et al., 2020). The most fundamental performance of the frame structure is structural stiffness. Therefore, stiffness maximization or compliance minimization is an important goal for obtaining a stable structure. On the other hand, having minimum weight is also one of the most critical factors in structural design.
For these two aims, the costing function in structural optimization is the sum of compliance C(x,d) and total member volume V(x,d), which is also the weighted sum method (Stadler, 1979), given as f ind = { 1 , 2 , … , } = { 1 , 2 , … , } min ( , , ) = 1 + 2 s.t. − = 0 (1) The displacement field u is the solution to a linear elastic problem with the stiffness matrix K and the loadings vector F. The design variable x is the nodal coordinate vector of M nodes. x k (k = 1,2,. . . ,M) is the coordinate vector of the kth node with an unknown location in the design domain with M as the number of such nodes, d is the vector of the cross-sectional diameters of the bars, and d j (j = 1,2,. . . ,N) is the coordinate vector of the jth frame element with N as the number of such elements. This multi-objective optimization problem uses the weighted sum method to solve this problem. The weighting factors w 1 and w 2 have to be positive (Marler & Arora, 2004).
The total volume and compliance are given, respectively, as where N is the number of elements, u e is the elemental displacement vector, x e is the elemental nodal coordinate vector, and d e is the elemental diameter. L e is the elemental length, and the elemental compliance C e can be defined as K e is the elemental stiffness matrix related to the nodal coordinate vector x e and the elemental cross-section diameter d e .
For a frame with M nodes and N elements, the number of the design variables is 2M + N in 2D because a node has two coordinating design variables, and a member has one thickness variable. In 3D, the number of design variables becomes 3M + N.
Because MMA (Svanberg, 1987) is used to update the design variables, the first-order derivative of the costing function needs to be derived. For simplification, the derivation process is focused on a 2D problem ( Figure 1). Note that the nodes could be shifted in two directions, namely, F I G U R E 1 Schematic of a set of p members sharing the kth node x and y. In addition, the cross-sectional diameters are also used as a design variable. Thus, the global sensitivity information should be the assembly of the sensitivity to three variables; that is, For simplicity, γ represents a specific design variable. The sensitivity of the objective function with respect to the design variable is The difficulty in calculating Equation (6) is the partial derivative of the displacement u about the design variables. Generally, the adjoint variable method is used to derive this partial derivative (Choi & Haug, 1983;Komkov et al., 1986). It has been used in node-movable optimization (Ding et al., 2017;H. Wang et al., 2021). To use this method, the governing equation, − = 0, weighted by a Lagranginia multiplier λ is added into the objective function as It is necessary to solve the derivation of the Lagrange function to a design variable γ because it is necessary to find a stationary point; at that point, the derivation is equal to 0, given as It is noted the partial derivative of the load F about the design variable is zero because of its invariance in opti-mization. The first-order derivation of the function is Rearranging the previous Equation (9) leads to factoring out the derivation of the displacement The adjoint term should be equal to zero to remove the derivation of displacement, which means that our added unknown variable has to satisfy Thus, The derivative of g of the entire structure for x k can be written as The derivative of g of the entire structure for d j can be written as where the subscript e represents the eth frame element. In that manner, the global sensitivity becomes

NODAL-BASED EVOLUTIONARY FRAME OPTIMIZATION
The key idea of the NEFO algorithm is to evolutionarily remove the most inefficient node during updating nodal coordinators and elemental diameters. As shown in its flow chart in Figure 2, each iteration contains three major steps, including meshing initial nodes, sensitivity analysis, and deleting the most inefficient node. In Figure 2a, the initial nodes are distributed in the rectangular design domain. Using the Delaunay triangulation algorithm, a network of frame members with fixed/pin nodes (gray in Figure 2b) and a loading node (white) could be obtained. The frame structure generated by the Delaunay triangulation algorithm can avoid extremely small internal angles or extremely short members (Lee & Schachter, 1980). Consider a frame structure generated by Delaunay triangulation has m members, n nodes, r reactions, and s hinges; the Maxwell criterion (3m + r) -3n -s ≥ 0 is satisfied because of m ≥ 2n -3 (Fenton et al., 2015), and the structure is geometrically stable. Thus, a statically determinate or indeterminate frame structure could be obtained when updating nodal connections in optimization using the Delaunay triangulation algorithm.
After the frame members are obtained (Figure 2b), the finite element analysis is performed to calculate the nodal displacement and sensitivity in Equation (15). Except for boundary nodes and loading nodes, the rest nodes named free nodes are subsequently moved to their new position as shown in Figure 2c using MMA. After that, the most inefficient node is deleted based on the average nodal strain energy, which is mathematically determined as where p is the number of members shared by node n. For instance, Node 1 (dark star nodes in Figure 2c) is shared by p = 5 members. The idea of node deletion is inspired by element removal in the well-known ESO method (Xie & Steven, 1993). The method iteratively analyses the structure using the finite element method and then removes "inefficient" elements. The efficiency of elements is usually evaluated as a function of their mechanical properties and responses, like strain energy, where elements with the lowest strain energy are considered inefficient. The node with the lowest average strain energy is regarded as inefficient in this proposed method because the frame elements connected to this node cannot effectively use materials to share the strain energy of the overall structure. Deleting nodes changes the topology of the structure, reduces design variables, and improves calculation speed. Thus, its average strain energy is the average compliance of its neighboring members. The displacement used in Equation (16) could be obtained from the finite element analysis in the sensitivity analysis. However, our program will perform finite element analysis again to obtain more accurate nodal strain energy for such an updated structure. In order to avoid an unstable structure, neither the loading nodes nor the boundary nodes can be deleted. Node removal naturally results in misfunction of related members, and therefore the structure could become F I G U R E 2 The flowchart of a representative iteration in the nodal-based evolutionary frame optimization algorithm. (a) The initial nodes in the design domain (gray nodes are pined/fixed, the loading is imposed on the white node); (b) the network of frames generated by the Delaunay triangulation algorithm; (c) nodes at their new positions and frames with new diameters updated by the method of moving asymptotes method; (d) the frame structure after nodal deletion; (e) the frame structure after Delaunay triangulation unstable in some cases. Thus, the NEFO method again used Delaunay triangulation to establish frame elements between these nodes to produce a stable system as shown in Figure 2e. However, using all nodes as the input of the Delaunay triangulation algorithm to get a new frame structure may greatly change the nodal connection, leading to discontinuities in optimization. Because Delaunay triangulation always generates a convex geometric shape, not all members are contained in the non-convex design domain. In order to avoid this problem, this method performs a local Delaunay triangulation that only considers the affected nodes (dark square nodes and gray square nodes in Figure 2c) that belong to the members sharing the most inefficient node in the same frame elements. For instance, if Node 1 in Figure 2c is the most inefficient node to be deleted, Nodes 2, 3, 4, 5, and 6 are affected. They are input into Delaunay triangulation to generate a local network as shown in Figure 2e. The diameter of the newly generated element is where V del is the total volume of the deleted frame elements, and L new is the total length of the new frame ele-ments. This strategy allows the algorithm to adapt to structures with non-convex design domains such as the Lshaped bracket. More importantly, it significantly saves computation time spent in the second Delaunay triangulation. This iteration will continue until the nodal number reaches the limit. The influence of the nodal number will be discussed in Section 5.1 because it is key to the optimal structure.

NUMERICAL IMPLEMENTATION
This section will introduce how to implement the NEFO algorithm from objective function normalization, sensitivity calculation, and design variable updating. In order to avoid too much difference between the volume and the compliance, it is necessary to normalize two items of the objective function; that is, where C 0 and L 0 are the total compliance and total volume of the initial structure, respectively. With the normalization process, the non-dimensional objective function obtains an upper limit of one and a lower limit of zero.

TA B L E 1 Pseudocode of sensitivity analysis
Since the objective function is normalized, the corresponding sensitivities also need to be transformed; that is, The degree of freedom (DOF) of the entire structure will also change with the nodal number, making it difficult to obtain an explicit objective function expression. However, the sensitivity of an arbitrary node only relates to itself and the affected nodes. Therefore, a similar concept of assembling the global stiffness matrix in finite element analysis is used to perform sensitivity analysis. For instance, the sensitivity of Node 1 in Figure 2b depends on the coordinates of itself and the affected nodes (Nodes 2∼6). The pseudocode to calculate the sensitivity is given in Table 1. In the previous section, the sensitivity to update the design variables is derived in Equation (15). After entering the properties of the frame, the frame element stiffness matrix is immediately assembled, the symbolic expression for the derivative of the design variables, ∕ is also "pre-assembled" as explicit expressions in MATLAB; that is, dgx. After applying boundary conditions and finite element analysis, the displacement u of each node can be calculated. The sensitivity of the coordinates of each node to the objective function is output in a matrix form dfdx, so the size of the matrix size is consistent with nodal coordinates xval (Line 1). This algorithm calculates one by one in the order of frame elements (Line 2). The algorithm first extracts the nodal coordinates x (Line 3) and the member node vector [p 1 , p 2 ] for the two nodes (Line 4) of the nth element. Each node has two coordinate values corresponding to two design variables. (Line 5), two translational DOFs and one rotational DOF, and its DOF vector df is obtained (Line 6). Through df, the displacement vector ue of the corresponding node can be extracted from the global displacement vector (Line 7). After that, the corresponding elemental sensitivity matrix W of size 5 × 1 is calculated through the built-in eval function of MATLAB. The evaluation of dgx requires nodal coordinates, elemental diameters, and displacement (Line 8). According to dv, this elemental sensitivity matrix is assembled to the corresponding position of the sensitivity matrix dfdx (Line 9). The algorithm will calculate the sensitivity of all elements. This "pre-assembled" explicit symbolic expression can significantly reduce repetitive calculations and improve computational efficiency.
This optimization is a discrete optimization and includes the process of node deletion. Provided that the allowed range of node movement is not too big, the error introduced by this method is acceptable. Ideally, the nodal position can vary freely in the entire design domain at each iteration, making the algorithm unstable for many variables. All nodes are only allowed to move within a small margin to avoid abrupt shape changes, namely, where the move limit δ is a positive value to be determined case by case and enforces small variations from one iteration to the next, which results in a more cautious progression toward the optimum. In order to ensure the convergence of the results and avoid the fluctuation of the results, when the nodal number is equal to the preset value, the proposed method adopts a gradually decreasing movement limit; that is, where α is the reduction factor. Generally, α is equal to 0.9. In some designs, the structure needs to maintain symmetry. However, due to numerical errors, symmetric nodes often deviate after one optimization iteration. This deviation may be gradually amplified in subsequent iterations and eventually make the structure lose symmetry. Therefore, symmetry constraints are necessary. After each optimization iteration, the node coordinates need to be corrected according to the symmetry plane. Furthermore, nodes and symmetry points will be deleted together in one iteration.

2D NUMERICAL EXAMPLES
In this section, the efficiency of the proposed algorithm is verified in a series of dimensionless examples in 2D. Without loss of generality, the initial cross-section of all members is circular with a diameter D 0 = 1. Also, Young's modulus is E = 1.

F I G U R E 3 (a)
The initial configuration of Michell arch problem with the load and supports; (b) the optimized structure in the 15th iteration; (c) the optimal structure in the 50th iteration

Michell arch
The classical Michell arch structure is often used to test the efficiency of the topology optimization method. The truss model has been used to study shape optimization based on the weight minimization standard. (Lamberti & Pappalettere, 2004;D. Wang et al., 2002). Herein, this benchmark example is used to investigate the effectiveness of this nodal-based optimization method that aims to minimize compliance and total volume with a specified nodal number. A downward unit loading P = 1 is imposed at the bottom-middle node as shown in Figure 3a. The initial 43 nodes are nearly homogeneously distributed in a rectangular design domain with a width of 1600 and a height of 400. (Figure 3a). Delaunay triangulation connects nodes in pairs to form 106 frame elements. The pin node at the bottom-left corner is represented in gray, and all three degrees θ = u x = u v = 0 (u = 0). The roller support at the bottom-right corner can move horizontally (u v = 0). Except for the loading, pin, and roller nodes, the rest free nodes can move within the design domain. However, the allowable moving range for an arbitrary free node is a square with a side of 10 for stability reasons. Such a constraint to the design variable has been used in the wellknown SIMP method to avoid abrupt density change in each optimization iteration (Andreassen et al., 2011). The allowable diameter range is 0.2∼2. During the optimization process, the symmetry constraint is activated, making the coordinates of the corresponding node symmetri-F I G U R E 4 The convergence of the normalized objective function, total volume, and compliance for the Michell arch cal about the y-axis. The weighting factors used in the case are w 1 = w 2 = 0.5.
At the beginning of optimization, corner nodes far from the pin position are deleted. The structure gradually transforms into a semi-ellipse shape with 16 nodes (Figure 3b) at the 15th iteration. The upper frame elements are thicker than the other frame elements. The five frame elements connected to the loading point are spread out radially to disperse the load. Then, the upper free nodes gradually approach each other, annihilating themselves to form a thick semi-elliptical arc at the 50th iteration. (Figure 3c). Such an optimal frame structure is similar to weight minimization criteria (Lamberti & Pappalettere, 2004;D. Wang et al., 2002). These optimization results of this type of structure are a closed semi-ellipse or semi-circle geometric shape composed of multiple triangles. If the design domain aspect ratio is close to 2:1, the shape of the optimal structure is close to a semi-circle. The aspect ratio of our design domain is 4:1, so the result is closer to a semi-circle ( Figure 3c). In the final result, the diameters of Frames 1, 2, 3, and 4 are 1.074, 1.294, 1.290, and 1.23, respectively, while the value of Frames 5, 6, and 7 are 1.882, 1.732, and 1.487, respectively, which shows that the frame elements in the upper part are thicker than those connecting to the loading point. The SIMP method also obtained an optimal semi-circle Michell arch structure (Browne, 2013). Since the SIMP method is a continuous topology optimization method, its structure has a smoother arc. The arc part is wider than other positions. Figure 4 indicates that the normalized objective function decreases from 1 to 0.590. Its objective function decreases in the first 10 iterations. The initial structure has more material, which undoubtedly means a stiffer structure. The iterative process resulted in material reduction and an increase in compliance. In addition, node deletion and remeshing cause a sudden change in the stress distribution, making the compliance higher than the initial F I G U R E 5 (a) The initial configuration of Michell arch with 27 nodes and the optimized result in the 20th iteration; (b) the initial configuration of Michell arch with 43 nodes and the optimal result in the 43rd iteration; (c) the initial configuration of Michell arch with 43 nodes and the optimal result in the 45th iteration value. Therefore, in the initial iteration process, its volume is drastically reduced, leading to increased compliance. Since this algorithm includes the remeshing of discrete elements, the value of the normalized objective function inevitably fluctuates and reaches the maximum value of 1.082. Afterward, it reaches the minimum in an oscillatory way, controlled by the size and shape optimization. In this final stage, the nodes move within the allowable range of movement, and the elements adjust the diameter to find a possible better structure. The allowable range is gradually reduced and eventually converges. With the node evolution and remeshing, the nodal number is finally reduced to 8. The normalized total volume is reduced to 0.565, and the normalized compliance drops to 0.615.
The proposed method only takes 50 iterations to reach convergence. In the work of Shojaee et al., a gradient-based method for the size and layout optimization of the truss structures was proposed to optimize the Michelle arch structure (Shojaee et al., 2013). Their algorithm took 60 iterations to approximate the optimal solution. Moreover, their method has fewer initial elements, and the topology of the optimal solution is consistent with the initial configuration. The proposed method is efficient, and if the number of initial structures is further reduced, the convergence efficiency can be further improved (Figure 5a).
The process of deleting nodes changes the topology of the structure so that the result is less dependent on the initial structure. This work began the optimization from three structures, as shown in the left column in Figure 5, to test the dependence of the optimal structure on the initial configurations. Because the initial structure in Figure 5a has fewer nodes to be removed, it therefore can find the optimal point with fewer steps (20 iterations). Although the initial configurations are not the same in Figure 5b-c, the proposed method takes almost the same iterations (43 and 45 iterations) to converge because both initials have 43 nodes. All cases in Figure 5 have the same optimal struc-TA B L E 2 The influence of the weight factor on the optimization result of Michell type structure tures. Therefore, it can conclude that the NEFO algorithm is initial-independent.
Unlike single-objective optimization problems, multiobjective optimization problems often do not have a single optimal solution and only a set of potential solutions called Pareto boundaries. There is no single solution that simultaneously achieves the optimal value for each criterion. The goal of the multi-objective optimization problem is to find a set of representative Pareto optimal solutions. Then, the user selects a solution that meets the needs through the trade-offs between the various criteria. The proposed method uses the linear scalar method to solve the optimization problem and calculate a set of Pareto optimal solutions by weighted summation of these two objectives. In Equation (1), the optimization tendency varies by changing w 1 and w 2 . When w 1 is larger than w 2 , the optimization tends to minimize the structure compliance. If the opposite is the case, optimization is more inclined to minimize the volume of the structure. Table 2 shows the evolution in the optimization structure under different weighting factors. The change of the weighting factors can significantly affect the final structure, and its objective values have also varied accordingly. Finally, an alternative optimal solution set is obtained. As w 1 increases, the upper free nodes gradually shift to the axis of symmetry to obtain a shorter volume. Two weighting factors are set to split the optimization problem into multiple sub-optimization problems to derive the Pareto optimal boundary. The solutions under different weighting factors and the intermediate solutions obtained during the process are drawn as curves (Figure 6).

F I G U R E 6 Pareto optimality for Michell arch structure
The selected optimization needs to obtain better results in the actual structural design according to specific requirements (such as fixed conditions, user needs, economy, and beauty). Figure 6 shows that the NEFO algorithm can adjust the weighting factor to deal with the engineering situation. Compared with other single-objective methods, this makes optimization more flexible and practical.

L-shaped bracket
The classic "L-shaped" bracket structure (Mela, 2014) is also discussed. The boundary conditions and dimensions are shown in Figure 7a. Since the design domain of the L-shaped bracket is non-convex, this design domain is divided into three convex sub-areas. Then, the frame elements are created to avoid the generation of elements outside the design domain. Two gray nodes remain fully fixed (u = 0), and a point load is applied on the right end (P = −1). Delaunay triangulation connects nodes in pairs to form 80 frame elements (Figure 7a). Suppose the nodes with applied boundary conditions have no state change capability. The allowable moving range for an arbitrary free node is a square with a side of 5. The allowable diameter range is 0.2∼2. However, the design domain is an L-shaped non-convex geometry, so the translation of the edge position of the node is partially restricted so that it cannot move to the outside. Thus, 30 free nodes can be deleted according to the rules during the evolution process. The maximum nodal number for this structure is 10. The weighting factors used in the case are w 1 = w 2 = 0.5. With the optimization iteration, its structure gradually changes toward the optimal solution, and the structure of the lower half gradually shows a semi-circle (Figure 7b). The optimal L-shaped bracket (Figure 7c) is very close to the analytical solution (Lewiński & Rozvany, 2008). The optimal structure comprises a quadrilateral and a sector in the analytical solution, and Delaunay triangulation causes F I G U R E 7 (a) The initial configuration of the L-shaped bracket problem with the load and supports; (b) the optimized structure in the 15th iteration; (c) the optimal structure in the 40th iteration the quadrilateral to be cut into three triangles. The load does not affect Frame 1, so its diameter continues to drop to the lower limit of 0.2. The diameters of Frames 2-4 also gradually decrease, similar to the structure obtained by SIMP (Bendsoe & Sigmund, 2013). The normalized objective value decreases from 1 to 0.624 (Figure 8). After the 23rd iteration, the objective function value increased to the maximum value of 1.101. Similar to the previous structural optimization results, since the coordinate points oscillate F I G U R E 8 The convergence of the normalized objective function, total volume, and compliance for the L-shaped bracket between the two values during the final convergence, the objective function fluctuates slightly. Similarly, with the node evolution and partial remeshing, nodes are finally reduced to 10. The normalized total volume is reduced to 0.565, and the normalized total compliance drops to 0.600.

Messerschmidt-Bölkow-Blohm (MBB) beam
The MBB beam has a rectangular design domain of length 2400 and height 400. Due to symmetry, the design domain can be reduced to half (Figure 9a). The algorithm first initializes the entire structure to get the truss-like initial structure (Figure 9a). A point load is applied (P = −1), and the roller boundary conditions are imposed on the right end. The algorithm first generates a series of nodes at equal intervals and then inserts a node in the four nodes of the square, generating 33 nodes. Delaunay triangulation connects nodes in pairs to form 80 frame elements. Suppose the nodes with applied boundary conditions have no state change capability. In contrast, other nodes (i.e., free nodes) can shift in evolution. The allowable moving range for an arbitrary free node is a square with a side of 15. The allowable diameter range is 0.2∼2. The weighting factors used in the case are w 1 = w 2 = 0.5.
The upper right corner nodes are deleted, and an arcshaped structure is gradually formed (Figure 9b). With the node evolution and partial remeshing, the nodal number is finally reduced to 6, as the optimal design for the MMB beam (Andreassen et al., 2011). The optimal MBB structure (Figure 9c) is slightly different from the MBB beam structure obtained by other methods. The MBB beam has an extra horizontal frame element along the mirror plane for the same reason. Frame 1 is deleted in other optimization methods, and its diameter reached the lower limit of 0.2 in the NEFO method. The diameters of Frame 2-4 gradually become thicker in the iterative process, which F I G U R E 9 (a) The initial configuration of the Messerschmidt-Bölkow-Blohm (MBB) beam problem with the load and supports; (b) the optimized structure in the 15th iteration; (c) the optimal structure in the 50th iteration F I G U R E 1 0 The convergence of the normalized objective function, total volume, and compliance for the MBB beam can be observed in the optimized structure of MBB beams obtained by the SIMP method. Figure 10 exhibits that the objective decreases from 1 to 0.604. After the 30th cycle, the value of the objective function increases sharply because of the node removal process. Since this is a discrete optimization, this sudden change is inevitable. Similar to the previous structural optimization results, since the coordinate points oscillate between the two values during the final convergence, the objective function fluctuates slightly. The normalized total volume is reduced to 0.598, and the normalized total compliance drops to 0.611. The 2D numerical experiments are carried out on a PC Intel Core i7 9700K with 3.6 GHz and 32 Gb of memory RAM. The algorithm is coded in MATLAB R2020a. The CPU time of the 2D case shown in Figures 3, 6, and 8 are shown in Table 3. The average CPU time of each iteration of these 2D calculation examples is less than 1 second. Node deletion reduces design variables and improves calculation speed. The "pre-assembled" explicit symbolic expression also improves the efficiency of the code.

The influence of the nodal number limitation
This algorithm can limit the nodal number of the final optimization result. When the nodal number is reduced to less than the preset value, the removal of nodes is stopped, and only the nodal position and member thickness are optimized. This section will discuss the relationship between the nodal number and the final optimization results. Table 4 shows the topology optimization results of the Michell arch under different numbers of nodes. When the TA B L E 5 Optimal structures of MBB beam with a different nodal number nodal number is 8 or 10, the optimal structure is close to other researchers' results (D. Wang et al., 2002). When the nodal number is 10, it has the lowest compliance. Furthermore, the structure is very close to the Michell circular arch beam (Gill et al., 2007). As the nodal number increases, its structure becomes more complex. However, its objective function does not increase too much with its complexity.
The influence of nodal number on the optimal MBB beam structure is summarized in Table 5. When the nodal number is 5 and 6, the algorithm can obtain the results from the SIMP (Andreassen et al., 2011). When the nodal number is 8 or 9, an arc-like structure appears in the upper right area of the MBB beam. Moreover, this arc-like feature also often appears in the optimal structure of MBB beams (Larsen et al., 2018;Martinez et al., 2007).
For the L-shaped bracket (Table 6), when the nodal number is 9 or 10, a result close to the analytical solution can be obtained (Lewiński & Rozvany, 2008). As the nodal number in the final iteration increases, the lower half of the L-shaped bracket gradually becomes a semi-circle. Furthermore, multiple frames connect the semi-circle and the center of the circle. As its shape gradually approaches the analytical solution, its stiffness also increases. The Lshaped bracket compliance with 10 nodes is 0.624, the lowest among several results.
Through the above three sets of examples, it can be found that the increase in the nodal number can increase the geometric complexity of the topology optimization results. However, the relationship between the nodal number and the final objective function value is inconsistent. More nodes sometimes lead to an increase in the objective function due to unnecessary structure. The final structure TA B L E 6 Optimal structures of the L-shaped bracket with a different nodal number stiffness depends on whether the nodal number is reasonable, whether close to the theoretical solution. Especially in the Michell arch and L-shaped bracket, the two structures with the nodal number of 10 show close to the known most solvable design, thus obtaining the best stiffness performance. Therefore, it is significant to set a reasonable limit on the nodal number.

SPACE FRAME OPTIMIZATION
This section shows the application potential of the NEFO method in the actual engineering design field. There have been some optimization methods based on nodal positions and cross-sectional diameters. They optimized the spatial truss element structure, including the boom and the electric tower (Cao et al., 2018;Mathakari et al., 2007). However, most of them use truss elements in this type of research, which does not reflect the actual bar properties. In order to check whether the frame fails due to yielding, von Mises stress of the frame element is introduced, which is where is the axial stress, ℎ is the shear stress, and ,max is the maximum bending stress that occurs at the extreme fiber. Suppose the yield strength of a specific location is greater than the material yield strength, in that case, the structure is considered to be unsafe.

Crane arm
The proposed method is used to design a spatial steel crane arm with an elastic modulus of 210 GPa and a Poisson's ratio of 0.25. Forty-eight nodes are distributed into a 50m cylinder with an equiangular triangle shape and connected to form 10 self-repeated sections with a pattern (Figure 11a). Moreover, their cross-sections are all circles with a radius of 0.05 m. All nodes on the left end are fixed, and a concentrated downward force equivalent to 1000 N is at the middle point of the bottom edge on the right surface. Except for the load and fixed nodes, the remaining nodes can move within a cube area with a side length of 0.25 m. The entire design domain is a cuboid, represented by a dashed line in Figure 11a. The allowable diameter range is 0.01∼0.1 m. During the optimization process, the symmetry constraint is activated. The weighting factors used in the case are w 1 = w 2 = 0.5, and the target nodal number is 35. In Figure 11b, the frame structure near the loading point is in the shape of a small pyramid. In contrast, the frame structures of other sections shrink in the width direction and show good periodicity (Figure 11c). At the four nodes of the fixed end, the four frame elements are not affected by the load, so their diameters are reduced to the lower limit of 0.01 m. The cross-sectional diameter of the frame element close to the load point is larger, with a maximum diameter of 0.06 m (Figure 11b). Nan et al. used a multi-objective evolutionary algorithm to optimize spatial truss arms and obtained similar periodic structures consisting of triangular pyramids (Nan et al., 2020). Stress analysis shows that the maximum von Mises stress is 2.6 × 10 5 Pa, which is much lower than the yield strength of steel (∼200 MPa). The convergence curve in Figure 12 shows that the objective function drops slightly in the early stages of evolution. The first few iterations of NEFO method optimization F I G U R E 1 1 (a) The initial geometry configuration of the crane arm in the perspective view, (b) the optimal design of the crane arm with the same viewpoint (the left-bottom inset denotes the zoomed-in part with maximum stress while the right-top inset is the distribution of member diameter), and (c) the side view of the optimal crane arm F I G U R E 1 2 The convergence of the normalized objective function, total volume, and compliance for the crane arm need to delete unnecessary nodes to reduce design variables and add new frame elements to ensure structural stability, which leads to discontinuities in the objective function. Especially in the iterations from the sixth iteration to the 14th iteration, the deletion of nodes caused the volume of the structure to decrease continuously, which caused fluctuations in compliance. In the 15th iteration, the nodal number reaches the limit, and the optimization process is controlled by shape and size optimization. No more node deletion and element increase, so the objective function gradually decreases steadily. In the 40th to 60th iteration, the objective function, volume, and compliance are gradually reduced and stable, reflecting the good convergence of the NEFO algorithm. The proposed algorithm converges to the optimal solution in 35 iterations.
Compared with the initial structure, the compliance and volume of the optimal tower crane boom have been reduced by nearly 27% and 36%, respectively, which shows that the NEFO algorithm uses less material to get a stronger structure.

Transmission tower
The NEFO method is used to optimize the transmission tower (Cretoff, 2014) to highlight its strengths. Unlike conventional approaches (Couceiro et al., 2016;Shea & Smith, 2006), this work uses frame elements other than truss elements to make the design approach to the actual scenario more closely. All frame elements are steel with an elastic modulus of 210 GPa and a Poisson's ratio of 0.25. And their cross-sections are all circles with a radius of 0.015 m. As shown in Figure 13a, two external loadings, namely, P 1 = 200 N and P 2 = 100 N, are exerted at some specified nodes at which electrical lines are connected. The gravity and wind load are also considered in the optimization. A uniform wind pressure of 21 N/m from the front is adopted for the wind load applied to the tower. To make an unchanged design domain, a cuboid with a size of 1.8 m×14.6 m×13.2 m, the fixed bottom nodes, the loading nodes cannot move in the optimization. Also, the top nodes cannot move in the vertical direction to make the structure has the same height as the initial structure. The allowable diameter range is 0.01∼0.05 m. The weighting factors used in the case are w 1 = w 2 = 0.5, and the symmetry constraint is activated during the optimization process. Figure 13a-d shows the optimal transmission tower. Compared with the original uniform diameter, the outer frame elements are thicker, and the inner frame elements are thinner. For example, their diameters are increased because the frame at positions #1, #2, #3, and #4 needs to bear the load. In particular, the diameters of the frames at position #2 have been increased to the upper limit of 0.05 m. Similar findings have been reported in (de Souza et al., 2016;Mathakari et al., 2007). Compared with the initial structure, the top part becomes thinner. Similar to position #5, the frame becomes thinner due to the small load and lighter weight. Such a geometrical feature is in line F I G U R E 1 3 (a) The optimization result of a transmission tower in the perspective view (the left inset denotes the zoomed-in part with maximum stress while the right inset is the distribution of member diameter), (b) the top view of the optimal tower, (c) the front view of the optimal tower, and (d) the side view of the optimal tower with the load-bearing characteristics of this type of structure. Because of their low strain energy, the NEFO algorithm deletes many nodes in the top parts. Stress analysis is performed on the optimal structure to verify structural safety. The distribution of von Mises stress indicates that high stress is more likely to occur near the junctions of the frame (Figure 13a). Although these joints bear the weight of the entire top, its von Mises stress is still only 73.0 MPa, which is much lower than the yield strength of common steel.
Because of the deletion of some critical nodes, the objective function fluctuates in the optimization process. However, its value is much lower than the initial structure, reducing from 1 to 0.652 in 50 iterations. As shown in Figure 14, it takes 20 iterations to drop the initial 311 nodes to the target of 287 nodes. The compliance and the normalized total volume are 0.654 and 0.655, approximately 34.6% and 34.5% less than the initial configuration. Obtaining lower compliance with less material justifies the effectiveness of the NEFO algorithm for designing engineering applications.

Aquatic dome
The third application example is the aquatic dome structure ( Figure 15a). The design is inspired by the steel dome structure of the London Aquatics Centre (Tarboush & Akçay, 2019). For the sake of simplification, it is assumed F I G U R E 1 4 The convergence of the normalized objective function, total volume, and compliance for the transmission tower design that nodes can only move on bi-layer surfaces defined as z = x 2 /50 − y 2 /200 and z = x 2 /50 -y 2 /200 + 1. A total of 3200 t weight is applied to each node on average to represent the weight of the entire dome. The Delaunay triangulation algorithm obtains the connections of the nodes. All frame elements are steel with an elastic modulus of 210 GPa and a Poisson's ratio of 0.23. Moreover, their crosssections are all circles with a diameter of 0.25 m. In the optimization process, the projection of the movable range of free nodes on the x-y plane is a square with a side length of 0.05 m. In order to ensure symmetry, the mirror plane is activated, and only half the structure needs to be F I G U R E 1 5 (a) The initial geometry configuration of the aquatic dome in the perspective view and (b) its optimization result with the same viewpoint (the left-bottom inset denotes the zoomed-in part with maximum stress while the right-bottom inset is the distribution of member diameter), (c) the front view of the optimal dome, (d) the top view of the optimal dome, and (e) the side view of the optimal dome calculated. The node with point load and the node used for fixation cannot be moved. The initial nodal number is 312, which is limited to 236. The weighting factors used in the case are w 1 = w 2 = 0.5. The allowable diameter range is 0.01∼0.5 m. The concept of node groups is introduced to speed up the optimization process. In the initial situation, nodes on the same surface with the same initial ycoordinate are regarded as a nodal group. The nodes in the node group will be deleted together, and the average strain energy of the node will be replaced with the average strain energy of the node group. Figure 15b-e shows the optimization result. Compared with the initial structure, the steel bars along the x-axis direction are thicker, but the bars in other directions inside have become thinner. This structure is conducive to load transfer, and a similar structure can be seen in the final design of the London Aquatics Centre (Tarboush & Akçay, 2019). The maximum von Mises stress is located at the backside joints, lower than the steel yield strength. The number of elements dropped from 1109 to 941. The compliance and the normalized total volume of the optimal structure are 0.601 and 0.604. The NEFO method can be applied to problems with many nodes and elements by deleting multiple nodes. Table 7 shows the number of iterations, total CPU time, and average CPU time per iteration for spatial examples of the proposed method and the 2D truss example of the GA (Dede et al., 2011). Because of the changed nodal number, only the number of freedoms of the initial TA B L E 7 Comparison of the CPU time between the nodal-based evolutionary frame optimization method and genetic algorithm (GA) iterations and 1052 s to get the optimal truss with 722 freedoms using the GAs with a population size of 80. It is noted that the average CPU time per iteration per freedom is 8.12 × 10 -3 s in the aquatic dome case, slightly lower than that (1.12 × 10 -2 s) of the GA. However, the average CPU time per freedom is 0.52 s, much faster than 1.46 s of the GA. Considering the difference in computer performance and the extra thickness design in the NEFO method, it is plausible to conclude that the proposed method is well competitive with the GA.

CONCLUSION
In this paper, a novel nodal-based evolutionary design optimization algorithm is proposed to minimize the volume and compliance of the frame structure. For the first time, the nodal position and member thickness appear in design variables for frame optimization. Also, the admissible solution space is extended significantly since the nodes can move freely within the design domain. Though the explicitly expressed sensitivity of the objective function is lengthy, it can efficiently get by assembling it in iterations, similar to the global stiffness matrix in finite element analysis. Like the ESO method, the method deletes the most inefficient node in each iteration until it reaches a prescribed limit. Besides, the Delaunay triangulation algorithm reconnects nodes to ensure structural stability in each iteration. The design of three 2D benchmark examples shows that the proposed method can get the optimal configurations within 50 iterations. Besides, it can generate a series of related solutions by adjusting weighting factors. The proposed method can yield elegant buildings by inhibiting nodes moving out from a surface pre-designed by an architect. It can also precisely control the nodal number by adjusting its limitation. Compared with the initial spatial shapes of the crane arm, transmission tower, and aquatic dome, the methods can reduce compliance by approximately 40% with less material (∼60%).
The present work is a preliminary exploration of designing real-world frame structures. In the future, it will include more key engineering factors like the crosssectional shape, safety, reliability, and costs.

A C K N O W L E D G M E N T
This project was supported by the Australian Research Council (DP200102190, FL190100014).