Реализация алгоритма ПЧЕЛИНОЙ КОЛОНИИ на Python

Искусственная пчелиная колония (ABC) — стохастический метод поиска, основанный на роевом интеллекте, который мимикрирует рой медоносных пчёл, собирающих пищу. В этом алгоритме каждое потенциальное решение представлено положением источника пищи в пространстве поиска и качеством порции нектара, использованного для оценки целевой функции потерь.

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

Пример на Python решения задачи оптимизации с помощью пчелиного алгоритма

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

Один из них:

x = 3; y = 2

3D-представление значений функции Химмельблау

Блок-схема выполнения пчелиного алгоритма

Этапы решения задачи

Инициализация

1. Функции имеет следующий вид:

F = (x^2 + y - 11)^2 + (x + y^2 - 7)^2

Зададим границы поиска:

-4 ≤ x, y ≤4

1. Задаём количество пчёл = 100. Количество источников пищи равно количеству пчёл.
2. Инициализируем популяцию источников пищи (решения задачи) случайными позициями в пространстве поиска.
3. Рассчитываем функции потерь начальной популяции.

Фаза рабочей пчелы

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

v_ij = x_ij + phi_{ij}(x_ij - x_kj)

где

• v_ij– новый кандидат решения
• x_ij – текущее решение
• x_kj – случайно выбранное решение
• phi_{ij} – случайное число в диапазоне[-1, 1]

2. Рассчитываем функцию потерь для нового решения.
3. Применяем процесс greedy для выбора. Если новое решение лучше или равно текущему, заменяем им текущее. Иначе – сохраняем имеющееся значение.
4. Рассчитываем вероятности выбора: нормируем текущую функцию потерь на сумму значений этой функции для всех пчёл.

Фаза наблюдателя

1. Для каждой пчелы-наблюдателя выбираем решение, с вероятностью выбора, рассчитанной на предыдущем шаге.
2. Генерируем новый источник пищи вблизи выбранного решения.
3. Рассчитываем функцию потерь для нового кандидата на решение.
4. Применяем процесс greedy для выбора.

Фаза разведчика

1. Ищем брошенные решения. Если находим, заменяем на новое решение, сгенерированное случайным образом.
2. Запоминаем лучшее решение.

Критерий окончания

1. В данном случае критерием является достижение заданного максимального количества итераций.

Реализация пчелиного алгоритма на Python

Реализация выполнена на языке Python 3.11 в IDE MS Visual Studio Code.

import numpy as np
import matplotlib.pyplot as plt

def CostFunction(input):
    # Функция Химмельблау
    x = input[0]
    y = input[1]

    a = (x**2 + y - 11)
    b = (x + y**2 - 7)
    result = a**2 + b**2
    return result

# Инициализировать параметры

# Число и диапазон параметров в функции стоимости
nVar = 2
VarMin = -4
VarMax = 4

# Основные параметры алгоритма роя
iter_max = 40
nPop = 100
nOnLooker = 100
L = np.around(0.6*nVar*nPop)
a = 1

# Создать матрицы
PopPosition = np.zeros([nPop,nVar])
PopCost = np.zeros([nPop,1])
Probability = np.zeros([nPop,1])
BestSol = np.zeros([iter_max+1,nVar])
BestCost = np.inf*np.ones([iter_max+1,1])
Mine = np.zeros([nPop,1])

# Инициализировать местоположение источника меда
PopPosition = 8*np.random.rand(nPop,nVar) - 4
for i in range(nPop):
    PopCost[i][0] = CostFunction(PopPosition[i])
    if PopCost[i][0] < BestCost[0][0]:
        BestCost[0][0] = PopCost[i][0]
        BestSol[0] = PopPosition[i]

for iter in range(iter_max):

    # Фаза рабочей пчелы

    # Найти следующий источник меда
    for i in range(nPop):
        while True:
            k = np.random.randint(0,nPop)
            if k != i:
                break
        phi = a*(-1+2*np.random.rand(2))
        NewPosition = PopPosition[i] + phi*(PopPosition[i]-PopPosition[k])

        # Greedy выбор
        NewCost = CostFunction(NewPosition)
        if NewCost < PopCost[i][0]:
            PopPosition[i] = NewPosition
            PopCost[i][0] = NewCost
        else:
            Mine[i][0] = Mine[i][0]+1

    # Фаза пчелы-наблюдателя

    # Рассчитать матрицу вероятности выбора
    Mean = np.mean(PopCost)
    for i in range(nPop):
        Probability[i][0] = np.exp(-PopCost[i][0]/Mean)
    Probability = Probability/np.sum(Probability)
    CumProb = np.cumsum(Probability)

    for k in range(nOnLooker):

        # Dыбор методом  рулетки
        m = 0
        for i in range(nPop):
            m = m + CumProb[i]
            if m >= np.random.rand(1):
                break
       
        while True:
            k = np.random.randint(0,nPop)
            if k != i:
                break
        phi = a*(-1+2*np.random.rand(2))
        NewPosition = PopPosition[i] + phi*(PopPosition[i]-PopPosition[k])

        # Greedy выбор
        NewCost = CostFunction(NewPosition)
        if NewCost < PopCost[i][0]:
            PopPosition[i] = NewPosition
            PopCost[i][0] = NewCost
        else:
            Mine[i][0] = Mine[i][0]+1

    # Фаза пчелы-разведчика
    for i in range(nPop):
        if Mine[i][0] >= L:
            PopPosition[i] = 8*np.random.rand(1,nVar) - 4
            PopCost[i][0] = CostFunction(PopPosition[i])
            Mine[i][0] = 0

    # Сохранить оптимальное решение
    for i in range(nPop):
        if PopCost[i][0] < BestCost[iter+1][0]:
            BestCost[iter+1][0] = PopCost[i][0]
            BestSol[iter+1] = PopPosition[i]

