一、带时间窗车辆路径优化问题（Vehicle Routing Problem with Time Window，VRPTW）

1.1 问题描述

VRPTW问题定义在一队车辆$\mathcal{V}$、一组客户$\mathcal{C}$和一个有向图$\mathcal{G}$上。图由$|\mathcal{C}|+2$个点组成，其中客户由$1,2,\cdots ,n$表示，车站由$0$（始发站）和$n+1$（终点站）表示。所有节点的集合可表示为 N = { 0 , 1 , 2 , ⋯ , n , n + 1 } \mathcal{N}=\{0,1,2,\cdots,n,n+1\} N={0,1,2,⋯,n,n+1}。所有边的集合$\mathcal{A}$表示了车站与客户之间，以及客户之间的有向连接。其中没有边从$0$点终止或者从$n+1$点开始。每一条弧$(i,j),i\ne j$都对应一个成本$c_{ij}$和时间$t_{ij}$，时间可以包括在弧上$(i,j)$的行驶时间和在$i$点的服务时间。所有车辆通常是同质化的，每辆车都存在容量上限$Q$，每一个客户点$i$都有需要被满足的需求$d_i$，并且需要在时间窗$[a_i,b_i]$之内被服务[1]。待优化的问题即为，如何决策车辆访问客户的路径，使得在满足一定约束的条件下，实现最小化总成本的目标。

1.2 模型构建

参数

• $\mathcal{V}$：车辆集合， V = { 0 , 1 , ⋯ , V } \mathcal{V}=\{0,1,\cdots, V\} V={0,1,⋯,V}；
• $\mathcal{C}$：客户集合， C = { 0 , 1 , ⋯ , n } \mathcal{C}=\{0,1,\cdots, n\} C={0,1,⋯,n}；
• $\mathcal{N}$：节点集合， N = { 0 , 1 , 2 , ⋯ , n , n + 1 } \mathcal{N}=\{0,1,2,\cdots,n,n+1\} N={0,1,2,⋯,n,n+1}；
• $c_{ij}$：从$i$到$j$的行驶距离；
• $d_i$：每一个客户点$i$都有需要被满足的需求量；
• $Q$：车辆容量；
• $t_{ij}$：从$i$点到$j$点的行驶时间，以及服务$i$点的时间之和；
• $[a_i,b_i]$：访问$i$点的时间窗；
• $M_{ij}$：足够大的正数。

决策变量

• $x_{ijk}$：车辆路径决策变量，$x_{ijk}=1$，车辆$k$经过弧$(i,j)$。
• $s_{ik}$：车辆$k$服务$i$点的开始时刻；

混合整数规划模型

