运筹优化

运筹优化

本文为转载文章，原文地址：

(本文内容主要参考《an introduction to management science(13e)》中的章节)

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

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

老板的预算是$40000，并且现在有200平方的空地；工厂老板应该怎样购买来让每天的利润最大？
如果我们设x 1 x_1x1​是购买多少台印刷机，x 2 x_2x2​是购买多少台车床，那么可以写出优化问题的模型

相比线性优化模型，上面的数学模型就多出必须为整数的约束。如果我们把整数约束去掉，那么可以得到最优解:

如果我们把这个解四舍五入，可以得到一个可行的整数解

但是这个整数解并不是实际的最优解，因为一旦加上整数约束，线性优化模型往往是非凸的，线性优化算法往往无能为力。而分支界定法则是把原问题进行松弛(即去掉整数约束)，以松弛线性模型为根节点不断地划分子问题和更新问题上下界，通过评判上下界来忽略无效子问题和停止算法。首先我们构造根节点，该节点就是上面模型去掉整数约束后的松弛模型，对于这个松弛线性最大化问题，我们构造一下上下界，用它的最优解(Z = 1055.56)作为问题的上界，用四舍五入后的整数解作为下界；上下界的意义在于，原始整数优化模型的解一定是在上下界之间的，所以对于所有最优解比下界更小的子问题，我们可以直接忽略；而子问题的上界一定不会大于父问题，通过不断分支，这个上下界间的范围将不断缩减。

因为根节点的最优解不是整数解，所以我们选择一个非整数解进行分支。这里我们采取选择离四舍五入值最远的变量，对于根节点，即选择了变量x 2 = 5.56。我们以变量x 2 为划分点向下分支出两个子问题，分别在根节点的松弛优化模型的基础上新增约束x 2 ≤ 5和x 2 ≥ 6 ，也就是说，对于左侧的子节点2，松弛优化模型为：

右侧子节点3的问题模型为:

对于节点2，可以求得最优解为x 1 = 2.5 , x 2 = 5 , Z = 1000 ，而节点3的最优解为x 1 = 1.33 , x 2 = 6 , Z = 1033.33 ，这两个解可分别作为这两个节点的上界。可以看到这两个上界都小于根节点的上界，因为每一次分支，都是在原松弛优化模型的基础上新增约束，它的最优解肯定不会好过原问题。再来分析下界，因为这两个子节点的解都不是完整整数解，所以这两个节点下界保持不变，仍是950

我们还没有找到最优可行解，所以要从结点2或者3继续分支。选择哪个结点呢，这里我们按上界更大优先的原则来选取结点，也就是从结点3进行分支。采用这种原则其实是一种间接地剪枝，如果某个结点的最优解是非整数解并且小于当前的全局下界(即上界比下界更小)，那么需要把这个结点及其所有子节点都从搜索数里去除，它们的解肯定不会是最优解了，通过上界越大越优先的原则来选取，这些结点实际上根本不会被选中。结点3的x 2已经是整数了，所以用x 1​进行划分：x 1 ≤ 1和x 1 ≥ 2 ，分出了结点4和结点5，结点4的问题模型是在3的基础上加上x 1 ≤ 1 

结点5则是在3的基础上加上x 1 ≥ 2 

对于结点4，其最优解是(x 1 = 1 , x 2 = 6.17 , Z = 1025 )；对于结点5，则无解。那就只需要对结点4进行定界，上界就是最优解1025.5，同时最优解也不是完整整数解，那下界仍保持不变。

我们还是没有找到最优整数解，需要继续分支。现在有3个叶节点：2、4和5，结点5已经无解，所以不需要考虑，而在2和4中，4的上界更大，所以从结点4开始分支。在结点4，x 1已经是整数，所以用x 2 ​进行划分：x 2 ≤ 6 和x 2 ≥ 7 。结点6的问题模型如下，其最优解为(x 1 = 1 , x 2 = 6 , Z = 1000)

结点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);}}}