Outline
- 前言
- Example-1
- Example-2
- 運行說明
00 前言
前面一篇文章我們講了branch and bound算法的相關概念。可能大家對精確算法實現的印象大概只有一個，調用求解器進行求解，當然這只是一部分。其實精確算法也好，啟發式算法也好，都是獨立的算法，可以不依賴求解器進行代碼實現的，只要過程符合算法框架即可。
只不過平常看到的大部分是精確算法在各種整數規划模型上的應用，為此難免脫離不了cplex等求解器。這里簡單提一下。今天給大家帶來的依然是branch and bound算法在整數規划中的應用的代碼實現，所以還是會用到部分求解器的。
注：本文代碼下載請關注公眾號【程序猿聲】，后台回復【bbcode】，不包括【】即可下載！
01 Example-1
首先來看第一個代碼實例，該代碼求解的是整數優化的模型，關於branch and bound求解整數規划的具體原理就不再概述了，和上一篇文章差不多但是有所區別。代碼文件層次如下：
其中branch and bound算法主要部分在BnB_Guide.java這個文件。ExampleProblem.java內置了三個整數規划模型的實例。調用的是scpsolver這個求解器的wrapper，實際調用的還是lpsolver這個求解器用以求解線性松弛模型。下面着重講講BnB_Guide.java這個文件。
public BnB_Guide(int demoProblem){
example = new ExampleProblem(demoProblem);
LinearProgram lp = new LinearProgram();
lp = example.getProblem().getLP();
solver = SolverFactory.newDefault();
double[] solution = solver.solve(lp); // Solution of the initial relaxation problem
int maxElement = getMax(solution); // Index of the maximum non-integer decision variable's value
if(maxElement == -1 ) // We only got integers as values, hence we have an optimal solution
verifyOptimalSolution(solution,lp);
else
this.solveChildProblems(lp, solution, maxElement); // create 2 child problems and solve them
printSolution();
}
該過程是算法主調用過程：
- 首先變量lp保存了整數規划的松弛問題。
- 在調用求解器求解松弛模型以后，判斷是否所有決策變量都是整數了，如果是，已經找到最優解。
- 如果不是，根據找出最大的非整數的決策變量，對該變量進行分支，solveChildProblems。
接着是分支子問題的求解過程solveChildProblems如下：
public void solveChildProblems(LinearProgram lp, double[] solution ,int maxElement){
searchDepth++;
LinearProgram lp1 = new LinearProgram(lp);
LinearProgram lp2 = new LinearProgram(lp);
String constr_name = "c" + (lp.getConstraints().size() + 1); // Name of the new constraint
double[] constr_val = new double[lp.getDimension()]; // The variables' values of the new constraint
for(int i=0;i<constr_val.length;i++){ // Populate the table
if(i == maxElement )
constr_val[i] = 1.0;
else
constr_val[i] = 0;
}
//Create 2 child problems: 1. x >= ceil(value), 2. x <= floor(value)
lp1.addConstraint(new LinearBiggerThanEqualsConstraint(constr_val, Math.ceil(solution[maxElement]), constr_name));
lp2.addConstraint(new LinearSmallerThanEqualsConstraint(constr_val, Math.floor(solution[maxElement]), constr_name));
solveProblem(lp1);
solveProblem(lp2);
}
具體的分支過程如下：
- 首先新建兩個線性的子問題。
- 兩個子問題分別添加需要分支的決策變量新約束：1. x >= ceil(value), 2. x <= floor(value)。
- 一切准備就緒以后，調用solveProblem求解兩個子問題。
而solveProblem的實現代碼如下：
private void solveProblem(LinearProgram lp) {
double[] sol = solver.solve(lp);
LPSolution lpsol = new LPSolution(sol, lp);
double objVal = lpsol.getObjectiveValue();
if(lp.isMinProblem()) {
if(objVal > MinimizeProblemOptimalSolution) {
System.out.println("cut >>> objVal = "+ objVal);
return;
}
}
else {
if(objVal < MaximizeProblemOptimalSolution) {
System.out.println("cut >>> objVal = "+ objVal);
return;
}
}
System.out.println("non cut >>> objVal = "+ objVal);
int maxElement = this.getMax(sol);
if(maxElement == -1 && lp.isFeasable(sol)){ //We found a solution
solutionFound = true;
verifyOptimalSolution(sol,lp);
}
else if(lp.isFeasable(sol) && !solutionFound) //Search for a solution in the child problems
this.solveChildProblems(lp, sol, maxElement);
}
該過程如下：
- 首先調用求解器求解傳入的線性模型。
- 然后實行定界剪支，如果子問題的objVal比當前最優解還要差，則剪掉。
- 如果不剪，則判斷是否所有決策變量都是整數以及解是否可行，如果是，找到新的解，更新當前最優解。
- 如果不是，根據找出最大的非整數的決策變量，對該變量再次進行分支，進入solveChildProblems。
從上面的邏輯過程可以看出，solveChildProblems和solveProblem兩個之間相互調用，其實這是一種遞歸。該實現方式進行的就是BFS廣度優先搜索的方式遍歷搜索樹。
02 Example-2
再來看看第二個實例：
input是模型的輸入，輸入的是一個整數規划的模型。由於輸入和建模過程有點繁瑣，這里就不多講了。挑一些重點講講具體是分支定界算法是怎么運行的就行。
首先該代碼用了stack的作為數據結構，遍歷搜索樹的方式是DFS即深度優先搜索，我們來看BNBSearch.java這個文件：
public class BNBSearch {
Deque<searchNode> searchStack = new ArrayDeque<searchNode>();
double bestVal = Double.MAX_VALUE;
searchNode currentBest = new searchNode();
IPInstance solveRel = new IPInstance();
Deque<searchNode> visited = new ArrayDeque<searchNode>();
public BNBSearch(IPInstance solveRel) {
this.solveRel = solveRel;
searchNode rootNode = new searchNode();
this.searchStack.push(rootNode);
};
BNBSearch 這個類是branch and bound算法的主要過程，成員變量如下：
- searchStack ：構造和遍歷生成樹用的，棧結構。
- bestVal：記錄當前最優解的值，由於求的最小化問題，一開始設置為正無窮。
- currentBest ：記錄當前最優解。
- solveRel ：整數規划模型。
- visited ：記錄此前走過的分支，避免重復。
然后在這里展開講一下searchNode就是構成搜索樹的節點是怎么定義的：
public class searchNode {
HashMap<Integer, Integer> partialAssigned = new HashMap<Integer, Integer>();
public searchNode() {
super();
}
public searchNode(searchNode makeCopy) {
for (int test: makeCopy.partialAssigned.keySet()) {
this.partialAssigned.put(test, makeCopy.partialAssigned.get(test));
}
}
}
其實非常簡單，partialAssigned 保存的是部分解的結構，就是一個HashMap，key保存的是決策變量，而value對應的是決策變量分支的取值（0-1）。舉上節課講過的例子：
比如：
- 節點1的partialAssigned == { {x3, 1} }。
- 節點2的partialAssigned == { {x3, 0} }。
- 節點3的partialAssigned == { {x3, 1}, {x2, 1} }。
- 節點4的partialAssigned == { {x3, 1}, {x2, 0} }。
- 節點7的partialAssigned == { {x3, 0}, {x1, 1}, {x2, 1}}。
……
想必各位已經明白得不能再明白了。
然后就可以開始BB過程了：
public int solveIP() throws IloException {
while (!this.searchStack.isEmpty()) {
searchNode branchNode = this.searchStack.pop();
boolean isVisited = false;
for (searchNode tempNode: this.visited) {
if (branchNode.partialAssigned.equals(tempNode.partialAssigned)){
isVisited = true;
break;
}
}
if (!isVisited) {
visited.add(new searchNode(branchNode));
double bound = solveRel.solve(branchNode);
if (bound > bestVal || bound == 0) {
//System.out.println(searchStack.size());
}
if (bound < bestVal && bound!=0) {
if (branchNode.partialAssigned.size() == solveRel.numTests) {
//分支到達低端，找到一個滿足整數約束的可行解，設置為當前最優解。
//System.out.println("YAY");
this.bestVal = bound;
this.currentBest = branchNode;
}
}
if (bound < bestVal && bound!=0) {
//如果還沒到達低端，找一個變量進行分支。
if (branchNode.partialAssigned.size() != solveRel.numTests) {
int varToSplit = getSplitVariable(branchNode);
if (varToSplit != -1) {
searchNode left = new searchNode(branchNode);
searchNode right = new searchNode(branchNode);
left.partialAssigned.put(varToSplit, 0);
right.partialAssigned.put(varToSplit, 1);
this.searchStack.push(left);
this.searchStack.push(right);
}
}
}
}
}
return (int) bestVal;
}
首先從搜索棧里面取出一個節點，判斷節點代表的分支是否此前已經走過了，重復的工作就不要做了嘛。
如果沒有走過，那么在該節點處進行定界操作，從該節點進入，根據partialAssigned 保存的部分解結構，添加約束，建立松弛模型，調用cplex求解。具體求解過程如下：
public double solve(searchNode node) throws IloException {
try {
cplex = new IloCplex();
cplex.setOut(null);
IloNumVarType [] switcher = new IloNumVarType[2];
switcher[0] = IloNumVarType.Int;
switcher[1] = IloNumVarType.Float;
int flag = 1;
IloNumVar[] testUsed = cplex.numVarArray(numTests, 0, 1, switcher[flag]);
IloNumExpr objectiveFunction = cplex.numExpr();
objectiveFunction = cplex.scalProd(testUsed, costOfTest);
cplex.addMinimize(objectiveFunction);
for (int j = 0; j < numDiseases*numDiseases; j++) {
if (j % numDiseases == j /numDiseases) {
continue;
}
IloNumExpr diffConstraint = cplex.numExpr();
for (int i = 0; i < numTests; i++) {
if (A[i][j/numDiseases] == A[i][j%numDiseases]) {
continue;
}
diffConstraint = cplex.sum(diffConstraint, testUsed[i]);
}
cplex.addGe(diffConstraint, 1);
diffConstraint = cplex.numExpr();
}
for (int test: node.partialAssigned.keySet()) {
cplex.addEq(testUsed[test], node.partialAssigned.get(test));
}
//System.out.println(cplex.getModel());
if(cplex.solve()) {
double objectiveValue = (cplex.getObjValue());
for (int i = 0; i < numTests; i ++) {
if (cplex.getValue(testUsed[i]) == 0) {
node.partialAssigned.put(i, 0);
}
else if (cplex.getValue(testUsed[i]) == 1) {
node.partialAssigned.put(i, 1);
}
}
//System.out.println("LOL"+node.partialAssigned.size());
return objectiveValue;
}
}
catch(IloException e) {
System.out.println("Error " + e);
}
return 0;
}
中間一大堆建模過程就不多講了，具體分支約束是這一句：
for (int test: node.partialAssigned.keySet()) {
cplex.addEq(testUsed[test], node.partialAssigned.get(test));
}
此后，求解完畢后，把得到整數解的決策變量放進partialAssigned，不是整數后續操作。然后返回目標值。
然后依舊回到solveIP里面，在進行求解以后，得到目標值，接下來就是定界操作了：
- if (bound > bestVal || bound == 0)：剪支。
- if (bound < bestVal && bound!=0)：判斷是否所有決策變量都為整數，如果是，找到一個可行解，更新當前最優解。如果不是，找一個小數的決策變量入棧，等待后續分支。
03 運行說明
Example-1：
運行說明，運行輸入參數1到3中的數字表示各個不同的模型，需要在32位JDK環境下才能運行，不然會報nullPointer的錯誤，這是那份求解器wrapper的鍋。怎么設置參數參考cplexTSP那篇，怎么設置JDK環境就不多說了。
然后需要把代碼文件夾下的幾個jar包給添加進去，再把lpsolve的dll給放到native library里面，具體做法還是參照cplexTSP那篇，重復的內容我就不多說了。
Example-2：
最后是運行說明：該實例運行調用了cplex求解器，所以需要配置cplex環境才能運行，具體怎么配置看之前的教程。JDK環境要求64位，無參數輸入。
代碼來源GitHub，經過部分修改。