S. Kok and J. A. Snyman

Abstract

A novel dynamic interacting particle swarm optimization algorithm (DYN-PSO) is proposed. The algorithm can be considered to be the synthesis of two established trajectory methods for unconstrained minimization. In the new method, the minimization of a function is achieved through the dynamic motion of a strongly interacting particle swarm, where each particle in the swarm is simultaneously attracted by all other particles located at positions of lower function value. The force of attraction experienced by a particle at higher function value due to a particle at a lower function value is equal to the difference between the respective function-values divided by their stochastically perturbed position difference. The resultant motion of the particles under the influence of the attracting forces is computed by solving the associated equations of motion numerically. An energy dissipation strategy is applied to each particle. The specific chosen force law and the dissipation strategy result in the rapid collapse (convergence) of the swarm to a stationary point. Numerical results show that, in comparison to the standard particle swarm algorithm, the proposed DYN-PSO algorithm is promising.

1. Introduction

A new direct search method using only function values is proposed for finding a local minimizer with associated function value of a real valued function , . The method proposed may be considered as the synthesis of two unconventional trajectory methods for the unconstrained minimization of a multivariable function. The first is the dynamic method of Snyman [1, 2], and the second method is the particle swarm optimization (PSO) method of Eberhardt and Kennedy [3]. In the dynamic method, also known as the leap-frog algorithm, the minimum of the function is sought by considering the dynamic motion of a single particle of unit mass in an -dimensional force field, where the potential energy of the particle is represented by the function to be minimized. In the computation of the numerical trajectory (by means of the leap-frog integration scheme of Greenspan [4]), an interfering strategy is applied to the motion of the particle by extracting kinetic energy whenever it moves “uphill" along its trajectory. In this way, the particle is forced to converge to a local minimum. This method requires the availability of the gradient vector of the function, denoted , the negative of which represents the force acting on the particle. The particle's acceleration is therefore proportional to . This is in contrast to classical gradient-based optimization, where position updates (or parts thereof) are proportional to . In the PSO method, the motion of a swarm of loosely interacting particles is considered. In this method, each particle is attracted to the best location (lowest function value position) along its path, as well as to the globally overall best position over all the particle trajectories to date. This method requires no gradient information and may therefore be considered a direct search method.

In the new method proposed here, the minimization of a function is achieved through the dynamic motion of a strongly interacting particle swarm, where each particle in the swarm is simultaneously attracted by all other particles located at positions of lower function value. The specific force law for the interaction between the individual particles within the swarm dictates that the force of attraction experienced by a particle at higher function value position due to a particle at a lower function value position be equal to the difference between the respective function values divided by their distance of separation. The resultant motion of the particles under the influence of the attracting forces is computed by solving the associated equations of motion numerically by using again the leap-frog numerical integration scheme. An energy dissipation strategy, similar to that used in the original dynamic method, is applied to each particle by extracting kinetic energy from a particle whenever it moves “uphill". The specific chosen force law and the dissipation strategy result in the rapid collapse (convergence) of the swarm to a stationary point. To prevent the collapse of the swarm to a minimum in a subspace of , the computed components of the respective attracting forces are individually stochastically perturbed in computing the respective particle trajectories. Because of the strong interaction, the resulting algorithm converges rapidly to a local minimum, as does the original dynamic algorithm, but now without the need of explicit gradient information. Note that the method proposed here is not similar to the method proposed by Engelbrecht [5], where the position of only the global best particle is adjusted using the dynamic method, with the remainder of the swarm using a standard PSO algorithm.

The proposed dynamic particle swarm optimization (DYN-PSO) algorithm promises to be an extremely reliable, robust, and easy-to-use method. As a preliminary test, it was decided to evaluate its performance against that of both trajectory methods it was constructed from, namely, a standard PSO algorithm [6] and the dynamic method of Snyman. The choice of a standard PSO algorithm is a deliberate one. The intention with this paper is to determine whether or not DYN-PSO deserves further investigation, not to suggest that this is a superior algorithm. A comparison to a standard PSO rather than a highly refined PSO variant serves this purpose better.

