数学建模算法之粒子群算法_粒子群模型与算法如何结合

数学建模之粒子群算法

粒子群优化算法(PSO)是一种进化计算技术(evolutionary computation)，1995 年由Eberhart 博士和kennedy 博士提出，源于对鸟群捕食的行为研究 。该算法最初是受到飞鸟集群活动的规律性启发，进而利用群体智能建立的一个简化模型。粒子群算法在对动物集群活动行为观察基础上，利用群体中的个体对信息的共享使整个群体的运动在问题求解空间中产生从无序到有序的演化过程，从而获得最优解。
粒子公式：

在找到这两个最优值时，粒子根据如下的公式来更新自己的速度和新的位置：
v[] = w * v[] + c1 * rand() * (pbest[] - present[]) + c2 * rand() * (gbest[] - present[]) (a)

present[] = present[] + v[] (b)

v[] 是粒子的速度, w是惯性权重,present[] 是当前粒子的位置. pbest[] and gbest[] 如前定义. rand () 是介于（0， 1）之间的随机数. c1, c2 是学习因子. 通常 c1 = c2 = 2.
这里我们主要讲解粒子群算法解决TSP问题
旅行商问题，即TSP问题（Traveling Salesman Problem）又译为旅行推销员问题、货郎担问题，是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。

例：根据city.xlsx文件，遍历51座城市找到最短路径最优值

解：算法的实现
粒子的表示：
1.TSP问题的一个解为一个序列，可以表示为一个粒子；
2.速度的表示：用一个序列的交换序列表示粒子的速度。
适应度函数的定义：
3.当前序列的路径长度即为适应度值，通过经纬度坐标计算。
惯性因子的定义：
4.自身的交换序列即惯性因子。
该MATLAB代码部分来自MATLAB智能算法分析
实现如下：

function indiFit=fitness(x,City ,Citydistance)% 该函数用于计算个体适应度值
m=size(x,1);
n=size(City,1);
indiFit=zeros(m,1);
for i=1:m
    for j=1:n-1
        indiFit(i)=indiFit(i)+ Citydistance(x(i,j),x(i,j+1));
    end
    indiFit(i)=indiFit(i)+ Citydistance(x(i,1),x(i,n));
end

function dist=dist(x,D)% 该函数用于计算总距离
n=size(x,2);
dist=0;
for i=1:n-1
    dist=dist+D(x(i),x(i+1));
end
dist=dist+D(x(1),x(n));

clc;clear

%导入数据 成3列矩阵，第一列为城市名，第二列为 x 坐标，第三列为 y坐标

cityCoor=[data(:,2) data(:,3)];%data为导入矩阵cityCoor城市坐标矩阵

figure
plot(cityCoor(:,1),cityCoor(:,2),'ms','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g')
legend('城市位置')
ylim([4 78])%限定y轴上限为78，下限为4


title('城市分布图','fontsize',12)
xlabel('km','fontsize',12)
ylabel('km','fontsize',12)
%ylim([min(cityCoor(:,2))-1 max(cityCoor(:,2))+1])

grid on %grid on是打开显示坐标轴网格线


n=size(cityCoor,1);            %城市数目
cityDist=zeros(n,n);           %城市距离矩阵
for i=1:n
    for j=1:n
        if i~=j
            cityDist(i,j)=((cityCoor(i,1)-cityCoor(j,1))^2+...
                (cityCoor(i,2)-cityCoor(j,2))^2)^0.5;%欧氏距离公式
        end
        cityDist(j,i)=cityDist(i,j);
    end
end
nMax=400;                      %进化次数
indiNumber=4000;               %个体数目
individual=zeros(indiNumber,n);
%^初始化粒子位置
for i=1:indiNumber
    individual(i,:)=randperm(n);    
end


indiFit=fitness(individual,cityCoor,cityDist);% 计算种群适应度
[value,index]=min(indiFit);
tourPbest=individual;                              %当前个体最优
tourGbest=individual(index,:) ;                    %当前全局最优
recordPbest=inf*ones(1,indiNumber);                %个体最优记录
recordGbest=indiFit(index);                        %群体最优记录
xnew1=individual;

