这一专栏将介绍电机控制算法,为接下来的课题研究提供参考。
本文将介绍群智能算法中的蚁群算法,它是历史最悠久的群智能算法之一。
本文首先结合旅行商问题介绍蚁群算法的原理,然后将蚁群算法应用于电机转速环PID控制器的参数优化。
蚁群算法概述
蚁群算法(Ant Colony Optimization, ACO)是一种模拟蚂蚁群体觅食行为的群智能优化算法,该算法采用的分布式并行计算机制,易于与其它方法结合,具有较强的鲁棒性,适用于组合优化问题的求解。
首先简单总结一下蚁群觅食的现象:
- 蚂蚁在运动时,会在路径上释放出一种称为信息素(Pheromone)的物质,并且能够感知这种物质的存在及其强弱,并根据信息素的浓度指导运动的方向。
- 当遇到没有走过的路口时,蚂蚁会随机挑选一个方向前进。
- 蚂蚁走的路径越长则在该路径上的单位长度上的信息素越小。
- 当后来的蚂蚁再次遇到这个路口时,选择信息素强度高的方向的概率就会增加,这样就形成了一个正反馈机制。
- 最短路径上的信息素的强度越来越大,而其它路径上的信息素的强度随时间的流逝而逐渐消失。最终寻找到最短的路径。
- 当原先的最短路径上出现障碍物时,蚂蚁能够重新寻找另外一条最短路径,这体现了蚂蚁对外部环境的适应性(鲁棒性)。
旅行商问题(Travelling Salesman Problem,TSP):
假设地图上有n个城市,n个城市间均有道路连接,不同城市间的道路长度不同,这使得不同城市间的交通成本不同。
一个商人需要从某个城市出发,按照一定的顺序访问这n个城市,且每个城市恰好访问一次,最终回到出发城市。
目标是求出一个最优的访问顺序,使得总交通成本最低。

