Python Optimization with PSO
Let me start with a little joke:
“Did you know that before the invention of clocks, people had to actively walk around and ask the time?”
This simple anecdote illustrates an important concept: information available to one member of a group can be spread to others. This idea has deep implications and has applications in many areas.
Let us consider self-organizing systems in nature, for example, flocks of birds or fish. Let us imagine such a system as a set of particles, where each individual is a separate particle. We can assume that the movement of each particle in space is determined by two main factors:
Individually optimal position: what the individual considers best for himself.
AND globally optimal position: determined by the collective interaction of particles, a kind of “instruction” received by an individual from the “group leader”.
This raises a natural question: what is considered “optimal” in nature? What is best for an individual and for the entire group? Not being a biologist, I cannot answer these questions. However, by observing such behavior in nature, we can develop an effective optimization algorithm. In other words, having defined the criteria for “optimality”, we can apply this evolutionary approach to optimize a given function.
This algorithm is known as Particle Swarm Optimization (PSO). This may sound a bit complicated. What is meant by “optimization”? What is the role of mathematics in this process? What exactly is being optimized? In this article, I will try to explain all these points in detail. Moreover, we will apply OOP in Python to create our own ParticleSwarmOptimizer() class. And thus, we will go from the theoretical foundations of PSO to their practical implementation.
So, let's get started! Enjoy reading.
0. Introduction to the concept of “optimization”
If you are interested in the “particle swarm optimization” method, you are probably already familiar with the general concept of optimization. However, the term “particle swarm” may be new to you. Therefore, in order not to overload you with unnecessary information, I will not dwell on the basic concepts of optimization in detail.
Let's clarify the problem statement. As we already mentioned, what is meant by “the best solution for the swarm” from a mathematical point of view?
Mathematically, the “swarm birds” (or “particles”) are points in a K-dimensional space. In the case of real physical space, our “search space” could be the familiar three-dimensional space, but in general the dimensions can be much larger. It is important to note that this method is applicable not only to the example with birds: we can use it to find the minimum or maximum of any given function.
Let's look at an example: let's imagine that the price of a house is determined automatically by an algorithm (although relying entirely on an algorithm in such a matter would be extremely risky):
Initially it may not be entirely clear what is meant by “Home”. In practice we are dealing with a list of characteristics such as location, size, etc.
Now the picture becomes clearer. Each house can be thought of as a point in a multidimensional feature space. In other words, each house is a point in a k-dimensional space, where k is the number of features (location, size, etc.). The process of converting these features into a house price can be thought of as a “black box”:
The key question we ask ourselves is:
“What values of features (x) should a house have for its cost to be minimal?”
This is precisely the essence of our optimization task.
Ultimately, the particle swarm optimization method will help us find the x coordinatesminthat is, those feature values that will provide the minimum cost. Next, we will consider in detail how this works.
1. Introduction to Particle Swarm Optimization
Before going into details, I recommend you read the original articleThis is one of the most accessible scientific articles I have come across: it is written in simple and understandable language, but at the same time it is informative and engaging.
So, what is the basic idea of the method? As mentioned, we create a set of particles that interact and exchange information to find the global minimum of the objective function. If we draw an analogy, it is similar to the TV series “Stranger Things”, where the characters unite and use all available means of communication to defeat a common enemy.
In terms of implementation, we start with a given number of particles (num_particles). Each particle is assigned a random initial position in the feature space. For each position, the value of the objective function L is calculated. For example, particle 1 will have coordinates x1 and the value of the function L(x1). This is the zero iteration, after which our system is ready for evolution.
How does this evolution happen? With a parameter called “velocity”. For each particle x, each dimension (i) and each iteration (j), we define the velocity as follows:
Or, in more compact vector form:
Now let's look at how the velocity vector v is determined. This is where the elegance of the particle swarm optimization method comes into play. The vector v consists of three main components:
Let's look at each component in more detail:
vinertial: This component represents the “memory” of the particle's previous velocity. The term is borrowed from physics, where an inertial frame of reference is a system in which a body, not subject to external forces, maintains a state of rest or uniform rectilinear motion.
where kinertial – this is a constant.
vcognitive: This component represents the speed offered by the particle itself (hence the name “cognitive”). It guides the particle towards its best known position, i.e. where the objective function value was minimal for that particular particle. To add randomness to the search process, we also use a random number r1:
kcognitive – is also a constant xbest– these are the coordinates of the best known position for particle x, that is, the coordinates where the value of the objective function was minimal for this particle.
vsocial: This component is responsible for taking into account information from all other particles in the swarm. It directs the particle towards the globally best position, that is, where the value of the objective function was minimal for the entire swarm.
ksocial – is a constant, r2– a random number, and xglobal best – the coordinates of the point in the feature space where the objective function reaches its global minimum (the smallest value for the entire swarm).
This approach is actually quite compact in my opinion. Let's summarize:
We start by generating a given number of particles and randomly placing them in the feature space.
Particles move in feature space based on their velocity vectors.
The speed of each particle is determined by three components: inertia, cognitive and social components.
We repeat the steps of movement and velocity calculation for a given number of iterations (num_iteration).
At the end of the iterations, we select the particle that achieved the best result (the minimum value of the objective function).
That's it for the theoretical part. Let's move on to practice and write some Python code!
2. Implementation of PSO in practice
To implement the PSO algorithm, we will create two classes. The first class will describe a single swarm particle, and the second, using the first as a basis, will be responsible for creating and managing the entire swarm.
First, let's create a file constants.py. In this file we will store all default parameters, such as the boundaries of the search area for particles, the number of algorithm iterations and the values of the k coefficients (inertia, cognitive and social components).
How to change these parameters? To do this, you need to create a file pso_config.json in JSON format and place it in a folder config. In this file you can set your own values for the parameters. If the file pso_config.json is missing, the script we will write for the classes will create it automatically using default parameters.
NUM_PARTICLES = 30
NUM_ITERATIONS = 100
INERTIA = 0.5
COGNITIVE = 2
SOCIAL = 2
BOUNDS = [[-10,10],[-10,10]]
CONFIG_FILENAME = 'config/pso_config.json'
DEFAULT_CONFIG_NAMES = ["num_particles","num_iterations","inertia","cognitive","social","bounds"]
DEFAULT_PARTICLE = [NUM_PARTICLES, NUM_ITERATIONS,INERTIA, COGNITIVE, SOCIAL, BOUNDS]Now everything we discussed is implemented in the code, which is located in the file particle_swarm_optimizer.py.
import numpy as np
import json
from constants import * # Импортируем константы из файла constants.py
import os
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np
class Particle:
def __init__(self, bounds):
# Инициализируем позицию частицы случайным образом с учетом границ
self.position = np.array([np.random.uniform(low, high) for low, high in bounds])
# Инициализируем скорость частицы нулями
self.velocity = np.zeros_like(self.position)
# Инициализируем лучшую позицию частицы текущей позицией
self.best_position = np.copy(self.position)
# Инициализируем лучшее значение целевой функции бесконечностью
self.best_value = float('inf')
def update_velocity(self, global_best_position, inertia, cognitive, social):
# Генерируем случайные числа для когнитивной и социальной составляющих скорости
r1 = np.random.random(len(self.position))
r2 = np.random.random(len(self.position))
# Вычисляем когнитивную составляющую скорости
cognitive_velocity = cognitive * r1 * (self.best_position - self.position)
# Вычисляем социальную составляющую скорости
social_velocity = social * r2 * (global_best_position - self.position)
# Обновляем скорость частицы
self.velocity = inertia * self.velocity + cognitive_velocity + social_velocity
def update_position(self, bounds):
# Обновляем позицию частицы
self.position += self.velocity
# Проверяем, не вышла ли частица за границы области поиска
for i in range(len(self.position)):
if self.position[i] < bounds[i][0]:
self.position[i] = bounds[i][0]
if self.position[i] > bounds[i][1]:
self.position[i] = bounds[i][1]
class ParticleSwarmOptimizer:
def __init__(self, particle_class, objective_function, config_filename = None):
self.particle_class = particle_class # Класс частицы, который будет использоваться
if config_filename is None:
self.config = self.load_or_create_config()
else:
self.config = self.load_or_create_config(config_filename)
self.objective_function = objective_function # Целевая функция, которую нужно оптимизировать
self.bounds = self.config['bounds'] # Границы области поиска
self.num_particles = self.config['num_particles'] # Количество частиц в рое
self.num_iterations = self.config['num_iterations'] # Количество итераций алгоритма
self.inertia = self.config['inertia'] # Коэффициент инерции
self.cognitive = self.config['cognitive'] # Коэффициент когнитивной составляющей
self.social = self.config['social'] # Коэффициент социальной составляющей
# Создаем рой частиц
self.swarm = [self.particle_class(self.bounds) for _ in range(self.num_particles)]
# Инициализируем глобально лучшую позицию и значение
self.global_best_position = np.copy(self.swarm[0].position)
self.global_best_value = float('inf')
def load_or_create_config(self,filename=CONFIG_FILENAME):
if os.path.exists(filename):
# Если файл конфигурации существует, загружаем его
with open(filename, 'r') as f:
config = json.load(f)
print("Loaded configuration from", filename)
else:
# Если файл конфигурации не найден, создаем его с значениями по умолчанию
with open(filename, 'w') as f:
default_config = {DEFAULT_CONFIG_NAMES[i]:DEFAULT_PARTICLE[i] for i in range(len(DEFAULT_PARTICLE))}
json.dump(default_config, f, indent=4)
config = default_config
print(f"Configuration file not found. Created default config at {filename}")
return config
def optimize(self, save_iterations=False, print_iterations=False):
# Инициализируем списки для хранения лучших значений и позиций на каждой итерации, если требуется
if save_iterations:
best_values_per_iteration = np.zeros(self.num_iterations)
best_positions_per_iteration = np.zeros((self.num_iterations, len(self.bounds)))
for iteration in range(self.num_iterations):
for particle in self.swarm:
# Вычисляем значение целевой функции для текущей позиции частицы
value = self.objective_function(particle.position)
# Обновляем личный лучший результат частицы
if value < particle.best_value:
particle.best_value = value
particle.best_position = np.copy(particle.position)
# Обновляем глобально лучший результат
if value < self.global_best_value:
self.global_best_value = value
self.global_best_position = np.copy(particle.position)
# Сохраняем текущее глобально лучшее значение и позицию для этой итерации, если требуется
if save_iterations:
best_values_per_iteration[iteration] = self.global_best_value
best_positions_per_iteration[iteration] = np.copy(self.global_best_position)
# Обновляем скорость и позицию каждой частицы
for particle in self.swarm:
particle.update_velocity(self.global_best_position, self.inertia, self.cognitive, self.social)
particle.update_position(self.bounds)
# Выводим информацию о текущей итерации, если требуется
if print_iterations:
print(f"Iteration {iteration+1}/{self.num_iterations}, Global Best Value: {self.global_best_value}")
# Возвращаем лучшие значения и позиции для каждой итерации, если save_iterations=True
if save_iterations:
self.best_values_per_iteration = best_values_per_iteration
self.best_positions_per_iteration = best_positions_per_iteration
return self.global_best_position, self.global_best_value, self.best_values_per_iteration, self.best_positions_per_iteration
else:
return self.global_best_position, self.global_best_value, None, None
def plot_pso_convergence(self, animated=False, gif_name="pso_convergence.gif"):
"""
Строит график сходимости алгоритма PSO с контурным графиком целевой функции.
Параметры:
- objective_function: целевая функция для построения контурного графика.
- animated: если True, создает анимированный GIF, показывающий процесс сходимости с движением частиц.
- gif_name: имя GIF-файла для сохранения, если animated=True.
"""
objective_function = self.objective_function
iteration_values = self.best_values_per_iteration
iteration_positions = self.best_positions_per_iteration
# Создаем сетку точек для вычисления целевой функции для построения контурного графика
x = np.linspace(self.bounds[0][0], self.bounds[0][1], 100)
y = np.linspace(self.bounds[1][0], self.bounds[1][1], 100)
X, Y = np.meshgrid(x, y)
Z = np.array([[objective_function([x_val, y_val]) for x_val in x] for y_val in y])
if animated:
# Создаем анимированный GIF с движением частиц и глобально лучшей позицией
fig, ax = plt.subplots()
# Строим контурный график целевой функции
contour = ax.contourf(X, Y, Z, levels=50, cmap='viridis')
plt.colorbar(contour, ax=ax)
# Инициализируем точечный график для позиций частиц
particles = ax.scatter([], [], c="blue", label="Particles", s=50)
global_best = ax.scatter([], [], c="red", label="Global Best", s=100, marker="x")
ax.set_title("PSO Particle Movement and Convergence")
ax.set_xlabel("X Position")
ax.set_ylabel("Y Position")
def init():
particles.set_offsets([])
global_best.set_offsets([])
return particles, global_best
def update(frame):
# Обновляем позиции частиц (в этом случае только глобально лучшую позицию)
particle_positions = iteration_positions[frame]
particles.set_offsets(particle_positions)
# Обновляем глобально лучшую позицию (красный крестик)
global_best_position = iteration_positions[frame]
global_best.set_offsets(global_best_position)
return particles, global_best
ani = animation.FuncAnimation(fig, update, frames=len(iteration_positions), init_func=init, blit=True)
# Сохраняем анимацию как GIF
ani.save(gif_name, writer="imagemagick", fps=5)
print(f"Animation saved as {gif_name}")
else:
# Статический график с контуром и сходимостью глобально лучшего значения по итерациям
fig, ax = plt.subplots()
# Строим контурный график целевой функции
contour = ax.contourf(X, Y, Z, levels=50, cmap='viridis')
plt.colorbar(contour, ax=ax)
# Строим график глобально лучших значений по итерациям
plt.plot(iteration_values, marker="o", linestyle="-", color="b")
plt.title("PSO Convergence Plot")
plt.xlabel("Iterations")
plt.ylabel("Global Best Value")
plt.grid(True)
plt.show()Let's take a quick look at the structure of the implemented code:
Particle() class:
This class is responsible for creating and describing individual swarm particles. Each particle is initialized with a random position in the feature space. The class also contains methods update_velocity() And update_position()which update the particle's velocity and position at each iteration of the algorithm according to the rules described above.
ParticleSwarmOptimizer() class:
This class is responsible for creating and managing the entire particle swarm. When creating an object of this class, two arguments must be passed:
Particle: Used to create particle objects. You can modify or extend the class
Particle()to change the behavior of the algorithm or add new features.objective_function: This is the function we are trying to optimize (find its minimum). It is a “black box” that takes as input the coordinates of a particle in the feature space and returns the value of the objective function at that point.
The main logic of the PSO algorithm is implemented in the class ParticleSwarmOptimizer(). Method optimize() starts the optimization process. You can control the output of information during optimization using the parameters print_iterations (whether to output information about each iteration) and save_iterations (whether to keep history of best values and positions at each iteration). Method plot_pso_convergence() is responsible for visualizing the algorithm convergence process. If you set the parameter animated=Truethe method will create a GIF animation that clearly demonstrates the movement of particles and the change in the position of the global minimum at each iteration (we will see an example shortly).
3. Practical examples
We have done all the necessary preparatory work and can now run the PSO algorithm with just a few lines of code.
Let's start with a simple quadratic function:
import numpy as np
import json
import os
from particle_swarm_optimizer import *
def quadratic_objective_function(x):
return x[0]**2 + x[1]**2 + 50 # Простая квадратичная функция
pso = ParticleSwarmOptimizer(Particle,quadratic_objective_function)
best_position, best_value, iteration_values,iteration_positions = pso.optimize(save_iterations=True,print_iterations=False)
print("Best Position:", best_position)
print("Best Value:", best_value)
pso.plot_pso_convergence(animated=True)Here is the gif we get:
The PSO algorithm copes well even with fairly complex objective functions, for example:
Running the same code as in the previous example, we get:
def complex_objective_function(x):
# Комбинация квадратичной функции, синусоиды и экспоненциального затухания
term1 = (x[0]**2 + x[1]**2) # Квадратичная функция
term2 = 10 * np.sin(x[0]) * np.sin(4*x[1]) # Синусоида, вносящая нелинейность
term3 = np.exp(-0.1 * (x[0]**2 + x[1]**2)) # Экспоненциальное затухание для создания локальных минимумов
return term1 - term2 + 5 * term3 + 20 # 20 просто сдвигает всю поверхность вверх
# Создаем объект оптимизатора PSO
pso = ParticleSwarmOptimizer(Particle,complex_objective_function)
# Запускаем оптимизацию и сохраняем значения на каждой итерации
best_position, best_value, iteration_values,iteration_positions = pso.optimize(save_iterations=True,print_iterations=False)
# Выводим результаты
print("Best Position:", best_position)
print("Best Value:", best_value)
# Строим график сходимости с анимацией
pso.plot_pso_convergence(animated=True)Our gif:
We have verified that the PSO algorithm not only efficiently finds the minimum for functions with one clearly defined minimum point, but is also capable of finding the global minimum (at least we can assume that it is the global minimum, since it is analytically difficult to verify) even for complex functions with many local minima.
4. Conclusions
Thank you for taking the time to read this article and diving into the fascinating world of PSO optimization with me. I hope I have managed to explain the main ideas of the algorithm in an accessible and understandable way.
In this article we:
We considered an approach to solving optimization problems based on the analysis of information received from a set of particles (agents).
We illustrated the problem of “optimization” using a simple and clear example of finding the minimum price for a house.
We got acquainted with the basic principles of particle swarm optimization (PSO) and its key idea.
We studied the mathematical foundations of the PSO algorithm and traced how it evolves from an initial random state of a swarm consisting of num_particles particles to finding the best estimate of the global minimum of the objective function.
Implemented the PSO algorithm from scratch using OOP in Python.
We tested our implementation on a simple quadratic function and on a more complex function with many local minima. In both cases, the algorithm demonstrated high efficiency and accuracy!
