particle swarm algorithm

1. Introduction to the algorithm

The idea of ​​particle swarm optimization originated from the research on the predation behavior of birds, simulating the behavior of bird flocks flying and foraging.
Imagine such a scenario, a group of birds are randomly searching for food, there is a piece of food in a certain area, all the birds do not know where the food is, but they can feel how far the current position is from the food, at this time What is the optimal strategy for finding food? The answer is: search the area around the bird closest to the food, and judge where the food is based on your own flying experience.

PSO is inspired by this model:
1. The solution to each optimization problem is imagined as a bird, called a 'particle', and all particles are searched in a D-dimensional space.
2. All particles are determined by a fitness function to determine the fitness value of the current position.
3. Each particle is given a memory function, which can remember the best position searched for
4. Each particle also has a velocity to determine the distance and direction of flight. This speed is dynamically adjusted according to its own flying experience and the flying experience of its companions.
5. It can be seen that the basis of particle swarm algorithm lies in the social sharing of information

In the D-dimensional space, there are N particles;
Position of particle i: xi=(xi1,xi2,...xiD), substitute xi into the fitness function f(xi) to obtain the fitness value;
Particle i speed: vi=(vi1,vi2,...viD)
The best position experienced by the individual particle i: pbesti=(pi1,pi2,...piD)
The best position the population has experienced: gbest=(g1,g2,…gD)

Usually, the position variation range in the d-th (1≤d≤D) dimension is limited to ( X m i n , d , X m a x , d X_{min,d},X_{max,d} Xmin,d​,Xmax,d​),
The speed variation range is limited to ( − V m a x , d , V m a x , d -V_{max,d},V_{max,d} −Vmax,d​,Vmax,d​) (that is, if in the iteration
beyond the boundary value, the velocity or position of the dimension is limited to the maximum velocity or boundary of the dimension
Location)
Update formulas for speed and position
With the above understanding, a very key question is how to update the speed and position of the particle?

The d-dimensional velocity update formula of particle i is:

The d-dimensional position update formula of particle i

v i d k v_{id}^k vidk​: The d-th dimension component of the flight velocity vector of the k-th iteration particle i
x i d k x_{id}^k xidk​: The d-th dimension component of the k-th iteration particle i position vector
c1,c2: acceleration constant, adjust the maximum learning step size (hyperparameter)
r1,r2: Two random functions, the value range [0,1], to increase randomness
w: inertia weight, non-negative, adjusts the search range of the solution space. (hyperparameter)
It can be seen that the pso algorithm contains a total of 3 hyperparameters.

For the particle velocity update formula, here are some more introductions:
Its velocity update formula consists of three parts:
The first part is the previous velocity of the particle
The second part is the "cognitive" part, which represents the thinking of the particle itself, which can be understood as the distance between the current position of particle i and its best position
The third part is the 'social part', which represents the information sharing and cooperation between particles, which can be understood as the distance between the current position of particle i and the best position of the group.
As shown in the figure below, the process of particle motion is affected by the above three aspects:

2. Algorithm process

  1. Initial:
    Initialize a particle swarm (swarm size n), including random positions and velocities.
  2. Evaluation:
    According to the fitness function, the fitness of each particle is evaluated.
  3. Find the Pbest:
    For each particle, its current fitness value is compared with the fitness value corresponding to its individual historical best position (pbest), and if the current fitness value is higher, the historical best position pbest will be updated with the current position.
  4. Find the Gbest:
    For each particle, compare its current fitness value with the fitness value corresponding to the global best position (gbest). If the current fitness value is higher, the global best position gbest will be updated with the current particle position.
  5. Update the Velocity:
    Update the velocity and position of each particle according to the formula.
  6. If the end condition is not met, go back to step 2
    Usually the algorithm stops when the maximum number of iterations is reached or the increment of the best fitness value is less than a given threshold.

The flow chart of the algorithm is as follows:

3. Algorithm example








4. Algorithm implementation

