参考资料:微信公众号“数据魔法师”;Column Generation
模型如下:
我 在华中科技大学黄楠博士的原作 (非常感谢原作者!)的基础上,做了小修,代码如下:
整体思路:
先设置问题的背景,数据的存放和组织,具体如:客户节点、车辆数、每个车辆的路径、访问客户的时间段和时刻,这些都在类Data中。实际读入数据之后(见 函数process_solomon),删除掉没用的边(有的边绕远路,有的端点组合的时间窗、距离明显不可搭配等等),这里所说“删除”是指data.arc[i][j]设为0,表示不可走。
然后,设置模型,具体见函数build_model;
调用cplex解模型。
最后,输出解(要考虑把每个车辆要访问的客户按顺序记录到ArrayList中等等),见class Solution.
// VRPTW.java
package HuangVRPTW;
import ilog.concert.*;
import ilog.cplex.*;
import java.io.BufferedReader;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.util.ArrayList;
import java.util.Scanner;
//import HuangVRPTW.Data.Solution;
public class VRPTW {
static double epsilon = 0.0001;
Data data;
IloCplex model;
public IloNumVar[][][] x;
public IloNumVar[][] w;
double cost; // obj value
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<>();
// initialize 'routes', 'servetimes'
for (int k = 0; k < data.vecnum; k++) {
ArrayList<Integer> r = new ArrayList<>();
ArrayList<Double> t = new ArrayList<>();
routes.add(r); // 一辆车有一条路径,有一条服务时间线
servetimes.add(t);
}
if(model.solve()==false){
System.out.println("The Problem cannot be solved!");
return;
}else{
// we found a solution, so print it!
for (int k = 0; k < data.vecnum; k++) {
boolean terminate = true;
int i=0;
routes.get(k).add(0); // depot 0
servetimes.get(k).add(0.0);
while(terminate){
for (int j = 0; j < data.vetexnum; j++) {
if(data.arcs[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; // this car continues to move foward the next customer
break;
}
}
if(i==data.vetexnum-1){
terminate=false;
}
}
}
}
solution = new Solution(data,routes,servetimes);
cost = model.getObjValue();
System.out.println("routes= "+solution.routes);
} // function solve ends
private void build_model() throws IloException {
// 流程:
// 首先 建一个model;
// 然后,向model里面添加决策变量;
// 接着,yongmodel建立一个表达式作为目标函数,设置目标;
// 最后,求解这个model
// model
model = new IloCplex();
model.setOut(null); // ?? what is set out
// decision vars
x = new IloNumVar[data.vetexnum][data.vetexnum][data.vecnum];
w = new IloNumVar[data.vetexnum][data.vecnum];
// more details about decision vars
for (int i = 0; i < data.vetexnum; i++) {
for (int k = 0; k < data.vecnum; k++) {
w[i][k] = model.numVar(0,1e5,IloNumVarType.Float,"w"+i+","+k);
}
for (int j = 0; j < data.vetexnum; j++) {
if(data.arcs[i][j]==0){
x[i][j]=null;
}else{
for (int k = 0; k < data.vecnum; k++) {
x[i][j][k] = model.numVar(0,1,IloNumVarType.Int,"x"+i+","+j+","+k);
}
}
}
}// end for i
// obj func
IloNumExpr obj = model.numExpr();
for (int i = 0; i < data.vetexnum; i++) {
for (int j = 0; j < data.vetexnum; j++) {
if(data.arcs[i][j]==1){
for (int k = 0; k < data.vecnum; k++) {
obj = model.sum(obj,model.prod(data.dist[i][j],x[i][j][k]));
}
}
}
} // end for i
model.addMinimize(obj);
// constraints
// constraint 1, every customer needs to be served
for (int i = 1; i < data.vetexnum-1; i++) {
IloNumExpr expr1 = model.numExpr();
for (int j = 1; j < data.vetexnum; j++) {
if(data.arcs[i][j]==1){
for (int k = 0; k < data.vecnum; k++) {
expr1 = model.sum(expr1,x[i][j][k]);
}
}
}
model.addEq(expr1,1);
}
// constraint 2, every cars start from depot 0
for (int k = 0; k < data.vecnum; k++) {
IloNumExpr expr2 = model.numExpr();
for (int j = 1; j < data.vetexnum; j++) {
expr2 = model.sum(expr2,x[0][j][k]);
}
model.addEq(expr2,1);
}
// constraint 3, for every car and every customer, #inFlow = #outFlow
for (int k = 0; k < data.vecnum; k++) {
for (int i = 1; i < data.vetexnum-1; i++) {
IloNumExpr subExpr1 = model.numExpr();
IloNumExpr subExpr2 = model.numExpr();
for (int j = 0; j < data.vetexnum; j++) {
if(data.arcs[i][j]==1){
subExpr1 = model.sum(subExpr1,x[i][j][k]);
}
if(data.arcs[j][i]==1){
subExpr2 = model.sum(subExpr2,x[j][i][k]);
}
}
model.addEq(subExpr1,subExpr2);
}
}
// constraint 4, every car goes back to the depot n-1
for (int k = 0; k < data.vecnum; k++) {
IloNumExpr expr4 = model.numExpr();
for(int i=0;i<data.vetexnum-1;i++){
if(data.arcs[i][data.vetexnum-1]==1){
expr4 = model.sum(expr4,x[i][data.vetexnum-1][k]);
}
}
model.addEq(expr4,1);
}
// constraint 5, time window constraint
double M=1e5;
for (int k = 0; k < data.vecnum; k++) {
for (int i = 0; i < data.vetexnum; i++) {
for (int j=0;j<data.vetexnum;j++){
if(data.arcs[i][j]==1){
IloNumExpr expr5 = model.numExpr();
IloNumExpr expr6 = model.numExpr();
// 注意 data.s[i]+data.dist[i][j] 两项都是 数值,所以可以直接用加法
expr5 = model.sum(w[i][k],data.s[i]+data.dist[i][j]); // s[i], the time period of serving customer i
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);
}
}
}
}
// constraint 6, every customer is served during his own time window
for (int i = 1; i < data.vetexnum-1; i++) {
for (int k = 0; k < data.vecnum; k++) {
IloNumExpr expr5 = model.numExpr();
expr5 = model.sum(expr5,w[i][k]);
model.addLe(data.a[i],expr5);
model.addLe(expr5,data.b[i]);
}
}
// constraint 7, car's capacity
for (int k = 0; k < data.vecnum; k++) {
IloNumExpr expr7 = model.numExpr();
for (int i=1;i<data.vetexnum-1;i++){
IloNumExpr expr8= model.numExpr();
for (int j = 0; j < data.vetexnum; j++) {
if(data.arcs[i][j]==1){
expr8 = model.sum(expr8,x[i][j][k]);
}
}
expr7 = model.sum(expr7,model.prod(data.demands[i],expr8));
}
model.addLe(expr7,data.cap);
}
}// end build_model
public static void process_solomon(String path, Data data, int vetexnum) throws FileNotFoundException {
// 把 path指定的文件 读入到 data中
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+");
// initialize params
data.vetexnum = vetexnum;
data.vecnum = Integer.parseInt(substr[1]);
data.cap = Integer.parseInt(substr[2]); //
data.vertexs = new int[data.vetexnum][2];
data.demands = new int[data.vetexnum];
data.vehicles = new int[data.vecnum];
data.a = new double[data.vetexnum];
data.b = new double[data.vetexnum];
data.s = new double[data.vetexnum]; // time period of customers
data.arcs = new int[data.vetexnum][data.vetexnum];
data.dist = new double[data.vetexnum][data.vetexnum];
for (int i = 0; i < 4; i++) {
line = cin.nextLine();
}
// read the content of (data.vetexnum-1) lines
for (int i = 0; i < data.vetexnum-1; i++) {
line = cin.nextLine();
line.trim();
substr = line.split("\\s+");
data.vertexs[i][0] = Integer.parseInt(substr[2]);
data.vertexs[i][1]= Integer.parseInt(substr[3]);
data.demands[i] = Integer.parseInt(substr[4]);
data.a[i] = Integer.parseInt(substr[5]);
data.b[i] = Integer.parseInt(substr[6]);
data.s[i] = Integer.parseInt(substr[7]);
}
cin.close();
// depot n-1
data.vertexs[vetexnum-1] = data.vertexs[0];
data.demands[vetexnum-1] = 0;
data.a[vetexnum-1] = data.a[0];
data.b[vetexnum-1] = data.b[0];
data.s[vetexnum-1] = 0; // dont need time to be served
data.E=data.a[0];
data.L = data.b[vetexnum-1];
double min1 = 1e15;
double min2 = 1e15;
// initialize the distance matrix
for (int i = 0; i < data.vetexnum; i++) {
for (int j = 0; j < data.vetexnum; j++) {
if(i==j){
data.dist[i][j] = 0;
continue;
}
data.dist[i][j] = Math.sqrt(
(data.vertexs[i][0]-data.vertexs[j][0])*(data.vertexs[i][0]-data.vertexs[j][0])
+ (data.vertexs[i][1]-data.vertexs[j][1])*(data.vertexs[i][1]-data.vertexs[j][1]));
data.dist[i][j] = data.double_truncate(data.dist[i][j]);
}
}
data.dist[0][data.vetexnum-1]=0;
data.dist[data.vetexnum-1][0]=0;
// triangular relationship ?? why do this // I guess it's because you don't need to take extra effort。 你没必要走弯路,如果可以抄近道的话
for (int k = 0; k < data.vetexnum; k++) {
for (int i = 0; i < data.vetexnum; i++) {
for (int j = 0; j < data.vetexnum; j++) {
if(data.dist[i][j]>data.dist[i][k]+data.dist[k][j]){ // 我后来想到这里是不是也应该用double_compare function
data.dist[i][j]=data.dist[i][k]+data.dist[k][j];
}
}
}
}
// set data.arcs
// step 1,
for (int i = 0; i < data.vetexnum; i++) {
for (int j = 0; j < data.vetexnum; j++) {
if(i!=j){
data.arcs[i][j] = 1;
}else {
data.arcs[i][j] = 0;
}
}
}
// step 2, considering time window and cap
for (int i = 0; i < data.vetexnum; i++) {
for (int j = 0; j < data.vetexnum; j++) {
if(data.arcs[i][j]==1){
// 这里是double型时刻数据在比较大小,我觉得应该同这套代码的其他地方一样,调用double_compare函数
// 这里也是 和原作不同的地方
if(double_compare( // make double_compare static
data.a[i]+data.s[i]+data.dist[i][j],data.b[j])>0||
data.demands[i]+data.demands[j]> data.cap){
data.arcs[i][j] = 0;
}
if(double_compare( // make double_compare static
data.a[0]+data.s[i]+data.dist[0][i]+data.dist[i][data.vetexnum-1],data.b[data.vetexnum-1])>0){
System.out.println("the current data is false");
}
}
}
}
// step 3: considering time window
for (int i = 1; i < data.vetexnum-1; i++) {
if(double_compare( data.b[i]-data.dist[0][i],min1)<0){
min1 = data.b[i]-data.dist[0][i];
}
if(double_compare(data.a[i]+data.s[i]+data.dist[i][data.vetexnum-1],min2)<0){
min2 = data.a[i]+data.s[i]+data.dist[i][data.vetexnum-1];
}
}
if(min1<data.E || min2>data.L){
System.out.println("Duration false");
}
// step 4 : considering depot 0 and depot n-1 as the start and the terminal
data.arcs[0][data.vetexnum-1]=1;
data.arcs[data.vetexnum-1][0]=0;
for (int i = 1; i < data.vetexnum-1; i++) { // 客户点 无法回到 分配点0
data.arcs[i][0] = 0;
}
for (int i = 1; i < data.vetexnum-1; i++) { // 分配点n-1 无法去到客户点
data.arcs[data.vetexnum-1][i]=0;
}
}
public static int double_compare(double a, double b){
if(a<b-epsilon){
return -1;
}
if(a>b+epsilon){
return 1;
}
return 0;
}
public static void main(String[] args) throws FileNotFoundException, IloException {
Data data = new Data();
int vetexnum=102; // 100+2 // 100 customers and 2 depots
String path = "D:\\VRP问题实践\\CGVRPTW_Model-master\\CGVRPTW_Seminar\\resource\\c101-2.txt";//算例地址
process_solomon(path,data,vetexnum);
System.out.println("Input Successfully");
System.out.println("Cplex procedure ....");
VRPTW model = new VRPTW(data);
model.build_model();
double time1 = System.nanoTime();
model.solve();
model.solution.feasible();
double time2 = System.nanoTime();
System.out.println("Runing time is " + (time2 - time1) / 1e9 + ", best cost is " + model.cost);
}
}
下面的文件同上面的文件在同一个目录下。
// Data_Solution.java
package HuangVRPTW;
import ilog.concert.*;
import ilog.cplex.*;
import java.util.ArrayList;
// this file defines two class, class Data and class Solution
// class Solution has a Data-class object
// class Data is used to read data from a given file
// class Solution is used to print the final solution
class Data {
int vetexnum;// #customers + 2 // 2 depots
double E; // start of the whole duration
double L;// end of the whole duration
int vecnum; // #vehicles
double cap;// vehicle capacity
int[][] vertexs;
int[] demands;
int[] vehicles; // carNo. 车辆编号
double[] a; // ready time of Customers
double[] b; // due time of Customers
double[] s; // service time period of customers
int[][] arcs; // arcs[i][j]=1 means this arc is feasible
double[][] dist;
public double double_truncate(double v){
int iv = (int)v;
if(iv+1 -v<=0.000000001){
return iv+1;
}
double dv = (v-iv)*10;
int idv = (int)dv;
double rv = iv + idv/10.0; // !!/10.0 , not /10
return rv;
}
}
class Solution{ // inner class of class Data
double epsilon = 0.0001;
Data data = new Data();
ArrayList<ArrayList<Integer>> routes = new ArrayList<>();// routes of all the vehicles
ArrayList<ArrayList<Double>> servetimes = new ArrayList<>();//the time point of every customer on the route
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 a, double b){ // Double 型数据比较大小要注意!
if(a<b-epsilon){
return -1;
}
if(a>b+epsilon){
return 1;
}
return 0;
}
public void feasible() {
if(routes.size()>data.vecnum){ // feasibility of the amount of cars
System.out.println("Error: vecnum!!");
System.exit(0);
}
// feasibility of capacity of cars
for (int k = 0; k < data.vecnum; k++) {
ArrayList<Integer> route = routes.get(k);
double trueCap = 0;
for(int i=0;i<route.size();i++){
trueCap += data.demands[route.get(i)];
}
if(trueCap>data.cap){
System.out.println("error: cap!!");
System.exit(0);
}
}
// feasibility of time window
for (int k = 0; k < data.vecnum; k++) {
ArrayList<Integer> route = new ArrayList<>();
for (int i = 0; i < route.size()-1; i++) {
int origin = route.get(i);
int destination = route.get(i+1);
double si = servetimes.get(k).get(i);
double sj = servetimes.get(k).get(i+1);
if(si<data.a[origin] || si>data.b[origin]){
System.out.println("Error: serveTime");
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 servetime");
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 serveTime");
System.exit(0);
}
}
}
}
}
运行结果:
Input Successfully
Cplex procedure ....
routes= [[0, 81, 78, 76, 71, 70, 73, 77, 79, 80, 101],
[0, 57, 55, 54, 53, 56, 58, 60, 59, 101],
[0, 98, 96, 95, 94, 92, 93, 97, 100, 99, 101],
[0, 90, 87, 86, 83, 82, 84, 85, 88, 89, 91, 101],
[0, 5, 3, 7, 8, 10, 11, 9, 6, 4, 2, 1, 75, 101],
[0, 101],
[0, 20, 24, 25, 27, 29, 30, 28, 26, 23, 22, 21, 101], [0, 32, 33, 31, 35, 37, 38, 39, 36, 34, 101], [0, 67, 65, 63, 62, 74, 72, 61, 64, 68, 66, 69, 101], [0, 13, 17, 18, 19, 15, 16, 14, 12, 101], [0, 101], [0, 43, 42, 41, 40, 44, 46, 45, 48, 51, 50, 52, 49, 47, 101]]
Runing time is 16.1905963, best cost is 827.3
算例文件中 限制车辆数是12,输出答案中所有车辆都有自己的路径,只是有的车不被使用,比如:
这里[0, 101]
表示对应的这辆车不被使用(因为depot 101 and depot 0 are the same one.)