1、前言

記得公衆號好久以前推出過一個branch and price的概念推文，後來小編找到了部分（不完整）的代碼，通過研究之後補齊了這部分代碼，可以運行之後也分享了給你們。詳情能夠看：java

乾貨 | 10分鐘帶你掌握branch and price（分支訂價）算法超詳細原理解析）node

乾貨 | Branch and Price算法求解VRPTW問題（附JAVA代碼分享）web

前陣子一個學姐問我這個代碼，我看了一下而後發現本身看不懂了……這幾天恰好要作相關的課題，又回來惡補了一番。把代碼的邏輯梳理了一遍，今天就寫寫，方便你們學習。算法

2、branch and price

branch and price實際上是column generation和branch and bound的結合。爲何要結合呢？前面的文章中介紹過，column generation是求解大規模線性規劃問題的，注意是線性規劃問題，不是整數或者混合整數規劃問題。編程

所以在用column generation求解VRPTW的Set-Partitioning Model的時候，須要先將決策變量由整數鬆弛爲實數，詳情參見：微信

乾貨 | 10分鐘教你使用Column Generation求解VRPTW的線性鬆弛模型app

那麼問題來了，既然column generation只能求解線性鬆弛模型，求得的解若是剛好是整數解那還好說。可是大部分狀況獲得的多是一個非整數解，這顯然不是咱們VRPTW問題的最終解。框架

這個時候就須要用到branch and bound去找整數解了。我記得前面的推文也提到過，branch and bound自己是一個比較通用的算法，只要設計好問題的branch和bound，便可用來求解該問題。所以branch and bound也能夠用來求解整數規劃問題，詳情能夠參見：less

乾貨 | 10分鐘帶你全面掌握branch and bound（分支定界）算法-概念篇編輯器

乾貨 | 10分鐘搞懂branch and bound算法的代碼實現附帶java代碼

branch and bound在求解整數規劃模型的時候經過對當前模型加約束的方式，來對決策變量進行分支，而支路的lower bound能夠經過求解該支路下整數規劃模型的線性鬆弛而得到。求解通常線性規劃問題能夠用單純形法，可是VRPTW的Set-Partitioning Model線性鬆弛後是超大規模的線性規劃問題，無法枚舉全部列。所以採用column generation來進行求解。

不知道你們暈了沒有呢！暈的話能夠把上面的推文好好看一遍，這些都是通過整理的，能節省大家再去瞎找資料亂學習所須要的時間成本。

3、branch and bound

關於VRPTW的Set-Partitioning Model是這樣的：

決策變量爲整數，我當時想着按照branch and bound求解整數規劃模型的思路看代碼，結果半天沒看出個因此然來。由於我沒找到的相關分支代碼。而後看着看着就懵逼了。最後我問了下學長怎麼進行分支，他一句話點醒了我：你能夠按照Set-Partitioning Model進行分支，也能夠按照以前的arc-flow model進行分支。

瞬間醒悟，代碼中的分支方式也是採用arc-flow model進行分支的，那麼這兩個是怎樣關聯起來的呢？實際上是靠"邊"，Set-Partitioning Model中的一列其實就表明一條路徑：

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;