带时间窗物流配送车辆路径问题

带时间窗物流配送车辆路径问题

摘要

本题是一个带有时间窗的车辆路径安排问题（VRPTW问题）。根据题目条件，本文建立了一个求解最小派送费用的VRPTW优化模型，采用遗传算法，给出了该模型的求解方法。然后，对一个实际问题进行求解，给出了一个比较好的路线安排方式。

模型一（见5.1.2）针对问题一，在需求量、接货时间段、各种费用消耗已知的情况下，决定采用规划模型，引入0-1变量，建立各个约束条件，包括车辆的容量限制，到达每个客户的车辆和离开每个客户的车辆均为1的限制，总车辆数的限制，目标函数为费用的最小化，费用包括车辆的行驶费用，车辆早到或晚到造成的损失。

模型一的求解采用遗传算法（见5.1.3），对题目给出的实际问题进行求解，

首先按照需求期望根据模型一得到一个比较好的方案，然后按照这一方案进行送货，在送货过程中，如果出现需求量过大的情况，允许车辆返回仓库进行补充。

模型一的思路清晰，考虑条件全面。但最优解解决起来困难，遗传算法只是一种相对好的解决方法，可以找出最优解的近似解。模型二的想法比较合理，易于实施，但还有待改进。

关键词：规划 时间窗 物流 车辆路径 遗传算法

一、 问题重述

一个中心仓库，拥有一定数量容量为Q 的车辆，负责对N 个客户进行货物派送工作，客户i 的货物需求量为i q ，且i q Q <，车辆必须在一定的时间范围[],i i a b 内到达，早于i a 到达将产生等待损失，迟于i b 到达将处以一定的惩罚，请解决如下问题：

（1）给出使派送费用最小的车辆行驶路径问题的数学模型及其求解算法。并具体求解以下算例：

客户总数N=8,每辆车的容量Q=8（吨/辆）, 各项任务的货运量i q （单位：吨）、装货（或卸货）时间i s （单位：小时）以及要求每项任务开始执行的时间范围[],i i a b 由附录1给出，车场0与各任务点以及各任务点间的距离（单位：公里）由附件二给出，这里假设车辆的行驶时间与距离成正比，每辆车的平均行驶速度为50公里/小时，问如何安排车辆的行驶路线使总运行距离最短； （2）进一步请讨论当客户i 的货物需求量i q 为随机参数时的数学模型及处理方法。

二、 问题分析

本题主要在两种不同情况下，研究使派送费用最小的车辆行驶路径问题。车辆行驶派送的费用主要包括运输成本、车辆在客户要求到达时间之前到达产生的等待损失和车辆在客户要求到达时间之后到达所受惩罚等等。为满足派送费用最小的需求，即要使所选行车路径产生的总费用最小，从而确定出最佳的车辆派送方案。

当客户i 的货物需求量i q 固定时，首先，我们根据题意，取若干辆车进行送货，然后，主要考虑每辆车各负责哪些客户的送货任务，我们可以给出满足题中限制条件的很多参考方案供选用，并考虑以所选行车路径产生的总费用最小为目标的情况下，建立最优化模型确定最佳的车辆派送方案。

进一步讨论，当客户i 的货物需求量i q 为随机参数时，我们首先可以简化随机模型，根据客户i 的货物需求量的期望与方差，确定每天应该运送给客户i 的货物量，即i q ，再根据第一题，确定最佳的车辆派送方案。

但考虑到客户的储存能力有限及货物在客户处的储存费用，客户不需要将一天的货物一次性接收完，只要满足缺货的情况出现的概率很低，客户可以让配送中心一天几次送货，这样可以得到很多满足约束的方案，考虑以单位时间的储存

费用最小为目标，建立最优化模型，确定配送中心给每位客户每次的配送量、配送周期与最有车辆行驶路径。

三、模型假设

（1）每个客户的需求只能由一辆配送车满足；

（2）每辆车送货时行驶的路程不超过它所能行驶的最远路程；

（3）中心仓库的车辆总数大于或等于当派送费用最小时所需的车辆数；

（4）从配送中心到各个用户、各个用户之间的运输距离已知；

（5）配送中心有足够的资源以供配送。

