VRPTW建模与求解—基于粒子群算法_vrptw数学模型

VRPTW建模与求解—基于粒子群算法

1 VRPTW简要描述

VRPTW（Vehicle Routing Problem with Time Windows）是指在经典VRP的前提上，给每个客户增添时间窗约束以进行运输服务，时间窗是一个顾客与厂商根据双方需求协商好的特定时间段，包含最早可到达时间和最晚必须到达时间。带时间窗约束使得 VRP 更加复杂性，通常可把时间窗分为硬时间窗（必须满足）、软时间窗（可以不满足，但会受到惩罚）。

2 课题场景设计

2.1 场景

单向：纯取货/纯送货；
单配送中心：只有一个配送中心/车场；
单车型：只考虑一种车型，
需求不可拆分：客户需求只能有一辆车满足；
车辆封闭：完成配送任务的车辆需回到配送中心；
车辆充足：不限制车辆数量，即配送车辆需求均能满足；
非满载：任意客户点的需求量小于车辆最大载重；
软时间窗：[最早可服务时间ET，最晚可服务时间LT]，早于ET会产生等待成本，晚于LT会有惩罚成本；

2.2 要求

优化目标：最小化车辆启动成本、车辆行驶成本、等待成本和惩罚成本之和；
约束条件：车辆行驶距离约束，重量约束，时间窗约束；
已知信息：配送中心位置、客户点位置、客户点需求、客户服务时间窗要求，客户服务时间，车辆最大载重、车辆最大行驶距离、车辆行驶速度、车辆启动成本、车辆单位距离行驶成本，等待成本，惩罚成本；

3 数学模型

3.1 符号定义：

3.2 数学模型

4 粒子群算法设计

4.1 算法设计

见【CVRP建模与求解-基于粒子群算法（python实现)】，算法设计一致。唯一区别在于在计算适应度时，VRPTW需要将车辆等待成本和惩罚成本记录，这两项惩罚成本的计算方法将见数学模型约束 f.惩罚成本。

4.2 python程序设计

# -*- coding: utf-8 -*-
import math
import random
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']  # 添加这条可以让图形显示中文


def calDistance(CityCoordinates):
    ''' 
    计算城市间距离
    输入：CityCoordinates-城市坐标；
    输出：城市间距离矩阵-dis_matrix
    '''
    dis_matrix = pd.DataFrame(data=None,columns=range(len(CityCoordinates)),index=range(len(CityCoordinates)))
    for i in range(len(CityCoordinates)):
        xi,yi = CityCoordinates[i][0],CityCoordinates[i][1]
        for j in range(len(CityCoordinates)):
            xj,yj = CityCoordinates[j][0],CityCoordinates[j][1]
            dis_matrix.iloc[i,j] = round(math.sqrt((xi-xj)**2+(yi-yj)**2),2)
    return dis_matrix


def greedy(CityCoordinates,dis_matrix):
    '''
    贪婪策略构造初始解
    输入：CityCoordinates-节点坐标,dis_matrix-距离矩阵
    输出：初始解-line
    '''
    #修改dis_matrix以适应求解需要
    dis_matrix = dis_matrix.astype('float64')
    for i in range(len(CityCoordinates)):dis_matrix.loc[i,i]=math.pow(10,10)
    dis_matrix.loc[:,0]=math.pow(10,10)#0不在编码内
    line = []#初始化
    now_city = random.randint(1,len(CityCoordinates)-1)#随机生成出发城市
    line.append(now_city)#添加当前城市到路径
    dis_matrix.loc[:,now_city] = math.pow(10,10)#更新距离矩阵，已经过城市不再被取出
    for i in range(1,len(CityCoordinates)-1):
        next_city = dis_matrix.loc[now_city,:].idxmin()#距离最近的城市
        line.append(next_city)#添加进路径
        dis_matrix.loc[:,next_city] = math.pow(10,10)#更新距离矩阵
        now_city = next_city#更新当前城市
    return line


