Начну с небольшой шутки:

"Знаете ли вы, что до изобретения часов людям приходилось активно ходить повсюду и спрашивать время?"

Этот незамысловатый анекдот иллюстрирует важную концепцию: информация, доступная одному члену группы, может распространяться среди остальных. Эта идея имеет глубокий смысл и находит применение во многих областях.

Рассмотрим самоорганизующиеся системы в природе, например, стаи птиц или рыб. Представим такую систему как совокупность частиц, где каждая особь – это отдельная частица. Можно предположить, что движение каждой частицы в пространстве определяется двумя основными факторами:

Индивидуально оптимальная позиция: то, что особь считает наилучшим для себя.

И глобально оптимальная позиция: определяемая коллективным взаимодействием частиц, своего рода "инструкция", получаемая особью от "лидера группы".

В связи с этим возникает естественный вопрос: что считать "оптимальным" в природе? Что является наилучшим для отдельной особи и для всей группы? Не являясь биологом, я не могу дать ответы на эти вопросы. Однако, наблюдая за подобным поведением в природе, мы можем разработать эффективный алгоритм оптимизации. Другими словами, определив критерии "оптимальности", мы можем применить этот эволюционный подход для оптимизации заданной функции.

Данный алгоритм известен как оптимизация роем частиц (Particle Swarm Optimization, PSO). Возможно, это звучит несколько сложно. Что подразумевается под "оптимизацией"? Какова роль математики в этом процессе? Что именно оптимизируется? В статье я постараюсь подробно разъяснить все эти моменты. Более того, мы применим ООП на Python для создания собственного класса ParticleSwarmOptimizer(). И таким образом, мы пройдем путь от теоретических основ PSO до их практической реализации.

Итак, приступим! Желаю приятного чтения.

0. Введение в понятие "оптимизация"

Если вас заинтересовал метод "оптимизации роем частиц", вероятно, вы уже знакомы с общим понятием оптимизации. Однако, возможно, термин "рой частиц" для вас нов. Поэтому, чтобы не перегружать вас лишней информацией, я не буду подробно останавливаться на базовых понятиях оптимизации.

Давайте все же уточним постановку задачи. Как мы уже упомянули, что с математической точки зрения подразумевается под "наилучшим решением для роя"?

С математической точки зрения, "птицы роя" (или "частицы") представляют собой точки в K-мерном пространстве. В случае реального физического пространства, нашей "областью поиска" может быть привычное нам трехмерное пространство, однако в общем случае размерность может быть значительно больше. Важно отметить, что данный метод применим не только к примеру с птицами: мы можем использовать его для нахождения минимума или максимума любой заданной функции.

Рассмотрим пример: представим, что цена дома определяется автоматически с помощью алгоритма (хотя полностью полагаться на алгоритм в таком вопросе было бы крайне рискованно):

Первоначально может быть не совсем понятно, что подразумевается под "Домом". На практике мы имеем дело со списком характеристик, таких как местоположение, размер и т.д.

Теперь картина становится яснее. Каждый дом можно представить как точку в многомерном пространстве признаков. Иными словами, каждый дом – это точка в k-мерном пространстве, где k – это количество характеристик (местоположение, размер и т.д.). Процесс преобразования этих характеристик в цену дома можно представить в виде "черного ящика":

Ключевой вопрос, который мы ставим перед собой:

"Какими значениями признаков (x) должен обладать дом, чтобы его стоимость была минимальной?"

Именно в этом и заключается суть нашей задачи оптимизации.

В конечном счете, метод оптимизации роем частиц поможет нам найти координаты x_min, то есть те значения признаков, которые обеспечат минимальную стоимость. Далее мы подробно рассмотрим, как это работает.

1. Знакомство с оптимизацией роем частиц

Прежде чем углубиться в детали, рекомендую ознакомиться с оригинальной статьей. Это одна из наиболее доступных научных статей, с которыми мне приходилось сталкиваться: она написана простым и понятным языком, но при этом информативна и увлекательна.

Итак, в чем заключается основная идея метода? Как уже упоминалось, мы создаем набор частиц, которые взаимодействуют и обмениваются информацией, чтобы найти глобальный минимум целевой функции. Если провести аналогию, это похоже на сериал "Очень странные дела", где герои объединяются и используют все доступные средства коммуникации, чтобы победить общего врага.

С точки зрения реализации, мы начинаем с заданного количества частиц (num_particles). Каждой частице присваивается случайное начальное положение в пространстве признаков. Для каждого положения рассчитывается значение целевой функции L. Например, частица 1 будет иметь координаты x₁ и значение функции L(x₁). Это нулевая итерация, после которой наша система готова к эволюции.

Как происходит эта эволюция? С помощью параметра, называемого "скоростью". Для каждой частицы x, каждого измерения (i) и каждой итерации (j) мы определяем скорость следующим образом:

Или, в более компактной векторной форме:

Теперь давайте разберемся, как определяется вектор скорости v. Именно здесь проявляется вся элегантность метода оптимизации роя частиц. Вектор v состоит из трех основных компонентов:

Рассмотрим каждый компонент подробнее:

v_inertial: Этот компонент представляет собой "память" о предыдущей скорости частицы. Термин заимствован из физики, где инерциальная система отсчета – это система, в которой тело, не подверженное воздействию внешних сил, сохраняет состояние покоя или равномерного прямолинейного движения.

где k_inertial – это константа.

