文章目录
前言： 目前已经做过几种类型的VRP问题了，从最初只会多种交通工具到后来学会如何加时间窗已经充电站，总算是学会了基本的求解。本来想把代码存到GitHub上面去的，结果网络实在不给力，总是上传失败，所以干脆再花点时间写一写求解的过程。
1.背景
假设在一个供求关系系统中，车辆从货源取货，配送到对应的若干配送点。车辆存在最大载货量，且配送可能有时间限制。需要合理安排取货时间，组织适当的行车路线，使用户需求得到满足，同时使某个代价函数最小，比如总工作时间最少、路径最短等。
可以看出TSP问题是VRP问题的一种简单特殊形式。因此，VRP也是一种NP hard 问题。
目前解决此种问题的办法有多种，主要以启发式算法为主。包括退火算法、遗传算法、蚁群算法、禁忌算法等。
本文介绍该问题的两种情形的求解：1、单种运输工具带时间窗；2、多种运输工具不带时间窗。
2. 单种运输工具带时间窗
该类问题的算法实现参照了高升的硕士学位论文《基于电动汽车的带时间窗的路径优化问题研究》。
2.1 带时间窗车辆路径问题的描述
带时间窗的车辆路径优化问题(Vehicle Routing Problems with Time Windows，VRPTW)是以经典的车辆路径优化问题为基础，新增了每个客户点都有物流配送时间的限制。我们称客户点的时间限制为时间窗。带时间窗的物流配送服务在实际的配送过程中会出现三种情况。
- 第一种情况，物流配送车辆早于客户要求的最早时间到达。这种情况下，配送车辆需要等待到客户要求的最早时间点才能进行卸货。
- 第二种情况，物流配送车辆正好在客户点规定的时间窗内到达，配送车辆到达后就可以直接进行卸货。
- 第三种情况，物流配送车辆晚于客户点规定的最晚时间到达客户点。在这种情况下会出现两种结果。
第一种结果，客户点拒绝接收货物。我们称这种情形下的时间窗为硬时间窗。硬时间窗的客户点要求物流配送车辆可以早于时间窗规定的最早服务时间到达，但是需要等待到客户点最早服务时间才能开始卸货，配送车辆也可以正好在客户规定的时间内到达，此时配送车辆到达客户点后可以直接进行卸货。
第二种结果，客户点同意接收配送的货物，但是物流配送中心需要接受一定的惩罚。我们称这种情形下的时间窗为软时间窗。软时间窗的客户点允许物流配送车辆在客户点规定的时间外到达。软时间窗与硬时间窗相同的是，如果配送车辆早于客户点规定的时间到达，需要等待到客户点的最早接受服务的时间才能开始卸货；如果配送车辆正好在客户规定的时间内到达，此时配送车辆到达客户点后可以直接进行卸货。软时间窗与硬时间窗不同的是，在软时间窗情形下，配送车辆晚于客户点规定的时间到达时，客户点接受货物；但是在硬时间窗的情形下，如果配送车辆晚于客户点规定的时间窗到达，客户点不会接受货物。
与构建普通的车辆路径优化问题模型不同，构建带时间窗的车辆路径问题模型需要考虑物流配送中心因为违反客户点规定的时间而遭受到的惩罚成本。硬时间窗情形下，惩罚成本只包含配送车辆早到而遭受到的惩罚。而软时间窗情形下，惩罚成本包括两部分，因为配送车辆早到遭受的惩罚成本和因为配送车辆晚到而遭受到来自客户点的惩罚成本。
通常情况下，带软时间窗的车辆路径优化问题更加符合实际情况，因此以软时间窗为例。
2.2 遗传算法的求解
2.2.1 编码方案的设计
根据论文中模型的特点，采用自然数编码方式。染色体长度为n+k+1，其中n表示客户点数量，k为配送车辆数量。配送中心的编码为
0；编码1,2,⋯,n 代表各个客户点被分配的自然数序号。
例如配送中心为8个客户进行配送服务，假如其中一条染色体为：0，3，2，9，6，0，7，10，1，5，0，8，4，0。表示的行驶路线为：第1辆电动汽车的行驶路径为0 3 2 6 0，代表这辆电动汽车从配送中心出发，先后经过客户点3、客户点2与客户点6返回配送中心；第2辆电动汽车的行驶路线为0 7 1 5 0，代表此车从配送中心出发经客户点7，然后经客户点1、客户点5返回配送中心：第3辆电动汽车的行驶路线为0 8 4 0，代表此车从配送中心出发经客户点8、客户点4返回配送中心。
2.2.2 种群初始化
首先将所有的客户点代码随机排成一列，qi′代表染色体中第i个客户点的货物需求量。如果Σi=1aqi′≤Q且Σi=1a+1qi′>Q，则在染色体第a位后面插入0；然后从插入0后开始重复计算直至插入n−1个0，再将染色体首位及最后一位分别插入一个0，最终形成一条初始染色体，反复上述过程产生N条个体构成初始种群。
MATLAB代码(chushihua.m)：
% 生成初始种群
% 函数说明：
% 染色体长度为m+k+1，其中m为客户点数目，k为车辆数目，加1是因为第一个位置为配送中心，数值为0
% 具体的染色体编码与解之间的关系见论文41页的5.2.1部分
% 种群初始化时需插入k个0，这样子可以满足式4.10和4.11，具体插入方式见42页5.2.2，5.2.2可以满足式4.9
function s = chushihua(NIND,m,k,q,Qk)
%% 生成NIND个配送中心，这是放在第一列的
center = zeros(NIND,1);
%% 生成长度为m+k的染色体，用以表示车辆配送路径设置
path_center = zeros(NIND,m+k); % 这样的好处一方面加速运行，另一方面保证每一行最后一个数一定为0
for i = 1:NIND
% 初始的，使用randperm可以保证每个客户点能且仅能被服务一次，即满足式4.8
path = randperm(m);
% 插入k个0
upload = 0; % 当前这辆车已载货物量
kk = 0; % 已使用的车辆数
for j = 1:m
if upload + q(path(j)+1) <= Qk % q的索引为path(j)+1是因为q的第一个数字为配送中心
path_center(i,j+kk) = path(j); % 第j个需求点应该放在第j-1+kk+1的位置上
upload = upload + q(path(j)+1);
else % 当当前车辆已经载货量已达到上限，需要更换一辆车了
path_center(i,j+kk) = 0; % 加kk的原因同上
path_center(i,j+kk+1) = path(j); % 这个再加1因为是放在0的后面
kk = kk+1;
upload = q(path(j)+1);
end
end
end
%% 合并，s的第一列和最后一列都是0，即满足式4.7
s = [center,path_center];
2.2.3 约束处理与适应度函数
常见的约束条件的处理方法有–平op。方法一是将问题的约束在染色体中表现出来的直接处理约束的方法。该方法的适用领域有限，设计专门的染色体和遗传算子较为困难。方法二是在编码的时候不考虑约束，而是在遗传算法的计算过程中对得到的染色体对应的解是否可行进行检测。此方法只适用于约束简单、可行解易于得到的优化问题。本文采用第三种方法一一惩罚函数的方法处理约束。
适应度函数要求越大越好，因此我们取上述目标函数的倒数为适应度函数，即适应度函数为：fit(i)=1/z形。
MATLAB代码(fitness.m):
% 适应度函数
% 目标函数 = 运输成本+惩罚成本+超载成本（最后一项为电量不足产生的成本，在这里不再计算）
% 其中超载出现的原因为交叉或变异，初始化生成的一定不会超载
% 以上内容参考自5.2.3
function [fit,ff] = fitness(s,m,Ck,C1,C2,LT,ET,Qk,q,k,speed,TT,D)
%% 初始化
M = 1000000; % M为很大的数，作为约束条件的惩罚因子
[NIND,~] = size(s); % 这个NIND可以不是种群的大小
d = zeros(NIND,1); % 运输距离
overload = zeros(NIND,k); % 是否超载
%% 适应度计算
% 运输总路程
for i = 1:NIND
for j = 1:m+k
d(i) = d(i) + D(s(i,j)+1,s(i,j+1)+1);
end
end
% 惩罚成本
pun = punish(s,C1,C2,LT,ET,m,k,speed,TT,D);
% k辆车的载重量
L = ostation(k,s);
for i = 1:NIND
for j = 1:k % 计算k辆车的运载量
% 这个计算起来比较复杂，一下子讲不清楚，所以得先自己理解一下
upload(i,j) = sum(q(s(i,L(i,j)+1:L(i,j+1)-1)+1)); % 计算第j辆车的运载量
% 判断是否超载，主要是方便后面计算目标函数
overload(i,j) = max(upload(i,j)-Qk,0);
end
end
% 总成本，最后一项为k辆车总的超重成本
f = Ck*d + pun + sum(M*overload,2); % 列向量求和运算
% 适应度
fit = 1./f;
% 真实的目标函数（不包含惩罚项）
ff = Ck*d + pun;
2.2.4 选择算子
采用轮盘赌的方法，MATLAB代码(sel.m)：
% 选择
% sel_num为选择的个体数目，该函数返回两个索引值
function seln = sel(m,Ck,C1,C2,LT,ET,s,Qk,q,k,speed,TT,D)
sel_num = 2; % 被选择的个体数目
seln = zeros(sel_num,1); % 从种群中选择sel_num个个体
%% 选择概率的计算
% Step1：计算适应度
[fit,~] = fitness(s,m,Ck,C1,C2,LT,ET,Qk,q,k,speed,TT,D);
% Step2-Step3:计算选择概率
fitsum = sum(fit);
p = fit/fitsum;
% Step4:计算累计概率
ps = cumsum(p); % cumsum用于求累计和
%% 个体的选择
for i = 1:sel_num
r = rand;
% 找到比随机数r大的ps中的最小数对应的下标
ind = min(find(ps > r));
% 如果ind还没有被选择过的话才可以被选
while length(find(seln == ind)) ~= 0 % seln中至少有一个值为ind，说明被选中过，那么要重新选
r = rand;
ind = min(find(ps > r));
end
seln(i) = ind;
end
2.2.5 交叉算子
为了改善算法进化后期的寻优能力，本文采用新的改进交叉算子，最大程度的保留父代的优秀子路经信息。具体步骤见论文的43页。
MATLAB代码(cro.m):
% 交叉
function scro = cro(s,m,pc,Ck,C1,C2,LT,ET,Qk,q,k,speed,TT,D)
%% 初始化
[NIND,~] = size(s);
scro=zeros(NIND,m+k+1);
%% 交叉
for i = 1:2:NIND-1
% 选择的两个父代染色体的索引
seln=sel(m,Ck,C1,C2,LT,ET,s,Qk,q,k,speed,TT,D);
% 两个父代染色体
scro(i,:) = s(seln(1),:);
scro(i+1,:) = s(seln(2),:);
% 两个父代染色体0的位置
L = ostation(k,scro(i:i+1,:));
if pc > rand
path = zeros(2,m+k+1);
%% Step1: 分别在两个父代染色体上随机选择一段子路经
path_num1 = round(1 + rand*(k-1)); % 生成1,2,3,...k中的某一个
path_num2 = round(1 + rand*(k-1));
path_elc1 = scro(i,L(1,path_num1):L(1,path_num1+1));
path_elc2 = scro(i+1,L(2,path_num2):L(2,path_num2+1));
%% Step2：被选择的子路段前置
% 备份scro
scro_copy1 = scro(i,:);
scro_copy2 = scro(i+1,:);
% 父代染色体1
path(1,1:length(path_elc1)) = path_elc1;
scro(i,1:length(path_elc1)) = path_elc1; % 被选中的段落移前面去
scro_copy1(L(1,path_num1):L(1,path_num1+1)-1) = []; % 删除前移部分，还要-1是因为得留下一个0
scro(i,length(path_elc1):end) = scro_copy1; % 剩下部分在后面接上去
% 父代染色体2
path(2,1:length(path_elc2)) = path_elc2;
scro(i+1,1:length(path_elc2)) = path_elc2;
scro_copy2(L(2,path_num2):L(2,path_num2+1)-1) = [];
scro(i+1,length(path_elc2):end) = scro_copy2;
%% Step3：参见论文43页最下方，以下步骤得到的path的后两个一定为0
% 补全子染色体1
for j = 1:m+k-1-length(path_elc1) % 后面还有两个0的位置，故为m+k+1-2
for ij = 2:m+k % 遍历子染色体2，从第2个遍历到倒数第二个
if length(find(path(1,:) == scro(i+1,ij))) == 0 % 如果被遍历的位置对应的需求点没在1当中
path(1,length(path_elc1)+j) = scro(i+1,ij);
break;
end
end
end
% 补全子染色体2
for j = 1:m+k-1-length(path_elc2)
for ij = 2:m+k % 遍历子染色体1
if length(find(path(2,:) == scro(i,ij))) == 0
path(2,length(path_elc2)+j) = scro(i,ij);
break;
end
end
end
%% Step4：选择一个最好的位置插入0
% 子染色体1
best_fit1 = 0;
best_s1 = zeros(1,m+k+1);
path_copy = path;
for j = 1:m+k-2-length(path_elc1) % 0不能连着存在，因此有3个位置不能插入0
% path_copy自插入位置(length(path_elc1)+1+j)开始，后面的那些非0数统一往后移1位
path_copy(1,length(path_elc1)+1+j+1:m+k) = path(1,length(path_elc1)+1+j:m+k-1);
% 插入位置为0
path_copy(1,length(path_elc1)+1+j) = 0;
% 插入位置之前的部分照常放
path_copy(1,1:length(path_elc1)+j) = path(1,1:length(path_elc1)+j);
% 求此时的适应度
[fit,~] = fitness(path_copy(1,:),m,Ck,C1,C2,LT,ET,Qk,q,k,speed,TT,D);
if fit > best_fit1 % 如果插更换入0的位置之后适应度更大，则更新最优解
best_fit1 = fit;
best_s1 = path_copy(1,:);
end
end
% 子染色体2
best_fit2 = 0;
best_s2 = zeros(1,m+k+1);
for j = 1:m+k-2-length(path_elc2) % 0不能连着存在，因此有3个位置不能插入0
ind = length(path_elc2)+1+j;% 记录插入位置
% path_copy自插入位置ind开始，后面的那些非0数统一往后移1位
path_copy(2,ind+1:m+k) = path(2,ind:m+k-1);
path_copy(2,ind) = 0;
path_copy(2,1:ind-1) = path(2,1:ind-1);
[fit,~] = fitness(path_copy(2,:),m,Ck,C1,C2,LT,ET,Qk,q,k,speed,TT,D);
if fit > best_fit2
best_fit2 = fit;
best_s2 = path_copy(2,:);
end
end
path = [best_s1;best_s2];
scro(i:i+1,:) = path;
end
end
end
2.2.6 变异算子
MATLAB代码(mut.m)：
% 变异操作
function smut = mut(s,m,k,Ck,C1,C2,LT,ET,Qk,q,speed,TT,pm,D)
[NIND,col] = size(s); % col = m+k+1
smut = s;
ss = zeros(6,col);
%% 变异
% 不妨记三个位置为1,2,3
for i = 1:NIND
if pm > rand
%% 找出三个交换的位置
r = randperm(m);
exchange_num = r(1:3); % 随机序列的前3个作为交换位置
ind_1 = find(s(i,:) == exchange_num(1)); % 1
ind_2 = find(s(i,:) == exchange_num(2)); % 2
ind_3 = find(s(i,:) == exchange_num(3)); % 3
ss(1,:) = s(i,:);
%% 5种情况
% 1,3,2
ss(2,:) = exchange(s(i,:),ind_2,ind_3);
% 2,1,3
ss(3,:) = exchange(s(i,:),ind_1,ind_2);
% 2,3,1
ss(4,:) = exchange(s(i,:),ind_1,ind_2);
ss(4,:) = exchange(ss(4,:),ind_2,ind_3);
% 3,1,2
ss(5,:) = exchange(s(i,:),ind_1,ind_3);
ss(5,:) = exchange(ss(5,:),ind_2,ind_3);
% 3,2,1
ss(6,:) = exchange(s(i,:),ind_1,ind_3);
%% 找到5种情况里最好的
[fit,~] = fitness(ss,m,Ck,C1,C2,LT,ET,Qk,q,k,speed,TT,D);
[~,ind_max] = max(fit);
smut(i,:) = ss(ind_max,:);
end
end
end
2.2.7 其他部分
- 精英个体选择(
reins.m)
% 保留父代精英个体
function [new_s,best_s,best_fit]=reins(s,smut,m,Ck,C1,C2,LT,ET,Qk,q,k,speed,TT,D,NIND)
ss = [s;smut];
[fit,ff] = fitness(ss,m,Ck,C1,C2,LT,ET,Qk,q,k,speed,TT,D);
% 取原种群与经交叉变异后得到的种群中最优的NIND个个体
[~,index] = sort(fit,'descend');
new_s = ss(index(1:NIND),:); % 按适应度从小到大排序的个体
new_fit = fit(index(1:NIND));
best_s = ss(index(1),:); % 这一代最佳染色体
best_fit = fit(index(1)); % 这一代最优适应度
end
- 惩罚函数(
punish.m)
% 惩罚函数
% 惩罚成本的计算公式见论文14页的2.9
function pun = punish(s,C1,C2,LT,ET,m,k,speed,TT,D)
[NIND,~] = size(s);
%% 参数设置
pun = zeros(NIND,1);
T = my_time(NIND,m,k,s,ET,speed,TT,D); % 计算到达每一个点的时刻
%% 惩罚计算
for i = 1:NIND
for j = 1:m % 这里j是对m各需求点按序号进行遍历
% 已有的惩罚加上当前点(第j个)产生的惩罚，而惩罚分早到与晚到两种，根据式4.21算得
pun(i) = pun(i) + C1*max(ET(j)-T(i,j),0) + C2*max(T(i,j)-LT(j),0);
end
end
- 查找染色体0的位置(
ostation.m)
% 函数功能：确定0的位置,把他赋值给L
% 如果知道每一条染色体0的位置，我们就可以将每一辆车的路径与染色体片段对应起来，方便求解目标函数,
% 在本程序中即更加方便求解适应度
function L=ostation(k,s)
[NIND,~] = size(s);
L = zeros(NIND,k+1); % 此处将4改为k+1，尽量避免硬编码
for i = 1:NIND
L(i,:) = find(s(i,:) == 0);
end
- 计算达到客户点的时间(
my_time.m)
% 该函数用于计算到达每一个客户点的时间
% 主要参考的内容在论文37页上半部分
function T = my_time(NIND,m,k,s,ET,speed,TT,D)
T = zeros(NIND,m); % 用以记录每个解中每个需求点的供货车辆到达时间
for i = 1:NIND
time_spend = 0; % 满足式4.12，从配送中心出发的时间是0时刻
for j = 1:m+k+1 % 遍历染色体的每个基因
if s(i,j) ~= 0
% 到达时间 = 已经花费的时间(离开上一个点的时间)+运输时长，满足4.15
T(i,s(i,j)) = time_spend + D(s(i,j-1)+1,s(i,j)+1)/speed; % 满足式4.13
% 离开目前这个点的时间=到达时间+早到停留时长+客户停留时长(TT)，满足4.14与4.16
time_spend = T(i,s(i,j)) + max(ET(s(i,j))-T(i,s(i,j)),0) + TT(s(i,j)) ;
else % 如果s(i,j)==0，说明换了辆车，已花费时间就要清零
time_spend = 0;
end
end
end
end
- 交换染色体元素(
exchange.m)
% 此函数主要用于交换s两个位置的元素
function s = exchange(s,ind_1,ind_2)
t = s(ind_1);
s(ind_1) = s(ind_2);
s(ind_2) = t;
end
- 欧氏距离(
Distance.m)
%计算两两需求点之间的欧式距离
function D = Distance(X)
row = size(X,1); % 城市个数（在这里是k+m）
D = zeros(row,row);
for i = 1:row
for j = i+1:row
D(i,j) = ((X(i,1) - X(j,1))^2 + (X(i,2) - X(j,2))^2)^0.5;
D(j,i) = D(i,j);
end
end
- 主函数(
main.m)
%% 注意
% 不用电动汽车，然后1个运输中心，20个客户，3辆车
% time 函数的子函数，就是时间窗怎么设置，我也需要知道。
% 重点讲解ostation、chushihua函数
% 需要的结果：
% ①：路径规划与每一条路径的配送成本（起点和终点都是运输中心）
% ②：绘出路径图
clear,clc
%% 参数设置
k=3; % 3辆车
m=20; % 20个客户点
% ch=2; % 2个充电站
% 各客户点需求量，其中配送中心需求量为0，其中q(1)为配送中心，以此类推
q = [0 0.5 0.8 0.4 0.5 0.7 0.7 0.6 0.2 0.2 0.4 0.1 0.1 0.2 0.5 0.2 0.7 0.2 0.7 0.1 0.5];
% 配送中心和客户点的坐标位置，配送中心在第1位
X=[ 56 56
40 48
32 80
16 69
88 96
48 96
32 104
80 56
48 40
24 16
48 8
16 32
8 48
32 64
24 96
72 104
72 32
72 16
88 8
104 56
104 32 ];
D = Distance(X); % 各个点之间的距离
% 时间参数
% 客户要求到货的时间始点（最早点）
ET=[3 0 7 1 4 1 3 0 2 2 7 6 7 1 1 8 6 7 6 4];
% 客户要求到货的时间终点（最晚点）
LT=[5 2 8 3 5 2 4 1 4 3 8 8 9 3 3 10 10 8 7 6];
% 客户点的停留时间
TT=[0.5 0.8 0.4 0.5 0.7 0.7 0.6 0.2 0.2 0.4 0.1 0.1 0.2 0.5 0.2 0.7 0.2 0.7 0.1 0.5];
% 成本参数
Ck=10; % 第k辆车运行单位距离的费用（运输成本）
C1=20; % 车辆在任务点处等待单位时间的机会成本（早到惩罚）
C2=30; % 车辆在要求时间之后到达单位时间所处的惩罚值（晚到惩罚）
% 汽车参数
Dis = 160; % 续驶里程
Qk = 4; % 车辆额定载重量
speed = 40; % 汽车的行驶速度
% 遗传算法参数
NIND = 100; % 种群大小，一般种群大小在100到200之间比较好
maxgen = 500; % 遗传代数
pc = 0.8; % 交叉概率
pm = 0.1; % 变异概率
% GGAP = 0.5; % 代购，决定了选择操作部分选择的个体数目
%% 模型求解
s = chushihua(NIND,m,k,q,Qk); % 初始化种群
gen=0; % 表示当前遗传代数
% 计算适应度和模型目标函数值
[~,ff] = fitness(s,m,Ck,C1,C2,LT,ET,Qk,q,k,speed,TT,D);
preff=min(ff); % 初始种群的最优目标函数值
sperfect=zeros(maxgen,m+k+1); % 每一代的最优个体
best_fit = zeros(maxgen,1); % 每一代最优适应度
while gen<maxgen
%% 遗传算子
% 交叉操作
scro = cro(s,m,pc,Ck,C1,C2,LT,ET,Qk,q,k,speed,TT,D);
% 变异操作
smut = mut(scro,m,k,Ck,C1,C2,LT,ET,Qk,q,speed,TT,pm,D);
%% 记录并更新全局最优解
% 重插入子代
[s,s_gen,best_fit_gen] = reins(s,smut,m,Ck,C1,C2,LT,ET,Qk,q,k,speed,TT,D,NIND);
% 更新迭代次数
gen=gen+1;
% 记录当前代最好的染色体
sperfect(gen,:)=s_gen;
% 记录这一代最好的适应度
best_fit(gen) = best_fit_gen;
% gen
end
%% 模型输出
best_path = sperfect(end,:);
% 输出最优解的路线和最优目标函数值
disp('最优解：')
disp(best_path) % 最优路线
disp(['目标函数值：',num2str(1/best_fit(end))])
%% 绘图
% 遗传图
figure;
hold on; box on;
plot(best_fit)
xlim([0,maxgen]); % 设置x坐标轴显示范围
title('优化过程');
xlabel('代数');
ylabel('适应度值');
% 配送路径图
figure
plot(X(1,1),X(1,2),'o','markersize',10) % 配送中心
hold on
plot(X(2:end,1),X(2:end,2),'.','markersize',10) % 需求点
for i = 1:m+k
line([X(best_path(i)+1,1),X(best_path(i+1)+1,1)],[X(best_path(i)+1,2),X(best_path(i+1)+1,2)]);
end
xlabel('横坐标X')
ylabel('纵坐标Y')
2.3 粒子群算法的求解
- 主函数(
main.m)
%%—————粒子群算法求解带时间窗的车辆路径规划问题——————%%%
clc
clear all
close all
%% 目标约束初始化
%同时需要给出距离矩阵（DS）--DS
DS=load('算法求解\DS.mat');
DS=struct2cell(DS);
DS=cell2mat(DS);
%车辆速度（speed)
speed=50;
%各任务的时间窗([ETi,LTi])---CT
CT=load('算法求解\CT.mat');
CT=struct2cell(CT);
CT=cell2mat(CT);
%服务时间（STi)--ST
ST=load('算法求解\ST.mat');
ST=struct2cell(ST);
ST=cell2mat(ST);
%车辆容量（W）
w=7.5;
%各场点的货运量（gi)--g
g=load('算法求解\g.mat');
g=struct2cell(g);
g=cell2mat(g);
%运输成本矩阵y（与运输距离成正比）
y=3;
%定义目标函数的罚金成本PE,PL
PE=50;%早到时间惩罚
PL=50; %到达加上该场点的服务时间达到LT的惩罚
%% 粒子初始化
n=40; %粒子规模
c1=1.49; %参数1
c2=2.33; %参数2
m=9; %%任务数
cdm=linspace(1,m,m); %%生成任务编号
vn=4; %%总的车辆数
Xv=zeros(n,m); %各任务对应的车辆编号，为0-vn之间的整数
Xr=zeros(n,m); %为各车辆进行各项任务的顺序
Vv=zeros(n,m); %为车辆编号的速度
Vr=zeros(n,m); %为任务顺序的速度
maxgen=100; %迭代最大次数
pbestXr=zeros(n,m); %%各粒子的在历史过程中的最优解
pbestXv=zeros(n,m); %%各粒子的在历史过程中的最优解
gbestXr=zeros(1,m); %%粒子种群在历史过程中的最优
gbestXr=zeros(1,m); %%粒子种群在历史过程中的最优
%% 算法实现
%i=1时
Xv=round((rand(n,m)*(vn-1)+1)); %为0-vn之间的整数
Xr=rand(n,m)*(m-1)+1; %为0-m之间的实数
Vv=round((rand(n,m)*(vn-3)+1)); %为0-(vn-1)之间的整数
Vr=(rand(n,m)*(m-1)+1); %为0-(m-1)之间的实数
[gbestXv gbestXr fav]=eval(Xv,Xr,y,PE,PL,CT,ST,g,vn,n,m,speed,DS,w); %%进行粒子评价并得到粒子的gbestXr,gbestXv
pbestXv=Xv;
pbestXr=Xr;
fav1=fav; %评估值进行存储，用于下一状态的评价
%当i>1时
for i=2:maxgen
Xv2=Xv;
Xr2=Xr;
for j=1:n
Xv1=Xv(j,:);
Xr1=Xr(j,:);
Vv1=Vv(j,:);
Vr1=Vr(j,:);
Vv(j,:)=round(Vv1+c1*rand*(pbestXv(j,:)-Xv1)+c2*rand*(gbestXv-Xv1)); %%各任务对应的车辆速度更新
Vr(j,:)=Vr1+c1*rand*(pbestXr(j,:)-Xr1)+c2*rand*(gbestXr-Xr1); %%车辆接受各任务的序号的速度更新
%%%速度校正
Vvsz=Vv(j,:);
Xvsz=Xv(j,:);
Vrsz=Vr(j,:);
Xrsz=Xr(j,:);
for j1=1:m
if Vvsz(j1)>vn-1 | Vrsz>m-1
Vvsz(j1)=vn;
Vrsz(j1)=m-1;
elseif Vvsz(j1)<1-vn | Vrsz<1-m
Vvsz(j1)=1-vn;
Vrsz(j1)=1-m;
end
end
%%%位置校正
Xv(j,:)=Xv(j,:)+Vv(j,:); %%各任务对应的车辆序号位置更新
Xr(j,:)=Xr(j,:)+Vr(j,:); %%各车辆进行各项任务顺序的位置更新
%------------------------------------%%
if Xvsz(j1)<=0 | Xrsz<1-m
Xvsz(j1)=1;
Xrsz(j1)=rand*(1-m);
elseif Xvsz(j1)>vn | Xrsz>m-1
Xvsz(j1)=round(rand*(vn)+1);
Xrsz(j1)=rand*(m-1);
end
Vv(j,:)=Vvsz;
Xv(j,:)= Xvsz;
Vr(j,:)=Vrsz;
Xr(j,:)= Xrsz;
end
[gbestXv gbestXr fav]=eval(Xv,Xr,y,PE,PL,CT,ST,g,vn,n,m,speed,DS,w); %%进行粒子评价
for k=1:n
if fav(k)<=fav1(k)
pbestXv(k,:)=Xv(k,:);
pbestXr(k,:)=Xr(k,:);
else
pbestXv(k,:)=Xv2(k,:);
pbestXr(k,:)=Xr2(k,:);
end
end
fav1=fav;
end
gbestXv;
min(fav);
gbestXr;
%% 最优粒子解码
zuiyou=zeros(vn+1,m);
zuiyou(1,:)=linspace(1,m,m);
xh=0;
for i=1:vn
colv1=find(gbestXv==i);
for j=1:length(colv1)
zuiyou(i+1,colv1(j))=gbestXr(colv1(j));
end
zy0=zuiyou(1,:)';
zy1=zuiyou(i+1,:)';
zy2=[zy0 zy1];
szy=sortrows(zy2,2);
szy=szy';
t=0;
for k=1:m
if szy(2,k)~=0
% xh
t=t+1;
szy(1,k);
xh(t)=szy(1,k);
end
end
xh;
st=0;
st=num2str(st);
%% 结果输出
for d=1:length(xh)
strr=num2str(xh(d));
strr=strcat(strr,'->');
st=[st,'->', strr];
end
st1=0;
st1=num2str(st1);
st=[st,st1];
disp(['第',num2str(i),'辆车','的运输路线为：',num2str(st)])
xh=0;
end
disp(['目标最优值为：',num2str(min(fav))])
- 例子评估函数(
eval.m)
%%----评估粒子---%%
%是否为可行解--判断依据：各车辆是否超重超重
%根据目标函数计算目标值，找出粒子群的gbest和各粒子的pbest
%i,j,k,c,t,j1,j2
function [gbestXv gbestXr fav]=eval(Xv,Xr,y,PE,PL,CT,ST,g,vn,n,m,speed,DS,w)
fav=zeros(n,1); %%%
bm=linspace(1,n,n);
ET=CT(:,1);
LT=CT(:,2);
g;
t1=0;
for i=1:n
Xvi=Xv(i,:); %获取任务对应车辆编号的粒子行
Xri=Xr(i,:); %获取车辆对应任务顺序行
st=zeros(vn+1,m);
st(1,:)=linspace(1,m,m);
t=0;
bmzys=zeros(1);
sg=zeros(1,vn); %存储没辆车的货物
ST1=zeros(1,vn);%存储每辆车的时间惩罚
juli=zeros(1,vn);
%对每个粒子的Xv和Xr进行解码
for j=1:vn
colv=find(Xvi==j);
for k=1:length(colv)
st(j+1,colv(k))=Xri(colv(k)); %获取每一辆车的Xr值
end
st0=st(1,:)';
stj=st(j+1,:)';
sst=[st0 stj];
pxst=sortrows(sst,2); %%进行排序,获取每辆车接送任务的顺序
for c=1:m
if pxst(c,2)~=0
t=t+1;
cxh(t)=pxst(c,1); %获取每辆车的进行各项任务的编号
end
end
for j1=1:t
sg(j)=sg(j)+g(cxh(j1)); %%每辆车所载货物
if j1~=t
juli(j)=juli(j)+DS(1,cxh(1))+DS(cxh(j1),cxh(j1+1)); %每辆车途经距离
else
juli(j)=juli(j)+DS(cxh(t),1); %每辆车途经距离
end
if t~=1
if j1~=t
cxh(j1+1);
cxh(j1);
ST0=DS(1,cxh(1))/speed+ST(cxh(1))+PE*max(ET(cxh(1))-ST(1),0)+PL*max(LT(cxh(1))-ST(1),0); %每辆车途经第一个节点时的时间（不包括时间惩罚）
ST1(j)=ST0+ST1(j)+DS(cxh(j1),cxh(j1+1))/speed+PE*max(ET(cxh(j1+1))-ST1(j),0)+PL*max(LT(cxh(j1+1))-ST1(j),0)+ST(cxh(j1));
else
ST1(j)=ST1(j)+DS(cxh(t),1)/speed;
end
else
ST0=DS(1,cxh(1))/speed+ST(cxh(1))+PE*max(ET(cxh(1))-ST(1),0)+PL*max(LT(cxh(1))-ST(1),0);
ST1(j)=ST0;
end
end
if sum(sg)>w %%约束条件判断
t1=t1+1;
bmzys(t1)=j;
end
end
%%%目标函数计算
fav(i)=y*sum(juli)+sum(ST1);
end
bmzyscl=find(bmzys~=0); %%读取不满足约束的粒子
for j3=1:length(bmzyscl)
bm(bmzys(bmzyscl(j3)))=0;
fav(bmzys(bmzyscl(j3)))=inf;
end
bm1=bm(bm~=0);
for j4=1:length(bmzyscl)
lh=length(bm1);
lh1=round(rand*(lh-2)+1);
Xv(j4,:)=Xv(lh1,:);
end
gbestj=find(fav==min(fav));
gbestXr=Xr(gbestj,:);
gbestXv=Xv(gbestj,:);
3. 多种运输工具不带时间窗
多种运输工具的编码等算子参考了李红蕾的硕士学位论文《立体交通条件下医药物流配送路径优化研究》
3.1 立体交通介绍
随着城市化进程与经济活动集约化的快速发展，城市内部的人口数量激增，人们自有汽车的数量也随之呈现出不断上涨的趋势，传统的地面交通系统已经不能很好地解决交通拥堵问题。为了缓解现有的交通压力，城市开始逐步开展立体交通运输体系建设，开始向地面、地下、空中的空间范围内发展。总的来说，立体交通概念的界定主要体现在以下几个方面：
- 向地下、空中的空间范围内发展；
- 包含多种不同类型的交通运输方式；
- 形成了立体化综合运输的网络格局。
简单来说，立体交通就是对地面、地下以及空中三个空间范围内的运输方式的综合利用，比如说汽车、地铁、轻轨、火车、飞机等。通过选用合理的运输工具来实现运输效率的最大化也是发展立体交通的根本目的。在城市内部之间选用轨道交通、快速公交、出租车等运输工具进行运输，在城际之间选用火车、飞机等进行对外交通联系。本文在进行立体交通条件下的 VRP 问题研究时，选用了汽车、火车、飞机三种不同空间范围内的不同类型的运输工具共同参与配送，目的是为了验证立体交通的发展有助于提高配送环节的效率、降低配送成本，有利于实现整个运输网络的优化。
3.2 遗传算法求解
3.2.1 染色体编码设计
本文研究的是立体交通条件下的 VRP 问题，通过第三章可以看出，本文在构建配送模型时，将多种不同类型的运输工具纳入了配送范围，因此在进行编码方式设计时，必须要将运输工具车型不同的因素考虑进去。与此同时，还要满足该运输工具的服务总量不超过其载重量的约束条件。
由于考虑到因素相对较多，因此本文选用实数编码的方式来设计染色体的编码方案。现假设只有一个配送中心，用数字 0 来表示；假设该配送中心拥有三种不同的类型的运输工具可供调用，用编号Ⅰ、Ⅱ、Ⅲ分别进行表示；假设此次配送方案中有 9 个客户点向配送中心提出了配送要求，用自然数 1~9 来进行表示。此时可以用如下的编码串表示配送方案：4 6 12 9 1 7 10 8 11 5 3 2。通过对上述编码串的解码可以得到一个对应的配送方案。解码过程如下：
首先要把 10 以下的整数找出来按上述编码串的顺序进行排列，形成一个将所有客户点包括在内的子集，此时子集内部的每一个数字都代表一个客户，数字排列的先后顺序与客户被服务的顺序是一致的，所以说这个可行解对应的客户被服务的顺序是 4 6 9 1 7 8 5 3 2。而编码方案中出现的 10 以上的数字排序则是对不同类型的运输工具的顺序安排。
对于客户服务线路的安排，由于运输工具有配载量的限制，因此在安排配送线路时必须要考虑配载量的限制。编号为 4 的客户作为第一个被服务的对象，加入到第一条配送路径中，此时需要将客户 4 的需求量与运输工具的载重量的大小进行比较，若需求量较小，则可以将编号为 6 的客户纳入第一条配送路径之中。此时仍需进行 4、6 两个客户的需求总量和载重量比较的过程。假设需求总量较小，继续按顺序将编号为 9 的客户纳入到第一条配送路径当中，再次判断需求总量与载重量的大小。若需求总量仍小于载重量，继续将编号为 1 的客户纳入同一配送线路。假设此时需求总量大于运输工具的载重量，则说明编号为 1 的客户不能加入到第一条配送路径中。此时得到的第一条配送路径为0-4-6-9-0，表示配送中心的运输工具从中心出发后，按照 4,6,9 的先后顺序对客户进行服务，最后返回原配送中心的一条路线。
对于运输工具类型的安排，在编码过程中，配送中心的Ⅰ、Ⅱ、Ⅲ类型的运输工具分别用实数 10、11、12 进行表示。在上边的例子中，三个 10 以上的实数排列顺序为 12、10、11，其对应代表的运输工具类型为Ⅲ、Ⅰ、Ⅱ。在安排运输工具时，首先安排进行配送任务的就是类型为Ⅲ的运输工具，也就是说第一条配送路径对应的的运输工具类型为Ⅲ；第二条配送路径对应的类型为Ⅰ；如果存在第三条配送路径，则对应的运输工具类型为Ⅱ。
采用上述编码方式，考虑到了立体交通条件下的多种运输工具参与配送的实际情况，通过这种新的编码方式，能够更直观地找出哪些客户都由哪种类型的运输工具进行服务的。在解决多类型运输工具配送的问题时，改进后的编码方式更为实用。
MATLAB代码(decode.m)
function [car_path,upload_num,dist_all] = decode(chrom,demand,N,max_vc,max_v,D)
% 该函数用于解码，传入的chrom为一个个体,N为结点个数,max_vc为最大载重
% max_v为每一种车辆的最大可用量,demand为每个点的需求量,D为距离矩阵
%可以得到：
% 1：个体的目标函数值
% 2：个体对应解的路径与其车辆安排
%% 参数设置
% 分开
renwu = chrom(find(chrom<=N)); % 分出来的任务
cheliang = chrom(find(chrom>N)); % 分出来的车辆
j = 1; % 用来记录第几辆车
num_cheliang = zeros(1,length(max_v)); % 每一种运输工具的使用量
k = 1; % 用来记录当前是第几条路径
m = 1; % 用来记录当前路径已经有多少个货物需求点了
upload_num = 0; % 所有路径上已经转载的货物量
dist_all = []; % 所有路径的总长度
%% 选择初始路径交通工具并确定下来最大转载量
car_path = []; % 这里面第一个数字为路径运输工具类型及路径(第一个为类型，1位货车，2为火车，3位飞机)
if cheliang(j)>N && cheliang(j) <= N+max_v(1)
num_cheliang(1) = num_cheliang(1)+1; % 使用第一种车辆的数目增加1
capacity = max_vc(1); % 这一条线路最大运输量
car_path(k,1) = 1;
else if cheliang(j) >N+max_v(1) && cheliang(j) <= N+max_v(1)+max_v(2)
num_cheliang(2) = num_cheliang(2)+1;
capacity = max_vc(2);
car_path(k,1) = 2;
else if cheliang(j) > N+max_v(1)+max_v(2) && cheliang(j) <= N+max_v(1)+max_v(2)+max_v(3)
num_cheliang(3) = num_cheliang(3)+1;
capacity = max_vc(3);
car_path(k,1) = 3;
end
end
end
%% 遍历任务
for i = 1:length(renwu)
if upload_num(k)+demand(renwu(i)+1,2) <= capacity % 可以再装当前货物需求点的货
car_path(k,m+1) = renwu(i); % 设置第k条路径第m个货物需求点
if m == 1
dist_all(k) = D(1,renwu(i)+1); % 供货点到第一个需求点的距离
else
dist_all(k) = dist_all(k)+D(renwu(i-1)+1,renwu(i)+1); % 前一个需求到到当前需求点的距离
end
m = m+1; % 更换得到下一个需求点
upload_num(k) = upload_num(k)+demand(renwu(i)+1,2); % 当前路径上的货物转载量
else
% 需要更换路径了
dist_all(k) = dist_all(k) + D(1,renwu(i-1)+1);
j = j+1;
% 需要换一个路径
k = k+1;
% 这条路径已经供货的需求点数量清为1
m = 1;
% 装载量清0
upload_num(k) = 0;
% 确认是哪种运输工具
if cheliang(j)>N && cheliang(j) <= N+max_v(1)
num_cheliang(1) = num_cheliang(1)+1; % 使用第一种车辆的数目增加1
capacity = max_vc(1); % 这一条线路最大运输量
car_path(k,1) = 1; % 这条路径的运输工具类型
else if cheliang(j) >N+max_v(1) && cheliang(j) <= N+max_v(1)+max_v(2)
num_cheliang(2) = num_cheliang(2)+1;
capacity = max_vc(2);
car_path(k,1) = 2;
else if cheliang(j) > N+max_v(1)+max_v(2) && cheliang(j) <= N+max_v(1)+max_v(2)+max_v(3)
num_cheliang(3) = num_cheliang(3)+1;
capacity = max_vc(3);
car_path(k,1) = 3;
end
end
end
% 给超了的任务布置路径
car_path(k,m+1) = renwu(i); % 设置第k条路径第m个货物需求点
if m == 1
dist_all(k) = D(1,renwu(i)+1); % 供货点到第一个需求点的距离
else
dist_all(k) = dist_all(k)+D(renwu(i-1)+1,renwu(i)+1); % 前一个需求到到当前需求点的距离
end
m = m+1; % 更换得到下一个需求点
upload_num(k) = upload_num(k)+demand(renwu(i)+1,2); % 当前路径上的货物转载量
end
end
3.2.2 适应度函数
MATLAB代码(fitness2.m)
function [all_fitness,best_fitness,all_obj,best_obj,best_chrom,best_car_path] = fitness2(chrom,demand,N,max_vc,max_v,cb,cv,D)
% chrom为整个种群,N为任务需求点个数,max_vc为最大载重
% max_v为每一种交通工具的最大可用量,demand为每个点的需求量
% cb为固定出车成本,cv为单位运输成本
%% 参数设置
num = size(chrom,1); % 个体数
best_fitness = 0;% 种群最佳适应度
best_obj = 0;
best_chrom = []; % 最优解对应的个体
best_car_path = []; % 最优解对应的路径以及交通工具安排情况
all_fitness = zeros(1,num);
all_obj = zeros(1,num);
%% 正式计算适应度
for i = 1:num
% 解码得到路径及交通工具安排情况
[car_path,upload_num,dist_all] = decode(chrom(i,:),demand,N,max_vc,max_v,D);
% 计算个体对应的目标函数值
obj = 0;
% 遍历路径计算总费用
for k = 1:size(car_path,1) % 遍历每一条路径
% 加上固定出车成本费,car_path(k,1)为第k条路径的运输工具序号
obj = obj+cb(car_path(k,1));
% 加上这一条路径的运输费用
% k,length(dist_all)
obj = obj + cv(car_path(k,1))*upload_num(k)*dist_all(k);
end
% 计算适应度
fit = 1/obj;
all_fitness(i) = fit;
all_obj(i) = obj;
if fit > best_fitness % 优于目前最优个体则更新最优个体
best_fitness = fit;
best_obj = obj;
best_chrom = chrom(i,:);
best_car_path = car_path;
end
end
end
3.2.3 选择
还是采用轮盘赌，MATLAB代码(Select.m)
%% 选择操作
%输入
%Chrom 种群
%FitnV 适应度值
%GGAP：代沟
%输出
%SelCh 被选择的个体
function SelCh=Select(Chrom,FitnV,GGAP)
NIND=size(Chrom,1); % 种群大小
NSel=max(floor(NIND*GGAP+.5),2); % 被选择的个体数目
ChrIx=Sus(FitnV,NSel); % ChrIx为被选择个体的索引号
SelCh=Chrom(ChrIx,:); % 根据索引选择个体
个体选择代码(Sus.m)
% 输入:
%FitnV 个体的适应度值
%Nsel 被选择个体的数目
% 输出:
%NewChrIx 被选择个体的索引号
function NewChrIx = Sus(FitnV,Nsel)
[Nind,~] = size(FitnV); % Nind为个体总数，也即索引最大值
cumfit = cumsum(FitnV); % 累计求和适应度
trials = cumfit(Nind) / Nsel * (rand + (0:Nsel-1)');
Mf = cumfit(:, ones(1, Nsel));
Mt = trials(:, ones(1, Nind))';
[NewChrIx, ~] = find(Mt < Mf & [ zeros(1, Nsel); Mf(1:Nind-1, :) ] <= Mt);
[~, shuf] = sort(rand(Nsel, 1));
NewChrIx = NewChrIx(shuf);
3.2.4 交叉
MATLAB代码(Recombin.m)
%% 交叉操作
% 输入
%SelCh 被选择的个体
%Pc 交叉概率
%输出：
% SelCh 交叉后的个体
function SelCh=Recombin(SelCh,Pc)
NSel=size(SelCh,1); % NSel 为被选择的个体数量
for i=1:2:NSel-mod(NSel,2) % 临近两个个体交叉
if Pc>=rand %交叉概率Pc
[SelCh(i,:),SelCh(i+1,:)]=intercross(SelCh(i,:),SelCh(i+1,:));
end
end
%% 交叉函数
function [a,b]=intercross(a,b)
L=length(a);
cross_p=randperm(L-2,1)+1; %交叉位置
% 以下为交叉步骤
c1=a(1:cross_p);
a(1:cross_p)=b(1:cross_p);
b(1:cross_p)=c1;
c2=a(cross_p+1:length(a));
a(cross_p+1:length(a))=b(cross_p+1:length(b));
b(cross_p+1:length(b))=c2;
3.2.5 变异
MATLAB代码(Mutate.m)：
%% 变异操作
%输入：
%SelCh 被选择的个体
%Pm 变异概率
%输出：
% SelCh 变异后的个体
%%
function SelCh=Mutate(SelCh,Pm)
[NSel,L]=size(SelCh);
for i=1:NSel % 遍历被选择的个体，逐次查看是否需要变异
if Pm>=rand
R=randperm(L);
SelCh(i,R(1:2))=SelCh(i,R(2:-1:1));
end
end
3.2.6 其他部分
- 主函数(
main.m)
clear,clc
load zuobiao_X %载入节点坐标:【节点编号，X,Y】（包含始节点）
load demand %载入需求矩阵：【节点编号，需求量】（包含始节点）
%% ---------------------运输工具灵敏度-----------------%%
num_truck_list = 0:6;
num_train_list = 0:4;
num_plane_list = 0:4;
num_truck = 6;
num_train = 3;
num_plane = 2;
% 货车数量改变
for i = 1:length(num_truck_list)
[best_truck(i)] = GA_main(num_truck_list(i),num_train,num_plane,demand,X);
end
% 火车数量改变
for i = 1:length(num_train_list)
[best_train(i)] = GA_main(num_truck,num_train_list(i),num_plane,demand,X);
end
% 飞机数量改变
for i = 1:length(num_plane_list)
[best_plane(i)] = GA_main(num_truck,num_train,num_plane_list(i),demand,X);
end
%% 绘制灵敏度分析图
% 货车
figure
subplot(1,3,1)
Q1 = plot(num_truck_list,best_truck);
set(Q1,'linewidth',1.1);
xlabel('货车数量','fontweight','bold');
ylabel('费用(单位：元)','fontweight','bold');
grid on
set(gca,'linewidth',1.1) % 设置坐标轴句柄属性
% 火车
subplot(1,3,2)
Q2 = plot(num_train_list,best_train);
set(Q2,'linewidth',1.1);
xlabel('火车数量','fontweight','bold');
ylabel('费用(单位：元)','fontweight','bold');
grid on
set(gca,'linewidth',1.1) % 设置坐标轴句柄属性
% 飞机
subplot(1,3,3)
Q3 = plot(num_plane_list,best_plane);
set(Q3,'linewidth',1.1);
xlabel('飞机数量','fontweight','bold');
ylabel('费用(单位：元)','fontweight','bold');
grid on
set(gca,'linewidth',1.1) % 设置坐标轴句柄属性
%% -----------------------需求量-----------------%%
% 选取6个需求点做灵敏度分析
demand_dot = [2 6 9 14 19 26]; % 一定要是偶数个
demand_range = 1:5;
demand_new = demand;
for i = 1:length(demand_dot)
for j = 1:length(demand_range)
demand_new(demand_dot(i)+1,2) = demand_range(j); % 修改新需求
[best_demad{i}(j)] = GA_main(num_truck,num_train,num_plane,demand_new,X);
end
end
% 绘图
figure
for i = 1:length(demand_dot)
subplot(2,length(demand_dot)/2,i)
Q = plot(demand_range,best_demad{i});
set(Q,'linewidth',1.1);
xlabel('需求量','fontweight','bold');
ylabel('费用(单位：元)','fontweight','bold');
grid on
set(gca,'linewidth',1.1) % 设置坐标轴句柄属性
end
- 遗传算法部分主函数(
GA_main.m)
function [best_all_and_all] = GA_main(num_truck,num_train,num_plane,demand,X)
%遗传算法VRP问题求解
%% 变量说明和场景描述
%% 输入数据：
% D：节点距离矩阵
% NIND： 种群规模
% X：节点坐标
% MAXGEN：最大迭代次数
% max_vc:最大车辆容量
% Demand:节点需求矩阵
% max_d:车辆运输最大距离
% v:车辆速度
% cv:车辆运输单价
% max_v:可供调配的最大车辆数
%% 优化目标：运输成本最小 ，车辆数最少
%% 约束条件：满足车辆容量限制，满足车辆运输距离限制，同一节点只接收一辆车服务,
%% 单个节点位置的需求小于车辆的最大载重
% clear
% clc
%% 初始化
% step1 加载车辆和服务信息数据
% load zuobiao_X %载入节点坐标:【节点编号，X,Y】（包含始节点）
% 计算结点距离矩阵
for i = 1:size(X,1)
for j = 1:size(X,1)
D(i,j) = sqrt((X(i,2)-X(j,2))^2+(X(i,3)-X(j,3))^2);
end
end
N=size(X,1)-1; %配送节点个数
gj_num = 3; % 运输工具类型数量
% load demand %载入需求矩阵：【节点编号，需求量】（包含始节点）
max_vc = [50,27,13]; %运输工具最大载重
cv = [30,20,10]; %单位运输成本
cb = [200,150,100]; % 固定出车成本
%max_d=40; %单车行驶最大距离
% max_v=[6,3,2]; %假设可调配车辆数的最大值
max_v = [num_truck,num_train,num_plane];
best_obj_list = [];
% step2 GA参数初始化
NIND=150; %种群大小
MAXGEN=2000; %最大迭代次数
Pc=0.9; %交叉概率
Pm=0.05; %变异概率
GGAP=0.9; %代沟
%% 开始迭代
% step1 舒适化种群
chrom=initpop(NIND,N+sum(max_v));
gen=0;
% % 迭代过程图
% figure
% hold on
% box on;
% xlim([0,MAXGEN])% x坐标轴上下限 x limit
% title('优化过程')
% xlabel('代数')
% ylabel('最优值')
% 计算第一代个体的适应度，并找出最优个体及对应的解与染色体
[all_fitness,best_fitness,all_obj,best_obj,best_chrom,best_car_path] = fitness2(chrom,demand,N,max_vc,max_v,cb,cv,D);
best_obj_list = [best_obj_list best_obj];
pre_obj = best_obj; % 上一代的最优目标函数值
for gen = 1:MAXGEN
% gen
%% 选择
FitnV = all_fitness;
SelCh=Select(chrom,FitnV,GGAP);
%% 交叉操作
SelCh=Recombin(SelCh,Pc);
%% 变异
SelCh=Mutate(SelCh,Pm);
%% 重插入子代的新种群
ObjV = all_obj;
chrom=Reins(chrom,SelCh,ObjV); % 合并成新种群
%% 计算新种群的适应度等
[all_fitness,best_fitness,all_obj,best_obj,best_chrom,best_car_path] = fitness2(chrom,demand,N,max_vc,max_v,cb,cv,D);
best_obj_list = [best_obj_list best_obj];
% %% 绘图
% line([gen-1,gen],[pre_obj,best_obj]);
% pre_obj = best_obj; % 上一代的最优目标函数值
end
best_all_and_all = best_obj_list(end);
% %% 结果分析
% % 最优目标函数值
% disp(['最优目标函数值:',num2str(best_obj_list(end))])
% % 最优时所用的交通工具数量（即路径数量）
% car_num = size(best_car_path,1);
% disp(['最优时总的交通工具使用量为：',num2str(car_num)])
% %% 绘制路径图
% lujing = cell(1,car_num);
% % 找出路径
% for i = 1:car_num
% lujing{i}(1,1:2) = X(1,2:3);
% store_num = length(best_car_path(i,find(best_car_path(i,:) ~= 0)))-1; % 需求点数
% for j = 1:store_num
% lujing{i}(j+1,1:2) = X(best_car_path(i,j+1)+1,2:3);
% end
%
% end
% % 绘图
% figure
% for i = 1:car_num
% plot(lujing{i}(:,1),lujing{i}(:,2),'-*','linewidth',1.1);
% hold on
% end
% xlabel('横坐标','fontweight','bold')
% ylabel('纵坐标','fontweight','bold')
% title('路径规划图','fontweight','bold')
% axis([-5 max(X(:,2))*1.1 0 max(X(:,3))*1.1]) % 设置坐标轴范围
% text(X(:,2),X(:,3),num2str(X(:,1))) % 标注文字
%
% trans_type = ['货车';'火车';'飞机'];
% legend_lujing = [];
% for i = 1:car_num
% lujing{i} = [lujing{i};X(1,2:3)];
% lujing_name = ['路径' num2str(i) '(' trans_type(best_car_path(i,1),:) ')'];
% legend_lujing = [legend_lujing;lujing_name ];
% end
% legend(legend_lujing)
% set(gca,'linewidth',1.1,'gridlinestyle','--')
% grid on
- 初始化(
initpop.m)
%% 初始化种群
%输入：
% NIND：种群大小
% N： 个体染色体长度（这里为城市的个数）
%输出：
%初始种群
function chrom=initpop(NIND,N)
chrom=zeros(NIND,N);%用于存储种群
for i=1:NIND
chrom(i,:)=randperm(N);%随机生成初始种群
end
- 精英个体的选择(
Reins.m)
%% 重插入子代的新种群
%输入：
%Chrom 父代的种群
%SelCh 子代种群
%ObjV 父代适应度
%输出
% Chrom 组合父代与子代后得到的新种群
function chrom=Reins(Chrom,SelCh,ObjV)
NIND=size(Chrom,1); % 原始种群大小
NSel=size(SelCh,1); % 被选择的种群大小
[~,index]=sort(ObjV);
% 将初始种群里面最优秀的NIND-NSel个个体放入新种群，再加上选择得到的个体合并成一个种群
chrom=[Chrom(index(1:NIND-NSel),:);SelCh];