四、符号说明

五、模型的建立和求解

5.1 问题一模型的建立及求解：

5.1.1问题的分析

中心仓库为了给N个客户派送货物，供发出m辆车，为了派货的节约和方便，每辆车载着适量的货物出发，可以给某一片的若干个满足约束条件的客户派送货物，见图一：

图一中心仓库派送货物图

中心仓如上图库派送货物时，必须满足约束条件：

（1）各个客户群的总需求小于或等于运输车的装载量；

（2）每个客户都必须且只能由一辆运输车运输所需货物；

（3）运输车为每位客户开始服务的时间必须尽可能在时间窗内。

根据如上的约束条件，我们可以得到很多可行解，但考虑到以所选行车路径产生的总费用最小为目标的情况下，我们可以建立最优化模型确定最佳的车辆派送方案，最优路径产生图如下：

图二 最优路径产生图

5.1.2模型的建立

（1）中心仓库使用车辆数量的确定 设配送中心需要向N 个客户送货，每个客户的货物需求量是gi （i ＝1，2，…..N ），每辆配送车的载重量是Q ，且gi

1/N

i i m g aQ ==∑

[ ]表示取整，a 为参数，0

（2）引入0—1变量：

1)ijs

x 表示车辆s 是否从客户i 行驶到客户j 。定义其为0—1变量，则

?

??=01

x ijs

否则行驶到客户从客户车辆j i s 1,,s m =

2)is y 表示客户i 的任务由车辆s 完成。同样定义其为0—1变量，则

???=0

1

is y

否则完成的任务由车辆客户s i 1,,s m =

（3） 非线性规划模型的建立： a ．目标函数的确定。

题目要求所选行车路径产生的总费用最小，我们确定总费用为目标函数，记为Z 。

总费用由运输成本A 、等待损失B 和迟到所收惩罚C 组成，根据题意有：

11

1

N

N m

ij ijs

i j s A c c x

====∑

∑∑

1*max{,0}N

i i B d D ==∑

1

*max{,0}N

i i C e X ==∑

所以，总费用Z 最小化为：

1

1

1

min *max{,0}*max{,0}N

N m

N N

ij ijs

i i i j s i i Z c c x

d D

e X ======++∑

∑∑∑∑

b ．约束条件的确定。 约束1：

车辆k 的运送总重量应不超过车辆的最大载重，即车辆有一定的运送能力，则可引入约束条件，

1

N

i is

i q y

Q =≤∑ (},,2,1{m k ∈?)

约束2：

每个客户只能由一辆车来配送，则可引入约束条件，

11

m

is s y m

=?=??∑ 1,2,3,...0i N i == 约束3：

保证到达一个客户的车辆也离开该客户，则可引入约束条件，

11

1m N

ijs

s i x

===∑∑ (1,2,3,

,;j N =)

11

1m N

ijk

s j x

===∑∑ (1,2,3,

,;i N =)

c ．变量之间关系的确定

由上可确定等待时间i D ，超时时间i X 为：

i i i

i i i

D a t X t b =-=- 1,2,....1,2,....i N

i N ==

车辆k 从客户i 到客户j 需经过两段时间ij t 为：

/ij ij t c v = ,1,2,....i j N =

设车辆为客户i 运送完货物后即为客户j 运送，则到达客户i 处时间i t 和到达客户j 处时间j t 之间的关系为：

11(max{,0})N

m

j ijs i i ij i i s t x t D t s ===+++∑∑

d ．此非线性规划模型为：

1

1

1

min *max{,0}*max{,0}N

N m

N N

ij ijs

i i i j s i i Z c c x

d D

e X ======++∑

∑∑∑∑

..t s

1

N

i is

i q y

Q =≤∑ m s ,,2,1 =

11

m

is

s y m

=?=??∑ 1,2,3,...0i N i == 11

1m N

ijs

s i x

===∑∑ 1,2,3,

,;j N =

11

1m

N ijk

s j x

===∑∑ 1,2,3,

,;i N =

i i i

i i i

D a t X t b =-=- 1,2,....1,2,....i N

i N ==

/ij ij t c v = ,1,2,....i j N =

11

(max{,0})N m

