一、带时间窗车辆路径问题描述
二、VRPTW数学模型
三、基于Cplex的VRPTW问题建模求解(JavaAPI)
package vrptw;
import java.io.BufferedReader;
import java.io.FileReader;
import java.util.Scanner;
import static net.mindview.util.Print.*;
//定义参数
class Data {
static int vertexNum; //所有点编号(包括配送中心和客户点,首尾(0和n)为配送中心)
static double E; //配送中心时间窗开始时间
static double L; //配送中心时间窗结束时间
static int vehNum; //车辆数
static double capacity; //车辆载荷
static int[][] xy; //所有点的坐标x,y
static int[] demand; //需求量
static double[] a; //时间窗开始时间【a[i],b[i]】
static double[] b; //时间窗结束时间【a[i],b[i]】
static double[] s; //客户点的服务时间
static int[][] arc; //arcs[i][j]表示i到j点的弧
static double[][] dist; //距离矩阵,满足三角关系,暂用距离表示花费 C[i][j]=dist[i][j]
//函数功能:从txt文件中读取数据并初始化参数
public static void readFile(String path, int vetexnum) throws Exception {
String line;
String[] substr;
Scanner cin = new Scanner(new BufferedReader(new FileReader(path))); //读取文件
for (int i = 0; i < 4; i++) {
line = cin.nextLine(); //读取一行
}
line = cin.nextLine();
line.trim(); //返回调用字符串对象的一个副本,删除起始和结尾的空格
substr = line.split(("\\s+")); //以空格为标志将字符串拆分
//初始化参数
vertexNum = vetexnum;
vehNum = Integer.parseInt(substr[1]);
capacity = Integer.parseInt(substr[2]);
xy = new int[vertexNum][2]; //所有点的坐标x,y
demand = new int[vertexNum]; //需求量
a = new double[vertexNum]; //时间窗开始时间
b = new double[vertexNum]; //时间窗结束时间
s = new double[vertexNum]; //服务时间
arc = new int[vertexNum][vertexNum];
dist = new double[vertexNum][vertexNum]; //距离矩阵,满足三角关系,用距离表示cost
for (int i = 0; i < 4; i++) {
line = cin.nextLine();
}
//读取vetexnum-1行数据
for (int i = 0; i < vertexNum - 1; i++) {
line = cin.nextLine();
line.trim();
substr = line.split("\\s+");
xy[i][0] = Integer.parseInt(substr[2]);
xy[i][1] = Integer.parseInt(substr[3]);
demand[i] = Integer.parseInt(substr[4]);
a[i] = Integer.parseInt(substr[5]);
b[i] = Integer.parseInt(substr[6]);
s[i] = Integer.parseInt(substr[7]);
}
cin.close();//关闭流
//初始化配送中心参数
xy[vertexNum - 1] = xy[0];
demand[vertexNum - 1] = 0;
a[vertexNum - 1] = a[0]; //时间窗开始时间【a[i],b[i]】
b[vertexNum - 1] = b[0]; //时间窗结束时间【a[i],b[i]】
E = a[0]; //配送中心时间窗开始时间
L = b[0]; //配送中心时间窗结束时间
s[vertexNum - 1] = 0;
}
public static void initialize() {
//距离矩阵初始化
for (int i = 0; i < vertexNum; i++) {
for (int j = 0; j < vertexNum; j++) {
if (i == j) {
dist[i][j] = 0;
continue;
}
dist[i][j] = Math.sqrt((xy[i][0] - xy[j][0]) * (xy[i][0] - xy[j][0]) +
(xy[i][1] - xy[j][1]) * (xy[i][1] - xy[j][1]));
dist[i][j] = UtilMeth.doubleTruncate(dist[i][j]);
}
}
dist[0][vertexNum - 1] = 0;
dist[vertexNum - 1][0] = 0;
//距离矩阵满足三角关系
for (int k = 0; k < vertexNum; k++)
for (int i = 0; i < vertexNum; i++)
for (int j = 0; j < vertexNum; j++)
if (dist[i][j] > dist[i][k] + dist[k][j])
dist[i][j] = dist[i][k] + dist[k][j];
double min1 = 1e15;
double min2 = 1e15;
//初始化为完全图
for (int i = 0; i < vertexNum; i++) {
for (int j = 0; j < vertexNum; j++) {
if (i != j) arc[i][j] = 1;
else arc[i][j] = 0;
}
}
//除去不符合时间窗和容量约束的边
for (int i = 0; i < vertexNum; i++) {
for (int j = 0; j < vertexNum; j++) {
if (i == j) {
continue;
}
if (a[i] + s[i] + dist[i][j] > b[j] || demand[i] + demand[j] > capacity) {
arc[i][j] = 0;
}
if (a[0] + dist[0][i] + s[i] + dist[i][vertexNum - 1] > b[vertexNum - 1]) {
print("the example is false");
}
}
}
for (int i = 1; i < vertexNum - 1; i++) {
if (b[i] - dist[0][i] < min1) {
min1 = b[i] - dist[0][i];
}
if (a[i] + s[i] + dist[i][vertexNum - 1] < min2) {
min2 = a[i] + s[i] + dist[i][vertexNum - 1];
}
}
if (E > min1 || L < min2) {
System.out.println("Duration false!");
System.exit(0);//终止程序
}
//初始化配送中心0,n+1两点的参数
arc[vertexNum - 1][0] = 0;
arc[0][vertexNum - 1] = 1;
for (int i = 1; i < vertexNum - 1; i++) {
arc[vertexNum - 1][i] = 0;//配送中心-顾客
}
for (int i = 1; i < vertexNum - 1; i++) {
arc[i][0] = 0;//顾客-配送中心
}
}
public static void display() {
printf("the number of customers:%10s\n", vertexNum);
printf("depot time windows:%12s[%s, %s]\n", " ", E, L);
printf("the number of vehicles:%10s\n", vehNum);
printf("vehicle capacity:%14s%.1f\n", " ", capacity);
for (int i = 0; i < xy.length; i++)
printf("vertex %03d:\t(\t%s,\t%s\t)\t%s\t(\t%s,\t%s\t)\t%s\n",
i, xy[i][0], xy[i][1], demand[i], a[i], b[i], s[i]);
}
}
package vrptw;
import java.util.ArrayList;
import ilog.concert.*;
import ilog.cplex.*;
import static net.mindview.util.Print.*;
//类功能:建立模型并求解
public class VRPTW {
IloCplex model; //定义cplex内部类的对象
public IloNumVar[][][] x; //x[i][j][k]表示弧arcs[i][j]被车辆k访问
public IloNumVar[][] w; //车辆访问所有点的时间矩阵
double cost; //目标值object
Solution solution;
//函数功能:根据VRPTW数学模型建立VRPTW的cplex模型
public void buildModel() throws IloException {
int vertexNum = Data.vertexNum;
int vehNum = Data.vehNum;
//model
model = new IloCplex();
model.setOut(null);
//variables
x = new IloNumVar[vertexNum][vertexNum][vehNum];
w = new IloNumVar[vertexNum][vehNum]; //车辆访问点的时间
//定义cplex变量x和w的数据类型及取值范围
for (int i = 0; i < vertexNum; i++) {
for (int k = 0; k < vehNum; k++) {
w[i][k] = model.numVar(0, 1e15, IloNumVarType.Float, "w" + i + "," + k);
}
for (int j = 0; j < vertexNum; j++) {
if (Data.arc[i][j] == 0) {
x[i][j] = null;
} else {
//Xijk,公式(10)-(11)
for (int k = 0; k < vehNum; k++) {
x[i][j][k] = model.numVar(0, 1, IloNumVarType.Int, "x" + i + "," + j + "," + k);
}
}
}
}
eq1(model, vertexNum, vehNum);
eq2(model, vertexNum, vehNum);
eq3(model, vertexNum, vehNum);
eq4(model, vertexNum, vehNum);
eq5(model, vertexNum, vehNum);
eq6(model, vertexNum, vehNum);
eq7(model, vertexNum, vehNum);
eq8(model, vertexNum, vehNum);
eq9(model, vertexNum, vehNum);
}
//函数功能:解模型,并生成车辆路径和得到目标值
public void solve() throws IloException {
ArrayList<ArrayList<Integer>> routes = new ArrayList<>(); //定义车辆路径链表
ArrayList<ArrayList<Double>> servetimes = new ArrayList<>(); //定义花费时间链表
//初始化车辆路径和花费时间链表,链表长度为车辆数k
for (int k = 0; k < Data.vehNum; k++) {
ArrayList<Integer> r = new ArrayList<>(); //定义一个对象为int型的链表
ArrayList<Double> t = new ArrayList<>(); //定义一个对象为double型的链表
routes.add(r); //将上述定义的链表加入到链表routes中
servetimes.add(t); //同上
}
//判断建立的模型是否可解
if (!model.solve()) {
//模型不可解
System.out.println("problem should not solve false!!!");
return; //若不可解,则直接跳出solve函数
} else {
//模型可解,生成车辆路径
for (int k = 0; k < Data.vehNum; k++) {
boolean terminate = true;
int i = 0;
routes.get(k).add(0);
servetimes.get(k).add(0.0);
while (terminate) {
for (int j = 0; j < Data.vertexNum; j++) {
if (Data.arc[i][j] >= 0.5 && model.getValue(x[i][j][k]) >= 0.5) {
routes.get(k).add(j);
servetimes.get(k).add(model.getValue(w[j][k]));
i = j;
break;
}
}
if (i == Data.vertexNum - 1) {
terminate = false;
}
}
}
}
solution = new Solution(routes, servetimes);
cost = model.getObjValue();
print("final solution:");
for (int i = 0; i < routes.size(); i++) {
print(solution.routes.get(i));
}
}
Objective
//Eq(1).
void eq1(IloCplex model, int vertexNum, int vehNum) throws IloException {
IloNumExpr obj = model.numExpr();
for (int i = 0; i < vertexNum; i++) {
for (int j = 0; j < vertexNum; j++) {
if (Data.arc[i][j] == 0)
continue;
for (int k = 0; k < vehNum; k++)
obj = model.sum(obj, model.prod(Data.dist[i][j], x[i][j][k]));
}
}
model.addMinimize(obj);
}
Constraint
//Eq(2).每个客户只能有一辆车服务一次
void eq2(IloCplex model, int vertexNum, int vehNum) throws IloException {
for (int i = 1; i < vertexNum - 1; i++) {
IloNumExpr cstr1 = model.numExpr();
for (int k = 0; k < vehNum; k++) {
for (int j = 1; j < vertexNum; j++) {
if (Data.arc[i][j] == 1)
cstr1 = model.sum(cstr1, x[i][j][k]);
}
}
model.addEq(cstr1, 1);
}
}
//Eq(3).车辆必须从配送中心0点出发
void eq3(IloCplex model, int vertexNum, int vehNum) throws IloException {
for (int k = 0; k < vehNum; k++) {
IloNumExpr expr2 = model.numExpr();
for (int j = 1; j < vertexNum; j++) {
if (Data.arc[0][j] == 1)
expr2 = model.sum(expr2, x[0][j][k]);
}
model.addEq(expr2, 1);
}
}
//Eq(4).第k辆车服务j点之后必须离开
void eq4(IloCplex model, int vertexNum, int vehNum) throws IloException {
for (int k = 0; k < vehNum; k++) {
for (int j = 1; j < vertexNum - 1; j++) {
IloNumExpr expr3;
IloNumExpr subExpr1 = model.numExpr();
IloNumExpr subExpr2 = model.numExpr();
for (int i = 0; i < vertexNum; i++) {
if (Data.arc[i][j] == 1) {
subExpr1 = model.sum(subExpr1, x[i][j][k]);
}
if (Data.arc[j][i] == 1) {
subExpr2 = model.sum(subExpr2, x[j][i][k]);
}
}
expr3 = model.sum(subExpr1, model.prod(-1, subExpr2));
model.addEq(expr3, 0);
}
}
}
//Eq(5).每辆车必须停留在配送中心n+1
void eq5(IloCplex model, int vertexNum, int vehNum) throws IloException {
for (int k = 0; k < vehNum; k++) {
IloNumExpr expr4 = model.numExpr();
for (int i = 0; i < vertexNum - 1; i++) {
if (Data.arc[i][vertexNum - 1] == 1) {
expr4 = model.sum(expr4, x[i][vertexNum - 1][k]);
}
}
model.addEq(expr4, 1);
}
}
void eq6(IloCplex model, int vertexNum, int vehNum) throws IloException {
//Eq(6)车辆时间约束
double M = 1e5;
for (int k = 0; k < vehNum; k++) {
for (int i = 0; i < vertexNum; i++) {
for (int j = 0; j < vertexNum; j++) {
if (Data.arc[i][j] == 1) {
IloNumExpr expr5;
IloNumExpr expr6;
expr5 = model.sum(w[i][k], Data.s[i] + Data.dist[i][j]);
expr5 = model.sum(expr5, model.prod(-1, w[j][k]));
expr6 = model.prod(M, model.sum(1, model.prod(-1, x[i][j][k])));
model.addLe(expr5, expr6);
}
}
}
}
}
void eq7(IloCplex model, int vertexNum, int vehNum) throws IloException {
//Eq(7)车辆时间约束
for (int k = 0; k < vehNum; k++) {
for (int i = 1; i < vertexNum - 1; i++) {
IloNumExpr expr7 = model.numExpr();
for (int j = 0; j < vertexNum; j++) {
if (Data.arc[i][j] == 1) {
expr7 = model.sum(expr7, x[i][j][k]);
}
}
model.addLe(model.prod(Data.a[i], expr7), w[i][k]);
model.addLe(w[i][k], model.prod(Data.b[i], expr7));
}
}
}
//Eq(8)车辆时间约束
void eq8(IloCplex model, int vertexNum, int vehNum) throws IloException {
for (int k = 0; k < vehNum; k++) {
model.addLe(Data.E, w[0][k]);
model.addLe(w[0][k], Data.L);
model.addLe(Data.E, w[vertexNum - 1][k]);
model.addLe(w[vertexNum - 1][k], Data.L);
}
}
//Eq(9)车辆容量约束
void eq9(IloCplex model, int vertexNum, int vehNum) throws IloException {
for (int k = 0; k < vehNum; k++) {
IloNumExpr expr8 = model.numExpr();
for (int i = 1; i < vertexNum - 1; i++) {
IloNumExpr expr9 = model.numExpr();
for (int j = 0; j < vertexNum; j++) {
if (Data.arc[i][j] == 1) {
expr9 = model.sum(expr9, x[i][j][k]);
}
}
expr8 = model.sum(expr8, model.prod(Data.demand[i], expr9));
}
model.addLe(expr8, Data.capacity);
}
}
}
package vrptw;
import java.util.ArrayList;
import static net.mindview.util.Print.*;
//类功能:解的可行性判断
class Solution {
ArrayList<ArrayList<Integer>> routes;
ArrayList<ArrayList<Double>> servetimes;
public Solution(ArrayList<ArrayList<Integer>> routes, ArrayList<ArrayList<Double>> servetimes) {
super();
this.routes = routes;
this.servetimes = servetimes;
}
//函数功能:解的可行性判断
public void fesible() {
int vehNum = Data.vehNum;
double capacity = Data.capacity; //车辆载荷
int[] demand = Data.demand; //需求量
double[] a = Data.a; //时间窗开始时间【a[i],b[i]】
double[] b = Data.b; //时间窗结束时间【a[i],b[i]】
double[][] dist = Data.dist; //距离矩阵,满足三角关系,暂用距离表示花费 C[i][j]=dist[i][j]
//车辆数量可行性判断
if (routes.size() > vehNum) {
System.err.println("error: vehicle number!!!");
System.exit(0);
}
//车辆载荷可行性判断
for (int k = 0; k < routes.size(); k++) {
ArrayList<Integer> route = routes.get(k);
double capasity = 0;
for (int i = 0; i < route.size(); i++) {
capasity += demand[route.get(i)];
}
if (capasity > capacity) {
System.err.println("error: vehicle capacity!!!");
System.exit(0);
}
}
//时间窗、车容量可行性判断
for (int k = 0; k < routes.size(); k++) {
ArrayList<Integer> route = routes.get(k);
ArrayList<Double> servertime = servetimes.get(k);
for (int i = 0; i < route.size() - 1; i++) {
int origin = route.get(i);
int destination = route.get(i + 1);
double si = servertime.get(i);
double sj = servertime.get(i + 1);
if (si < a[origin] && si > b[origin]) {
System.out.println("error: servertime!");
System.exit(0);
}
if (UtilMeth.doubleCompare(si + dist[origin][destination], b[destination]) > 0) {
print(origin + ": [" + a[origin] + "," + b[origin] + "]" + " " + si);
print(destination + ": [" + a[destination] + "," + b[destination] + "]" + " " + sj);
print(dist[origin][destination]);
print(destination + ":");
print("error: forward servertime!");
System.exit(0);
}
if (UtilMeth.doubleCompare(sj - dist[origin][destination], a[origin]) < 0) {
print(origin + ": [" + a[origin] + "," + b[origin] + "]" + " " + si);
print(destination + ": [" + a[destination] + "," + b[destination] + "]" + " " + sj);
print(dist[origin][destination]);
print(destination + ":");
print("error: backward servertime!");
System.exit(0);
}
}
}
}
}
package vrptw;
import static net.mindview.util.Print.printf;
public class Test {
public static void main(String[] args) throws Exception {
int vetexnum = 102;//所有点个数,包括0,n+1两个配送中心点
//读入不同的文件前要手动修改vetexnum参数,参数值等于所有点个数,包括配送中心
String path = "D:\\Users\\36297\\JavaWorkspace\\algori\\src\\vrptw\\c101.txt";//算例地址
Data.readFile(path, vetexnum);
Data.initialize();
Data.display();
VRPTW vrptw = new VRPTW();
vrptw.buildModel();
long t1 = System.currentTimeMillis();
vrptw.solve();
vrptw.solution.fesible();
long t2 = System.currentTimeMillis();
printf("求解时间:\t%ss\n目标函数值:\t%s", (t2 - t1) / 1000, vrptw.cost);
}
}
参考
- https://mp.weixin.qq.com/s/OdzX_b9sLJ2PZ4WoDM5NnQ