智能优化系列|第五期：禁忌搜索算法求解VRPTW问题（带时间窗的车辆路径规划问题）- 附python代码

车辆路径问题（VRP）最早是由 Dantzig 和 Ramser 于1959年首次提出，它是指一定数量的客户，各自有不同数量的货物需求，配送中心向客户提供货物，由一个车队负责分送货物，组织适当的行车路线，目标是使得客户的需求得到满足，并能在一定的约束下，达到诸如路程最短、成本最小、耗费时间最少等目的。

Miracle_lucy

2570人浏览 · 2024-09-02 09:00:00

Miracle_lucy · 2024-09-02 09:00:00 发布

经过上一篇无形忍者-也许你和禁忌搜索算法的距离就差这一文(附python案例代码) 的概念介绍和算法核心步骤介绍，我们对禁忌搜索已经有了一定的了解。下面我们来一起复习下禁忌搜索算法的核心步骤：

1.禁忌搜索算法核心步骤

① 初始化

利用贪婪算法等局部搜索算法生成一个初始解，清空禁忌表，设置禁忌长度。

② 邻域搜索产生候选解

根据步骤①产生的初始解，通过搜索算子（search operators），如relocation、exchange、2-opt等，生成候选解(candidate solution)，并计算各个候选解的适应值（即解对应的目标函数值）。

③ 选择最好的候选解

从步骤②产生的所有候选解中选出这一轮适应值最好的候选解，将其与全局最优解进行比较，

如果优于全局最优解，那么就不考虑其是否被禁忌，用这个最好的候选解来更新全局最优解，并且作为下一个迭代的当前解，然后将对应的操作加入禁忌表；

如果不优于全局最优解，就从所有候选解中选出不在禁忌状态下的最好解作为新的当前解，并且作为下一个迭代的当前解，然后将对应的操作加入禁忌表；这时候不更新全局最优解。

④ 藐视准则

如步骤3中介绍的：如果优于全局最优解，那么就不考虑其是否被禁忌，用这个最好的候选解来更新全局最优解，并且作为下一个迭代的当前解，然后将对应的操作加入禁忌表；

⑤ 判断终止条件

若满足终止条件，则立即停止并输出当前最好解；否则继续搜索。一般终止条件为是否到达一定的迭代次数或者达到了一个时间限制。

2.带时间窗车辆路径问题（VRPTW）简介

车辆路径问题（VRP）最早是由 Dantzig 和 Ramser 于1959年首次提出，它是指一定数量的客户，各自有不同数量的货物需求，配送中心向客户提供货物，由一个车队负责分送货物，组织适当的行车路线，目标是使得客户的需求得到满足，并能在一定的约束下，达到诸如路程最短、成本最小、耗费时间最少等目的。

由于VRP问题的持续发展，考虑需求点对于车辆到达的时间有所要求之下，在车辆途程问题之中加入时窗的限制，便成为带时间窗车辆路径问题（VRP with Time Windows, VRPTW）。带时间窗车辆路径问题（VRPTW）是在VRP上加上了客户的被访问的时间窗约束。

3.VRPTW问题的建模实例

3.1问题背景

VRPTW 问题是不仅考虑CVRP的所有约束，还需要考虑时间窗约束，也就是每个顾客对应一个时间窗 $[a_i,b_i]$ ,其中和 $a_i$ 和 $b_i$ 分别代表该点的最早到达时间和最晚到达时间。

3.2约束条件

①每辆车必须从仓库出发，最终返回仓库；

②每个客户只被1辆车服务；

③每辆车的载重不超过总容量Q；

④每辆车的总时间不超过Dk；

⑤需要在时间窗 $[a_i,b_i]$ 内服务客户i,并且需要在 $[a_0,b_0]$ 时间窗内返回仓库点

3.3目标

每辆车辆的总时间最短（由于车辆速度相同，时间最短相当于路程最短）的路线。（允许不使用某些车辆）

4.利用禁忌搜索解决VRPTW问题

算法的核心步骤

首先，我们怎么去产生一个初始解？这是利用的贪心的思想为每个客户指派路径：

尝试将客户i加入到某条路径k中，
判断是否满足路径k的容量约束条件，如果不满足则新建一条路径，k=k+1;
然后再路径k中，任意找节点j,j+1;如果客户i满足时间窗满足 $a_j \leq a_i \leq a_{j+1}$ 。
我们就将客户i分派到这条路径k,并放置在节点j,j+1之间，同时累加该路径的总长度。

