Популяционные алгоритмы оптимизации: Оптимизация Стаей Серых Волков (Grey Wolf Optimizer - GWO)

Andrey Dik | 2 декабря, 2022

Содержание:

1. Введение
2. Описание алгоритма
3. Тестовые функции
4. Результаты тестов

1. Введение

Алгоритм серых волков — метаэвристический стохастический алгоритм роевого интеллекта, разработанный в 2014 году, идея которого построена на основе модели охоты стаи серых волков. Выделяются четыре типа волков: альфа, бета, дельта и омега. Альфа имеет наибольший «вес» в принятии решений и управлении стаей. Далее идут бета и затем дельта, которые подчиняются альфе и имеют власть над остальными волками. Волк омега всегда подчиняется остальным доминирующим волкам.

В математической модели иерархии волков Альфа-α-волк считается доминирующим волком в стае, и его приказы должны выполняться членами стаи. Бета-β-подчиненные волки, которые помогают альфе в принятии решений и считаются лучшим кандидатом на роль альфы. Волки дельты δ должны подчиняться альфе и бета, но они доминируют над омегой. Волки Омега ω считаются козлами отпущения в стае, являются наименее важными особями в стае, и им разрешается есть только в конце. Альфа считается наиболее благоприятным решением.

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

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

2. Описание алгоритма

Иерархия в стае представлена схематично рисунком 1.  Распределение ролей между членами стаи, доминирующую роль играет alpha.

Рисунок 1. Социальная иерархия в стае волков.

Математическая модель и алгоритм
Социальная иерархия:

• Лучшее решение в виде альфа-волка (α).
• Второе лучшее решение в роли бета-волка (β).
• Третье лучшее решение в качестве дельта-волка (δ).
• Остальные возможные решения как волки Омега (ω).

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

Рисунок 2. Стадии охоты: поиск, окружение, нападение.

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

После этапа "поиск" определяется значение фитнес функции и только после этого можно переходить к этапу "окружение". На этом этапе коэффициент а больше 1, это означает, что альфа бета и дельта отдаляются от своих предыдущих позиций, позволяя тем самым уточнить положение предполагаемой добычи. Когда коэффициент а становится равным 1, тогда начинается этап "нападение" и до окончания итераций коэффициент стремится к 0. Это приводит к приближению волков к добыче, предполагая, что лучшее её положение уже найдено. Хотя, если на этом этапе один из волков обнаружит более лучшее решение, то положение добычи будет обновлено и иерархия волков обновится, но коэффициент все равно будет продолжать стремится к 0. Процесс изменения а представлен нелинейной функцией, этапы схематично показаны на рисунке 2.

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

Рисунок 3. Схема движения омеги по отношению к альфа, бета и дельта.

Псевдокод алгоритма GWO заключается в следующем:

1) Случайным образом инициализировать популяцию серых волков.
2) Рассчитать пригодность каждого члена популяции.
3) Лидеры стаи:
-α = член с лучшим значением пригодности
-β = второй лучший участник (с точки зрения пригодности)
-δ = третий лучший участник (с точки зрения значения пригодности)
Обновить положение всех омега-волков по формулам в зависимости от α, β, δ
4) Рассчитать пригодность каждого члена популяции.
5) повторить шаг 3.

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

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

//——————————————————————————————————————————————————————————————————————————————
struct S_Wolf
{
  double c []; //coordinates
  double p;    //prey
};
//——————————————————————————————————————————————————————————————————————————————

Волчья стая представлена компактным и понятным классом. Здесь объявлены диапазоны и шаг оптимизируемых параметров, лучшее положение добычи, лучшее значение решения и вспомогательные функции.

//——————————————————————————————————————————————————————————————————————————————
class C_AO_GWO //wolfpack
{
  //============================================================================
  public: double rangeMax  []; //maximum search range
  public: double rangeMin  []; //manimum search range
  public: double rangeStep []; //step search
  public: S_Wolf wolves    []; //wolves of the pack
  public: double cB        []; //best prey coordinates
  public: double pB;           //best prey