j ijs i i ij i i s t x t D t s ===+++∑∑

5.1.3模型的求解

我们采用遗传算法解决上面的问题： 1．编码

采用自然数编码，即序数编码。货物运输路线可以编成长度为N+m 的染色体

11121s 21210,,,,0,,,0,,0,,,)t m mw i i i i i i i (，，其中，ik i 表示第ik i 项任务。0

表示车场，m 表示完成任务所需的车辆数。 2.出生初始群体

初始群体随机产生，即产生N 项货物运输任务点的全排列，如12,,,N i i i ，

如果

1

1

s ij

j q

Q -=≤∑，且1

s

ij j q Q =>∑，将s 至N 的数向后移动一位，将0插入第s 位。

接着，继续上述操作，直到m 个0全部插入为止。这样就构成了一条初始染色体。用这种方法构造一个群体的染色体。如：82576314，该编码插零之后变成

0825*******。它代表着需要三辆车运输货物。其中，第1辆车行走路线为

08250，即从仓库出发到依次到8、2、5商店再回到仓库。第2辆车行走路线

为07630，第3辆车行走路线为0140。 3.适应度函数

适应度函数取'

k k

bz f z =

，其中k f 为染色体k v 的适应度，b 为常数，'z 为初始种群中最好的染色体的运输成本，k z 为染色体k v 对应的运输成本。

4.遗传算子

选取最佳保留的轮盘赌复制法进行染色体的复制。变异算子采用反转变异。交叉算子用最大保留交叉，其操作过程为：

a) 若染色体交叉点处的两个基因都为0，则直接进行顺序交叉运算； b) 若染色体交叉点处的基因不全为0，则将交叉点左移（右移），直到左右

两个交叉处的基因都为0，再进行顺序交叉运算。

5.算法的实现步骤

Step1：采用自然数编码的方式，构造表示可行车路线的染色体；

Step2：设置控制参数，包括交叉率0.7c p =、变异率0.1m p =、群体规模10n =； Step3：初始化，令0d =，随机产生初始群体(0)p ，群体中包括n 个染色体，

每个染色体代表一条行车线路；

Step4：令1i =；

Step5：将群体()p d 中的第i 个染色体译为线路长度； Step6：计算适应度；

Step7：若满足算法终止条件，则停止，否则继续； Step8：1i

i =+；

Step9：若i n ≤，回到step5，否则，转step10；

Step11：进行最大保留交叉、基于位的变异和倒位操作； Step12: 1d

d =+；

Step13：若满足算法终止条件，则停止，否则转step4。

运用matlab 软件编写程序得到在车辆总行程最短的条件小的派送方案为：

5.1.4结果分析

从上面解出的派送方案可以看出，上面的每条路线在车辆送完货物后，直接从最后一个客户处返回中心仓库，我们通过求中心仓库到各个客户的最短路径可以发现，上面的返回路线不是最优的，返回路线可以经过某些客户使得所走路线最短。

公里。

5.2 问题二模型的建立及求解：

5.2.1问题的分析与假设

在问题一中，根据已知的各个商店的需求分布，根据遗传算法求解，预先确定一条总费用最小的路径（或者虽不是最优路径，但是此结果能够接受的）。车辆沿该路径服务商店，因此服务商店的顺序固定不变。

而问题二中，在车辆服务商店的过程中，事先并不知道商店的需求量。而这个不确定信息要随着车辆的服务逐步确定，商店具体的需求只有车辆到达用户后才确定。这样第一问的求解方法得到的路径并不适用于第二问的求解。 假设：

（1）商店的需求变化符合正态分布，第i 个商店的需求量的期望和方差分别为i μ和2i σ。

（2）商店的供货可以分为多次补充但在每次供货中最大程度满足用户需求，即只有出现当前车载余量小于用户需求量时才出现下一次的供货。 5.2.2模型的建立

基于第一问解决了在已知用户需求概率情况下，确定一个服务方案，满足所有n 个商店的需求并且使车辆期望行程费用最小这个问题。我们由假设可知，第

i 个商店的需求量的期望为i μ，则由需求量的期望得到一个使车辆期望行程费用

最小的服务方案，称该方案为A 。

