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

import numpy as np
np.set_printoptions(precision=2, linewidth=140)
from numpy import random
import math
from dataclasses import dataclass
import matplotlib.pyplot as plt

Рассматривается статья метаэвристическое решения задачи о булевом рюкзаке, разновидностью гибридной комбинации пары алгоритмов «роевого интеллекта»:

• «Алгоритм кукушки»

  • Кукушка закладывает по одному яйцу за раз, и сбрасывает его в случайно выбранном гнезде;

  • «Лучшие гнезда с высоким качеством яиц» (решения) переносятся на следующие поколения;

  • Количество доступных хозяйских гнезд фиксировано, и хозяин может обнаружить инородное яйцо с вероятностью Pa ∈ [0, 1]. В этом случае птица-хозяин может либо выбросить яйцо из гнезда, либо покинуть гнездо, чтобы построить совершенно новое гнездо в новом месте.

• «Алгоритм перетасованных лягушачьих прыжков»

  • Существует популяция лягушек - P. Она состоит из непересекающихся подгрупп лягушек называемых "мемеплекс". Всего мемеплексов M штук. В каждом мемеплексе N лягушек.

  • У каждой лягушки в популяции есть своя позиция (задаваемая вектором x). Процедура эволюции (один шаг алгоритма) состоит из трех прыжков каждой лягушки (с каждым прыжком лягушка меняет свое положение, то есть меняется вектор x). Эти три прыжка рассмотрены ниже.

frog-alg.png

Нужно применить «Алгоритм перетасованных лягушачьих прыжков» для задачи о рюкзаке.

Задача о рюкзаке

• Пусть имеется набор предметов, каждый из которых имеет два параметра — масса и ценность. Также имеется рюкзак определённой грузоподъёмности. Задача заключается в том, чтобы собрать рюкзак с максимальной ценностью предметов внутри, соблюдая при этом ограничение рюкзака на суммарную массу.

• Математическое описание задачи приведено на картинке

knapsack.png

Эксперимент

В алгоритме перетасованных лягушачьих прыжков (SFLA) используется положение лягушки при помощи вектора x. Применяя алгоритм к задаче о рюкзаке, вектор x будет задаваться: vector_x.png, где d — кол-во ценных вещей, то есть размерность задачи о рюкзаке.

Чтобы преобразовать этот вектор в бинарный вектор b = vector_b.png (который определяет кладется ли вещь в рюкзак) используется сигмоидная функция:

sigmoid.png

def sigmoid(x):
    return 1 / (1 + math.e ** -x)

x = np.linspace(start=-3, stop=3, num=100)
plt.plot(x, sigmoid(x))
plt.show()

Для всей популяции инициализируется матрицы X и B, состоящие из векторов для каждой лягушки <$X_i$, $B_i$>.

Подготовка данных

Некоррелированный случай

uncorrelated.png

d = 1000 # размерность задачи о рюкзаке
w = random.randint(10, 100, d) # вектор весов
p = random.randint(10, 100, d) # вектор ценности

c = 0.75 * w.sum() # ограничение рюкзака по весу

Популяция лягушек

@dataclass
class Population:
    L: float # Минимально разрешенное изменение позиции
    U: float # Максимально разрешенное изменение позиции
    M: int # число мемеплексов («супергенов»)
    N: int # число лягушек в каждом из мемеплексов → Каждый «мемеплекс лягушек» как бы использует свою стратегию поиска оптимальных решений.

#     def __repr__(self):
#        'Если хотим хитро печатать объект, можно использовать такого типа функцию, но может хватит и просто датаклассов.'    
#        return "Популяция размера {d}. {L}←→{U} {w} {p} {c} {M}×{N} ".format(**vars(self))

    def __post_init__(self):
        self.P = self.M * self.N # общее число лягушек в популяции
        self.X = self.initialize_population((self.P, d))

    def initialize_population(self, size = None):
        if size is None:
            size = d
        X = random.uniform(self.L, self.U, size)
        return X

    def reset_population(self):
        self.X = self.initialize_population((self.P, d))

