蚁群算法研究、基础代码实现及在路由中的应用

前言

为何选此题目

笔者简单浏览了一篇关于资源受限无线传感器网络中的多约束路由问题的论文，文中提到，多约束路由问题常常用计算智能方法来解决。而计算智能方法包括神经网络、机器学习、遗传算法、模糊计算、蚁群算法、人工鱼群算法、粒子群算法、免疫算法、禁忌搜索、进化算法、启发式算法、模拟退火算法、混合智能算法等等。笔者本次选用常常在论文和本子中看到的蚁群算法做一个基础的研究。

多目标优化的概念

在某个情景中在需要达到多个目标时，由于容易存在目标间的内在冲突，一个目标的优化是以其他目标劣化为代价，因此很难出现唯一最优解，取而代之的是在他们中间做出协调和折衷处理，使总体的目标尽可能的达到最优。
多目标优化问题可以描述为：

其中，k ≥ 2 ：优化的目标函数的个数；S：系统决策变量的可行域。

算法原理

蚂蚁在寻找食物源时，能在其走过的路上释放信息素，使得一定范围内的其他蚂蚁能够察觉到并由此影响它们以后的行为。同时，信息素还会随着时间逐渐回挥发。
当一些路径上通过的蚂蚁越来越多时，其留下的信息素也越来越多，以致信息素强度增大，所以蚂蚁选择该路径的概率也越高，从而更增大了该路径的信息素浓度。这其实是一种自催化的正反馈机制。

数学模型

1. 初始化
   包括信息素初始化，启发信息初始化，以及种群规模、信息素挥发率等参数初值的设置。
2. 解空间构建 （在问题空间依据状态转移规则如何构建候选解？）
   地点x和地点y之间的距离为$d_{xy}$。${\tau }_{xy}$为xy路径上的信息素浓度，$P_{xy}^k$表示蚂蚁k从x到y的概率。蚂蚁去到哪个地点，受到两方面因素影响：地点吸引力与路径信息素水平。
   
   其中${\eta }_{xy}$启发函数，表示蚂蚁从x到y的期望程度，通常取xy距离的倒数。
3. 信息素更新
   信息素更新包括两个环节：
   • 信息素挥发
     设$\rho$为挥发系数，用于降低路径上的信息素，减小信息素对未来蚂蚁行为的影响，增加算法的探索能力；
   • 信息素释放，蚂蚁在其所经过的路径上释放信息素，加强对蚂蚁未来选择该路径的影响，增强算法的开发能力。
     
     
     $L_k$：第k只蚂蚁本次旅行的成本(通常是长度)；
     $Q$：蚂蚁一次旅行释放的信息素总量，为常数。

算法代码

MATLAB

来源于https://www.bilibili.com/video/BV1ZA411v7pC?from=search&seid=6994083714146285315&spm_id_from=333.337.0.0，部分注释系笔者自己添加

%% I. 清空环境
clc
clear all
%% II. 符号说明
% C                         -- n个城市的坐标
% NC_max                    -- 最大迭代次数
% m                         -- 蚁群中蚂蚁的数量，一般设置为城市的1.5倍
% D(i, j)                   -- 两城市i和之间的距离
% Eta(i, j) = 1 ./ D(i, j)  -- 启发函数
% alpha                     -- 表征信息素重要程度的参数
% beta                      -- 表征启发函数重要程度的参数
% rho                       -- 信息素挥发因子
% Q                         --
% rBest                     -- 各次迭代最佳的路线
% lBest                     -- 各次迭代最佳路线的长度
% lAverage                  --各次迭代的平均长度
 
%% III. 导入城市位置数据
citys = [18.4700   95.1000
     16.4700   94.6400
     20.0900   94.5400
     14.3900   93.3700
     25.2300   97.2400
     22.0000   93.0500
     23.4700   92.0200
     16.2000   96.2900
     17.3000   97.3800
     13.0500   98.1200
     15.5300   97.3800
     24.5200   95.5900
     16.4100   97.1300
     15.0900   92.5500];
 
%% IV. 计算距离矩阵
D = Distance(citys);                         % 计算距离矩阵
n = size(D, 1);                              % 城市的个数
  