当用户的需求未知（当车辆服务商店时需求变成已知量）时，由于用户的需求量在区间(3)i i μσ±的概率是96%，而在区间外的事件可以看成小概率事件，由小概率事件定理可知，在一次试验中，小概率事件可以看成不可能事件。由此可知，用户需求量就在区间(3)i i μσ±里。

用户需求在区间(3)i i μσ±里，而需求决定服务方案。由上面可知，服务方案在A 方案附近变化，而变化的幅度由方差2i σ决定。当2i σ越小时（说明需求量接近一个常数期望i μ），最优方案（或满意方案）与A 方案越接近（即在A 方案上面稍作改动即可）。反之，A 方案需作较大的方案。

由假设（2）知，车辆只要满足当前用户部分需求，就服务该用户，用户未满足的部分以后再服务。在服务用户后，车辆根据当前的位置信息、车载余量以及需要服务用户的信息，决定下一个服务用户(包括当前还未满足需求的用户)或回库房装载货物。

先按A 方案进行分配车辆路线（若增加车辆数目虽可以更好满足条件，但会增加多于开支）。假设m 辆车都从仓库0出发，按A 方案中的预定路线进行服务，当服务完第一个商店时，再判断剩余的货物量。此时货物量为'i Q Q q =-（其中Q 为货车满载量，i q 为车辆到达商店i 时确定了i 个商店的需求量）。然后判断该剩余量是否在大于下一个商店需求量(3)i i μσ+，若大于则进行下一个商店服

务，若小于则判断下个商店是否满足大于(3)i i μσ+这个条件，若满足，则对其进行服务，否则继续下个判断，直至其所有负责区域遍历完一遍。

当判断其负责区域所有的商店不满足条件时再判断该剩余量是否大于下一个商店需求量(3)i i μσ-，若大于则进行下一个商店服务，若小于则判断下个商店是否满足大于(3)i i μσ-这个条件，若满足，则对其进行服务，否则继续下个判断，直至其所有负责区域遍历完一遍。

若所有商店的需求量(3)i i μσ-大于车辆货物剩余量，则说明此车辆的剩余量不能满足其所负责的区域，因此该车辆需要回仓库进行货物补充。当货物补充完之后进行判断剩余未服务商店的时间窗口和路程距离进行判断（产生方法同于第二问A 方案的产生方法），然后再进行服务。服务商店之后再进行前面的判断，直至其所有负责商店都被服务完回到仓库。

六、 模型的评价和推广

6.1 模型的评价

由5.1和5.2建立的模型，我们得到了有软时间限制的物流配送车辆路径问题的解决办法。首先我们对问题进行了合理的假设，模型简化了实际中的复杂问题，考虑了主要的约束条件与目标，思路清晰，结果也基本符合实际。

但是，模型对实际进行简化的同时忽略了实际中很多次要的条件，由累加效果可能会产生很大误差，同时，我们进行模型的求解时只是简单使用了遗传算法，没有对其进行改进。

6.2 模型的推广

1．模型中不合实际的假设：

（1） 在问题一和问题二总费用的组成上，我们没有考虑组织送货的费用，即每组织一辆车进行送货都需要一定的费用，我们设为S 元/（次、辆）； （2） 在问题一中，我们假设每辆车送货时行驶的路程不超过它所能行驶的最远路程，但是实际上，考虑到车辆行驶中的耗油等因素，每辆车的行驶最大距离都有限制，我们设车辆行驶的最大距离为L ，我们可以得到约束条件：

00

N N

ijs ij

i j x

c L ==≤∑∑ 1,2,...,s m =

（3） 在实际中，为了公平，每位送货车辆司机行驶的路程相差必须控制在一定范围之内，我们取这个差值的最大允许值为M ，则有：

00

00

N N

N N

ijs ij

ijk ij i j i j x

c x c M ====-≤∑∑∑∑ ,1,2,...,s k m =

2. 模型的推广

目标函数：

1

1

1

min *max{,0}*max{,0}N

N m

N N

ij ijs

i i i j s i i Z c c x

d D

e X mS

======+++∑

∑∑∑∑

1

N

i is

i q y

Q =≤∑ m s ,,2,1 =

