这个内容是转载自华中科技大学——数据魔法师团队的原创代码(Java调用cplex求解VRPTW),增加了一些注解和心得,相信站在巨人的肩膀上可以看得更远学得更快,谢谢原作者们辛苦的汗水。
首先,列出VRPTW模型,如下:
其次,列出算法源码,分三个板块,具体如下:
具体Java怎么调用cplex可以下载我的方法文档(链接:https://download.csdn.net/download/weixin_39106536/11290822),下面就是具体的代码了,我在原有代码中增加了自己的理解
import java.io.BufferedReader;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.util.Scanner;
public class Data {
public static int pointnum=102;//所有点集合n(包括配送中心和客户点,首尾(0和n)为配送中心)
public static int carnum;//车辆数
public static int cap;//车辆载荷
public static double E;//配送中心时间窗开始时间
public static double L;//配送中心时间窗结束时间
public static int point[][]=new int[pointnum][2];
public static int s[]=new int[pointnum];//客户点的服务时间
public static double a[]=new double[pointnum];//时间窗结束时间【a[i],b[i]】
public static double b[]=new double[pointnum];//时间窗结束时间【a[i],b[i]】
public static int demand[]=new int[pointnum];//需求量
public static int car[]=new int[pointnum];//车辆编号
public static double dist[][]=new double[pointnum][pointnum];//距离矩阵,满足三角关系,暂用距离表示花费 C[i][j]=dist[i][j]
public static double arcs[][]=new double[pointnum][pointnum];
//截断小数3.26434-->3.2
public static double double_truncate(double v){//取前几位数字
int iv = (int) v;
if(iv+1 - v <= 0.000000000001)
return iv+1;
double dv = (v - iv) * 10;
int idv = (int) dv;
double rv = iv + idv / 10.0;
return rv;
}
//读取数据并且初始化参数
public static void read_data(String path,Data data) throws Exception {
String line=null;
String substr[]=null;
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+");
carnum=Integer.parseInt(substr[1]);
cap=Integer.parseInt(substr[2]);
for (int i = 0; i < 4; i++) {
line=cin.nextLine();
}
for (int i = 0; i < pointnum-1; i++) {
line=cin.nextLine();
line.trim();//去除空格
substr=line.split("\\s+");//以空格分开
point[i][0]=Integer.parseInt(substr[2]);//初始化定义为整数
point[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();
//终点数据初始化
point[pointnum-1]=point[0];
demand[pointnum-1]=demand[0];
a[pointnum-1]=a[0];
b[pointnum-1]=b[0];
s[pointnum-1]=s[0];
E=a[0];
L=b[0];
double min1=1e15;
double min2=1e15;
//距离矩阵初始化
for (int i = 0; i < pointnum; i++) {
for (int j = 0; j < pointnum; j++) {
if(i==j) {
dist[i][j]=0;
continue;
}else {
dist[i][j]=Math.sqrt(Math.pow(point[i][0]-point[j][0], 2)+Math.pow(point[i][1]-point[j][1], 2));
dist[i][j]=double_truncate(dist[i][j]);
}
}
}
dist[pointnum-1][0]=0;
dist[0][pointnum-1]=0;//起点和终点同一位置
//距离矩阵满足三角关系,满足最小行走规则,不绕路
for (int k = 0; k < pointnum; k++) {
for (int i = 0; i < pointnum; i++) {
for (int j = 0; j < pointnum; j++) {
if(dist[i][j]>dist[i][k]+dist[k][j]) {
dist[i][j]=dist[i][k]+dist[k][j];
}
}
}
}
//初始化为完全图,先不管约束,先把所有的弧连起来,除了自身不能原地打转,后续再消除
for (int i = 0; i < pointnum; i++) {
for (int j = 0; j < pointnum; j++) {
if(i==j) {
arcs[i][j]=0;//不能原地打转
}else {
arcs[i][j]=1;
}
}
}
//除去不符合时间窗和容量约束的边
for (int i = 0; i < pointnum; i++) {
for (int j = 0; j < pointnum; j++) {
if(i==j) {
continue;
}
if(a[i]+s[i]+dist[i][j]>b[j]||demand[i]+demand[j]>cap) {//初始化必须满足容量和服务时间,以消除多余的边
arcs[i][j]=0;
}
if(a[0]+s[i]+dist[0][i]+dist[i][pointnum-1]>b[pointnum-1]) {
System.out.println("计算错误");
}
}
}
for (int i = 1; i < pointnum-1; i++) {
if(b[i]-dist[0][i]<min1) {
min1=b[i]-dist[0][i];
}
if(a[i]+s[i]+dist[i][pointnum-1]<min2) {
min2=a[i]+s[i]+dist[i][pointnum-1];
}
}
if(E>min1||L<min2) {
System.out.println("计算错误");
System.exit(0);
}
//初始化配送中心0,n+1的arcs
arcs[pointnum-1][0]=0;
arcs[0][pointnum-1]=1;
for (int i = 1; i < pointnum-1; i++) {
arcs[pointnum-1][i]=0;
arcs[i][0]=1;
}
}
}
import java.util.ArrayList;
import ilog.concert.IloException;
public class Solution {
double epsilon = 0.0001;
Data data = new Data();
ArrayList<ArrayList<Integer>> routes = new ArrayList<>();
ArrayList<ArrayList<Double>> servetimes = new ArrayList<>();
public Solution(Data data, ArrayList<ArrayList<Integer>> routes, ArrayList<ArrayList<Double>> servetimes) {
super();
this.data = data;
this.routes = routes;
this.servetimes = servetimes;
}
//函数功能:比较两个数的大小
public int double_compare(double v1,double v2) {
if (v1 < v2 - epsilon) {
return -1;
}
if (v1 > v2 + epsilon) {//[-epsilon,+epsilon]
return 1;
}
return 0;
}
//函数功能:解的可行性判断
public void fesible() throws IloException {
//车辆数量可行性判断
if (routes.size() > data.carnum) {
System.out.println("error: carnum!!!");
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 += data.demand[route.get(i)];
}
if (capasity > data.cap) {
System.out.println("error: cap!!!");//过程中装货容量检查,溢出报错
System.exit(0);
}
}
//时间窗、车容量可行性判断
for (int k = 0; k < routes.size(); k++) {
ArrayList<Integer> route = routes.get(k);
ArrayList<Double> servertime = servetimes.get(k);//定义一个序列列表
double capasity = 0;
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 < data.a[origin] && si > data.b[origin]) {
System.out.println("error: servertime!");
System.exit(0);
}
if (double_compare(si + data.dist[origin][destination],data.b[destination]) > 0) {
System.out.println(origin + ": [" + data.a[origin] + ","+data.b[origin]+"]"+ " "+ si);
System.out.println(destination + ": [" + data.a[destination] + ","+data.b[destination]+"]"+ " "+ sj);
System.out.println(data.dist[origin][destination]);
System.out.println(destination + ":" );
System.out.println("error: forward servertime!");//来早了,晚点来
System.exit(0);
}
if (double_compare(sj - data.dist[origin][destination],data.a[origin]) < 0) {
System.out.println(origin + ": [" + data.a[origin] + ","+data.b[origin]+"]"+ " "+ si);
System.out.println(destination + ": [" + data.a[destination] + ","+data.b[destination]+"]"+ " "+ sj);
System.out.println(data.dist[origin][destination]);
System.out.println(destination + ":" );
System.out.println("error: backward servertime!");//来晚了,早点来
System.exit(0);
}
}
if (capasity > data.cap) {
System.out.println("error: cap!!!");//容量溢出了,容量为0的有限排除
System.exit(0);
}
}
}
}
import java.util.ArrayList;
import ilog.concert.IloException;
import ilog.concert.IloNumExpr;
import ilog.concert.IloNumVar;
import ilog.concert.IloNumVarType;
import ilog.cplex.IloCplex;
public class Vrptw {
Data data; //定义类Data的对象
IloCplex cplex; //定义cplex内部类的对象
public IloNumVar[][][] x; //x[i][j][k]表示弧arcs[i][j]被车辆k访问
public IloNumVar[][] w; //车辆访问所有点的时间矩阵
double cost; //目标值object
Solution solution;
public Vrptw(Data data) {
this.data = data;
}
//函数功能:解模型,并生成车辆路径和得到目标值
public void solve() throws IloException {
ArrayList<ArrayList<Integer>> routes = new ArrayList<>(); //定义车辆路径链表
ArrayList<ArrayList<Double>> servetimes = new ArrayList<>(); //定义花费时间链表
//初始化车辆路径和花费时间链表,链表长度为车辆数k
for (int k = 0; k < data.carnum; k++) {
ArrayList<Integer> r = new ArrayList<>(); //定义一个对象为int型的链表
ArrayList<Double> t = new ArrayList<>(); //定义一个对象为double型的链表
routes.add(r); //将上述定义的链表加入到链表routes中
servetimes.add(t); //同上
}
//判断建立的模型是否可解
if(cplex.solve() == false){
//模型不可解
System.out.println("problem should not solve false!!!");
return; //若不可解,则直接跳出solve函数
}
else{
//模型可解,生成车辆路径
for(int k = 0; k < data.carnum; 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.pointnum; j++) {
if (data.arcs[i][j]>=0.5 && cplex.getValue(x[i][j][k])>=0.5) {
routes.get(k).add(j);
servetimes.get(k).add(cplex.getValue(w[j][k]));
i = j;
break;
}
}
if (i == data.pointnum-1) {
terminate = false;
}
}
}
}
solution = new Solution(data,routes,servetimes);
cost = cplex.getObjValue();
System.out.println("routes="+solution.routes);
}
//函数功能:根据VRPTW数学模型建立VRPTW的cplex模型
private void build_cplex() throws IloException {
//model
cplex = new IloCplex();
cplex.setOut(null);
//variables
x = new IloNumVar[data.pointnum][data.pointnum][data.carnum];
w = new IloNumVar[data.pointnum][data.carnum]; //车辆访问点的时间
//定义cplex变量x和w的数据类型及取值范围
for (int i = 0; i < data.pointnum; i++) {
for (int k = 0; k < data.carnum; k++) {
w[i][k] = cplex.numVar(0, 1e15, IloNumVarType.Float, "w" + i + "," + k);
}
for (int j = 0; j < data.pointnum; j++) {
if (data.arcs[i][j]==0) {
x[i][j] = null;
}
else{
//Xijk,公式(10)-(11)
for (int k = 0; k < data.carnum; k++) {
x[i][j][k] = cplex.numVar(0, 1, IloNumVarType.Int, "x" + i + "," + j + "," + k);
}
}
}
}
//加入目标函数
//公式(1)
IloNumExpr obj = cplex.numExpr();
for(int i = 0; i < data.pointnum; i++){
for(int j = 0; j < data.pointnum; j++){
if (data.arcs[i][j]==0) {
continue;//不计算不连通的弧,也即没连起来的节点不计算
}
for(int k = 0; k < data.carnum; k++){
obj = cplex.sum(obj, cplex.prod(data.dist[i][j], x[i][j][k]));
}
}
}
cplex.addMinimize(obj);//最小化目标函数
//加入约束
//公式(2)
for(int i= 1; i < data.pointnum-1;i++){
IloNumExpr expr1 = cplex.numExpr();
for (int k = 0; k < data.carnum; k++) {
for (int j = 1; j < data.pointnum; j++) {
if (data.arcs[i][j]==1) {
expr1 = cplex.sum(expr1, x[i][j][k]);//每个节点在所有车辆当中,求入度和
}
}
}
cplex.addEq(expr1, 1);//求入度和与1进行对比,也即每个节点最多有一辆车辆服务
}
//公式(3)
for (int k = 0; k < data.carnum; k++) {
IloNumExpr expr2 = cplex.numExpr();
for (int j = 1; j < data.pointnum; j++) {
if (data.arcs[0][j]==1) {
expr2 = cplex.sum(expr2, x[0][j][k]);//所有从配送中心出发的出度求和
}
}
cplex.addEq(expr2, 1);//求和与1进行对比,也即对每一辆车必须从配送中心出发
}
//公式(4)
for (int k = 0; k < data.carnum; k++) {
for (int j = 1; j < data.pointnum-1; j++) {
IloNumExpr expr3 = cplex.numExpr();
IloNumExpr subExpr1 = cplex.numExpr();
IloNumExpr subExpr2 = cplex.numExpr();
for (int i = 0; i < data.pointnum; i++) {
if (data.arcs[i][j]==1) {
subExpr1 = cplex.sum(subExpr1,x[i][j][k]);//入度+1
}
if (data.arcs[j][i]==1) {
subExpr2 = cplex.sum(subExpr2,x[j][i][k]);//如果i,j两节点之间连通,出度+1
}
}
expr3 = cplex.sum(subExpr1,cplex.prod(-1, subExpr2));//两个加和相减,表示车辆服务该节点必须离开
cplex.addEq(expr3, 0);//之差和0进行对比,也即对于每个节点都满足:出度=入度
}
}
//公式(5)
for (int k = 0; k < data.carnum; k++) {
IloNumExpr expr4 = cplex.numExpr();
for (int i = 0; i < data.pointnum-1; i++) {
if (data.arcs[i][data.pointnum-1]==1) {
expr4 = cplex.sum(expr4,x[i][data.pointnum-1][k]);//每辆车最后都必须回到配送中心,也即无论那一节点i都必须指向终点
}
}
cplex.addEq(expr4, 1);//与0进行比较
}
//公式(6)
double M = 1e5;//首先定义一个超级大的数,也可以定义其它无穷的数
for (int k = 0; k < data.carnum; k++) {
for (int i = 0; i < data.pointnum; i++) {
for (int j = 0; j < data.pointnum; j++) {
if (data.arcs[i][j] == 1) {
IloNumExpr expr5 = cplex.numExpr();
IloNumExpr expr6 = cplex.numExpr();
expr5 = cplex.sum(w[i][k], data.s[i]+data.dist[i][j]);//wik+si+tij
expr5 = cplex.sum(expr5,cplex.prod(-1, w[j][k]));//更新为wik+si+tij-wjk
expr6 = cplex.prod(M,cplex.sum(1,cplex.prod(-1, x[i][j][k])));//公式(6)的右项,M(1-xijk)
cplex.addLe(expr5, expr6);//满足左右项约束成立
}
}
}
}
//公式(7)
for (int k = 0; k < data.carnum; k++) {
for (int i = 1; i < data.pointnum-1; i++) {
IloNumExpr expr7 = cplex.numExpr();
for (int j = 0; j < data.pointnum; j++) {
if (data.arcs[i][j] == 1) {
expr7 = cplex.sum(expr7,x[i][j][k]);
}
}
cplex.addLe(cplex.prod(data.a[i], expr7), w[i][k]);
cplex.addLe(w[i][k], cplex.prod(data.b[i], expr7));
}
}
//公式(8)
for (int k = 0; k < data.carnum; k++) {
cplex.addLe(data.E, w[0][k]);//服务时间必须在服务窗内
cplex.addLe(data.E, w[data.pointnum-1][k]);
cplex.addLe(w[0][k], data.L);
cplex.addLe(w[data.pointnum-1][k], data.L);
}
//公式(9)
for (int k = 0; k < data.carnum; k++) {
IloNumExpr expr8 = cplex.numExpr();//定义式子
for (int i = 1; i < data.pointnum-1; i++) {
IloNumExpr expr9 = cplex.numExpr();
for (int j = 0; j < data.pointnum; j++) {
if (data.arcs[i][j] == 1) {
expr9=cplex.sum(expr9,x[i][j][k]);//所有节点和i节点连通求和,对每辆车,也即第i个节点路过
}
}
expr8 = cplex.sum(expr8,cplex.prod(data.demand[i],expr9));//所有路过的节点需求求和,对应需求相乘
}
cplex.addLe(expr8, data.cap);//小于最大容量
}
}
public static void main(String[] args) throws Exception {
Data data=new Data();//创建一个数据对象
String path="data/c102.txt";//数据的路径
data.read_data(path, data);//导入数据
System.out.println("输入成功");//看看有没有导入进来
System.out.println("程序执行中**************************");
Vrptw cplex=new Vrptw(data);//将数据转换成VRPTW的框架
cplex.build_cplex();//搭建计算对象
double cplex_time1=System.nanoTime();//起始时间
cplex.solve();//进行求解
cplex.solution.fesible();//输出可行解
double cplex_time2=System.nanoTime();//结束时间
double cplex_time=( cplex_time2- cplex_time1)/1e9;//求解时间
System.out.println("cplex_time:"+cplex_time);//输出时间
System.out.println("bestcost:"+cplex.cost);//输出最优目标值
}
}
最后的计算结果如下:
心得和理解:总体而言,上面的代码可以看作是求解一般运筹规划问题的框架了,可以在其基础上进一步修改以适应自己的研究问题,关键问题——根据自己的问题增加变量和添加新的约束,看懂了就很简单了。优点:资料比较全面,扩展性极强。缺点:有些地方注释不够,对新手不是很友好。再次感谢这个优秀团队的贡献。
参考来源:https://paste.ubuntu.com/25476905/和https://blog.csdn.net/cc098818/article/details/99624281?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromBaidu-4.not_use_machine_learn_pai&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromBaidu-4.not_use_machine_learn_pai。