%% V. 初始化参数
NC_max = 200;                                % 最大迭代次数，取100~500之间
m = 22;                                      % 蚂蚁的个数，一般设为城市数量的1.5倍
alpha = 1;                                   % α 选择[1, 4]比较合适
beta = 4;                                    % β 选择[3 4 5]比较合适
rho = 0.2;                                   % ρ 选择[0.1, 0.2, 0.5]比较合适
Q = 20;
NC = 1;                                      % 迭代次数，一开始为1
 
Eta = 1 ./ D;                                % η = 1 / D(i, j) ,这里是矩阵
Tau = ones(n, n);                            % Tau(i, j)表示边(i, j)的信息素量，一开始都为1
Table = zeros(m, n);                         % 路径记录表
 
rBest = zeros(NC_max, n);                    % 记录各代的最佳路线
lBest = inf .* ones(NC_max, 1);              % 记录各代的最佳路线的总长度
lAverage = zeros(NC_max, 1);                 % 记录各代路线的平均长度
 
%% VI. 迭代寻找最佳路径
while NC <= NC_max
    % 第1步：随机产生各个蚂蚁的起点城市
    start = zeros(m, 1);        % m行1列的零向量
    for i = 1: m
        temp = randperm(n);     % 生成随机含有1-n这n个元素的行向量temp
        start(i) = temp(1);
    end
    Table(:, 1) = start;        % Table路径记录表的第一列即是所有蚂蚁的起点城市，共m行，对应m只蚂蚁的起点城市
    citys_index = 1: n;         % 所有城市索引的一个集合
    % 第2步：逐个蚂蚁路径选择
    for i = 1: m
        % 逐个城市路径选择
        for j = 2: n
            tabu = Table(i, 1: (j - 1));                    % 蚂蚁i已经访问的城市集合（称禁忌表）
            allow_index = ~ismember(citys_index, tabu);     % citys_index在集合tabu中相应的位置为0，其余为1，生成的allow_index为逻辑数组
            Allow = citys_index(allow_index);               % 用逻辑数组allow_index来剔除citys_index中与allow_index中为0的对应位置上的元素。Allow表：存放待访问的城市
            P = Allow;
             
            % 计算从城市j到剩下未访问的城市的转移概率
            for k = 1: size(Allow, 2)                       % 待访问的城市数量
                P(k) = Tau(tabu(end), Allow(k))^alpha * Eta(tabu(end), Allow(k))^beta;
            end
            P = P / sum(P);                                 % 归一化
             
            % 轮盘赌法选择下一个访问城市（为了增加随机性）
            Pc = cumsum(P);    % Pc为P的累计和
            target_index = find(Pc >= rand);    % rand为生成0~1间的随机数，find返回的是Pc中所有>=随机数的位置的索引值
            target = Allow(target_index(1));    % 返回所有>=随机数的位置当中的第一个位置
            Table(i, j) = target;               % 路径表中第i行第j列就代表第i只蚂蚁第j个到达的城市是哪个 
        end
    end
     
    % 第3步：计算各个蚂蚁的路径距离
    length = zeros(m, 1);
    for i = 1: m
        Route = Table(i, :);
        for j = 1: (n - 1)
            length(i) = length(i) + D(Route(j), Route(j + 1));
        end
        length(i) = length(i) + D(Route(n), Route(1));
    end
     
    % 第4步：计算最短路径距离及平均距离
    if NC == 1
        [min_Length, min_index] = min(length);
        lBest(NC) = min_Length;
        lAverage(NC) = mean(length);
        rBest(NC, :) = Table(min_index, :);
    else
        [min_Length, min_index] = min(length);
        lBest(NC) = min(lBest(NC - 1), min_Length);  % 若本次迭代的最短路径长度比上一次长，则还是取上一次为最短路径，若比上一次短，则取本次
        lAverage(NC) = mean(length);
        if lBest(NC) == min_Length                   % 某次迭代的最短路径和已经得到的最短路径长度相同，即这次迭代的最短路径长度比之前任何一次的都短，则取此次迭代的路径为最短路径
            rBest(NC, :) = Table(min_index, :);
        else
            rBest(NC, :) = rBest((NC - 1), :);       % 否则，取上一次的路径为最短路径
        end
    end
    % 第5步：更新信息素
    Delta_tau = zeros(n, n);
    for i = 1: m         
        for j = 1: (n - 1)
            Delta_tau(Table(i, j), Table(i, j + 1)) = Delta_tau(Table(i, j), Table(i, j + 1)) + Q / length(i);
        end
        Delta_tau(Table(i, n), Table(i, 1)) = Delta_tau(Table(i, n), Table(i, 1)) + Q / length(i);
    end
    Tau = (1 - rho) .* Tau + Delta_tau;
 
    % 第6步：迭代次数加1，并且清空路径记录表
    NC = NC + 1;
    Table = zeros(m, n);
