智能算法系列之蚁群算法(详细补充版)

  本博客封面由ChatGPT + DALL·E 2共同创作而成。

前言

  本篇是智能算法(Python复现)专栏的第六篇文章，主要是对蚁群算法的补充。
  在上篇博客中对蚁群算法的原理进行了简单的介绍，并用python实现蚁群算法在求函数极值方面的应用，但仔细研读一番发现重点没有突出来，而且与原始TSP问题中的蚁群算法又较大区别，感觉只涉及到蚁群算法的外壳，并没有触及灵魂。
  因此，在本篇中对蚁群算法在求函数极值方面的应用进行一个详细补充，并尽可能的触及蚁群算法的灵魂。

  论文参考：应用蚁群算法求解函数所有极值(知网)

1. 改进方法

  设一维函数$y=f(x)$的定义域为$[a,b]$。
  (1) 将区间$[a,b]$划分成若干个长度相等的子区间，这些子区间记为 { I 1 , I 2 , . . . , I n } \{I_1,I_2,...,I_n\} {I1​,I2​,...,In​}，取区间$I_i$的中点记为$x_i$；
  (2) 设与区间$I_i$相邻的区间为$I_{i+1}$，假设$I_i$和$I_{i+1}$之间有一条虚拟的边$e(I_{i+1},I_{i+1})$，该边的权重(虚拟距离)和$f(x_i)-f(x_{i+1})$的大小有关，$f(x_i)-f(x_{i+1})$越大，蚂蚁从区间$I_i$转移到区间$I_{i+1}$的可能性就越大；
  (3) 蚂蚁从区间$I_i$转移到区间$I_{i+1}$后，要留下信息素，所留信息素的大小是一个与$f(x_i)-f(x_{i+1})$有关的量，区间$I_{i+1}$的信息素越多，就越吸引相邻区间蚂蚁向其转移；
  (4) 蚂蚁经过许多次转移之后，有的区间含有许多蚂蚁，而有的区间则不含蚂蚁。那些含有蚂蚁的区间正是包含极值点的区间，不含蚂蚁的区间，不大可能包含极值点；
  (5) 取出包含蚂蚁的区间，并将它们重新细化，重复上述搜索过程，直到细化后的区间足够小。
  最后，蚂蚁都停留在了极值点附近，蚂蚁所在的区间的中点位置正是极值点的位置。

2. 蚁群初始化

  首先将问题的定义域$[a,b]$进行$n$等分，等分后的各区间长度为$$\sigma =\frac{b-a}{n}$$各区间记为 { I 1 , I 2 , . . . , I n } \{I_1,I_2,...,I_n\} {I1​,I2​,...,In​}，)，其中$I_i(i=1,2,..,n)$表示第$i$个区间，则$$I_i=[a+(i-1)\sigma ,a+i\sigma ]$$，表示区间$I_i$的左右两个端点，区间$I_i$的中点位置记为$x_i$，则$$x_i=a+(i-\frac{1}{2})\sigma$$  假设蚁群中共有$m$只蚂蚁，记为 { a 1 , a 2 , . . . , a n } \{a_1,a_2,...,a_n\} {a1​,a2​,...,an​}一般情况下，$m\ge n$。在初始时，为每只蚂蚁随机分配一个区间，蚂蚁的位置即为区间中点的位置。${\tau }_i(t)$表示$t$时刻子区间$i$上的信息素，初始时各个子区间的信息素浓度相同，默认为1.0；$\Delta {\tau }_i(t)$表示$t$时刻子区间$i$上的信息素增量，初始时各个子区间的信息素增量也相同，默认为0。