Taking the above example as an example, the implementation of the algorithm is as follows. If you need to optimize other problems, you only need to adjust the fitness function.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def fit_fun(x):  # fitness function
    return sum(100.0 * (x[0][1:] - x[0][:-1] ** 2.0) ** 2.0 + (1 - x[0][:-1]) ** 2.0)


class Particle:
    # initialization
    def __init__(self, x_max, max_vel, dim):
        self.__pos = np.random.uniform(-x_max, x_max, (1, dim))  # the position of the particle
        self.__vel = np.random.uniform(-max_vel, max_vel, (1, dim))  # the speed of the particles
        self.__bestPos = np.zeros((1, dim))  # best location for particles
        self.__fitnessValue = fit_fun(self.__pos)  # fitness function value

    def set_pos(self, value):
        self.__pos = value

    def get_pos(self):
        return self.__pos

    def set_best_pos(self, value):
        self.__bestPos = value

    def get_best_pos(self):
        return self.__bestPos

    def set_vel(self, value):
        self.__vel = value

    def get_vel(self):
        return self.__vel

    def set_fitness_value(self, value):
        self.__fitnessValue = value

    def get_fitness_value(self):
        return self.__fitnessValue


class PSO:
    def __init__(self, dim, size, iter_num, x_max, max_vel, tol, best_fitness_value=float('Inf'), C1=2, C2=2, W=1):
        self.C1 = C1
        self.C2 = C2
        self.W = W
        self.dim = dim  # particle size
        self.size = size  # Number of particles
        self.iter_num = iter_num  # number of iterations
        self.x_max = x_max
        self.max_vel = max_vel  # maximum particle speed
        self.tol = tol  # As of the conditions
        self.best_fitness_value = best_fitness_value
        self.best_position = np.zeros((1, dim))  # optimal location of the population
        self.fitness_val_list = []  # Optimal fitness value for each iteration

        # Initialize the population
        self.Particle_list = [Particle(self.x_max, self.max_vel, self.dim) for i in range(self.size)]

    def set_bestFitnessValue(self, value):
        self.best_fitness_value = value

    def get_bestFitnessValue(self):
        return self.best_fitness_value

    def set_bestPosition(self, value):
        self.best_position = value

    def get_bestPosition(self):
        return self.best_position

    # refresh rate
    def update_vel(self, part):
        vel_value = self.W * part.get_vel() + self.C1 * np.random.rand() * (part.get_best_pos() - part.get_pos()) \
                    + self.C2 * np.random.rand() * (self.get_bestPosition() - part.get_pos())
        vel_value[vel_value > self.max_vel] = self.max_vel
        vel_value[vel_value < -self.max_vel] = -self.max_vel
        part.set_vel(vel_value)

    # update location
    def update_pos(self, part):
        pos_value = part.get_pos() + part.get_vel()
        part.set_pos(pos_value)
        value = fit_fun(part.get_pos())
        if value < part.get_fitness_value():
            part.set_fitness_value(value)
            part.set_best_pos(pos_value)
        if value < self.get_bestFitnessValue():
            self.set_bestFitnessValue(value)
            self.set_bestPosition(pos_value)

    def update_ndim(self):

        for i in range(self.iter_num):
            for part in self.Particle_list:
                self.update_vel(part)  # refresh rate
                self.update_pos(part)  # update location
            self.fitness_val_list.append(self.get_bestFitnessValue())  # After each iteration, save the current optimal fitness to the list
            print('the first{}The next best fit is{}'.format(i, self.get_bestFitnessValue()))
            if self.get_bestFitnessValue() < self.tol:
                break

        return self.fitness_val_list, self.get_bestPosition()

if __name__ == '__main__':
    # test banana function
    pso = PSO(4, 5, 10000, 30, 60, 1e-4, C1=2, C2=2, W=1)
    fit_var_list, best_pos = pso.update_ndim()
    print("optimal location:" + str(best_pos))
    print("Optimal solution:" + str(fit_var_list[-1]))
    plt.plot(range(len(fit_var_list)), fit_var_list, alpha=0.5)