节点与路径:假设蚂蚁由起点S出发,中途经过节点A、B、C,最后到达终点D,那么在这一过程中的SA、AB、BC、CD均为所经过的路径(路径是相邻两点之间构成的线段)。
在旅行商问题中,各个城市是节点,城市间的道路是路径。
信息素:蚂蚁在经过路径之时,会在路径上留下信息素,信息素浓度越高,意味着蚂蚁越愿意选择这一路径(形成正反馈)。
可以用\(\tau_{ij}\)表示节点\(i\)到节点\(j\)之间的信息素。
信息素的更新:每一轮迭代后都会更新信息素,其表达式为:
\[\tau_{ij}(t+1)=(1-\rho)\tau_{ij}(t)+\Delta
\tau_{ij}\] 式中\(\rho \in
(0,1)\)为信息素的挥发系数,这意味着若没有新的蚂蚁走过此路径,则该路径上的信息素浓度会不断降低。
\(\Delta \tau_{ij}=\sum\limits_{k = 1}^m
\Delta
\tau_{ij}^k\),为信息素在本次迭代的增量,是各个蚂蚁在本次迭代中释放的信息素之和。
单个蚂蚁释放的信息素:主要分为三种模型。
首先是ant-cycle模型,其公式为\(\Delta
\tau_{ij}^k = \left\{ {\begin{array}{*{20}{c}} \frac{Q}{L_k} &&
若第k只蚂蚁在本次迭代经过了路径ij \\ 0 && 其他 \end{array}}
\right.\)
式中\(Q\)是一个正数,代表信息素的浓度,\(L_k\)为第\(k\)只蚂蚁在本轮迭代所经过路径的总长度。
(总长度越低,意味着这条路越好走,释放的信息素就越多)
将\(\frac{Q}{L_k}\)换成\(\frac{Q}{d_{ij}}\)就是ant-quantity模型(\(d_{ij}\)是路径\(ij\)的长度),再换成\(Q\)就是ant-density模型。
ant-cycle模型利用全局信息更新路径上的信息素量,其优化性能要好于采用局部信息的ant-quantity模型与ant-density模型。
状态转移概率公式:在第\(t\)次迭代,蚂蚁位于\(i\)节点时,选择\(j\)节点的概率如下 \[p_{ij}^k(t)= \left\{
{\begin{array}{*{20}{c}} \frac{[\tau_{ij}(t)]^\alpha \cdot
[\eta_{ij}(t)]^\beta}{\sum\limits_{s\in J_k(i)}[\tau_{is}]^\alpha \cdot
[\eta_{is}]^\beta} && 当j \in J_k(i)时 \\ 0 &&
其他 \end{array}} \right.\] 式中\(J_k(i)\)表示蚂蚁下一步可以选择的所有节点,在程序运行过程中,需要定义一个禁忌表\(tabu_k\),该表存储当前蚂蚁已经走过的节点,\(J_k(i)\)是所有节点去除\(tabu_k\)后剩余的部分。
\(\eta_{ij}\)是启发式因子,通常取\(\eta_{ij}=\frac{1}{d_{ij}}\)。
\(\alpha\)与\(\beta\)是信息素与启发式因子的权重系数,体现了信息素和启发式因子的相对重要程度。
该公式采用的是轮盘赌算法,分子对应单个节点,分母对应所有节点的和。路径上信息素浓度越大或路径越短时,节点被选中的概率越大。
蚁群算法应用于旅行商问题
下面是将蚁群算法应用于解决旅行商问题的matlab实例,代码来自智能优化算法及其MATLAB实例(第2版)例5.1。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119%%%%%%%%%%%%%%%%%%%%蚁群算法解决TSP问题%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; %清除所有变量
close all; %清图
clc; %清屏
m=50; %蚂蚁个数
Alpha=1; %信息素重要程度参数
Beta=5; %启发式因子重要程度参数
Rho=0.1; %信息素蒸发系数
G_max=200; %最大迭代次数
Q=100; %信息素增加强度系数
C=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;...
3238 1229;4196 1044;4312 790;4386 570;3007 1970;2562 1756;...
2788 1491;2381 1676;1332 695;3715 1678;3918 2179;4061 2370;...
3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;...
3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;...
2370 2975]; %31个省会城市坐标
%%%%%%%%%%%%%%%%%%%%%%%%第一步:变量初始化%%%%%%%%%%%%%%%%%%%%%%%%
n=size(C,1); %n表示问题的规模(城市个数)
D=zeros(n,n); %D表示两个城市距离间隔矩阵
for i=1:n
for j=1:n
if i~=j
D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;
else
D(i,j)=eps;
end
D(j,i)=D(i,j);
end
end
Eta=1./D; %Eta为启发因子,这里设为距离的倒数
Tau=ones(n,n); %Tau为信息素矩阵
Tabu=zeros(m,n); %存储并记录路径的生成
NC=1; %迭代计数器
R_best=zeros(G_max,n); %各代最佳路线
L_best=inf.*ones(G_max,1); %各代最佳路线的长度
figure(1);%优化解
while NC<=G_max
%%%%%%%%%%%%%%%%%%第二步:将m只蚂蚁放到n个城市上%%%%%%%%%%%%%%%%
Randpos=[];
for i=1:(ceil(m/n))
Randpos=[Randpos,randperm(n)];
end
Tabu(:,1)=(Randpos(1,1:m))';
%%%%%第三步:m只蚂蚁按概率函数选择下一座城市,完成各自的周游%%%%%%
for j=2:n
for i=1:m
visited=Tabu(i,1:(j-1)); %已访问的城市
J=zeros(1,(n-j+1)); %待访问的城市
P=J; %待访问城市的选择概率分布
Jc=1;
for k=1:n
if length(find(visited==k))==0
J(Jc)=k;
Jc=Jc+1;
end
end
%%%%%%%%%%%%%%%%%%计算待选城市的概率分布%%%%%%%%%%%%%%%%
for k=1:length(J)
P(k)=(Tau(visited(end),J(k))^Alpha)...
*(Eta(visited(end),J(k))^Beta);
end
P=P/(sum(P));
%%%%%%%%%%%%%%%%按概率原则选取下一个城市%%%%%%%%%%%%%%%%
Pcum=cumsum(P);
Select=find(Pcum>=rand);
to_visit=J(Select(1));
Tabu(i,j)=to_visit;
end
end
if NC>=2
Tabu(1,:)=R_best(NC-1,:);
end
%%%%%%%%%%%%%%%%%%%第四步:记录本次迭代最佳路线%%%%%%%%%%%%%%%%%%
L=zeros(m,1);
for i=1:m
R=Tabu(i,:);
for j=1:(n-1)
L(i)=L(i)+D(R(j),R(j+1));
end
L(i)=L(i)+D(R(1),R(n));
end
[L_best(NC),pos]=min(L);
R_best(NC,:)=Tabu(pos(1),:);
%%%%%%%%%%%%%%%%%%%%%%%%%第五步:更新信息素%%%%%%%%%%%%%%%%%%%%%%
Delta_Tau=zeros(n,n);
for i=1:m
for j=1:(n-1)
Delta_Tau(Tabu(i,j),Tabu(i,j+1))=...
Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i);
end
Delta_Tau(Tabu(i,n),Tabu(i,1))=...
Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i);
end
Tau=(1-Rho).*Tau+Delta_Tau;
%%%%%%%%%%%%%%%%%%%%%%%第六步:禁忌表清零%%%%%%%%%%%%%%%%%%%%%%
Tabu=zeros(m,n);
%%%%%%%%%%%%%%%%%%%%%%%%%历代最优路线%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n-1
plot([ C(R_best(NC,i),1), C(R_best(NC,i+1),1)],...
[C(R_best(NC,i),2), C(R_best(NC,i+1),2)],'bo-');
hold on;
end
plot([C(R_best(NC,n),1), C(R_best(NC,1),1)],...
[C(R_best(NC,n),2), C(R_best(NC,1),2)],'ro-');
title(['优化最短距离:',num2str(L_best(NC))]);
hold off;
pause(0.005);
NC=NC+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%第七步:输出结果%%%%%%%%%%%%%%%%%%%%%%%%%%
Pos=find(L_best==min(L_best));
Shortest_Route=R_best(Pos(1),:); %最佳路线
Shortest_Length=L_best(Pos(1)); %最佳路线长度
figure(2),
plot(L_best)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')
参数初始化:
定义蚁群个体数量\(m=50\),迭代次数\(G_{max}=200\)(迭代的循环变量是NC
)。
定义单个蚂蚁释放信息素多少相关的\(Q\),状态转移概率公式中的信息素与启发式因子的权重系数\(\alpha\)与\(\beta\),信息素更新公式中的挥发系数\(\rho\)。
\(C\)是一个\(31 \times
2\)矩阵,用于存储平面上所有节点的坐标。
n=size(C,1);
将矩阵\(C\)的行数赋值给变量\(n\),\(n\)表示节点(城市)的数量。
接下来建立了一个\(31 \times
31\)的矩阵\(D\)来存储节点之间所构成的路径的距离,用到了距离公式\(d=\sqrt{(x_1-x_2)^2+(y_1-y_2)^2}\)。
这里使用eps
来定义同一个节点间的距离(\(i=j\)时的\(D(i,j)\)),这是个matlab函数,返回值为一个很小的正数,其值为\(2^{-52}\)。
另外定义\(D(j,i)=D(i,j)\),说明这是一个对称的TSP问题,两个节点在两个方向(从\(i\)到\(j\)与从\(j\)到\(i\))的距离相等。若是非对称TSP问题,则会有\(D(j,i) \ne D(i,j)\)。
定义启发因子\(\eta_{ij}=\frac{1}{d_{ij}}\)。
定义矩阵\(\tau_{n \times
n}\),存储各个路径的信息素。
定义矩阵\(tabu_{m \times
n}\),也就是禁忌表,存储一次迭代过程中\(m\)只蚂蚁各自所选择的路线(经过节点的顺序)。
定义矩阵\(R_{best}\),储存各次迭代的最优路线。
定义矩阵\(L_{best}\),储存各次迭代最优线路的总长度。定义中用到了inf
函数,使得其默认值为无穷大。
自此初始化工作全部结束,下面进行的是循环迭代过程。
将m只蚂蚁放到n个城市上:
首先确定蚂蚁的初始位置,定义Randpos
用于储存随机生成的数据,该数组的元素个数并不固定,因此默认设为空。
然后进入一个for循环,cell
函数是向上取整,这里\(m=50,n=31\),因此循环次数为2。
randperm(n)
函数会返回一个包含从1到\(n\)的随机序列,根据代码所写,每次循环都会使Randpos
增加\(n\)个互不相同的元素。
对于本程序,Randpos
的最终会有62个元素,前31个互不相同,后31个也互不相同。
最后Randpos
中的前\(m\)个元素会赋值给禁忌表\(tabu\)的第一列。
相较于直接生成1到\(n\)之间的随机数,这样可以确保当\(m \ge n\)时,每个城市上至少有一只蚂蚁。
蚂蚁按照状态转移概率公式经过各个节点:
外循环\(j\)从2到\(n\),共\((n-1)\)次,每次循环对应一次节点的选择,对应起点之外的另外\((n-1)\)个节点。
内循环\(i\)对应每一只蚂蚁。
首先将这个蚂蚁已经走过的节点存储于数组visited
中,共有\((j-1)\)个节点。
剩下的\((n-j+1)\)个节点就是未走过的节点,建立一个长度为\((n-j+1)\)的全零数组\(J\),用于存储这些节点。
然后将\(J\)赋值给\(P\)(使二者拥有相同的元素个数),这个\(P\)在后面将用于存储蚂蚁下一步要走的各个节点的概率。
接下来又进入一个循环\(k\),这个循环用于将所有未走过的节点填充到数组\(J\)中。
这里根据matlab给出的提示,将原文if(length(find(visited==k))==0)
修改为if(isempty(find(visited==k, 1)))
,以提高代码运行效率。
这个循环依次对各个节点进行判断,find(visited==k, 1)
会返回visited
数组中值为\(k\)(代表第\(k\)个节点)的第1个元素的元素下标(这个下标是几并不重要,重要的是能不能找到这个元素)。
find
函数的返回值是数组,若visited
数组中找不到值为\(k\)的元素,则会返回一个空数组,isempty
的返回值为1(真)。若找到了,isempty
的返回值为0(假)。
所以说,这个循环的作用是依次判断蚂蚁是否经过第\(k\)个节点,若未经过,则将这个节点存储到数组\(J\)中。
跳出这个小循环后,进入另一个小循环,该循环的变量也是\(k\)。这个循环按照状态转移概率公式得到蚂蚁下一步要走的各个节点的概率,最后储存于数组\(P\)。
最后一部分代码用于确定要去的下一个节点,cumsum
是累计求和函数,所返回的数组与输入的数组具有相同的元素个数,返回数组的第1个元素与输入相同,第2个元素是输入的前2个之和,第3个元素的输入的前3个之和,以此类推,最后一个元素是输入数组的所有元素值之和。
所以得到的数组\(P_{cum}\)的值会单调递增,最后一个元素的值为1。
rand
函数会生成一个在0到1之间的随机数,正好可以与\(P_{cum}\)的范围重合。
接下来的find(Pcum>=rand)
会寻找\(P_{cum}\)数组中大于0到1之间的随机数的所有元素,存储于数组\(Select\)中,而\(Select\)的第一个元素就是下一次迭代要经过的节点。
这里举一个例子,假设还剩下3个待走过的节点,根据状态转移概率公式计算得到的概率\(P=[0.1,0.3,0.6]\),那么通过累计求和后得到\(P_{cum}=[0.1,0.4,1]\)。
若\(0 < rand \le 0.1\),则数组\(Select\)将存储全部三个节点,最后要去的就是第一个节点;若\(0.1 <rand \le 0.4\),则数组\(Select\)将存储第二个与第三个节点,最后要去的是第二个节点;若\(0.4 <rand < 1\),则数组\(Select\)只存储第三个节点,最后要去的就是第三个节点。
接下来是一个判断,只有迭代到第二轮及以后才执行下一行代码。下一行代码将上一次迭代得到的最优路径赋值给禁忌表的第一行。
(这样写意味着从第二轮开始,蚁群中的1号蚂蚁的路径会一直被上一轮迭代产生的最优路径所替代)
在接下来的路径对比中,这个最优路径会与蚁群迭代产生的新路径进行对比,若新路径的总距离均大于这个最优路径,则这个最优路径依旧是最优路径;若新路径的总距离有小于这个最优路径的,新路径就会取代它成为新的最优路径。
记录本次迭代最佳路线: 列向量\(L\)存储本次迭代各个蚂蚁走过的总距离。
代码比较好理解,调用的是禁忌表\(tabu\)的数据,首先计算第一个节点到第二个节点,第二个节点到第三个节点,一直到最后一个节点间的距离,然后还要加上最后一个节点到第一个节点的距离(题目要求最后回到出发节点)。
然后将\(L\)中的最小值幅值给\(L_{best}\),作为这一次迭代蚂蚁走过总距离的最小值,并获得最小值对应的下标。
(这里将原文L_best(NC)=min(L);pos=find(L==L_best(NC));
修改为[L_best(NC),pos]=min(L);
)
最后将总距离最小的那只蚂蚁的路径赋值给\(R_{best}\),作为这一轮迭代的最优路径。
更新信息素:
按照上文所述公式更新路径间的信息素。
首先求单个蚂蚁释放的信息素,采用的是ant-cycle模型,外循环对应\(m\)个蚂蚁,内循环对应一只蚂蚁在本次迭代所经过的路径。
然后将其代入到信息素更新公式。
然后将禁忌表清零,为下一轮迭代做准备。
代码的最后一部分就不必多说了,绘制迭代过程中的图形,求出最优路径与最小的总距离,得到最小总距离随迭代次数变化的曲线,如下图所示。

