【路径规划】基于matlab一种带交叉因子的双向寻优粒子群栅格地图路径规划【含Matlab源码 117期】

一、双向寻优粒子群及栅格地图简介

1 双向寻优粒子群简介
针对传统粒子群算法易早熟、精度低、后期收敛速度慢等问题,结合反向学习理论,提出了一种基于交叉因子的双向寻优粒子群优化算法(CBMPSO)。该算法使初始种群在搜索区域均匀分布,计算粒子及其反向粒子的适应值,取最优作为初始种群;迭代过程增加对全局最差粒子的跟踪,随机开启基于交叉因子的双向学习机制。对几种典型函数的测试结果表明,CBMPSO算法的寻优能力及收敛速度有了显著提高,并且能够有效避免早熟收敛问题。

2 栅格地图简介

2 栅格地图
2.1 栅格法应用背景
路径规划时首先要获取环境信息，建立环境地图，合理的环境表示有利于建立规划方法和选择合适的搜索算法，最终实现较少的时间开销而规划出较为满意的路径。一般使用栅格法在静态环境下建立环境地图。
2.2 栅格法实质
将AGV的工作环境进行单元分割，将其用大小相等的方块表示出来，这样栅格大小的选取是影响规划算法性能的一个很重要的因素。栅格较小的话，由栅格地图所表示的环境信息将会非常清晰，但由于需要存储较多的信息，会增大存储开销，同时干扰信号也会随之增加，规划速度会相应降低，实时性得不到保证；反之，由于信息存储量少，抗干扰能力有所增强，规划速随之增快，但环境信息划分会变得较为模糊，不利于有效路径的规划。在描述环境信息时障碍物所在区域在栅格地图中呈现为黑色，地图矩阵中标为1，可自由通行区域在栅格地图中呈现为白色，地图矩阵中标为0。路径规划的目的就是在建立好的环境地图中找到一条最优的可通行路径，所以使用栅格法建立环境地图时，栅格大小的合理设定非常关键。
2.3 10乘10的静态环境地图

10乘10的静态环境地图代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%建立环境地图%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function DrawMap(map)
n = size(map);
step = 1;
a = 0 : step :n(1);
b = 0 : step :n(2);
figure(1)
axis([0 n(2) 0 n(1)]); %设置地图横纵尺寸
set(gca,'xtick',b,'ytick',a,'GridLineStyle','-',...
'xGrid','on','yGrid','on');
hold on
r = 1;
for(i=1:n(1))         %设置障碍物的左下角点的x,y坐标for(j=1:n(2))if(map(i,j)==1)p(r,1)=j-1;p(r,2)=i-1;fill([p(r,1) p(r,1) + step p(r,1) + step p(r,1)],...[p(r,2) p(r,2) p(r,2) + step p(r,2) + step ],'k');r=r+1;hold onendend
end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%栅格数字标识%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_text = 1:1:n(1)*n(2); %产生所需数值.
for i = 1:1:n(1)*n(2)[row,col] = ind2sub([n(2),n(1)],i);text(row-0.9,col-0.5,num2str(x_text(i)),'FontSize',8,'Color','0.7 0.7 0.7');
end
hold on
axis square

建立环境矩阵，1代表黑色栅格，0代表白色栅格，调用以上程序，即可得到上述环境地图。

map=[0 0 0 1 0 0 1 0 0 0;1 0 0 0 0 1 1 0 0 0;0 0 1 0 0 0 1 1 0 0;0 0 0 0 0 0 0 0 0 0;0 0 0 0 0 1 0 0 1 0;1 0 0 0 0 1 1 0 0 0;0 0 0 1 0 0 0 0 0 0;1 1 1 0 0 0 1 0 0 0;0 0 0 0 0 1 1 0 0 0;0 0 0 0 0 1 1 0 0 0;];DrawMap(map);         %得到环境地图

2.4 栅格地图中障碍栅格处路径约束
移动体栅格环境中多采用八方向的移动方式，此移动方式在完全可通行区域不存在运行安全问题，当
移动体周围存在障碍栅格时此移动方式可能会发生与障碍物栅格的碰撞问题，为解决此问题加入约束
条件，当在分别与障碍物栅格水平方向和垂直方向的可行栅格两栅格之间通行时，禁止移动体采用对
角式移动方式。

