【TWVRP】基于matlab智能水滴算法带时间窗的车辆路径与调度优化问题【含Matlab源码 1228期】

一、智能水滴算法简介

1 概述
在大自然中，河道是由流动的水滴组成的巨大群体移动创造的，相对于水滴群体本身来说，它们流经的河道周围就是它们所处在的环境，而这个环境正在被水滴群体的流动方向剧烈的改变着，同时周围的环境也在影响着水滴群体的流向。例如，环境当中的坚便的巨石泥土阳碍水滴的流动，而柔软松弛的泥土对于河水的流动的阻力就要小得多，实际上，大自然中的河流样貌就是由丁水滴群体和周围环境长期相互作用的结果而形成的。
依据人们的常识，大自然当中大部分的河流水道都是曲曲折折的，然而河流中的水清并不是像人类一样依靠视觉才到湖中或是海里的，而是，水滴在地球重力的作用下不断向前流动的。我们都知道，地球重力的方向指向地心，所以如果周围环境当中没有障碍物的话，水滴应该会沿着重力的方向向地心运动：由于在地球重力的作用下的，水滴在流动的过程当中会获得一定的加速度，因此随着运动其速度也会慢慢地不断增大。但是，在现实情况下，并不存在这样理想的状态也不存在从水清源头到地心的理想直线路径，甚至由于障砖物的影响，水消运动的实际路径与理想路径是远远不同的，因此才会有今天我们看到的河流样貌.而经研究发现，水滴在移动情况下建立起米的河道是最优的。河流当中的水滴具有两种属性：其本身具有一定的速度和携带一定量的泥土。在这样的情况下，水滴在运动时，就能够将泥上从一个地方向搬运到另一个地方。据我们所知，泥土会被水流从速度较高的地方携带到水流速度较低的地方，原因是当水滴速度较高时其具有的动能较大，能够克服地球重力对其所抛带泥土的作用，当水流具有的速度较低时，水滴携带的泥土会在地球重力的作用在水滴行进的路径中沉积下来。由于这个原因，泥土就会在水流速度较高的区域被搬移的越来越多，河道相应变得越来越深，也就吸引了越来越多的水流流向此处，进而可以容纳更多的水，同时水滴携带的泥土会在水流速度变慢时由于地球重力的作用
卸下来并沉积在河床上。
当一个水滴在河流中从一个地方流动另一个地方时，水滴具有的速度、水滴携带的泥土量、两个地方之间的泥土量三者都会发生明显的变化：
1)水滴流动时具有的速度会越来越大；
2)水滴携带的泥土量会越来越多：
3)位于这两点之间的河床上的泥土量会越来越少。

事实上，正如前文所描述，水清在从一个地方移动另一个地方的过程中会从河床上携带走一定量的泥土，河道上的泥土量会越来越少相应地水滴木身携带的泥土量会越来越多。在这个移动的过程中，水滴流动时具有的速度也会不断的变化。速度较大的水滴在流动时会从河床上携带走更多的泥上，除此之外，由于水滴流动时具有速度的快慢还和它周用所处的环境有关，周围的环境对水滴的移动起着阻碍的作用，所以当周围的泥土量较大时，水清在此段路径上流动时速度增长的就会越慢，当周围环境的泥土最较少时，水滴在此段路径上获得的速度就相应的会越来越快，向前流动就会更加容易：内此我们得到自然界中移动时的水滴具备以下属性：
1)流动速度高的水滴比流动速度低的水滴在路径中移动时能够携带走更多的泥土：
2)水滴能够在泥土量较小的路径中获取的速度增量大于比在泥土最较多的路径中获取的速度增量：
3)水滴在选择路径流动时会以较大的可能性朝着泥土量较少的路径移动.
3 智能水滴算法的基本原理
前文已经介绍了大自然中水滴具有某些典型的特点，根据上文中描述的这些特性，根据人们的想象，可以虚拟建立人工模型，称为智能水滴(IWD) ，它具备两个重要的属性：我们一般用soil(IWD) 表示水滴木身携带的泥土量，用veloci tv(/WD) 表示水清向前流动时本身具有的速度口，在工程问题巾，智能水滴流动时周围的环境就代表着我们需要解决的工程问题，而由于相互作用智能水滴寻找的最短路径就是我们解决工程问题找到的最优路径，在子找最优路径的过程中，智能水滴具有的两项重要属性即soil(JWD) 和
velocity(IWD) 是不断变化的。
我们假设某些智能水滴在某种环境下流动，从始点流向某个终点，在这种情况下，从始点到达终点的路径可能不止一条.这里所指的环境就是实践当中我们需要解决的工程类问题，可以分为两种问题：
其一：在终点已知的情况下，解决此类问题需要按照某个标准(一般情况下是距离)从始点到终点找到最佳路径.
其二：在终点的位置事先不知道的情况下，解决此类问题需要按照代价或者其他标准找到终点最优的位置。

