- Email: [email protected]
Efficient force method for the analysis of finite element models comprising of triangular elements using ant colony optimization
Finite Elements in Analysis and Design 45 (2009) 710 -- 720
Contents lists available at ScienceDirect
Finite Elements in Analysis and Design journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / f i n e l
Efficient force method for the analysis of finite element models comprising of triangular elements using ant colony optimization A. Kaveh a,b, ∗ , M. Daei a a b
Department of Civil Engineering, Iran University of Science and Technology, Narmak, Tehran 16, Iran Institute for Mechanics of Materials and Structures, Vienna University of Technology, Karlsplatz 13, A-1040 Wien, Austria
A R T I C L E
I N F O
Article history: Received 27 January 2009 Received in revised form 8 June 2009 Accepted 8 June 2009 Available online 5 July 2009 Keywords: Finite elements Triangular element Force method Null basis Flexibility matrix Sparsity Ant colony system
A B S T R A C T
An efficient algorithm is presented for the formation of null basis of triangular plane stress and plane strain finite element models, corresponding to highly sparse flexibility matrices. This is achieved by applying a modified ant colony system (ACS). An integer linear programming formulation is also presented to evaluate the quality of the results obtained by the proposed ant colony system algorithm. The efficiency of the present algorithm is illustrated through some examples. © 2009 Elsevier B.V. All rights reserved.
1. Introduction The force method of structural analysis, in which the member forces are used as unknowns, is appealing to engineers, since the properties of members of a structure most often depend on the member forces rather than joint displacements. Four different approaches are adopted for the force method of structural analysis, which can be classified as: (1) topological force methods, (2) algebraic force methods, (3) mixed algebraic–topological force methods, and (4) integrated force method. Topological methods have been developed by Henderson and Maunder [1] and Maunder [2] for rigid-jointed skeletal structures using manual selection of cycle bases. Methods suitable for computer programming are due to Kaveh [3,4]. Algebraic methods have been developed by Denke [5], Robinson [6], Topçu [7], Kaneko et al. [8], and mixed algebraic–topological methods have been used by Gilbert and Heath [9], Coleman and Pothen [10,11]. The integrated force method has been developed by Patnaik [12] and Patnaik et al. [13], in which the equilibrium equations and compatibility conditions are satisfied simultaneously in terms of the force variables.
∗ Corresponding author at: Department of Civil Engineering, Iran University of Science and Technology, Narmak, Tehran 16, Iran. Tel.: +98 21 44202710; fax: +98 21 77240398. E-mail address: [email protected] (A. Kaveh). 0168-874X/$ - see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2009.06.005
The force method of structural analysis requires the formation of a maximal set of independent self-equilibrating stress systems (SESs), known as a null basis [14,15]. The elements of this basis form the columns of an m×(S) matrix, B1 , known as the self-stress matrix. The main problem in the application of the force method is the formation of a self-stress matrix corresponding to a sparse flexibility matrix G = Bt1 Fm B1 , where Fm contains the flexibility matrices of the individual members of the structure in a block diagonal form. The combinatorial methods for the force method are very efficient for skeletal structures and, in particular, for rigid-jointed frames. For a general structure, the underlying graph or hypergraph of a SES has not yet been properly defined, and further research is needed. Algebraic methods, on the other hand, can be formulated in a more general form to cover different types of structures such as skeletal structures and finite element models (FEM). The main drawbacks of these methods are the large storage requirements and the higher number of operations. Heuristic algorithms, such as ant colony algorithms, have found many applications in optimization problems in the last decade. The essence of these algorithms lies in the fact that their capability to converge to a good solution does not depend on the specific search space to which they are applied. In this paper, the ant colony system (ACS) which is a variation of the ant colony optimization (ACO) is applied to the formation of null bases of triangular plane stress and plane strain finite element models corresponding to highly sparse and banded flexibility matrices. An integer linear programming formulation is presented to evaluate the quality of the results obtained
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
by the proposed ACS algorithm. The efficiency of the present method is illustrated through simple examples. 2. Formulation of the force method Consider a structure S, which is (S) times statically indeterminate. (S) independent unknown forces are selected as redundants. These unknown forces can be selected from external reactions and or internal forces of the structure. These redundants are denoted by a vector as q = {q1 ,q2 , . . . ,q (S)} t . In order to obtain a statically determinate structure, known as the basic (released or primary) structure of S, the constraints corresponding to redundants are removed. Consider the joint loads vector as p = {p1 ,p2 , . . . ,pn }t , where n is number of entries of the applied nodal load vector and let r denote the mdimensional vector of generalized independent element forces. The equilibrium conditions of the structure can then be expressed as Ar = p
(1)
where A is an n×m equilibrium matrix. The element forces can be written as r = B0 p + B1 q
(2)
where B0 and B1 are rectangular matrices each having m rows, and n and (S) columns, respectively. B1 is called a self-stress matrix as well as null basis matrix. Each column of B1 is known as a null vector. Minimizing the complementary potential energy requires that r minimize the quadratic form, 12 rt Fm r subjected to the constraint of the equilibrium conditions as in Eq. (1). Fm is a m×m block diagonal element flexibility matrix. Using Eq. (2), it can be seen that q must satisfy the following equation: (Bt1 Fm B1 )q = −Bt1 Fm B0 p
(3)
tF
where B1 m B1 = G is the overall flexibility matrix of the structure. Computing the redundant forces q from Eq. (3), r can be found using Eq. (2), i.e. r = [B0 − B1 (BT1 Fm B1 )−1 BT1 Fm B0 ]p
(4)
The structure of the matrix G is again important, and its sparsity, bandwidth and conditioning govern the efficiency of the force method. For the sparsity of G one can search for a sparse B1 matrix, which is often referred to as the sparse null basis problem. 3. Constant stress triangular element For this element, the element forces, Fi = {Fi , Fi , Fi }, are taken
as the natural forces acting along the side of the triangles, as shown in Fig. 1. The corresponding displacements are denoted by i = {i , i , i }. In a global coordinate system, the nodal forces for each element have six components and the nodal forces and element forces can be related by projection. The simple flexibility matrix is as fi =
1 l l tA ct
(5)
where l = {L , L , L }, A is the area of the element, t is the thickness, and
⎡⎡
1 cos2 1 ⎣⎣ 2 ct = cos 1 2G cos2 cos2
⎡ ⎤ 1 cos2 ⎣ 2 ⎦ 1 cos + 1+ 1 1
1 1 1
⎤⎤ 1 1 ⎦⎦ 1
4. Mathematical modeling for optimization problem Since the overall flexibility matrix of the structure G is Bt1 Fm B1 , for the sparsity of G one should select a null basis corresponding to sparse B1 matrix, which is often referred to as the sparse null basis problem. The main objective of this paper is to find sparse selfstress matrices to simplify the solution and to ensure the formation of well-conditioned flexibility matrices. For a SES (null vector), no applied load is required, thus the equilibrium conditions can be expressed as AB1 = 0
(7)
This equation shows that the columns of the matrix A, which is an n×m matrix with rank of n are linearly dependent. There are m−n = t independent columns of B1 which will satisfy this equation, thus forming a set of SESs as a basis. A fact to emphasize is that there are many sets of SESs, which have independent columns and satisfy the above equation. However, the problem is to find a set corresponding to highly sparse B1 matrix. Let us denote the columns of matrix B1 by Si as B1 = [S1 , S2 , . . . , Sg , . . . , St ]
(8)
Suppose the first null vector S1 is found, then it can be normalized by the following equation: et1 S1 = 1
(9)
where e1 = {1 0 . . . 0 . . . 0} is an m×1 unit vector with 1 in the first entry position. The second column S2 can be normalized and must be independent of S1 and these conditions are expressed as et1 S2 = 0
(10)
et2 S2 = 1
(11)
where e2 = {0 1 0 . . . 0 . . . 0} is an m×1 unit vector with 1 in the second entry position. It is obvious that the conditions analogous to these relationships can be formed for the subsequent null vectors. In this section, first the mathematical programming is employed for selecting the column S2 and then extended for the formation of the complete set of the SESs. The first null vector S1 , is arbitrary. Now we find the second null vector S2 , satisfying the following equations: AS2 = 0
(12)
et1 S2 = 0
(13)
et2 S2 = 1
(14)
or more concisely A 0 S2 = 0 e2 I2
(15)
where e2 = {01} is a 2×1 unit vector, with 1 in the gth position which minimizes the function Z = |S2 |. Here, |S2 | denotes the cardinality of S2 and it is equal to the number of non-zero entries of S2 . This can be generalized for the gth null vector Sg , after all the previous null vectors up to g−1 have been obtained. The problem can now be stated as follows: Minimize the objective function of the form Z = |Sg | satisfying
A Ig
(6)
711
0
Sg =
0 eg
(16)
where e¯ g = {0 0 . . . 0 . . . 1} is a g×1 unit vector, with 1 in the gth entry position.
712
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
Fig. 1. A triangular element, natural forces acting along the side of the element.
Fig. 2. Ant technique to find an optimum solution.
By performing the above series of operations, t null vectors B1 = [S1 , S2 , . . . , Sg , . . . , St ] are generated consecutively one after another. It should be noted that for the last null vector (the tth system (t = m−n)), there is no choice. This is because the number of equations is equal to the number of variables, i.e. there are n original equations in the n×m matrix A, and m−n = t orthogonalising equations, thus forming n+t = m equations (with the number of variables being equal to m), leading to a unique solution for the last null vector. A point to notice is the numbering of the members in the structure. The importance of numbering pattern can be recognized by considering the additional equations used in the normalization and orthogonalization. Here, the ant colony system will be applied to choose the members such that the resulting null vectors lead to sparse B1 matrices. 5. Optimization by ant colony systems A meta-heuristic algorithm based on the ants behavior was developed in early 1990s by Dorigo and Gambardella [16]. This algorithm was called ant colony optimization because it was motivated by social behavior of ants. Ant colony system is a variation of the ACO which has proven to behave more robustly and provide far better results for certain problems. In this work, ACS is chosen as a suitable tool for finding sparse null vectors. A brief description of ACO is given in the next section, when describing the process of adapting ACS to the problem of finding sparse null basis. The building blocks of these algorithms are cooperative agents called ants. These agents have simple capabilities, which make their behavior similar to real ants. Real ants are capable of finding the shortest path from food source to their nest or vice versa by smelling pheromones which are chemical substances they leave on the ground while walking. Each ant probabilistically prefers to follow a direction rich in pheromone. Since pheromones do evaporate and lose strength over time, the final result is that more ants tend to pass over the
shortest path and this path is visited more often as the amount of pheromone being laid increases. As an illustrative example, consider the sketch shown in Fig. 2. The number of dashed lines in Fig. 2(c) is approximately proportional to the amount of pheromone deposited by ants. 6. The effect of generator sequence and edge ordering on the sparsity of null basis According to the proposed mathematical modeling, the numbering the members affect the results of the selected null vectors. This can be found out by considering the additional equations which are used in the process of normalization and orthogonalisation. Taking et1 S1 = 1, it is required to have a force equal to unity in the first entry of S1 vector. This entry is called the generator of S1 . Since et1 S2 = 0 and et2 S2 = 1, therefore the first entry in S2 vector, which is the generator of S1 , must be zero, while the second entry must be equal to one. This second entry is known as the generator for the second column in null basis matrix, i.e. the second null vector S2 . Therefore, for the gth null vector, Sg , the forces in the previous generators are zero while in its generator position it is equal to one. As an example, consider a finite element model as shown in Fig. 3. This model is divided into 12 elements and its degree of static indeterminacy is equal to 15. The numbering of the edges for the triangular elements is shown in the discretized structure in Fig. 4. First, the set of generators are chosen as (2 → 3 → 4 → 8 → 9 → 10 → 14 → 15 → 20 → 22 → 26 → 28 → 32 → 7 → 12) . In this set, the first 13 generators are selected from double edges, therefore the resulted null vectors have only two non-zero entries as shown in Fig. 5, and the last two null vectors are illustrated in Fig. 6. In this figure for each null vector, the bold line shows the corresponding generator.
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
713
(2 → 3 → 12 → 8 → 9 → 10 → 14 → 15 → 20 → 22 → 26 → 28 → 32 → 7 → 4) is is used, then the 3rd and 15th null vectors will be as those illustrated in Fig. 8, which correspond to a better null basis than the first one. The remaining 13 null vectors are related to double edges which are unaltered. Therefore, in the next section each sequence of edges as generators is considered as a tour for an ant travel, and the best ant search for the generators is the one leading to the sparsest possible null basis. An efficient algorithm based on the ant colony system is presented in the following section for finding an optimum solution. 7. ACS algorithm for the formation of sparse null basis
Fig. 3. A finite element model with 12 elements and 15 degrees of static indeterminacy.
Fig. 4. The numbering of edges for the elements of the discretized model.
Fig. 5. Typical null vectors resulted for the first 13 generators.
In order to show the effect of choosing different sets of generators, the generator for the 14th SES is replaced by 5 in the previous set. Therefore the sequence of generators is changed to (2 → 3 → 4 → 8 → 9 → 10 → 14 → 15 → 20 → 22 → 26 → 28 → 32 → 5 → 12) in the second attempt. The reason for performing this change is to obtain better results for the 14th and 15th null vectors, while the first 13 null vectors are kept identical to the previous set (Fig. 5). The new results for the 14th and 15th null vectors are shown in Fig. 7. Obviously, this system is sparser than the previous null basis, resulting in a more sparse flexibility matrix. Selection of a different set of generators as well as using a different order of edges in this sequence can alter the results. In the first set of generators, if the position of the third generator (4) is exchanged by the position of last generator (12), i.e. if the sequence
In order to apply the ACO algorithm to a specific problem, it is necessary to represent it as a set of different paths for ants to travel. In the problem of finding sparse null basis, different sequence of generators is considered as a tour for an ant to travel, therefore the cooperative ant agents search to find the best generator sequence resulting in a sparse null basis. Since both the edge numbering, and its order in the generator sequence are important, therefore the pheromone amount is specified by two indices (ij ), where the index i is the generator order in the set of generators, and the index j shows the edge number. As an example, 25 shows the amount of pheromone for selection of the edge number 5 as the 2nd generator in the generators set. In our algorithm first m artificial ants are initially positioned on m edge of elements as primary generators, and then ACS algorithm is applied as follows: An ant k chooses the rth generator by applying the rule of the following equation:
if q ⱕ q0 arg maxu∈Lk (r) ru · ru (17) j= J Otherwise where q is a random number uniformly distributed in [0..1], q0 is a parameter 0 ⱕ q0 ⱕ 1, and J is a random variable selected according to the probability distribution given in the following equation: ⎧ rs ⎪ ⎨ rs if s ∈ Lk (r) k Prs = u∈Lk (r) ru ru (18) ⎪ ⎩ 0 Otherwise Lk (r) is the set of generators that remain to be chosen by ant k as the rth generator and rs is the amount of pheromone deposited on the generator number s as a candidate for being the rth generator. It is assumed that there is an equal amount of pheromone 0 , deposited initially on each generator. rs is the corresponding heuristic value which remains constant throughout the iterations and unlike pheromone amount this is not modified. Moreover, is a parameter controlling the relative importance between and . After an ant chooses one as a generator, the local updating rule on that chosen generator is performed in order to shuffle the solution and prevent focusing on a specific solution. The local updating rule modifies the amount of pheromone by
rs ← (1 − )rs + 0
(19)
where 0 < < 1 is a parameter adjusting the pheromone previously deposited on rs . Once all the ants complete their own tour, the pheromone will be updated for all the edges according to the global updating rule. This pheromone updating is intended to allocate a greater amount of pheromone to shorter tours. The rule is given by the following equation:
rs ← (1 − )rs + rs
(20)
714
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
Fig. 6. The last two null vectors with the edges 7 and 12 as the generators: (a) 14th SES and (b) 15th SES.
Fig. 7. The new result for the last two null vectors with the edge 5 and 12 as their generators: (a) 14th SES and (b) 15th SES.
Fig. 8. The new result for the 3rd and 15th null vectors: (a) 3rd SES and (b) 15th SES.
where
rs =
−1
(Dgb ) 0
if (r, s) ∈ global best tour otherwise
(21)
where Dgb is the sparsity coefficient of the globally best tour (number of non-zero elements in the selected null basis) and 0 < < 1 is the pheromone decay parameter. The best ant tries to find the sparse null basis, and then to reduce the bandwidth corresponding to this null basis, the numbering pattern that is proposed in [17] is applied to the final result. 8. A lower bound for optimal null basis selection Heuristic optimization algorithms like ant colony optimization, seek good feasible solutions to optimization problems in circumstances where the complexities of the problem or the limited time available for solution do not allow exact solution, although these algorithms are often capable of leading to the optimal solutions.
In this section, an integer linear programming formulation is presented to evaluate the efficiency of the suboptimal solution which is obtained by the proposed ACS algorithm. In this modeling, the graph-theoretical approaches developed by Kaveh et al. [17] are applied. In order to transfer the topological property of a finite element model into the connectivity of a graph, the interface graph is used. The interface graph of a FEM, denoted by IG(FEM), is constructed by considering the skeleton of the FEM and adding one edge to the interface edge of each pair of adjacent elements. Thus, a typical overlap of two elements in FEM is represented by double edges in IG(FEM). As an example, the interface graph IG(FEM) of the FEM shown in Fig. 9(a) is illustrated in Fig. 9(b). First, each double edge of the interface graph is selected as the underlying subgraph of a self-stress system. In other words, one of these edges is taken as generator and its corresponding row in the null vector is assigned +1, and for the other edge, we have −1 in the null vector. This process is repeated for all the double edges of the
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
715
Fig. 9. A FEM and the corresponding IG(FEM): (a) a finite element model (FEM) and (b) the corresponding interface graph IG(FEM).
Fig. 10. A finite element model with 24 elements and its disconnected elements.
IG(FEM). Since for the remaining null vectors, the forces in the previous generators should be zero, thus after removing the generating edges of all the double edges from the interface graph, it is changed to a simple graph with no multiple edges. The remaining complementary SESs are chosen from subgraphs of this simple graph by using a graph-theoretical approach as follows: As before, the columns of the matrix B1 is shown by Si , B1 = [S1 , S2 , . . . , Sg , . . . , S ]
(22)
For each column of this matrix, the non-zero entries form a subgraph which corresponds to a null vector. The related subgraph to the kth column of this matrix is denoted by Ck . A sequence of expansions is considered, where in each step one subgraph is selected until the entire simple graph is re-formed as C ∗1 → C ∗2 → C ∗3 → · · · → C ∗(S)
In this sequence, where in the kth step a subgraph Cs is added to C ∗(k−1) , Cs is called admissible if the following equation holds:
(C ∗k ) = (C ∗(k−1) ∪ Cs ) = (C ∗(k−1) ) + 1
(24)
This means that the degree of static indeterminacy must be increased only by unity in each step of the expansion. In what follows, the different constants and variables which are employed in the proposed integer linear programming formulation are defined and then the mathematical model is presented. Consider a resulted simple graph model G = (V,E) with N nodes (vertices) and M members (edges) together with a relation of incidence which associates with each edge a pair of nodes, called its ends, where V = {1, 2, . . . , v, . . . , N}
(23)
E = {1, 2, . . . , e, . . . , M}
(25)
716
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
For finding the kth null vector, which is added to C*(k−1) in the kth step of the expansion process of C ∗1 → C ∗2 → C ∗3 → · · · → C ∗(S) two parameters are defined as follows: ∀i ∈ V : MSG (i) = 1 if edge i is in subgraph C ∗(k−1) and 0 otherwise ∀i ∈ E : NSG (i) = 1 if vertex i is in subgraph C ∗(k−1) and 0 otherwise
The two parameters MSG and NSG , determine those nodes and edges which belong to the previously selected null vectors, i.e. the nodes and edges which belong to the subgraph C*(k−1) . The total number of edges and nodes in this subgraph are denoted by Mt and Nt , respectively. The nodes and edges which belong to the kth null vector, Ck , are determined by Mc and Nc parameters as follows: ∀i ∈ V : Mc (i) = 1 if edge i is in Ck and 0 otherwise ∀i ∈ E : Nc (i) = 1 if vertex i is in Ck and 0 otherwise The following is an integer linear programming formulation for the formation of the kth null vector with minimum non-zero entries. This procedure is repeated for finding all the other null vectors. The first null vector is arbitrary and therefore should be selected as simple as possible. In this formulation, Deg(i) denotes the degree of vertex i and Adj(i,j) is the jth adjacent node of vertex i min :
M
Mc (i)
(26)
i=1
s.t. : [A][B1 ] = 0
(27)
−U × Mc (i) ⱕ S(i) ⱕ U × Mc (i) Nc (i) ⱕ
Deg(i)
∀i ∈ E
(28)
Mc (Adj(i, j)) ⱕ Deg(i) × Nc (i)
∀i ∈ V
(29)
Fig. 11. Pattern of the null basis matrix resulted by the best ant.
j=1
Mt +
M
⎛ Mc (i) [1 − MSG (i)] − 2 ⎝Nt +
i=1
i=1
Mc (i) − 2
⎞ Nc (i) [1 − NSG (i)]⎠
9. Numerical results
i=1
+3=k M
N
(30) N
Nc (i) + 3 = 1
(31)
i=1
The equilibrium condition is expressed by the constraint (27). According to the constraint (28), if S(i), which shows the force in edge i, is non-zero, then its corresponding parameter in the subgraph, Mc (i), will be equal to 1. In this formulation, U is an upper bound for the forces in the structure, which is assumed as a big number. Constraint (29) ensures that if a node is in the selected SES, then its corresponding parameter in the subgraph, Nc (i), will be equal to 1. The independency control is insured by constraint (30), and constraint (31) makes it certain that each selected SES has a static indeterminacy of unit value. Based on this formulation, for a resulted simple graph (the graph with the generators of the double edges being removed) with M edges and N nodes, there is (2M+4N+2) constraints and (2M+N) variables, and the computational time obtained with LINGO solver, is extremely high for large models. However, in this paper the result of this method is used only to evaluate the quality of the solutions which are obtained by the proposed ACS algorithm.
In this section, four examples are presented to illustrate the performance of the proposed ACS algorithm, and to provide a measure of its efficiency. This algorithm is coded by MATLAB, and is run on a personal computer Pentium䉸 4 CPU 3.40 GHz, 1.00 GB of RAM. Example 1. Consider a finite element model as shown in Fig. 10(a). The numbering of the edges for the elements of the model is shown on the disconnected elements in Fig. 10(b). This model is divided into 24 elements, and its degree of static indeterminacy is equal to 35. Twenty-nine generators are selected from the double edges, having the null vectors with only two non-zero entries. For the remaining 6 null vectors, the set of generators is chosen by the best ant as (51→7→60→66→13→23). The total number of non-zero entries of the null basis matrix is 130, having the pattern as shown in Fig. 11. Example 2. A finite element model with two cut-outs is considered as shown in Fig. 12. This model consists of 36 elements and has 47 degrees of static indeterminacy. There are 39 double edges which result in 39 null vectors with two non-zero entries. The other 8 complementary null vectors are
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
717
Fig. 12. A finite element model with two cut-outs.
Fig. 13. Eight complementary null vectors selected by ACS.
illustrated in Fig. 13. For each SES, the corresponding generator is shown in bold line. The total number of non-zero entries of the null basis matrix is 222, having the pattern as shown in Fig. 14. The elapsed run time is less than 1 s. The number of non-zero entries of the null basis matrix for the LU factorization method, the graph-theoretical method [17] and the present approach are depicted in Table 1 for comparison.
Example 3. Consider a beam model, divided into 400 triangular elements as shown in Fig. 15. This beam has the following properties: thickness = 0.05 m,
E = 2e + 8 kN/m2 ,
= 0.3,
DSI = 711.
The proposed algorithm selected 555 null vectors with 2 entries (double edges), and 156 complementary null vectors with 12 entries with a typical pattern shown in Fig. 16.
718
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
Fig. 16. Typical pattern of the selected null vector with 12 entries.
Fig. 14. The pattern of null basis matrix resulted by the best ant.
Table 1 Comparison of the results. Method
Number of non-zero entries
Present method Graph-theoretical method LU factorization method
222 234 461 Fig. 17. Pattern of the null basis matrix resulted by the best ant.
Fig. 15. A beam model with 400 triangular elements.
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
719
Fig. 19. An H-shaped FEM with 312 triangular elements.
Fig. 18. Pattern of the resulted flexibility matrix.
Table 2 Comparison of the results. Method
Number of non-zero entries
CPU time (s)
Present method Graph-theoretical method LU factorization method
2982 2982 28,632
23.041 14.547 28.094
For this sparse basis, the total number of non-zero entries is 2982. The pattern of sparse null basis and the resulted flexibility matrix are shown in Figs. 17 and 18, respectively. The CPU time and the number of non-zero entries for the LU factorization method, the graphtheoretical method [17] and the present approach are depicted in Table 2 for comparison. For this sparse basis, the total number of non-zero entries is 2982. The pattern of sparse null basis and the resulted flexibility matrix are shown in Figs. 17 and 18, respectively. Example 4. An H-shaped FEM with 312 triangular elements is considered as shown in Fig. 19. This model has the same properties as that of Example 3. The present method leads to the formation of 430 null vectors with two entries (double edges) and 119 complementary null vectors with 12 entries, corresponding to highly sparse and banded flexibility matrix as shown in Fig. 20. The results obtained for all 4 examples are identical to those obtained by the integer linear programming formulation of Section 8. However, the only difference is the computational time which is much lower for the ACS algorithm. 10. Conclusions In this paper, an ant colony system is developed for the formation of sparse null basis leading to sparse self-stress matrices, and correspondingly highly sparse flexibility matrices for triangular plane stress and plane strain finite element models. An integer linear programming formulation is also presented to evaluate the quality of the solutions obtained by the proposed ACS algorithm. The method
Fig. 20. Pattern of the resulted flexibility matrix.
presented in this paper is not confined to plane stress and plane strain triangular elements and it can easily be extended to other finite elements. Acknowledgement The first author is grateful to Iran National Science Foundation for the support. References [1] J.C.de.C. Henderson, E.A.W. Maunder, A problem in applied topology, J. Inst. Math. Appl. 5 (1969) 254–269. [2] E.A.W. Maunder, Topological and linear analysis of skeletal structures, Ph.D. Thesis, London University, Imperial College, March 1971.
720
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
[3] A. Kaveh, Application of topology and matroid theory to the analysis of structures, Ph.D. Thesis, London University, Imperial College, March 1974. [4] A. Kaveh, Improved cycle bases for the flexibility analysis of structures, Comput. Methods Appl. Mech. Eng. 9 (1976) 267–272. [5] P.H. Denke, A general digital computer analysis of statically indeterminate structures, NASA-TD-D-1666, 1962. [6] J. Robinson, Integrated Theory of Finite Element Methods, Wiley, New York, 1973. [7] A. Topçu, A contribution to the systematic analysis of finite element structures using the force method, Doctoral Dissertation, Essen University, 1979 (in German). [8] I. Kaneko, M. Lawo, G. Thierauf, On computational procedures for the force methods, Int. J. Numer. Methods Eng. 18 (1982) 1469–1495. [9] J.R. Gilbert, M.T. Heath, Computing a sparse basis for the null space, SIAM J. Alg. Disc. Methods 8 (1987) 446–459. [10] T.F. Coleman, A. Pothen, The null space problem I; complexity, SIAM J. Alg. Disc. Methods 7 (4) (1986) 527–537.
[11] T.F. Coleman, A. Pothen, The null space problem II; Algorithms, SIAM J. Alg. Disc. Methods 8 (4) (1987) 544–561. [12] S.N. Patnaik, An integrated force method for discrete analysis, Int. J. Numer. Methods Eng. 6 (2) (1973) 237–251. [13] S.N. Patnaik, L. Berke, R.H. Gallagher, Integrated force method versus displacement method for finite element analysis, Comput. Struct. 38 (4) (1991) 377–407. [14] A. Kaveh, A combinatorial optimization problem; optimal generalized cycle bases, Comput. Methods Appl. Mech. Eng. 20 (3) (1979) 9–52. [15] A. Kaveh, Structural Mechanics: Graph and Matrix Methods, third ed, Research Studies Press (Wiley), Somerset, UK, 2004. [16] M. Dorigo, L.M. Gambardella, Ant colony system: a cooperative learning approach to the traveling salesman problem, IEEE Trans. Evol. Comput. 1 (1) (1997) 53–66. [17] A. Kaveh, K. Koohestani, N. Taghizadieh, Efficient finite element analysis by graph-theoretical force method, Finite Elem. Anal. Des. 43 (2007) 543–554.
Contents lists available at ScienceDirect
Finite Elements in Analysis and Design journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / f i n e l
Efficient force method for the analysis of finite element models comprising of triangular elements using ant colony optimization A. Kaveh a,b, ∗ , M. Daei a a b
Department of Civil Engineering, Iran University of Science and Technology, Narmak, Tehran 16, Iran Institute for Mechanics of Materials and Structures, Vienna University of Technology, Karlsplatz 13, A-1040 Wien, Austria
A R T I C L E
I N F O
Article history: Received 27 January 2009 Received in revised form 8 June 2009 Accepted 8 June 2009 Available online 5 July 2009 Keywords: Finite elements Triangular element Force method Null basis Flexibility matrix Sparsity Ant colony system
A B S T R A C T
An efficient algorithm is presented for the formation of null basis of triangular plane stress and plane strain finite element models, corresponding to highly sparse flexibility matrices. This is achieved by applying a modified ant colony system (ACS). An integer linear programming formulation is also presented to evaluate the quality of the results obtained by the proposed ant colony system algorithm. The efficiency of the present algorithm is illustrated through some examples. © 2009 Elsevier B.V. All rights reserved.
1. Introduction The force method of structural analysis, in which the member forces are used as unknowns, is appealing to engineers, since the properties of members of a structure most often depend on the member forces rather than joint displacements. Four different approaches are adopted for the force method of structural analysis, which can be classified as: (1) topological force methods, (2) algebraic force methods, (3) mixed algebraic–topological force methods, and (4) integrated force method. Topological methods have been developed by Henderson and Maunder [1] and Maunder [2] for rigid-jointed skeletal structures using manual selection of cycle bases. Methods suitable for computer programming are due to Kaveh [3,4]. Algebraic methods have been developed by Denke [5], Robinson [6], Topçu [7], Kaneko et al. [8], and mixed algebraic–topological methods have been used by Gilbert and Heath [9], Coleman and Pothen [10,11]. The integrated force method has been developed by Patnaik [12] and Patnaik et al. [13], in which the equilibrium equations and compatibility conditions are satisfied simultaneously in terms of the force variables.
∗ Corresponding author at: Department of Civil Engineering, Iran University of Science and Technology, Narmak, Tehran 16, Iran. Tel.: +98 21 44202710; fax: +98 21 77240398. E-mail address: [email protected] (A. Kaveh). 0168-874X/$ - see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2009.06.005
The force method of structural analysis requires the formation of a maximal set of independent self-equilibrating stress systems (SESs), known as a null basis [14,15]. The elements of this basis form the columns of an m×(S) matrix, B1 , known as the self-stress matrix. The main problem in the application of the force method is the formation of a self-stress matrix corresponding to a sparse flexibility matrix G = Bt1 Fm B1 , where Fm contains the flexibility matrices of the individual members of the structure in a block diagonal form. The combinatorial methods for the force method are very efficient for skeletal structures and, in particular, for rigid-jointed frames. For a general structure, the underlying graph or hypergraph of a SES has not yet been properly defined, and further research is needed. Algebraic methods, on the other hand, can be formulated in a more general form to cover different types of structures such as skeletal structures and finite element models (FEM). The main drawbacks of these methods are the large storage requirements and the higher number of operations. Heuristic algorithms, such as ant colony algorithms, have found many applications in optimization problems in the last decade. The essence of these algorithms lies in the fact that their capability to converge to a good solution does not depend on the specific search space to which they are applied. In this paper, the ant colony system (ACS) which is a variation of the ant colony optimization (ACO) is applied to the formation of null bases of triangular plane stress and plane strain finite element models corresponding to highly sparse and banded flexibility matrices. An integer linear programming formulation is presented to evaluate the quality of the results obtained
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
by the proposed ACS algorithm. The efficiency of the present method is illustrated through simple examples. 2. Formulation of the force method Consider a structure S, which is (S) times statically indeterminate. (S) independent unknown forces are selected as redundants. These unknown forces can be selected from external reactions and or internal forces of the structure. These redundants are denoted by a vector as q = {q1 ,q2 , . . . ,q (S)} t . In order to obtain a statically determinate structure, known as the basic (released or primary) structure of S, the constraints corresponding to redundants are removed. Consider the joint loads vector as p = {p1 ,p2 , . . . ,pn }t , where n is number of entries of the applied nodal load vector and let r denote the mdimensional vector of generalized independent element forces. The equilibrium conditions of the structure can then be expressed as Ar = p
(1)
where A is an n×m equilibrium matrix. The element forces can be written as r = B0 p + B1 q
(2)
where B0 and B1 are rectangular matrices each having m rows, and n and (S) columns, respectively. B1 is called a self-stress matrix as well as null basis matrix. Each column of B1 is known as a null vector. Minimizing the complementary potential energy requires that r minimize the quadratic form, 12 rt Fm r subjected to the constraint of the equilibrium conditions as in Eq. (1). Fm is a m×m block diagonal element flexibility matrix. Using Eq. (2), it can be seen that q must satisfy the following equation: (Bt1 Fm B1 )q = −Bt1 Fm B0 p
(3)
tF
where B1 m B1 = G is the overall flexibility matrix of the structure. Computing the redundant forces q from Eq. (3), r can be found using Eq. (2), i.e. r = [B0 − B1 (BT1 Fm B1 )−1 BT1 Fm B0 ]p
(4)
The structure of the matrix G is again important, and its sparsity, bandwidth and conditioning govern the efficiency of the force method. For the sparsity of G one can search for a sparse B1 matrix, which is often referred to as the sparse null basis problem. 3. Constant stress triangular element For this element, the element forces, Fi = {Fi , Fi , Fi }, are taken
as the natural forces acting along the side of the triangles, as shown in Fig. 1. The corresponding displacements are denoted by i = {i , i , i }. In a global coordinate system, the nodal forces for each element have six components and the nodal forces and element forces can be related by projection. The simple flexibility matrix is as fi =
1 l l tA ct
(5)
where l = {L , L , L }, A is the area of the element, t is the thickness, and
⎡⎡
1 cos2 1 ⎣⎣ 2 ct = cos 1 2G cos2 cos2
⎡ ⎤ 1 cos2 ⎣ 2 ⎦ 1 cos + 1+ 1 1
1 1 1
⎤⎤ 1 1 ⎦⎦ 1
4. Mathematical modeling for optimization problem Since the overall flexibility matrix of the structure G is Bt1 Fm B1 , for the sparsity of G one should select a null basis corresponding to sparse B1 matrix, which is often referred to as the sparse null basis problem. The main objective of this paper is to find sparse selfstress matrices to simplify the solution and to ensure the formation of well-conditioned flexibility matrices. For a SES (null vector), no applied load is required, thus the equilibrium conditions can be expressed as AB1 = 0
(7)
This equation shows that the columns of the matrix A, which is an n×m matrix with rank of n are linearly dependent. There are m−n = t independent columns of B1 which will satisfy this equation, thus forming a set of SESs as a basis. A fact to emphasize is that there are many sets of SESs, which have independent columns and satisfy the above equation. However, the problem is to find a set corresponding to highly sparse B1 matrix. Let us denote the columns of matrix B1 by Si as B1 = [S1 , S2 , . . . , Sg , . . . , St ]
(8)
Suppose the first null vector S1 is found, then it can be normalized by the following equation: et1 S1 = 1
(9)
where e1 = {1 0 . . . 0 . . . 0} is an m×1 unit vector with 1 in the first entry position. The second column S2 can be normalized and must be independent of S1 and these conditions are expressed as et1 S2 = 0
(10)
et2 S2 = 1
(11)
where e2 = {0 1 0 . . . 0 . . . 0} is an m×1 unit vector with 1 in the second entry position. It is obvious that the conditions analogous to these relationships can be formed for the subsequent null vectors. In this section, first the mathematical programming is employed for selecting the column S2 and then extended for the formation of the complete set of the SESs. The first null vector S1 , is arbitrary. Now we find the second null vector S2 , satisfying the following equations: AS2 = 0
(12)
et1 S2 = 0
(13)
et2 S2 = 1
(14)
or more concisely A 0 S2 = 0 e2 I2
(15)
where e2 = {01} is a 2×1 unit vector, with 1 in the gth position which minimizes the function Z = |S2 |. Here, |S2 | denotes the cardinality of S2 and it is equal to the number of non-zero entries of S2 . This can be generalized for the gth null vector Sg , after all the previous null vectors up to g−1 have been obtained. The problem can now be stated as follows: Minimize the objective function of the form Z = |Sg | satisfying
A Ig
(6)
711
0
Sg =
0 eg
(16)
where e¯ g = {0 0 . . . 0 . . . 1} is a g×1 unit vector, with 1 in the gth entry position.
712
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
Fig. 1. A triangular element, natural forces acting along the side of the element.
Fig. 2. Ant technique to find an optimum solution.
By performing the above series of operations, t null vectors B1 = [S1 , S2 , . . . , Sg , . . . , St ] are generated consecutively one after another. It should be noted that for the last null vector (the tth system (t = m−n)), there is no choice. This is because the number of equations is equal to the number of variables, i.e. there are n original equations in the n×m matrix A, and m−n = t orthogonalising equations, thus forming n+t = m equations (with the number of variables being equal to m), leading to a unique solution for the last null vector. A point to notice is the numbering of the members in the structure. The importance of numbering pattern can be recognized by considering the additional equations used in the normalization and orthogonalization. Here, the ant colony system will be applied to choose the members such that the resulting null vectors lead to sparse B1 matrices. 5. Optimization by ant colony systems A meta-heuristic algorithm based on the ants behavior was developed in early 1990s by Dorigo and Gambardella [16]. This algorithm was called ant colony optimization because it was motivated by social behavior of ants. Ant colony system is a variation of the ACO which has proven to behave more robustly and provide far better results for certain problems. In this work, ACS is chosen as a suitable tool for finding sparse null vectors. A brief description of ACO is given in the next section, when describing the process of adapting ACS to the problem of finding sparse null basis. The building blocks of these algorithms are cooperative agents called ants. These agents have simple capabilities, which make their behavior similar to real ants. Real ants are capable of finding the shortest path from food source to their nest or vice versa by smelling pheromones which are chemical substances they leave on the ground while walking. Each ant probabilistically prefers to follow a direction rich in pheromone. Since pheromones do evaporate and lose strength over time, the final result is that more ants tend to pass over the
shortest path and this path is visited more often as the amount of pheromone being laid increases. As an illustrative example, consider the sketch shown in Fig. 2. The number of dashed lines in Fig. 2(c) is approximately proportional to the amount of pheromone deposited by ants. 6. The effect of generator sequence and edge ordering on the sparsity of null basis According to the proposed mathematical modeling, the numbering the members affect the results of the selected null vectors. This can be found out by considering the additional equations which are used in the process of normalization and orthogonalisation. Taking et1 S1 = 1, it is required to have a force equal to unity in the first entry of S1 vector. This entry is called the generator of S1 . Since et1 S2 = 0 and et2 S2 = 1, therefore the first entry in S2 vector, which is the generator of S1 , must be zero, while the second entry must be equal to one. This second entry is known as the generator for the second column in null basis matrix, i.e. the second null vector S2 . Therefore, for the gth null vector, Sg , the forces in the previous generators are zero while in its generator position it is equal to one. As an example, consider a finite element model as shown in Fig. 3. This model is divided into 12 elements and its degree of static indeterminacy is equal to 15. The numbering of the edges for the triangular elements is shown in the discretized structure in Fig. 4. First, the set of generators are chosen as (2 → 3 → 4 → 8 → 9 → 10 → 14 → 15 → 20 → 22 → 26 → 28 → 32 → 7 → 12) . In this set, the first 13 generators are selected from double edges, therefore the resulted null vectors have only two non-zero entries as shown in Fig. 5, and the last two null vectors are illustrated in Fig. 6. In this figure for each null vector, the bold line shows the corresponding generator.
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
713
(2 → 3 → 12 → 8 → 9 → 10 → 14 → 15 → 20 → 22 → 26 → 28 → 32 → 7 → 4) is is used, then the 3rd and 15th null vectors will be as those illustrated in Fig. 8, which correspond to a better null basis than the first one. The remaining 13 null vectors are related to double edges which are unaltered. Therefore, in the next section each sequence of edges as generators is considered as a tour for an ant travel, and the best ant search for the generators is the one leading to the sparsest possible null basis. An efficient algorithm based on the ant colony system is presented in the following section for finding an optimum solution. 7. ACS algorithm for the formation of sparse null basis
Fig. 3. A finite element model with 12 elements and 15 degrees of static indeterminacy.
Fig. 4. The numbering of edges for the elements of the discretized model.
Fig. 5. Typical null vectors resulted for the first 13 generators.
In order to show the effect of choosing different sets of generators, the generator for the 14th SES is replaced by 5 in the previous set. Therefore the sequence of generators is changed to (2 → 3 → 4 → 8 → 9 → 10 → 14 → 15 → 20 → 22 → 26 → 28 → 32 → 5 → 12) in the second attempt. The reason for performing this change is to obtain better results for the 14th and 15th null vectors, while the first 13 null vectors are kept identical to the previous set (Fig. 5). The new results for the 14th and 15th null vectors are shown in Fig. 7. Obviously, this system is sparser than the previous null basis, resulting in a more sparse flexibility matrix. Selection of a different set of generators as well as using a different order of edges in this sequence can alter the results. In the first set of generators, if the position of the third generator (4) is exchanged by the position of last generator (12), i.e. if the sequence
In order to apply the ACO algorithm to a specific problem, it is necessary to represent it as a set of different paths for ants to travel. In the problem of finding sparse null basis, different sequence of generators is considered as a tour for an ant to travel, therefore the cooperative ant agents search to find the best generator sequence resulting in a sparse null basis. Since both the edge numbering, and its order in the generator sequence are important, therefore the pheromone amount is specified by two indices (ij ), where the index i is the generator order in the set of generators, and the index j shows the edge number. As an example, 25 shows the amount of pheromone for selection of the edge number 5 as the 2nd generator in the generators set. In our algorithm first m artificial ants are initially positioned on m edge of elements as primary generators, and then ACS algorithm is applied as follows: An ant k chooses the rth generator by applying the rule of the following equation:
if q ⱕ q0 arg maxu∈Lk (r) ru · ru (17) j= J Otherwise where q is a random number uniformly distributed in [0..1], q0 is a parameter 0 ⱕ q0 ⱕ 1, and J is a random variable selected according to the probability distribution given in the following equation: ⎧ rs ⎪ ⎨ rs if s ∈ Lk (r) k Prs = u∈Lk (r) ru ru (18) ⎪ ⎩ 0 Otherwise Lk (r) is the set of generators that remain to be chosen by ant k as the rth generator and rs is the amount of pheromone deposited on the generator number s as a candidate for being the rth generator. It is assumed that there is an equal amount of pheromone 0 , deposited initially on each generator. rs is the corresponding heuristic value which remains constant throughout the iterations and unlike pheromone amount this is not modified. Moreover, is a parameter controlling the relative importance between and . After an ant chooses one as a generator, the local updating rule on that chosen generator is performed in order to shuffle the solution and prevent focusing on a specific solution. The local updating rule modifies the amount of pheromone by
rs ← (1 − )rs + 0
(19)
where 0 < < 1 is a parameter adjusting the pheromone previously deposited on rs . Once all the ants complete their own tour, the pheromone will be updated for all the edges according to the global updating rule. This pheromone updating is intended to allocate a greater amount of pheromone to shorter tours. The rule is given by the following equation:
rs ← (1 − )rs + rs
(20)
714
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
Fig. 6. The last two null vectors with the edges 7 and 12 as the generators: (a) 14th SES and (b) 15th SES.
Fig. 7. The new result for the last two null vectors with the edge 5 and 12 as their generators: (a) 14th SES and (b) 15th SES.
Fig. 8. The new result for the 3rd and 15th null vectors: (a) 3rd SES and (b) 15th SES.
where
rs =
−1
(Dgb ) 0
if (r, s) ∈ global best tour otherwise
(21)
where Dgb is the sparsity coefficient of the globally best tour (number of non-zero elements in the selected null basis) and 0 < < 1 is the pheromone decay parameter. The best ant tries to find the sparse null basis, and then to reduce the bandwidth corresponding to this null basis, the numbering pattern that is proposed in [17] is applied to the final result. 8. A lower bound for optimal null basis selection Heuristic optimization algorithms like ant colony optimization, seek good feasible solutions to optimization problems in circumstances where the complexities of the problem or the limited time available for solution do not allow exact solution, although these algorithms are often capable of leading to the optimal solutions.
In this section, an integer linear programming formulation is presented to evaluate the efficiency of the suboptimal solution which is obtained by the proposed ACS algorithm. In this modeling, the graph-theoretical approaches developed by Kaveh et al. [17] are applied. In order to transfer the topological property of a finite element model into the connectivity of a graph, the interface graph is used. The interface graph of a FEM, denoted by IG(FEM), is constructed by considering the skeleton of the FEM and adding one edge to the interface edge of each pair of adjacent elements. Thus, a typical overlap of two elements in FEM is represented by double edges in IG(FEM). As an example, the interface graph IG(FEM) of the FEM shown in Fig. 9(a) is illustrated in Fig. 9(b). First, each double edge of the interface graph is selected as the underlying subgraph of a self-stress system. In other words, one of these edges is taken as generator and its corresponding row in the null vector is assigned +1, and for the other edge, we have −1 in the null vector. This process is repeated for all the double edges of the
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
715
Fig. 9. A FEM and the corresponding IG(FEM): (a) a finite element model (FEM) and (b) the corresponding interface graph IG(FEM).
Fig. 10. A finite element model with 24 elements and its disconnected elements.
IG(FEM). Since for the remaining null vectors, the forces in the previous generators should be zero, thus after removing the generating edges of all the double edges from the interface graph, it is changed to a simple graph with no multiple edges. The remaining complementary SESs are chosen from subgraphs of this simple graph by using a graph-theoretical approach as follows: As before, the columns of the matrix B1 is shown by Si , B1 = [S1 , S2 , . . . , Sg , . . . , S ]
(22)
For each column of this matrix, the non-zero entries form a subgraph which corresponds to a null vector. The related subgraph to the kth column of this matrix is denoted by Ck . A sequence of expansions is considered, where in each step one subgraph is selected until the entire simple graph is re-formed as C ∗1 → C ∗2 → C ∗3 → · · · → C ∗(S)
In this sequence, where in the kth step a subgraph Cs is added to C ∗(k−1) , Cs is called admissible if the following equation holds:
(C ∗k ) = (C ∗(k−1) ∪ Cs ) = (C ∗(k−1) ) + 1
(24)
This means that the degree of static indeterminacy must be increased only by unity in each step of the expansion. In what follows, the different constants and variables which are employed in the proposed integer linear programming formulation are defined and then the mathematical model is presented. Consider a resulted simple graph model G = (V,E) with N nodes (vertices) and M members (edges) together with a relation of incidence which associates with each edge a pair of nodes, called its ends, where V = {1, 2, . . . , v, . . . , N}
(23)
E = {1, 2, . . . , e, . . . , M}
(25)
716
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
For finding the kth null vector, which is added to C*(k−1) in the kth step of the expansion process of C ∗1 → C ∗2 → C ∗3 → · · · → C ∗(S) two parameters are defined as follows: ∀i ∈ V : MSG (i) = 1 if edge i is in subgraph C ∗(k−1) and 0 otherwise ∀i ∈ E : NSG (i) = 1 if vertex i is in subgraph C ∗(k−1) and 0 otherwise
The two parameters MSG and NSG , determine those nodes and edges which belong to the previously selected null vectors, i.e. the nodes and edges which belong to the subgraph C*(k−1) . The total number of edges and nodes in this subgraph are denoted by Mt and Nt , respectively. The nodes and edges which belong to the kth null vector, Ck , are determined by Mc and Nc parameters as follows: ∀i ∈ V : Mc (i) = 1 if edge i is in Ck and 0 otherwise ∀i ∈ E : Nc (i) = 1 if vertex i is in Ck and 0 otherwise The following is an integer linear programming formulation for the formation of the kth null vector with minimum non-zero entries. This procedure is repeated for finding all the other null vectors. The first null vector is arbitrary and therefore should be selected as simple as possible. In this formulation, Deg(i) denotes the degree of vertex i and Adj(i,j) is the jth adjacent node of vertex i min :
M
Mc (i)
(26)
i=1
s.t. : [A][B1 ] = 0
(27)
−U × Mc (i) ⱕ S(i) ⱕ U × Mc (i) Nc (i) ⱕ
Deg(i)
∀i ∈ E
(28)
Mc (Adj(i, j)) ⱕ Deg(i) × Nc (i)
∀i ∈ V
(29)
Fig. 11. Pattern of the null basis matrix resulted by the best ant.
j=1
Mt +
M
⎛ Mc (i) [1 − MSG (i)] − 2 ⎝Nt +
i=1
i=1
Mc (i) − 2
⎞ Nc (i) [1 − NSG (i)]⎠
9. Numerical results
i=1
+3=k M
N
(30) N
Nc (i) + 3 = 1
(31)
i=1
The equilibrium condition is expressed by the constraint (27). According to the constraint (28), if S(i), which shows the force in edge i, is non-zero, then its corresponding parameter in the subgraph, Mc (i), will be equal to 1. In this formulation, U is an upper bound for the forces in the structure, which is assumed as a big number. Constraint (29) ensures that if a node is in the selected SES, then its corresponding parameter in the subgraph, Nc (i), will be equal to 1. The independency control is insured by constraint (30), and constraint (31) makes it certain that each selected SES has a static indeterminacy of unit value. Based on this formulation, for a resulted simple graph (the graph with the generators of the double edges being removed) with M edges and N nodes, there is (2M+4N+2) constraints and (2M+N) variables, and the computational time obtained with LINGO solver, is extremely high for large models. However, in this paper the result of this method is used only to evaluate the quality of the solutions which are obtained by the proposed ACS algorithm.
In this section, four examples are presented to illustrate the performance of the proposed ACS algorithm, and to provide a measure of its efficiency. This algorithm is coded by MATLAB, and is run on a personal computer Pentium䉸 4 CPU 3.40 GHz, 1.00 GB of RAM. Example 1. Consider a finite element model as shown in Fig. 10(a). The numbering of the edges for the elements of the model is shown on the disconnected elements in Fig. 10(b). This model is divided into 24 elements, and its degree of static indeterminacy is equal to 35. Twenty-nine generators are selected from the double edges, having the null vectors with only two non-zero entries. For the remaining 6 null vectors, the set of generators is chosen by the best ant as (51→7→60→66→13→23). The total number of non-zero entries of the null basis matrix is 130, having the pattern as shown in Fig. 11. Example 2. A finite element model with two cut-outs is considered as shown in Fig. 12. This model consists of 36 elements and has 47 degrees of static indeterminacy. There are 39 double edges which result in 39 null vectors with two non-zero entries. The other 8 complementary null vectors are
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
717
Fig. 12. A finite element model with two cut-outs.
Fig. 13. Eight complementary null vectors selected by ACS.
illustrated in Fig. 13. For each SES, the corresponding generator is shown in bold line. The total number of non-zero entries of the null basis matrix is 222, having the pattern as shown in Fig. 14. The elapsed run time is less than 1 s. The number of non-zero entries of the null basis matrix for the LU factorization method, the graph-theoretical method [17] and the present approach are depicted in Table 1 for comparison.
Example 3. Consider a beam model, divided into 400 triangular elements as shown in Fig. 15. This beam has the following properties: thickness = 0.05 m,
E = 2e + 8 kN/m2 ,
= 0.3,
DSI = 711.
The proposed algorithm selected 555 null vectors with 2 entries (double edges), and 156 complementary null vectors with 12 entries with a typical pattern shown in Fig. 16.
718
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
Fig. 16. Typical pattern of the selected null vector with 12 entries.
Fig. 14. The pattern of null basis matrix resulted by the best ant.
Table 1 Comparison of the results. Method
Number of non-zero entries
Present method Graph-theoretical method LU factorization method
222 234 461 Fig. 17. Pattern of the null basis matrix resulted by the best ant.
Fig. 15. A beam model with 400 triangular elements.
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
719
Fig. 19. An H-shaped FEM with 312 triangular elements.
Fig. 18. Pattern of the resulted flexibility matrix.
Table 2 Comparison of the results. Method
Number of non-zero entries
CPU time (s)
Present method Graph-theoretical method LU factorization method
2982 2982 28,632
23.041 14.547 28.094
For this sparse basis, the total number of non-zero entries is 2982. The pattern of sparse null basis and the resulted flexibility matrix are shown in Figs. 17 and 18, respectively. The CPU time and the number of non-zero entries for the LU factorization method, the graphtheoretical method [17] and the present approach are depicted in Table 2 for comparison. For this sparse basis, the total number of non-zero entries is 2982. The pattern of sparse null basis and the resulted flexibility matrix are shown in Figs. 17 and 18, respectively. Example 4. An H-shaped FEM with 312 triangular elements is considered as shown in Fig. 19. This model has the same properties as that of Example 3. The present method leads to the formation of 430 null vectors with two entries (double edges) and 119 complementary null vectors with 12 entries, corresponding to highly sparse and banded flexibility matrix as shown in Fig. 20. The results obtained for all 4 examples are identical to those obtained by the integer linear programming formulation of Section 8. However, the only difference is the computational time which is much lower for the ACS algorithm. 10. Conclusions In this paper, an ant colony system is developed for the formation of sparse null basis leading to sparse self-stress matrices, and correspondingly highly sparse flexibility matrices for triangular plane stress and plane strain finite element models. An integer linear programming formulation is also presented to evaluate the quality of the solutions obtained by the proposed ACS algorithm. The method
Fig. 20. Pattern of the resulted flexibility matrix.
presented in this paper is not confined to plane stress and plane strain triangular elements and it can easily be extended to other finite elements. Acknowledgement The first author is grateful to Iran National Science Foundation for the support. References [1] J.C.de.C. Henderson, E.A.W. Maunder, A problem in applied topology, J. Inst. Math. Appl. 5 (1969) 254–269. [2] E.A.W. Maunder, Topological and linear analysis of skeletal structures, Ph.D. Thesis, London University, Imperial College, March 1971.
720
A. Kaveh, M. Daei / Finite Elements in Analysis and Design 45 (2009) 710 -- 720
[3] A. Kaveh, Application of topology and matroid theory to the analysis of structures, Ph.D. Thesis, London University, Imperial College, March 1974. [4] A. Kaveh, Improved cycle bases for the flexibility analysis of structures, Comput. Methods Appl. Mech. Eng. 9 (1976) 267–272. [5] P.H. Denke, A general digital computer analysis of statically indeterminate structures, NASA-TD-D-1666, 1962. [6] J. Robinson, Integrated Theory of Finite Element Methods, Wiley, New York, 1973. [7] A. Topçu, A contribution to the systematic analysis of finite element structures using the force method, Doctoral Dissertation, Essen University, 1979 (in German). [8] I. Kaneko, M. Lawo, G. Thierauf, On computational procedures for the force methods, Int. J. Numer. Methods Eng. 18 (1982) 1469–1495. [9] J.R. Gilbert, M.T. Heath, Computing a sparse basis for the null space, SIAM J. Alg. Disc. Methods 8 (1987) 446–459. [10] T.F. Coleman, A. Pothen, The null space problem I; complexity, SIAM J. Alg. Disc. Methods 7 (4) (1986) 527–537.
[11] T.F. Coleman, A. Pothen, The null space problem II; Algorithms, SIAM J. Alg. Disc. Methods 8 (4) (1987) 544–561. [12] S.N. Patnaik, An integrated force method for discrete analysis, Int. J. Numer. Methods Eng. 6 (2) (1973) 237–251. [13] S.N. Patnaik, L. Berke, R.H. Gallagher, Integrated force method versus displacement method for finite element analysis, Comput. Struct. 38 (4) (1991) 377–407. [14] A. Kaveh, A combinatorial optimization problem; optimal generalized cycle bases, Comput. Methods Appl. Mech. Eng. 20 (3) (1979) 9–52. [15] A. Kaveh, Structural Mechanics: Graph and Matrix Methods, third ed, Research Studies Press (Wiley), Somerset, UK, 2004. [16] M. Dorigo, L.M. Gambardella, Ant colony system: a cooperative learning approach to the traveling salesman problem, IEEE Trans. Evol. Comput. 1 (1) (1997) 53–66. [17] A. Kaveh, K. Koohestani, N. Taghizadieh, Efficient finite element analysis by graph-theoretical force method, Finite Elem. Anal. Des. 43 (2007) 543–554.