【運籌學教學?勤實踐】-2 分支定界法解帶時間窗的車輛路徑規劃問題（附代碼及注釋）

05-22

來自專欄數據魔術師

歷盡千辛萬苦，外加外援幫助，本辣雞小編終於搞定了這個大坑-用分支定界法(Branch and bound, B&B)解帶時間窗的車輛路徑規劃問題(VRPTW)。

預備知識

前面的推文中有提到過，分支定界法是一種精確解演算法，之前推文「運籌學教學|分枝定界求解旅行商問題」中對於分支定界的基本思想進行了詳細的闡述，有不記得的小夥伴可以點擊上面的鏈接傳送到之前推文。

帶時間窗的車輛路徑規劃問題(下簡稱：VRPTW)在之前的推文中已經被詳細的介紹過了，為了方便讀者的閱讀，我們在這裡給出傳送門

乾貨|十分鐘快速掌握CPLEX求解VRPTW數學模型（附JAVA代碼及CPLEX安裝流程）

除此之外還要先學習一種數據結構叫做優先隊列。優先隊列(priority queue)是一種常用的數據結構，在這種數據結構中，隊頭永遠是存儲優先順序最高的元素，取隊頭和插入元素的操作的時間複雜度都是O(logn)。在JAVA和C++中都內置了這一種數據結構，因此，親愛的讀者們不要害怕。當然如果有代碼實現能力強的讀者想要手工實現優先隊列，也是可以的，想要學習優先隊列以事先學習堆(heap)這種數據結構，可以完美的實現優先隊列的功能。

當你仔細的閱讀了上面兩篇推文並理解了優先隊列的原理之後，小編相信聰明的你一定不會對於接下來要講的內容感到陌生。

代碼以及解釋

代碼共分為4個類包括：

BaB_Vrptw ：主類，用於建模以及分支定界求解VRPTW。
Check ：解的可行性判斷
Data ：定義參數
Node ：定義分支定界的節點

Data 類

Data類的作用就是讀入數據以及數據預處理，在這裡我們便不做過多的解釋，為了方便後面的閱讀以及篇幅限制，我們在這裡便不對其進行展開描述，代碼中的注釋對於各個變數含義有較為詳細的介紹。但是由於之後的程序會調用這些變數，我們便首先講解這個類。

