其他分享
首页 > 其他分享> > branch and price求解VRPTW问题代码详解

branch and price求解VRPTW问题代码详解

作者:互联网

一、前言

记得公众号很久之前推出过一个branch and price的概念推文,后来小编找到了部分(不完整)的代码,经过研究以后补齐了这部分代码,能够运行以后也分享了给大家。详情可以看:

干货 | 10分钟带你掌握branch and price(分支定价)算法超详细原理解析)

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

前阵子一个学姐问我这个代码,我看了一下然后发现自己看不懂了……这几天刚好要做相关的课题,又回来恶补了一番。把代码的逻辑梳理了一遍,今天就写写,方便大家学习。

二、branch and price

branch and price其实是column generation和branch and bound的结合。为什么要结合呢?前面的文章中介绍过,column generation是求解大规模线性规划问题的,注意是线性规划问题,不是整数或者混合整数规划问题。

因此在用column generation求解VRPTW的Set-Partitioning Model的时候,需要先将决策变量由整数松弛为实数,详情参见:

干货 | 10分钟教你使用Column Generation求解VRPTW的线性松弛模型

那么问题来了,既然column generation只能求解线性松弛模型,求得的解如果恰好是整数解那还好说。但是大部分情况得到的可能是一个非整数解,这显然不是我们VRPTW问题的最终解。

这个时候就需要用到branch and bound去找整数解了。我记得前面的推文也提到过,branch and bound本身是一个比较通用的算法,只要设计好问题的branch和bound,即可用来求解该问题。因此branch and bound也可以用来求解整数规划问题,详情可以参见:

干货 | 10分钟带你全面掌握branch and bound(分支定界)算法-概念篇

干货 | 10分钟搞懂branch and bound算法的代码实现附带java代码

branch and bound在求解整数规划模型的时候通过对当前模型加约束的方式,来对决策变量进行分支,而支路的lower bound可以通过求解该支路下整数规划模型的线性松弛而获得。求解一般线性规划问题可以用单纯形法,但是VRPTW的Set-Partitioning Model线性松弛后是超大规模的线性规划问题,没法枚举所有列。因此采用column generation来进行求解。

不知道大家晕了没有呢!晕的话可以把上面的推文好好看一遍,这些都是经过整理的,能节省你们再去瞎找资料乱学习所需要的时间成本。

三、branch and bound

关于VRPTW的Set-Partitioning Model是这样的:

watermark,size_16,text_QDUxQ1RP5Y2a5a6i,color_FFFFFF,t_100,g_se,x_10,y_10,shadow_90,type_ZmFuZ3poZW5naGVpdGk=

决策变量为整数,我当时想着按照branch and bound求解整数规划模型的思路看代码,结果半天没看出个所以然来。因为我没找到的相关分支代码。然后看着看着就懵逼了。最后我问了下学长怎么进行分支,他一句话点醒了我:你可以按照Set-Partitioning Model进行分支,也可以按照之前的arc-flow model进行分支。

瞬间醒悟,代码中的分支方式也是采用arc-flow model进行分支的,那么这两个是怎样关联起来的呢?其实是靠"边",Set-Partitioning Model中的一列其实就代表一条路径:

watermark,size_16,text_QDUxQ1RP5Y2a5a6i,color_FFFFFF,t_100,g_se,x_10,y_10,shadow_90,type_ZmFuZ3poZW5naGVpdGk=

for example,在上面的模型中就包含了4列,每一列就代表一条路径:

r1 : 1->4

r2 : 2->3

r3 : 1->3

r4 : 1->2->3->4

arc-flow model中的“边”是不是包含在上面的路径中了呢,当Set-Partitioning Model中的(k=1,2,3,4)不为整数时,那么这些路径中相应的“边”在arc-flow model也不是整数了。这样两个模型就可以关联起来了。因此我们在branch and bound的时候也可以按照“边”来进行分支。

好了,讲了这么多,我们来看看代码是怎么操作的吧!当然了我这里也只是讲讲大概的思路,详细的还是自己慢慢去研读代码哈。

这里分支树的搜索采用的是递归的方式,遍历方式是有一个公式进行下一条“边”的选择的:

(1) 首先判断上下界是否达到指定的gap,如果是,则认为已经完成了搜索,跳出递归。

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


(2) 分支节点为null时,说明为初始状态,那么分配根节点,开始进行分支。

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

(3) branchValue<1时,该边被禁止,否则该边被选择(一定要经过)。

// 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");
}

(4) 利用column generation计算该节点的lower bound,注意上面有些边被禁止了,因此这些边在column generation中是无法被访问或者强制要经过的。

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

CGobj = CG.computeColGen(userParam, routes);

(5) 新的lower bound出来后,看看要不要更新全局的lower bound(其实这里维护一个全局最小的lower bound也没啥用,就是拿来输出看看gap的):

// 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 = Math.min(branching.lowestValue, branching.father.son0.lowestValue);
        branching.toplevel = true;
    } else if (branching.father == null) {
        // root node
        lowerbound = CGobj;
    }

(6) 接下来是常规操作,判断lower bound和upper bound的大小,进行剪枝或者再次分支。如果该分支的lower bound > upperbound,那么cut掉:

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
}

(7) 否则,先判断column generation找到的解是不是整数解(所有的边都不能有小数):

// ///////////////////////////////////////////////////////////////////////////
// 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 = Math.min(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;
            }
        }
    }
}

(8) 如果是整数解,那么就找到了一个新的可行解,判断是否更新upper bound:

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

(9) 否则,找一条边继续进行分支(这条边具体的选择也会影响分支的速度,选择看步骤(7)中的 bestEdge1和 bestEdge2),EdgesBasedOnBranching函数的作用是通过设置距离矩阵各边的距离,来禁止或者指定选择一些边(比如将一条边的距离设置为正无穷,那么该边就无法访问了)。这里是先分左支,即该边被禁止,要移除column generation的RLMP中包含该边的所有路径:

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

(10) 然后是分右支,该支限定该边一定要经过,因此要要移除column generation的RLMP中不包含该边的所有路径:

// 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);
}
//reinitialize了因此需要recur递归一下
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 = Math.min(newNode1.lowestValue, newNode2.lowestValue);

return ok;

branch and bound的代码就到这了,是不是非常简单呢!细节上碍于篇幅我就不多讲了。大家可以慢慢看代码,不懂在留言区提出来就好了。

四、column generation

这部分分为Master problem和pricing problem,这两部分的内容公众号已经介绍过了。Master problem的代码基本上是固定的框架,直接拿过来用就好了。而pricing problem的代码用的是labeling的算法。

注意pricing problem找路径时,是要结合之前branch and bound禁止or已经选择的边进行最短路的寻找,关于这部分的内容可以参照:

干货 | VRPTW子问题ESPPRC的介绍及其求解算法的C++代码

最短路问题与标号算法(label correcting algorithm)研究(1) - 开篇介绍

最短路问题与标号算法(label correcting algorithm)研究(2) - 最短路径问题简介

最短路问题与标号算法(label correcting algorithm)研究(3)

最短路问题与标号算法(label correcting algorithm)研究(4)

五、代码下载

参见这篇推文:

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

这篇推文包含的内容不多,但是如果你刚好在学习branch and price算法,不妨看看这个,看看代码,相信对你会有所帮助的哦。国内这块的资料和代码都太少了,大佬们的代码又长又臭,压根看不下去。

watermark,size_16,text_QDUxQ1RP5Y2a5a6i,color_FFFFFF,t_100,g_se,x_10,y_10,shadow_90,type_ZmFuZ3poZW5naGVpdGk=

 

看在我写了这么多的份上,能不能帮我点个在看呀~

推荐阅读:

干货 | 想学习优化算法,不知从何学起?

干货 | 运筹学从何学起?如何快速入门运筹学算法?

干货 | 学习算法,你需要掌握这些编程基础(包含JAVA和C++)

干货 | 算法学习必备诀窍:算法可视化解密

干货 | 模拟退火、禁忌搜索、迭代局部搜索求解TSP问题Python代码分享

标签:branching,price,bound,VRPTW,upperbound,branch,routes,代码
来源: https://blog.51cto.com/u_14328065/2919972