  public: void InitPack (const int    coordinatesP,   //number of opt. parameters
                         const int    wolvesNumberP,  //wolves number
                         const int    alphaNumberP,   //alpha beta delta number
                         const int    epochCountP);   //epochs number

  public: void TasksForWolves      (int epochNow);
  public: void RevisionAlphaStatus ();


  //============================================================================
  private: void   ReturnToRange (S_Wolf &wolf);
  private: void   SortingWolves ();
  private: double SeInDiSp      (double In, double InMin, double InMax, double Step);
  private: double RNDfromCI     (double Min, double Max);

  private: int    coordinates;     //coordinates number
  private: int    wolvesNumber;    //the number of all wolves
  private: int    alphaNumber;     //Alpha beta delta number of all wolves
  private: int    epochCount;

  private: S_Wolf wolvesT    [];   //temporary, for sorting
  private: int    ind        [];   //array for indexes when sorting
  private: double val        [];   //array for sorting

  private: bool   searching;       //searching flag
};
//——————————————————————————————————————————————————————————————————————————————

Традиционно после объявления класса в описании алгоритма идет инициализация. Здесь мы сбрасываем до минимального значения double приспособленности волков, распределяем размер массивов.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_GWO::InitPack (const int    coordinatesP,   //number of opt. parameters
                         const int    wolvesNumberP,  //wolves number
                         const int    alphaNumberP,   //alpha beta delta number
                         const int    epochCountP)    //epochs number
{
  MathSrand (GetTickCount ());
  searching = false;
  pB        = -DBL_MAX;

  coordinates  = coordinatesP;
  wolvesNumber = wolvesNumberP;
  alphaNumber  = alphaNumberP;
  epochCount   = epochCountP;

  ArrayResize (rangeMax,  coordinates);
  ArrayResize (rangeMin,  coordinates);
  ArrayResize (rangeStep, coordinates);
  ArrayResize (cB,        coordinates);

  ArrayResize (ind, wolvesNumber);
  ArrayResize (val, wolvesNumber);

  ArrayResize (wolves,  wolvesNumber);
  ArrayResize (wolvesT, wolvesNumber);

  for (int i = 0; i < wolvesNumber; i++)
  {
    ArrayResize (wolves  [i].c, coordinates);
    ArrayResize (wolvesT [i].c, coordinates);
    wolves  [i].p = -DBL_MAX;
    wolvesT [i].p = -DBL_MAX;
  }
}
//——————————————————————————————————————————————————————————————————————————————

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

//----------------------------------------------------------------------------
//space has not been explored yet, then send the wolf in a random direction
if (!searching)
{
  for (int w = 0; w < wolvesNumber; w++)
  {
    for (int c = 0; c < coordinates; c++)
    {
      wolves [w].c [c] = RNDfromCI (rangeMin [c], rangeMax [c]);
      wolves [w].c [c] = SeInDiSp  (wolves [w].c [c], rangeMin [c], rangeMax [c], rangeStep [c]);
    }
  }
   
  searching = true;
  return;
}

В канонической версии описания алгоритма присутствуют формулы, оперирующие векторами, но в виде кода они гораздо понятнее. Расчет omega волков идет до расчета alpha, beta и delta волков, потому что нужно использовать предыдущие значения лидеров.

Основной компонент, обеспечивающий три этапа охоты (поиск, окружение, нападение) - коэффициент a. Он представляет собой нелинейную зависимость от текущей итерации и общего количества итераций и стремится к 0.
Следующие компоненты формулы Ai и Сi:

• Ai = 2.0 * a * r1 - a;
• Ci = 2.0 * r2;