约束条件的加入，实质是改变栅格地图的邻接矩阵，将障碍栅格(数字为“1”的矩阵元素)的对角栅格
设为不可达，即将对角栅格的距离值改为无穷大。其实现MATLAB代码如下：
代码：

%约束移动体在障碍栅格对角运动
%通过优化邻接矩阵实现
%%%%%%%%%%%%%%%%%% 约束移动体移动方式 %%%%%%%%%%%%%%%%%
function W=OPW(map,W)
% map 地图矩阵  % W 邻接矩阵
n = size(map);
num = n(1)*n(2);
for(j=1:n(1))for(z=1:n(2))if(map(j,z)==1)if(j==1)                  %若障碍物在第一行if(z==1)               %若障碍物为第一行的第一个W(j+1,j+n(2)*j)=Inf;W(j+n(2)*j,j+1)=Inf;elseif(z==n(2))         %若障碍物为第一行的最后一个W(n(2)-1,n(2)+n(1)*j)=Inf;W(n(2)+n(1)*j,n(2)-1)=Inf;else                %若障碍物为第一行的其他W(z-1,z+j*n(2))=Inf;W(z+j*n(2),z-1)=Inf;W(z+1,z+j*n(2))=Inf;W(z+j*n(2),z+1)=Inf;endendendif(j==n(1))               %若障碍物在最后一行if(z==1)               %若障碍物为最后一行的第一个W(z+n(2)*(j-2),z+n(2)*(j-1)+1)=Inf;W(z+n(2)*(j-1)+1,z+n(2)*(j-2))=Inf;elseif(z==n(2))            %若障碍物为最后一行的最后一个W(n(1)*n(2)-1,(n(1)-1)*n(2))=Inf;W((n(1)-1)*n(2),n(1)*n(2)-1)=Inf;else                   %若障碍物为最后一行的其他W((j-2)*n(2)+z,(j-1)*n(2)+z-1)=Inf;W((j-1)*n(2)+z-1,(j-2)*n(2)+z)=Inf;W((j-2)*n(2)+z,(j-1)*n(2)+z+1)=Inf;W((j-1)*n(2)+z+1,(j-2)*n(2)+z)=Inf;endendendif(z==1)              if(j~=1&&j~=n(1))       %若障碍物在第一列非边缘位置 W(z+(j-2)*n(2),z+1+(j-1)*n(2))=Inf;W(z+1+(j-1)*n(2),z+(j-2)*n(2))=Inf;W(z+1+(j-1)*n(2),z+j*n(2))=Inf;W(z+j*n(2),z+1+(j-1)*n(2))=Inf;endendif(z==n(2))if(j~=1&&j~=n(1))         %若障碍物在最后一列非边缘位置 W((j+1)*n(2),j*n(2)-1)=Inf;W(j*n(2)-1,(j+1)*n(2))=Inf;W(j*n(2)-1,(j-1)*n(2))=Inf;W((j-1)*n(2),j*n(2)-1)=Inf;endendif(j~=1&&j~=n(1)&&z~=1&&z~=n(2))   %若障碍物在非边缘位置W(z+(j-1)*n(2)-1,z+j*n(2))=Inf;W(z+j*n(2),z+(j-1)*n(2)-1)=Inf;W(z+j*n(2),z+(j-1)*n(2)+1)=Inf;W(z+(j-1)*n(2)+1,z+j*n(2))=Inf;W(z+(j-1)*n(2)-1,z+(j-2)*n(2))=Inf;W(z+(j-2)*n(2),z+(j-1)*n(2)-1)=Inf;W(z+(j-2)*n(2),z+(j-1)*n(2)+1)=Inf;W(z+(j-1)*n(2)+1,z+(j-2)*n(2))=Inf;endendendend
end

2.5 栅格法案例
下面以Djkstra算法为例，其实现如下：