二、部分源代码

clear all;
clc;
format short g;
global D;
global q;
global q1;
global ss;
global E;
global L;
global ELL;
% test=xlsread('test.xlsx',1);
%
% position=test(:,2:3);% data=[
%     100 100 0   0    0   40
%     60 100 10 0.2  0   3
%     20 110 15 0.3  0 4
%     160 150 21 0.3  1 6
%     160 110 16 0.4  0 4
%     30 20 25 0.3  1 4
%     100 60 22 0.2  0 6
%     100 170 15 0.1  0 5
%     30 10 12 0.4  0 6
%     60 50 15 0.3  1 7
%     160 160 20 0.2  0 6
%     50 140 23 0.3  0 6];
position=[18.70 15.2916.47 8.4520.07 10.1419.39 13.3725.27 14.2422.00 10.0425.47 17.0215.79 15.1016.60 12.3814.05 18.1217.53 17.3823.52 13.4519.41 18.1322.11 12.5111.25 11.0414.17 9.7624.00 19.8912.21 14.50];
% 配送中心及各个需求点之间的距离矩阵D
D=squareform(pdist(position(:,1:2),'euclidean'));
% 车辆的单位行驶成本: h
% 动用每辆车的固定成本：R.
% 车辆的行驶速度v;
% 车辆的最大装载量 Qmax。
h=1;R=10;v=1;Qmax=200;epsilon=0.001;
R0=2;R1=0.8;
% 车辆在配送中心及各个需求点之间行驶的时间矩阵 T
T=D/v;% 最大迭代次数
Iter=30;
% 水滴个数N
N=size(position,1);
% 速度更新参数：
av=1;bv=0.01;cv=1;
% 泥土量更新参数：
as=1;bs=0.01;cs=1;
% 局部泥土更新权系数
alpha=0.9;
% 全局泥土更新权系数
beta=0.9;
% 任意两点间的初始泥土量 Initsoil
Initsoil=1000;
% 初始泥土量矩阵W ???
W=ones(N)*1000;
% 每个水滴的初始速度InitVel；
% InitVel=100;
InitVel=randperm(100,1);
% 全局最优目标函数值
TotalZ=1000000;
% 全局最优路径
TotalRoute=[];
% 主程序
t=0;
WaterDrop(1,N)=struct('SW',[],... % 水滴对应的初始泥土量矩阵'Source',[],... % 水滴的出发点为1; 即配送中心'Target',[],... %'VisitNode',[],... % 已访问过的点序列 (访问路径)'UnvisitNode',[],... % 未访问过的点的集合'Vel',[],... % 水滴的初始速度'Soil',[],... % 水滴携带的初始泥土量'Q',[],... % 水滴出发时的装载量'S',[],... % 水滴到达source(k)点的时刻'ZZ',[],... %水滴对应的目标函数初始值'FK',[]); % 从source(k)出发可以去的下一个需求点的集合 FV
%  q=test(:,4); % 每个需求点的需求量
% q1=TestData(:,5);
q=[0 6 5 11 6 3 8 5 6 4 5 7 6 10 9 4 7 8 ]';
% q1=[0 3.6 2  4.6 3.6 2.4 4.8 3 3.6 2.4 3 4.2 3.6 4 5.4 2.6 2.4 3]';
%  ss=test(:,7); % 每个需求点的服务时间
ss=[0  1.8  1.0  2.3  1.8  1.2  2.4  1.5  1.8  1.2  1.5  2.1  1.8  2.0  2.7  1.3  1.2  1.5]';
E=[0  5.0  4  1  2.0  5  2  1 3  1  2  2.0  2  3  2  3  2  1]';
L=[40  20  15  20  20  15  18  24  27  20  16  20  10  25  28  24.0  20  23]';
%
% E=test(:,5); % 每个需求点的时间窗下限
%  L=test(:,6); % 每个需求点的时间窗上限
ELL=L-E;
% EP=test(:,8);
% LP=test(:,9);
I=1:N;
while t < Iterclcfprintf('第%d次进化\n',t+1) ;% 设置初始动态变量：% 本次迭代的最优目标函数值Z=1000000;% 本次迭代的最优路径Route=[];for k=1:NWaterDrop(k).SW=W;WaterDrop(k).Source=1;WaterDrop(k).VisitNode=[1];WaterDrop(k).UnvisitNode=I;  WaterDrop(k).Vel= InitVel;% 水滴携带的初始泥土量WaterDrop(k).Soil=0;% 水滴出发时的装载量WaterDrop(k).Q(WaterDrop(k).Source)=0;% 水滴到达source(k)点的时刻；WaterDrop(k).S(WaterDrop(k).Source)=0;% 第k个水滴对应的目标函数初始值WaterDrop(k).ZZ=0;        endfor k=1:Nwhile WaterDrop(k).Source~=1 || ~isequal(WaterDrop(k).UnvisitNode,[1])WaterDrop(k).FK=[];if WaterDrop(k).Source==1alt=2;elsealt=1;end             for i=alt:length(WaterDrop(k).UnvisitNode)WaterDrop(k).Q(WaterDrop(k).UnvisitNode(i))=WaterDrop(k).Q(WaterDrop(k).Source)+q(WaterDrop(k).UnvisitNode(i));% 判断i点的载重量是否小于车辆最大载重且到达i点的时间是否在i点所要求的时间窗内if WaterDrop(k).Q(WaterDrop(k).UnvisitNode(i))<Qmax; WaterDrop(k).FK=[WaterDrop(k).FK,WaterDrop(k).UnvisitNode(i)];endend
%             for i=alt:N
%                 if ismember(i,WaterDrop(k).UnvisitNode)
%                % if sum(ismember(WaterDrop(k).UnvisitNode,i))
%                     WaterDrop(k).Q(i)=WaterDrop(k).Q(WaterDrop(k).Source)+q(i);
%                     WaterDrop(k).S(i)=WaterDrop(k).S(WaterDrop(k).Source)+ss(WaterDrop(k).Source)+T(WaterDrop(k).Source,i);
%                    % 判断i点的载重量是否小于车辆最大载重且到达i点的时间是否在i点所要求的时间窗内
%                     if WaterDrop(k).Q(i)<Qmax && WaterDrop(k).S(i)>=E(i) && WaterDrop(k).S(i)<=L(i)
%                         WaterDrop(k).FK=[WaterDrop(k).FK,i];
%                     end
%                 end
%             endif isempty(WaterDrop(k).FK)WaterDrop(k).Target=1;else% 计算去下一个可以服务的需求点的概率% 判断从source(k)到下一个可以服务的需求点的路径上泥土量的最小值是否小于0Minsoil=0;for u=1:length(WaterDrop(k).FK)if WaterDrop(k).SW(WaterDrop(k).Source,WaterDrop(k).FK(u))<MinsoilMinsoil=WaterDrop(k).SW(WaterDrop(k).Source,WaterDrop(k).FK(u));endend% 求下一个可以服务需求点对应的函数f之和（可以增加改进，调整节点选择概率）SumF=0;for u=1:length(WaterDrop(k).FK)g(WaterDrop(k).Source,WaterDrop(k).FK(u))=WaterDrop(k).SW(WaterDrop(k).Source,WaterDrop(k).FK(u))-Minsoil;f(WaterDrop(k).Source,WaterDrop(k).FK(u))=1/(0.01+g(WaterDrop(k).Source,WaterDrop(k).FK(u)));SumF=SumF+f(WaterDrop(k).Source,WaterDrop(k).FK(u));end% 求可以服务需求点对应的概率P=[];U=[];for u=1:length(WaterDrop(k).FK)P=[P,f(WaterDrop(k).Source,WaterDrop(k).FK(u))/SumF];U=[U,WaterDrop(k).FK(u)];end% 按照概率P(source(k),u), 用轮盘赌法从FV(k) 集合中随机选择一个点WaterDrop(k).Target=U(RoutedGame(P));          end% 更新水滴的速度WaterDrop(k).Vel=WaterDrop(k).Vel+av/(bv+cv*WaterDrop(k).SW(WaterDrop(k).Source,WaterDrop(k).Target)^2);Maxtime=max(epsilon,WaterDrop(k).Vel);% 计算水滴流动到目标点需要的时间TT(WaterDrop(k).Source,WaterDrop(k).Target)=D(WaterDrop(k).Source,WaterDrop(k).Target)/Maxtime;% 计算水滴携带泥土的增量DeltaSW(WaterDrop(k).Source,WaterDrop(k).Target)= as/(bs+cs*TT(WaterDrop(k).Source,WaterDrop(k).Target)^2);% 更新水滴流过的路径上的泥土量WaterDrop(k).SW(WaterDrop(k).Source,WaterDrop(k).Target)=...(1-alpha)*WaterDrop(k).SW(WaterDrop(k).Source,WaterDrop(k).Target)-alpha*DeltaSW(WaterDrop(k).Source,WaterDrop(k).Target);% 更新水滴携带的泥土量WaterDrop(k).Soil=WaterDrop(k).Soil+DeltaSW(WaterDrop(k).Source,WaterDrop(k).Target);
%             %早到晚到惩罚系数
%             cfe1=0;
%             cfe2=0;
%             cfe0=0;WaterDrop(k).S(WaterDrop(k).Target)=WaterDrop(k).S(WaterDrop(k).Source)+ss(WaterDrop(k).Source)+T(WaterDrop(k).Source,WaterDrop(k).Target);
%             if WaterDrop(k).S(WaterDrop(k).Target)<E
%                 cfe1=EP(WaterDrop(k).Target)*abs(E-WaterDrop(k).Target);
%             elseif WaterDrop(k).S(WaterDrop(k).Target)>L
%                 cfe1=LP(WaterDrop(k).Target)*abs(WaterDrop(k).Target-L);
%             end
%             cfe0=cfe1+cfe2;if  WaterDrop(k).Target==1% （有序集合）配送车辆返回配送中心WaterDrop(k).VisitNode=[WaterDrop(k).VisitNode,WaterDrop(k).Target];% 目标函数中增加动用一辆车的固定成本（表示形成了一辆车的回路）WaterDrop(k).ZZ=WaterDrop(k).ZZ+h*D(WaterDrop(k).Source,WaterDrop(k).Target)+D(WaterDrop(k).Source,WaterDrop(k).Target)*(R0+R1*(q(WaterDrop(k).Target)))+R;WaterDrop(k).Source=WaterDrop(k).Target;WaterDrop(k).Q(WaterDrop(k).Source)=0;WaterDrop(k).S(WaterDrop(k).Source)=0;else% 访问Target(k)点，已访问点增加1个WaterDrop(k).VisitNode=[WaterDrop(k).VisitNode,WaterDrop(k).Target];% 未访问点减少一个WaterDrop(k).UnvisitNode(WaterDrop(k).UnvisitNode==WaterDrop(k).Target)=[];  %%%%%？？？？？% 目标函数只增加车辆运行成本WaterDrop(k).ZZ=WaterDrop(k).ZZ+h*D(WaterDrop(k).Source,WaterDrop(k).Target)+D(WaterDrop(k).Source,WaterDrop(k).Target)*(R0+R1*(q(WaterDrop(k).Target)));% 计算第k辆车到达下一个点时的装载量WaterDrop(k).Q(WaterDrop(k).Target)=WaterDrop(k).Q(WaterDrop(k).Source)+q(WaterDrop(k).Target);% 计算第k辆车到达下一个点的时刻WaterDrop(k).S(WaterDrop(k).Target)=WaterDrop(k).S(WaterDrop(k).Source)+...ss(WaterDrop(k).Source)+T(WaterDrop(k).Source,WaterDrop(k).Target);% 将下一个点作为出发点WaterDrop(k).Source=WaterDrop(k).Target;end %如果第k个水滴的路径已经完成（即所有需求点都被服务完毕）endif WaterDrop(k).ZZ<ZZ=WaterDrop(k).ZZ;Route=WaterDrop(k).VisitNode;Soil=WaterDrop(k).Soil;endend% 对本次循环得到的最优解对应的路径Route上的泥土量进行更新M=size(Route,2);% 路径以外的边上泥土量不变。for i=1:M-1W(Route(i),Route(i+1))=(1+beta)* W(Route(i),Route(i+1))-beta*Soil/(M-1);end%更新迭代次数t=t+1;
end
TotalZ
TotalRoute
LL
DrawPath(TotalRoute,position);
%迭代图
figure
x=1:1:Iter;
y=bestZ;
plot(x,y,'r--')
%plot(x,y);
hold on
%plot(x,y1,'r--');%适应值平均数
title('优化过程')
xlabel('迭代次数')
ylabel('最优适应值')
axis([0,Iter,500,1200])

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

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