5. Algorithm application

Note that the algorithm needs to pay attention to the following points when solving specific problems:
1. Population size m
If m is small, it is easy to fall into local optimum. If m is large, the optimization ability of pso is very good. When the number of population increases to a certain level, further growth will no longer have a significant effect.
2. Weighting factor
Three parts for particle velocity update:
a. The inertia factor w=1 means the basic particle swarm algorithm, and w=0 means losing the speed memory of the particle itself.
b. The learning factor c1=0 of the self-awareness part means that the selfless particle swarm algorithm has only society and no self, which will make the group lose diversity, which will easily lead to falling into local optimum and unable to jump out.
c. The learning factor c2=0 of the social experience part represents the self-type particle swarm algorithm, only the self has no society, which leads to social sharing without information, and the algorithm converges slowly.
The choice of these three parameters is very important. How to adjust these three parameters so that the algorithm can avoid premature and converge relatively quickly is of great significance for solving practical problems.
3. Maximum speed
The role of the speed limit is to maintain the balance between the exploration ability and the development ability of the algorithm.
V m V_m When Vm​is large, the exploration ability is strong, but the particles are easy to fly over the optimal solution
V m V_m When Vm is small, the development ability is strong, but it is easy to fall into the local optimal solution
V m V_m Vm​ is generally set to 10%~20% of the variation range of each dimension variable
4. Stop Guidelines
a. Maximum number of iterations
b. Acceptable satisfactory solution (judging whether it is satisfactory or not through fitness function)
5. Initialization of particle space
A better choice of the initialization space of the particles will greatly shorten the convergence time. The initialization space varies according to the specific problem, and is set according to the specific problem. The setting of the few key parameters of the algorithm has a great effect on the accuracy and efficiency of the algorithm.
Significant impact.
6. Linearly decreasing weights (untested)

w m a x w_{max} wmax​Maximum inertia weight, w m i n w_{min} wmin​Minimum inertia weight, the current number of iterations of run, r u n m a x run_{max} runmax​ is the total number of algorithm iterations

Larger w has better global convergence ability, and smaller w has stronger local convergence ability. Therefore, with the increase of the number of iterations, the inertia weight w should be continuously reduced, so that the particle swarm algorithm has a strong global convergence ability in the early stage, and a strong local convergence ability in the late stage.
7. Shrinkage factor method (untested)

In 1999, Clerc introduced a shrinkage factor to ensure the convergence of the algorithm, that is, the speed update formula was changed to the above formula, where the shrinkage factor K is w limited by φ1 φ2 . φ1 and φ2 are model parameters that need to be set in advance.

The shrinkage factor method controls the final convergence of the system behavior, and can effectively search different regions, and this method can obtain higher quality solutions.
8. Neighborhood topology of the algorithm (untested)
The neighborhood topology structure of particle swarm optimization includes two kinds, one is that all individuals in the group are regarded as the neighborhood of the particle, and the other is that only some individuals in the group are regarded as the neighborhood of the particle. The neighborhood topology determines the historical optimal position of the group, so the algorithm is divided into global particle swarm optimization and local particle swarm optimization. What I implemented above is the global particle swarm optimization.
Global Particle Swarm Optimization
1. The particle's own historical optimal value
2. Global Optimum of Particle Swarm
local particle swarm algorithm
1. The historical optimal value of the particle itself
2. The optimal value of the particle in the particle neighborhood
The neighborhood grows linearly with the number of iterations, and finally the neighborhood expands to the entire particle swarm.
Practice has proved that the global version of the particle swarm algorithm has a fast convergence speed, but it is easy to fall into the local optimum. The local version of the particle swarm algorithm has a slow convergence speed, but it is difficult to fall into a local optimum. Most of the current particle swarm optimization algorithms focus on the two aspects of convergence speed and getting rid of local optimum. In fact, these two aspects are contradictory. See how a better compromise is made.

Posted by jkarr on Tue, 10 May 2022 02:33:07 +0300