def calFitness(birdPop,Demand,dis_matrix,CAPACITY,DISTABCE,C0,C1,C2,C3,time,V):
    '''
    贪婪策略分配车辆（解码），计算路径距离（评价函数）
    输入：birdPop-路径，Demand-客户需求,dis_matrix-城市间距离矩阵，CAPACITY-车辆最大载重,DISTABCE-车辆最大行驶距离,C0-车辆启动成本,
    C1-车辆单位距离行驶成本,C2-等待成本，C3-惩罚成本，time-服务时间窗和服务时间；
    输出：birdPop_car-分车后路径,fits-适应度
    '''
    birdPop_car,fits = [],[]#初始化
    for i in range(len(birdPop)):
        bird = birdPop[i]
        lines = []#存储线路分车
        line = [0]#每辆车服务客户点
        dis_sum = 0#线路距离
        dis,d = 0,0#当前客户距离前一个客户的距离、当前客户需求量
        i = 0#指向配送中心
        time_point = 0#
        wait = 0
        late = 0
        
        while i < len(bird):
            if line == [0]:#车辆未分配客户点
                dis += dis_matrix.loc[0,bird[i]]#记录距离
                line.append(bird[i])#为客户点分车
                d += Demand[bird[i]]#记录需求量
                time_point += dis_matrix.loc[0,bird[i]]/V
                if time_point < time[bird[i]][0]:
                    wait = time[bird[i]][0] - time_point
                    time_point = time_point + wait + time[bird[i]][2]
                elif time_point > time[bird[i]][1]:
                    late = time_point - time[bird[i]][1]
                    time_point = time_point + time[bird[i]][2]
                else:
                    time_point = time_point + time[bird[i]][2]
                
                i += 1#指向下一个客户点
            else:#已分配客户点则需判断车辆载重和行驶距离
                if (dis_matrix.loc[line[-1],bird[i]]+dis_matrix.loc[bird[i],0]+ dis <= DISTABCE) & (d + Demand[bird[i]]<=CAPACITY ) :
                    dis += dis_matrix.loc[line[-1],bird[i]]
                    time_point += dis_matrix.loc[line[-1],bird[i]]/V
                    if time_point < time[bird[i]][0]:
                        wait = time[bird[i]][0] - time_point
                        time_point = time_point + wait + time[bird[i]][2]
                    elif time_point > time[bird[i]][1]:
                        late = time_point - time[bird[i]][1]
                        time_point += time[bird[i]][2]
                    else:
                        time_point = time_point + time[bird[i]][2]
                    
                    line.append(bird[i])
                    d += Demand[bird[i]]
                    i += 1
                else:
                    dis += dis_matrix.loc[line[-1],0]#当前车辆装满
                    line.append(0)
                    dis_sum += dis
                    lines.append(line)
                    #下一辆车
                    dis,d = 0,0
                    line = [0]
                    time_point = 0
        
        #最后一辆车
        dis += dis_matrix.loc[line[-1],0]
        line.append(0)
        dis_sum += dis
        lines.append(line)
        birdPop_car.append(lines)
        fits.append(round(C1*dis_sum+C0*len(lines)+C2*wait+C3*late,1))
        
    return birdPop_car,fits


def crossover(bird,pLine,gLine,w,c1,c2):
    '''
    采用顺序交叉方式；交叉的parent1为粒子本身，分别以w/(w+c1+c2),c1/(w+c1+c2),c2/(w+c1+c2)
    的概率接受粒子本身逆序、当前最优解、全局最优解作为parent2,只选择其中一个作为parent2；
    输入：bird-粒子,pLine-当前最优解,gLine-全局最优解,w-惯性因子,c1-自我认知因子,c2-社会认知因子；
    输出：交叉后的粒子-croBird；
    '''
    croBird = [None]*len(bird)#初始化
    parent1 = bird#选择parent1
    #选择parent2（轮盘赌操作）
    randNum = random.uniform(0, sum([w,c1,c2]))
    if randNum <= w:
        parent2 = [bird[i] for i in range(len(bird)-1,-1,-1)]#bird的逆序
    elif randNum <= w+c1:
        parent2 = pLine
    else:
        parent2 = gLine
    
    #parent1-> croBird
    start_pos = random.randint(0,len(parent1)-1)
    end_pos = random.randint(0,len(parent1)-1)
    if start_pos>end_pos:start_pos,end_pos = end_pos,start_pos
    croBird[start_pos:end_pos+1] = parent1[start_pos:end_pos+1].copy()
    
    # parent2 -> croBird
    list2 = list(range(0,start_pos))
    list1 = list(range(end_pos+1,len(parent2)))
    list_index = list1+list2#croBird从后往前填充
    j = -1
    for i in list_index:
        for j in range(j+1,len(parent2)+1):
            if parent2[j] not in croBird:
                croBird[i] = parent2[j]
                break 
                    
    return croBird


def draw_path(car_routes,CityCoordinates):
    '''
    #画路径图
    输入：line-路径，CityCoordinates-城市坐标；
    输出：路径图
    '''
    for route in car_routes:
        x,y= [],[]
        for i in route:
            Coordinate = CityCoordinates[i]
            x.append(Coordinate[0])
            y.append(Coordinate[1])
        x.append(x[0])
        y.append(y[0])
        plt.plot(x, y,'o-', alpha=0.8, linewidth=0.8)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.show()


