前言
00
前面一篇文章我们讲了branch and bound算法的相关概念。可能大家对精确算法实现的印象大概只有一个，调用求解器进行求解，当然这只是一部分。
其实精确算法也好，启发式算法也好，都是独立的算法，可以不依赖求解器进行代码实现的，只要过程符合算法框架即可。
只不过平常看到的大部分是精确算法在各种整数规划模型上的应用，为此难免脱离不了cplex等求解器。这里简单提一下。
今天给大家带来的依然是branch and bound算法在整数规划中的应用的代码实现，所以还是会用到部分求解器的。
注：本文代码下载请移步留言区。
Example-1
01
首先来看第一个代码实例，该代码求解的是整数优化的模型，关于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(); }
该过程是算法主调用过程：
1. 首先变量lp保存了整数规划的松弛问题。
2. 在调用求解器求解松弛模型以后，判断是否所有决策变量都是整数了，如果是，已经找到最优解。
3. 如果不是，根据找出最大的非整数的决策变量，对该变量进行分支，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. 首先新建两个线性的子问题。
2. 两个子问题分别添加需要分支的决策变量新约束：1. x >= ceil(value), 2. x <= floor(value)。
3. 一切准备就绪以后，调用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); }
该过程如下：
1. 首先调用求解器求解传入的线性模型。
2. 然后实行定界剪支，如果子问题的objVal比当前最优解还要差，则剪掉。
3. 如果不剪，则判断是否所有决策变量都是整数以及解是否可行，如果是，找到新的解，更新当前最优解。
4. 如果不是，根据找出最大的非整数的决策变量，对该变量再次进行分支，进入solveChildProblems。
从上面的逻辑过程可以看出，solveChildProblems和solveProblem两个之间相互调用，其实这是一种递归。
该实现方式进行的就是BFS广度优先搜索的方式遍历搜索树。
Example-2
02
再来看看第二个实例：
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; }
中间一大堆建模过程就不多讲了，具体分支约束是这一句：
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; }
此后，求解完毕后，把得到整数解的决策变量放进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，经过部分修改。
