一、引言
上篇关于使用matlab编写遗传算法求解车辆路径问题写完后，我发现南柯一梦那篇文章的参考文献应该是一篇《中国管理科学》的文章《用混合遗传算法求解物流配送路径优化问题的研究》，果然期刊水平高，文章的质量还是更有保障。这篇文章提出用爬山算法对遗传算法进行改进，本文将对基于爬山算法的遗传算法进行复现，并用于求解原文中的车辆路径优化问题。问题描述在上一篇文章中，这里不再赘述。
二、算法思路及GA代码
2.0算法思路
本文针对VRP问题所编写代码的核心过程为车辆路径的分配过程，即把哪些顾客安排给第几辆车。不同于上一篇文章，本文代码中的种群个体没有包含虚拟配送中心，只包含了需要配送的需求点。所以首先要做的就是为顾客安排车辆，具体的思路如下。以6-9-4-2-1-3-8-7-5为例，首先将第一个顾客6安排给第一辆车（这里的隐含假设是，不存在需求量或配送距离大于单个车辆限额的顾客），接着计算将第二个顾客9加入第一个路径后的配送距离及载重量，如果都没超过，则将9安排给第一辆车，否则将9安排给下一个车辆，以此类推。
对于一个安排好车辆的个体，如果安排的车辆总数小于等于配送中心的车辆限额，则为可行解。如果安排的车辆总数大于配送中心得车辆限额，则对该个体进行惩罚。具体操作就是给该个体的总路径增加一个较大的惩罚值。接着进行的算法步骤为：选择-交叉-变异-爬山-重插入。爬山操作的思想是，对于通过遗传操作形成的每代群体中的最优个体，通过领域搜索实施爬山操作。具体做法是：（1）在个体中随机选择两个基因，交换其位置；（2）判断基因换位后适应值是否增加，若增加，则以换位后的个体取代原个体；（3）重复操作（1）（2）直到达到指定次数为止。
2.1车辆路径分配及个体总路径距离计算
%% LengthInd函数用于计算个体的总路径长度% 输入：个体a,车辆数carNum,客户之间的距离里矩阵DL,% 距离限制disMax,客户重量需求DW,车辆重量限制capMax,客户到需求点距离X% 输出：个体的总路径长度%%% 整体思路是给个体中每个客户赋予对应的车辆编号，从小到大，% 如果新增加的个体超过当前车辆的距离或重量限制，则分配给下一个车辆，% 如果被分配的车辆总数超过了给定车辆数，对个体进行惩罚function lengthInd = LengthInd(a,carNum,DL,disMax,DW,capMax,X)custNum = length(a); % 记录顾客的数量matchA = zeros(custNum,4); % 记录顾客和车辆的匹配情况，第一列为顾客的编号，第二列为匹配的车辆编号% 第三列为子路径总距离只记录在一个路径最后一个顾客对应的位置，第四列为需求重量matchA(:,1) = a';matchA(:,4) = DW(a)';matchA(1,2) = 1; % 第一个顾客一定在第一个路径中k = 1; % 车辆的编号for i = 2:custNum% 对于第i个顾客，首先要计算的是在当前路径加入该顾客后的路径总距离以及总重量custTemp = find(matchA(:,2)==k); % 先找到第k个路径中包含的顾客custK = [matchA(custTemp,1)' matchA(i,1)]; % 在第k个路径中加入第i个顾客后的总顾客情况i1 = custK(1:end-1);i2 = custK(2:end);custKL = sum(DL((i1-1)*custNum+i2));custKL = custKL + X(custK(1)) + X(custK(end)); % 计算第k辆车的总路径custKW = sum(DW(custK)); % 计算第k辆车的总载重if (custKL > disMax || custKW > capMax)k = k + 1;matchA(i,2) = k;% 此时判断出了上一辆车的最后一个需求点，因此需要记录上一个路径的距离custK(end) = [];i1 = custK(1:end-1);i2 = custK(2:end);custKL = sum(DL((i1-1)*custNum+i2));custKL = custKL + X(custK(1)) + X(custK(end));matchA(i-1,3) = custKL;elsematchA(i,2) = k;endif i == custNumi1 = custK(1:end-1);i2 = custK(2:end);custKL = sum(DL((i1-1)*custNum+i2));custKL = custKL + X(custK(1)) + X(custK(end));matchA(i-1,3) = custKL;endendaLength = sum(matchA(:,3)); % 计算a的总路径if matchA(end,2) > carNumlengthInd = aLength + 1000; % 车辆数超过限制，给予惩罚elselengthInd = aLength;end
2.2选择
本文中选择操作采用的是轮盘赌方式，具体代码如下：
%% Select函数采用轮盘赌方式对种群中的个体进行选择% 输入：种群population,车辆数carNum,客户之间的距离里矩阵DL,% 距离限制disMax,客户重量需求DW,车辆重量限制capMax,客户到需求点的距离X% 输出：选择后的种群newPop%%function newPop = Select(population,carNum,DL,disMax,DW,capMax,X)% 首先对种群的适应度进行计算popNum = size(population,1); % 记录种群的个数fitInd = zeros(popNum,1);% 记录个体的适应度newPop = population;% 记录选择后的个体for i = 1:popNumfitInd(i) = 1/LengthInd(population(i,:),carNum,DL,disMax,DW,capMax,X);endprobInd = fitInd./sum(fitInd); % 每个个体被选择的概率probIndAcu = cumsum(probInd);% probIndAcu是probInd的累加概率，这一步为轮盘构建for i = 1:popNumnewPop(i) = population(popNum - length(find(rand<=probIndAcu)) + 1); % 根据轮盘赌概率选择个体end
2.3交叉
首先是配对函数
function [new_pop_intercross]=Mating_pool(population_num,population,Pc)%%%输入:population,population_num,Pc%输出：1.new_popopulation_intercross%2.c3，配对池：随机将种群population两两配对%3.pool%%pl=randperm(population_num); % 将种群打散用于配对num=population_num/2;c3=zeros(2,num);pool=[];new_pop_intercross=population;for kj=1:numc3(1,kj)=pl(2*kj-1);c3(2,kj)=pl(2*kj);end %生成"配对池c3 %%判断“配对池c3”每一对个体的随机数是否小于交叉概率Pcrd=rand(1,num);for kj=1:numif rd(kj)<Pcpool=[pool,c3(:,kj)];endend%%判断配对池每一对个体的随机数是否小于交叉概率Pc,若小于，保存到“产子池pool”pool_num=size(pool,2);for kj=1:pool_numc1=population(pool(1,kj),:); % 产子池第kj列的第一个个体c2=population(pool(2,kj),:); % 产子池第kj列的第二个个体[new_c1,new_c2]=cross(c1,c2); % 对两个个体进行交叉操作new_pop_intercross(pool(1,kj),:)=new_c1;new_pop_intercross(pool(2,kj),:)=new_c2;endend
其次是两个体的交叉函数
% 因为VRP问题中的一个重要约束是要遍历需求点，且每个点只访问一次，所以A，B两个染色体配对成功后% 首先要确定交叉的基因段，选择好后先对染色体A进行交叉操作，清楚A中交叉基因段，在B中找到A中被% 清除的需求点的位置，去除B中其余的需求点，将剩余的需求点直接插入A中的交叉基因段，这样就保证了% 染色体A仍然遍历所有的需求点且每个需求点只经过一次function [A,B]=cross(A,B)A_1=A;B_1=B;r=randperm(length(A));c=min(r(1,1:2));d=max(r(1,1:2));for i=c:dA(i)=0;endB_1(ismember(B_1,A))=[];A(1,c:d)=B_1;for i=c:dB(i)=0;endA_1(ismember(A_1,B))=[];B(1,c:d)=A_1;end
2.4变异
function [Mut_Pop]=Mutation(Cross_Pop,Pm)Mut_Pop=Cross_Pop;Cross_Pop_num=size(Cross_Pop,1);for j=1:Cross_Pop_numA=Cross_Pop(j,:); % 要进行变异的染色体A_1=A; % 对A进行记录，之后交换变异需求点的时候要用到n=size(A,2);r=rand(1,n);Pe=find(r<Pm); % 这句代码意思是染色体A中的每个基因都可能变异，返回一个变异基因位置的数组sum_Pe=size(Pe,2);for i=1:sum_Pec=A(Pe(i));A_1(Pe(i))=A_1(find(r==max(r))); % 找到r中最大值对应的位置，这个基因大概率不会变异A_1(find(r==max(r)))=c;% 接着两个基因进行交换Mut_Pop=[Mut_Pop;A_1]; % 这里是把变异后的个体直接加在交叉后的种群后面endend
2.5爬山
%% ClimbMoun函数对变异后的最优个体进行爬山寻优操作% 输入：个体a,车辆数carNum,客户之间的距离里矩阵DL,% 距离限制disMax,客户重量需求DW,车辆重量限制capMax,客户到需求点距离X% 输出：爬山寻优后的个体function aClimb = ClimbMoun(a,carNum,DL,disMax,DW,capMax,X)aNum = length(a);% 记录a的长度lengthOrigina = LengthInd(a,carNum,DL,disMax,DW,capMax,X);for i = 1:100 % 爬山寻优次数temp = a;r1 = randsrc(1,1,1:aNum);r2 = randsrc(1,1,1:aNum); % 在1-aNum中找出两个随机数while r1 ==r2r2 = randsrc(1,1,1:aNum);endtempr1 = temp(r1);temp(r1) = temp(r2);temp(r2) = tempr1;lengthClimba = LengthInd(temp,carNum,DL,disMax,DW,capMax,X);if lengthOrigina > lengthClimba % 寻得更优个体aClimb = temp;break;endendif lengthOrigina <= lengthClimba% 未寻得更优个体aClimb = a;end
2.6主函数
clcclearclose all;tic;%% 算法参数popNum=80; % 种群规模Max_gen=300; % 迭代次数Pc=0.8;% 交叉概率Pm=0.05;% 变异概率%% 问题参数carNum=5; % 车辆最大数量custNum=20; % 顾客人数capMax=8; % 每辆车的载重限制disMax=50; % 每辆车的总路径限制trace = zeros(2,Max_gen); % 用于记录总路径的进化过程load data_DW% 加载需求点的位置和需求量数据 %%% 种群初始化，并计算需求点之间的距离矩阵DL，需求量DW，配送中心到需求点距离XDL = zeros(custNum,custNum); % 建立需求点之间的距离矩阵X = zeros(1,custNum); % 建立需求点到配送中心的距离for i = 1:custNumfor j = i:custNumDL(i,j) = ((data_DW(i,1)-data_DW(j,1))^2+(data_DW(i,2)-data_DW(j,2))^2)^0.5;DL(j,i) = DL(i,j);endX(i) = ((data_DW(i,1)-data_DW(custNum + 1,1))^2+(data_DW(i,2)-data_DW(custNum + 1,2))^2)^0.5;endDW = data_DW(:,3)'; % 建立需求点的重量population=zeros(popNum,custNum);for i=1:popNumpopulation(i,:)=randperm(custNum); % 产生含有虚拟配送中心的随机个体end%% 基于爬山算法的遗传算法y=1;%循环计数器Total_Dis = zeros(popNum,1);for i = 1:popNumTotal_Dis(i) = LengthInd(population(i,:),carNum,DL,disMax,DW,capMax,X);endtrace(2,1) = min(Total_Dis);while y<Max_gen% 选择，根据个体的适应度函数，采用轮盘赌的方式对个体进行选择popSelect = Select(population,carNum,DL,disMax,DW,capMax,X);%交叉 交叉后种群的数量不变[new_pop_intercross]=Mating_pool(popNum,population,Pc);%变异 变异后的种群数量会增加[new_pop_mutation]=Mutation(new_pop_intercross,Pm);% 计算变异后种群中个体的总距离mutation_num=size(new_pop_mutation,1);Total_Dis=zeros(mutation_num,1);for i = 1:mutation_numTotal_Dis(i) = LengthInd(new_pop_mutation(i,:),carNum,DL,disMax,DW,capMax,X);endminDisInd = find(Total_Dis==min(Total_Dis)); % 找到距离最短的个体% 对最优个体进行爬山操作for i = 1:length(minDisInd)new_pop_mutation(minDisInd(i),:) = ClimbMoun(new_pop_mutation(minDisInd(i),:),carNum,DL,disMax,DW,capMax,X);end%更新种群 根据目标函数保留变异后的前population_num个染色体new_pop_new=zeros(popNum,custNum);[~, index] = sort(Total_Dis);for k=1:popNumnew_pop_new(k,:)=new_pop_mutation(index(k),:);endpopulation=new_pop_new;%迭代次数加一y=y+1;trace(1,y) = y;trace(2,y) = min(Total_Dis);endX_Best=new_pop_mutation(minDisInd(1),:);customer_car=CarAllocate(new_pop_mutation(minDisInd(1),:),DL,disMax,DW,capMax,X);Y_Obj=Total_Dis(minDisInd(1),1);t=toc;
2.7得出车辆分配结果
这个函数和路径计算函数差不多，但我对函数的运用还不太熟练，所以又增加了一个函数对其进行求解
%% CarAllocate函数用于得到个体的车辆安排% 输入：个体a,车辆数carNum,客户之间的距离里矩阵DL,% 距离限制disMax,客户重量需求DW,车辆重量限制capMax,客户到需求点距离X% 输出：个体的车辆安排%%% 整体思路是给个体中每个客户赋予对应的车辆编号，从小到大，% 如果新增加的个体超过当前车辆的距离或重量限制，则分配给下一个车辆，% 如果被分配的车辆总数超过了给定车辆数，对个体进行惩罚function carAllocate = CarAllocate(a,DL,disMax,DW,capMax,X)custNum = length(a); % 记录顾客的数量matchA = zeros(custNum,4); % 记录顾客和车辆的匹配情况，第一列为顾客的编号，第二列为匹配的车辆编号% 第三列为子路径总距离只记录在一个路径最后一个顾客对应的位置，第四列为需求重量matchA(:,1) = a';matchA(:,4) = DW(a)';matchA(1,2) = 1; % 第一个顾客一定在第一个路径中k = 1; % 车辆的编号for i = 2:custNum% 对于第i个顾客，首先要计算的是在当前路径加入该顾客后的路径总距离以及总重量custTemp = find(matchA(:,2)==k); % 先找到第k个路径中包含的顾客custK = [matchA(custTemp,1)' matchA(i,1)]; % 在第k个路径中加入第i个顾客后的总顾客情况i1 = custK(1:end-1);i2 = custK(2:end);custKL = sum(DL((i1-1)*custNum+i2));custKL = custKL + X(custK(1)) + X(custK(end)); % 计算第k辆车的总路径custKW = sum(DW(custK)); % 计算第k辆车的总载重if (custKL > disMax || custKW > capMax)k = k + 1;matchA(i,2) = k;% 此时判断出了上一辆车的最后一个需求点，因此需要记录上一个路径的距离custK(end) = [];i1 = custK(1:end-1);i2 = custK(2:end);custKL = sum(DL((i1-1)*custNum+i2));custKL = custKL + X(custK(1)) + X(custK(end));matchA(i-1,3) = custKL;elsematchA(i,2) = k;endendcarAllocate = matchA(:,2)';
三、算例展示
这里的算例选用的是文章中的实例2，首先将需要载入的数据给出，及配送中心和需求点的坐标及需求量
12.80000000000008.500000000000000.10000000000000018.40000000000003.400000000000000.40000000000000015.400000000000016.60000000000001.2000000000000018.900000000000015.20000000000001.5000000000000015.500000000000011.60000000000000.8000000000000003.9000000000000010.60000000000001.3000000000000010.60000000000007.600000000000001.700000000000008.600000000000008.400000000000000.60000000000000012.50000000000002.100000000000001.2000000000000013.80000000000005.200000000000000.4000000000000006.7000000000000016.90000000000000.90000000000000014.80000000000002.600000000000001.300000000000001.800000000000008.700000000000001.3000000000000017.1000000000000111.900000000000007.4000000000000011.700000000000000.2000000000000002.800000000000001.1000000000000011.900000000000019.80000000000001.5000000000000013.200000000000015.10000000000001.600000000000006.400000000000005.600000000000001.700000000000009.6000000000000014.80000000000001.5000000000000014.5000000000000130
多次实验得到的最优解为98.3Km，具体的路径安排如下， 0-14-4-3-17-18-0, 0-20-11-6-13-19-8-0, 0-7-15-9-12-2-10-1-5, 0-4-0，进化情况如图所示。10次实验所得到的平均最短路径长度为109.27Km，结果均由于原文的结果。
参考文献
[1] 郎茂祥,胡思继.用混合遗传算法求解物流配送路径优化问题的研究[J].中国管理科学,2002(05):52-57.DOI:10.16381/ki.issn1003-207x.2002.05.011.
[2] 遗传算法求解车辆路径问题 - 知乎
