基于NSGA2算法的VRP模型求解及可视化

由于该问题涉及到多个约束和目标函数，需要使用多目标优化算法，本文选用NSGA2算法来解决该问题。以下是完整代码：

import numpy as np
import geatpy as ea
import matplotlib.pyplot as plt
import seaborn as sns

# 定义目标函数类
class VRPProblem(ea.Problem):             
    def __init__(self):
        name = 'VRPProblem'            # 初始化函数名称
        M = 2                          # 两个目标函数
        maxormins = [1, 1]             # 两个目标函数均为最小化问题
        Dim = 100 * 2 + 10 + 10        # 10辆车，每辆车最多服务10个客户，每个客户有横纵坐标和需求量
        varTypes = [1] * Dim           # 所有变量均为实数变量
        lb = [0] * Dim                 # 所有变量的下界为0
        ub = [50] * 200 + [10] * 10    # 前200个变量为横纵坐标，后面20个变量为需求量和预期送达时间
        lbin = [1] * Dim               # 所有变量的取值范围均为闭区间
        ubin = [1] * Dim
        self.customer_num = 100        # 客户数量
        self.vehicle_num = 10          # 车辆数量
        self.capacity = 100            # 车辆容量
        self.speed = 40                # 车辆速度
        self.start_time = 8            # 回到仓库的最迟时间
        self.fixed_cost = 500          # 每辆车的固定成本
        self.distance_cost = 0.05      # 路径行驶成本与行驶距离成正比的比例系数
        self.customers = np.zeros((self.customer_num, 3))  # 客户信息（坐标、需求量、预期送达时间）
        self.customers[:, :2] = np.random.rand(self.customer_num, 2) * 50
        self.customers[:, 2] = np.random.randint(2, 11, self.customer_num)
        self.customers[:, 3] = np.random.rand(self.customer_num) * 8
        self.customers = self.customers.astype(int)
        ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)  # 调用父类构造方法

    # 计算每个客户的满意度
    def calc_satisfaction(self, routes):
        satisfaction = 0
        for i, route in enumerate(routes):
            time = 0
            for j in range(len(route)):
                customerId = route[j]
                customer = self.customers[customerId]
                distance = np.sqrt((customer[0] - time) ** 2 + (customer[1] - time) ** 2)  # 客户到达时间
                time += distance / self.speed  # 更新时间
                satisfaction += (customer[3] - time) * (customer[3] > time) * 0.8 + (time - customer[3]) * (time >= customer[3]) * 1.3
            satisfaction += (self.start_time - time) * (self.start_time > time) * 0.8 + (time - self.start_time) * (time >= self.start_time) * 1.3
        return satisfaction

    # 计算成本
    def calc_cost(self, routes):
        cost = 0
        for i, route in enumerate(routes):
            distance = 0
            time = 0
            for j in range(len(route)):
                customerId = route[j]
                customer = self.customers[customerId]
                distance += np.sqrt((customer[0] - time) ** 2 + (customer[1] - time) ** 2)  # 计算客户到达时间
                time += distance / self.speed  # 更新时间
            distance += np.sqrt((0 - time) ** 2 + (0 - time) ** 2)  # 计算回到仓库的时间
            cost += self.fixed_cost + distance * self.distance_cost  # 计算成本
        return cost

    # 目标函数
    def aimFunc(self, pop):
        x = pop.Phen  # 取出种群所有个体的决策变量
        routes = [[] for i in range(self.vehicle_num)]  # 路径列表
        loads = [0] * self.vehicle_num  # 车辆负载列表
        for i in range(self.customer_num):
            vehicle = np.argmin(loads)  # 找到负载最小的车辆
            if loads[vehicle] + self.customers[i][2] <= self.capacity:  # 判断是否能够装下该客户的货物
                routes[vehicle].append(i)
                loads[vehicle] += self.customers[i][2]
        satisfaction = self.calc_satisfaction(routes)  # 计算客户满意度
        cost = self.calc_cost(routes)  # 计算成本
        pop.ObjV[:, 0] = cost  # 第一个目标函数为总成本最小
        pop.ObjV[:, 1] = -satisfaction  # 第二个目标函数为总满意度最大（注意要取负数）

接下来，我们需要使用NSGA2算法来求解该问题，代码如下：

# 实例化NSGA2算法对象
algorithm = ea.moea.NSGA2_templet(problem=VRPProblem(), population_size=100, offspring_population_size=100, max_iteration=100)

# 运行算法
algorithm.run()

# 输出结果
best_x = algorithm.best_gen.Phen[0]
routes = [[] for i in range(10)]
loads = [0] * 10
for i in range(100):
    vehicle = np.argmin(loads)
    if loads[vehicle] + VRPProblem().customers[i][2] <= 100:
        routes[vehicle].append(i)
        loads[vehicle] += VRPProblem().customers[i][2]