end
%% VII. 结果显示
[shortest_Length, shortest_index] = min(lBest);
shortest_Route = rBest(shortest_index, :);
disp(['最短距离:' num2str(shortest_Length)]);
disp(['最短路径:' num2str([shortest_Route shortest_Route(1)])]);   %是个循环，最终要回到起点
 
%% VIII. 绘图
figure(1)
plot([citys(shortest_Route,1); citys(shortest_Route(1),1)],...
     [citys(shortest_Route,2); citys(shortest_Route(1),2)],'o-');
grid on
for i = 1: size(citys, 1)
    text(citys(i, 1), citys(i, 2), ['   ' num2str(i)]);
end
text(citys(shortest_Route(1), 1), citys(shortest_Route(1), 2), '       起点');
text(citys(shortest_Route(end), 1), citys(shortest_Route(end), 2), '       终点');
xlabel('城市位置横坐标')
ylabel('城市位置纵坐标')
title(['蚁群算法优化路径(最短距离: ' num2str(shortest_Length) ')'])
figure(2)
plot(1: NC_max, lBest, 'b', 1: NC_max, lAverage, 'r:')
legend('最短距离', '平均距离')
xlabel('迭代次数')
ylabel('距离')
title('各代最短距离与平均距离对比')

其中用到的Distance函数（Distance.m）如下：

function D = Distance(citys)
%% 计算两两城市之间的距离
% 输入：各城市的位置坐标（citys）
% 输出：两两之间的距离(D)

n = size(citys,1);  % n为citys的第一个维度的长度，也就是行数，也就是城市的个数
D = zeros(n,n);
for i=1:n
    for j = i + 1 : n
        D(i,j)=sqrt((citys(i,1)-citys(j,1))^2+(citys(i,2)-citys(j,2))^2);
        D(j,i) = D(i,j);
    end
    D(i,i)=1e-4;     %近似于零，因为后面启发因子beta要取倒数，因此用一个很小的数字代替0
end

运行结果

最短距离:35.6082
最短路径:2 1 3 6 7 12 5 9 13 8 11 10 4 14 2

Python

来源于https://www.cnblogs.com/xxhbdk/p/9177423.html，部分注释系笔者自己添加

import numpy as np
import matplotlib.pyplot as plt


# 建立“蚂蚁”类
class Ant(object):
    def __init__(self, path):
        self.path = path                       # 蚂蚁当前迭代整体路径
        self.length = self.calc_length(path)   # 蚂蚁当前迭代整体路径长度

    def calc_length(self, path_):              # path=[A, B, C, D, A]注意路径闭环
        length_ = 0
        for i in range(len(path_)-1):
            delta = (path_[i].x - path_[i+1].x, path_[i].y - path_[i+1].y)
            length_ += np.linalg.norm(delta)   # 求delta的二范数
        return length_

    @staticmethod
    def calc_len(A, B):                        # 静态方法，计算城市A与城市B之间的距离
        return np.linalg.norm((A.x - B.x, A.y - B.y))   # ✓[(A.x-B.x)^2+(A.y - B.y)^2]，求直线距离


# 建立“城市”类
class City(object):
    def __init__(self, x, y):
        self.x = x
        self.y = y


# 建立“路径”类
class Path(object):
    def __init__(self, A):                     # A为起始城市
        self.path = [A, A]

    def add_path(self, B):                     # 追加路径信息，方便计算整体路径长度
        self.path.append(B)
        self.path[-1], self.path[-2] = self.path[-2], self.path[-1]