11m

is s y m

=?=??∑ 1,2,3,...0i N

i == 11

1m N

ijs

s i x

===∑∑ 1,2,3,,;j N =

11

1m N

ijk

s j x

===∑∑ 1,2,3,,;i N =

约束条件

i i i D a t =- 1,2,....i N = i i i X t b =- 1,2,....i N =

/ij ij t c v = ,1,2,....i j N =

11(max{,0})N

m

j ijs i i ij i i s t x t D t s ===+++∑∑

00N N

ijs ij

i j x

c L ==≤∑∑ 1,2,...,s m = 00

00N

N N

N

ijs ij

ijk ij i j i j x

c x c M ====-≤∑∑∑∑ ,1,2,...,s k m =

根据以上目标函数与约束条件，带入实际数据，由遗传算法即可求解出总费

用最小的车辆行驶路径方案。

七、参考文献

[1]李军，谢秉磊，郭耀煌，非满载车辆调度问题的遗传算法。系统工程理论与实践，2000，20（3）：235—239；

[2]邢文训，谢金星编著现代优化计算方法。北京：清华大学出版社，1999.140；

[3]魏明高成修胡润洲，一种带时间窗和容量约束的车辆路线问题及其Tabu Search 算法，运筹与管理，Vol. 11 ,No. 3：P49-P53，2002；

[4]胡大伟朱志强胡勇，车辆路径问题的模拟退火算法，中国公路学报，Vol.1 9 No.4：P123-P126,2006

[5]阎庆鲍远律，新型遗传模拟退火算法求解物流配送路径问题，https://www.360docs.net/doc/8411631075.html,/grid2008/detail.aspx?filename=HNSZ200804014&dbname=CJF

D2008，2009/8/16

[6]张志霞邵必林，基于改进蚁群算法的运输调度规划，公路交通科技，Vo1．25 No．4：P137-P140，2008；

[7]杜文衰庆达周再玲，一类随机库存/运输联合优化问题求解过程分析，中国公路学报，Vol. 17 No．1：P114-P118，2004.

八、附录

9.1附录一

表1 任务的特征及其要求

表2 点对之间的距离

9.2附录二

Matlab 7.0 程序：

function gatsp(s)

sum=0;

M=1000000; %无穷大数

inn=100;

citynum=8;

K=2;

gnmax=1000; %最大代数

pm=0.8; %变异概率

pc=0.2;

c=[0,40,60,75,90,200,100,160,80;40,0,65,40,100,50,75,110,100; 60,65,0,75,100,100,75,75,75;75,40,75,0,100,50,90,90,150; 90,100,100,100,0,100,75,75,100;200,50,100,50,100,0,70,90,75;

100,75,75,90,75,70,0,70,100;160,110,75,90,75,90,70,0,100; 80,100,75,150,100,75,100,100,0];

q=[2 1.5 4.5 3 1.5 4 2.5 3];

s=[1 2 1 3 2 2.5 3 0.8];

a=[1 4 1 4 3 2 5 1];

b=[4 6 2 7 5.5 5 8 4];

%-------------------------------------------------------------------------- %产生初始种群

m=zeros(1,inn);

m=m';

KM=citynum+K-1;

s=zeros(inn,citynum+K-1);

for i=1:1:inn

s(i,:)=randperm(KM);

end

s=[m s];

for i=1:inn

for j=1:KM-1

if s(i,j)>citynum

s(i,j)=0;

end

end

end

%-------------------------------------------------------------------------- %主程序

function ga

[f,p]=objf(s)

gn=1;

while gn

for j=1:2:inn

seln=sell(s,ps); %选择操作

scro=cross(s,seln,pc); %交叉操作

scnew(j,:)=scross(1,:);

scnew(j+1,:)=scross(2,:);

smnew(j,:)=chang(scnew(j,:),pm); %变异操作

smnew(j+1,:)=chang(scnew(j+1,:),pm);

end

s=smnew; %产生了新的种群

[f,p]=objf(s,dislist); %计算新种群的适应度

%记录当前代最好和平均的适应度

[fmax,nmax]=max(f);

ymean(gn)=1/mean(f);

ymax(gn)=1/fmax;