где r1 и r2 случайные числа в диапазоне [0.0;1.0].
В выражении координаты волков корректируются по усреднённым значениям волков лидеров. Поскольку в алгоритме можно задавать любое количество лидеров, то суммирование координат производим в цикле. После этого полученную сумму делим на количество лидеров. Данную операцию проводим для каждой координаты по отдельности, генерируя каждый раз новое значение r1 и r2. Новое положение омега волков, как видим, корректируется по положениям лидеров с учетом собственного текущего положения.

//----------------------------------------------------------------------------
double a  = sqrt (2.0 * (1.0 - (epochNow / epochCount)));
double r1 = 0.0;
double r2 = 0.0;

double Ai = 0.0;
double Ci = 0.0;
double Xn = 0.0;

double min = 0.0;
double max = 1.0;

//omega-----------------------------------------------------------------------
for (int w = alphaNumber; w < wolvesNumber; w++)
{
  Xn = 0.0;

  for (int c = 0; c < coordinates; c++)
  {
    for (int abd = 0; abd < alphaNumber; abd++)
    {
      r1 = RNDfromCI (min, max);
      r2 = RNDfromCI (min, max);
      Ai = 2.0 * a * r1 - a;
      Ci = 2.0 * r2;
      Xn += wolves [abd].c [c] - Ai * (Ci * wolves [abd].c [c] - wolves [w].c [c]);
    }

    wolves [w].c [c] = Xn /= (double)alphaNumber;
  }

  ReturnToRange (wolves [w]);
}

А вот и расчет лидеров. Для них так же рассчитываются коэффициенты a, Ai и Ci для каждой координаты. Отличие заключается лишь в том, что положение лидеров изменяется в отношении к лучшим координатам добычи на текущий момент и собственного положения. Лидеры кружат вокруг добычи, то приближаясь, то отдаляясь, и управляют второстепенными волками в нападении.

//alpha, beta, delta----------------------------------------------------------
for (int w = 0; w < alphaNumber; w++)
{
  for (int c = 0; c < coordinates; c++)
  {
    r1 = RNDfromCI (min, max);
    r2 = RNDfromCI (min, max);

    Ai = 2.0 * a * r1 - a;
    Ci = 2.0 * r2;

    wolves [w].c [c] = cB [c] - Ai * (Ci * cB [c] - wolves [w].c [c]);
  }

  ReturnToRange (wolves [w]);
}

Это второй публичный метод, вызываемый на каждой итерации. Тут обновляется статус лидеров в стае, а по сути - производится сортировка волков по значению приспособленности. Если обнаружены координаты добычи лучше, чем сохранённые по всей стае, то обновляем значения.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_GWO::RevisionAlphaStatus ()
{
  SortingWolves ();

  if (wolves [0].p > pB)
  {
    pB = wolves [0].p;
    ArrayCopy (cB, wolves [0].c, 0, 0, WHOLE_ARRAY);
  }
}
//——————————————————————————————————————————————————————————————————————————————

3. Тестовые функции

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

1. Глобальный экстремум не должен находиться на границах диапазона. В случае, если в алгоритме отсутствует проверка выхода за диапазон, то возможны ситуации, когда алгоритм будет показывать превосходные результаты. Это связано с тем, что благодаря внутреннему дефекту, значения будут находиться на границах.
2. Глобальный экстремум не должен находиться в центре координат диапазона. В этом случае будет учтена ситуация, когда алгоритм генерирует значения усредненные в диапазоне.
3. Глобальный минимум должен находиться в центре координат. Это необходимо, чтобы заведомо исключить ситуации, описанные в п.2.
4. Подсчет результатов тестовой функции должен учитывать момент, когда случайно сгенерированные числа по всей области определения функции (когда функция многопеременная) дадут усредненный результат около 50% от максимума, хотя по факту данные результаты получены случайно.
   

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

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

Рисунок 4. Тестовая функция Skin.

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

Рисунок 5. Тестовая функция Forest.

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