map=[0 0 0 1 0 0 1 0 0 0;1 0 0 0 0 1 1 0 0 0;0 0 1 0 0 0 1 1 0 0;0 0 0 0 0 0 0 0 0 0;0 0 0 0 0 1 0 0 1 0;1 0 0 0 0 1 1 0 0 0;0 0 0 1 0 0 0 0 0 0;1 1 1 0 0 0 1 0 0 0;0 0 0 0 0 1 1 0 0 0;0 0 0 0 0 1 1 0 0 0;];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%建立环境矩阵map%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DrawMap(map); %得到环境地图
W=G2D(map);   %得到环境地图的邻接矩阵
W(W==0)=Inf;  %邻接矩阵数值处理
W=OPW(map,W); %优化邻接矩阵
[distance,path]=dijkstra(W,1,100);%设置起始栅格，得到最短路径距离以及栅格路径
[x,y]=Get_xy(distance,path,map);   %得到栅格相应的x，y坐标
Plot(distance,x,y);   %画出路径

运行结果如下：

其中函数程序：
DrawMap(map) 详见建立栅格地图
W=G2D(map) ；详见建立邻接矩阵
[distance， path] =dijkstra(W， 1， 100) 详见Djk stra算法
[x， y] =Get_xy(distance， path， map) ；
Plot(distance， x， y) ；

二、部分源代码