%记录当前代的最佳个体

x=s(nmax,:);

gn=gn+1;

%pause;

end

gn=gn-1;

figure(2);

plot(ymax,'r'); hold on;

plot(ymean,'b');grid;

title('搜索过程');

legend('最优解','平均解');

function pcc=pro(pc);

test(1:100)=0;

l=round(100*pc);

test(1:l)=1;

n=round(rand*99)+1;

pcc=test(n);

%-------------------------------------------------------------------------- %计算适应度

function [f,p]=objf(s);

y=zeros(citynum+1,citynum+1);

for i=1:inn-1

a=s(i,:);

for j=1:KM-1

m=a(j);

n=a(j+1);

m=m+1;

n=n+1;

end

y(m,n)=1;

y=y';

for i=1:citynum

for j=1:citynum

mubiaob=c(i,j)*y(i,:);

end

end

xuq1=0;

for i=1:citynum

for j=1:citynum

xuq1=xuq1+s(i)*y(i,:)-q(i);

end

xuqiu=max((xuq1),0)*M;

end

end

shij1=0;

shij2=0;

for i=1:citynum

for j=1:citynum

for l=1:citynum

shij1=shij1+t(i)-a(i);

shij2=shij12+b(i)-t(i);

end

shij3=max((shij1),0);

shij4=max((shij2),0);

shijian=M*shij3+M*shij4;

end

end

f=mubiao+xuqiu+shijian;

f=1/f;

end

%-------------------------------------------------------------------------- %计算选择概率

fsum=0;

for i=1:inn

fsum=fsum+f(i);

end

for i=1:inn

ps(i)=f(i)/fsum;

end

%计算累积概率

p(1)=ps(1);

for i=2:inn

p(i)=p(i-1)+ps(i);

end

p=p';

p

%-------------------------------------------------------------------------- %“选择”操作

%从种群中选择两个个体

function seln=sell(s,p)

inn=size(p,1);

for i=1:2

r=rand; %产生一个随机数

prand=p-r;

j=1;

while prand(j)<0

j=j+1;

end

seln(i)=j; %选中个体的序号

end

sel1=seln(1);

sel2=seln(2);

%-------------------------------------------------------------------------- %“交叉”操作

function snew=cross(A,B,pc)

A=s(sel1,:);

B=s(sel2,:);

c=find(A==0);

d=find(B==0);

a=sym(A);

b=sym(B);

k=size(a,2);

for i=1:size(a,2)

for j=1:k

e(i,c(k))=a(i+k-1);

end

end

for i=1:size(a,2)

for j=1:k

e(i,d(k))=b(i+k-1);

end

end

c0=round(rand*(k-1))+1;

c1=round(rand*(k-1))+1;

a=[f(:,c0),a];

b=[e(:,c1),b];

for i=1:size(a,2);

j=1:size(e,2)

if a(i)==e(j)

a(i)==[];

end

end

for i=1:size(b,2);

j=1:size(f,2)

if b(i)==f(j)

b(i)==[];

end

end

a=double(a);

b=double(b);

g=zeros(size(A));

for i=1:size(a)

for j=1:size(c)

g(i+j)=a(i);

end

end

h=zeros(size(A));

for i=1:size(b)

for j=1:size(d)

h(i+j)=b(i);

end

end

g=g';

h=h';

snew=[g h];

%--------------------------------------------------------------------------

%变异

function snew=chang(snew,pm)

bn=size(snew,2);

snnew=snew;

c2=round(rand*(bn-2))+1; %在[1,bn-1]范围内随机产生一个变异位

c3=round(rand*(bn-2))+1;

chb1=min(c2,c2);

chb2=max(c3,c2);

x=snew(chb1+1:chb2);

snnew(chb1+1:chb2)=fliplr(x);

pmm=pro(pm); %根据变异概率决定是否进行变异操作，1则是，0则否if pmm==1

c2=round(rand*(bn-2))+1; %在[1,bn-1]范围内随机产生一个变异位c3=round(rand*(bn-2))+1;

chb1=min(c2,c3);

chb2=max(c2,c3);

x=snew(chb1+1:chb2);

snnew(chb1+1:chb2)=fliplr(x);

end