min ⁡ ∑ k ∈ V ∑ i ∈ N ∑ j ∈ N c i j x i j k subject to ∑ i ∈ C d i ∑ j ∈ N x i j k ⩽ Q , ∀ k ∈ V , ∑ k ∈ V ∑ j ∈ N x i j k = 1 , ∀ i ∈ C ∑ j ∈ N x 0 j k = 1 , ∀ k ∈ V , ∑ i ∈ N x i h k − ∑ j ∈ N x h j k = 0 , ∀ h ∈ C , ∀ k ∈ V , ∑ i ∈ N x i , n + 1 , k = 1 ∀ k ∈ V , s i k + t i j − M i j ( 1 − x i j k ) ⩽ s j k , ∀ i , j ∈ N , ∀ k ∈ V , a i ⩽ s i k ⩽ b i , ∀ i ∈ N , ∀ k ∈ V ,  x i j k ∈ { 0 , 1 } , ∀ i , j ∈ N , ∀ k ∈ V s i k ≥ 0 , ∀ i ∈ N , ∀ k ∈ V \begin{align} \min \quad &\sum_{k \in \mathcal{V}} \sum_{i \in \mathcal{N}} \sum_{j \in \mathcal{N}} c_{i j} x_{i j k} \\ \text {subject to} \quad &\sum_{i \in \mathcal{C}} d_{i} \sum_{j \in \mathcal{N}} x_{i j k} \leqslant Q, \quad \forall k \in \mathcal{V}, \\ &\sum_{k \in \mathcal{V}} \sum_{j \in \mathcal{N}} x_{i j k}=1, \quad \forall i \in \mathcal{C} \\ &\sum_{j \in \mathcal{N}} x_{0 j k}=1, \quad \forall k \in \mathcal{V}, \\ &\sum_{i \in \mathcal{N}} x_{i h k}-\sum_{j \in \mathcal{N}} x_{h j k}=0 , \quad \forall h \in \mathcal{C}, \forall k \in \mathcal{V}, \\ &\sum_{i \in \mathcal{N}} x_{i, n+1, k}=1 \quad \forall k \in \mathcal{V}, \\ &s_{i k}+t_{i j}-M_{i j}\left(1-x_{i j k}\right) \leqslant s_{j k}, \quad \forall i, j \in \mathcal{N}, \forall k \in \mathcal{V}, \\ &a_{i} \leqslant s_{i k} \leqslant b_{i} , \quad \forall i \in \mathcal{N}, \forall k \in \mathcal{V} \text {, } \\ &x_{i j k} \in\{0,1\}, \quad \forall i, j \in \mathcal{N}, \forall k \in \mathcal{V} \\ &s_{i k} \geq 0, \quad \forall i \in \mathcal{N}, \forall k \in \mathcal{V} \\ \end{align} minsubject to​k∈V∑​i∈N∑​j∈N∑​cij​xijk​i∈C∑​di​j∈N∑​xijk​⩽Q,∀k∈V,k∈V∑​j∈N∑​xijk​=1,∀i∈Cj∈N∑​x0jk​=1,∀k∈V,i∈N∑​xihk​−j∈N∑​xhjk​=0,∀h∈C,∀k∈V,i∈N∑​xi,n+1,k​=1∀k∈V,sik​+tij​−Mij​(1−xijk​)⩽sjk​,∀i,j∈N,∀k∈V,ai​⩽sik​⩽bi​,∀i∈N,∀k∈V, xijk​∈{0,1},∀i,j∈N,∀k∈Vsik​≥0,∀i∈N,∀k∈V​​

下面是各个约束的具体含义：

• 目标函数(1)最小化总体行驶距离；
• 约束(2)说明车辆的载重量不能超过其容量上限；
• 约束(3)保证了每个客户点都被访问了一次；
• 约束(4-6)分别保证了每辆车必须从始发站出发，到达并离开每个客户点，并最终回到终点站；
• 约束(7)规定了车辆开始服务点的时刻，不能早于开始服务点的时刻，加上从点到 - 点的行驶时间以及服务点的时间之和，其中的取值下界为；
• 约束(8)使得开始服务点的时刻是在规定的时间窗范围内；
• 约束(9-10)是对于决策变量的约束。

二、使用Java调用Cplex求解VRPTW问题

在pom.xml中先导入依赖：

<dependency><groupId>commons-logging</groupId><artifactId>commons-logging</artifactId><version>1.2</version>
</dependency><dependency><groupId>org.apache.commons</groupId><artifactId>commons-lang3</artifactId><version>3.8</version>
</dependency><dependency><groupId>org.python</groupId><artifactId>jython-slim</artifactId><version>2.7.2</version>
</dependency>

2.1 完整代码

