2. it技术
  3. 分支界定算法实践(一)求解线性整数优化问题

分支界定算法实践(一)求解线性整数优化问题

it2024-11-04  7

分支界定算法的主要思想是把问题的可行解划分为多个子集，对子集进行系统化的评估直到最优解找到。分支界定只是一种算法框架，因为如何构造子集、在每个子集如何搜索最优解可以按不同情况有不同的具体实现，当用于求解线性整数优化问题时，就需要结合一般的线性优化算法或者求解器来实现。下面先用一个实例来说明算法的思路，然后用代码实现一下。

一家工厂打算购买两种设备，印刷机和车床。老板估计每台印刷机每台可以增加$100的利润，每台车床每天可以增加$150的利润；但同时购买数量也会收到预算和场地面积的约束，每台印刷机或车床的花费和占地如下表所示：

MachineRequired Floor SpacePurchase Price印刷机15$8000车床30$4000

老板的预算是$40000，并且现在有200平方的空地；工厂老板应该怎样购买来让每天的利润最大？ 如果我们设 x 1 x_1 x1​是购买多少台印刷机， x 2 x_2 x2​是购买多少台车床，那么可以写出优化问题的模型 m a x i m i e Z = 100 x 1 + 150 x 2 s u b j e c t   t o 8000 x 1 + 4000 x 2 ≤ 40000 15 x 1 + 30 x 2 ≤ 200 x 1 , x 2 ≥ 0   a n d   i n t e g e r \begin{aligned} maximie\quad& Z=100x_1+150x_2\\ subject\ to\quad& 8000x_1+4000x_2\leq 40000 \\ & 15x_1+30x_2\leq 200 \\ & x_1,x_2\geq 0\ and\ integer \end{aligned} maximiesubject to​Z=100x1​+150x2​8000x1​+4000x2​≤4000015x1​+30x2​≤200x1​,x2​≥0 and integer​ 相比线性优化模型，上面的数学模型就多出必须为整数的约束。如果我们把整数约束去掉，那么可以得到最优解: x 1 = 2.22 , x 2 = 5.56 ,   a n d   Z = 1055.56 x_1=2.22,x_2=5.56,\ and\ Z=1055.56 x1​=2.22,x2​=5.56, and Z=1055.56 如果我们把这个解四舍五入，可以得到一个可行的整数解 x 1 = 2 , x 2 = 5 ,   a n d   Z = 950 x_1=2,x_2=5,\ and\ Z=950 x1​=2,x2​=5, and Z=950 但是这个整数解并不是实际的最优解，因为一旦加上整数约束，线性优化模型往往是非凸的，线性优化算法往往无能为力。而分支界定法则是把原问题进行松弛(即去掉整数约束)，以松弛线性模型为根节点不断地划分子问题和更新问题上下界，通过评判上下界来忽略无效子问题和停止算法。首先我们构造根节点，该节点就是上面模型去掉整数约束后的松弛模型，对于这个松弛线性最大化问题，我们构造一下上下界，用它的最优解( Z = 1055.56 Z=1055.56 Z=1055.56)作为问题的上界，用四舍五入后的整数解( x 1 = 2 , x 2 = 5 , Z = 950 x_1=2,x_2=5,Z=950 x1​=2,x2​=5,Z=950)作为下界；上下界的意义在于，原始整数优化模型的解一定是在上下界之间的，所以对于所有最优解比下界更小的子问题，我们可以直接忽略；而子问题的上界一定不会大于父问题，通过不断分支，这个上下界间的范围将不断缩减。 因为根节点的最优解不是整数解，所以我们选择一个非整数解进行分支。这里我们采取选择离四舍五入值最远的变量，对于根节点，即选择了变量 x 2 = 5.56 x_2=5.56 x2​=5.56。我们以变量 x 2 x_2 x2​为划分点向下分支出两个子问题，分别在根节点的松弛优化模型的基础上新增约束 x 2 ≤ 5 x_2\leq5 x2​≤5和 x 2 ≥ 6 x_2\geq6 x2​≥6，也就是说，对于左侧的子节点2，松弛优化模型为： m a x i m i e Z = 100 x 1 + 150 x 2 s u b j e c t   t o 8000 x 1 + 4000 x 2 ≤ 40000 15 x 1 + 30 x 2 ≤ 200 x 1 , x 2 ≥ 0 x 2 ≤ 5 \begin{aligned} maximie\quad& Z=100x_1+150x_2\\ subject\ to\quad& 8000x_1+4000x_2\leq 40000 \\ & 15x_1+30x_2\leq 200 \\ & x_1,x_2\geq 0\\ & x_2 \leq5 \end{aligned} maximiesubject to​Z=100x1​+150x2​8000x1​+4000x2​≤4000015x1​+30x2​≤200x1​,x2​≥0x2​≤5​ 右侧子节点3的问题模型为: m a x i m i e Z = 100 x 1 + 150 x 2 s u b j e c t   t o 8000 x 1 + 4000 x 2 ≤ 40000 15 x 1 + 30 x 2 ≤ 200 x 1 , x 2 ≥ 0 x 2 ≥ 6 \begin{aligned} maximie\quad& Z=100x_1+150x_2\\ subject\ to\quad& 8000x_1+4000x_2\leq 40000 \\ & 15x_1+30x_2\leq 200 \\ & x_1,x_2\geq 0\\ & x_2 \geq6 \end{aligned} maximiesubject to​Z=100x1​+150x2​8000x1​+4000x2​≤4000015x1​+30x2​≤200x1​,x2​≥0x2​≥6​ 对于节点2，可以求得最优解为 x 1 = 2.5 , x 2 = 5 , Z = 1000 x_1=2.5,x_2=5,Z=1000 x1​=2.5,x2​=5,Z=1000，而节点3的最优解为 x 1 = 1.33 , x 2 = 6 , Z = 1033.33 x_1=1.33,x_2=6,Z=1033.33 x1​=1.33,x2​=6,Z=1033.33，这两个解可分别作为这两个节点的上界。可以看到这两个上界都小于根节点的上界，因为每一次分支，都是在原松弛优化模型的基础上新增约束，它的最优解肯定不会好过原问题。再来分析下界，因为这两个子节点的解都不是完整整数解，所以这两个节点下界保持不变，仍是950