Before proceeding, the constrained molecular dynamics PSO by Poli and Stephens [7] deserves to be mentioned. In essence, their method is similar is spirit, but the details differ. They embed the cost function (the height above the search space) as an artificial coordinate, and then constrain the particles “to be on the fitness landscape” via an equality constraint. Their general formulation allows for different types of forces, one of which is a gravity-like force. In our implementation, this gravity-like force is the only external force. Poli and Stephens also consider particle-interaction forces. In particular, they discuss the case of particles connected by springs and energy dissipation that occur via friction/viscous damping. They then also numerically integrate the resulting differential equation, but elect to use the forward Euler method. Their method still requires cost function gradients, as well as second-order derivatives, which they approximate numerically. Limited numerical testing was performed.

2. Description of the Dynamic-PSO Method

2.1. Computation of Particle Trajectories

The DYN-PSO method is started by generating at time , a swarm of particles, each of unit mass, with random initial positions denoted by , within the region (“box”) of interest in . Initially, at , these particles all have zero velocities, that is, . We now postulate that at time each particle experiences a force , which is to be a function of the positions and corresponding function values of all the particles at time . The explicit analytical form of the force law giving will be discussed in the next subsection. Thus, the trajectories of the particles are given by the solution to the system of initial value problems(1) with initial conditions (2) for

In practice, these equations are solved numerically by discretizing the time interval into time steps , and computing for approximations to and to at discrete time mesh points by some suitable numerical integration scheme. Here, we use the simple leap-frog numerical integration scheme of Greenspan [4].

Given and for , then for iterations compute for : (3)where denotes the resultant force on particle due to the individual forces of all the other particles at respective positions , , (see next subsection).

This scheme has been found to be stable for sufficiently small time steps . It is also approximately energy conserving in the absence of energy dissipating forces, and was successfully used in the original single-particle dynamic method of Snyman [1].

2.2. The Interacting Force Law

In the original dynamic method of Snyman [1, 2], a single particle is considered. The force acting on this particle is equal to the negative of the function gradient. Since we propose a direct search method, we no longer have gradient information available. Rather, we will use information available in the swarm to generate the particle forces.

We now postulate that at iteration each particle in the swarm is simultaneously attracted by all other particles located at current computed positions of lower function value. No force is exerted on particle by particles at higher function value positions. The explicit force law that we assume here dictates that the force of attraction experienced by a particle at higher function value due to a particle at a lower function value position is equal to the difference between the respective function values, divided by their distance of separation. If the particle positions are randomly distributed, this force law can be viewed as a coarse finite difference computation, with the finite difference step equal to the distance of separation. As the particle positions become biased, the analogy between this proposed force law and a finite difference computation breaks down. This breakdown occurs since multiple forces in the same direction add up, opposed to standard finite difference perturbations which form an orthogonal set. Even considering generalized finite differences [8], the gradient components due to multiple perturbations in the same direction are averaged, rather than added together as in our proposed force law. A stochastic element is introduced into the proposed force law by randomly perturbing the direction of action of the attracting force. More explicitly and precisely, at each iteration the resultant force exerted on particle at current position with function value by the other particles at current positions with corresponding function values is given by the force with components(4)for . Here and is, for each , an independent random number in [0, 1]. Note that the multiplication of a uniform random number in [0,1] by 2 results in a scheme that on average assigns unit weight to each of the computed components. If this stochastic element is omitted, the proposed algorithm suffers from premature convergence due to collapse of the swarm to a subspace in . Similar observations have been made for the linear PSO [9].

This proposed force law is in contrast to the standard PSO, where particles are attracted to historic positions (personal and global bests) rather than the current position of other particles. In addition, this force law automatically results in a dynamic neighborhood, where it is possible that completely different particles exert forces on particle in iteration , as compared to iteration . In a given iteration , the complete range of interaction is also covered. One extreme is the particle with the current worst function value, which experiences forces from all other particles. The other extreme is the particle with the current best function value, which experiences no force and travels at constant velocity.