蚁群算法应用于PID控制器参数优化
为了使得蚁群算法应用于PID控制器的参数优化,需要对算法本身进行一定的调整。对于旅行商问题,优化的目的是找到一条最短的路径以经过所有的节点。而对于电机PID控制器,其优化目标是寻找出一组最优的\(K_p,K_i,K_d\)参数,使得电机的控制性能最好。
信息素并不释放在路径上,而是在节点上,而节点就是待优化的\(K_p,K_i,K_d\)参数。
另外,需要定义与电机控制性能相关的评价函数,这相当于旅行商问题中的路径长度。

蚁群算法优化PID控制参数的路径图如上图所示,每一只蚂蚁需要依次经过\(K_p,K_i,K_d\)这三组节点,每一组节点选择其中的一个节点经过,最后迭代得到最优的一组\(K_p,K_i,K_d\)参数。
下面简要介绍所编写程序的思路:
首先是节点的取值,三个参数各有100个元素,构成的是等差数列。初步确定\(K_p\)为1:1:100,\(K_i\)为0.1:0.1:10,\(K_p\)为10:10:1000(首项:等差:末项)。
然后是信息素矩阵,由于信息素释放在节点上,因此可以定义一个\(100 \times 3\)矩阵。
为确保初始蚁群能够经过所有的节点,蚁群个体数量\(m=100\),初始蚁群的\(K_p,K_i,K_d\)均用randperm(m)
生成。
也就是说,在程序中节点取值个数与蚂蚁个数相等。
此外,该代码对状态转移概率公式做出调整,公式中只含\(\alpha\)项,没有\(\beta\)项,这是因为信息素全部在节点上,节点之间没有信息素。
1 | %蚁群算法优化PID控制器参数 |
蚁群算法的改进方向
该程序接下来的改进思路如下:
1.对于群智能算法,通常在迭代初期需要有尽可能大的搜索范围,避免陷入局部最优解;而迭代后期则需要减小搜索范围,以获得尽可能精确的解。可以根据迭代次数调整信息素浓度\(Q\),使之先大后小。或者限制信息素\(\tau\)的最大值与最小值。
2.增加信息素对周围节点的影响,目前信息素仅释放在蚂蚁走过的节点上,可以将释放的范围扩大至临近节点。比如蚂蚁走过\(K_p=50\)对应的节点时,同时将信息素释放到\(K_p=47 \sim
53\)范围内的其余节点,释放的量随距离的增大而减小。
3.PID参数是连续性的参数,而蚁群算法主要针对离散参数的优化,因此这里将PID参数离散化成等差数列,这样就降低了参数优化的精确性。
To Be Countinued未完待续......