我们还没有找到最优可行解，所以要从结点2或者3继续分支。选择哪个结点呢，这里我们按上界更大优先的原则来选取结点，也就是从结点3进行分支。采用这种原则其实是一种间接地剪枝，如果某个结点的最优解是非整数解并且小于当前的全局下界(即上界比下界更小)，那么需要把这个结点及其所有子节点都从搜索数里去除，它们的解肯定不会是最优解了，通过上界越大越优先的原则来选取，这些结点实际上根本不会被选中。结点3的 x 2 x_2 x2​已经是整数了，所以用 x 1 x_1 x1​进行划分： x 1 ≤ 1 x_1\leq1 x1​≤1和 x 1 ≥ 2 x_1\geq2 x1​≥2，分出了结点4和结点5，结点4的问题模型是在3的基础上加上 x 1 ≤ 1 x_1\leq1 x1​≤1 m a x i m i e Z = 100 x 1 + 150 x 2 s u b j e c t   t o 8000 x 1 + 4000 x 2 ≤ 40000 15 x 1 + 30 x 2 ≤ 200 x 1 , x 2 ≥ 0 x 2 ≥ 6 x 1 ≤ 1 \begin{aligned} maximie\quad& Z=100x_1+150x_2\\ subject\ to\quad& 8000x_1+4000x_2\leq 40000 \\ & 15x_1+30x_2\leq 200 \\ & x_1,x_2\geq 0\\ & x_2 \geq6 \\ & x_1 \leq1 \end{aligned} maximiesubject to​Z=100x1​+150x2​8000x1​+4000x2​≤4000015x1​+30x2​≤200x1​,x2​≥0x2​≥6x1​≤1​ 结点5则是在3的基础上加上 x 1 ≥ 2 x_1\geq2 x1​≥2 m a x i m i e Z = 100 x 1 + 150 x 2 s u b j e c t   t o 8000 x 1 + 4000 x 2 ≤ 40000 15 x 1 + 30 x 2 ≤ 200 x 1 , x 2 ≥ 0 x 2 ≥ 6 x 1 ≥ 2 \begin{aligned} maximie\quad& Z=100x_1+150x_2\\ subject\ to\quad& 8000x_1+4000x_2\leq 40000 \\ & 15x_1+30x_2\leq 200 \\ & x_1,x_2\geq 0\\ & x_2 \geq6 \\ & x_1 \geq2 \end{aligned} maximiesubject to​Z=100x1​+150x2​8000x1​+4000x2​≤4000015x1​+30x2​≤200x1​,x2​≥0x2​≥6x1​≥2​ 对于结点4，其最优解是( x 1 = 1 , x 2 = 6.17 , Z = 1025 x_1=1,x_2=6.17,Z=1025 x1​=1,x2​=6.17,Z=1025)；对于结点5，则无解。那就只需要对结点4进行定界，上界就是最优解1025.5，同时最优解也不是完整整数解，那下界仍保持不变。