package org.example;import ilog.concert.*;
import ilog.cplex.IloCplex;import org.apache.commons.lang3.time.StopWatch;
import org.apache.commons.lang3.tuple.ImmutableTriple;
import org.apache.commons.logging.Log;
import org.apache.commons.logging.LogFactory;
import org.apache.commons.lang3.tuple.Triple;import java.io.*;
import java.util.*;public class Main {private static final Log log = LogFactory.getLog(Main.class);//【注】M=1e8合适public static final double M = 1e8;// public static final double M = Double.POSITIVE_INFINITY;// M = Double.MAX_VALUE;// M = Integer.MAX_VALUE;public static void main(String[] args) {String pathname = "src/main/resources/Solomon/100_customer/c101.txt";VRP problem = new VRP(pathname);System.out.println(problem);try {StopWatch watch = StopWatch.createStarted();//声明cplex优化模型IloCplex model = new IloCplex();model.setOut(null);// 定义决策变量x—ijkIloIntVar[][][] x = new IloIntVar[problem.nodeNum][problem.nodeNum][problem.vehicleNum];for (int i : problem.nodes.keySet()) {for (int j : problem.nodes.keySet()) {if (i != j && i != problem.endDepotId && j != 0) {for (int k = 0; k < problem.vehicleNum; k++) {x[i][j][k] = model.boolVar();}}}}// 定义决策变量s-ikIloNumVar[][] s = new IloNumVar[problem.nodeNum][problem.vehicleNum];for (int i : problem.nodes.keySet()) {for (int k = 0; k < problem.vehicleNum; k++) {s[i][k] = model.numVar(problem.nodes.get(problem.startDepotId).readyTime, problem.nodes.get(problem.endDepotId).dueDate, IloNumVarType.Float);}}// 目标函数（1）最小化车辆行驶距离IloLinearNumExpr objExpr = model.linearNumExpr();for (int i : problem.nodes.keySet()) {for (int j : problem.nodes.keySet()) {if (i != j && i != problem.endDepotId && j != 0) {for (int k = 0; k < problem.vehicleNum; k++) {objExpr.addTerm(problem.distance[i][j], x[i][j][k]);}}}}model.addMinimize(objExpr);// 约束（2）说明车辆的载重量不能超过其容量上限；int Q = problem.vehicleCapacity;for (int k = 0; k < problem.vehicleNum; k++) {IloLinearIntExpr expr = model.linearIntExpr();for (int i : problem.customers.keySet()) {int d = problem.customers.get(i).demand;for (int j : problem.nodes.keySet()) {if (i != j && i != problem.endDepotId && j != 0) {expr.addTerm(d, x[i][j][k]);}}}model.addLe(expr, Q);}// 约束（3）每个客户只能有一辆车服务一次（每个客户只能被访问一次）for (int i : problem.customers.keySet()) {IloLinearIntExpr expr = model.linearIntExpr();for (int j : problem.nodes.keySet()) {if (i != j && j != 0) {for (int k = 0; k < problem.vehicleNum; k++) {expr.addTerm(1, x[i][j][k]);}}}model.addEq(expr, 1);}// 约束（4）每辆车必须从depot出发（0）for (int k = 0; k < problem.vehicleNum; k++) {IloLinearIntExpr expr = model.linearIntExpr();for (int j : problem.nodes.keySet()) {if (j != 0) {expr.addTerm(1, x[0][j][k]);}}model.addEq(expr, 1);}// 约束（5）车辆服务j点后必须离开for (int k = 0; k < problem.vehicleNum; k++) {for (int h : problem.customers.keySet()) {IloLinearIntExpr expr1 = model.linearIntExpr();for (int i : problem.nodes.keySet()) {if (h != i && i != problem.endDepotId) {expr1.addTerm(1, x[i][h][k]);}}IloLinearIntExpr expr2 = model.linearIntExpr();for (int j : problem.nodes.keySet()) {if (h != j && j != 0) {expr2.addTerm(1, x[h][j][k]);}}model.addEq(expr1, expr2);}}// 约束（6）每辆车必须返回配送中心（n+1）for (int k = 0; k < problem.vehicleNum; k++) {IloLinearIntExpr expr = model.linearIntExpr();for (int i : problem.nodes.keySet()) {if (i != problem.endDepotId) {expr.addTerm(1, x[i][problem.endDepotId][k]);}}model.addEq(expr, 1);}// 约束（7）规定了车辆开始服务顾客的时刻for (int k = 0; k < problem.vehicleNum; k++) {for (int i : problem.nodes.keySet()) {for (int j : problem.nodes.keySet()) {if (i != j && i != problem.endDepotId && j != 0) {IloLinearNumExpr expr5 = model.linearNumExpr();double tij = problem.distance[i][j] / problem.vehicleVelocity;expr5.addTerm(1, s[i][k]);expr5.addTerm(-1, s[j][k]);expr5.addTerm(M, x[i][j][k]);model.addLe(expr5, M - tij - problem.nodes.get(i).serviceTime);}}}}// 约束（8）使得开始服务点的时刻是在规定的时间窗范围内；for (int k = 0; k < problem.vehicleNum; k++) {for (int i : problem.customers.keySet()) {IloLinearNumExpr expr = model.linearNumExpr();expr.addTerm(1, s[i][k]);model.addLe(expr, problem.customers.get(i).dueDate);model.addGe(expr, problem.customers.get(i).readyTime);}}//            model.exportModel("VRPTW.lp");if (!model.solve()) {System.out.println("求解失败");System.out.println("求解状态：" + model.getCplexStatus());} else {System.out.println("求解成功");System.out.println("求解状态：" + model.getCplexStatus());System.out.printf("目标函数值：%.2f\n", model.getObjValue());Map<Integer, ArrayList<Integer>> vehicle2path = new HashMap<>();for (int k = 0; k < problem.vehicleNum; k++) {vehicle2path.put(k, new ArrayList<>());}//模型可解，得到车辆路径for (int k = 0; k < problem.vehicleNum; k++) {boolean terminate = false;int i = problem.startDepotId;vehicle2path.get(k).add(i);while (!terminate) {for (int j : problem.nodes.keySet()) {if (i != j && j != 0 && model.getValue(x[i][j][k]) >= 0.5) {vehicle2path.get(k).add(j);i = j;break;}}if (i == problem.endDepotId) {terminate = true;}}if (vehicle2path.get(k).size() == 2) {vehicle2path.remove(k);}}for (Integer vehicle : vehicle2path.keySet()) {System.out.println(vehicle2path.get(vehicle));}System.out.println("使用车辆：" + vehicle2path.size());System.out.printf("程序运行时间：%ds", watch.getTime() / 1000);}} catch (Exception e) {//noinspection CallToPrintStackTracee.printStackTrace();log.error(e);}}
}class VRP {public final int vehicleNum;public final int vehicleCapacity;public final Map<Integer, Node> nodes;public final Map<Integer, Node> customers;public final Node startDepot;public final Node endDepot;public final int startDepotId;public final int endDepotId;public final int nodeNum;public final double[][] distance;public final double vehicleVelocity = 5.;public VRP(String pathname) {Triple<Integer, Integer, Map<Integer, Node>> txt = readFile(pathname);this.vehicleNum = txt.getLeft();this.vehicleCapacity = txt.getMiddle();this.nodes = new HashMap<>(txt.getRight());startDepotId = 0;endDepotId = nodes.size();startDepot = this.nodes.get(startDepotId);endDepot = new Node(startDepot.x, startDepot.y, startDepot.demand, startDepot.readyTime, startDepot.dueDate, startDepot.serviceTime);this.nodes.put(this.nodes.size(), endDepot);this.customers = new HashMap<>();for (int i = 1; i < this.nodes.size() - 1; i++) {customers.put(i, nodes.get(i));}this.nodeNum = this.nodes.size();distance = calcDistanceMatrix();}public double[][] calcDistanceMatrix() {double[][] distance = new double[this.nodes.size()][this.nodes.size()];for (int i = 0; i < this.nodes.size(); i++) {int xi = this.nodes.get(i).x;int yi = this.nodes.get(i).y;for (int j = 0; j < this.nodes.size(); j++) {int xj = this.nodes.get(j).x;int yj = this.nodes.get(j).y;distance[i][j] = Math.sqrt(Math.pow(xi - xj, 2) + Math.pow(yi - yj, 2));}}//距离矩阵满足三角关系for (int k = 0; k < nodeNum; k++) {for (int i = 0; i < nodeNum; i++) {for (int j = 0; j < nodeNum; j++) {if (distance[i][j] > distance[i][k] + distance[k][j]) {distance[i][j] = distance[i][k] + distance[k][j];}}}}return distance;}/*** @param pathname C101.txt文件路径* @return 车辆数量，车辆容量，顶点集合*/public static Triple<Integer, Integer, Map<Integer, Node>> readFile(String pathname) {Map<Integer, Node> nodeList = new HashMap<>();int vehicleNum;int vehicleCapacity;try (Scanner scanner = new Scanner(new FileReader(pathname))) {// 跳过4行for (int i = 0; i < 4; i++) {scanner.nextLine();}String line = scanner.nextLine().trim();String[] data = line.split("\\s+");vehicleNum = Integer.parseInt(data[0]);vehicleCapacity = Integer.parseInt(data[1]);// 跳过4行for (int i = 0; i < 4; i++) {scanner.nextLine();}while (scanner.hasNextLine()) {  //按行读取字符串line = scanner.nextLine();line = line.trim();data = line.split("\\s+");int id = Integer.parseInt(data[0]);int x = Integer.parseInt(data[1]);int y = Integer.parseInt(data[2]);int demand = Integer.parseInt(data[3]);int readyTime = Integer.parseInt(data[4]);int deliveryTime = Integer.parseInt(data[5]);int serviceTime = Integer.parseInt(data[6]);nodeList.put(id, new Node(x, y, demand, readyTime, deliveryTime, serviceTime));}} catch (FileNotFoundException e) {throw new RuntimeException(e);}return new ImmutableTriple<>(vehicleNum, vehicleCapacity, nodeList);}@Overridepublic String toString() {return "节点数=" + nodes.size() +", 顾客数=" + customers.size() +", 车辆数：" + vehicleNum +", 车辆容量：" + vehicleCapacity;}
}final class Node {public final int x;public final int y;public final int demand;public final int readyTime;public final int dueDate;public final int serviceTime;public Node(int x, int y, int demand, int readyTime, int dueDate, int serviceTime) {this.x = x;this.y = y;this.demand = demand;this.readyTime = readyTime;this.dueDate = dueDate;this.serviceTime = serviceTime;}@Overridepublic String toString() {return "Node[" +"x=" + x + ", " +"y=" + y + ", " +"demand=" + demand + ", " +"readyTime=" + readyTime + ", " +"dueDate=" + dueDate + ", " +"serviceTime=" + serviceTime + ']';}@Overridepublic boolean equals(Object obj) {if (obj == this) return true;if (obj == null || obj.getClass() != this.getClass()) return false;var that = (Node) obj;return this.x == that.x &&this.y == that.y &&this.demand == that.demand &&this.readyTime == that.readyTime &&this.dueDate == that.dueDate &&this.serviceTime == that.serviceTime;}@Overridepublic int hashCode() {return Objects.hash(x, y, demand, readyTime, dueDate, serviceTime);}}

2.2 求解结果

采用Solomon数据集c101.txt进行求解，

C101VEHICLE
NUMBER     CAPACITY25         200CUSTOMER
CUST NO.  XCOORD.   YCOORD.    DEMAND   READY TIME  DUE DATE   SERVICE   TIME0      40         50          0          0       1236          0   1      45         68         10        912        967         90   2      45         70         30        825        870         90   3      42         66         10         65        146         90   4      42         68         10        727        782         90   5      42         65         10         15         67         90   6      40         69         20        621        702         90   7      40         66         20        170        225         90   8      38         68         20        255        324         90   9      38         70         10        534        605         90   10      35         66         10        357        410         90   11      35         69         10        448        505         90   12      25         85         20        652        721         90   13      22         75         30         30         92         90   14      22         85         10        567        620         90   15      20         80         40        384        429         90   16      20         85         40        475        528         90   17      18         75         20         99        148         90   18      15         75         20        179        254         90   19      15         80         10        278        345         90   20      30         50         10         10         73         90   21      30         52         20        914        965         90   22      28         52         20        812        883         90   23      28         55         10        732        777         90   24      25         50         10         65        144         90   25      25         52         40        169        224         90   26      25         55         10        622        701         90   27      23         52         10        261        316         90   28      23         55         20        546        593         90   29      20         50         10        358        405         90   30      20         55         10        449        504         90   31      10         35         20        200        237         90   32      10         40         30         31        100         90   33       8         40         40         87        158         90   34       8         45         20        751        816         90   35       5         35         10        283        344         90   36       5         45         10        665        716         90   37       2         40         20        383        434         90   38       0         40         30        479        522         90   39       0         45         20        567        624         90   40      35         30         10        264        321         90   41      35         32         10        166        235         90   42      33         32         20         68        149         90   43      33         35         10         16         80         90   44      32         30         10        359        412         90   45      30         30         10        541        600         90   46      30         32         30        448        509         90   47      30         35         10       1054       1127         90   48      28         30         10        632        693         90   49      28         35         10       1001       1066         90   50      26         32         10        815        880         90   51      25         30         10        725        786         90   52      25         35         10        912        969         90   53      44          5         20        286        347         90   54      42         10         40        186        257         90   55      42         15         10         95        158         90   56      40          5         30        385        436         90   57      40         15         40         35         87         90   58      38          5         30        471        534         90   59      38         15         10        651        740         90   60      35          5         20        562        629         90   61      50         30         10        531        610         90   62      50         35         20        262        317         90   63      50         40         50        171        218         90   64      48         30         10        632        693         90   65      48         40         10         76        129         90   66      47         35         10        826        875         90   67      47         40         10         12         77         90   68      45         30         10        734        777         90   69      45         35         10        916        969         90   70      95         30         30        387        456         90   71      95         35         20        293        360         90   72      53         30         10        450        505         90   73      92         30         10        478        551         90   74      53         35         50        353        412         90   75      45         65         20        997       1068         90   76      90         35         10        203        260         90   77      88         30         10        574        643         90   78      88         35         20        109        170         90   79      87         30         10        668        731         90   80      85         25         10        769        820         90   81      85         35         30         47        124         90   82      75         55         20        369        420         90   83      72         55         10        265        338         90   84      70         58         20        458        523         90   85      68         60         30        555        612         90   86      66         55         10        173        238         90   87      65         55         20         85        144         90   88      65         60         30        645        708         90   89      63         58         10        737        802         90   90      60         55         10         20         84         90   91      60         60         10        836        889         90   92      67         85         20        368        441         90   93      65         85         40        475        518         90   94      65         82         10        285        336         90   95      62         80         30        196        239         90   96      60         80         10         95        156         90   97      60         85         30        561        622         90   98      58         75         20         30         84         90   99      55         80         10        743        820         90   100      55         85         20        647        726         90

求解结果如下：

"C:\Program Files (x86)\openjdk-17.0.2_windows-x64_bin\jdk-17.0.2\bin\java.exe" "-javaagent:C:\Program Files\JetBrains\IntelliJ IDEA Community Edition 2023.2\lib\idea_rt.jar=50414:C:\Program Files\JetBrains\IntelliJ IDEA Community Edition 2023.2\bin" -Dfile.encoding=UTF-8 -classpath "C:\Users\pengkangzhen\IdeaProjects\cplex-vrp\target\classes;C:\Users\pengkangzhen\.m2\repository\commons-logging\commons-logging\1.2\commons-logging-1.2.jar;C:\Users\pengkangzhen\.m2\repository\org\apache\commons\commons-lang3\3.8\commons-lang3-3.8.jar;C:\Program Files\IBM\ILOG\CPLEX_Studio1263\cplex\lib\cplex.jar" org.example.Main
节点数=102, 顾客数=100, 车辆数：25, 车辆容量：200
求解成功
求解状态：Optimal
目标函数值：828.94
[0, 98, 96, 95, 94, 92, 93, 97, 100, 99, 101]
[0, 67, 65, 63, 62, 74, 72, 61, 64, 68, 66, 69, 101]
[0, 5, 3, 7, 8, 10, 11, 9, 6, 4, 2, 1, 75, 101]
[0, 32, 33, 31, 35, 37, 38, 39, 36, 34, 101]
[0, 20, 24, 25, 27, 29, 30, 28, 26, 23, 22, 21, 101]
[0, 57, 55, 54, 53, 56, 58, 60, 59, 101]
[0, 81, 78, 76, 71, 70, 73, 77, 79, 80, 101]
[0, 43, 42, 41, 40, 44, 46, 45, 48, 51, 50, 52, 49, 47, 101]
[0, 90, 87, 86, 83, 82, 84, 85, 88, 89, 91, 101]
[0, 13, 17, 18, 19, 15, 16, 14, 12, 101]
使用车辆：10
程序运行时间：57s
Process finished with exit code 0

三、求解过程中踩的坑

在使用大M法的时候，务必注意M的取值：不能取太小，也不能取太大！取太小可能导致出现不可行解，取太大可能会因为计算机的精度问题导致约束失效。如果Double.MAX_VALUE或Double.POSITIVE_INFINITY会导致求解出现问题，得不到最优解。

Double.MAX_VALUE是double可以表示的最大值(大约$1.7\times 10^{308}$)。如果尝试减去数据类型的最大可能值，则这将导致一些计算问题。

public static void main(String[] args) {System.out.println(Double.MAX_VALUE);System.out.println(Double.POSITIVE_INFINITY);
}

1.7976931348623157E308
Infinity

参考

• VRPTW问题Solomon标准测试数据集
• https://blog.csdn.net/m0_45008035/article/details/134127916
• https://cloud.tencent.com/developer/article/1103529