L_best=zeros(1,nMax);% 循环寻找最优路径
for N=1:nMax
    
    %计算适应度值
    indiFit=fitness(individual,cityCoor,cityDist);
    
    %更新当前最优和历史最优
    for i=1:indiNumber
        if indiFit(i)<recordPbest(i)
            recordPbest(i)=indiFit(i);
            tourPbest(i,:)=individual(i,:);
        end
        if indiFit(i)<recordGbest
            recordGbest=indiFit(i);
            tourGbest=individual(i,:);
        end
    end
    
    [value,index]=min(recordPbest);
    recordGbest(N)=recordPbest(index);
    
  
    for i=1:indiNumber% 与个体最优进行交叉
       
        c1=unidrnd(n-1); %产生交叉位
        c2=unidrnd(n-1); %产生交叉位
        while c1==c2
            c1=round(rand*(n-2))+1;
            c2=round(rand*(n-2))+1;
        end
        chb1=min(c1,c2);
        chb2=max(c1,c2);
        cros=tourPbest(i,chb1:chb2);
        ncros=size(cros,2);      
        %删除与交叉区域相同元素
        for j=1:ncros
            for k=1:n
                if xnew1(i,k)==cros(j)
                    xnew1(i,k)=0;
                    for t=1:n-k
                        temp=xnew1(i,k+t-1);
                        xnew1(i,k+t-1)=xnew1(i,k+t);
                        xnew1(i,k+t)=temp;
                    end
                end
            end
        end
        %插入交叉区域
        xnew1(i,n-ncros+1:n)=cros;
        %新路径长度变短则接受
        dist=0;
        for j=1:n-1
            dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));
        end
        dist=dist+cityDist(xnew1(i,1),xnew1(i,n));
        if indiFit(i)>dist
            individual(i,:)=xnew1(i,:);
        end
        
        % 与全体最优进行交叉
        c1=round(rand*(n-2))+1;  %产生交叉位
        c2=round(rand*(n-2))+1;  %产生交叉位
        while c1==c2
            c1=round(rand*(n-2))+1;
            c2=round(rand*(n-2))+1;
        end
        chb1=min(c1,c2);
        chb2=max(c1,c2);
        cros=tourGbest(chb1:chb2); 
        ncros=size(cros,2);      
        %删除与交叉区域相同元素
        for j=1:ncros
            for k=1:n
                if xnew1(i,k)==cros(j)
                    xnew1(i,k)=0;
                    for t=1:n-k
                        temp=xnew1(i,k+t-1);
                        xnew1(i,k+t-1)=xnew1(i,k+t);
                        xnew1(i,k+t)=temp;
                    end
                end
            end
        end
        %插入交叉区域
        xnew1(i,n-ncros+1:n)=cros;
        %新路径长度变短则接受
        dist=0;
        for j=1:n-1
            dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));
        end
        dist=dist+cityDist(xnew1(i,1),xnew1(i,n));
        if indiFit(i)>dist
            individual(i,:)=xnew1(i,:);
        end
        
       
        c1=round(rand*(n-1))+1;   % 变异操作,产生变异位
        c2=round(rand*(n-1))+1;   %产生变异位
        while c1==c2
            c1=round(rand*(n-2))+1;
            c2=round(rand*(n-2))+1;
        end
        temp=xnew1(i,c1);
        xnew1(i,c1)=xnew1(i,c2);
        xnew1(i,c2)=temp;
        
        %新路径长度变短则接受
        dist=0;
        for j=1:n-1
            dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));
        end
        dist=dist+cityDist(xnew1(i,1),xnew1(i,n));
        if indiFit(i)>dist
            individual(i,:)=xnew1(i,:);
        end
    end

    [value,index]=min(indiFit);
    L_best(N)=indiFit(index);
    tourGbest=individual(index,:); 
    
end

figure
plot(L_best)
title('算法训练过程')
xlabel('迭代次数')
ylabel('适应度值')
grid on

figure
hold on
plot([cityCoor(tourGbest(1),1),cityCoor(tourGbest(n),1)],[cityCoor(tourGbest(1),2),...
    cityCoor(tourGbest(n),2)],'ms-','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g')
hold on
for i=2:n
    plot([cityCoor(tourGbest(i-1),1),cityCoor(tourGbest(i),1)],[cityCoor(tourGbest(i-1),2),...
        cityCoor(tourGbest(i),2)],'ms-','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g')
    hold on
end
legend('规划路径')
scatter(cityCoor(:,1),cityCoor(:,2));
title('规划路径','fontsize',10)
xlabel('km','fontsize',10)
ylabel('km','fontsize',10)

grid on
ylim([4 80])
% % 注意：代码仅供参考，不可直接用于自己的数模论文中
% % 国赛对于论文的查重要求非常严格，代码雷同也算作抄袭
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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217

该算法求解的答案不一定为最优解，需多次计算。
目前最短路径为
14 24 43 23 7 26 8 48 6 27 51 46 12 47 18 4 17 37 5 38 49 9 50 16 2 11 32 1 22 31 28 3 36 35 20 29 21 34 30 10 39 33 45 15 44 42 19 40 41 13 25
该路径长度为432.5745.