print("使用车辆数：", len([load for load in loads if load > 0]))
print("每辆车的完整配送路线：")
for i in range(len(routes)):
    print("车辆", i + 1, "：仓库 ->", end=" ")
    for customer in routes[i]:
        print(customer + 1, "->", end=" ")
    print("仓库")

# 可视化帕累托前沿
data = algorithm.optima_history
sns.scatterplot(x=data[:, 0], y=-data[:, 1])
plt.xlabel("Total cost")
plt.ylabel("Total satisfaction")
plt.show()

# 可视化车辆配送路径
plt.figure(figsize=(8, 8))
plt.scatter([0], [0], c="blue", s=200, marker="s")
for i in range(10):
    route = routes[i]
    x = [VRPProblem().customers[j][0] for j in route]
    y = [VRPProblem().customers[j][1] for j in route]
    plt.scatter(x, y, c="orange")
    plt.plot(x + [0], y + [0], c="orange")
plt.xlim([0, 50])
plt.ylim([0, 50])
plt.xlabel("X")
plt.ylabel("Y")
plt.show()

运行结果如下：

using pygmo

# 先定义一个可行性函数
def feasible(x):
    loads = [0] * 10
    for i in range(100):
        vehicle = np.argmin(loads)
        if loads[vehicle] + x[i * 2 + 2] <= 100:
            loads[vehicle] += x[i * 2 + 2]
        else:
            return False
    return True

# 然后定义目标函数
def objective(x):
    routes = [[] for i in range(10)]
    loads = [0] * 10
    for i in range(100):
        vehicle = np.argmin(loads)
        if loads[vehicle] + x[i * 2 + 2] <= 100:
            routes[vehicle].append(i)
            loads[vehicle] += x[i * 2 + 2]
    satisfaction = 0
    for i, route in enumerate(routes):
        time = 0
        for j in range(len(route)):
            customerId = route[j]
            customer = customers[customerId]
            distance = np.sqrt((customer[0] - time) ** 2 + (customer[1] - time) ** 2)
            time += distance / 40
            satisfaction += (customer[3] - time) * (customer[3] > time) * 0.8 + (time - customer[3]) * (time >= customer[3]) * 1.3
        satisfaction += (8 - time) * (8 > time) * 0.8 + (time - 8) * (time >= 8) * 1.3
    cost = sum([500 + np.sqrt((customers[route[0]][0] - 0) ** 2 + (customers[route[0]][1] - 0) ** 2) * 0.05 +             sum([np.sqrt((customers[route[j]][0] - customers[route[j - 1]][0]) ** 2 + (customers[route[j]][1] - customers[route[j - 1]][1]) ** 2) * 0.05 for j in range(1, len(route))]) +             np.sqrt((0 - customers[route[-1]][0]) ** 2 + (0 - customers[route[-1]][1]) ** 2) * 0.05 for route in routes])
    return [cost, -satisfaction]

# 然后使用pygmo来求解
algo = pg.algorithm(pg.nsga2(gen=1000))
prob = pg.problem(pg.real(220))
prob.set_bounds(0, 50)
prob.set_bounds(200, 220, 2, 10)
prob.set_bounds(220, 240, 0, 8)
prob.set_feasibility(feasible)
prob.set_evaluate(objective)
pop = pg.population(prob, 100)
pop = algo.evolve(pop)
data = np.array(pop.get_f())
sns.scatterplot(x=data[:, 0], y=-data[:, 1])
plt.xlabel("Total cost")
plt.ylabel("Total satisfaction")
plt.show()

best_x = pop.get_x()[0]
routes = [[] for i in range(10)]
loads = [0] * 10
for i in range(100):
    vehicle = np.argmin(loads)
    if loads[vehicle] + best_x[i * 2 + 2] <= 100:
        routes[vehicle].append(i)
        loads[vehicle] += best_x[i * 2 + 2]

print("使用车辆数：", len([load for load in loads if load > 0]))
print("每辆车的完整配送路线：")
for i in range(len(routes)):
    print("车辆", i + 1, "：仓库 ->", end=" ")
    for customer in routes[i]:
        print(customer + 1, "->", end=" ")
    print("仓库")

plt.figure(figsize=(8, 8))
plt.scatter([0], [0], c="blue", s=200, marker="s")
for i in range(10):
    route = routes[i]
    x = [customers[j][0] for j in route]
    y = [customers[j][1] for j in route]
    plt.scatter(x, y, c="orange")
    plt.plot(x + [0], y + [0], c="orange")
plt.xlim([0, 50])
plt.ylim([0, 50])
plt.xlabel("X")
plt.ylabel("Y")
plt.show()