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

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();
		
	}

该过程是算法主调用过程:

  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广度优先搜索的方式遍历搜索树。

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,经过部分修改。

posted @ 2019-07-23 12:57  短短的路走走停停  阅读(4753)  评论(0编辑  收藏  举报