A NEW ALGORITHM IN NONLINEAR ANALYSIS OF STRUCTURES USING PARTICLE SWARM OPTIMIZATION

Solving systems of nonlinear equations is a difficult problem in numerical computation. Probably the best known and most widely used algorithm to solve a system of nonlinear equations is the Newton-Raphson method. A significant shortcoming of this method becomes apparent when attempting to solve problems with limit points. Once a fixed load is defined in the first step, there is no way to modify the load vector should a limit point occur within the increment. To overcome this defect, displacement control methods for passing limit points can be used. In the displacement control method, the load ratio in the first step of an increment is defined so that a particular key displacement component will change only by a prescribed amount. In this paper, the load ratio is obtained using a Particle Swarm Optimization (PSO) algorithm so that the complex behavior of structures can be followed, automatically. The design variable is load ratio, and its unbalanced force is also considered an objective function in the optimization process. Numerical results are performed under geometrical nonlinear analysis, elastic post-buckling analysis and inelastic post-buckling analysis. The efficiency and accuracy of proposed method are demonstrated by solving numerical examples. ABSTRAK: Menyelesaikan sistem persamaan tak linear adalah masalah yang sukar dalam pengiraan numerik. Mungkin algoritma yang paling terkenal dan paling banyak digunakan untuk menyelesaikan sistem persamaan tak linear adalah kaedah NewtonRaphson. Satu kekurangan besar kaedah ini menjadi jelas apabila mencuba menyelesaikan masalah dengan titik had. Setelah beban tetap ditakrifkan dalam langkah pertama, tidak ada cara untuk mengubahsuai vektor beban jika titik had berlaku dalam pertambahan. Untuk mengatasi kecacatan ini, kaedah kawalan sesaran untuk melepasi titik had boleh digunakan. Dalam kaedah kawalan sesaran, nisbah beban dalam langkah pertama peningkatan akan ditakrifkan supaya komponen sesaran utama tertentu akan berubah hanya dengan jumlah yang ditetapkan. Dalam kertas ini, nisbah beban diperolehi dengan menggunakan algoritma pengoptimuman zarah swarm (PSO) supaya tingkah laku yang kompleks daripada struktur boleh diikuti, secara automatik. Pembolehubah reka bentuk adalah nisbah beban dan daya tidak seimbang yang juga dianggap sebagai fungsi objektif dalam proses pengoptimuman. Keputusan numerik dilakukan di bawah analisis geometri tak linear, analisis pasca-lengkokan anjal dan analisis pasca-lengkokan tidak anjal. Kecekapan dan ketepatan kaedah yang dicadangkan akan dibuktikan dengan menyelesaikan contoh numerik.