我们还是没有找到最优整数解，需要继续分支。现在有3个叶节点：2、4和5，结点5已经无解，所以不需要考虑，而在2和4中，4的上界更大，所以从结点4开始分支。在结点4， x 1 x_1 x1​已经是整数，所以用 x 2 x_2 x2​进行划分： x 2 ≤ 6 x_2\leq6 x2​≤6和 x 2 ≥ 7 x_2\geq7 x2​≥7。结点6的问题模型如下，其最优解为( x 1 = 1 , x 2 = 6 , Z = 1000 x_1=1,x_2=6,Z=1000 x1​=1,x2​=6,Z=1000) m a x i m i e Z = 100 x 1 + 150 x 2 s u b j e c t   t o 8000 x 1 + 4000 x 2 ≤ 40000 15 x 1 + 30 x 2 ≤ 200 x 1 , x 2 ≥ 0 x 2 ≥ 6 x 1 ≥ 2 x 2 ≤ 6 \begin{aligned} maximie\quad& Z=100x_1+150x_2\\ subject\ to\quad& 8000x_1+4000x_2\leq 40000 \\ & 15x_1+30x_2\leq 200 \\ & x_1,x_2\geq 0\\ & x_2 \geq6 \\ & x_1 \geq2 \\ & x_2\leq6 \end{aligned} maximiesubject to​Z=100x1​+150x2​8000x1​+4000x2​≤4000015x1​+30x2​≤200x1​,x2​≥0x2​≥6x1​≥2x2​≤6​ 结点7的问题模型如下，它同样无解 m a x i m i e Z = 100 x 1 + 150 x 2 s u b j e c t   t o 8000 x 1 + 4000 x 2 ≤ 40000 15 x 1 + 30 x 2 ≤ 200 x 1 , x 2 ≥ 0 x 2 ≥ 6 x 1 ≥ 2 x 2 ≥ 7 \begin{aligned} maximie\quad& Z=100x_1+150x_2\\ subject\ to\quad& 8000x_1+4000x_2\leq 40000 \\ & 15x_1+30x_2\leq 200 \\ & x_1,x_2\geq 0\\ & x_2 \geq6 \\ & x_1 \geq2 \\ & x_2\geq7 \end{aligned} maximiesubject to​Z=100x1​+150x2​8000x1​+4000x2​≤4000015x1​+30x2​≤200x1​,x2​≥0x2​≥6x1​≥2x2​≥7​ 现在搜索树如下图所示

在结点6，因为出现了新的整数解，并且该整数解的目标值大于之前的下界950，所以结点6的下界更新为1000，这时结点6的上下界已经相等。现在可比较叶节点是2和6(排除掉无解的5和7)，可以发现结点6是完整的整数解并且上界是叶节点中最大的，那就说明结点6的解就是原始整数规划的最优解。为什么呢，结点6的上界比其他叶节点的上界大，而从叶节点继续向下分支的上界总是比叶节点的小，所以结点6的上界肯定是比其他所有子节点的上界更大。由此我们找到了原问题的解，即购买1台印刷机和6台车床可以达到最大日利润$1000

我们来归纳一下分支界定算法的一般步骤。我注意到不同的文章里对算法步骤的描述都有所区别，不过核心的几个要素是一致的。

Step 1: 计算原始整数规划松弛化后的最优值，如果无解，则原问题就是无解，退出算法；如果最优解就是整数解，那么直接得到最终解，退出算法Step 2: 建立根节点，在根节点上，上界是原始整数规划松弛后的最优解，下界则需要预先给出一个可行解来得到Step 3: 评估当前所有的叶节点，不过不需要加入无解的结点。如果叶节点某个结点的最优解是完整整数解，并且其上界比其他叶节点都要大，则该结点的解就是全局最优解，可以退出算法；否则继续Step 4；另外如果目前没有任何可行的叶节点，那么就表明除了初始给出的可行解，不存在任何其他可行解，也直接退出算法Step 4: 从最大上界的可行叶节点继续分支。选择该结点的一个非整数解，产生两个新的 ≤ \leq ≤和 ≥ \geq ≥约束，分别加到该节点的问题模型上产生两个子节点。计算更新两个子节点的上下界，上界值等于子节点的最优解，下界等于全局的最大下界。Step 5: 回到Step 3

最后我用C#代码实现一下上面问题求解。首先需要说明的是，求解松弛问题，我使用了Google OR-Tools的线性求解器。定义了一个LooseLinearConstraintProblem类，来表示基于上面例子的松弛优化问题模型，在内部定义了最基本的约束，并且定义了additionalConstraints来表示相对父问题的额外约束，如果要从一个父松弛问题生成一个子松弛问题，可以通过拷贝构造函数来生成