Рисунок 6. Тестовая функция Megacity.

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

4. Результаты тестов

В связи с внесенными изменениями в тестовые функции обновился и тестовый стенд. В правой части экрана стенда по-прежнему можно наблюдать графики сходимости алгоритмов оптимизации. Зеленая линия - результаты сходимости на функции с двумя переменными. Голубая - функции с 40 переменными. Красная - функции с 1000 переменными. Черный кружок большего размера указывает на положение глобального максимума функции. Черный кружок меньшего размера указывает на положение текущего значения решения алгоритма оптимизации. Перекрестие белых линий указывает на геометрический центр тестовых функций и соответствует глобальному минимуму, это введено для лучшего визуального восприятия поведения тестируемых алгоритмов. Белыми точками обозначены усредненные промежуточные решения. Цветными мелкими точками обозначены пары координат соответствующей размерности, цвет указывает на порядковое положение размерности тестовой функции.

Обновлённые результаты тестирования алгоритмов оптимизации, рассмотренных в предыдущих статьях, на новом стенде, можно увидеть в обновлённой таблице. Из таблицы убрана строка о скорости сходимости — это визуально можно определить на анимации стенда, а таблицу неоправданно загружает для восприятия. Добавлена колонка с описанием алгоритма.

Результаты испытаний ACOm (Ant Colony Optimization):

Результаты испытаний ABCm (Artificial Bee Colony):

Результаты испытаний ABC (Artificial Bee Colony):

Результаты испытаний PSO (Particle Swarm Optimization)

Результаты испытаний RND (Random):

  GWO на тестовой функции Skin.

GWO на тестовой функции Forest.

  GWO на тестовой функции Megacity.

Результаты испытаний GWO.

Алгоритм оптимизации "Серых волков", GWO — это один из недавних биоинспирированных алгоритмов оптимизации, основанный на имитации загонной охоты стаи серых волков. В среднем алгоритм очень хорошо зарекомендовал себя на различных типах функций, как по точности поиска экстремума, так и по скорости сходимости. В некоторых тестах он является лидером в рейтинговой таблице. Основные показатели эффективности алгоритма оптимизации «Серых волков» лучше, чем у алгоритма оптимизации роем частиц, который в классе биоинсперированых алгоритмов оптимизации является "классическим".

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

Причем невысокая точность найденного экстремума проявилась на всех размерностях функции Forest, и результаты худшие среди всех участников таблицы. Хотя на гладкой функции Skin алгоритм проявил себя превосходно, являясь лидером, особенно с функцией Skin большой размерности. Так же GWO является третьим в таблице, кто смог добиться 100% попадания в глобальный максимум на функции Megacity.

AO	Description	Skin			Forest			Megacity (discrete)			Final result
		2 params (1 F)	40 params (20 F)	1000 params (500 F)	2 params (1 F)	40 params (20 F)	1000 params (500 F)	2 params (1 F)	40 params (20 F)	1000 params (500 F)
ACOm	ant colony optimization	0,98229	0,79108	0,12602	1,00000	0,62077	0,11521	0,38333	0,44000	0,02377	0,49805222
ABCm	artificial bee colony M	1,00000	0,63922	0,08076	0,99908	0,20112	0,03785	1,00000	0,16333	0,02823	0,46106556
ABC	artificial bee colony	0,99339	0,73381	0,11118	0,99934	0,21437	0,04215	0,85000	0,16833	0,03130	0,46043000
GWO	grey wolf optimizer	0,99900	0,48033	0,18924	0,83844	0,08755	0,02555	1,00000	0,10000	0,02187	0,41577556
PSO	particle swarm optimisation	0,99627	0,38080	0,05089	0,93772	0,14540	0,04856	1,00000	0,09333	0,02233	0,40836667
RND	random	0,99932	0,44276	0,06827	0,83126	0,11524	0,03048	0,83333	0,09000	0,02403	0,38163222

Выводы: