SDPTWVRP基于matlab头脑风暴算法求解带时间窗和同时取送货车辆路径问题含Matlab源码 1990期

Posted 2022-07-24 佛怒火莲

一、头脑风暴优化算法（BSO）简介

头脑风暴(BSO)是一种将搜索空间不断缩减的算法。此算法通过不断迭代,最终将局部最优解慢慢精确至全局最优解。其基本过程如图1所示。

图1 头脑风暴算法流程

1 种群初始化
随机产生n个个体,每个个体代表一种配送方案。这里应用Matlab仿真软件进行仿真,故个体采用二维编码的方式进行表示。例如,若用数字0表示配送中心,数字1-10表示顾客,数字11-13表示充电站,则个体随机产生的过程可表示如下:首先,得到顾客经过顺序,如[1 2 3 4 5 6 7 8 9 10];其次,鉴于顾客货物需求以及电动货车容量进行分组,若顾客需求皆为3,电动货车容量为10,则上述顺序可更新为[0 1 2 3 0 4 5 6 0 7 8 9 0 10 0],每两个数字0之间的数字即为一辆电动货车的配送路线,此顺序即采用了4辆电动货车;再次,根据上述顺序中每辆电动货车的配送顺序求解耗电量、剩余电量,根据剩余电量是否能支撑电动货车到达下一个顾客点最近充电站为判定依据进行充电行为,充电量约束满足式(11),最高电量为充满,最低电量为刚好到达下一充电站或配送中心。例如[0 1 2 11 3 0 4 5 6 0 7 12 8 9 0 10 0],此顺序中就表示第一辆与第三辆电动货车需要在配送过程中充电,而此顺序即为个体表示。最后,将个体存放至二维数组中,每一行代表一个个体。

2 聚 类
将n个个体进行分类,聚成m个类,分别计算每个个体的适应度值并选出每一类的最优个体作为聚类中心。此适应度值为耗电量与顾客满意度倒数之和,适应度值最小者即为最优个体。耗电量与顾客满意度的求解方式见2.1节、2.2节。

3 产生新个体
如公式所示:XNEW=XSELECT+δ(μ,σ),其中XNEW为新个体,XSELECT为选择的个体,δ(μ,σ)为以μ为均值、σ为方差的高斯随机函数。经过一系列的变异、交叉操作,产生新个体。其中,用于新个体产生公式中的个体选择操作参照图1中的选择过程。通过选择条件P1、P2、P3、P4、P5来决定选择对象,然后根据选择的对象进行交叉、变异操作。该研究中的交叉操作采用两点交叉、变异也采用两点变异。两点交叉即为通过选择操作P5选择的两个对象作为交叉父本,通过设置两个交叉点来进行交叉互换,然后通过映射消除重复位置的因素,最终得到交叉后的个体。例如,现有两个选择对象:[1 2 3 4 5 6]和[6 5 4 3 2 1],交叉点为第3、第5个位置,那么保持交叉点中间位置不变其余位置互换得到:[6 5 3 4 5 1]和[1 2 4 3 2 6],经过映射消除得到:[6 2 3 4 5 1]和[1 5 4 3 2 6],即为交叉后结果。而变异操作则通过对选择操作P4选择的对象进行两点调换。例如,现有一选择对象:[1 2 3 4 5 6],变异点同样为第3、第5个位置,经过变异操作后为:[1 2 5 4 3 6]。

4 更新个体
将新个体与原个体相比较,保留适应度值更好的个体。将每一代的最优个体依次存在二维数组中,每一行代表此迭代的最优结果。

5 输出最优解
当达到最大迭代次数后,输出适应度值最优的个体。

二、部分源代码