然后开始禁忌搜索，在禁忌搜索的每一轮迭代的过程中，我们都需要找出一个最好的候选解，这个候选解与全局最优解进行比较，如果优于全局最优解，那么就不考虑其是否被禁忌，用这个最好的候选解来更新全局最优解，并且作为下一个迭代的当前解，然后将对应的操作加入禁忌表；如果不优于全局最优解，就从所有候选解中选出不在禁忌状态下的最好解作为新的当前解，并且作为下一个迭代的当前解，然后将对应的操作加入禁忌表；这时候不更新全局最优解。

那么我们怎么搜索最好的候选解呢？

首先我们假设每一轮迭代中的当前解为bestv，当前最好客户为bestc，当前最好路径为bestr，当前最好位置为bestp。然后，我们对每一个客户节点i进行遍历，先找出它在当前路径中的位置，即为P，接着我们将客户i从原路径的第P个位置中移除，然后遍历所有路径的所有节点，取任意路径k,将客户i插入路径k的任意位置x，计算路径k加入新节点i后的评价函数值，记为tempv。我们在每轮迭代中找到最好的候选解，只符合①-③的限制条件，所以我们在每一轮找到bestv后，需要check是否满足条件④-⑤，如果也满足④-⑤ ，才更新全局最佳解为Ans。

4.1插入算子

具体的插入操作为:

1. 将客户i从原路径k0的第P个位置中移除

2. 选择路径k(k≠ k0), 确保加入节点 i 后路径k不会超过容量约束的上限。

3. 在路径 k中找到点x,使x 满足 $a_x \leq a_i \leq a_{x+1}$ 。

4. 往路径 k中插入点i,计算路径k加入新节点i后的评价函数值。

5.恢复路径k0，恢复路径k。

4.2评价函数

我们引入目标函数主要由三个部分组成：d路径总长度（优化目标），q超出容量约束总量，t超出时间窗约束总量
目标函数结构为 f(r) = d + alpha * q + beta* t,

第一项为问题最小化目标，后两项为惩罚部分，其中alpha与beta为变量，分别根据当前解是否满足两个约束进行变化，根据每轮迭代得到的解在check函数中更新：分别根据约束满足的情况和控制系数sita更新惩罚系数alpha和beta值。新路径满足条件，惩罚系数减小；新路径违反条件，惩罚系数加大。

4.3禁忌表

我们这里用了2个禁忌表：禁忌表tabu是一个二维表,禁忌将节点i插入路线k的操作; 禁忌表tabu_create是一个一维表,禁忌使用新车辆(开始新路线) k'的操作。

5.python代码实现

代码主要分为以下几个类：

5.1代码模块说明

代码一共分为2个模块：实体模块和算法模块。

实体模块中2个类：CustomerType和RouteType；

算法模块TabuSearch中5个主要的方法：read_data()读取初始文件信息；init_solve()构造初始解；check()检验是否满足容量约束和时间窗约束;evaluate()计算路径规划r的目标函数值；tabu_search()完成核心的禁忌搜索流程，算法中采取插入算子，即从一条路径中选择一点首先从原路径所在的位置中移除，遍历插入到任意车辆路径的任意位置，在该操作下形成的邻域中选取使目标函数最小化的解，然后恢复源路径，继续遍历重新选择。

5.2文本数据输入格式说明

采取读取文本格式来作为数据输入的形式，具体格式见下表：

序号	横坐标	纵坐标	服务量	时间窗起始	时间窗结束	服务时间
0	35.00	35.00	0.00	0.00	230．00	0.00
1	41.00	49.00	10.00	161.00	171.00	10.00

由上表可知，输入数据为一行七列的形式，依次为对应点的序号、横坐标、纵坐标、所需服务量、服务时间窗起始时间点、服务时间窗结束时间点和服务所需时间。由此可见，节点0为仓库（depot）,其他节点为待服务点。

5.3部分代码展示

import copy
import time
from math import sqrt
import random
from entity_vrptw import CustomerType, RouteType