v_cognitive: Этот компонент представляет собой скорость, предлагаемую самой частицей (отсюда и название "когнитивный"). Он направляет частицу в сторону ее лучшего известного положения, то есть туда, где значение целевой функции было минимальным для этой конкретной частицы. Для добавления случайности в процесс поиска мы также используем случайное число r₁:

k_cognitive – это также константа. x_best– это координаты лучшего известного положения для частицы x, то есть координаты, где значение целевой функции было минимальным для этой частицы.

v_social: Этот компонент отвечает за учет информации от всех остальных частиц роя. Он направляет частицу в сторону глобально лучшего положения, то есть туда, где значение целевой функции было минимальным для всего роя.

k_social – это константа, r₂– случайное число, а x_{global best} – координаты точки в пространстве признаков, где целевая функция достигает своего глобального минимума (наименьшего значения для всего роя).

Этот подход, на мой взгляд, действительно довольно компактен. Давайте подведем итог:

• Мы начинаем с создания заданного количества частиц и случайного размещения их в пространстве признаков.

• Частицы перемещаются в пространстве признаков, основываясь на своих векторах скорости.

• Скорость каждой частицы определяется тремя составляющими: инерцией, когнитивной и социальной составляющими.

• Мы повторяем шаги движения и расчета скорости на протяжении заданного количества итераций (num_iteration).

• По окончании итераций мы выбираем частицу, которая достигла наилучшего результата (минимального значения целевой функции).

На этом теоретическая часть завершена. Давайте перейдем к практике и напишем код на Python!

2. Реализация PSO на практике

Для реализации алгоритма PSO мы создадим два класса. Первый класс будет описывать отдельную частицу роя, а второй, используя первый как основу, будет отвечать за создание и управление всем роем.

Для начала создадим файл constants.py. В этом файле мы будем хранить все параметры по умолчанию, такие как границы области поиска для частиц, количество итераций алгоритма и значения коэффициентов k (инерции, когнитивной и социальной составляющих).

Как изменить эти параметры? Для этого необходимо создать файл pso_config.json в формате JSON и поместить его в папку config. В этом файле вы сможете задать собственные значения для параметров. Если файл pso_config.json отсутствует, скрипт, который мы напишем для классов, создаст его автоматически, используя параметры по умолчанию.

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]

Теперь все, что мы обсуждали, реализовано в коде, который находится в файле 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()

Давайте кратко рассмотрим структуру реализованного кода:

Класс Particle():

Этот класс отвечает за создание и описание отдельных частиц роя. Каждая частица инициализируется случайным положением в пространстве признаков. Класс также содержит методы update_velocity() и update_position(), которые обновляют скорость и положение частицы на каждой итерации алгоритма в соответствии с описанными выше правилами.

Класс ParticleSwarmOptimizer():

Этот класс отвечает за создание и управление всем роем частиц. При создании объекта этого класса необходимо передать два аргумента:

1. Particle: Используется для создания объектов-частиц. Вы можете модифицировать или расширить класс Particle(), чтобы изменить поведение алгоритма или добавить новые функции.

2. objective_function: Это функция, которую мы стремимся оптимизировать (найти ее минимум). Она представляет собой "черный ящик", который принимает на вход координаты частицы в пространстве признаков и возвращает значение целевой функции в этой точке.

Основная логика алгоритма PSO реализована в классе ParticleSwarmOptimizer(). Метод optimize() запускает процесс оптимизации. Вы можете управлять выводом информации во время оптимизации с помощью параметров print_iterations (выводить ли информацию о каждой итерации) и save_iterations (сохранять ли историю лучших значений и позиций на каждой итерации). Метод plot_pso_convergence() отвечает за визуализацию процесса сходимости алгоритма. Если установить параметр animated=True, метод создаст GIF-анимацию, наглядно демонстрирующую движение частиц и изменение положения глобального минимума на каждой итерации (вскоре мы увидим пример).

3. Практические примеры

Мы проделали всю необходимую подготовительную работу и теперь можем запустить алгоритм PSO всего несколькими строками кода.

Для начала рассмотрим простую квадратичную функцию:

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)

Вот гиф, которую мы получим:

Алгоритм PSO отлично справляется даже с достаточно сложными целевыми функциями, например:

Запустив тот же код, что и в предыдущем примере, мы получим:

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)

Наша гифка:

Мы убедились, что алгоритм PSO не только эффективно находит минимум для функций с одной ярко выраженной точкой минимума, но и способен находить глобальный минимум (по крайней мере, мы можем предположить, что это именно глобальный минимум, так как аналитически это сложно проверить) даже для сложных функций с множеством локальных минимумов.

4. Выводы

Благодарю вас за то, что уделили время прочтению этой статьи и вместе со мной погрузились в увлекательный мир оптимизации с помощью PSO. Надеюсь, мне удалось доступно и понятно изложить основные идеи алгоритма.

В рамках статьи мы:

• Рассмотрели подход к решению задач оптимизации, основанный на анализе информации, получаемой от множества частиц (агентов).

• Проиллюстрировали задачу "оптимизации" на простом и наглядном примере поиска минимальной цены дома.

• Ознакомились с основными принципами оптимизации роем частиц (PSO) и ее ключевой идеей.

• Изучили математические основы алгоритма PSO и проследили, как он эволюционирует от начального случайного состояния роя, состоящего из num_particles частиц, к нахождению наилучшей оценки глобального минимума целевой функции.

• Реализовали алгоритм PSO с нуля, используя ООП на Python.

• Протестировали нашу реализацию на простой квадратичной функции и на более сложной функции с множеством локальных минимумов. В обоих случаях алгоритм продемонстрировал высокую эффективность и точность!