class Population(Population):
    def B(self, X = None):
        '''
          Бинаризация вектора или матрицы X по сигмоиде
          Если вектор или матрица X не подается на вход, берется текущее значение положений лягушек для данной популяции
        '''
        if X is None:
            X = self.X
        return sigmoid(X) >= 0.5

    def f(self, X = None, B = None):
        '''
          Значение функции ценности f
          Если вектор или матрица B не подается на вход, B считается по сигмоиде от матрицы или вектора X
        '''
        if B is None:
            B = self.B(X)
        return (p * B).sum(len(B.shape) - 1)

    def isLegal(self, x):
        '''
        Возвращает булево значение вместимости вектора x в рюкзак, а также вес рюкзака
        '''
        weight = (w * x).sum(len(x.shape) - 1)
        return weight <= c, weight

    def B_best(self, X = None):
        '''
        Возвращает наилучшее текущее решение по всей популяции
        '''
        B = self.B(X)
        return B[np.argmax(self.f(X))]

    def B_k_best(self, X = None):
        '''
        Возвращает наилучшее текущее решение для каждого мемеплекса
        '''
        B = self.B(X)
        return [B[np.argmax(self.f(X)[k::self.M])] for k in range(self.M)]

Улучшенный алгоритм SFLA

Ранее были представлены 3 прыжка лягушки в алгоритме SFLA:

3steps_1.png | 3steps_2.png

Теперь мы улучшаем его:

improvedsfla.png

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

Также первые два прыжка лягушки (4)-(5) заменяются на один прыжок (8)

Псевдокод улучшенного алгоритма SFLA

sfla.png

Лягушачьи параметры

M — число мемплексов («супергенов»), в каждом из которых по N «лягушек»

Каждый «мемплекс лягушек» как бы использует свою стратегию поиска оптимальных решений

M = 4
N = 10

L = -3 # Минимально разрешенное изменение позиции
U =  3 # Максимально разрешенное изменение позиции

$p_m$ — это вероятность случайной инициализации нового положения лягушки.

p_m = 0.15

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

def p1_random(i, maxval):
    '''
    На вход подается текущая лягушка i и диапазон случайного выбора другой лягушки
    На выходе индекс лягушки не равный i
    '''
    p1 = random.randint(maxval)
    if p1 != i:
        return p1
    else:
        return p1_random(i, maxval)

Инициализируем популяцию лягушек

population = Population(L, U, M, N)
population

np.set_printoptions(precision=1)
population.X # случайные положения

population.B()

plt.bar(np.arange(population.M * population.N), population.f())
plt.title("Значение функции ценности для каждой лягушки")
plt.show()

np.set_printoptions(precision=2)

def frog_leaping_algorithm(population, p_m, d):
    for i in range(population.M * population.N):
        k = i % population.M # номер мемеплекса
        p1 = p1_random(i, population.M * population.N) # случайный выбор другой лягушки
        r1 = np.array([1 if (random >= 0.5) else -1 for random in random.rand(d)]) # прибавляем или вычитаем
        r2 = random.rand(d)
        Y = population.B_best() + r1 * r2 * (population.B_k_best()[k] - population.X[p1])

        if (population.f(Y) > population.f(population.X[i])):
            population.X[i] = Y
        else:
            if (random.rand() <= p_m):
                population.X[i] = population.initialize_population()

old_population_f = population.f().copy()

frog_leaping_algorithm(population, p_m, d) # одна итерация алгоритма

bw = 0.3
plt.bar(np.arange(population.M * population.N), old_population_f, bw, color = 'r')
plt.bar(np.arange(population.M * population.N) + bw, population.f(), bw, color = 'g')
plt.title("Значение функции ценности для старого и нового положения лягушки (зеленые - новые положения)")
plt.show()

Стоит отметить, что получившиеся вектора b (принадлежность предмета рюкзаку) может выходить за лимит по весу рюкзака:

print("Ограничение по весу = ", c, " (выделено красным)")
plt.plot(np.arange(population.M * population.N), np.ones(population.M * population.N) * c, color = 'r')
plt.bar(np.arange(population.M * population.N), population.isLegal(population.B())[1])
plt.title("Значение весов для каждой лягушки")
plt.show()

Как мы видим, в некоторых случаях так и есть. Чтобы этого избежать, необходимо использовать восстанавливающий GTM алгоритм.