%% 用xlsread函数来读取xlsx文件
data=xlsread(实例验证数据.xlsx,转换后数据,A2:H17);
cap=150;                            %车辆最大装载量
v=30/60;                            %车辆行驶速度=30km/h=30/60km/min
%% 提取数据信息
E=data(1,6);                        %配送中心时间窗开始时间
L=data(1,7);                        %配送中心时间窗结束时间
vertexs=data(:,2:3);                %所有点的坐标x和y
customer=vertexs(2:end,:);          %顾客坐标
cusnum=size(customer,1);            %顾客数
v_num=10;                           %车辆最大允许使用数目
demands=data(2:end,4);              %需求量
pdemands=data(2:end,5);             %回收量
a=data(2:end,6);                    %顾客时间窗开始时间[a[i],b[i]]
b=data(2:end,7);                    %顾客时间窗结束时间[a[i],b[i]]
s=data(2:end,8);                    %客户点的服务时间
h=pdist(vertexs);
dist=squareform(h);                 %距离矩阵
N=cusnum+v_num-1;                   %解长度=顾客数目+车辆最多使用数目-1
%% 参数初始化
alpha=10;                           %违反的容量约束的惩罚函数系数
belta=100;                          %违反时间窗约束的惩罚函数系数
MAXGEN=150;                         %最大迭代次数
NIND=50;                            %种群数目
cluster_num=5;                      %聚类数目
p_replace=0.1;                      %用随机解替换一个聚类中心的概率
p_one=0.5;                          %选择1个聚类的概率
p_two=1-p_one;                      %选择2个聚类的概率，p_two=1-p_one
p_one_center=0.3;                   %选择1个聚类中聚类中心的概率
p_two_center=0.2;                   %选择2个聚类中聚类中心的概率
%% 种群初始化
Population=InitPop(NIND,N);
%% 主循环
gen=1;                              %计数器初始化
bestInd=Population(1,:);            %初始化全局最优个体
bestObj=ObjFunction(bestInd,v_num,cusnum,cap,demands,pdemands,a,b,s,L,dist,v,alpha,belta);  %初始全局最优个体的目标函数值
BestPop=zeros(MAXGEN,N);            %记录每次迭代过程中全局最优个体
BestObj=zeros(MAXGEN,1);            %记录每次迭代过程中全局最优个体的目标函数值 
BestTD=zeros(MAXGEN,1);            %记录每次迭代过程中全局最优个体的总距离
while gen<=MAXGEN
    %% 计算目标函数值
    Obj=ObjFunction(Population,v_num,cusnum,cap,demands,pdemands,a,b,s,L,dist,v,alpha,belta);
    %% K-means聚类
    Idx=kmeans(Obj,cluster_num,Distance,cityblock,Replicates,2);
    cluster=cell(cluster_num,2);                    %将解储存在每一个聚类中
    order_cluster=cell(cluster_num,2);              %将储存在每一个聚类中的个体按照目标函数值排序
    for i=1:cluster_num
        clusteri,1=Population(Idx==i,:);          %将个体按照所处的聚类编号储存到对应的聚类中
        cluster_row(i)=size(clusteri,1,1);           %计算当前聚类中个体数目
        for j=1:cluster_row(i)
            Individual=clusteri,1(j,:);           %当前聚类中第j个个体
            %计算当前聚类中第j个个体的目标函数值
            clusteri,2(j,1)=ObjFunction(Individual,v_num,cusnum,cap,demands,pdemands,a,b,s,L,dist,v,alpha,belta);
        end
        [order_clusteri,2,order_index]=sort(clusteri,2) ; %将当前聚类中的所有个体按照目标函数值从小到大的顺序进行排序
        order_clusteri,1=clusteri,1(order_index,:); %将当前聚类中的所有个体按照排序结果重新排列
        order_index=0;                                  %重置排序序号
    end
    cluster_fit=cell2mat(order_cluster);                %将聚类的元胞数组转换为矩阵，最后一列为个体的目标函数值
    %% 以一定的概率随机从m个聚类中心中选择出一个聚类中心，并用一个新产生的随机解更新这个被选中的聚类中心
    order_cluster=replace_center(p_replace,cluster_num,N,order_cluster,...
    v_num,cusnum,cap,demands,pdemands,a,b,s,L,dist,v,alpha,belta);
    %% 更新这n个个体
    Population=update_Population(Population,cluster_num,cluster_row,Idx,order_cluster,p_one,p_one_center,p_two_center...
    ,v_num,cusnum,cap,demands,pdemands,a,b,s,L,dist,v,alpha,belta);
    %% 计算原始Population的目标函数值
    Obj=ObjFunction(Population,v_num,cusnum,cap,demands,pdemands,a,b,s,L,dist,v,alpha,belta);
    %% 从更新后的Population中选择目标函数值在前50%的个体
    offspring=select(Population,v_num,cusnum,cap,demands,pdemands,a,b,s,L,dist,v,alpha,belta);
    %% 对选择出的个体进行局部搜索操作
    offspring=LocalSearch(offspring,v_num,cusnum,a,b,s,L,dist,demands,pdemands,cap,alpha,belta,v);
    %% 将局部搜索后的offspring与原来的Population进行合并
    Population=merge(Population,offspring,Obj);
    %% 计算合并后的Population的目标函数值
    mObj=ObjFunction(Population,v_num,cusnum,cap,demands,pdemands,a,b,s,L,dist,v,alpha,belta);
    %% 找出合并后的Population中的最优个体
    [min_len,min_index]=min(mObj);               %当前种群中最优个体以及所对应的序号
    %如果当前迭代最优个体目标函数值小于全局最优目标函数值，则更新全局最优个体
    if min_len<bestObj
        bestObj=min_len;
        bestInd=Population(min_index,:);
    end
    %% 打印各代最优解
    disp([第,num2str(gen),代最优解:])
    [bestVC,bestNV,bestTD,best_vionum,best_viocus]=decode(bestInd,v_num,cusnum,cap,demands,pdemands,a,b,s,L,dist,v);
     disp([目标函数值：,num2str(bestObj),，车辆使用数目：,num2str(bestNV),，车辆行驶总距离：,num2str(bestTD),，违反约束路径数目：,num2str(best_vionum),，违反约束顾客数目：,num2str(best_viocus)]);
    fprintf(\\n)
    %% 距离全局最优个体
    BestObj(gen,1)=bestObj;                     %记录全局最优个体的目标函数值 
    BestPop(gen,:)=bestInd;                     %记录全局最优个体
    BestTD(gen,1)=bestTD;                       %记录全局最优个体的总距离 
    %% 计数器加1
    gen=gen+1;
