基於branch and bound插入的large neighborhood search
一、前言
今年開年那會還在做一個課題的實驗，那時候想用large neighborhood search來做一個問題，但是後來發現常規的一些repair、destroy運算元效果並不是很好。後來才知道，large neighborhood search以及它的衍生演算法，這類框架給人一種非常通用的感覺，就是無論啥問題都能往裡面套。
往往的結果是套進去效果也是一般。這也是很多剛入行的小夥伴經常喜歡乾的事吧，各種演算法框架套一個問題，發現結果不好了就感覺換下一個。最後復現了N多個演算法發現依然no process，這時候就會懷疑人生了。其實要想取得好的performance，肯定還是要推導一些問題特性，設計相應的運算元也好，鄰域結構也好。
好了，回到正題。當時我試了好幾個large neighborhood search運算元，發現沒啥效果的時候，心裡難受得很。那幾天晚上基本上是轉輾反側，難以入眠，當然了是在思考問題。然後一個idea突然浮現在我的腦瓜子裡，常規的repair運算元難以在問題中取得好的performance，是因為約束太多了，插入的時候很容易違背約束。在不違背約束的條件下又難以提升解的質量，我就想能不能插入的啥時候採取branch and bound。遍歷所有的可能插入方式，然後記錄過程中的一個upper bound用來刪掉一些分支。
感覺是有搞頭的，後來想想，這個branch的方法以及bound的方法似乎是有點難設計。然後又擱置了幾天，最後沒進展的時候突然找了一篇論文，是好多年前的一篇文章了。裡面詳細講解了large neighborhood search中如何利用branch and bound進行插入，後來實現了以下感覺還可以。感覺這個方法還是有一定的參考價值的，因此今天就來寫寫（其實當時就想寫了，只不過一直拖到了現在。。。）
二、large neighborhood search
關於這個演算法，我在此前的推文中已經有過相應的介紹，詳情小夥伴們可以戳這篇的連結進行檢視：
自適應大鄰域搜尋(Adaptive Large Neighborhood Search)入門到精通超詳細解析-概念篇
我把其中的一段話摘出來：
大多數鄰域搜尋演算法都明確定義它們的鄰域。 在LNS中，鄰域是由\(destroy\)和\(repair\)方法隱式定義的。\(destroy\)方法會破壞當前解的一部分，而後\(repair\)方法會對被破壞的解進行重建。\(destroy\)方法通常包含隨機性的元素，以便在每次呼叫\(destroy\)
那麼，解\(x\)的鄰域\(N(x)\)就可以定義為：首先通過利用\(destroy\)方法破壞解\(x\)，然後利用\(repair\)方法重建解\(x\)，從而得到的一系列解的集合。LNS演算法框架如下：
有關該演算法更詳細的介紹可以參考Handbook Of Metaheuristics這本書2019版本中的Chapter 4
Large Neighborhood Search（David Pisinger and Stefan Ropke），文末我會放出下載的連結。
關於destroy運算元呢，有很多種，比如隨機移除幾個點，貪心移除一些比較差的點，或者基於後悔值排序移除一些點等，這裡我給出文獻中的一種移除方式，Shaw (1998)提出的基於\(relateness\)進行移除：
假設需要從解中所有的\(Customers\)移除\(n\)個，它首先隨機選擇一個\(Customer\)放進\(S\)（已經移除的\(Customer\)列表）（第1行），然後迭代地（3–6行）移除剩下的\(n-1\)個\(Customer\)。 每次這樣的迭代都會先從\(S\)中隨機選擇一個\(Customer\)，並根據相關標準對其餘未移除的\(Customer\)進行排序（第3-4行）。 在第5行中計算要插入的新\(Customer\)的下標，然後插入到\(S\)中（第6行），直到迭代結束。 關聯度的定義如Shaw（1998）所述：
\[relateness(i, j) = \frac{1}{c'_{ij}+v_{ij}} \]其中，customer \(i\)和 \(j\)在不同的路徑中時\(v_{ij}=1\)，否則為0。
三、branch and bound
上面講了Large Neighborhood Search以及介紹了一個\(destroy\)方法，下面就是重頭戲，如何利用branch and bound進行插入了。
3.1 branch
其實插入的分支方式還是挺好設計的，這玩意兒呢我將也比較難講清楚，我就畫圖好了，還是基於VRP問題示例，其他問題類似，假如我們現在有這樣一個解\(s\)：
為了演示我就不畫太多點太多路徑了，免得大家看得心累。
紅色箭頭就是能夠插入的位置。現在，假如我們插入\(2\)（由於branch and bound是需要遍歷所有可能的插入組合，因此先插入哪個後插入哪個都是可以的，但是分支定界的速度可能會受到很大的影響，這個咱們暫時不討論）：
為了讓大家看得更加清楚，我把\(2\)的位置用粉紅色給標記出來了，一共有3條分支，有個候選的位置就有多少條分支。
現在，還剩下\(4\)，插入\(4\)的時候，我們需要繼續進行分支：
55~畫分支樹真是畫死我啦（大家一定要給個贊，點個在看呀~），可以看到，最後每條路徑就是一個完成的解。在兩個點的\(partial\, solution\)插入兩個點最後分支完成的\(completed\, solution\)居然有12個！！！大家可以自行腦補下，在90個點的\(partial\, solution\)中插入10個點最終形成的分支會有多少。毫無疑問會很多很多，多到你無法想象。下面是DFS搜尋分支樹的過程：
如果要插入的客戶組為空，則可以認為所有客戶已經插入到solution中，形成了一個\(completed\, solution\)，因此判斷找到的一個upper bound是否比最優的upper bound還要好，是的話就對upper bound進行更新。 否則，它會選擇插入效果最好的客戶，這會使目標函式下降得最大（Shaw 1998中也使用了這種啟發式方法）。 然後，對所有插入客戶後形成的分支按照lower bound進行排序，從lower bound低的分支開始繼續往下分支（可以算是一種加速的策略）。 同樣，請注意，該演算法僅探索其lower bound比upper bound更好的分支。
3.2 bound
開始之前大家想想bound的難點在哪裡呢？首先想想bound中最重要的兩個界：upper bound和lower bound：
- lower bound是指搜尋過程中一個partial solution（比如上圖插入\(2\)後形成的3個\(partial\, solution\)）的目標值，因為\(partial\, solution\)並不能算完整的一個解，繼續往下的時候只可能增加（最小化問題）或者減少（最大化問題），因此它的意思是說當前支路的最終形成解的目標值下界（最終目標值不可能比這個lower bound更好）。
- upper bound是指搜尋過程中找到的一個feasible solution（比如上圖插入\(4\)後形成的12個\(completed\, solution\)中滿足所有約束的就是\(feasible\, solution\)）的目標值，如果存在某支路的lower bound比upper bound還要差，那麼該支路顯然是沒有繼續往下的必要了，可以剪去。
顯然可以使用LNS在destroy之前的解的目標值作為upper bound，因為我們總是期望找到比當前解更好的解，才會去進行destroy和repair。現在的問題是如何對一個\(partial\, solution\)的lower bound應該怎樣計算。下面講講幾種思路：
(1) 文獻中給出的思路，利用最小生成樹：
這個方案我試了，但是找到的lower bound實在是太低了，這個lower bound只考慮了距離因素，但問題中往往還存在時間窗等約束。因此這個方法在我當時做的問題中只能說聊勝於無。
(2) 按照greedy的方法將所有未插入的Customer插入到他們最好的位置上，形成一個\(completed\, solution\)，然後該\(completed\, solution\)的目標值作為lower bound。
但是這個lower bound是有缺陷的，因為很難保證不會錯過某些比較有潛力的分支。
(3) 直接利用當前的\(partial\, solution\)的目標值作為lower bound，也比較合理。但是該值往往太低了，這可能會導致要遍歷更多的分支，消耗更多時間。
以上就是一些思路，至於有沒有更好的bound方法，我後面也沒有往下深究了。當時實現出來以後效果是有的，就是時間太長了，然後也放棄了。
當然這篇paper後面也給了一個利用LDS進行搜尋以加快演算法的速度，這裡就不展開了，有空再說。感興趣的小夥伴可以去看看原paper，我會放到留言區的。
四、程式碼環節
程式碼實現放兩個，一個是我當時寫的一個DFSEXPLORER，採用的是思路2作為bound的，（程式碼僅僅提供思路）如下：
private void DFSEXPLORER5(LNSSolution node, LNSSolution upperBound, int dep) {
Queue<LNSSolution> queue = new LinkedList<LNSSolution>();
LNSSolution s_c_ = node;
queue.add(s_c_);
int es = 1;
while (!queue.isEmpty()) {
s_c_ = queue.remove();
//v是一個完整的解
if(s_c_.removalCustomers.isEmpty()) {
if(s_c_.cost < upperBound.cost && Math.abs(s_c_.cost-upperBound.cost)>0.001) {
//System.out.println("new found > "+s_c_.cost+" feasible = "+s_c_.feasible());
upperBound.cost = s_c_.cost;
upperBound.routes = s_c_.routes;
}
}else {
//System.out.println("l > "+s_c_.removalCustomers.size() + " cost = "+s_c_.cost);
double minIDelta = Double.POSITIVE_INFINITY;
int minIndex = -1;
Customer c=null;
for(int i = 0; i < s_c_.removalCustomers.size(); ++i) {
Customer cu = s_c_.removalCustomers.get(i);
double d1 = s_c_.minInsertionDeltas[cu.getCustomerNo()];
if(minIDelta > d1) {
minIDelta = d1;
c = cu;
minIndex = i;
}
}
ArrayList<LNSSolution> neighborI_c = new ArrayList<LNSSolution>();
for( int i = 0; i < s_c_.routes.length; ++i) {
Route route = s_c_.routes[i];
if(!MyUtil.checkCompatibility(c, route.getAssignedVehicle())) {
continue;
}
for (int j = 0; j <= route.getCustomersLength(); j++) {
LNSSolution s_i = s_c_.solClones();
s_i.insertCustomer(s_i.routes[i], s_i.removalCustomers.get(minIndex), j, minIndex);
//updateIDAfterOneInserted(s_i, s_i.routes[i]);
//s_i.calcLowerBound();
double o_c = s_i.lb;
updateInsertionDelta(s_i);
double n_c = s_i.lb;
//if(o_c != n_c)System.out.println("o = "+o_c+" n = "+n_c);
neighborI_c.add(s_i);
}
}
Collections.sort(neighborI_c);
for(LNSSolution s:neighborI_c) {
//System.out.println("lBound "+s.lb+" upperBound = "+upperBound.cost);
//updateInsertionDelta(s);
//s.calcLowerBound();
if(s.lb < upperBound.cost /*&& dep > 0*/) {
//System.out.println("lBound "+s.lb+" upperBound = "+upperBound.cost);
//System.out.println(s.removalCustomers.size());
queue.add(s);
es++;
dep--;
}
}
}
}
//System.out.println(es);
}
第二個是GitHub上找到的一個人復現的，我已經fork到我的倉庫中了：
https://github.com/dengfaheng/vrp
這個思路bound的思路呢沒有按照paper中的，應該還是用的貪心進行bound。看起來在R和RC系列的算例中效果其實也一般般，因為用了LDS吧可能。下面是執行的c1_2_1的截圖：
匯入idea或者eclipse後等他安裝完依賴，執行下面的檔案即可，更改算例的位置如圖所示：
這個思路是直到借鑑的，大家在用LNS的時候也可以想想有什麼更好的bound方法。
欲下載本文相關的完整程式碼及算例，請關注公眾號【程式猿聲】，檢視相應推文留言區即可