repair1.png | repair2.png

Как описано, данный восстанавливающий алгоритм содержит в себе 2 фазы. Первая фаза восстановления проверяет каждое решение, причем идет в порядке уменьшения полезности предмета (уменьшение p/w), и подтверждает принадлежность предмета если ограничение по весу не нарушается. Вторая фаза, названная фазой оптимизации, меняет оставшиеся нулевые значения вектора решения (то есть те предметы, которые пока не в рюкзаке) на единицу, если при этом ограничение по весу не нарушается.

Для начала отсортируем индексы в порядке убывания p/w

s = np.argsort((p / w))[::-1]

def GTM_repair(population, x):
    '''
    На вход подается вектор положения лягушек
    На выходе получаем восстановленный вектор положения лягушек
    '''
    b = population.B(x)
    weight = 0
    for index in s: # Фаза восстановления
        if (b[index]):
            if (weight:= weight + w[index]) > c:
                weight -= w[index]
                x[index] = -1.0
                b[index] = 0
    rest = c - weight
    for index in s: # Фаза оптимизации
        if (w[index] <= rest and not b[index]):
            x[index] = 0.001
            b[index] = 1
            rest -= w[index]
            weight += w[index]
    return x

print("Ограничение по весу = ", c)
print("Старый вес десятой лягушки = ", population.isLegal(population.B()[10])[1])
print("Вес после GTM восстанавливающего алгоритма = ", population.isLegal(population.B(GTM_repair(population, population.X[10])))[1])

CSISFLA algorithm

Полеты Леви

В природе животные ищут пищу случайным образом.

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

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

      Xc(t+1)  = Xc(t) + α ⊕ Levy(λ)

где,

• t – поколение или номер итерации,

• a – размер шага, который зависит от масштаба задачи (обычно принимают равной единице),

• ⊕ — произведение Адамара.

• Levy — берется из распределения Леви, Levy (𝜆) ∼ 𝑢 = 𝑡^{−𝜆}

| d.png | levy-principles.png |

Полет Леви используется в алгоритме CSISFLA для прыжка лягушки. Реализуем это распределение.

# Параметр распределения:

lamda = 1.5

def Levy():
    # Генерация прыжка используя распределение Леви
    sigma1 = np.power((math.gamma(1 + lamda) * np.sin((np.pi * lamda) / 2)) \
                      / math.gamma((1 + lamda) / 2) * np.power(2, (lamda - 1) / 2), 1 / lamda)
    sigma2 = 1
    u = np.random.normal(0, sigma1, size = d)
    v = np.random.normal(0, sigma2, size = d)
    step = u / np.power(np.fabs(v), 1 / lamda)

    return step # вектор изменения положения лягушки

Псевдокод улучшенного алгоритма SFLA

csisfla.png

@dataclass
class CSISFLA:
    population: Population # популяция
    p_m: float # вероятность случайной инициализации нового положения лягушки.
    alpha: float # размер шага для полета Леви
    s: np.array # 1 шаг: Сортировка индексов
    amount_of_iter: int # Количество итераций главного алгоритма

    def __post_init__(self):
        self.population.reset_population() # 2 шаг: Инициализация популяции
        self.G = 0 # итерации алгоритма
        self.best_B = self.population.B_best()
        self.best_B_value = self.population.f(B = self.best_B)
        self.history_best_value = [self.best_B_value]
        self.history_weight = [self.population.isLegal(self.best_B)[1]]

    def main_algorithm(self):
        while self.stopping_criteria(): # 3 шаг: Главный алгоритм с использованием полетов Леви
            self.shuffled_frog_leaping_algorithm()
            B_best_now = self.population.B_best()
            if (self.population.f(B = B_best_now) > self.best_B_value):
                self.best_B_value = self.population.f(B = B_best_now)
                self.best_B = B_best_now
            self.history_best_value.append(self.best_B_value)
            self.history_weight.append(self.population.isLegal(self.best_B)[1])
            self.G += 1
            self.shuffle_memeplexes() # 4 шаг: перемешать ляшушек в каждом мемеплексе

    def shuffle_memeplexes(self):
        for k in range(self.population.M):
            random.shuffle(population.X[k::self.population.M])

    def shuffled_frog_leaping_algorithm(self):
        '''
        Улучшенный SFLA с добавлением полетов Леви
        '''
        for i in range(self.population.P):
            k = i % self.population.M
            p1 = p1_random(i, self.population.P)
            self.population.X[i] += self.alpha * Levy()
            r1 = np.array([1 if (random >= 0.5) else -1 for random in random.rand(d)])
            r2 = random.rand(d)
            Y = self.population.B_best() + r1 * r2 * (self.population.B_k_best()[k] - self.population.X[p1])

            if (self.population.f(Y) > self.population.f(self.population.X[i])):
                self.population.X[i] = Y
            else:
                if (random.rand() <= self.p_m):
                    self.population.X[i] = self.population.initialize_population()
            self.population.X[i] = GTM_repair(population, self.population.X[i])
            #print("profits ", self.population.f(self.population.X[i]))

    def stopping_criteria(self):
        '''
        Устанавливаем критерий остановки алгоритма
        '''
        return self.G <= self.amount_of_iter