# 构建“蚁群算法”的主体
class ACO(object):
    def __init__(self, ant_num=50, maxIter=300, alpha=1, beta=5, rho=0.1, Q=1):
        self.ants_num = ant_num   # 蚂蚁个数
        self.maxIter = maxIter    # 蚁群最大迭代次数
        self.alpha = alpha        # 信息启发式因子
        self.beta = beta          # 期望启发式因子
        self.rho = rho            # 信息素挥发速度
        self.Q = Q                # 信息素强度
        ###########################
        self.deal_data('coordinates.dat')                         # 提取所有城市的坐标信息
        ###########################
        self.path_seed = np.zeros(self.ants_num).astype(int)      # 记录一次迭代过程中每个蚂蚁的初始城市下标
        self.ants_info = np.zeros((self.maxIter, self.ants_num))  # 记录每次迭代后所有蚂蚁的路径长度信息
        self.best_path = np.zeros(self.maxIter)                   # 记录每次迭代后整个蚁群的“历史”最短路径长度
        ###########################
        self.solve()              # 完成算法的迭代更新
        self.display()            # 数据可视化展示

    def deal_data(self, filename):
        with open(filename, 'rt') as f:
            temp_list = list(line.split() for line in f)
            # 临时存储提取出来的坐标信息，每行是列表temp_list中的一个元素
        self.cities_num = len(temp_list)                                                   # 1. 获取城市个数
        self.cities = list(City(float(item[0]), float(item[1])) for item in temp_list)     # 2. 构建城市列表
        self.city_dist_mat = np.zeros((self.cities_num, self.cities_num))                  # 3. 构建城市距离矩阵（初始化为全零）
        # 构建城市距离矩阵过程：
        for i in range(self.cities_num):
            A = self.cities[i]
            for j in range(i, self.cities_num):
                B = self.cities[j]
                self.city_dist_mat[i][j] = self.city_dist_mat[j][i] = Ant.calc_len(A, B)
        self.phero_mat = np.ones((self.cities_num, self.cities_num))                       # 4. 初始化信息素矩阵
        # self.phero_upper_bound = self.phero_mat.max() * 1.2                              ###信息素浓度上限
        self.eta_mat = 1/(self.city_dist_mat + np.diag([np.inf]*self.cities_num))          # 5. 初始化启发函数矩阵

    def solve(self):
        iterNum = 0                                                            # 当前迭代次数
        while iterNum < self.maxIter:
            self.random_seed()                                                 # 使整个蚁群产生随机的起始点
            delta_phero_mat = np.zeros((self.cities_num, self.cities_num))     # 初始化每次迭代后信息素矩阵的增量
            ##########################################################################
            for i in range(self.ants_num):
                city_index1 = self.path_seed[i]                                # 每只蚂蚁访问的第一个城市下标
                ant_path = Path(self.cities[city_index1])                      # 记录每只蚂蚁访问过的城市
                tabu = [city_index1]                                           # 记录每只蚂蚁访问过的城市下标，禁忌城市下标列表
                non_tabu = list(set(range(self.cities_num)) - set(tabu))
                for j in range(self.cities_num-1):                             # 对余下的城市进行访问
                    up_proba = np.zeros(self.cities_num-len(tabu))             # 初始化状态迁移概率的分子
                    for k in range(self.cities_num-len(tabu)):
                        up_proba[k] = np.power(self.phero_mat[city_index1][non_tabu[k]], self.alpha) * \
                        np.power(self.eta_mat[city_index1][non_tabu[k]], self.beta)
                    proba = up_proba/sum(up_proba)                             # 每条可能子路径上的状态迁移概率
                    while True:                                                # 提取出下一个城市的下标
                        random_num = np.random.rand()
                        index_need = np.where(proba > random_num)[0]
                        if len(index_need) > 0:
                            city_index2 = non_tabu[index_need[0]]
                            break
                    ant_path.add_path(self.cities[city_index2])
                    tabu.append(city_index2)
                    non_tabu = list(set(range(self.cities_num)) - set(tabu))
                    city_index1 = city_index2
                self.ants_info[iterNum][i] = Ant(ant_path.path).length
                if iterNum == 0 and i == 0:                                    # 完成对最佳路径城市的记录
                    self.best_cities = ant_path.path
                else:
                    if self.ants_info[iterNum][i] < Ant(self.best_cities).length: self.best_cities = ant_path.path
                tabu.append(tabu[0])                                           # 每次迭代完成后，使禁忌城市下标列表形成完整闭环
                for l in range(self.cities_num):
                    delta_phero_mat[tabu[l]][tabu[l+1]] += self.Q/self.ants_info[iterNum][i]

            self.best_path[iterNum] = Ant(self.best_cities).length

            self.update_phero_mat(delta_phero_mat)                             # 更新信息素矩阵
            iterNum += 1

    def update_phero_mat(self, delta):
        self.phero_mat = (1 - self.rho) * self.phero_mat + delta
        # self.phero_mat = np.where(self.phero_mat > self.phero_upper_bound, self.phero_upper_bound, self.phero_mat) # 判断是否超过浓度上限

    def random_seed(self):                                                     # 产生随机的起始点下表，尽量保证所有蚂蚁的起始点不同
        if self.ants_num <= self.cities_num:                                   # 蚂蚁数 <= 城市数
            self.path_seed[:] = np.random.permutation(range(self.cities_num))[:self.ants_num]
        #     np.random.permutation表示随机排序
        else:                                                                  # 蚂蚁数 > 城市数
            self.path_seed[:self.cities_num] = np.random.permutation(range(self.cities_num))
            temp_index = self.cities_num
            while temp_index + self.cities_num <= self.ants_num:
                self.path_seed[temp_index:temp_index + self.cities_num] = np.random.permutation(range(self.cities_num))
                temp_index += self.cities_num
            temp_left = self.ants_num % self.cities_num
            if temp_left != 0:
                self.path_seed[temp_index:] = np.random.permutation(range(self.cities_num))[:temp_left]

    def display(self):                                                         # 数据可视化展示
        plt.figure(figsize=(6, 10))
        plt.subplot(211)
        plt.plot(self.ants_info, 'g.')
        plt.plot(self.best_path, 'r-', label='history_best')
        plt.xlabel('Iteration')
        plt.ylabel('length')
        plt.legend()
        plt.subplot(212)
        plt.plot(list(city.x for city in self.best_cities), list(city.y for city in self.best_cities), 'g-')
        plt.plot(list(city.x for city in self.best_cities), list(city.y for city in self.best_cities), 'r.')
        plt.xlabel('x')
        plt.ylabel('y')
        plt.savefig('ACO.png', dpi=500)
        plt.show()
        plt.close()
        