if __name__ == '__main__':
    #车辆参数
    CAPACITY = 8#车辆最大容量
    DISTABCE = 1000#车辆最大行驶距离
    V = 40#速度，km/h
    C0 = 100#启动成本
    C1 = 2#行驶成本/km
    C2 = 10#等待成本/h
    C3 = 40#惩罚成本/h
    
    #PSO参数
    birdNum = 50#粒子数量
    w = 0.2#惯性因子
    c1 = 0.4#自我认知因子
    c2 = 0.4#社会认知因子
    pBest,pLine =0,[]#当前最优值、当前最优解，（自我认知部分）
    gBest,gLine = 0,[]#全局最优值、全局最优解，（社会认知部分）
    
    #其他参数
    iterMax = 1000#迭代次数
    iterI = 1#当前迭代次数
    bestfit = [] #记录每代最优值
    
    #读入数据
    Customer = [(70,70),(107,77),(109,139),(120,22),(48,47),(116,22),(12,138),(86,40),(121,124),(61,57),
                (40,113),(129,24),(12,84),(44,116),(102,52),(41,36),(132,138),(104,139),(104,54),(22,104),(46,133)]
    Demand = [0,3.4,0.8,3.9,1.9,3.2,1.4,2.2,2.1,3.5,2.3,1.8,1.6,2.7,1.5,1.3,2.4,2.9,1.3,1.1,0.7]
    time = [(0,10,0),(0.5,4.5,0.2),(2,6.5,0.2),(1,6,0.2),(0.5,6,0.4),(1,6.5,0.2),(3,9,0.5),(0.5,4,0.4),(1.5,6,0.2),(3.5,9,0.2),
                   (1,4.5,0.2),(1.5,6.5,0.4),(1,5.5,0.4),(2,7,0.2),(2.5,6.5,0.2),(3,8,0.2),(2,7,0.4),(2,6,0.2),(0,4.5,0.4),(1,5.5,0.4),(1,5.5,0.4)]
    dis_matrix = calDistance(Customer)#计算城市间距离
    
    
    birdPop = [greedy(Customer,dis_matrix) for i in range(birdNum)]#贪婪算法构造初始解
    # birdPop = [random.sample(range(1,len(Customer)),len(Customer)-1) for i in range(birdNum)]#客户点编码，随机初始化生成种群
    
    birdPop_car,fits = calFitness(birdPop,Demand,dis_matrix,CAPACITY,DISTABCE,C0,C1,C2,C3,time,V)#分配车辆，计算种群适应度
    
    gBest = pBest = min(fits)#全局最优值、当前最优值
    gLine = pLine = birdPop[fits.index(min(fits))]#全局最优解、当前最优解
    gLine_car = pLine_car = birdPop_car[fits.index(min(fits))]
    bestfit.append(gBest)
    
    
    while iterI <= iterMax:#迭代开始
        for i in range(birdNum):
            birdPop[i] = crossover(birdPop[i],pLine,gLine,w,c1,c2)
        
        birdPop_car,fits = calFitness(birdPop,Demand,dis_matrix,CAPACITY,DISTABCE,C0,C1,C2,C3,time,V)#分配车辆，计算种群适应度
        pBest,pLine,pLine_car =  min(fits),birdPop[fits.index(min(fits))],birdPop_car[fits.index(min(fits))]
        if min(fits) <= gBest:
            gBest,gLine,gLine_car =  min(fits),birdPop[fits.index(min(fits))],birdPop_car[fits.index(min(fits))]
            
        bestfit.append(gBest)
        print(iterI,gBest)#打印当前代数和最佳适应度值
        iterI += 1#迭代计数加一
        
    print(gLine_car)#路径顺序
    draw_path(gLine_car,Customer)#画路径图
    
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235

4.3 例子求解结果

本文的测试例子借鉴期刊【混合蚁群算法求解双目标时间窗VRP】，具体数据如下:

车辆启动成本 C0 取100，车辆单位距离行驶成本C1 取2，等待成本取10，惩罚成本取40，车辆最长行驶距离为1000，车辆速度取40，算法运算结果为2616.5,，路径为[[0, 4, 15, 9, 0], [0, 1, 16, 8, 0], [0, 12, 19, 10, 13, 0], [0, 18, 14, 5, 0], [0, 7, 3, 11, 0], [0, 2, 17, 20, 6, 0]] ，路径图如下：

记录学习过程，欢迎指正