class Data{

int vertex_num; //所有點集合n（包括配送中心和客戶點，首尾（0和n）為配送中心）

double E; //配送中心時間窗開始時間

double L; //配送中心時間窗結束時間

int veh_num; //車輛數

double cap; //車輛載荷

int[][] vertexs; //所有點的坐標x,y

int[] demands; //需求量

int[] vehicles; //車輛編號

double[] a; //時間窗開始時間【a[i],b[i]】

double[] b; //時間窗結束時間【a[i],b[i]】

double[] s; //客戶點的服務時間

int[][] arcs; //arcs[i][j]表示i到j點的弧

double[][] dist; //距離矩陣，滿足三角關係,暫用距離表示花費 C[i][j]=dist[i][j]

double gap= 1e-6; // 一個小數，表示精讀

double big_num = 100000; // 無窮大

//截斷小數3.26434-->3.2

public double double_truncate(double v){

int iv = (int) v;

if(iv+1 - v <= gap)

return iv+1;

double dv = (v - iv) * 10;

int idv = (int) dv;

double rv = iv + idv / 10.0;

return rv;

}

public Data() {

super();

}

//函數功能：從txt文件中讀取數據並初始化參數

public void Read_data(String path,Data data,int vertexnum) 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+")); //以空格為標誌將字元串拆分

//初始化參數

data.vertex_num = vertexnum;

data.veh_num = Integer.parseInt(substr[1]);

data.cap = Integer.parseInt(substr[2]);

data.vertexs =new int[data.vertex_num][2]; //所有點的坐標x,y

data.demands = new int[data.vertex_num]; //需求量

data.vehicles = new int[data.veh_num]; //車輛編號

data.a = new double[data.vertex_num]; //時間窗開始時間

data.b = new double[data.vertex_num]; //時間窗結束時間

data.s = new double[data.vertex_num]; //服務時間

data.arcs = new int[data.vertex_num][data.vertex_num];

//距離矩陣,滿足三角關係,用距離表示cost

data.dist = new double[data.vertex_num][data.vertex_num];

for(int i =0; i < 4;i++){

line = cin.nextLine();

}

//讀取vetexnum-1行數據

for (int i = 0; i < data.vertex_num - 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();//關閉流

//初始化配送中心參數

data.vertexs[data.vertex_num-1] = data.vertexs[0];

data.demands[data.vertex_num-1] = 0;

data.a[data.vertex_num-1] = data.a[0];

data.b[data.vertex_num-1] = data.b[0];

data.E = data.a[0];

data.L = data.b[0];

data.s[data.vertex_num-1] = 0;

double min1 = 1e15;

double min2 = 1e15;

//距離矩陣初始化

for (int i = 0; i < data.vertex_num; i++) {

for (int j = 0; j < data.vertex_num; 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.vertex_num-1] = 0;

data.dist[data.vertex_num-1][0] = 0;

//距離矩陣滿足三角關係

for (int k = 0; k < data.vertex_num; k++) {

for (int i = 0; i < data.vertex_num; i++) {

for (int j = 0; j < data.vertex_num; j++) {

if (data.dist[i][j] > data.dist[i][k] + data.dist[k][j]) {

data.dist[i][j] = data.dist[i][k] + data.dist[k][j];

}

//初始化為完全圖

for (int i = 0; i < data.vertex_num; i++) {

for (int j = 0; j < data.vertex_num; j++) {

if (i != j) {

data.arcs[i][j] = 1;

}

else {

data.arcs[i][j] = 0;

}

//除去不符合時間窗和容量約束的邊

for (int i = 0; i < data.vertex_num; i++) {

for (int j = 0; j < data.vertex_num; j++) {

if (i == j) {

continue;

}

if (data.a[i]+data.s[i]+data.dist[i][j]>data.b[j] ||

data.demands[i]+data.demands[j]>data.cap) {

data.arcs[i][j] = 0;

}

if (data.a[0]+data.s[i]+data.dist[0][i]+data.dist[i][data.vertex_num-1]>

data.b[data.vertex_num-1]) {

System.out.println("the calculating example is false");

}

for (int i = 1; i < data.vertex_num-1; i++) {

if (data.b[i] - data.dist[0][i] < min1) {

min1 = data.b[i] - data.dist[0][i];

}

if (data.a[i] + data.s[i] + data.dist[i][data.vertex_num-1] < min2) {

min2 = data.a[i] + data.s[i] + data.dist[i][data.vertex_num-1];

}

if (data.E > min1 || data.L < min2) {

System.out.println("Duration false!");

System.exit(0);//終止程序

}

//初始化配送中心0，n+1兩點的參數

data.arcs[data.vertex_num-1][0] = 0;

data.arcs[0][data.vertex_num-1] = 1;

for (int i = 1; i < data.vertex_num-1; i++) {

data.arcs[data.vertex_num-1][i] = 0;

}

for (int i = 1; i < data.vertex_num-1; i++) {

data.arcs[i][0] = 0;

}

Node類

Node類的主要作用是記錄分支節點，下面一段代碼是Node類定義的對象

Data data;

int d;

double node_cost; //目標值object

double[][][]lp_x;//記錄lp解

int[][][] node_x_map;//node_xij=1時,node_x_mapijk=1表示必須訪問，node_x_mapijk=0表示不能訪問

int[][] node_x;//0表示弧可以訪問，1表示必須訪問，-1表示不能訪問

ArrayList<ArrayList<Integer>> node_routes; //定義車輛路徑鏈表

ArrayList<ArrayList<Double>> node_servetimes; //定義花費時間鏈表

Node類的初始化如下，注意新生成的node_cost 的初始值是無窮大，因為在沒有操作的情況下，這是一個非法解。

public Node(Data data) {

super();

this.data = data;

node_cost = data.big_num;

lp_x = new double [data.vertex_num][data.vertex_num][data.veh_num];

node_x_map = new int[data.vertex_num][data.vertex_num][data.veh_num];

node_x = new int[data.vertex_num][data.vertex_num];

node_routes = new ArrayList<ArrayList<Integer>>();

node_servetimes = new ArrayList<ArrayList<Double>>();

}

由於要進行多次的生成節點，為了方便，我們設置了一個函數note_copy()來完成這項操作以及兩個節點比較大小的函數。

public Node note_copy() {

Node new_node = new Node(data);

new_node.d = d;

new_node.node_cost = node_cost;

for (int i = 0; i < lp_x.length; i++) {

for (int j = 0; j < lp_x[i].length; j++) {

new_node.lp_x[i][j] = lp_x[i][j].clone();

}

for (int i = 0; i < node_x.length; i++) {

new_node.node_x[i] = node_x[i].clone();

}

for (int i = 0; i < node_x_map.length; i++) {

for (int j = 0; j < node_x_map[i].length; j++) {

new_node.node_x_map[i][j] = node_x_map[i][j].clone();

}

for (int i = 0; i < node_routes.size(); i++) {

new_node.node_routes.add((ArrayList<Integer>) node_routes.get(i).clone());

}

for (int i = 0; i < node_servetimes.size(); i++) {

new_node.node_servetimes.add((ArrayList<Double>) node_servetimes.get(i).clone());

}

return new_node;

}

public int compareTo(Object o){

Node node = (Node) o;

if(node_cost < node.node_cost)

return -1;

else if(node_cost == node.node_cost)

return 0;

else

return 1;

}

BaB_Vrptw類

這是整個程序中最重要的一個部分，因此本文將花費大篇幅來講解這個問題(????)??嗨。看程序先看什麼？答案是-主函數。

public static void main(String[] args) throws Exception {

Data data = new Data();

int vetexnum = 102;//所有點個數，包括0，n+1兩個配送中心點

//讀入不同的文件前要手動修改vetexnum參數，參數值等於所有點個數,包括配送中心

String path = "data/c102.txt";//算例地址

data.Read_data(path,data,vetexnum);

System.out.println("input succesfully");

System.out.println("cplex procedure###########################");

BaB_Vrptw lp = new BaB_Vrptw(data);

double cplex_time1 = System.nanoTime();

//刪除未用的車輛，縮小解空間

lp=lp.init(lp);

System.out.println(": "+lp.data.veh_num);

lp.branch_and_bound(lp);

Check check = new Check(lp);

check.fesible();

double cplex_time2 = System.nanoTime();

double cplex_time = (cplex_time2 - cplex_time1) / 1e9;//求解時間，單位s

System.out.println("cplex_time " + cplex_time + " bestcost " + lp.cur_best);

for (int i = 0; i < lp.best_note.node_routes.size(); i++) {

ArrayList<Integer> r = lp.best_note.node_routes.get(i);

System.out.println();

for (int j = 0; j < r.size(); j++) {

System.out.print(r.get(j)+" ");

}

上面的函數就是主函數，前面的11行都是數據讀入的內容，相信大家都能看懂，在這裡就不做贅述，遇到的第一個操作init，這個函數的作用是確定有合法解的最小車輛數量，由於直接求解，解空間太大，且有很多車輛不能使用，因此，我們刪去無用的車輛，來縮小解空間（這是一個小優化，能夠加快程序速度）

public BaB_Vrptw init(BaB_Vrptw lp) throws IloException {

lp.build_model();

if (lp.model.solve()) {

lp.get_value();

int aa=0;

for (int i = 0; i < lp.routes.size(); i++) {

if (lp.routes.get(i).size()==2) {

aa++;

}

System.out.println(aa);

if (aa==0) {

data.veh_num -=1;

lp.model.clearModel();

lp = new BaB_Vrptw(data);

return init(lp);

}else {

data.veh_num -=aa;

lp.model.clearModel();

lp = new BaB_Vrptw(data);

return init(lp);

}

}else {

data.veh_num +=1;

System.out.println("vehicle number: "+data.veh_num);

lp.model.clearModel();

lp = new BaB_Vrptw(data);

lp.build_model();

if (lp.model.solve()) {

lp.get_value();

return lp;

}else {

System.out.println("error init");

return null;

}

如上面的程序所示，具體的做法就是建立一個鬆弛了的cplex模型，並計算使用的車輛數，如果有aa輛未使用車輛就減少aa輛可用車輛，否則減少一輛直到沒有可行解。當然，最後我們可使用的車輛是最少的車輛啦~

鬆弛的模型代碼如下，這就是之前「乾貨|十分鐘快速掌握CPLEX求解VRPTW數學模型（附JAVA代碼及CPLEX安裝流程）」中的模型把x_ijk的整數約束去掉得到的。

private void build_model() throws IloException {

//model

model = new IloCplex();

model.setOut(null);

// model.setParam(IloCplex.DoubleParam.EpOpt, 1e-9);

// model.setParam(IloCplex.DoubleParam.EpGap, 1e-9);

//variables

x = new IloNumVar[data.vertex_num][data.vertex_num][data.veh_num];

w = new IloNumVar[data.vertex_num][data.veh_num]; //車輛訪問點的時間

//定義cplex變數x和w的數據類型及取值範圍

for (int i = 0; i < data.vertex_num; i++) {

for (int k = 0; k < data.veh_num; k++) {

w[i][k] = model.numVar(0, 1e15, IloNumVarType.Float, "w" + i + "," + k);

}

for (int j = 0; j < data.vertex_num; j++) {

if (data.arcs[i][j]==0) {

x[i][j] = null;

}

else{

//Xijk,公式(10)-(11)

for (int k = 0; k < data.veh_num; k++) {

x[i][j][k] = model.numVar(0, 1, IloNumVarType.Float, "x" + i + "," + j + "," + k);

}

//加入目標函數

//公式(1)

IloNumExpr obj = model.numExpr();

for(int i = 0; i < data.vertex_num; i++){

for(int j = 0; j < data.vertex_num; j++){

if (data.arcs[i][j]==0) {

continue;

}

for(int k = 0; k < data.veh_num; k++){

obj = model.sum(obj, model.prod(data.dist[i][j], x[i][j][k]));

}

model.addMinimize(obj);

//加入約束

//公式(2)

for(int i= 1; i < data.vertex_num-1;i++){

IloNumExpr expr1 = model.numExpr();

for (int k = 0; k < data.veh_num; k++) {

for (int j = 1; j < data.vertex_num; j++) {

if (data.arcs[i][j]==1) {

expr1 = model.sum(expr1, x[i][j][k]);

}

model.addEq(expr1, 1);

}

//公式(3)

for (int k = 0; k < data.veh_num; k++) {

IloNumExpr expr2 = model.numExpr();

for (int j = 1; j < data.vertex_num; j++) {

if (data.arcs[0][j]==1) {

expr2 = model.sum(expr2, x[0][j][k]);

}

model.addEq(expr2, 1);

}

//公式(4)

for (int k = 0; k < data.veh_num; k++) {

for (int j = 1; j < data.vertex_num-1; j++) {

IloNumExpr expr3 = model.numExpr();

IloNumExpr subExpr1 = model.numExpr();

IloNumExpr subExpr2 = model.numExpr();

for (int i = 0; i < data.vertex_num; i++) {

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]);

}

expr3 = model.sum(subExpr1,model.prod(-1, subExpr2));

model.addEq(expr3, 0);

}

//公式(5)

for (int k = 0; k < data.veh_num; k++) {

IloNumExpr expr4 = model.numExpr();

for (int i = 0; i < data.vertex_num-1; i++) {

if (data.arcs[i][data.vertex_num-1]==1) {

expr4 = model.sum(expr4,x[i][data.vertex_num-1][k]);

}

model.addEq(expr4, 1);

}

//公式(6)

double M = 1e5;

for (int k = 0; k < data.veh_num; k++) {

for (int i = 0; i < data.vertex_num; i++) {

for (int j = 0; j < data.vertex_num; j++) {

if (data.arcs[i][j] == 1) {

IloNumExpr expr5 = model.numExpr();

IloNumExpr expr6 = model.numExpr();

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);

}

//公式(7)

for (int k = 0; k < data.veh_num; k++) {

for (int i = 1; i < data.vertex_num-1; i++) {

IloNumExpr expr7 = model.numExpr();

for (int j = 0; j < data.vertex_num; j++) {

if (data.arcs[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));

}

//公式(8)

for (int k = 0; k < data.veh_num; k++) {

model.addLe(data.E, w[0][k]);

model.addLe(data.E, w[data.vertex_num-1][k]);

model.addLe(w[0][k], data.L);

model.addLe(w[data.vertex_num-1][k], data.L);

}

//公式(9)

for (int k = 0; k < data.veh_num; k++) {

IloNumExpr expr8 = model.numExpr();

for (int i = 1; i < data.vertex_num-1; i++) {

IloNumExpr expr9 = model.numExpr();

for (int j = 0; j < data.vertex_num; j++) {

if (data.arcs[i][j] == 1) {

expr9=model.sum(expr9,x[i][j][k]);

}

expr8 = model.sum(expr8,model.prod(data.demands[i],expr9));

}

model.addLe(expr8, data.cap);

}

之後就是branch and bound過程，在這裡，就是最重點的環節了，先說一下我們的定界方法，把VRPTW的數學模型鬆弛的成一個線性規劃問題可以求解出VRPTW問題的一個下界，分支的原則就是對於一個選定的x_ijk，且0<x_ijk<1，那麼，利用這個x_ijk進行分成兩支，左支是不能夠走弧ij，右支是必須走弧ij且必須由車輛k經過。即左支對於任意的t，x_ijt = 0。右邊則是x_ijk = 1。（關於x_ijk的含義請參考「乾貨|十分鐘快速掌握CPLEX求解VRPTW數學模型（附JAVA代碼及CPLEX安裝流程）」）增加上述約束後，再進行求解，進行定界。找到要分支的弧的代碼如下。

// 找到要分支的弧

public int[] find_arc(double[][][] x) {

int record[] = new int[3];//記錄分支頂點

for (int i = 0; i <data.vertex_num; i++) {

for (int j = 0; j < data.vertex_num; j++) {

if (data.arcs[i][j]>0.5) {

for (int k = 0; k <data.veh_num; k++) {

//若該弧值為0或1，則繼續

if (is_one_zero(x[i][j][k])) {

continue;

}

// cur_dif = get_dif(x[i][j][k]);

record[0] = i;

record[1] = j;

record[2] = k;

return record;

}

record[0] = -1;

record[1] = -1;

record[2] = -1;

return record;

}

分支定界的流程是：

確定一個下界（初始解LB），上界定為無窮大UB。
把初始問題構建一個節點加入優先隊列（因為是優先隊列，所以使用best first sloution，也就是每一次最好的目標值最前搜索）。
判斷隊列是否為空，如果為空跳轉至7，否則取出並彈出隊首元素，計算該節點的目標值P。
如果P > UB，返回3。否則判斷當前節點是否是合法解（對於任意i,j,k,x_ijk均為整數），如果是，跳轉5否則跳轉6。
如果P < UB, 記錄UB = P，當前節點為當前最優解BS。返回3.
設置兩個子節點L, R。L，R的建立方式如上，如果L的目標值L.P <= UB，把L加入隊列，如果R的目標值R.P <= UB，把R加入隊列。返回3.
結束，返回記錄的最優節點BS。如果BS為空則無解。

public void branch_and_bound(BaB_Vrptw lp) throws IloException {

cur_best = 3000;//設置上界

deep=0;

record_arc = new int[3];

node1 = new Node(data);

best_note = null;

queue = new PriorityQueue<Node>();

//初始解（非法解）

for (int i = 0; i < lp.routes.size(); i++) {

ArrayList<Integer> r = lp.routes.get(i);

System.out.println();

for (int j = 0; j < r.size(); j++) {

System.out.print(r.get(j)+" ");

}

lp.copy_lp_to_node(lp, node1);

// node1.node_cost = lp.cost;

// node1.lp_x = lp.x_map.clone();

// node1.node_routes =lp.routes;

// node1.node_servetimes = lp.servetimes;

node2 = node1.note_copy();

deep=0;

node1.d=deep;

queue.add(node1);

//branch and bound過程

while (!queue.isEmpty()) {

Node node = queue.poll();

//某支最優解大於當前最好可行解，刪除

if (doubleCompare(node.node_cost, cur_best)>0) {

continue;

}else {

record_arc = lp.find_arc(node.lp_x);

//某支的合法解,0,1組合的解,當前分支最好解

if (record_arc[0]==-1) {

//比當前最好解好，更新當前解

if (doubleCompare(node.node_cost, cur_best)==-1) {

lp.cur_best = node.node_cost;

System.out.println(node.d+" cur_best:"+cur_best);

lp.best_note = node;

}

continue;

}else {//可以分支

node1 = lp.branch_left_arc(lp, node, record_arc);//左支

node2 = lp.branch_right_arc(lp, node, record_arc);//右支

if (node1!=null && doubleCompare(node1.node_cost, cur_best)<=0) {

queue.add(node1);

}

if (node2!=null && doubleCompare(node2.node_cost, cur_best)<=0) {

queue.add(node2);

}

//設置左支

public Node branch_left_arc(BaB_Vrptw lp,Node father_node,int[] record) throws IloException {

if (record[0] == -1) {

return null;

}

Node new_node = new Node(data);

new_node = father_node.note_copy();

new_node.node_x[record[0]][record[1]] = -1;

for (int k = 0; k < data.veh_num; k++) {

new_node.node_x_map[record[0]][record[1]][k]=0;

}

// new_node.node_x_map[record[0]][record[1]][record[2]]=-1;

//設置左支

lp.set_bound(new_node);

if (lp.model.solve()) {

lp.get_value();

deep++;

new_node.d=deep;

lp.copy_lp_to_node(lp, new_node);

System.out.println(new_node.d+" "+lp.cost);

}else {

new_node.node_cost = data.big_num;

}

return new_node;

}

//設置右支

public Node branch_right_arc(BaB_Vrptw lp,Node father_node,int[] record) throws IloException {

if (record[0] == -1) {

return null;

}

Node new_node = new Node(data);

new_node = father_node.note_copy();

new_node.node_x[record[0]][record[1]] = 1;

// new_node.node_x_map[record[0]][record[1]][record[2]]=1;

for (int k = 0; k < data.veh_num; k++) {

if (k==record[2]) {

new_node.node_x_map[record[0]][record[1]][k]=1;

}else {

new_node.node_x_map[record[0]][record[1]][k]=0;

}

//設置右支

lp.set_bound(new_node);

if (lp.model.solve()) {

lp.get_value();

deep++;

new_node.d=deep;

System.out.println(new_node.d+" right: "+lp.cost);

lp.copy_lp_to_node(lp, new_node);

}else {

new_node.node_cost = data.big_num;

}

return new_node;

}

//找到需要分支的支點位置

這樣就完美的利用branch and bound解決了VRPTW。誒，等等，完美么？是不是忘了點什麼？解的合法性有沒有檢驗呢？

為了檢驗我們所求的解是不是合法的，我們利用遲遲沒出面的Check類來檢查這個問題。

Check類

Check類存在的目的，主要是檢驗解的可行性，包括解是否滿足車輛數量約束，是否滿足容量約束，時間窗約束等等。

包括函數

double_compare(v1, v2): 比較兩個數大小 v1 < v2 – eps 返回 -1， v1 > v2 + eps 返回1，兩數相等返回0。

fesible()：判斷解的可行性，包括車輛數量可行性，車輛載荷可行性，時間窗、車容量可行性判斷。

class Check{

double epsilon = 0.0001;

Data data = new Data();

ArrayList<ArrayList<Integer>> routes = new ArrayList<>();

ArrayList<ArrayList<Double>> servetimes = new ArrayList<>();

public Check(BaB_Vrptw lp) {

super();

this.data = lp.data;

this.routes = lp.routes;

this.servetimes = lp.servetimes;

}

//函數功能：比較兩個數的大小

public int double_compare(double v1,double v2) {

if (v1 < v2 - epsilon) {

return -1;

}

if (v1 > v2 + epsilon) {

return 1;

}

return 0;

}

//函數功能：解的可行性判斷

public void fesible() throws IloException {

//車輛數量可行性判斷

if (routes.size() > data.veh_num) {

System.out.println("error: vecnum!!!");

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.demands[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!!!");

System.exit(0);

}

運算結果

以Solomon測試算例為例，我們對其進行了測試，輸入數據如下：

輸出結果如下：其中第一列代表顧客編號，第二列和第三列分別代表顧客的橫縱坐標，第四列代表需求，第五列第六列代表時間窗，第七列代表服務時間。車輛數量20，容量200。

time 25.245537525 bestcost 827.3

0 13 17 18 19 15 16 14 12 101

0 43 42 41 40 44 46 45 48 51 50 52 49 47 101

0 5 3 7 8 10 11 9 6 4 2 1 75 101

0 67 65 63 62 74 72 61 64 68 66 69 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 57 55 54 53 56 58 60 59 101

0 81 78 76 71 70 73 77 79 80 101

0 90 87 86 83 82 84 85 88 89 91 101

0 98 96 95 94 92 93 97 100 99 101

第一行是運算時間和最優目標值，第二行到最後一行分別表示每輛車的運行路線。

終於搞完了，是不是覺得很舒爽呢？再來張老婆的圖片犒勞一下吧~

-The End-

文案 / 賀興（大三）

排版 / 賀興（大三）

代碼 / 黃楠（博二）（華中科技大學管理學院博士二年級、huangnanhust@163.com)

審核 / 賀興（大三）（華中科技大學管理學院本科三年級、hexing15@gmail.com)

指導老師 / 秦時明岳（professor.qin@qq.com）

如有疑問，歡迎諮詢~

更多最新數據分析相關原創文章，請關注微信公眾號：數據魔術師