INTRODUCTION
The structural problems with geometrically nonlinear features can often be analyzed by solving a system of nonlinear equations (algebraic or differential) to determine the path of nonlinear load-displacement.The literature on analysis of these problems is quite extensive and provides several approaches, most notably the incremental stiffness procedure [1], the perturbation method, the Newton-Raphson method and its modified variations, [1,2], the initial value approach [3], and the self-correcting incremental procedure [1,2], through all of which nonlinear equilibrium equations could be solved in an efficient manner.
Literature also provides several analytical approaches specifically developed for assessing the nonlinear behavior exhibited by structural members such as trusses; these include the approach suggested in a study by Wempner [4], where complete sets of equilibrium paths of a number of Mises trusses have been provided, the method of Saffari et al. [5], who have used the Newton-Raphson iterative algorithm in conjunction with the flow path normal in nonlinear static analysis, and the method of Saffari et al. [5], who have analyzed the nonlinear behavior of trusses via a modified normal flow algorithm.Other notable works on this subject include -but are not limited to -the study by Papadrakakis and Gantes [6], where the optimization of geometrically nonlinear shallow trusses with stability-related constrains has been examined, the work of Greco et al. [7], where space trusses has been analyzed by a new geometric nonlinear formulation based on nodal positions (instead of nodal displacements), and the article of Papadrakakis [6], where trusses have undergone a second-order and large deflection analysis using the dynamic relaxation scheme.There has also been much progress as the result of the study conducted by Huang and Vahidi [8], where the snap-through buckling of two sample trusses has been analyzed by the elastic theory of prismatic bars, the article of Kassimali and Parsi-Feraidoonian [9], which has been focused on nonlinear behavior of pre-stressed cable trusses, the work of Bellini [10], where the snap-buckling problem has been solved via a novel mathematical model, the article of Ramesh and Krishnamoorthy [11], where the dynamic relaxation method has been used to analyze the inelastic post-buckling of trusses, and the work of Thai and Kim [12], where space truss structures with geometric and material nonlinearities have been analyzed for their inelastic large deflections.Although solving a system of nonlinear equations could prove to be very difficult, literature has provided a number of methods for this purpose, nevertheless they do little in improving the computation time required by this procedure.The Newton-Raphson method, which attacks the problem via an iterative procedure of solving linear algebraic equations (corresponding to their non-linear counterparts), is perhaps one of the best techniques for solving nonlinear equations since it provides the quadratic convergence characteristics well suited for this application.
There are a number of iterative methods specifically developed to improve the performance of the Newton-Raphson technique [13][14][15], but they often have questionable practicality mainly because they still follow an approach that requires the computation of second or higher derivatives, which can be very time-consuming and computationintensive.There are some techniques that can omit the second derivatives (Hessian matrix) from this process [16,17], and these are classified into two groups: one-step techniques, and two-step techniques [14].Two-step techniques, like the one used in this paper, are in fact the combinations of the Newton-Raphson method and another one-step method.This approach, constructed in conformity with predictor-corrector methods, and so taking most frequently Newton's method in the first iteration, is used in this paper.The common disadvantages of one-point methods, most notably their high computational intensity and their inadequate convergence, can be sidestepped using a series of very effective techniques known as multipoint iterative solvers [18].Rezaiee-Pajand et al. [19,20] the efficiencies and capabilities of residual load minimization, normal plane, updated normal plane, cylindrical arc length, work control, residual displacement minimization, generalized displacement control, and modified normal flow were evaluated.Rezaiee-Pajand and Naserian [21] employed incremental-iterative methods for the base of iteration steps by setting each area to zero, and minimizing its perimeter separately.Application of cubic spline on large deformation analysis of structures has been proposed by Saffari et al. [22].A Dynamic Relaxation (DR) method, suitable for nonlinear structural analysis, has been suggested by Rezaiee-Pajand and Estitri [23].
The weakness of the Newton-Raphson method is that, in the course of computing structural response, having a solution point sufficiently close to the limit point will cause the method to exhibit significant divergence, possibly leading to significant errors in computation of failure loads.To deal with this problem, one can use a number of methods provided by literature, most important of which are classified as displacement control methods [24].In the present study, to compensate for this inability of the Newton-Raphson method, a normal flow algorithm is used.
The objective of this study is to assess the optimization of the structures that must undertake large deformations.A Eulerian formulation is used to model the geometrically nonlinear structural behaviors and then a normal flow algorithm is used, which is a displacement-control iterative method to determine the point of equilibrium.Finally, Particle Swarm Optimization (PSO) is used with the objective of effectively reducing the computational intensity of the analysis.

GEOMETRICAL NONLINEAR ANALYSIS OF SPACE TRUSSES
The system equilibrium equations are generally defined as: where {f} is the resultant of the nodal internal loads and {P} is the external nodal force or load.According to equation expressing the member force deformation, {f} is a function of {δ} and is predominantly non-linear.Equation 1, in its differential form, is expressed as: In the above equation, {∆δ} is the increment of displacement, {ΔP} is the increment of load, and [τ] denotes the tangent stiffness matrix.

Tangent Stiffness Matrix of the Structural Member
The end displacements have an incremental relationship with end member forces, which is expressed by the following equation [5]: In the above equation, [T] is the tangent stiffness matrix of the member.The following formula provides a method to calculate this parameter: The structure of geometric matrix [g] is as follows [5]:

Newton-Raphson Technique
Of several methods available to deal with non-linear problems, the Newton-Raphson technique is certainly one of the best and is generally the basis of other iterative algorithms developed to deal with this problem.While this section presents a brief introduction to this technique, the more detailed information regarding this subject can be found in [2].This technique starts its first iteration with a procedure similar to the one followed in the linear incremental method, with the only difference being the calculation of member forces and their transformation into the global coordinates at the end of procedure.The next step of this technique is the calculation of imbalance between the external loads and the internal nodal forces.This procedure of iteration continues until it reaches a predefined termination condition, often based on a convergence criterion.

Modified Newton-Raphson Technique
The difference between the normal Newton-Raphson technique and the modified Newton-Raphson technique is the frequency of reformation of the tangent stiffness, which will be repeated for every iteration of the normal technique, but will be done only once in the first iteration of the modified one.Omitting the reformation process increases the number of iterations required to reach convergence, but reduces the computational intensity of each single iteration, leading to a shorter run time for the entire process.