# Выходной результат
y = np.zeros(iter_max+1)

for i in range(iter_max):
#    if i % 5 == 0:
    print(f'Итерация {i+1} Позиция: {BestSol[i]}   Ошибка: {BestCost[i]}')
    y[i] = BestCost[i][0]

print(f'Лучшее Решение:  {BestSol[iter_max-1]}')
x = [i for i in range(iter_max+1)]
plt.plot(x,y)
plt.show()

Результат выполнения программы

Итерация 1 Позиция: [-3.6317594  -3.08142279]   Ошибка: [2.08705542]
Итерация 2 Позиция: [2.89598819 1.93412458]   Ошибка: [0.59310992]
Итерация 3 Позиция: [3.04432297 1.99628027]   Ошибка: [0.07066022]
Итерация 4 Позиция: [3.01764824 1.96506991]   Ошибка: [0.01968474]
Итерация 5 Позиция: [2.98874841 2.01226043]   Ошибка: [0.00447797]
Итерация 6 Позиция: [2.98874841 2.01226043]   Ошибка: [0.00447797]
Итерация 7 Позиция: [2.98874841 2.01226043]   Ошибка: [0.00447797]
Итерация 8 Позиция: [2.98874841 2.01226043]   Ошибка: [0.00447797]
Итерация 9 Позиция: [2.98874841 2.01226043]   Ошибка: [0.00447797]
Итерация 10 Позиция: [2.98874841 2.01226043]   Ошибка: [0.00447797]
Итерация 11 Позиция: [3.00559953 1.98379558]   Ошибка: [0.00377936]
Итерация 12 Позиция: [3.00559953 1.98379558]   Ошибка: [0.00377936]
Итерация 13 Позиция: [2.99650443 1.98814883]   Ошибка: [0.00365321]
Итерация 14 Позиция: [2.99650443 1.98814883]   Ошибка: [0.00365321]
Итерация 15 Позиция: [-2.812734    3.12982526]   Ошибка: [0.00199206]
Итерация 16 Позиция: [2.99939245 2.00857787]   Ошибка: [0.00116525]
Итерация 17 Позиция: [-2.812734    3.12982526]   Ошибка: [0.00199206]
Итерация 18 Позиция: [2.98939398 2.02250803]   Ошибка: [0.00807152]
Итерация 19 Позиция: [2.98939398 2.02250803]   Ошибка: [0.00807152]
Итерация 20 Позиция: [2.98939398 2.02250803]   Ошибка: [0.00807152]
Итерация 21 Позиция: [2.98939398 2.02250803]   Ошибка: [0.00807152]
Итерация 22 Позиция: [3.0010804  2.01053109]   Ошибка: [0.00216574]
Итерация 23 Позиция: [3.0010804  2.01053109]   Ошибка: [0.00216574]
Итерация 24 Позиция: [3.0010804  2.01053109]   Ошибка: [0.00216574]
Итерация 25 Позиция: [3.01306056 1.99177696]   Ошибка: [0.00533424]
Итерация 26 Позиция: [3.00839208 1.98353693]   Ошибка: [0.00442387]
Итерация 27 Позиция: [-3.77812268 -3.28543931]   Ошибка: [0.00038172]
Итерация 28 Позиция: [-3.77812268 -3.28543931]   Ошибка: [0.00038172]
Итерация 29 Позиция: [-3.77812268 -3.28543931]   Ошибка: [0.00038172]
Итерация 30 Позиция: [-3.77812268 -3.28543931]   Ошибка: [0.00038172]
Итерация 31 Позиция: [3.00275904 2.00160502]   Ошибка: [0.00041434]
Итерация 32 Позиция: [3.00275904 2.00160502]   Ошибка: [0.00041434]
Итерация 33 Позиция: [2.99654367 2.00451058]   Ошибка: [0.00047629]
Итерация 34 Позиция: [2.99654367 2.00451058]   Ошибка: [0.00047629]
Итерация 35 Позиция: [2.99654367 2.00451058]   Ошибка: [0.00047629]
Итерация 36 Позиция: [2.99654367 2.00451058]   Ошибка: [0.00047629]
Итерация 37 Позиция: [2.999689   1.99640036]   Ошибка: [0.00024586]
Итерация 38 Позиция: [2.999689   1.99640036]   Ошибка: [0.00024586]
Итерация 39 Позиция: [2.99733045 2.00134727]   Ошибка: [0.00022241]
Итерация 40 Позиция: [2.99733045 2.00134727]   Ошибка: [0.00022241]

Лучшее Решение:  [2.99733045 2.00134727]

График уменьшения ошибки решения

Итог выполнения программы нахождения одного из минимумов функции Химмельблау

Как видно из журнала выполнения программы и графика падения значения ошибки, алгоритм быстро сходится к значениям Х = 3 и Y = 2, что является одним из решений для функции Химмельблау.