ACO()

在本主程序同一文件夹下添加数据文件coordinates.dat，记录每个点的坐标位置，内容如下：

565.0   575.0
25.0   185.0
345.0   750.0
945.0   685.0
845.0   655.0
880.0   660.0
25.0    230.0
525.0   1000.0
580.0   1175.0
650.0   1130.0
1605.0   620.0
1220.0   580.0
1465.0   200.0
1530.0   5.0
845.0   680.0
725.0   370.0
145.0   665.0
415.0   635.0
510.0   875.0
560.0   365.0
300.0   465.0
520.0   585.0
480.0   415.0
835.0   625.0
975.0   580.0
1215.0   245.0
1320.0   315.0
1250.0   400.0
660.0   180.0
410.0   250.0
420.0   555.0
575.0   665.0
1150.0   1160.0
700.0   580.0
685.0   595.0
685.0   610.0
770.0   610.0
795.0   645.0
720.0   635.0
760.0   650.0
475.0   960.0
95.0   260.0
875.0   920.0
700.0   500.0
555.0   815.0
830.0   485.0
1170.0   65.0
830.0   610.0
605.0   625.0
595.0   360.0
1340.0   725.0
1740.0   245.0

运行结果

算法在路由中的应用

1. 初始化各通信节点之间的路径，通过节点广播信息完成，广播数据帧格式如下：
   
   • ID：每个节点被分配不同的ID
   • Positionx，Positiony：节点位置坐标（x，y）
   • 源节点：若该节点有数据发送需求，则该字段置1，同时将目的节点字段置为目的节点ID，并在发送序列字段对发送序列进行编号；若没有数据发送需求，则该字段置0
2. 将若干蚂蚁放在不同的通信节点上，每个通信节点维护自身的信息素列表，如下图所示：
   
   此表描述了当前节点的信息素含量和此刻与其他节点的距离。当节点位置移动后，此表中数据也进行更新。
3. 每只蚂蚁根据各节点至目的节点的距离d和信息素水平${\tau }_{ij}(t)$选择下一跳的节点，同时修改禁忌表tabu。
4. 所有蚂蚁都完成循环后，更新信息素列表中的信息素水平和节点位置信息。
5. 返回步骤二，迭代次数+1，直至满足结束条件，找到源节点与目的节点之间的最优路径（即路程最短），摒弃其他不必要的路径信息