2.3. Energy Dissipation Strategy

In computing the trajectories of the particles , for the function values at are monitored at each iteration so that the best (lowest) function value and the corresponding best position along each trajectory are recorded. The current overall globally best function value and the corresponding position are also recorded.

The following energy dissipation strategy is now applied to ensure local descent of a particle and the overall collapse of the swarm to a local minimum. Whenever a particle moves “uphill” at iteration , that is, when , then set and set . Notice that the current velocity is recomputed to be half the average velocity over the past two time steps. This interference will normally result in a decrease in kinetic energy and together with the “backward” adjustment of the position will initiate controlled motion of the particle towards a position of lower function value. The trajectory of particle is assumed to have converged if the relative function value difference from iteration to , that is, , is less than some prescribed tolerance . The computation of the trajectories are continued until a sufficient number () of the particles have converged to a stationary point. In practice, we choose .

A formal presentation of the basic DYN-PSO algorithm is presented in Algorithm 1. For the sake of clarity and simplicity of presentation of the algorithm, the function value computation and monitoring procedure, and the recording of best local and global function values and corresponding best local and global positions, and are not explicitly listed, but are implicitly assumed to be done in the execution of Algorithm 1. Also the computation of the forces , according to (4), is assumed to have been done as the need for the arises in Algorithm 1.

Algorithm 1: Basic DYN-PSO algorithm.

2.4. Selection of Suitable Integration Time Step

An outstanding matter is the selection of an appropriate time step to be used in the leap-frog integration scheme. This value can be chosen arbitrarily, but if chosen too large may result in individually erratic and unstable trajectories that fail to converge because of very large zigzagging steps being taken in space. On the other hand, if is too small, the steps in space will be correspondingly small and the collapse of the swarm may be very slow, requiring an excessively large number of iterations for convergence.

Many different schemes may be proposed to automatically select and control the time step so that acceptable convergence rates are obtained. Here, we firstly select an initial time step which guarantees sufficiently large trajectory steps in space for all the particles. Secondly, after the trajectories are initiated, for each iteration and for each particle , the magnitude of the actual step taken in space is monitored; and if larger than some specified step limit, the time step associated with the particular trajectory is reduced, that is, we now allow for different time steps , for the different particles. The details of these additional automatic time step control procedures follow below.

In Step 1 of Algorithm 1, after the generation of the random particles, compute the average magnitude of the forces acting on the particles . An associated average computed initial step size follows from the leap-frog scheme as . Requiring initially, on average, a step size of where is the diameter of the initial variable “box", we initially select as sufficiently large initial time steps , which is now inserted in Step 1 just after the generation of the initial random particle positions. However, for a particular particle, this initial choice for the time step may still be too large and result in erratic behavior along the computed trajectory. Thus, as a further control measure, the magnitude of the step taken by each particle is monitored at each iteration, and if it exceeds a specified limit , the current velocity is scaled down by a factor , the particle's time step for subsequent steps is reduced by a factor , and the step to is recomputed using the rescaled velocity and new time step. More specifically, we introduce at the end of Step 2 of Algorithm 1, for each iteration and for each particle , the additional procedure given in Algorithm 2.

Algorithm 2: Adjustment of time step .

In practice, good choices for the additional parameters introduced here are and . These values are used in the numerical tests that are reported here.

One complication of the above strategy arises if the nature of the cost function changes dramatically as the search proceeds. The initial time steps are appropriate during the initial stages of the search. However, the average force acting on the particles often decreases over time. Hence, the average particle step size also decreases over time. This seems appropriate for swarm convergence, but the step size decrease may be excessive if the average magnitudes of the forces computed by (4) decrease dramatically. This typically happens in the neighborhood of a local minimum of highly nonlinear functions. To overcome this drawback, we simply recompute the time steps , whenever appropriate, using the recomputed average force . For all the results presented in this paper, the time steps are recomputed every 100 iterations.

