一、问题描述
1.1 VRP（车辆路径规划）问题
在已知需求点位置（多为坐标）的情况下，通过计算获得到达这些需求点最短路径或是最小成本的线路。
1.2 VRPTW（带时间窗的路径规划）问题
图的左侧，有一个Depot（仓库、车场）节点和部分需求节点。其中：
Depot： 仓库（车场）、中心节点
N（Vehicle）：车辆数
C（Capacity）：车辆容量
需求点1：
（45,68）：需求点经纬度位置
[912,967]：需求点时间窗（最早到达时间、最晚到达时间）
S（90）：需求点服务时间
D（10）：需求点需求量
VRPTW问题是需要在仓库点Depot发出n辆车，车辆容量不得超出C（Capacity）的限制，在满足每一个需求点时间窗需求的前提下，找到最佳的运输车辆数以及最佳的运输路线（通常以距离最短作为目标，当然也有时间最短）。因此，我们可以得到相对应的目标函数和约束函数。
决策变量：两点之间的弧为0，1变量，判断i点是否到达j点
目标函数：成本最低（用距离替代，即距离最短）
约束条件：
- 容量约束
- 时间约束，包括时间窗约束和两点之间前后到达的时间关系约束（防止子回路出现的约束已在其中）
- 每个点只能访问一次的约束（每个点必须被访问）
- 起点约束，必须从起点出发（此处弧为单向）
- 终点约束，必须停在终点（此处弧为单向）
- 客户点不能停留的约束，即进入i点之后必须从i点出去（此处i点为客户点）（与点只能访问一次的约束相辅相成，即可成为每个客户点只能经过一次的约束）
解释：
wik表示第k辆车到达i点的时刻
约束（2）还表示每个点必须被访问（因为等于1，而不是小于等于1）
约束（6）表示前后到达时间之间的关系，同时消除了子回路（例如，0，1，0的子回路，则会有w0k<w1k<w0k，产生矛盾）
约束（7）表示i点的时间窗约束
约束（8）表示起始点与终点的时间窗约束
1.3 ESPPRC（带资源约束的基本最短路径规划）问题
ESPPRC问题是VRPTW问题的子问题，即VRPTW是多辆车的路径规划，而ESPPRC是里面一辆车时的路径规划，找出一辆车的可行路径。
此时cij为距离dij，也就是在容量约束、时间窗约束、每个客户点只能访问一次约束下求得起点到终点的最短路径。
其数学模型如下：
二、算法介绍
2.1 优化算法简介
我们常见的求解优化算法的方法有三类：
- 精确解算法
- 启发式算法
- 优化算法求解器
精确解算法： 指可求出最优解的算法。已提出的精确算法种类较多，有分支定界法、割平面法、整数规划算法和动态规划算法等。精确解可以求得优化算法的最优解，但是相对应其时间成本较大，当算例的规模较大时求解速率极低。
启发式算法： 通过某些特定的原理在一个解集空间中去搜寻相对最优的解，常见的启发式算法有遗传算法、粒子群算法等等。本文采用的模拟退火算法也是启发式算法中的一种。
优化算法求解器： 优化求解器产品是求解优化问题的专业设计软件，常见的求解器有Gurobi、CPLEX等。
2.2 模拟退火算法（SA）
2.2.1 简介
模拟退火算法(SA)来源于固体退火原理，是一种基于概率的算法。将固体加温至充分高的温度，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，分子和原子越不稳定。而徐徐冷却时粒子渐趋有序，能量减少，原子越稳定。在冷却（降温）过程中，固体在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。
模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。
2.2.2 原理
模拟退火算法包含两个部分即Metropolis准则和退火过程，分别对应内循环和外循环。外循环就是退火过程，将固体达到较高的温度（初始温度T（0）），然后按照降温系数alpha使温度按照一定的比例下降，当达到终止温度Tf时，冷却结束，即退火过程结束。
Metropolis准则是内循环，即在每次温度下，迭代L次，寻找在该温度下能量的最小值（即最优解）。下图中所示即为在一次温度下，跌代L次，固体能量发生的变化。在该温度下，整个迭代过程中温度不发生变化，能量发生变化，当前一个状态x(n)的能量大于后一个状态x(n+1)的能量时，状态x(n)的解没有状态x(n+1)的解好，所以接受状态x(n+1)。但是如果下一状态的能量比前一个状态的能量高时，该不该接受下一状态呢？在这里设置一个接受概率P,即如果下一状态的能量比前一个状态的能量高，则接受下一状态的概率为P,下面具体讲一下如何接受下一个状态。
Metropolis准则就是如何在局部最优解的情况下让其跳出来（如图中B、C、E为局部最优），是退火的基础。1953年Metropolis提出重要性采样方法，即以概率来接受新状态，而不是使用完全确定的规则，称为Metropolis准则，计算量较低。
假设开始状态在A，多次迭代之后更新到B的局部最优解，这时发现更新到B时，能力比A要低，则说明接近最优解了，因此百分百转移，状态到达B后，发现下一步能量上升了，如果是梯度下降则是不允许继续向前的，而这里会以一定的概率跳出这个坑，这各概率和当前的状态、能量等都有关系。所以说这个概率的设计是很重要的，下面从数学方面进行解释。
假设前一个状态为x(n),系统根据某一指标（梯度下降，上节的能量），状态变为x(n+1),相应的，系统的能量由E(n)变为E(n+1),定义系统由x(n)变为x(n+1)的接受概率P为：
从上式我们可以看到，如果能量减小了，那么这种转移就被接受（概率为1），如果能量增大了，就说明系统偏离全局最优值位置更远了，此时算法不会立刻将其抛弃，而是进行概率操作：首先在区间[0,1]产生一个均匀分布的随机数ϵ ，如果ϵ<P，则此种转移接受，否则拒绝转移，进入下一步，往复循环。其中P以能量的变化量和T进行决定概率P的大小，所以这个值是动态的。
用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件Tf。而温度的作用就是来计算转移概率P的。当温度每次下降后，转移概率也发生变化，因此在所有温度下迭代L次的结果也都是不相同的。在每个温度下迭代L次来寻找当前温度下的最优解，然后降低温度继续寻找，直到到达终止温度，即转移概率P接近于0。
接受状态的三条原则：
- 在固定温度下，接受使目标函数下降的候选解的概率要大于使目标函数上升的候选解概率；
- 随着温度的下降，接受使目标函数上升的解的概率要逐渐减小；
- 当温度趋于零时，只能接受目标函数下降的解。
2.2.3 退火过程中的参数控制
- 初始的温度T(0)应选的足够高，使的所有转移状态都被接受。初始温度越高，获得高质量的解的概率越大，耗费的时间越长。
- 退火速率，即温度下降，最简单的下降方式是指数式下降：
T(n) = α T(n) ,n =1,2,3,…
其中α 是小于1的正数，一般取值为0.8到0.99之间。使的对每一温度，有足够的转移尝试，指数式下降的收敛速度比较慢。 - 终止温度
如果温度下降到终止温度或者达到用户设定的阈值，则退火完成。
2.2.4 算法步骤
1.模拟退火的基本思想:
(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)，每个T值的迭代次数L
(2) 对k=1, …, L做第(3)至第6步：
(3) 产生新解S′
(4) 计算增量ΔT=C(S′)-C(S)，其中C(S)为目标函数,C(S)相当于能量
(5) 若ΔT<0则接受S′作为新的当前解，否则以概率exp(-ΔT/T)接受S′作为新的当前解.
(6) 如果满足终止条件则输出当前解作为最优解，结束程序。
(7) T逐渐减少，且T->0，然后转第2步。
2.模拟退火算法新解的产生和接受可分为如下四个步骤：
第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若ΔT<0则接受S′作为新的当前解S，否则以概率P接受S′作为新的当前解S。
第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。
三、问题实现
3.1 程序结构
3.1.1 主方法运行
- 计算程序运行时间
- 读取文件数据
- 初始化参数及问题初始条件
- 测试数据集
- 模拟退火算法得解
- 检验解结果
3.1.2 读取文件数据
IO流
3.1.3 初始化参数及问题初始条件
总点数、坐标、时间窗、需求量、服务时间数据、距离矩阵、无向图构建
3.1.4 测试数据（可有可无）
判断时间与需求数据等是否有误
3.1.5 模拟退火算法
- 生成初始解
- 运用算子产生新解
- 评价函数计算新解的目标函数值
- Metropolis准则判断是否取该新解
- 退火过程，循环终止得到循环过程中所产生的最优解
（启发式算法流程很简单，重要的是评价函数和生成新解的算子）
3.1.6 检验解结果
测试解的某些数据是否对应算例的数据，该解测试其车容量、时间窗是否满满足要求
注意两个整数值转成double值进行比较大小时，有精度误差，需要加上或者减去一定大的小数
3.2 算例（Solomon经典算例C101(100))
VEHICLE
NUMBER CAPACITY
25 200
CUSTOMER
CUST NO. XCOORD. YCOORD. DEMAND READY TIME DUE DATE SERVICE TIME
0 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
3.3 java程序源码
package test;
import java.io.File;
import java.io.FileInputStream;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;
import java.util.Random;
/**
* ClassName: SA_ESPPRC
* Package: test
* Description:模拟退火算法求解ESPPRC问题
*
* @Author:
* @Create: 2023/6/13 - 10:35
* @Version:
*/
public class SA_ESPPRC {
int locationNum;//地点总数
int vehicleNum;//车子总数
double[] locationData;//地点数据[x,y,d,a,b,s]，坐标+需求+时间窗+服务时间（读取数据用）
List<double[]> locationList;//存储点数据的集合（读取数据用）
double[][] location;//存储每个点的坐标[x,y]
int[][] arcs;//弧,arcs[i][j]表示i到j点的弧，值为1则此弧存在，值为0则此弧不存在
double[][] distance;//两点之间的距离矩阵，满足三角关系，暂用距离表示花费c[i][j]=distance[i][j]
double[] a;//时间窗[ai,bi]
double[] b;//时间窗[ai,bi]
double E;//仓库点（起始点、终点）时间窗[E,l]
double L;//仓库点（起始点、终点）时间窗[E,l]
double C;//车子的总容量
double CMin;//车子最低载货量
double[] d;//每个点的需求di
double[] s;//每个点的服务时间
int MAX_GEN;//模拟退火算法内循环最大迭代次数，提高该次数可以提高解的质量，即更好地找到全局最优解，但会增加求解时间
List<Integer> temGh;//存放临时编码，即存放临时的一条路径
List<Integer> currGh;//存放当前编码，即存放当前所得的一条路径
List<Integer> bestGh;//最优编码，即最优路径解
double tempEvaluation;//临时解
double currEvaluation;//当前解
double bestEvaluation;//最优解
double T;//模拟退火时的初始温度
double Tf;//模拟退火时的外循环终止温度
double α;//降温速度，速度越快，越不容易摆脱局部最优解，速度慢则运算时间长
Random random;//随机函数对象
//读取数据
public List<double[]> dataInput(String path) {
File file;
FileInputStream fileInputStream = null;
StringBuilder sb = null;
try {
file = new File(path);//创建文件对象
fileInputStream = new FileInputStream(file);//创建字节输入流对象，导入文件
//读入数据
int len;
sb = new StringBuilder();//创建该对象，将文本内容转为字符串存储到该对象中
while ((len = fileInputStream.read()) != -1) {
sb.append((char) len);//append()方法将一个个(char)len字符存储字符串中
}
} catch (IOException e) {
e.printStackTrace();
} finally {
try {
if (fileInputStream != null) {
fileInputStream.close();//关闭数据
}
} catch (IOException e) {
e.printStackTrace();
}
}
//将读取到的数据传到List<double[]> locationList中
String[] data = sb.toString().split("\n");//将文本内容按行存储进字符串数组中
String[] data1 = data[2].split("\\s+");//将第三行的数据按空格分开，存储进字符串数组中
vehicleNum = Integer.parseInt(data1[1]);
C = Integer.parseInt(data1[2]);//每辆车总容量
locationData = new double[6];//数据数组初始化
locationList = new ArrayList<>();//数据集合初始化
for (int i = 7; i < data.length; i++) {//将起始点与100个客户点数据传输到数组和集合里
String[] data2 = data[i].split("\\s+");
locationData[0] = Double.parseDouble(data2[2]);//坐标x
locationData[1] = Double.parseDouble(data2[3]);//坐标y
locationData[2] = Double.parseDouble(data2[4]);//需求d
locationData[3] = Double.parseDouble(data2[5]);//时间窗a
locationData[4] = Double.parseDouble(data2[6]);//时间窗b
locationData[5] = Double.parseDouble(data2[7]);//服务时间s
locationList.add(locationData.clone());
}
//将终点传输到数组和集合里，终点数据即起点数据，只是编号不一样而已
String[] data2 = data[7].split("\\s+");
locationData[0] = Double.parseDouble(data2[2]);//坐标x
locationData[1] = Double.parseDouble(data2[3]);//坐标y
locationData[2] = Double.parseDouble(data2[4]);//需求d
locationData[3] = Double.parseDouble(data2[5]);//时间窗a
locationData[4] = Double.parseDouble(data2[6]);//时间窗b
locationData[5] = Double.parseDouble(data2[7]);//服务时间s
locationList.add(locationData.clone());
return locationList;//返回数据集合，包含起点、终点、客户点，共n+2各点
}
//初始化参数及问题初始条件
public void init() {
locationNum = locationList.size();//总点数
//坐标、时间窗、需求量、服务时间数据
location = new double[locationNum][2];//数组初始化
a = new double[locationNum];
b = new double[locationNum];
d = new double[locationNum];
s = new double[locationNum];
for (int i = 0; i < locationNum; i++) {//取数据
location[i][0] = locationList.get(i)[0];//x坐标
location[i][1] = locationList.get(i)[1];//y坐标
a[i] = locationList.get(i)[3];
b[i] = locationList.get(i)[4];
d[i] = locationList.get(i)[2];
s[i] = locationList.get(i)[5];
}
//仓库点的时间窗，即起点和终点的时间窗
E = a[0];
L = b[0];
//车子最低载货量
double demandSum = 0.0;
for (int i = 0; i < locationNum; i++) {
demandSum = demandSum + d[i];
}
CMin = demandSum / vehicleNum;
//两点之间距离的获取
distance = new double[locationNum][locationNum];
for (int i = 0; i < locationNum; i++) {
for (int j = i; j < locationNum; j++) {
if (i == j) {
distance[i][j] = Double.MAX_VALUE;//不能自己到达自己，设置距离为无限大
} else {//欧式距离
distance[i][j] = Math.sqrt(Math.pow(location[i][0] - location[j][0], 2) + Math.pow(location[i][1] - location[j][1], 2));
distance[j][i] = distance[i][j];
}
}
}
//三点之间的距离满足三角关系，即a+b<c
for (int k = 0; k < locationNum; k++) {
for (int i = 0; i < locationNum; i++) {
for (int j = 0; j < locationNum; j++) {
if (distance[i][j] > distance[i][k] + distance[k][j]) {
distance[i][j] = distance[i][k] + distance[k][j];
}
}
}
}
//G=(V,A)无向图的构建，去除掉某些不符合原始数据条件的弧和规定了起点和终点是单向弧后所得的图
arcs = new int[locationNum][locationNum];//无向图中的弧
for (int i = 0; i < locationNum; i++) {//初始化完全图
for (int j = 0; j < locationNum; j++) {
if (i != j) {//其余i==j的弧值默认为0（数组默认值）
arcs[i][j] = 1;
}
}
}
for (int i = 0; i < locationNum; i++) {//除掉某些不符合时间数据和容量数据的弧
for (int j = 0; j < locationNum; j++) {
if (i != j) {
if (a[i] + s[i] + distance[i][j] > b[j] + 0.0001) {//不能晚到
arcs[i][j] = 0;
}
if (d[i] + d[j] > C) {
arcs[i][j] = 0;
}
}
}
}
//起点和终点是单向弧，起点只能出不能进，终点只能进不能出
for (int i = 0; i < locationNum; i++) {
arcs[i][0] = 0;
arcs[locationNum - 1][i] = 0;
}
arcs[0][locationNum - 1] = 0;//起点不能直接到达终点
//参数初始化
MAX_GEN = 100;
T = 100;
Tf = 0.01;
α = 0.9;
random = new Random(System.currentTimeMillis());//new Random(随机数种子)，种子不同则在一次运行时所得的随机数不同，种子相同时则在一次运行中的几个随机数会出现相同现象。
temGh = new ArrayList<>();
currGh = new ArrayList<>();
bestGh = new ArrayList<>();
//打印参数
System.out.println("算例客户点总数：" + (locationNum - 2));
}
//测试数据集是否有误
public void testData() {
for (int i = 0; i < locationNum; i++) {
for (int j = 0; j < locationNum; j++) {
if (i != j) {
if (a[0] + distance[0][i] + s[i] + distance[i][locationNum - 1] > b[locationNum - 1]) {
//从起始点出发，最早到达i点，再到达终点时，判断其所花费的时间是否满足整个过程的时间窗，不满足则该数据集有误
System.out.println("The calculating example is false!");
}
}
}
}
//最小值min1=最晚到达i点的时间（即时间窗bi) - 起始点到i点所花费的时间（距离），如果min1比整个过程的时间窗E还小，则数据集有误
//最小值min2=最早到达i点的时间（即时间窗ai）+ i点的服务时间 + i点到终点所花费的时间（距离），如果min2比整个过程的时间窗L还大，则数据集有误
double min1 = Double.MAX_VALUE;
double min2 = Double.MAX_VALUE;
for (int i = 1; i < locationNum - 1; i++) {
if (b[i] - distance[0][i] < min1) {
min1 = b[i] - distance[0][i];
}
if (a[i] + s[i] + distance[i][locationNum - 1] < min2) {
min2 = a[i] + s[i] + distance[i][locationNum - 1];
}
}
if (E > min1 || L < min2) {
System.out.println("Duration false!");
System.exit(0);
}
//如果一个地点的需求量超过一辆车的容量，则数据集有误
for (int i = 0; i < locationNum; i++) {
if (d[i] > C) {
System.out.println("Demand false!");
}
}
}
//评价函数（计算目标函数值）
public double evaluate(List<Integer> path) {
double pathLen = 0.0;//路径长度
for (int i = 0; i < path.size() - 1; i++) {
pathLen = pathLen + distance[path.get(i)][path.get(i + 1)];//整条路径的路程累加和
}
return pathLen;
}
//检验是否满足约束条件(路径弧、容量和时间窗）
public boolean checkConstrains(List<Integer> currPath) {
//判断生成的路径各点之间是否可连接
for (int i = 0; i < currPath.size() - 1; i++) {
if (arcs[currPath.get(i)][currPath.get(i + 1)] == 0) {
return false;
}
}
//车子容量约束
double capacitySum = 0.0;//车子载货容量累加和
for (int i : currPath) {
capacitySum = capacitySum + d[i];
}
if (capacitySum > C + 0.0001) {
return false;
}
if (currPath.contains(locationNum - 1)) {//一整条路径的车载量不小于CMin
if (capacitySum < CMin - 0.0001) {
return false;
}
}
//车子时间窗约束
double[] arriveTime = new double[currPath.size()];//车子到达j点的时刻
arriveTime[0] = 0.0;//车子在起点的出发时刻设定为0
for (int i = 1; i < currPath.size(); i++) {
arriveTime[i] = arriveTime[i - 1] + s[currPath.get(i - 1)] + distance[currPath.get(i - 1)][currPath.get(i)];
if (arriveTime[i] < a[currPath.get(i)]) {
arriveTime[i] = a[currPath.get(i)];
}
if (arriveTime[i] > b[currPath.get(i)] + 0.0001) {//不符合硬时间窗约束：可以早到等待，但不能迟到
return false;
}
}
return true;
}
// //针对初始解的产生而设置，只检验弧约束和时间窗约束（编码方式2所用）
// public boolean checkInitialConstrains(List<Integer> currPath) {
// //判断生成的路径各点之间是否可连接
// for (int i = 0; i < currPath.size() - 1; i++) {
// if (arcs[currPath.get(i)][currPath.get(i + 1)] == 0) {
// return false;
// }
// }
// //车子时间窗约束
// double[] arriveTime = new double[currPath.size()];//车子到达j点的时刻
// arriveTime[0] = 0.0;//车子在起点的出发时刻设定为0
// for (int i = 1; i < currPath.size(); i++) {
// arriveTime[i] = arriveTime[i - 1] + s[currPath.get(i - 1)] + distance[currPath.get(i - 1)][currPath.get(i)];
// if (arriveTime[i] < a[currPath.get(i)]) {
// arriveTime[i] = a[currPath.get(i)];
// }
// if (arriveTime[i] > b[currPath.get(i)] + 0.0001) {//不符合硬时间窗约束：可以早到等待，但不能迟到
// return false;
// }
// }
// return true;
// }
//随机生成一条除起点和终点外的序列，序列中没有重复的点
public List<Integer> initialSequence(){
List<Integer> list = new ArrayList<>();
while (true){
int customer = random.nextInt(locationNum - 2) + 1;
if (!list.contains(customer)){//序列里不重复的客户点
list.add(customer);
}
if (list.size() == locationNum - 2){//序列长度为客户点总数
break;
}
}
return list;
}
// //通过位置元素交换的方式，生成一条新的序列
// public List<Integer> newSequence(List<Integer> oldList){
// List<Integer> newList = new ArrayList<>();
// newList.addAll(oldList);
// int index1 = random.nextInt(newList.size());//随机生成序列中的一个位置索引
// int index2 = random.nextInt(newList.size());//随机生成序列中的一个位置索引
// //位置交换
// int temp = newList.get(index1);
// newList.set(index1,newList.get(index2));
// newList.set(index2,temp);
//
// return newList;
// }
//随机生成一条序列，在序列中取点进行编码，以此生成初始解
public void initialSolution() {
List<Integer> sequence = initialSequence();//随机生成除起点和终点外的初始序列
//编码方式1：直接取序列的点然后做约束判断，不满足约束则产生新序列（方式简单粗暴，代码简单，但是程序运行时间时长时短，依解长度的随机数而定，不稳定）
while (true){
int solutionLength = random.nextInt(locationNum - 2) + 1;
currGh.add(0);
for (int i = 0; i < solutionLength; i++) {
currGh.add(sequence.get(i));
}
// //编码方式2：在序列中依次取满足约束条件的点（会相对复杂些，但是减少了程序总的循环次数，使得程序运行时间较短且稳定
// while (true){
// currGh.add(0);//先加入起点
// int index = 1;//路径解的索引
// boolean flag1;//用于弧约束和时间窗约束的判断
boolean flag2;//用于弧、容量、时间窗约束的判断
// for (int i = 0; i < sequence.size(); i++) {
// //解路径依次加入序列中的客户点，不满足弧和时间窗约束的点则被跳过，直到满足容量约束为止，若整条序列找不到一个解路径，则产生新的序列继续生成
// currGh.add(sequence.get(i));
// flag1 = checkInitialConstrains(currGh);
// if (flag1){
// index++;//满足弧和时间窗约束则加入到解路径中
// //车子容量约束
// double capacitySum = 0.0;//车子载货容量累加和
// for (int j : currGh) {
// capacitySum = capacitySum + d[j];
// }
// if (capacitySum <= C + 0.0001 && capacitySum >= CMin - 0.0001) {
// break;
// }
// } else {
// currGh.remove(index);//不满足约束则删除该点
// }
// }
currGh.add(locationNum - 1);//加入终点
flag2 = checkConstrains(currGh);//遍历完整条序列后可能找不到一条合法的解路径，需要进行校验
if (flag2){
break;//满足约束则退出整个循环
} else {
currGh.removeAll(currGh);//不满足约束则删除之前加入的客户点
//生成序列有两种方式，第一种是重新随机生成序列，另一种是运用元素交换法在原始序列的基础上得出新序列
//随机方法使得程序运行时间时长时短，但代码简单粗暴；元素交换法能够减少循环次数，降低程序运行时间
sequence = initialSequence();//重新生成一条序列
// sequence = newSequence(sequence);//元素交换法生成一条新的序列
}
}
//复制当前路径解
temGh.addAll(currGh);
bestGh.addAll(currGh);
currEvaluation = evaluate(currGh);
tempEvaluation = currEvaluation;
bestEvaluation = currEvaluation;
//打印初始值
System.out.println("初始解为：" + currGh);
double demandSum = 0.0;
for (int i : currGh) {
demandSum = demandSum + d[i];
}
System.out.println("初始解车子载货量为：" + demandSum);
System.out.println("初始解路径长度为：" + currEvaluation);
}
//算子，remove算子、change算子和insert算子，轮盘赌法选择算子
//remove算子：随机选择解路径中一个点，移出该点
//change算子：随机选择一个客户点，随机选择解路径中一个点进行更换
//insert算子：随机选择一个客户点，并随机选择解路径中的一个点，客户点随机插入其左边或者右边
public List<Integer> generateNewGh(List<Integer> currPath) {//输入当前解
List<Integer> newPath = new ArrayList<>();//新解
newPath.addAll(currPath);
while (true) {//一定产生一个新的满足约束条件的解出来
boolean flag;
while (true){//一定产生一个新的解出来
double r = random.nextDouble();//产生[0,1)范围的随机数
int index = random.nextInt(newPath.size() - 2) + 1;//随机数，排除掉路径的起点和终点
if (r < 0.33) {//remove算子
if (currPath.size() > 3) {
newPath.remove(index);//解路径中起点和终点是能够保证不被删掉的，但是删除直到只剩下三个时，则不可再删除
break;
}
} else if (0.33 <= r && r < 0.66){//change算子
int customer = random.nextInt(locationNum - 2) + 1;//随机客户点，排除掉起点和终点
if (!currPath.contains(customer)){
newPath.set(index,customer);
break;
}
} else {//insert算子
double r1 = random.nextDouble();
int customer = random.nextInt(locationNum - 2) + 1;//随机客户点，排除掉起点和终点
if (!currPath.contains(customer)){
if (r1 < 0.5){//产生随机数，如果随机数小于0.5，则插入点的左边
newPath.add(index,customer);
break;
} else {//大于0.5时，插入点的右边
newPath.add(index + 1,customer);
break;
}
}
}
}
flag = checkConstrains(newPath);
if (flag) {
break;
} else {
newPath.removeAll(newPath);//此时集合的拷贝不能直接=，这样子会使两个集合同一地址，同时改变
newPath.addAll(currPath);
}
}
return newPath;//返回新解
}
//模拟退火核心算法
public void solve () {
initialSolution();//生成初始解
while (T > Tf) {//外层循环，根据温度进行退火，直到一定温度要求则停止
for (int i = 0; i < MAX_GEN; i++) {//内层循环，Metropolis准则
temGh.removeAll(temGh);
temGh.addAll(generateNewGh(currGh));//生成新解
tempEvaluation = evaluate(temGh);
if (tempEvaluation < bestEvaluation) {//新解优于最优解
bestEvaluation = tempEvaluation;
bestGh.removeAll(bestGh);
bestGh.addAll(temGh);
} else if (tempEvaluation < currEvaluation) {//新解比最优解差，但是比当前解好，则保留下来，用此新解产生另一个新解
currEvaluation = tempEvaluation;
currGh.removeAll(currGh);
currGh.addAll(temGh);
} else {//如果新解比最优解和当前解都要差，则采用Metropolis准则决定是否保留
double r = random.nextDouble();
if (r <= Math.exp(-1 * (Math.abs(tempEvaluation - currEvaluation) / T))) {
currEvaluation = tempEvaluation;
currGh.removeAll(currGh);
currGh.addAll(temGh);
}
}
}
T = T * α;//降温操作
}
//打印结果
double demandSum = 0.0;
for (int i : bestGh) {
demandSum = demandSum + d[i];
}
System.out.println("最短路径为：" + bestGh);
System.out.println("最短路径车子载货量为：" + demandSum);
System.out.println("最短路径长度为：" + bestEvaluation);
}
//检验解的结果
public void checkSolution () {
//解中车的车辆载荷是否超过车容量的判断
double capacity = 0;
for (int i : bestGh) {
capacity += d[i];
}
if (capacity > C) {
System.out.println("error:Capacity!!!");
System.exit(0);
}
if (capacity < CMin) {
System.out.println("可行路径，该路径的车载量小于车载量平均水平");
}
//解中时间窗可行性判断
double[] arriveTime = new double[bestGh.size()];//车子到达j点的时刻
arriveTime[0] = 0.0;//车子在起点的出发时刻设定为0
for (int i = 1; i < bestGh.size(); i++) {
arriveTime[i] = arriveTime[i - 1] + s[bestGh.get(i - 1)] + distance[bestGh.get(i - 1)][bestGh.get(i)];
if (arriveTime[i] < a[bestGh.get(i)]) {//早到则等待
arriveTime[i] = a[bestGh.get(i)];
}
}
for (int i = 0; i < bestGh.size() - 1; i++) {
int origin = bestGh.get(i);
int destination = bestGh.get(i + 1);
double originTime = arriveTime[i];
double destinationTime = arriveTime[i + 1];
if (originTime < a[origin] && destinationTime > b[destination]) {//到达时间不在时间窗范围内
System.out.println("Error:arriveTime!!!");
System.exit(0);
}
if (originTime + s[origin] + distance[origin][destination] > b[destination] + 0.0001) {//从i点到j点，到达j点时超过了j点的时间窗
System.out.println(origin + "->" + destination);
System.out.println("Error:be late");
System.exit(0);
}
if (destinationTime - s[origin] - distance[origin][destination] < a[origin] - 0.0001) {//double精度问题，两个整数值相等，换算成double值时不一定相等，容易造成误差，故在两个double值比较大小时，得加或者减去一定大的小数以消除精度差
double temp = destinationTime - s[origin] - distance[origin][destination];
System.out.println(origin + "->" + destination + ":" + temp);
System.out.println("Error");
System.exit(0);
}
}
}
//主方法运行
public static void main (String[]args){
long startTime = System.currentTimeMillis();//开始时间（毫秒），该时间返回的是long值
SA_ESPPRC sa_espprc = new SA_ESPPRC();
sa_espprc.dataInput("D:\\IDEA2022\\IDEA projects\\algorithm implementation\\CplexTest\\VRPTWDATAC101(100).txt");
sa_espprc.init();//初始化参数
sa_espprc.testData();//数据检测
sa_espprc.solve();//模拟退火算法
sa_espprc.checkSolution();//检查求解结果
System.out.println("程序运行时间：" + (System.currentTimeMillis() - startTime) / 1000d + "s");//所花费的时间，1000d指1000是double值，也就是1000.0
}
}
3.4 结果展示
1376
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