end
%% 画出最优配送路线图
draw_Best(bestVC,vertexs);
%% 打印每次迭代的全局最优个体的总距离变化趋势图
figure;
plot(BestTD,LineWidth,1);
title(优化过程)
xlabel(迭代次数);
ylabel(配送方案行驶总距离);
toc  
%% 将当前个体解码为配送方案
%输入Individual：  当前个体
%输入v_num：       车辆最大允许使用数目
%输入cusnum：      顾客数目
%输入cap：         货车最大装载量
%输入demands：     顾客需求量
%输入pdemands：    顾客回收量
%输入a：           顾客时间窗开始时间[a[i],b[i]]
%输入b：           顾客时间窗结束时间[a[i],b[i]]
%输入L：           配送中心时间窗结束时间
%输入s：           客户点的服务时间
%输入dist：        距离矩阵
%输入v：           车辆行驶速度
%输出VC：          配送方案，即每辆车所经过的顾客
%输出NV：          车辆使用数目
%输出TD：          车辆行驶总距离
%输出violate_num： 违反约束路径数目
%输出violate_cus： 违反约束顾客数目
function [VC,NV,TD,violate_num,violate_cus]=decode(Individual,v_num,cusnum,cap,demands,pdemands,a,b,s,L,dist,v)
violate_num=0;                                      %违反约束路径数目
violate_cus=0;                                      %违反约束顾客数目
VC=cell(v_num,1);                                   %每辆车所经过的顾客
count=1;                                            %车辆计数器，表示当前车辆使用数目
location0=find(Individual>cusnum);                  %找出个体中配送中心的位置
for i=1:length(location0)
    if i==1                                              %第1个配送中心的位置
        route=Individual(1:location0(i));                %提取两个配送中心之间的路径
        route(route==Individual(location0(i)))=[];       %删除路径中配送中心序号
    else
        route=Individual(location0(i-1):location0(i));   %提取两个配送中心之间的路径
        route(route==Individual(location0(i-1)))=[];     %删除路径中配送中心序号
        route(route==Individual(location0(i)))=[];       %删除路径中配送中心序号
    end
    VCcount=route;                                %更新配送方案
    count=count+1;                                  %车辆使用数目
end
route=Individual(location0(end):end);               %最后一条路径       
route(route==Individual(location0(end)))=[];        %删除路径中配送中心序号
VCcount=route;                                    %更新配送方案
[VC,NV]=deal_vehicles_customer(VC);                 %将VC中空的数组移除
for j=1:NV
    route=VCj;
    %判断一条配送路线上的各个点是否都满足装载量约束和时间窗约束，1表示满足，0表示不满足
    flag=JudgeRoute(route,demands,pdemands,cap,a,b,s,L,dist,v); 
    if flag==0
        violate_cus=violate_cus+length(route);      %如果这条路径不满足约束，则违反约束顾客数目加该条路径顾客数目
        violate_num=violate_num+1;                  %如果这条路径不满足约束，则违反约束路径数目加1
    end
end
TD=travel_distance(VC,dist);                        %该方案车辆行驶总距离
end

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 齐元豪,王凯,付亚平.基于头脑风暴算法的电动货车路径优化问题[J].计算机技术与发展. 2020,30(04)

3 备注
简介此部分摘自互联网，仅供参考，若侵权，联系删除