class TabuSearch:
    def __init__(self):
        self.inf = float('inf')
        self.customer_number = 100  # 算例中除仓库以外的顾客节点个数
        self.vehicle_number = 25
        self.depot_number = 1  # 算例中仓库节点个数
        self.capacity = 200  # 车辆最大容量
        self.iter_max = 300  # 搜索迭代次数

        self.tabu_tenure = 20  # 禁忌步长
        self.tabu = [[0 for _ in range(self.vehicle_number)] for _ in
                     range(self.customer_number + self.depot_number)]  # 禁忌表用于禁忌节点插入操作:[i][j]将节点i插入路径j中
        self.tabu_create = [0 for _ in range(self.customer_number + self.depot_number)]  # 禁忌表用于紧急拓展新路径或使用新车辆

        self.ans = 0.0  # 最优解距离
        self.alpha, self.beta, self.sita = 1, 1, 0.5  # alpha，beta为系数，计算目标函数值；sita控制系数改变的速度
        self.dis_matrix = [[0.0 for _ in range(self.customer_number + self.depot_number)] for _ in
                           range(self.customer_number + self.depot_number)]  # 距离矩阵
        self.customers = [None] * (self.customer_number + self.depot_number)  # 存储客户数据
        self.routes = [None] * self.vehicle_number  # 存储当前解路线数据
        self.route_ans = [None] * self.vehicle_number  # 存储最优解路线数据

    # 检查是否满足约束条件
    def check(self, r):
        q = 0
        t = 0

        # 检查是否满足容量约束
        for i in range(0, self.vehicle_number):
            if len(r[i].v) > 2 and r[i].load > self.capacity:  # 对有客户且超过容量约束的路径，记录超过部分
                q += r[i].load - self.capacity

        # 检查是否满足时间窗约束
        for i in range(0, self.vehicle_number):
            t += r[i].subt

        # 分别根据约束满足的情况和控制系数sita更新alpha和beta值
        # 新路径满足条件，惩罚系数减小，
        # 新路径违反条件，惩罚系数加大。
        if q == 0 and self.alpha >= 0.001:
            self.alpha /= (1 + self.sita)
        elif q != 0 and self.alpha <= 2000:
            self.alpha *= (1 + self.sita)

        if t == 0 and self.beta >= 0.001:
            self.beta /= (1 + self.sita)
        elif t != 0 and self.beta <= 2000:
            self.beta *= (1 + self.sita)

        if t == 0 and q == 0:
            return True
        else:
            return False

    # 更新路径r对时间窗的违反量
    def update_subt(self, r):
        arrive_time = 0
        for j in range(1, len(r.v)):
            arrive_time = arrive_time + r.v[j - 1].service + self.dis_matrix[r.v[j - 1].number][r.v[j].number]
            if arrive_time > r.v[j].end:
                r.subt = r.subt + arrive_time - r.v[j].end
            elif arrive_time < r.v[j].begin:
                arrive_time = r.v[j].begin

    # 计算路径规划r的目标函数值，通过该目标函数判断解是否较优
    def evaluate(self, r):
        # 目标函数主要由三个部分组成：d路径总长度（优化目标），q超出容量约束总量，t超出时间窗约束总量
        # 目标函数结构为 f(r) = d + Alpha * q + Beta * t, 第一项为问题最小化目标，后两项为惩罚部分
        # 其中alpha与beta为变量，分别根据当前解是否满足两个约束进行变化，根据每轮迭代得到的解在check函数中更新
        q = 0
        t = 0
        d = 0

        # 计算单条路径超出容量约束的总量
        for i in range(0, self.vehicle_number):
            if len(r[i].v) > 2 and r[i].load > self.capacity:
                q += r[i].load - self.capacity

        # 计算总超出时间
        for i in range(0, self.vehicle_number):
            t += r[i].subt

        # 计算路径总长度
        for i in range(0, self.vehicle_number):
            d += r[i].dis

        return d + self.alpha * q + self.beta * t  # 返回目标函数值
    
    @staticmethod
    def dis_cal(c1, c2):  # 计算图上各节点间的距离
        return sqrt((c1.x - c2.x) ** 2 + (c1.y - c2.y) ** 2)

    def read_data(self):

        with open("r101.txt", 'r') as file:  # 更换了算例
            for i in range(0, self.customer_number + self.depot_number):  # 把客户数据存在字典customers中，字典的长度是101
                line = file.readline().split()
                c = CustomerType()
                c.number = int(line[0])
                c.x = float(line[1])
                c.y = float(line[2])
                c.demand = float(line[3])
                c.begin = float(line[4])
                c.end = float(line[5])
                c.service = float(line[6])
                self.customers[i] = c
        # 初始化每条路径，默认路径收尾为仓库，
        # 首仓库最早最晚时间均为原仓库最早时间，
        # 尾仓库最早最晚时间均为原仓库最晚时间
        for i in range(0, self.vehicle_number):  # 把路由数据存在字典routes中
            r = RouteType()
            r.v.append(copy.deepcopy(self.customers[0]))
            r.v.append(copy.deepcopy(self.customers[0]))
            r.v[0].end = r.v[0].begin
            r.v[1].begin = r.v[1].end
            r.load = 0
            self.routes[i] = r
        self.ans = self.inf
        for i in range(0, self.customer_number + 1):
            for j in range(0, self.customer_number + 1):
                self.dis_matrix[i][j] = self.dis_cal(self.customers[i], self.customers[j])

    # 构造初始解
    def init_solve(self):
        customer_set = [i for i in range(0, self.customer_number + self.depot_number)] # 存的客户节点的编码:0,1,2,3,4,5,6...100
        sizeof_customer_set = self.customer_number
        current_route = 0

        while sizeof_customer_set > 0:
            k = random.randint(1, sizeof_customer_set) #随机挑选一个节点的index值，令其为K，K的范围是[1-101)左闭右开
            c = customer_set[k]
            customer_set[k] = customer_set[sizeof_customer_set]
            sizeof_customer_set -= 1

            if self.routes[current_route].load + self.customers[c].demand > self.capacity:
                current_route += 1
            for i in range(len(self.routes[current_route].v) - 1):
                if (self.routes[current_route].v[i].begin <= self.customers[c].begin) and (
                        self.customers[c].begin <= self.routes[current_route].v[i + 1].begin):
                    self.routes[current_route].v.insert(i + 1, copy.deepcopy(self.customers[c]))
                    self.routes[current_route].load += self.customers[c].demand
                    self.customers[c].r = current_route
                    break

        for i in range(0, self.vehicle_number):
            self.routes[i].subt = 0
            self.routes[i].dis = 0

            for j in range(1, len(self.routes[i].v)):
                self.routes[i].dis += self.dis_matrix[self.routes[i].v[j - 1].number][
                    self.routes[i].v[j].number]

            self.update_subt(self.routes[i])

    def print_ans(self, routes):
        print("************************************************************")
        print("The Minimum Total Distance = ", self.ans)
        print("Concrete Schedule of Each Route as Following : ")

        m = 0
        for i in range(0, self.vehicle_number):
            if len(routes[i].v) > 2:
                m += 1
                print("No.", m, " : ", end="")

                for j in range(len(routes[i].v)):
                    if j < len(routes[i].v) - 1:
                        print(routes[i].v[j].number, "->", end="")
                    else:
                        print(routes[i].v[j].number)

        print("************************************************************")

    def check_ans(self, routes):
        # 检验距离计算是否正确
        check_ans = 0
        for i in range(0, self.vehicle_number):
            for j in range(1, len(routes[i].v)):
                check_ans += self.dis_matrix[routes[i].v[j - 1].number][
                     routes[i].v[j].number]

        print("Check_Ans =", check_ans)

        # 检验是否满足时间窗约束
        flag = True
        for i in range(0, self.vehicle_number):
            self.update_subt(routes[i])
            if routes[i].subt > 0:
                flag = False

        if flag:
            print("Solution satisfies time windows construction")
        else:
            print("Solution does not satisfy time windows construction")

    def add_node(self, r, pos, cus):
        # 更新在路径r中加上节点cus的需求
        self.routes[r].load += self.customers[cus].demand

        # 更新在路径r中插入节点cus后所组成的路径距离
        self.routes[r].dis = self.routes[r].dis - \
                             self.dis_matrix[self.routes[r].v[pos - 1].number][self.routes[r].v[pos].number] + \
                             self.dis_matrix[self.routes[r].v[pos - 1].number][self.customers[cus].number] + \
                             self.dis_matrix[self.routes[r].v[pos].number][self.customers[cus].number]

        # 在路径r中插入节点cus
        self.routes[r].v.insert(pos, copy.deepcopy(self.customers[cus]))

    def remove_node(self, r, pos, cus):
        # 更新在路径r中去除节点cus的需求
        self.routes[r].load -= self.customers[cus].demand

        # 更新在路径r中去除节点cus后所组成的路径的距离
        self.routes[r].dis = self.routes[r].dis - \
                             self.dis_matrix[self.routes[r].v[pos - 1].number][self.routes[r].v[pos].number] - \
                             self.dis_matrix[self.routes[r].v[pos].number][self.routes[r].v[pos + 1].number] + \
                             self.dis_matrix[self.routes[r].v[pos - 1].number][self.routes[r].v[pos + 1].number]

        # 从路径r中去除节点cus
        del self.routes[r].v[pos]