3. Numerical Performance of the DYN-PSO Method

3.1. Illustrative Two-Dimensional Trajectories

To illustrate the mechanics of the DYN-PSO algorithm, we compute the trajectories for with 3 particles, that is, the case with . As starting points, we select the vertices . Figure 1 depicts the three computed trajectories (using diameter ) up to the iteration at which point , with .

Figure 1: Example of dynamic PSO particle trajectories for the 2D function .

3.2. Choice of Number of Particles

Throughout the numerical experiments done here, the number of particle trajectories required to converge before termination of the algorithm is taken as . It is now required to obtain an indication of what is a good or an optimum choice for the number of particles to be used for “normal" problems where one expects a single unique global minimum in the region of interest. To get an indication of what it should be, some experiments were performed (with function value tolerance ) on the extended homogeneous quadratic test function (see list of test problems) for and , using different values for and determining the average number of functions evaluations (over 100 independent runs) required for successful convergence in each case. The variation of the number of function evaluations against is shown in Figure 2. Note that for , the number of function evaluations appears to be almost insensitive to over the range to , with a variation of less than 4 percent being recorded. The corresponding variation of average function value at convergence with is depicted in Figure 3. The results appear to indicate that in general a choice of (i.e., with the positions of the particles defining the vertices of an -dimensional simplex in ) is probably a good one. Consequently, the choice is used throughout the experiments performed in the next subsection. Since this guideline is based on a very narrow test, some future effort should be directed towards improved guidelines to determine the number of particles. It is anticipated that highly multimodal problems might benefit from an increased number of particles.

Figure 2: Average number of function evaluations at convergence versus number of particles.

Figure 3: Average function value at convergence versus number of particles.

3.3. Performance on a Set of Test Functions

The newly proposed DYN-PSO algorithm is tested on the following eight test problems.

(i) Homogeneous quadratic (unimodal)(5)subject to for ; , at

(ii) Oren's power function (unimodal) [1](6)subject to for ; , at .

(iii) Extended Rosenbrock (multimodal) [10](7)subject to for ; , at .

(iv) Neumaier 3 (multimodal) [10](8)subject to for ; , at

(v) Extended Manevich (unimodal) [11](9)subject to for , at .

(vi) Zakharov (unimodal) [12](10)subject to for ; , at .

(vii) Griewank (multimodal) [10](11)subject to for ; , at .

(viii) Rastrigin (multimodal) [10](12)subject to for ; , at .