clc;
close all
clear
load('data4.mat')
S=(S_coo(2)-0.5)*num_shange+(S_coo(1)+0.5);%起点对应的编号
E=(E_coo(2)-0.5)*num_shange+(E_coo(1)+0.5);%终点对应的编号PopSize=20;%种群大小
OldBestFitness=0;%旧的最优适应度值
gen=0;%迭代次数
maxgen =100;%最大迭代次数c1=0.5;%认知系数
c2=0.7;%社会学习系数
c3=0.2;%反向因子
w=0.96;%惯性系数
%%
%初始化路径
w_min=0.5;
w_max=1;
Group=ones(num_point,PopSize);  %种群初始化
flag=1;
%% 初始化粒子群位置
for i=1:PopSizep_lin=randperm(num_point)';%随机生成1*400不重复的行向量%% 将起点编号放在首位index=find(p_lin==S);lin=p_lin(1);p_lin(1)=p_lin(index);p_lin(index)=lin;Group(:,i)=p_lin;%%将每个个体进行合理化处理[Group(:,i),flag]=deal_fun(Group(:,i),num_point,liantong_point,E,num_shange);fangxiang_Group(:,i)=fangxiang(Group(:,i),c3);%方向粒子数量while flag==1%如处理不成功，则初始化个体，重新处理%% 将起点编号放在首位index=find(p_lin==S);lin=p_lin(1);p_lin(1)=p_lin(index);p_lin(index)=lin;Group(:,i)=p_lin;fangxiang_Group(:,i)=p_lin;%%将每个个体进行合理化处理[Group(:,i),flag]=deal_fun(Group(:,i),num_point,liantong_point,E,num_shange);[fangxiang_Group(:,i),flag]=deal_fun(fangxiang_Group(:,i),num_point,liantong_point,E,num_shange);endend%初始化粒子速度（即交换序）
Velocity =zeros(num_point,PopSize);
for i=1:PopSizeVelocity(:,i)=round(rand(1,num_point)'*num_point/10); %round取整
end%计算每个个体对应路径的距离
for i=1:PopSizeEachPathDis(i) = PathDistance(Group(:,i)',E,num_shange);EachPathDis_fangxiang(i) = PathDistance(fangxiang_Group(:,i)',E,num_shange);
endIndivdualBest=Group;%记录各粒子的个体极值点位置,即个体找到的最短路径
IndivdualBestFitness=EachPathDis;%记录最佳适应度值,即个体找到的最短路径的长度
if min(EachPathDis)<min(EachPathDis_fangxiang)
[GlobalBestFitness,index]=min(EachPathDis);%找出全局最优值和相应序号
else[GlobalBestFitness,index]=min(EachPathDis_fangxiang);%找出全局最优值和相应序号
end
%寻优
while gen < maxgenw=w_max-(w_max-w_min)*gen/maxgen;%自适应权重%迭代次数递增gen = gen +1%更新全局极值点位置,这里指路径for i=1:PopSizeif min(EachPathDis)<min(EachPathDis_fangxiang)GlobalBest(:,i) = Group(:,index);elseGlobalBest(:,i) = fangxiang_Group(:,index);endend%求pij-xij ,pgj-xij交换序，并以概率c1，c2的保留交换序pij_xij=GenerateChangeNums(Group,IndivdualBest);  %根据个体最优解求交换序pij_xij=HoldByOdds(pij_xij,c1);%以概率c1保留交换序pgj_xij=GenerateChangeNums(Group,GlobalBest);%根据全局最优解求交换序pgj_xij=HoldByOdds(pgj_xij,c2);%以概率c2保留交换序pfj_xij=GenerateChangeNums(Group,fangxiang_Group);%根据反向求交换序pfj_xij=HoldByOdds(pfj_xij,c3);%以概率c3保留交换序%以概率w保留上一代交换序Velocity=HoldByOdds(Velocity,w);Group = PathExchange(Group,pfj_xij);%根据反向粒子位置进行交换Group = PathExchange(Group,Velocity); %根据交换序进行路径交换Group = PathExchange(Group,pij_xij);%粒子位置变换通过速度、全局性适应度和个体适应度对比来交换来实现，完成自我学习和社会学习Group = PathExchange(Group,pgj_xij);for i = 1:PopSize[Group(:,i),flag]=deal_fun(Group(:,i),num_point,liantong_point,E,num_shange);while flag==1p_lin=randperm(num_point)';index=find(p_lin==S);lin=p_lin(1);p_lin(1)=p_lin(index);p_lin(index)=lin;Group(:,i)=p_lin;[Group(:,i),flag]=deal_fun(Group(:,i),num_point,liantong_point,E,num_shange);endendfor i = 1:PopSize    % 更新各路径总距离EachPathDis(i) = PathDistance(Group(:,i)',E,num_shange);endIsChange = EachPathDis<IndivdualBestFitness;%更新后的距离优于更新前的，记录序号IndivdualBest(:, find(IsChange)) = Group(:, find(IsChange));%更新个体最佳路径IndivdualBestFitness = IndivdualBestFitness.*( ~IsChange) + EachPathDis.*IsChange;%更新个体最佳路径距离[GlobalBestFitness, index] = min(IndivdualBestFitness);%更新全局最佳路径,记录相应的序号if GlobalBestFitness~=OldBestFitness %比较更新前和更新后的适应度值;OldBestFitness=GlobalBestFitness;%不相等时更新适应度值best_route=IndivdualBest(:,index)';endBestFitness(gen) =GlobalBestFitness;%每一代的最优适应度
end
%最优解
index1=find(best_route==E);
route_lin=best_route(1:index1);%最优解
figure(3)
hold on
for i=1:num_shangefor j=1:num_shangeif sign(i,j)==1y=[i-1,i-1,i,i];x=[j-1,j,j,j-1];h=fill(x,y,'k');set(h,'facealpha',0.5)ends=(num2str((i-1)*num_shange+j));text(j-0.95,i-0.5,s,'fontsize',6)end
end
axis([0 num_shange 0 num_shange])%限制图的边界
plot(S_coo(2),S_coo(1), 'p','markersize', 10,'markerfacecolor','b','MarkerEdgeColor', 'm')%画起点
plot(E_coo(2),E_coo(1),'o','markersize', 10,'markerfacecolor','g','MarkerEdgeColor', 'c')%画终点
set(gca,'YDir','reverse');%图像翻转
for i=1:num_shangeplot([0 num_shange],[i-1 i-1],'k-');plot([i i],[0 num_shange],'k-');%画网格线
end
for i=2:index1Q1=[mod(route_lin(i-1)-1,num_shange)+1-0.5,ceil(route_lin(i-1)/num_shange)-0.5];Q2=[mod(route_lin(i)-1,num_shange)+1-0.5,ceil(route_lin(i)/num_shange)-0.5];plot([Q1(1),Q2(1)],[Q1(2),Q2(2)],'r','LineWidth',3)
end
title('粒子群算法-最优路线');%进化曲线
figure(4);
plot(BestFitness);
xlabel('迭代次数')
ylabel('适应度值')
grid on;
title('进化曲线');
disp('粒子群算法-最优路线方案：')
disp(num2str(route_lin))
disp(['起点到终点的距离:',num2str(BestFitness(end))]);
figure(5);
plot(BestFitness*100);
xlabel('迭代次数')
ylabel('适应度值')
grid on;
title('最佳个体适应度值变化趋势');

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.