def main():
    begin_time = time.time()
    ts = TabuSearch()
    ts.read_data()
    ts.init_solve()
    ts.print_ans(ts.routes)
    ts.check_ans(ts.routes)

    ts.tabu_search()
    ts.print_ans(ts.route_ans)
    ts.check_ans(ts.route_ans)

    end_time = time.time()
    used_time = end_time - begin_time
    print()
    print("程序耗时：" + str(used_time) + "s")


if __name__ == "__main__":
    main()

5.4结果输出

************************************************************
The Minimum Total Distance =  inf
Concrete Schedule of Each Route as Following : 
No. 1  : 0 ->65 ->44 ->30 ->71 ->7 ->86 ->53 ->6 ->57 ->46 ->37 ->54 ->0
No. 2  : 0 ->72 ->76 ->85 ->8 ->84 ->49 ->3 ->56 ->68 ->74 ->91 ->0
No. 3  : 0 ->63 ->36 ->47 ->21 ->38 ->87 ->50 ->43 ->26 ->35 ->48 ->77 ->80 ->0
No. 4  : 0 ->45 ->95 ->52 ->29 ->64 ->88 ->61 ->73 ->90 ->94 ->10 ->20 ->55 ->32 ->17 ->93 ->0
No. 5  : 0 ->5 ->33 ->28 ->82 ->98 ->12 ->19 ->22 ->97 ->96 ->25 ->89 ->70 ->100 ->0
No. 6  : 0 ->42 ->14 ->69 ->23 ->75 ->79 ->78 ->66 ->4 ->24 ->1 ->60 ->58 ->0
No. 7  : 0 ->92 ->59 ->39 ->83 ->15 ->11 ->99 ->40 ->18 ->81 ->34 ->13 ->0
No. 8  : 0 ->27 ->31 ->2 ->62 ->16 ->67 ->51 ->41 ->9 ->0
************************************************************
Check_Ans = 3523.7571791004384
Solution does not satisfy time windows construction
************************************************************
The Minimum Total Distance =  1774.7772089653195
Concrete Schedule of Each Route as Following : 
No. 1  : 0 ->40 ->53 ->0
No. 2  : 0 ->95 ->98 ->16 ->85 ->96 ->0
No. 3  : 0 ->36 ->47 ->19 ->49 ->48 ->0
No. 4  : 0 ->52 ->88 ->10 ->0
No. 5  : 0 ->6 ->0
No. 6  : 0 ->39 ->75 ->41 ->56 ->4 ->0
No. 7  : 0 ->82 ->18 ->84 ->17 ->93 ->0
No. 8  : 0 ->31 ->30 ->51 ->50 ->0
No. 9  : 0 ->63 ->64 ->90 ->66 ->1 ->0
No. 10  : 0 ->11 ->8 ->46 ->60 ->89 ->0
No. 11  : 0 ->14 ->44 ->38 ->43 ->58 ->0
No. 12  : 0 ->27 ->69 ->76 ->78 ->55 ->25 ->0
No. 13  : 0 ->59 ->5 ->61 ->86 ->37 ->91 ->100 ->0
No. 14  : 0 ->2 ->21 ->73 ->22 ->74 ->0
No. 15  : 0 ->28 ->12 ->26 ->0
No. 16  : 0 ->65 ->71 ->9 ->34 ->35 ->77 ->0
No. 17  : 0 ->62 ->7 ->0
No. 18  : 0 ->72 ->23 ->67 ->54 ->24 ->80 ->0
No. 19  : 0 ->92 ->42 ->15 ->87 ->57 ->97 ->13 ->0
No. 20  : 0 ->33 ->81 ->20 ->32 ->70 ->0
No. 21  : 0 ->45 ->83 ->99 ->94 ->0
No. 22  : 0 ->29 ->79 ->3 ->68 ->0
************************************************************
Check_Ans = 1774.7772089651903
Solution satisfies time windows construction

程序耗时：72.50237035751343s

Process finished with exit code 0

禁忌搜索算法（ Tabu Search ）吴浩博