运筹学教学|分支定界法解带时间窗的车辆路径规划问题(附代码及详细注释)

历尽千辛万苦,外加外援帮助,本辣鸡小编终于搞定了这个大坑-用分支定界法(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 : 定义分支定界的节点

01

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

02

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

03

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

分支定界的流程是:

  1. 确定一个下界(初始解LB),上界定为无穷大UB。
  2. 把初始问题构建一个节点加入优先队列(因为是优先队列,所以使用best first sloution,也就是每一次最好的目标值最前搜索)。
  3. 判断队列是否为空,如果为空跳转至7,否则取出并弹出队首元素,计算该节点的目标值P。
  4. 如果P > UB,返回3。否则判断当前节点是否是合法解(对于任意i,j,k,x_ijk均为整数),如果是,跳转5否则跳转6。
  5. 如果P < UB, 记录UB = P,当前节点为当前最优解BS。返回3.
  6. 设置两个子节点L, R。L,R的建立方式如上,如果L的目标值L.P <= UB,把L加入队列,如果R的目标值R.P <= UB,把R加入队列。返回3.
  7. 结束,返回记录的最优节点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类来检查这个问题。

01

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

第一行是运算时间和最优目标值,第二行到最后一行分别表示每辆车的运行路线。

原文发布于微信公众号 - 数据魔术师(data-magician)

原文发表时间:2018-05-10

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏Java成长之路

Java内存模型

多任务处理在现代计算机操作系统中几乎已经是一项必备的功能了。计算机cpu的运算速度与它的存储和通信子系统速度的差距太大,大量的时间都花费在磁盘I/O、网络通信或...

571
来自专栏Lambda

JavaScript排序算法详解

JS家的排序算法 引子 有句话怎么说来着: 雷锋推倒雷峰塔,Java implements JavaScript. 当年,想凭借抱Java大腿火一把...

2187
来自专栏Golang语言社区

[转载]Go JSON 技巧

相对于很多的语言来说, Go 的 JSON 解析可谓简单至极. 问题 通常情况下, 我们在 Go 中经常这样进行 JSON 的解码: package main ...

2913
来自专栏阿凯的Excel

文本数字拆分技巧(第二弹!)

上期刚刚分享了简单的通过智能填充和Len与LenB函数实现的文本数字拆分! 感兴趣可以点我先看上一期的! 本期难度较上期略有提高,和您分享新的技巧。 ? 没...

2797
来自专栏Golang语言社区

Golang-简洁的并发

转载原文:http://www.yankay.com/go-clear-concurreny/ 多核处理器越来越普及。有没有一种简单的办法,能够让我们写的软件释...

3364
来自专栏生信技能树

可能只是一个函数,却要耗费你大半天

好像不少人问过我一个聚类后的树如何根据肉眼观察到的cluster情况来提前指定的树的子集,有点类似于WGCNA分析把几千个基因划分成若干个module后能提取各...

1073
来自专栏C/C++基础

进程与线程的区别

在开发工作中,尤其是对负载较大的服务端程序的开发,为充分发挥处理器多核性能,提高硬件资源利用率,增加系统吞吐量,少不了并发编程。并发编程一般通过多进程和多线程的...

673
来自专栏潇涧技术专栏

Python Algorithms - C6 Divide and Combine and Conquer

Python算法设计篇(6) Chapter 6: Divide and Combine and Conquer

982
来自专栏C/C++基础

C/C++获取本地时间常见方法

(1)UTC (Coordinated Universal Time):协调世界时,又称世界标准时间。曾由格林威治平均时间(Greenwich Mean Tim...

591
来自专栏冰霜之地

如何设计并实现一个线程安全的 Map ?(上篇)

Map 是一种很常见的数据结构,用于存储一些无序的键值对。在主流的编程语言中,默认就自带它的实现。C、C++ 中的 STL 就实现了 Map,JavaScrip...

762

扫码关注云+社区