3. 蚂蚁移动策略

  假设$neighbor(I_i)$表示与区间$I_i$邻近的区间集合，则对于一维函数 n e i g h b o r ( I i ) = { { I i + 1 } , i = 1 { I i − 1 , I i + 1 } , i = 2 , 3 , . . . , n − 1 { I i − 1 } , i = n neighbor(I_i) = \begin{cases} \{I_{i+1}\}, & i=1\\ \{I_{i-1}, I_{i+1}\}, & i=2,3,...,n-1\\ \{I_{i-1}\}, & i=n \end{cases} neighbor(Ii​)=⎩ ⎨ ⎧​{Ii+1​},{Ii−1​,Ii+1​},{Ii−1​},​i=1i=2,3,...,n−1i=n​  位于区间$I_i$的蚂蚁向其邻近区间$I_j$进行转移，假设$I_i$和$I_j$之间有一条虚拟的边$e(I_i,I_j)$，该边的权重为$|f(x_i)-f(x_j)|$，则启发函数$${\eta }_{ij}=|f(x_i)-f(x_j)|$$  设蚂蚁$a_k$当前处于区间$I_i$，若$f(x_i)-f(x_j)>0$，表示当前区间$I_i$上的值较大，蚂蚁$a_k$就可以转移到其邻近区间$I_j$；否则，蚂蚁不向其转移。设用$allowe{d}_k$表示蚂蚁$a_k$下一步可以转移的子区间的集合。
  用$p_{ij}^k(t)$表示第$t$次循环蚂蚁$a_k$从子区间$I_i$转移到子区间$I_j$的概率，模仿基本蚁群算法，定义蚂蚁$a_k$的转移概率如下：$$p_{ij}^k(t)=\{\begin{array}{ll}\frac{{\tau }_j^{\alpha }(t){\eta }_{ij}^{\beta }}{{\sum }_{h\in allowe{d}_k}{\tau }_j^{\alpha }(t){\eta }_{ih}^{\beta }},& j\in allowe{d}_k\\ & \\ 0,& j\notin allowe{d}_k\end{array}$$  其中，$\alpha$为信息启发式因子(信息素重要程度因子)，$\beta$为期望启发式因子(启发函数重要程度因子)，${\tau }_j$为子区间$I_i$的信息素浓度。

4. 信息素更新

  如果蚂蚁$a_k$依概率从子区间$I_i$转移到了子区间$I_j$，那么蚂蚁$a_k$就会在区间留下信息素，所留信息素的量用$\Delta {\tau }_j^k(t)$表示，则$$\Delta {\tau }_j^k(t)=C(f(x_i)-f(x_j))$$  其中，$C$表示一个常量，$f(x_i)-f(x_j)$越大，表示蚂蚁$a_k$留下的信息素就越多，就越吸引其他蚂蚁向区间$I_j$转移。
  所有转移到区间$I_j$的蚂蚁都要留下相应的信息素，假设在第$t$次循环，有$b$只蚂蚁转移到了区间$I_j$，则这$b$只蚂蚁在区间$I_j$上所留信息素总和为$$\Delta {\tau }_j(t)=\sum _{c=1}^b\Delta {\tau }_j^c(t)$$  所有蚂蚁完成一次邻近区间转移后，要对信息素进行更新处
理。区间$I_j$的信息素浓度${\tau }_j$更新为$${\tau }_j(t+1)=(1-1\rho ){\tau }_j(t)+\Delta {\tau }_j(t)$$  其中，$\rho$表示信息素挥发系数。

5. 缩小蚁群搜索空间

  函数值较小的区间含有的信息素会较多，更容易吸引蚂蚁向其转移。经过一些循环，当所有蚂蚁都停止转移的时候，蚁群分布就会出现这样的特点：所有蚂蚁都分布在极值点较小的区间，而其他区间则没有蚂蚁，即含有蚂蚁的区间包含极小值点。
  取出这些包含蚂蚁的区间，并将这些区间重新细化，在下一次循环时，让蚁群在这些区间重新搜索极小值点。如此循环下去，蚁群的搜索范围就会越来越小，当细化后的区间足够小的时候，所有蚂蚁就会停留在极值点附近，蚂蚁所停留的区间的中点位置就是极值点的位置。
  论文中设了一个阈值，毕竟不可能无限细分下去，当区间步长小于阈值时，即停止细分区间，搜索结束。在本文的实现中并没有采取论文中的操作，而是在初始时直接采取一个较小的划分区间步长，这样搜索出的最优解很接近于理论最优值。

6. 代码实现

# -*- coding:utf-8 -*-
# Author:   xiayouran
# Email:    youran.xia@foxmail.com
# Datetime: 2023/5/6 14:59
# Filename: ant_colony_optimization.py
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from base_algorithm import BaseAlgorithm


__all__ = ['AntColonyOptimization']


class Ant:
    def __init__(self):
        self.id = None          # 蚂蚁所在区间的id
        self.eta = None         # 启发式信息, 存放本区间与其左右两个区间的函数差值
        self.tran_prob = None   # 状态转移概率


class Interval:
    def __init__(self):
        self.position = None    # 该区间上的中点
        self.tau = 1.           # 信息素浓度值
        self.delta_tau = 0.     # 信息素浓度增量
        self.count = 0          # 此区间上蚂蚁的数量


class AntColonyOptimization(BaseAlgorithm):
    def __init__(self, population_size=100, max_iter=200, alpha=1.5, beta=0.8, rho=0.3, epsilon=1e-4,
                 q=1.0, step=0.01, x_range=(0, 5), seed=10086):
        super(AntColonyOptimization, self).__init__()
        self.__population_size = population_size  # 蚂蚁种群大小
        self.__max_iter = max_iter  # 最大迭代次数
        self.__alpha = alpha        # 信息素重要程度因子
        self.__beta = beta          # 启发函数重要程度因子
        self.__rho = rho            # 信息素蒸发系数
        self.__epsilon = epsilon    # 停止搜索门限
        self.__q = q                # 信息素释放增量系数, 常量C
        self.__step = step          # 子区间间隔
        self.__x_range = x_range    # 变量x的定义域
        self.__population = []      # 蚁群
        self.__interval_set = []    # 区间集合(对应TSP问题中的city)
        self.__tabu = []            # 禁忌表
        self.__seed = seed
        self.optimal_solution = None

        np.random.seed(seed)

    def init_interval_set(self):
        for left_point in np.arange(*self.__x_range, self.__step):
            interval = Interval()
            interval.position = left_point + self.__step / 2  # 区间中点
            self.__interval_set.append(interval)

        if len(self.__interval_set) > self.__population_size:
            tmp_size = self.__population_size
            self.__population_size = int(np.ceil(len(self.__interval_set) / self.__population_size) * self.__population_size)
            print("Suggest a larger value for population_size, the value for population_size has been "
                  "changed from {} to {}".format(tmp_size, self.__population_size))

    def init_tabu(self):
        self.__tabu = np.zeros(shape=(len(self.__interval_set), 2))   # 初始禁忌表, 当前点与左右两坐标
        # TSP中的禁忌表的shape是(m, n), 其中m为蚂蚁数量, n为城市数量, 确保每只蚂蚁仅能访问城市一次

    def update_eta(self, ant):
        index = ant.id
        interval = self.__interval_set[index]
        ant.eta = []

        if index == 0:
            ant.eta.append(0)
            ant.eta.append(self.problem_function(interval.position) - self.problem_function(
                interval.position + self.__step))
        elif index == len(self.__interval_set) - 1:
            ant.eta.append(self.problem_function(interval.position) - self.problem_function(
                interval.position - self.__step))
            ant.eta.append(0)
        else:
            ant.eta.append(self.problem_function(interval.position) - self.problem_function(
                interval.position - self.__step))  # 当前区间(中点)与左邻居区间(中点)的差值
            ant.eta.append(self.problem_function(interval.position) - self.problem_function(
                interval.position + self.__step))  # 当前区间(中点)与右邻居区间(中点)的差值

    def update_tabu(self, ant):
        index = ant.id
        if ant.eta[0] > 0:
            self.__tabu[index, 0] = 1  # 表示左子区间值较小, 可以跳转
        if ant.eta[1] > 0:
            self.__tabu[index, 1] = 1  # 表示右子区间值较小, 可以跳转

    def init_population(self):
        # 初始化区间集合
        self.init_interval_set()

        # 初始化禁忌表
        self.init_tabu()

        for i in range(self.__population_size):
            index = np.random.choice(range(len(self.__interval_set)))   # 随机选择一个区间
            # index = i
            interval = self.__interval_set[index]
            interval.count += 1

            ant = Ant()
            ant.id = index

            # 更新eta
            self.update_eta(ant)

            # 更新禁忌表tabu
            self.update_tabu(ant)

            # 更新蚂蚁的状态转移概率
            ant.tran_prob = interval.tau * ant.eta[0] / (interval.tau * ant.eta[0] + interval.tau * ant.eta[1])  # 蚂蚁向左区间跳转概率

            self.__population.append(ant)

    def skip_left(self, ant):
        index = ant.id

        interval = self.__interval_set[index]
        interval.count -= 1  # 此区间蚂蚁数-1
        left_interval = self.__interval_set[index - 1]
        left_interval.count += 1  # 左区间蚂蚁数+1
        ant.id -= 1  # 蚂蚁跳转至左区间
        left_interval.delta_tau = left_interval.delta_tau + self.__q * ant.eta[0]

    def skip_right(self, ant):
        index = ant.id

        interval = self.__interval_set[index]
        interval.count -= 1  # 此区间蚂蚁数-1
        right_interval = self.__interval_set[index + 1]
        right_interval.count += 1  # 右区间蚂蚁数+1
        ant.id += 1  # 蚂蚁跳转至右区间
        right_interval.delta_tau = right_interval.delta_tau + self.__q * ant.eta[1]

    def search_local_optimal_solution(self):
        flag = np.ones(self.__population_size)
        while np.sum(flag):
            for i, ant in enumerate(self.__population):
                index = ant.id
                if self.__tabu[index, 0] and not self.__tabu[index, 1]:
                    # 蚂蚁可以向左区间跳转
                    self.skip_left(ant)
                elif not self.__tabu[index, 0] and self.__tabu[index, 1]:
                    # 蚂蚁可以向右区间跳转
                    self.skip_right(ant)
                elif self.__tabu[index, 0] and self.__tabu[index, 1]:
                    # 两个区间都可以跳转, 计算一下蚂蚁的状态转移概率
                    if ant.tran_prob > np.random.rand():
                        # 蚂蚁向左区间跳转
                        self.skip_left(ant)
                    else:
                        # 蚂蚁向右区间跳转
                        self.skip_right(ant)
                else:
                    flag[i] = 0     # 表示此蚂蚁不再进行跳转了

                self.update_eta(ant)    # 更新eta
                self.update_tabu(ant)   # 更新禁忌表

            for interval in self.__interval_set:
                # 更新区间上的信息素
                interval.tau = (1 - self.__rho) * interval.tau + interval.delta_tau

    def print_local_optimal_solution(self):
        print('local optimal solution:')
        local_optimal_solution = {}
        best_point = np.inf
        for ant in self.__population:
            index = ant.id
            if not local_optimal_solution.get(index, ''):
                local_optimal_solution[index] = (self.__interval_set[index].position,
                                                 self.problem_function(self.__interval_set[index].position))
                print(local_optimal_solution[index])
                if best_point > local_optimal_solution[index][1]:
                    best_point = local_optimal_solution[index][1]
                    self.optimal_solution = local_optimal_solution[index]

    def solution(self):
        self.init_population()
        self.search_local_optimal_solution()
        self.print_local_optimal_solution()

        print('the optimal solution is', self.optimal_solution)


if __name__ == "__main__":
	algo = AntColonyOptimization()
	algo.solution()

  蚁群算法可视化效果如下：

代码仓库：IALib[GitHub]

  本篇代码已同步至【智能算法(Python复现)】专栏专属仓库：IALib
  运行IALib库中的ACO算法：

git clone git@github.com:xiayouran/IALib.git
cd examples
python main.py -algo aco