Tracing periodic orbits in 3D galactic potentials using the particle swarm optimization method
Ch. Skokos, Konstantinos E. Parsopoulos, Panos A. Patsis, Michael N. Vrahatis · 2003
The class of Swarm Intelligence algorithms consists of stochastic optimization methods that exploit a population of interacting individuals to probe the search space simultaneously. The ability to work on nondifferentiable and discontinuous functions using only function value information renders these algorithms a useful tool, especially in cases where classical optimization methods fail. In this contribution, we apply Particle Swarm Optimization for locating periodic orbits in a 3D Ferrers bar model. An appropriate scheme that transforms the problem of finding periodic orbits to the corresponding problem of detecting the global minimizers of a function defined on the Poincare Surface of Section of the Hamiltonian system is employed. We succeeded in tracing systematically several periodic orbits of the system, a large fraction of which is reported for the first time. In particular, we found families of 2D and 3D periodic orbits, associated with inner resonance higher than the 8:1 resonance, and not belonging to the x1-tree. We were also able to locate a plethora of pperiodic orbits with p>1. 1 DESCRIPTION OF THE ALGORITHM 1.1 The Particle Swarm Optimization algorithm Particle Swarm Optimization (PSO) belongs to the category of Swarm Intelligence methods. The ideas that underlie PSO are inspired by the social dynamics of flocking organisms, which are governed by fundamental rules like nearest-neighbor velocity matching and information sharing. Assume the problem of finding a global minimizer of an n-dimensional function, f, which is called the objective function, defined in an n-dimensional search space, S O R. PSO is a population based algorithm, i.e., it exploits a population of individuals to probe promising regions of the search space, simultaneously. In this context, the population is called a swarm and the individuals (i.e., the search points) are called particles. Each particle moves with an adaptable velocity within the search space, and retains a memory of the best position it ever encountered. In minimization problems, the best positions are characterized by lower function values. In the local variant of PSO, each particle is assigned to a neighborhood that consists of a pre-specified number of particles, and the best position ever attained by the particles that comprise the neighborhood is communicated among them. On the other hand, in the global variant of PSO, the whole swarm is considered as the neighborhood of each particle. The local variant is characterized by better exploration behavior, i.e., more thorough search, while the global variant has better exploitation capabilities, i.e., faster convergence towards the best solutions detected so far. In the current work, we employed the local variant, due to its better exploration capabilities in multimodal problems, such as the detection of periodic orbits. Consider a swarm consisting of N particles. Each particle is an n-dimensional vector in S, Charalampos D. Skokos, Konstantinos E. Parsopoulos, Panos A. Patsis, and Michael N. Vrahatis. 1293 T i i1 i2 in X = (x , x , , x ) S, i=1,2, ,N ∈ ... ... , (1) where () denotes the transpose of a matrix. The velocities of the particles are also n-dimensional vectors, T i i1 i2 in V = (u , u , , u ) , i=1,2, ,N ... ... . (2) The best previous position encountered by the i-th particle is a point in S, denoted by T i i1 i2 in P = (p , p , , p ) S ∈ ... . (3) The positions, Xi, as well as the velocities, Vi, are initialized randomly and uniformly within S. The best positions, Pi, are initially set equal to Xi. Each particle is evaluated according to the objective function, f, i.e., the value f(Xi) is computed for all particles. Obviously, at the initialization phase, it holds that f(Pi) = f(Xi). Let Ni = (Xi-r, ..., Xi-1, Xi, Xi+1, ..., Xi+r), be a neighborhood of radius r of the i-th particle, Xi. Then gi is defined as the index of the best particle in the neighborhood of Xi, i.e., i g j f(P ) f(P ), j=i-r, , i+r ≤ ... . (4) The neighborhood's topology is usually cyclic, i.e., the first particle, X1, is assumed to follow after the last particle, XN. For example, a neighborhood of radius 3 of the particle X2 consists of the particles XN-1, XN, X1, X2, X3, X4 and X5. The particles are moving in S according to the equations: ( ) (q 1) (q) (q) (q) (q) (q) i i 1 1 i i 2 2 gi i V χ V c r (P X ) c r (P X ) + = + + , (5) (q 1) (q) (q 1) i i i X X V + + = + , (6) where i = 1, 2, ..., N; χ is a parameter called constriction factor; c1 and c2 are two fixed, positive parameters called cognitive and social parameter, respectively; r1, r2 are random numbers uniformly distributed in [0,1]; and q stands for the counter of iterations. The objective function, f, is computed again at the new positions Xi of the particles, and the best positions, Pi, i=1,2,...,N, are updated along with the new indices gi. Then, Eqs. (5) and (6) are applied again to proceed to the next iteration of the algorithm. The procedure stops when a sufficiently good solution is detected or a maximum number of iterations is reached. Let us now discuss the role of the various parameters that appear in equations (5) and (6). The constriction factor, χ, is used as a mechanism to control and adjust the magnitude of the velocities. Cognitive parameter c1 controls the effect of the knowledge of the best previous position, Pi, of i-th particle on its velocity, while social parameter c2 plays a similar role but it concerns the best previous position, Pgi, attained by any particle in the neighborhood. The size, N, of the swarm as well as the neighborhood radius, r, are usually set arbitrarily. However, it is a common belief in evolutionary computation that a population size equal from 2 to 10 times the dimension of the problem at hand is a good initial guess. Moreover, the neighborhood's size shall be problemdependent. In easy problems, larger neighborhoods result in faster convergence without loss of the algorithm's efficiency, while, in harder problems with a plethora of local minima, smaller neighborhoods are considered a better starting choice. 1.2 Detecting further minimizers thought Deflection PSO is able to detect one, in general arbitrary, minimizer of the objective function, per run. However, in some applications, several minimizers of the objective function are required. Restarting the algorithm does not guarantee the detection of a different minimizer. In such cases, the deflection technique can be used. This technique consists of a transformation of the objective function, f, once a minimizer Xi, i=1,...,nmin, has been detected, * -1 i i F(X) = [tanh(λ ||X-X ||)] f(X) , (7) Charalampos D. Skokos, Konstantinos E. Parsopoulos, Panos A. Patsis, and Michael N. Vrahatis. 1294 Figure 1. The effect of the deflection procedure on the function f(x) = cosx + 0.1, at the point x = π/2, for λ=1 (a), and λ=0.5 (b). where λi, i = 1,...,nmin, are nonnegative relaxation parameters, and nmin is the number of the detected minimizers. The transformed function, F, has exactly the same minimizers with f, with the exception of Xi. Alternative configurations of the parameter λ result in different shapes of the transformed function. For larger values of λ, the impact of the deflection technique on the objective function is relatively mild. On the other hand, using 0 0 is a constant, instead of f. The function f possesses all the information regarding the minimizers of f, but its global minimum is increased from zero to c. The value of c does not affect the performance of the algorithm and, thus, if there is no information regarding the global minimum of f, it can be selected arbitrarily large. The effect of the deflection procedure on the function f(x) = cosx+0.1, at the point x=π/2, is illustrated in Fig. 1. 2 THE GALACTIC POTENTIAL The 3D galactic bar model that is used in our study is described in detail by Skokos et al.. It consists of a Miyamoto disk, a Plummer bulge and a Ferrers bar. The potential of the Miyamoto disk is given by the formula, ( ) D D 2 2 2 2 2 G M V x y A B z = + + + + , (8) where MD represents the total mass of the disk, A and B are the horizontal and vertical scale lengths, respectively, and G is the gravitational constant. The bulge is a Plummer sphere, i.e., its potential is given by S S 2 2 2 2 s G M V x y z e = + + + , (9) where es is the bulge scale length and MS is its total mass. Finally, the bar is a triaxial Ferrers bar with density defined by Charalampos D. Skokos, Konstantinos E. Parsopoulos, Panos A. Patsis, and Michael N. Vrahatis. 1295 2 2 B 105 M (1-m ) for m 1 32 π abc ρ(m) 0 for m 1 ⎧ ≤ ⎪ ⎪ = ⎨ ⎪ ≥ ⎪ ⎩ (10)