Java调用Cplex解VRPTW问题(代码非原创)

参考资料:微信公众号“数据魔法师”;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.)

评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值