Each test problem is solved 100 times for problem dimension and 30. The DYN-PSO algorithm is compared to the dynamic method of Snyman [1, 2] (STD-DYN) and a standard PSO algorithm (STD-PSO) (http://www.particleswarm.info/Standard_PSO_2006.c), making use of settings proposed by Eberhart and Shi [6] (, constant inertia weight of 0.729 and swarm size of 20). The gradients required by the dynamic method are computed using a forward difference method, using a perturbation size of . These additional function evaluations are taken into consideration in the results of STD-DYN presented below. The convergence histories of the STD-DYN, DYN-PSO, and STD-PSO algorithms are depicted in Figures 4 to 11, for the eight test functions. In these graphs, at each iteration the absolute relative error(13)averaged over the 100 runs, is plotted against the corresponding average number of function evaluations.

Figure 4: Convergence history on the homogeneous quadratic test function.

Figure 5: Convergence history on Oren's power test function.

Two common scenarios exist to terminate optimization algorithms. First, some convergence criteria are satisfied, which indicates that no substantial improvement is likely and search can terminate. Alternatively, a maximum number of function evaluations may be specified. For the test problems considered here, both these scenarios occur. The many local minima present in the multimodal Griewank and Rastrigin test functions prevent the majority of particle trajectories to converge, hence these problems are terminated based on a maximum allowable number of function evaluations. The proposed convergence criteria work well for all the remaining test functions, where we used a convergence tolerance of .

The STD-DYN and DYN-PSO algorithms always locate the global minimum for the homogeneous quadratic, Oren's power function, Neumaier 3, Manevich and Zakharov problems. The STD PSO fails to do so for the 30D Neumaier problem. Also note that out of these 5 test problems, the STD-DYN algorithm is most efficient on the homogeneous quadratic, Oren's power function, and the Zakharov function, all of which are unimodal. The DYN-PSO algorithm performs best on the Neumaier 3 problem and the 30D Manevich function. STD PSO is only more efficient than the DYN-PSO algorithm for Oren's power function.

In the case of the Rosenbrock test function, the DYN-PSO algorithm does not always locate the global minimum. In those cases in which the global minimum is not found, the algorithm converges to the only other local minimum, as reported by Shang and Qiu [13]. In the Rosenbrock experiments, 89 and 96 of the 100 runs converged to the global minimum, for and 30, respectively. The STD-DYN method locates the global minimum 83 and 88 times, for and 30, respectively. Note however that the number of function evaluations is more than a factor 10 less compared to the DYN-PSO algorithm. In the case of the STD-PSO algorithm, the majority of runs do not converge to either the global or local minimum. The number of runs that has a global best less than 1 after function evaluations is 98 and 4, respectively, for n=10 and 30. The absolute relative error in Figure 6 is computed by only averaging over those runs that converge to the global minimum (DYN PSO and STD DYN) or those that have a global best less than unity after function evaluations (STD PSO).

Figure 6: Convergence history on the Rosenbrock test function.

Figure 7: Convergence history on the Neumaier 3 test function.

Figure 8: Convergence history on the Manevich test function.

Figure 9: Convergence history on the Zakharov test function.

Figure 10: Convergence history on the Griewank test function.

Figure 11: Convergence history on the Rastrigin test function.

Finally, the DYN-PSO algorithm never locates the global minimum for the multimodal Griewank and Rastrigin test functions. The dynamic method, which is the backbone of the DYN-PSO algorithm, was developed as a local minimizer. However, the energy dissipation strategy we opted in the DYN-PSO algorithm is unchanged from that of the original dynamic method. This strategy is quite severe on “uphill” moves and hence energy is lost quickly if many local minimizers are present, such as in the Griewank and Rastrigin test functions. Nevertheless, compared to the performance of the STD-DYN method, the DYN PSO is vastly superior on the Rastrigin test function. The STD-DYN method simply locates the first strong local minimum and cannot escape it. On the Griewank function, DYN PSO is superior to STD DYN for 10D, and vice versa for 30D. The STD PSO however far outperforms both STD DYN and DYN PSO for the Griewank function. However, in our experience, the performance of the DYN-PSO algorithm is still comparable to that of more traditional global optimization algorithms. The current algorithm does demonstrate the ability to escape some local minima, similar to the standard PSO algorithm. The performance on the 30D Rastrigin function is especially noteworthy, with a mean function value less than 20 after only 5 000 function evaluations. This is obtained with the standard settings proposed here, with no tuning of parameters to suit the problem.

In summary, the DYN-PSO algorithm seems to inherit the desired properties from both its ancestors. It can efficiently solve unimodal functions, sometimes even more efficiently than the gradient-based local minimizer it is based on, for example, the 30D Manevich problem. This is achieved without making use of gradient information. Therefore, the DYN-PSO method might even work for discontinuous problems. This efficient local character is blended with nonlocal behavior, where the swarm provides sufficient information to solve multimodal problems, illustrated best on the Rastrigin function.

4. Conclusions

We have proposed a novel dynamic interacting particle swarm optimization algorithm. The algorithm compares well to a standard PSO implementation, especially in terms of efficient solution of high-dimensional problems containing few local minimizers. Based on these promising results, the DYN-PSO algorithm deserves further development.

A number of outstanding issues remain. The importance of recomputing appropriate time-step sizes is already recognized, but more refined criteria should be developed. A guideline for the number of particles for efficient search is already proposed, but since it is based on a very narrow test, additional experiments are required. Also, the convergence criterion can be refined in order to also work for functions containing a very large number of local minimizers. The mechanism that probably requires the most attention is the energy dissipation scheme, which can be modified to increase the probability of convergence to the global minimum in the case of multiple local minima. This could be achieved by less aggressive energy dissipation during “uphill” moves, but will necessarily retard convergence.

The initial results indicate that the proposed DYN-PSO algorithm shows much promise as an alternative direct search method for solving large scale unconstrained optimization problems. The algorithm seems capable of solving unimodal problems economically, and it also has competitive performance on functions containing many local minima.

References

1. K. W. Kolodziej and J. Hjelm, Local Positioning Systems: LBS Applications and Services, CRC Press, Boca Raton, Fla, USA, 2006.
2. A. LaMarca, Y. Chawathe, S. Consolvo, et al., “Place lab: device positioning using radio beacons in the wild,” in Proceedings of the 3rd International Conference on Pervasive Computing (PERVASIVE '05), vol. 3468, pp. 116–133, Munich, Germany, May 2005.
3. J. Letchner, D. Fox, and A. LaMarca, “Large-scale localization from wireless signal strength,” in Proceedings of the 20th National Conference on Artificial Intelligence and the 17th Innovative Applications of Artificial Intelligence Conference (AAAI '05), pp. 15–20, Pittsburgh, Pa, USA, July 2005.
4. H. Lim, L.-C. Kung, J. C. Hou, and H. Luo, “Zero-configuration, robust indoor localization: theory and experimentation,” in Proceedings of the 25th IEEE International Conference on Computer Communications (INFOCOM '06), pp. 1–12, Barcelona, Spain, April 2006.
5. P. Bahl and V. N. Padmanabhan, “RADAR: an in-building RF-based user location and tracking system,” in Proceedings of the 19th Annual Joint Conference of the IEEE Computer and Communications Societies (INFOCOM '00), vol. 2, pp. 775–784, Tel Aviv, Israel, March 2000.
6. K. Kaemarungsi and P. Krishnamurthy, “Properties of indoor received signal strength for WLAN location fingerprinting,” in Proceedings of the 1st Annual International Conference on Mobile and Ubiquitous Systems: Networking and Services (MOBIQUITOUS '04), pp. 14–23, Boston, Mass, USA, August 2004.
7. D. Niculescu and B. Nath, “Ad hoc positioning system (APS),” in Proceedings of the IEEE Global Telecommunications Conference (GLOBECOM '01), vol. 5, pp. 2926–2931, San Antonio, Tex, USA, November 2001.
8. D. Niculescu and B. Nath, “Ad hoc positioning system (APS) using AOA,” in Proceedings of the 22nd Annual Joint Conference of the IEEE Computer and Communications Societies (INFOCOM '03), vol. 3, pp. 1734–1743, San Francisco, Calif, USA, March-April 2003.
9. R. Peng and M. L. Sichitiu, “Angle of arrival localization for wireless sensor networks,” in Proceedings of the 3rd Annual IEEE Communications Society Conference on Sensor, Mesh and Ad Hoc Communications and Networks (SECON '06), vol. 1, pp. 374–382, Reston, Va, USA, September 2006.
10. R. L. Moses, D. Krishnamurthy, and R. M. Patterson, “A self-localization method for wireless sensor networks,” EURASIP Journal on Applied Signal Processing, vol. 2003, no. 4, pp. 348–358, 2003.
11. D. Niculescu and B. Nath, “Position and orientation in ad hoc networks,” Ad Hoc Networks, vol. 2, no. 2, pp. 133–151, 2004.
12. N. Bulusu, J. Heidemann, and D. Estrin, “GPS-less low-cost outdoor localization for very small devices,” IEEE Personal Communications, vol. 7, no. 5, pp. 28–34, 2000.
13. C. Savarese, J. M. Rabaey, and K. Langendoen, “Robust positioning algorithms for distributed Ad-Hoc wireless sensor networks,” in Proceedings of USENIX Annual Technical Conference, pp. 317–327, Monterey, Calif, USA, June 2002.