干货 | Branch and Price算法求解VRPTW问题（附JAVA代码分享）

  public boolean BBnode(paramsVRP userParam, ArrayList<route> routes,
      treeBB branching, ArrayList<route> bestRoutes, int depth)
      throws IOException {
    // userParam (input) : all the parameters provided by the users (cities,
    // roads...)
    // routes (input) : all (but we could decide to keep only a subset) the
    // routes considered up to now (to initialize the Column generation process)
    // branching (input): BB branching context information for the current node
    // to process (branching edge var, branching value, branching from...)
    // bestRoutes (output): best solution encountered
    int i, j, bestEdge1, bestEdge2, prevcity, city, bestVal;
    double coef, bestObj, change, CGobj;
    boolean feasible;

    try {

      // check first that we need to solve this node. Not the case if we have
      // already found a solution within the gap precision
      if ((upperbound - lowerbound) / upperbound < userParam.gap)
        return true;

      // init
      if (branching == null) { // root node - first call
        // first call - root node
        treeBB newnode = new treeBB();
        newnode.father = null;
        newnode.toplevel = true;
        newnode.branchFrom = -1;
        newnode.branchTo = -1;
        newnode.branchValue = -1;
        newnode.son0 = null;
        branching = newnode;
      }

      // display some local info
      if (branching.branchValue < 1)
        System.out.println("\nEdge from " + branching.branchFrom + " to "
            + branching.branchTo + ": forbid");
      else
        System.out.println("\nEdge from " + branching.branchFrom + " to "
            + branching.branchTo + ": set");
      int mb = 1024 * 1024;
      Runtime runtime = Runtime.getRuntime();
      System.out.print("Java Memory=> Total:" + (runtime.totalMemory() / mb)
          + " Max:" + (runtime.maxMemory() / mb) + " Used:"
          + ((runtime.totalMemory() - runtime.freeMemory()) / mb) + " Free: "
          + runtime.freeMemory() / mb);

      // Compute a solution for this node using Column generation
      columngen CG = new columngen();

      CGobj = CG.computeColGen(userParam, routes);
      // feasible ? Does a solution exist?
      if ((CGobj > 2 * userParam.maxlength) || (CGobj < -1e-6)) { 
        // can only be true when the routes in the solution include forbidden edges (can happen when the BB set branching values)
        System.out.println("RELAX INFEASIBLE | Lower bound: " + lowerbound
            + " | Upper bound: " + upperbound + " | Gap: "
            + ((upperbound - lowerbound) / upperbound) + " | BB Depth: "
            + depth + " | " + routes.size() + " routes");
        return true; // stop this branch
      }
      branching.lowestValue = CGobj;

      // update the global lowerbound when required
      if ((branching.father != null) && (branching.father.son0 != null)
          && branching.father.toplevel) {
        // all nodes above and on the left have been processed=> we can compute
        // a new lowerbound
        lowerbound = (branching.lowestValue > branching.father.son0.lowestValue) ? branching.father.son0.lowestValue
            : branching.lowestValue;
        branching.toplevel = true;
      } else if (branching.father == null) // root node
        lowerbound = CGobj;

      if (branching.lowestValue > upperbound) {
        CG = null;
        System.out.println("CUT | Lower bound: " + lowerbound
            + " | Upper bound: " + upperbound + " | Gap: "
            + ((upperbound - lowerbound) / upperbound) + " | BB Depth: "
            + depth + " | Local CG cost: " + CGobj + " | " + routes.size()
            + " routes");
        return true; // cut this useless branch
      } else {
        // ///////////////////////////////////////////////////////////////////////////
        // check the (integer) feasibility. Otherwise search for a branching
        // variable
        feasible = true;
        bestEdge1 = -1;
        bestEdge2 = -1;
        bestObj = -1.0;
        bestVal = 0;

        // transform the path variable (of the CG model) into edges variables
        for (i = 0; i < userParam.nbclients + 2; i++)
          java.util.Arrays.fill(userParam.edges[i], 0.0);
        for (route r : routes) {
          if (r.getQ() > 1e-6) { // we consider only the routes in the current
                                 // local solution
            ArrayList<Integer> path = r.getpath(); // get back the sequence of
                                                   // cities (path for this
                                                   // route)
            prevcity = 0;
            for (i = 1; i < path.size(); i++) {
              city = path.get(i);
              userParam.edges[prevcity][city] += r.getQ(); // convert into edges
              prevcity = city;
            }
          }
        }

        // find a fractional edge
        for (i = 0; i < userParam.nbclients + 2; i++) {
          for (j = 0; j < userParam.nbclients + 2; j++) {
            coef = userParam.edges[i][j];
            if ((coef > 1e-6)
                && ((coef < 0.9999999999) || (coef > 1.0000000001))) {
              // this route has a fractional coefficient in the solution =>
              // should we branch on this one?
              feasible = false;
              // what if we impose this route in the solution? Q=1
              // keep the ref of the edge which should lead to the largest
              // change
              change = (coef < Math.abs(1.0 - coef)) ? coef : Math.abs(1.0 - coef);
              change *= routes.get(i).getcost();
              if (change > bestObj) {
                bestEdge1 = i;
                bestEdge2 = j;
                bestObj = change;
                bestVal = (Math.abs(1.0 - coef) > coef) ? 0 : 1;
              }
            }
          }
        }
        CG = null;

        if (feasible) {
          if (branching.lowestValue < upperbound) { // new incumbant feasible solution!
            upperbound = branching.lowestValue;
            bestRoutes.clear();
            for (route r : routes) {
              if (r.getQ() > 1e-6) {
                route optim = new route();
                optim.setcost(r.getcost());
                optim.path = r.getpath();
                optim.setQ(r.getQ());
                bestRoutes.add(optim);
              }
            }
            System.out.println("OPT | Lower bound: " + lowerbound
                + " | Upper bound: " + upperbound + " | Gap: "
                + ((upperbound - lowerbound) / upperbound) + " | BB Depth: "
                + depth + " | Local CG cost: " + CGobj + " | " + routes.size()
                + " routes");
            System.out.flush();
          } else
            System.out.println("FEAS | Lower bound: " + lowerbound
                + " | Upper bound: " + upperbound + " | Gap: "
                + ((upperbound - lowerbound) / upperbound) + " | BB Depth: "
                + depth + " | Local CG cost: " + CGobj + " | " + routes.size()
                + " routes");
          return feasible;
        } else {
          System.out.println("INTEG INFEAS | Lower bound: " + lowerbound
              + " | Upper bound: " + upperbound + " | Gap: "
              + ((upperbound - lowerbound) / upperbound) + " | BB Depth: "
              + depth + " | Local CG cost: " + CGobj + " | " + routes.size()
              + " routes");
          System.out.flush();
          // ///////////////////////////////////////////////////////////
          // branching (diving strategy)

          // first branch -> set edges[bestEdge1][bestEdge2]=0
          // record the branching information in a tree list
          treeBB newnode1 = new treeBB();
          newnode1.father = branching;
          newnode1.branchFrom = bestEdge1;
          newnode1.branchTo = bestEdge2;
          newnode1.branchValue = bestVal; // first version was not with bestVal
                                          // but with 0
          newnode1.lowestValue = -1E10;
          newnode1.son0 = null;

          // branching on edges[bestEdge1][bestEdge2]=0
          EdgesBasedOnBranching(userParam, newnode1, false);

          // the initial lp for the CG contains all the routes of the previous
          // solution less the routes containing this arc
          ArrayList<route> nodeRoutes = new ArrayList<route>();
          for (route r : routes) {
            ArrayList<Integer> path = r.getpath();
            boolean accept = true;
            if (path.size() > 3) { // we must keep trivial routes
                                   // Depot-City-Depot in the set to ensure
                                   // feasibility of the CG
              prevcity = 0;
              for (j = 1; accept && (j < path.size()); j++) {
                city = path.get(j);
                if ((prevcity == bestEdge1) && (city == bestEdge2))
                  accept = false;
                prevcity = city;
              }
            }
            if (accept)
              nodeRoutes.add(r);
          }

          boolean ok;
          ok = BBnode(userParam, nodeRoutes, newnode1, bestRoutes, depth + 1);
          nodeRoutes = null; // free memory
          if (!ok) {
            return false;
          }

          branching.son0 = newnode1;

          // second branch -> set edges[bestEdge1][bestEdge2]=1
          // record the branching information in a tree list
          treeBB newnode2 = new treeBB();
          newnode2.father = branching;
          newnode2.branchFrom = bestEdge1;
          newnode2.branchTo = bestEdge2;
          newnode2.branchValue = 1 - bestVal; // first version: always 1
          newnode2.lowestValue = -1E10;
          newnode2.son0 = null;

          // branching on edges[bestEdge1][bestEdge2]=1
          // second branching=>need to reinitialize the dist matrix
          for (i = 0; i < userParam.nbclients + 2; i++)
            System.arraycopy(userParam.distBase[i], 0, userParam.dist[i], 0,
                userParam.nbclients + 2);
          EdgesBasedOnBranching(userParam, newnode2, true);
          // the initial lp for the CG contains all the routes of the previous
          // solution less the routes incompatible with this arc
          ArrayList<route> nodeRoutes2 = new ArrayList<route>();
          for (route r : routes) {
            ArrayList<Integer> path = r.getpath();
            boolean accept = true;
            if (path.size() > 3) { // we must keep trivial routes
                                   // Depot-City-Depot in the set to ensure
                                   // feasibility of the CG
              prevcity = 0;
              for (i = 1; accept && (i < path.size()); i++) {
                city = path.get(i);
                if (userParam.dist[prevcity][city] >= userParam.verybig - 1E-6)
                  accept = false;
                prevcity = city;
              }
            }
            if (accept)
              nodeRoutes2.add(r);

          }
          ok = BBnode(userParam, nodeRoutes2, newnode2, bestRoutes, depth + 1);
          nodeRoutes2 = null;

          // update lowest feasible value of this node
          branching.lowestValue = (newnode1.lowestValue < newnode2.lowestValue) ? newnode1.lowestValue
              : newnode2.lowestValue;

          return ok;

        }
      }

    } catch (IOException e) {
      System.err.println("Error: " + e);
    }
    return false;
  }

  public double computeColGen(paramsVRP userParam, ArrayList<route> routes)
      throws IOException {
    int i, j, prevcity, city;
    double cost, obj;
    double[] pi;
    boolean oncemore;

    try {

      // ---------------------------------------------------------
      // construct the model for the Restricted Master Problem
      // ---------------------------------------------------------
      // warning: for clarity, we create a new cplex env each time we start a
      // Column Generation
      // this class contains (nearly) everything about CG and could be used
      // independently
      // However, since the final goal is to encompass it inside 锟� Branch and
      // Bound (BB),
      // it would (probably) be better to create only once the CPlex env when we
      // initiate the BB and to work with the same (but adjusted) lp matrix each
      // time
      IloCplex cplex = new IloCplex();

      IloObjective objfunc = cplex.addMinimize();

      // for each vertex/client, one constraint (chapter 3, 3.23 )
      IloRange[] lpmatrix = new IloRange[userParam.nbclients];
      for (i = 0; i < userParam.nbclients; i++)
        lpmatrix[i] = cplex.addRange(1.0, Double.MAX_VALUE); 
      // for each constraint, right member >=1
      // lpmatrix[i] = cplex.addRange(1.0, 1.0); 
      // or for each constraint, right member=1 ... what is the best?

      // Declaration of the variables
      IloNumVarArray y = new IloNumVarArray(); // y_p to define whether a path p
                                               // is used

      // Populate the lp matrix and the objective function
      // first with the routes provided by the argument 'routes' of the function
      // (in the context of the Branch and Bound, it would be a pity to start
      // again the CG from scratch at each node of the BB!)
      // (we should reuse parts of the previous solution(s))
      for (route r : routes) {
        int v;
        cost = 0.0;
        prevcity = 0;
        for (i = 1; i < r.getpath().size(); i++) {
          city = r.getpath().get(i);
          cost += userParam.dist[prevcity][city];
          prevcity = city;
        }

        r.setcost(cost);
        IloColumn column = cplex.column(objfunc, r.getcost()); 
        // obj coefficient
        for (i = 1; i < r.getpath().size() - 1; i++) {
          v = r.getpath().get(i) - 1;
          column = column.and(cplex.column(lpmatrix[v], 1.0)); 
          // coefficient of y_i in (3.23) => 0 for the other y_p
        }
        y.add(cplex.numVar(column, 0.0, Double.MAX_VALUE)); 
        // creation of the variable y_i
      }
      // complete the lp with basic route to ensure feasibility
      if (routes.size() < userParam.nbclients) { // a priori true only the first time
        for (i = 0; i < userParam.nbclients; i++) {
          cost = userParam.dist[0][i + 1]
              + userParam.dist[i + 1][userParam.nbclients + 1];
          IloColumn column = cplex.column(objfunc, cost); // obj coefficient
          column = column.and(cplex.column(lpmatrix[i], 1.0)); // coefficient of
                                                               // y_i in (3.23)
                                                               // => 0 for the
                                                               // other y_p
          y.add(cplex.numVar(column, 0.0, Double.MAX_VALUE)); // creation of the
                                                              // variable y_i
          route newroute = new route();
          newroute.addcity(0);
          newroute.addcity(i + 1);
          newroute.addcity(userParam.nbclients + 1);
          newroute.setcost(cost);
          routes.add(newroute);
        }
      }

      // cplex.exportModel("model.lp");

      // CPlex params
      cplex.setParam(IloCplex.IntParam.RootAlg, IloCplex.Algorithm.Primal);
      cplex.setOut(null);
      // cplex.setParam(IloCplex.DoubleParam.TiLim,30); // max number of
      // seconds: 2h=7200 24h=86400

      // ---------------------------------------------------------
      // column generation process
      // ---------------------------------------------------------
      DecimalFormat df = new DecimalFormat("#0000.00");
      oncemore = true;
      double[] prevobj = new double[100];
      int nbroute;
      int previ = -1;
      while (oncemore) {

        oncemore = false;
        // ---------------------------------------s------------------
        // solve the current RMP
        // ---------------------------------------------------------
        if (!cplex.solve()) {
          System.out.println("CG: relaxation infeasible!");
          return 1E10;
        }
        prevobj[(++previ) % 100] = cplex.getObjValue();
        // store the 30 last obj values to check stability afterwards

        // System.out.println(cplex.getStatus());
        // cplex.exportModel("model.lp");

        // ---------------------------------------------------------
        // solve the subproblem to find new columns (if any)
        // ---------------------------------------------------------
        // first define the new costs for the subproblem objective function
        // (SPPRC)
        pi = cplex.getDuals(lpmatrix);
        for (i = 1; i < userParam.nbclients + 1; i++)
          for (j = 0; j < userParam.nbclients + 2; j++)
            userParam.cost[i][j] = userParam.dist[i][j] - pi[i - 1];

        // start dynamic programming
        SPPRC sp = new SPPRC();
        ArrayList<route> routesSPPRC = new ArrayList<route>();

        nbroute = userParam.nbclients; // arbitrarily limit to the 5 first
                                       // shortest paths with negative cost
        // if ((previ>100) &&
        // (prevobj[(previ-3)%100]-prevobj[previ%100]<0.0003*Math.abs((prevobj[(previ-99)%100]-prevobj[previ%100]))))
        // {
        // System.out.print("/");
        // complete=true; // it the convergence is too slow, start a "complete"
        // shortestpast
        // }
        sp.shortestPath(userParam, routesSPPRC, nbroute);
        sp = null;

        // /////////////////////////////
        // parameter here
        if (routesSPPRC.size() > 0) {
          for (route r : routesSPPRC) {
//             if (userParam.debug) {
//             System.out.println(" "+r.getcost());
//             }
            ArrayList<Integer> rout = r.getpath();
            prevcity = rout.get(1);
            cost = userParam.dist[0][prevcity];
            IloColumn column = cplex.column(lpmatrix[rout.get(1) - 1], 1.0);
            for (i = 2; i < rout.size() - 1; i++) {
              city = rout.get(i);
              cost += userParam.dist[prevcity][city];
              prevcity = city;
              column = column.and(cplex.column(lpmatrix[rout.get(i) - 1], 1.0)); 
              // coefficient of y_i in (3.23) => 0 for the other y_p
            }
            cost += userParam.dist[prevcity][userParam.nbclients + 1];
            column = column.and(cplex.column(objfunc, cost));
            y.add(cplex.numVar(column, 0.0, Double.MAX_VALUE,
                "P" + routes.size())); // creation of the variable y_i
            r.setcost(cost);
            routes.add(r);

            oncemore = true;
          }
          System.out.print("\nCG Iter " + previ + " Current cost: "
                  + df.format(prevobj[previ % 100]) + " " + routes.size()
                  + " routes");
            System.out.flush();
        }
        //if (previ % 50 == 0)
        routesSPPRC = null;
      }

      System.out.println();

      for (i = 0; i < y.getSize(); i++)
        routes.get(i).setQ(cplex.getValue(y.getElement(i)));
      obj = cplex.getObjValue(); // mmmmhhh: to check. To be entirely safe, we
                                 // should recompute the obj using the distBase
                                 // matrix instead of the dist matrix

      cplex.end();
      return obj;
    } catch (IloException e) {
      System.err.println("Concert exception caught '" + e + "' caught");
    }
    return 1E10;
  }

  public void shortestPath(paramsVRP userParamArg, ArrayList<route> routes,int nbroute) {
    label current;
    int i,j,idx,nbsol,maxsol;
    double d,d2;
    int[]  checkDom;
    float tt,tt2;
    Integer currentidx;

    this.userParam=userParamArg;
    // unprocessed labels list => ordered TreeSet List (?optimal:  need to be sorted like this?)
    TreeSet<Integer> U = new TreeSet<Integer>(new MyLabelComparator());   // unprocessed labels list

    // processed labels list => ordered TreeSet List 
    TreeSet<Integer> P = new TreeSet<Integer>(new MyLabelComparator());   // unprocessed labels list

    // array of labels
    labels = new ArrayList<label>(2*userParam.nbclients); // initial size at least larger than nb clients
    boolean[] cust= new boolean[userParam.nbclients+2];
    cust[0]=true;
    for (i=1;i<userParam.nbclients+2;i++)
      cust[i]=false;
    labels.add(new label(0,-1,0.0,0,0,false,cust));  // first label: start from depot (client 0)
    U.add(0);

    // for each city, an array with the index of the corresponding labels (for dominance)
    checkDom = new int[userParam.nbclients+2];
    ArrayList<Integer>[] city2labels = new ArrayList[userParam.nbclients+2];
    for (i=0;i<userParam.nbclients+2;i++) {
      city2labels[i]=new ArrayList<Integer>();
      checkDom[i]=0;  // index of the first label in city2labels that needs to be checked for dominance (last labels added)
    }
    city2labels[0].add(0);

    nbsol = 0;
    maxsol = 2 * nbroute;
    while ((U.size()>0) && (nbsol<maxsol)) {  
      // second term if we want to limit to the first solutions encountered to speed up the SPPRC (perhaps not the BP)
      // remark: we'll keep only nbroute, but we compute 2xnbroute!  It makes a huge difference=>we'll keep the most negative ones
      // this is something to analyze further!  how many solutions to keep and which ones?
      // process one label => get the index AND remove it from U
      currentidx = U.pollFirst();
      current = labels.get(currentidx);

      // check for dominance
      // code not fully optimized: 
      int l1,l2;
      boolean pathdom;
      label la1,la2;
      ArrayList<Integer> cleaning = new ArrayList<Integer>();
      for (i = checkDom[current.city]; i < city2labels[current.city].size(); i++) {  // check for dominance between the labels added since the last time we came here with this city and all the other ones
        for (j = 0; j < i; j++) {
          l1 = city2labels[current.city].get(i);
          l2 = city2labels[current.city].get(j);
          la1 = labels.get(l1);
          la2 = labels.get(l2);
          if (!(la1.dominated || la2.dominated)) {  // could happen since we clean 'city2labels' thanks to 'cleaning' only after the double loop
            pathdom = true;
            for (int k = 1; pathdom && (k < userParam.nbclients+2); k++) 
              pathdom=(!la1.vertexVisited[k] || la2.vertexVisited[k]);
            if (pathdom && (la1.cost<=la2.cost) && (la1.ttime<=la2.ttime) && (la1.demand<=la2.demand)) {
              labels.get(l2).dominated = true;
              U.remove((Integer) l2);
              cleaning.add(l2);
              pathdom = false;
              //System.out.print(" ###Remove"+l2);
            } 
            pathdom = true;
            for (int k = 1; pathdom && (k < userParam.nbclients + 2); k++) 
              pathdom = (!la2.vertexVisited[k] || la1.vertexVisited[k]);

            if (pathdom && (la2.cost<=la1.cost) && (la2.ttime<=la1.ttime) && (la2.demand<=la1.demand)) {
              labels.get(l1).dominated = true;
              U.remove(l1);
              cleaning.add(l1);
              //System.out.print(" ###Remove"+l1);
              j = city2labels[current.city].size();
            }
          }
        }
      }

      for (Integer c : cleaning) 
        city2labels[current.city].remove((Integer) c);   // a little bit confusing but ok since c is an Integer and not an int!

      cleaning = null;
      checkDom[current.city] = city2labels[current.city].size();  // update checkDom: all labels currently in city2labels were checked for dom.

      // expand REF
      if (!current.dominated){ 
        //System.out.println("Label "+current.city+" "+current.indexPrevLabel+" "+current.cost+" "+current.ttime+" "+current.dominated);
        if (current.city == userParam.nbclients + 1) { // shortest path candidate to the depot!
          if (current.cost<-1e-7)  {        // SP candidate for the column generation
            P.add(currentidx);
            nbsol=0;
            for (Integer labi : P) {
              label s = labels.get(labi);
              if (!s.dominated)
                nbsol++;
            }
          }
        } else {  // if not the depot, we can consider extensions of the path
          for (i = 0; i < userParam.nbclients + 2; i++) {
            if ((!current.vertexVisited[i]) && (userParam.dist[current.city][i] < userParam.verybig-1e-6)) {  // don't go back to a vertex already visited or along a forbidden edge
              // ttime
              tt = (float) (current.ttime + userParam.ttime[current.city][i] + userParam.s[current.city]);
              if (tt < userParam.a[i])
                tt = userParam.a[i];
              // demand
              d = current.demand + userParam.d[i];
              //System.out.println("  -- "+i+" d:"+d+" t:"+tt);

              // is feasible?
              if ((tt <= userParam.b[i]) && (d <= userParam.capacity)) {
                idx = labels.size();
                boolean[] newcust = new boolean[userParam.nbclients + 2];
                System.arraycopy(current.vertexVisited, 0, newcust, 0, userParam.nbclients + 2);
                newcust[i] = true;
                //speedup: third technique - Feillet 2004 as mentioned in Laporte's paper
                for (j=1;j<=userParam.nbclients;j++)
                  if (!newcust[j]) {
                    tt2=(float) (tt+userParam.ttime[i][j]+userParam.s[i]);
                    d2=d+userParam.d[j];
                    if ((tt2>userParam.b[j]) || (d2>userParam.capacity))
                      newcust[j]=true;  // useless to visit this client
                  }

                labels.add(new label(i, currentidx, current.cost+userParam.cost[current.city][i], tt, d, false, newcust));  // first label: start from depot (client 0)
                if (!U.add((Integer) idx)) {  
                  // only happens if there exists already a label at this vertex with the same cost, time and demand and visiting the same cities before
                  // It can happen with some paths where the order of the cities is permuted
                  labels.get(idx).dominated = true; // => we can forget this label and keep only the other one
                } else 
                  city2labels[i].add(idx);

              }
            }
          }
        }
      }
    }
    // clean
    checkDom = null;

    // filtering: find the path from depot to the destination
    Integer lab;
    i = 0;
    while ((i < nbroute) && ((lab = P.pollFirst()) != null)) {
      label s = labels.get(lab);
      if (!s.dominated) {
        if (/*(i < nbroute / 2) ||*/ (s.cost < -1e-4)) {
          // System.out.println(s.cost);
//          if(s.cost > 0) {
//            System.out.println("warning >>>>>>>>>>>>>>>>>>>>");
//          }
          route newroute = new route();
          newroute.setcost(s.cost);
          newroute.addcity(s.city);
          int path = s.indexPrevLabel;
          while (path >= 0) {
            newroute.addcity(labels.get(path).city);
            path = labels.get(path).indexPrevLabel;
          }
          newroute.switchpath();
          routes.add(newroute);
          i++;
        }
      }

    }
  }