public class LooseLinearConstraintProblem { private Solver solver; private List<Variable> variables; public List<Variable> Variables { get { return variables; } } //private List<Constraint> constraints; private Objective objective; //item1: 第几个变量 //item2: 大于等于多少 //item3: 小于等于多少 private List<Tuple<int, double, double>> additionalConstraints; public List<Tuple<int, double, double>> AdditionalConstraints { set { this.additionalConstraints = value; } get { return this.additionalConstraints; } } private void BasicInit() { this.solver = Solver.CreateSolver("GLOP"); this.variables = new List<Variable>(); Variable x1 = solver.MakeNumVar(0.0, double.PositiveInfinity, "x1"); Variable x2 = solver.MakeNumVar(0.0, double.PositiveInfinity, "x2"); this.variables.Add(x1); this.variables.Add(x2); Constraint c0 = solver.MakeConstraint(double.NegativeInfinity, 40000); c0.SetCoefficient(x1, 8000); c0.SetCoefficient(x2, 4000); Constraint c1 = solver.MakeConstraint(double.NegativeInfinity, 200); c1.SetCoefficient(x1, 15); c1.SetCoefficient(x2, 30); this.objective = solver.Objective(); this.objective.SetCoefficient(x1, 100); this.objective.SetCoefficient(x2, 150); this.objective.SetMaximization(); } public LooseLinearConstraintProblem(List<Tuple<int, double, double>> additionalConstraints) { BasicInit(); if (additionalConstraints != null) { this.additionalConstraints = additionalConstraints; foreach(var additionalConstraint in additionalConstraints) { Constraint newConstraint = this.solver.MakeConstraint(additionalConstraint.Item2, additionalConstraint.Item3); for (int j = 0; j < this.variables.Count; j++) { if (j != additionalConstraint.Item1) { newConstraint.SetCoefficient(this.variables[j], 0); } else { newConstraint.SetCoefficient(this.variables[j], 1); } } } } } public LooseLinearConstraintProblem(LooseLinearConstraintProblem parentProblem,List<Tuple<int, double, double>> additionalConstraints) { BasicInit(); this.additionalConstraints = new List<Tuple<int, double, double>>(); if (parentProblem.AdditionalConstraints != null) { foreach (var additionalConstraint in parentProblem.AdditionalConstraints) { this.additionalConstraints.Add(additionalConstraint); Constraint newConstraint = this.solver.MakeConstraint(additionalConstraint.Item2, additionalConstraint.Item3); for (int j = 0; j < this.variables.Count; j++) { if (j != additionalConstraint.Item1) { newConstraint.SetCoefficient(this.variables[j], 0); } else { newConstraint.SetCoefficient(this.variables[j], 1); } } } } if (additionalConstraints != null) { foreach (var additionalConstraint in additionalConstraints) { this.additionalConstraints.Add(additionalConstraint); Constraint newConstraint = this.solver.MakeConstraint(additionalConstraint.Item2, additionalConstraint.Item3); for (int j = 0; j < this.variables.Count; j++) { if (j != additionalConstraint.Item1) { newConstraint.SetCoefficient(this.variables[j], 0); } else { newConstraint.SetCoefficient(this.variables[j], 1); } } } } } public Tuple<List<double>, double> Solve() { var status = solver.Solve(); if (status != Solver.ResultStatus.OPTIMAL) { return null; } else { List<double> vResults = new List<double>(); foreach(var v in this.variables) { vResults.Add(v.SolutionValue()); } double oResult = solver.Objective().Value(); return new Tuple<List<double>, double>(vResults, oResult); } } }

顶一个Node类来表示分支界定中的每个结点

//每个问题节点包括: 松弛线性优化子问题模型, 当前节点的问题上界, 当前节点的问题下界, 当前问题的求解情况(无完整整数解, 完整整数解, 无解) public class Node { private LooseLinearConstraintProblem problem; public LooseLinearConstraintProblem Problem { get { return this.problem; } } public double UpperBound { get;set; } public double LowerBound { get; set; } private SolveStatus solveStatus; public SolveStatus SolveStatus { get { return this.solveStatus; } } private Tuple<List<double>, double> solution; public Node(LooseLinearConstraintProblem problem) { this.problem = problem; this.solveStatus = SolveStatus.HaveNotSolve; } public Tuple<List<double>, double> Solve() { if (this.solveStatus != SolveStatus.HaveNotSolve) { return this.solution; } Tuple<List<double>, double> result = this.problem.Solve(); if (result == null) { this.solveStatus = SolveStatus.Infeasible; } else if (getIntegerIndex(result.Item1).Count != result.Item1.Count) { this.solveStatus = SolveStatus.NoneCompleteIntegerSolution; } else { this.solveStatus = SolveStatus.CompleteIntegerSolution; } this.solution = result; return this.solution; } }

Node类用到自定义的结构体SolveStatus来表示当前结点的求解状态，

public enum SolveStatus { HaveNotSolve, //还没求解该结点 NoneCompleteIntegerSolution, //非完整整数解 CompleteIntegerSolution, //完整整数解 Infeasible //无解 }

两个辅助方法，一个来判断一个double变量是否是整数，一个找出一个double变量序列中的所有非整数变量

//从一个double数列中找出所有为整数的元素的下标 public static List<int> getIntegerIndex(List<double> values) { List<int> results = new List<int>(); for(int i = 0; i < values.Count; i++) { if (isInteger(values[i])) { results.Add(i); } } return results; } public static bool isInteger(double n) { string s = Convert.ToString(n); return s.Contains(".") ? false : true; }

最后是算法的主流程

public void Test(List<double> initVars, double initObj) { //初始的松弛线性优化问题 LooseLinearConstraintProblem baseProblem = new LooseLinearConstraintProblem(null); var baseProblemResult = baseProblem.Solve(); //如果root问题就无解，那么原问题就是无解 if (baseProblemResult == null) { throw new Exception("No solution for this problem"); } List<int> integerIndexList = getIntegerIndex(baseProblemResult.Item1); //如果root问题的解全是整数，那么已经找到原问题的最优解 if (integerIndexList.Count == baseProblemResult.Item1.Count) { Console.WriteLine("Solution:"); for(int i = 0; i < baseProblem.Variables.Count; i++) { Console.WriteLine(" x"+i+" = " + baseProblemResult.Item1[i]); } Console.WriteLine("Optimal objective value = " + baseProblemResult.Item2); return; } //开始进行分支界定 double globalUpperBound = baseProblemResult.Item2;//初始情况下问题的上界就是原始松弛问题的最优解 double globalLowerBound = initObj;//初始的问题下界手动给定, x1=2,x2=5 List<Node> endNodes = new List<Node>();//缓存所有叶节点 Node rootNode = new Node(baseProblem); rootNode.Solve(); rootNode.UpperBound = globalUpperBound; rootNode.LowerBound = globalLowerBound; endNodes.Add(rootNode); Console.WriteLine("Start Branch and Bound: "); Console.WriteLine($"Root Node, UpperBound={rootNode.UpperBound}, LowerBound={rootNode.LowerBound}"); while (true) { //如果叶节点序列为空, 那就表明除了我们初始给定的整数解没有其他可行解了 if (endNodes.Count == 0) { Console.WriteLine("Solution:"); for (int i = 0; i < baseProblem.Variables.Count; i++) { Console.WriteLine(" x" + i + " = " + baseProblemResult.Item1[i]); } Console.WriteLine("Optimal objective value = " + baseProblemResult.Item2); break; } //搜索的停止条件是存在一个完整整数解, 并且该整数解的目标值是所有叶节点(非无解节点)中最大的. //因为子节点的目标值肯定是比父节点小的(额外的约束导致), 所以如果一个整数解的目标值比其他所有 //叶节点大, 那么这些叶节点的子节点不可能比这个解更好 double maxUpperBound = double.NegativeInfinity; int maxUpperBoundNodeIndex = 0; for(int i = 0; i < endNodes.Count; i++) { if (endNodes[i].UpperBound >= maxUpperBound) { maxUpperBound = endNodes[i].UpperBound; maxUpperBoundNodeIndex = i; } } if (endNodes[maxUpperBoundNodeIndex].SolveStatus == SolveStatus.CompleteIntegerSolution) { Console.WriteLine("Solution:"); for (int i = 0; i < endNodes[maxUpperBoundNodeIndex].Problem.Variables.Count; i++) { Console.WriteLine(" x" + i + " = " + endNodes[maxUpperBoundNodeIndex].Solve().Item1[i]); } Console.WriteLine("Optimal objective value = " + endNodes[maxUpperBoundNodeIndex].Solve().Item2); break; } //用最大上边界的节点来分支, 选择一个非整数解新增约束 Node parentNode = endNodes[maxUpperBoundNodeIndex]; endNodes.RemoveAt(maxUpperBoundNodeIndex); Tuple<List<double>, double> parentSolution = parentNode.Solve(); //选择一个离四舍五入整数解最远的非整数解来划分 int selectVariableIndex = 0; double distance = 0; for(int i = 0; i < parentSolution.Item1.Count; i++) { if (!isInteger(parentSolution.Item1[i]) && Math.Abs(Math.Round(parentSolution.Item1[i])-parentSolution.Item1[i])>distance) { selectVariableIndex = i; distance = Math.Abs(Math.Round(parentSolution.Item1[i]) - parentSolution.Item1[i]); } } double leftVar = Math.Floor(parentSolution.Item1[selectVariableIndex]); double rightVar = Math.Ceiling(parentSolution.Item1[selectVariableIndex]); //新分出两个节点 LooseLinearConstraintProblem leftSubProblem = new LooseLinearConstraintProblem(parentNode.Problem, new List<Tuple<int, double, double>> { new Tuple<int, double, double>(selectVariableIndex, double.NegativeInfinity, leftVar) }); Node leftNode = new Node(leftSubProblem); LooseLinearConstraintProblem rightSubProblem = new LooseLinearConstraintProblem(parentNode.Problem, new List<Tuple<int, double, double>> { new Tuple<int, double, double>(selectVariableIndex, rightVar,double.PositiveInfinity) }); Node rightNode = new Node(rightSubProblem); //分别计算两个子节点的结果 leftNode.Solve(); rightNode.Solve(); //更新当前的叶节点序列, 无解的节点不需要加入 if (leftNode.SolveStatus != SolveStatus.Infeasible) { endNodes.Add(leftNode); } if (rightNode.SolveStatus != SolveStatus.Infeasible) { endNodes.Add(rightNode); } //更新当前所有叶节点的上下界 List<int> needRemoveIndexList = new List<int>(); for(int i=0;i<endNodes.Count;i++) { var node = endNodes[i]; node.UpperBound = node.Solve().Item2; //如果该节点没有完整整数解, 那么它的下边界就是全局下界 if (node.SolveStatus == SolveStatus.NoneCompleteIntegerSolution) { node.LowerBound = globalLowerBound; } // 如果是一个完整整数解, 那么该节点的上下边界相等, 进一步, 如果这个目标值大于当前的 // 全局下界, 则更新全局下界 else { node.LowerBound = node.UpperBound; if (node.LowerBound >= globalLowerBound) { globalLowerBound = node.LowerBound; } } //如果在该节点下界大于上界，把这个结点从叶节点序列中去除，即剪枝 if (node.LowerBound > node.UpperBound) { needRemoveIndexList.Add(i); } } foreach(var needRemove in needRemoveIndexList) { endNodes.RemoveAt(needRemove); } Console.WriteLine($"LeftNode, additional constraint: x{selectVariableIndex + 1}<={leftVar}"); if (leftNode.SolveStatus == SolveStatus.Infeasible) { Console.WriteLine($"LeftNode is infeasible"); } else { Console.WriteLine($"LeftNode, UpperBound={leftNode.UpperBound}, LowerBound={leftNode.LowerBound}"); } Console.WriteLine($"RightNode, additional constraint: {rightVar}<=x{selectVariableIndex + 1}"); if (rightNode.SolveStatus == SolveStatus.Infeasible) { Console.WriteLine($"RightNode is infeasible"); } else { Console.WriteLine($"RightNode, UpperBound={rightNode.UpperBound}, LowerBound={rightNode.LowerBound}"); } } }

完整的代码:

using Google.OrTools.LinearSolver; using System; using System.Collections.Generic; using System.Text; namespace BranchAndBound { public class Simple { public enum SolveStatus { HaveNotSolve, NoneCompleteIntegerSolution, CompleteIntegerSolution, Infeasible } //每个问题节点包括: 松弛线性优化子问题模型, 当前节点的问题上界, 当前节点的问题下界, 当前问题的求解情况(无完整整数解, 完整整数解, 无解) public class Node { private LooseLinearConstraintProblem problem; public LooseLinearConstraintProblem Problem { get { return this.problem; } } public double UpperBound { get;set; } public double LowerBound { get; set; } private SolveStatus solveStatus; public SolveStatus SolveStatus { get { return this.solveStatus; } } private Tuple<List<double>, double> solution; public Node(LooseLinearConstraintProblem problem) { this.problem = problem; this.solveStatus = SolveStatus.HaveNotSolve; } public Tuple<List<double>, double> Solve() { if (this.solveStatus != SolveStatus.HaveNotSolve) { return this.solution; } Tuple<List<double>, double> result = this.problem.Solve(); if (result == null) { this.solveStatus = SolveStatus.Infeasible; } else if (getIntegerIndex(result.Item1).Count != result.Item1.Count) { this.solveStatus = SolveStatus.NoneCompleteIntegerSolution; } else { this.solveStatus = SolveStatus.CompleteIntegerSolution; } this.solution = result; return this.solution; } } public void Test(List<double> initVars, double initObj) { //初始的松弛线性优化问题 LooseLinearConstraintProblem baseProblem = new LooseLinearConstraintProblem(null); var baseProblemResult = baseProblem.Solve(); //如果root问题就无解，那么原问题就是无解 if (baseProblemResult == null) { throw new Exception("No solution for this problem"); } List<int> integerIndexList = getIntegerIndex(baseProblemResult.Item1); //如果root问题的解全是整数，那么已经找到原问题的最优解 if (integerIndexList.Count == baseProblemResult.Item1.Count) { Console.WriteLine("Solution:"); for(int i = 0; i < baseProblem.Variables.Count; i++) { Console.WriteLine(" x"+i+" = " + baseProblemResult.Item1[i]); } Console.WriteLine("Optimal objective value = " + baseProblemResult.Item2); return; } //开始进行分支界定 double globalUpperBound = baseProblemResult.Item2;//初始情况下问题的上界就是原始松弛问题的最优解 double globalLowerBound = initObj;//初始的问题下界手动给定, x1=2,x2=5 List<Node> endNodes = new List<Node>();//缓存所有叶节点 Node rootNode = new Node(baseProblem); rootNode.Solve(); rootNode.UpperBound = globalUpperBound; rootNode.LowerBound = globalLowerBound; endNodes.Add(rootNode); Console.WriteLine("Start Branch and Bound: "); Console.WriteLine($"Root Node, UpperBound={rootNode.UpperBound}, LowerBound={rootNode.LowerBound}"); while (true) { //如果叶节点序列为空, 那就表明除了我们初始给定的整数解没有其他可行解了 if (endNodes.Count == 0) { Console.WriteLine("Solution:"); for (int i = 0; i < baseProblem.Variables.Count; i++) { Console.WriteLine(" x" + i + " = " + baseProblemResult.Item1[i]); } Console.WriteLine("Optimal objective value = " + baseProblemResult.Item2); break; } //搜索的停止条件是存在一个完整整数解, 并且该整数解的目标值是所有叶节点(非无解节点)中最大的. //因为子节点的目标值肯定是比父节点小的(额外的约束导致), 所以如果一个整数解的目标值比其他所有 //叶节点大, 那么这些叶节点的子节点不可能比这个解更好 double maxUpperBound = double.NegativeInfinity; int maxUpperBoundNodeIndex = 0; for(int i = 0; i < endNodes.Count; i++) { if (endNodes[i].UpperBound >= maxUpperBound) { maxUpperBound = endNodes[i].UpperBound; maxUpperBoundNodeIndex = i; } } if (endNodes[maxUpperBoundNodeIndex].SolveStatus == SolveStatus.CompleteIntegerSolution) { Console.WriteLine("Solution:"); for (int i = 0; i < endNodes[maxUpperBoundNodeIndex].Problem.Variables.Count; i++) { Console.WriteLine(" x" + i + " = " + endNodes[maxUpperBoundNodeIndex].Solve().Item1[i]); } Console.WriteLine("Optimal objective value = " + endNodes[maxUpperBoundNodeIndex].Solve().Item2); break; } //用最大上边界的节点来分支, 选择一个非整数解新增约束 Node parentNode = endNodes[maxUpperBoundNodeIndex]; endNodes.RemoveAt(maxUpperBoundNodeIndex); Tuple<List<double>, double> parentSolution = parentNode.Solve(); //选择一个离四舍五入整数解最远的非整数解来划分 int selectVariableIndex = 0; double distance = 0; for(int i = 0; i < parentSolution.Item1.Count; i++) { if (!isInteger(parentSolution.Item1[i]) && Math.Abs(Math.Round(parentSolution.Item1[i])-parentSolution.Item1[i])>distance) { selectVariableIndex = i; distance = Math.Abs(Math.Round(parentSolution.Item1[i]) - parentSolution.Item1[i]); } } double leftVar = Math.Floor(parentSolution.Item1[selectVariableIndex]); double rightVar = Math.Ceiling(parentSolution.Item1[selectVariableIndex]); //新分出两个节点 LooseLinearConstraintProblem leftSubProblem = new LooseLinearConstraintProblem(parentNode.Problem, new List<Tuple<int, double, double>> { new Tuple<int, double, double>(selectVariableIndex, double.NegativeInfinity, leftVar) }); Node leftNode = new Node(leftSubProblem); LooseLinearConstraintProblem rightSubProblem = new LooseLinearConstraintProblem(parentNode.Problem, new List<Tuple<int, double, double>> { new Tuple<int, double, double>(selectVariableIndex, rightVar,double.PositiveInfinity) }); Node rightNode = new Node(rightSubProblem); //分别计算两个子节点的结果 leftNode.Solve(); rightNode.Solve(); //更新当前的叶节点序列, 无解的节点不需要加入 if (leftNode.SolveStatus != SolveStatus.Infeasible) { endNodes.Add(leftNode); } if (rightNode.SolveStatus != SolveStatus.Infeasible) { endNodes.Add(rightNode); } //更新当前所有叶节点的上下界 List<int> needRemoveIndexList = new List<int>(); for(int i=0;i<endNodes.Count;i++) { var node = endNodes[i]; node.UpperBound = node.Solve().Item2; //如果该节点没有完整整数解, 那么它的下边界就是全局下界 if (node.SolveStatus == SolveStatus.NoneCompleteIntegerSolution) { node.LowerBound = globalLowerBound; } // 如果是一个完整整数解, 那么该节点的上下边界相等, 进一步, 如果这个目标值大于当前的 // 全局下界, 则更新全局下界 else { node.LowerBound = node.UpperBound; if (node.LowerBound >= globalLowerBound) { globalLowerBound = node.LowerBound; } } //如果在该节点下界大于上界，把这个结点从叶节点序列中去除，即剪枝 if (node.LowerBound > node.UpperBound) { needRemoveIndexList.Add(i); } } foreach(var needRemove in needRemoveIndexList) { endNodes.RemoveAt(needRemove); } Console.WriteLine($"LeftNode, additional constraint: x{selectVariableIndex + 1}<={leftVar}"); if (leftNode.SolveStatus == SolveStatus.Infeasible) { Console.WriteLine($"LeftNode is infeasible"); } else { Console.WriteLine($"LeftNode, UpperBound={leftNode.UpperBound}, LowerBound={leftNode.LowerBound}"); } Console.WriteLine($"RightNode, additional constraint: {rightVar}<=x{selectVariableIndex + 1}"); if (rightNode.SolveStatus == SolveStatus.Infeasible) { Console.WriteLine($"RightNode is infeasible"); } else { Console.WriteLine($"RightNode, UpperBound={rightNode.UpperBound}, LowerBound={rightNode.LowerBound}"); } } } //从一个double数列中找出所有为整数的元素的下标 public static List<int> getIntegerIndex(List<double> values) { List<int> results = new List<int>(); for(int i = 0; i < values.Count; i++) { if (isInteger(values[i])) { results.Add(i); } } return results; } public static bool isInteger(double n) { string s = Convert.ToString(n); return s.Contains(".") ? false : true; } public class LooseLinearConstraintProblem { private Solver solver; private List<Variable> variables; public List<Variable> Variables { get { return variables; } } //private List<Constraint> constraints; private Objective objective; //item1: 第几个变量 //item2: 大于等于多少 //item3: 小于等于多少 private List<Tuple<int, double, double>> additionalConstraints; public List<Tuple<int, double, double>> AdditionalConstraints { set { this.additionalConstraints = value; } get { return this.additionalConstraints; } } private void BasicInit() { this.solver = Solver.CreateSolver("GLOP"); this.variables = new List<Variable>(); Variable x1 = solver.MakeNumVar(0.0, double.PositiveInfinity, "x1"); Variable x2 = solver.MakeNumVar(0.0, double.PositiveInfinity, "x2"); this.variables.Add(x1); this.variables.Add(x2); Constraint c0 = solver.MakeConstraint(double.NegativeInfinity, 40000); c0.SetCoefficient(x1, 8000); c0.SetCoefficient(x2, 4000); Constraint c1 = solver.MakeConstraint(double.NegativeInfinity, 200); c1.SetCoefficient(x1, 15); c1.SetCoefficient(x2, 30); this.objective = solver.Objective(); this.objective.SetCoefficient(x1, 100); this.objective.SetCoefficient(x2, 150); this.objective.SetMaximization(); } public LooseLinearConstraintProblem(List<Tuple<int, double, double>> additionalConstraints) { BasicInit(); if (additionalConstraints != null) { this.additionalConstraints = additionalConstraints; foreach(var additionalConstraint in additionalConstraints) { Constraint newConstraint = this.solver.MakeConstraint(additionalConstraint.Item2, additionalConstraint.Item3); for (int j = 0; j < this.variables.Count; j++) { if (j != additionalConstraint.Item1) { newConstraint.SetCoefficient(this.variables[j], 0); } else { newConstraint.SetCoefficient(this.variables[j], 1); } } } } } public LooseLinearConstraintProblem(LooseLinearConstraintProblem parentProblem,List<Tuple<int, double, double>> additionalConstraints) { BasicInit(); this.additionalConstraints = new List<Tuple<int, double, double>>(); if (parentProblem.AdditionalConstraints != null) { foreach (var additionalConstraint in parentProblem.AdditionalConstraints) { this.additionalConstraints.Add(additionalConstraint); Constraint newConstraint = this.solver.MakeConstraint(additionalConstraint.Item2, additionalConstraint.Item3); for (int j = 0; j < this.variables.Count; j++) { if (j != additionalConstraint.Item1) { newConstraint.SetCoefficient(this.variables[j], 0); } else { newConstraint.SetCoefficient(this.variables[j], 1); } } } } if (additionalConstraints != null) { foreach (var additionalConstraint in additionalConstraints) { this.additionalConstraints.Add(additionalConstraint); Constraint newConstraint = this.solver.MakeConstraint(additionalConstraint.Item2, additionalConstraint.Item3); for (int j = 0; j < this.variables.Count; j++) { if (j != additionalConstraint.Item1) { newConstraint.SetCoefficient(this.variables[j], 0); } else { newConstraint.SetCoefficient(this.variables[j], 1); } } } } } public Tuple<List<double>, double> Solve() { var status = solver.Solve(); if (status != Solver.ResultStatus.OPTIMAL) { return null; } else { List<double> vResults = new List<double>(); foreach(var v in this.variables) { vResults.Add(v.SolutionValue()); } double oResult = solver.Objective().Value(); return new Tuple<List<double>, double>(vResults, oResult); } } } } }
最新回复(0)