amount_of_iter = 10

alpha = 0.01 # размер шага для полета Леви

csisfla = CSISFLA(population, p_m, alpha, s, amount_of_iter)

# Старт алгоритма
csisfla.main_algorithm()

plt.plot(np.arange(len(csisfla.history_best_value)), csisfla.history_best_value, color = 'g')
plt.title("Изменение ценности функции по итерациям")
plt.show()

csisfla.history_best_value

print("Ограничение по весу = ", c, " (выделено красным)")
plt.plot(np.arange(len(csisfla.history_weight)), np.ones(len(csisfla.history_weight)) * c, color = 'r')
plt.plot(np.arange(len(csisfla.history_weight)), csisfla.history_weight, color = 'blue')
plt.title("Изменение веса по итерациям")
plt.show()

csisfla.history_weight

print('Итоговая ценность =', csisfla.best_B_value)
print('Итоговый вес:', population.isLegal(csisfla.best_B)[1])
print('Ограничение по весу: ', c)

Решение с помощью OR-tools

from ortools.algorithms import pywrapknapsack_solver

solver = pywrapknapsack_solver.KnapsackSolver(
    pywrapknapsack_solver.KnapsackSolver.
    KNAPSACK_MULTIDIMENSION_BRANCH_AND_BOUND_SOLVER, 'KnapsackExample')

values = list(p)
weights = [list(w)]
capacities = [c]
solver.Init(values, weights, capacities)
computed_value = solver.Solve()

packed_items = []
packed_weights = []
total_weight = 0
print('Итоговая ценность =', computed_value)
for i in range(len(values)):
    if solver.BestSolutionContains(i):
        packed_items.append(i)
        packed_weights.append(weights[0][i])
        total_weight += weights[0][i]
print('Итоговый вес:', total_weight)

print('Итоговая ценность csisfla =', csisfla.best_B_value)
print('Итоговая ценность OR-tools =', computed_value)
print('Итоговый вес csisfla:', population.isLegal(csisfla.best_B)[1])
print('Итоговый вес OR-tools:', total_weight)
print('Ограничение по весу: ', c)

График сходимости алгоритма CSISLFA к оптимуму

plt.plot(np.arange(len(csisfla.history_best_value[:4])), csisfla.history_best_value[:4], color = 'g', label = 'csislfa')
plt.plot(np.arange(len(csisfla.history_best_value[:4])), computed_value * np.ones(len(csisfla.history_best_value[:4])), color = 'r', label = 'OR-tools')
plt.legend()
plt.title("Изменение ценности функции по итерациям")
plt.xlabel("итерации", fontsize = 12)
plt.ylabel("ценность", fontsize = 12)
plt.show()

Вывод

Как мы можем видеть, алгоритм быстро сходится, начиная с 3-4 итерации уже достигается максимальный результат функции f.

Еще как мне кажется, сильно помогает оптимизировать решения алгоритм GTM, который восстанавливает положения лягушек и оптимизирует их.

По сравнению с решением из OR-tools алгоритм csisfla проигрывает совсем чуть чуть