NORMAL FLOW ALGORITHM
This paper uses the normal flow algorithm, which is an effective approach for tracing the complete equilibrium path for this very same purpose.Now, let i be the number of the steps, j be the number of the modifying iterations, and   j i P be the total load on the structure.In that case: where  is the product of a total ratio and {Pref} is a reference external load exerted via a series of load increments.Figure 1 shows the schematic representation of this method.More detailed information regarding this method can be found in [5].The following formulation is the general form of a nonlinear system of equilibrium equations: In this formulation, ) , ( represents the Jacobian matrix of order N×(N+1), and j i S denotes the step size of Newton-Raphson technique.This step size can be calculated via Eq.( 8): The step size of Newton-Raphson technique for a normal flow algorithm can be determined by finding the minimum of an infinite number of solutions obtainable for Eq. ( 9).The procedure stated below must be followed to find this minimum solution: First, the following equation should be used to calculate an initial solution {V}: where: In the above equation,   1   j i Q denotes the vector of the unbalanced forces,   j i I  represents the vector of tangential displacement at the point of convergence; and   j i R   denotes the vector of unbalanced displacements such that, where   j i F denotes the vector representing the resultant of internal forces at the nodes.
Next, the following equation should be used to obtain the vector of unbalanced force: And finally, the following equation should be used to obtain the minimum solution of the norm: This means that the process of calculating the unknown vectors {V} and   j i R   by solving the system of equations should be repeated at each iteration.

PARTICLE SWARM OPTIMIZATION (PSO)
This study uses the PSO method to optimize its objective.Authors of [25] have reported that the use of Binary PSO (BPSO) for structural optimization has resulted in significant improvement in solutions.First proposed by Eberhart et al. [26] in the 1990s, PSO has been inspired by the social and group behavior of animals such as fish, insects, and birds.The mentioned authors were seeking to perform a socio-cognitive study on the group behavior of bird swarms and to achieve that purpose, they developed this amazingly effective optimization algorithm.This algorithm starts with a number of particles (or birds), which are spread randomly in the search space of the objective function (swarm of birds in a habitat).Here, particles represent the potential solutions to the problem.At each iteration, all particles move around the search space based on the best positions (solutions) obtained at the previous iteration.This algorithm obtains this best solution by calculating the objective function and the fitness of all particles in the search space [26].Algorithm moves the position of particles based on the following equations: Where, k is the number of iterations, X is the current position of i-th particle and is the velocity of this particle; Pi is the best position previously attained by the i-th particle (called pbest) and Pg is the best position attained globally by all particles (called gbest); r1 and r2 are two uniform random sequences generated from interval [0, 1]; w k denotes the inertia weight that will be used to modify (reduce) the previous velocity of the preserved particle.According to Shi and Eberhart [25], the cognitive and social scaling parameters c1 and c2, should be selected such that c1 = c2 = 2.0, which consequently allows the product c1r1 or c2r2 to have a mean of 1.The minimum and maximum possible velocities are also incorporated into this design and are denoted by max i V and min i V (Fig. 2).This paper uses a modified Euclidian criterion in the following form to express the displacement control: In the above relationship, ε is the defined error of calculation.
In the updating scheme used in this paper, the number of performed iterations affects the load increment parameter.Meanwhile, this algorithm uses the following formula to determine the sign of the determinant of the tangential stiffness matrix resulting from the previous step: In the above formula, JD is the predefined number of iterations, JM is current count of performed iterations minus one, and γ is a constant.

NUMERICAL RESULTS
In this section, the results of a computer program developed to implement the algorithm described in previous section are presented.The duration of calculation, or total runtime, can be obtained by entering the corresponding command into the program.This program is run on a 64-bit PC with a Core™ i7-6700T Processor, 8M Cache, and uses similar parameters for all problem instances.These problem instances act as a measure to gauge the performance of the proposed method.Here, the objective is to perform elastic analysis, elastic post-buckling (EPB) analysis, and inelastic post-buckling (IPB) analysis on typical truss structures using the proposed method, the conventional Newton-Raphson technique and other sub-stepping algorithms.In all tests, r = 0.4, X1 = 50 and X2 = 100, and the maximum error is defined as ε = 10-5 [12].It should also be noted that all algorithms are tested with the smallest possible steps to ensure a uniform and accurate solution.

Problem Instance 1
The following geodesic dome truss has 156 members with identical cross-sections A = 6.5 cm 2 , and I = 1 cm 4 , and 61 nodes with outer ones attached by pin supports (Ramesh and Krishnamoorthy [11]).
This truss is built with a material with an elasticity modulus of E = 6895 kN/cm 2 and yield strength of Fy = 400 kN/cm 2 .Figure 4 shows the load-displacement curves obtained as a result of elastic analysis, post-buckling (EPB) analysis, and inelastic post-buckling (IPB) analysis conducted on this truss structure.

Elastic Analysis EPB Analysis IPB Analysis
The run times achieved in each mode of analysis by the proposed method and the conventional Newton-Raphson technique are shown in Table 1.This result demonstrates the superiority of the proposed method in all modes of analysis.

Problem Instance 2
The semi-spherical truss shown in Fig. 5 is composed of 168 elements with identical cross-sections with A = 50.431cm 2 , and I = 52.94cm 4 , and 73 nodes with total freedom degree of 147 [12].In this problem, a vertical load P = 500 kN is applied on the apex of this truss.The pin supports defined at the outer nodes of the truss restrict the off-plane displacements.
Fig. 5: Semi-spherical truss from problem instance 2 (all dimensions in cm) [12].This truss is built with a material with elasticity modulus of E = 2.04×104 kN/cm 2 and yield strength of Fy = 25 kN/cm 2 .In this problem, the parameters of Eq. ( 17) are 01 .0  2.These runtimes also show that the proposed method outperforms other algorithms and achieves the quickest convergence.Fig. 6: Load-displacement curves for node 2 of the semi-spherical truss.

Problem Instance 3
In this problem, a Schewdeler's truss is assessed with properties stated as EA = 640×103 kN, Fy = 25 kN/cm 2 , I = 30.04cm 4 built in the form of 264 elements and 97 nodes.Here, the outer nodes are again attached with pin supports [7].

Fig. 4 :
Fig. 4: Load-displacement curves for the apex of a geodesic dome truss.
displacement curves obtained via elastic analysis, post-buckling (EPB) analysis, and inelastic post-buckling (IPB) analysis conducted on this truss are shown in Fig.6.The run times achieved in each mode of analysis by the proposed method and other tested technique are shown in Table

Table 1 :
The runtime achieved for problem instance 1 (sec)

Table 2 :
The runtime achieved for problem instance 2 (sec)