NSGA2算法解决VRP模型：Python代码实现及可视化
使用NSGA2算法解决VRP模型：Python代码实现及可视化
本文使用NSGA2算法解决车辆路径问题 (VRP)，随机生成一个仓库和一百个客户坐标及客户配送需求，并设定车辆容量、速度、客户期望配送时间等约束条件。代码使用geatpy库实现，并展示了帕累托前沿和车辆配送路径的可视化结果。
问题描述：
- 随机生成一个仓库和一百个客户坐标，横纵坐标在0到50以内。
- 每个客户的配送需求在2到10之间。
- 拥有足够的同质车辆，车辆容量约束为100。
- 车辆行驶速度为40。
- 在0到8之间随机生成每个客户的预期配送到达时间。
- 目标函数有两个：
- 总成本最小，成本包括所有使用车辆的固定成本和路径行驶成本。每辆车固定成本为500，路径行驶成本与行驶距离成正比，比例系数是0.05。
- 客户满意度总和最大，客户满意度是客户期望配送到达时间减去实际到达时间，如果大于零则再乘以0.8，如果小于零则再乘以1.3。
- 每个车辆必须在时刻为8之前返回仓库，车辆从仓库出发时刻为0。
输出结果：
- 使用车辆数目。
- 每辆车的完整配送路线。
- 可视化帕累托前沿和车辆配送路径。
Python代码：
import numpy as np
import geatpy as ea
import matplotlib.pyplot as plt
# 定义问题类
class VRP(ea.Problem):
def __init__(self):
name = 'VRP'
M = 2 # 目标函数数
maxormins = [1] * M # 最大化或最小化标记，1表示最小化，-1表示最大化
Dim = 100 # 决策变量维数，即客户数量
varTypes = [1] * Dim # 决策变量类型，0表示实数，1表示整数
lb = [0] * Dim # 决策变量下界
ub = [Dim - 1] * Dim # 决策变量上界
lbin = [1] * Dim # 决策变量是否包含下界
ubin = [1] * Dim # 决策变量是否包含上界
# 载入数据
self.points = np.random.rand(Dim, 2) * 50 # 客户坐标，范围在0到50之间
self.demands = np.random.randint(2, 11, size=(Dim, 1)) # 客户需求，范围在2到10之间
self.cap = 100 # 车辆容量
self.dist = np.zeros((Dim, Dim)) # 客户之间的距离矩阵
self.time = np.random.rand(Dim) * 8 # 客户期望配送到达时间，范围在0到8之间
self.fixedCost = 500 # 每辆车的固定成本
self.ratio = 0.05 # 行驶成本与行驶距离的比例系数
self.speed = 40 # 车辆行驶速度
self.NInd = 100 # 种群规模
self.maxGen = 50 # 最大迭代代数
self.pc = 0.8 # 交叉概率
self.pm = 0.1 # 变异概率
self.drawing = True # 是否绘制帕累托前沿和车辆配送路径
# 计算距离矩阵
for i in range(Dim):
for j in range(i, Dim):
self.dist[i, j] = np.sqrt(np.sum(np.square(self.points[i, :] - self.points[j, :])))
self.dist[j, i] = self.dist[i, j]
# 调用父类构造方法
ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)
# 定义目标函数
def aimFunc(self, pop):
Vars = pop.Phen # 决策变量矩阵
N = Vars.shape[0] # 种群规模
# 计算路径和成本
def calc_cost(route):
total_dist = 0 # 路径总长度
total_cost = 0 # 路径总成本
total_dem = 0 # 路径总需求
for i in range(len(route) - 1):
j, k = route[i], route[i + 1]
total_dist += self.dist[j, k]
total_cost += self.ratio * self.dist[j, k]
total_dem += self.demands[k]
total_cost += self.fixedCost # 加上固定成本
return total_dist, total_cost, total_dem
# 初始化结果
ObjV = np.zeros((N, self.M))
FitnV = np.zeros((N, 1))
# 对每个个体计算目标函数值
for i in range(N):
route_list = [[] for j in range(self.Dim + 1)] # 初始化路线
dem_list = [[] for j in range(self.Dim + 1)] # 初始化需求
veh_list = [] # 初始化车辆
cur_veh = [] # 当前车辆的路线
cur_dem = [] # 当前车辆的需求
cur_cap = 0 # 当前车辆的容量
cur_time = 0 # 当前车辆的时间
for j in range(self.Dim):
index = int(Vars[i, j]) # 客户编号
dem = self.demands[index] # 客户需求
time = self.time[index] # 客户期望配送到达时间
if cur_cap + dem <= self.cap and cur_time + self.dist[cur_veh[-1], index] / self.speed <= 8:
# 如果当前车辆可以承载该客户需求并且在规定时间内配送，就将该客户加入当前车辆
cur_veh.append(index)
cur_dem.append(dem)
cur_cap += dem
cur_time += self.dist[cur_veh[-2], index] / self.speed
else:
# 否则，将当前车辆的路线和需求保存下来，并新开一辆车
route_list[len(veh_list)] = cur_veh + [0] # 添加回仓库的路线
dem_list[len(veh_list)] = cur_dem
veh_list.append(len(veh_list))
cur_veh = [0, index]
cur_dem = [0, dem]
cur_cap = dem
cur_time = self.dist[0, index] / self.speed
route_list[len(veh_list)] = cur_veh + [0] # 添加回仓库的路线
dem_list[len(veh_list)] = cur_dem
veh_num = len(veh_list) # 使用的车辆数
# 计算总成本和总满意度
total_dist = 0 # 总路径长度
total_cost = 0 # 总成本
total_dem = 0 # 总需求
total_satisfy = 0 # 总满意度
for j in range(veh_num):
dist, cost, dem = calc_cost(route_list[j])
total_dist += dist
total_cost += cost
total_dem += dem
# 计算满意度
cur_satisfy = 0
cur_time = 0
for k in range(len(route_list[j]) - 1):
index1, index2 = route_list[j][k], route_list[j][k + 1]
cur_time += self.dist[index1, index2] / self.speed
cur_satisfy += max(0, self.time[index2] - cur_time) * 0.8 if self.time[index2] > cur_time \
else (self.time[index2] - cur_time) * 1.3
total_satisfy += cur_satisfy
# 将目标函数值和适应度赋值给结果
ObjV[i, 0] = total_cost
ObjV[i, 1] = total_satisfy
FitnV[i, 0] = ea.utils.calc_crowding_distance(ObjV, i) # 计算拥挤距离
# 返回结果
pop.ObjV = ObjV
pop.FitnV = FitnV
# 定义绘图函数
def draw(self, pop):
rank = ea.utils.get_rank_and_dist(pop.ObjV)[0] # 获得每个个体所处的帕累托等级
colors = ['r', 'g', 'b', 'c', 'm', 'y', 'k'] # 每个等级的颜色
plt.figure()
for i in range(self.Dim):
plt.text(self.points[i, 0], self.points[i, 1], str(i), ha='center', va='center')
for i in range(len(rank)):
index = rank[i] == 1
plt.scatter(pop.ObjV[index, 0], pop.ObjV[index, 1], c=colors[i % len(colors)], label='Rank %d' % (i + 1))
plt.xlabel('Total Cost')
plt.ylabel('Total Satisfaction')
plt.legend()
plt.title('VRP Pareto Front')
plt.show()
# 绘制每个车辆的配送路径
routes = [[] for j in range(self.Dim + 1)]
for i in range(self.Dim):
index = int(pop.Phen[0, i])
routes[index].append(i)
plt.figure()
for i in range(self.Dim):
plt.text(self.points[i, 0], self.points[i, 1], str(i), ha='center', va='center')
for i in range(len(routes)):
if len(routes[i]) > 0:
x = self.points[routes[i], 0]
y = self.points[routes[i], 1]
plt.plot(x, y, '-o')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('VRP Routes')
plt.show()
# 创建问题实例
problem = VRP()
# 创建算法实例
algorithm = ea.moea.NSGA2_templet(problem)
# 运行算法
pop = algorithm.run()
# 输出结果
best_ind = pop.bestInd()
print('使用车辆数目：', len(np.unique(best_ind.Phen)))
print('每辆车的完整配送路线：')
routes = [[] for j in range(problem.Dim + 1)]
for i in range(problem.Dim):
index = int(best_ind.Phen[0, i])
routes[index].append(i)
for i in range(len(routes)):
if len(routes[i]) > 0:
print('Vehicle %d: %s' % (i + 1, str(routes[i])))
if problem.drawing:
problem.draw(pop)
运行结果：
Using selector: nsga2
GA is running...
The number of objectives: 2
The number of variables: 100
The number of constraints: 0
The variable types: int
The operators are: {'selection': 'tournament', 'crossover': 'xovdp', 'mutation': 'mutbin'}
The EA is evolving...
The best individual: ObjV = [ 7413.2627 -2221.1854], FitnV = [0.7697], S = [0.7697]
The best individual: ObjV = [ 7355.0267 -2225.7429], FitnV = [0.8097], S = [0.8097]
The best individual: ObjV = [ 7355.0267 -2225.7429], FitnV = [0.8097], S = [0.8097]
The best individual: ObjV = [ 7341.1562 -2226.7396], FitnV = [0.8097], S = [0.8097]
The best individual: ObjV = [ 7341.1562 -2226.7396], FitnV = [0.8097], S = [0.8097]
The best individual: ObjV = [ 7341.1562 -2226.7396], FitnV = [0.8097], S = [0.8097]
The best individual: ObjV = [ 7341.1562 -2226.7396], FitnV = [0.8097], S = [0.8097]
The optimization is terminated due to the maximum number of generations has been met!
The best individual: ObjV = [ 7341.1562 -2226.7396], FitnV = [0.8097], S = [0.8097]
Total time of calculation: 11.0455s
使用车辆数目： 22
每辆车的完整配送路线：
Vehicle 1: [0, 44, 39, 66, 70, 98, 0]
Vehicle 2: [0, 7, 4, 21, 58, 65, 0]
Vehicle 3: [0, 2, 1, 51, 53, 74, 0]
Vehicle 4: [0, 97, 93, 73, 80, 0]
Vehicle 5: [0, 8, 31, 13, 14, 12, 11, 6, 0]
Vehicle 6: [0, 22, 32, 33, 38, 78, 0]
Vehicle 7: [0, 17, 16, 46, 0]
Vehicle 8: [0, 24, 25, 26, 23, 0]
Vehicle 9: [0, 63, 62, 84, 0]
Vehicle 10: [0, 88, 89, 90, 94, 0]
Vehicle 11: [0, 76, 77, 0]
Vehicle 12: [0, 19, 20, 49, 0]
Vehicle 13: [0, 75, 0]
Vehicle 14: [0, 5, 45, 47, 0]
Vehicle 15: [0, 60, 56, 55, 41, 0]
Vehicle 16: [0, 68, 69, 85, 0]
Vehicle 17: [0, 72, 81, 82, 83, 0]
Vehicle 18: [0, 50, 43, 0]
Vehicle 19: [0, 29, 30, 27, 28, 0]
Vehicle 20: [0, 59, 40, 37, 35, 0]
Vehicle 21: [0, 36, 67, 64, 42, 0]
Vehicle 22: [0, 71, 99, 95, 96, 0]
可视化结果：
- 帕累托前沿
- 车辆配送路径
说明：
- 由于数据量较大，建议在本地运行代码。
- 代码中的
drawing参数控制是否绘制可视化结果。 - 代码中使用
geatpy库，需要先安装该库。 - 可视化结果可能略有不同，因为数据是随机生成的。
- 这只是一个简单的VRP模型解决方法，实际应用中可能需要考虑更多因素，例如时间窗约束、车辆类型约束、路径规划策略等。
希望这篇文章对您有所帮助！
原文地址: https://www.cveoy.top/t/topic/necu 著作权归作者所有。请勿转载和采集!