干货 | 10分钟掌握branch and cut算法原理附带C++求解TSP问题代码

简介: 干货 | 10分钟掌握branch and cut算法原理附带C++求解TSP问题代码

branch and cut其实还是和branch and bound脱离不了干系的。所以,在开始本节的学习之前,请大家还是要务必掌握branch and bound算法的原理

微信图片_20220422143035.gif


01 应用背景


Branch and cut is a method of combinatorial optimization for solving integer linear programs (ILPs), that is, linear programming (LP) problems where some or all the unknowns are restricted to integer values.
Branch and cut involves running a branch and bound algorithm and using cutting planes to tighten the linear programming relaxations.
Note that if cuts are only used to tighten the initial LP relaxation, the algorithm is called cut and branch.[1]


02 总体描述


前面说过,branch and cut其实还是和branch and bound脱离不了干系。其实是有很大干系的。在应用branch and bound求解整数规划问题的时候,如下图(好好复习一下该过程):

微信图片_20220422143042.png


假如,我们现在求一个整数规划最大化问题,在分支定界过程中,求解整数规划模型的LP松弛模型得到的非整数解作为上界,而此前找到的整数解作为下界。如果出现某个节点upper bound低于现有的lower bound,则可以剪掉该节点。否则,如果不是整数解继续分支。


此外,在求解整数规划模型的LP松弛时,If cutting planes are used to tighten LP relaxations。那么这时候的branch and bound就变成了branch and cut。


那么,什么是cutting planes呢?如下图:


微信图片_20220422143045.png


红色部分是整数规划的可行解空间。

蓝色部分是整数规划的LP松弛可行解空间。

在求解LP松弛时,加入橙色的cut,缩小解空间,同时又不影响整数解的解空间,可使解收敛得更快。


这就是branch and cut的过程了。比branch and bound高明之处就在于多了一个cutting planes,可能使branch and bound的效率变得更高。至于cutting planes是什么,等下一篇推文吧~


03 举个例子


为了让大家更好了解到branch and cut的精髓,必须得举一个简单的例子。对于同一个问题:

微信图片_20220422143047.png

branch and cut(左)和branch and bound(右)求解过程如下:


微信图片_20220422143050.jpg


可以看到,两者的不同之处就在子问题P2的处理上。


  • 对于branch and bound来说,求解线性松弛得到的Z = -29.5 < Z = -28。可知该支是可能隐含有更优解的,于是二话不说分支。无奈,分了两支以后发现居然没更优解,这种付出了却没有回报的感觉就像是受到了欺骗一样。

微信图片_20220422143053.jpg


  • 对于branch and cut来说,在求解线性松弛得到的Z = -29.5 < Z = -28时,并没有被兴奋冲昏头脑,它尝试着在线性松弛的解空间上砍下一块,但又不能影响到整数解的解空间范围。琢磨半天终于找到一块能砍的,于是Add cut: 2x1 + x2 <= 7。砍下来以后,形成新的子问题P3,赶紧看看P3的最优解是多少。在P3中,Z=-27.8 > -28,这一支果然不可取。


从上面的算法过程我们可以看到,求解同一个问题,branch and cut只用了3步,而branch and bound却用了4步。


There are many methods to solve the mixed-integer linear programming. Gomory Cutting Planes is fast, but unreliable. Branch and Bound is reliable but slow. The Branch and cut combine the advantages from these two methods and improve the defects.


It has proven to be a very successful approach for solving a wide variety of integer programming problems.


We can solve the MILP by taking some cutting planes before apply whole system to the branch and bound, Branch and cut is not only reliable, but faster than branch and bound alone. Finally, we understand that using branch and cut is more efficient than using branch and bound.[2]


04 算法过程


关于branch and cut的过程,可以总结如下:[1]


微信图片_20220422143055.png


相比branch and bound,其多了一个Cutting Planes的过程,先用Cutting Planes tighten LP relaxations,然后求解LP relaxations再判断是否有分支的必要。


其伪代码如下:

// ILP branch and cut solution pseudocode, assuming objective is to be maximized
ILP_solution branch_and_cut_ILP(IntegerLinearProgram initial_problem) {
    queue active_list; // L, above
    active_list.enqueue(initial_problem); // step 1
    // step 2
    ILP_solution optimal_solution; // this will hold x* above
    double best_objective = -std::numeric_limits<double>::infinity; // will hold v* above
    while (!active_list.empty()) { // step 3 above
        LinearProgram& curr_prob = active_list.dequeue(); // step 3.1
        do { // steps 3.2-3.7
            RelaxedLinearProgram& relaxed_prob = LP_relax(curr_prob); // step 3.2
            LP_solution curr_relaxed_soln = LP_solve(relaxed_prob); // this is x above
            bool cutting_planes_found = false;
            if (!curr_relaxed_soln.is_feasible()) { // step 3.3
                continue; // try another solution; continues at step 3
            }
            double current_objective_value = curr_relaxed_soln.value(); // v above
            if (current_objective_value <= best_objective) { // step 3.4
                continue; // try another solution; continues at step 3
            }
            if (curr_relaxed_soln.is_integer()) { // step 3.5
                best_objective = current_objective_value;
                optimal_solution = cast_as_ILP_solution(curr_relaxed_soln);
                continue; // continues at step 3
            }
            // current relaxed solution isn't integral
            if (hunting_for_cutting_planes) { // step 3.6
                violated_cutting_planes = search_for_violated_cutting_planes(curr_relaxed_soln);
                if (!violated_cutting_planes.empty()) { // step 3.6
                    cutting_planes_found = true; // will continue at step 3.2
                    for (auto&& cutting_plane : violated_cutting_planes) {
                        active_list.enqueue(LP_relax(curr_prob, cutting_plane));
                    }
                    continue; // continues at step 3.2
                }
            }
            // step 3.7: either violated cutting planes not found, or we weren't looking for them
            auto&& branched_problems = branch_partition(curr_prob);
            for (auto&& branch : branched_problems) {
                active_list.enqueue(branch);
            }
            continue; // continues at step 3
        } while (hunting_for_cutting_planes /* parameter of the algorithm; see 3.6 */
               && cutting_planes_found);
        // end step 3.2 do-while loop
    } // end step 3 while loop
    return optimal_solution; // step 4
}


相关代码


关于branch and cut 求解TSP问题的代码,请关注公众号【程序猿声】,在后台回复【bctsp】不包括【】即可获取。代码用了Gurobi,编译前配置好,不要出问题满世界跑火车说代码有问题。至于Gurobi,有时间再出教程吧。


注:对文中或者代码有疑问可联系小编,可提供有偿辅导服务。

// tsp.cc - traveling salesman code based on Gurobi using branch and cut
#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include <chrono>
#include <vector>
#include <deque>
#include <limits>
// Gurobi
#include "gurobi_c++.h"
// Data structure to represent an edge of the input graph
struct Edge
{
  int end1;
  int end2;
  double weight;
};
// If memory is limited, a soft limit for the maximum number of LPs in the
// queue can be set for the branch and cut algorithm.
// This is effectively disabled by default.
const int lp_soft_limit = 1000000;
// Tolerance to determine whether or not a number is integral
const double tol = 1.0e-9;
// Auxiliary function to decide whether or not a number is integral
inline bool is_integral(double d)
{
  return std::fabs(d - std::round(d)) < tol;
}
void find_components(int                       num_nodes,
                     int                       num_edges,
                     std::vector<Edge>         edges,
                     const std::vector<double> &x,
                     int                       &num_components,
                     std::vector<int>          &components)
{
  // Mark all unassigned nodes with -1.
  components.resize(num_nodes);
  for (int n = 0; n < num_nodes; n++) {
    components[n] = -1;
  }
  // Component index
  int c = 0;
  // Node indices
  int n1, n2;
  // Find all connected components.
  while (true) {
    // Find an unassigned node.
    n1 = -1;
    for (int n = 0; n < num_nodes; n++) {
      if (components[n] == -1) {
        n1 = n;
        break;
      }
    }
    if (n1 == -1) {
      // All nodes have been assigned.
      break;
    }
    // Assign node to current component.
    components[n1] = c;
    // Mark the entire connected component.
    while (true) {
      // Find an unassigned connected node.
      n2 = -1;
      for (int e = 0; e < num_edges; e++) {
        // Skip edges that are not used in the current solution.
        if (x[e] < tol) continue;
        if (edges[e].end1 == n1 && components[edges[e].end2] == -1) {
          n2 = edges[e].end2;
          break;
        }
        if (edges[e].end2 == n1 && components[edges[e].end1] == -1) {
          n2 = edges[e].end1;
          break;
        }
      }
      // No connected node found. Continue with next component.
      if (n2 == -1) break;
      // Assign the connected node to the current component.
      components[n2] = c;
      // Merge the two nodes.
      for (int e = 0; e < num_edges; e++) {
        if (edges[e].end1 == n2) {
          edges[e].end1 = n1;
        }
        if (edges[e].end2 == n2) {
          edges[e].end2 = n1;
        }
      }
    }
    c++;
  }
  num_components = c;
}
void solve_TSP(int                     num_nodes,
               int                     num_edges,
               const std::vector<Edge> &edges,
               bool                    integral_weights,
               std::vector<double>     &x_opt,
               int                     &lp_solves,
               int                     &subtour_constraints)
{
  // Allocate memory for the solution.
  std::vector<double> x;
  x.resize(num_edges);
  // Set up environment.
  GRBEnv env;
  // Create initial model.
  GRBModel initial_model(env);
  // Add variables.
  initial_model.addVars(num_edges, GRB_CONTINUOUS);
  initial_model.update();
  GRBVar *vars = initial_model.getVars();
  // Set up objective function.
  GRBLinExpr obj;
  for (int e = 0; e < num_edges; e++) {
    obj += GRBLinExpr(vars[e], edges[e].weight);
  }
  initial_model.setObjective(obj, GRB_MINIMIZE);
  // Add initial constraints.
  for (int n = 0; n < num_nodes; n++) {
    GRBLinExpr lhs;
    // Add all edges that are adjacent to the current node.
    for (int e = 0; e < num_edges; e++) {
      if (edges[e].end1 == n || edges[e].end2 == n) {
        lhs += GRBLinExpr(vars[e], 1.0);
      }
    }
    initial_model.addConstr(lhs, GRB_EQUAL, GRBLinExpr(2.0));
  }
  for (int e = 0; e < num_edges; e++) {
    initial_model.addConstr(GRBLinExpr(vars[e], 1.0),
                            GRB_LESS_EQUAL,
                            GRBLinExpr(1.0));
  }
  initial_model.update();
  lp_solves = 0;
  subtour_constraints = 0;
  // Branch and cut.
  double cost;
  double cost_opt = std::numeric_limits<double>::infinity();
  int num_components;
  std::vector<int> components;
  std::deque<GRBModel> problems;
  problems.push_back(initial_model);
  while (problems.size() > 0) {
    // Get next problem in queue.
    GRBModel model = problems.front();
    problems.pop_front();
    // In this loop, the LP is solved repeatedly until a solution without
    // subtours is found.
    bool skipped = false;
    while (true) {
      // Solve current model.
      model.optimize();
      lp_solves++;
      if (model.get(GRB_IntAttr_Status) != GRB_OPTIMAL) {
        // Do not continue branch if problem is infeasible.
        skipped = true;
        break;
      }
      // Check cost. If it is too high, stop following this branch.
      cost = model.get(GRB_DoubleAttr_ObjVal);
      // We can cut off branches more aggressively when the weights,
      // and thus the optimal cost, are integral.
      if (integral_weights && cost > cost_opt - 1.0 + tol) {
          skipped = true;
          break;
      }
      if (cost > cost_opt) {
        skipped = true;
        break;
      }
      // Get current solution.
      delete[] vars;
      vars = model.getVars();
      for (int e = 0; e < num_edges; e++) {
        x[e] = vars[e].get(GRB_DoubleAttr_X);
      }
      // Find connected components of the solution and eliminate subtours.
      find_components(num_nodes,
                      num_edges,
                      edges,
                      x,
                      num_components,
                      components);
      if (num_components == 1) {
        // There are no more subtours that could be eliminated.
        break;
      }
      // We will add one constraint per connected component.
      std::vector<GRBLinExpr> lhs;
      lhs.resize(num_components);
      for (int e = 0; e < num_edges; e++) {
        // Identify the component this edge belongs to.
        const int c1 = components[edges[e].end1];
        const int c2 = components[edges[e].end2];
        // Skip edges that connect two different components.
        if (c1 != c2) continue;
        // Add edge to the subtour elimination constraint.
        lhs[c1] += GRBLinExpr(vars[e], 1.0);
      }
      // Compute the size of each component.
      // This is required for the right-hand side of the constraints.
      std::vector<int> component_sizes;
      component_sizes.resize(num_components);
      for (int c = 0; c < num_components; c++) {
        component_sizes[c] = 0;
      }
      for (int n = 0; n < num_nodes; n++) {
        component_sizes[components[n]]++;
      }
      // Add constraints to model.
      for (int c = 0; c < num_components; c++) {
        model.addConstr(lhs[c],
                        GRB_LESS_EQUAL,
                        GRBLinExpr(component_sizes[c] - 1.0));
        subtour_constraints++;
      }
      std::cout << "Added " << num_components
                << " subtour elimination constraints." << std::endl;
      model.update();
    }
    if (skipped) continue;
    // Branch using a fractional variable.
    bool integral_sol = true;
    for (int e = 0; e < num_edges; e++) {
      if (!is_integral(x[e])) {
        integral_sol = false;
        // Add model with <= constraint for fractional variable.
        GRBModel model_le(model);
        model_le.addConstr(GRBLinExpr(vars[e], 1.0),
                           GRB_LESS_EQUAL,
                           GRBLinExpr(std::floor(x[e])));
        model_le.update();
        // Add model with >= constraint for fractional variable.
        GRBModel model_ge(model);
        model_ge.addConstr(GRBLinExpr(vars[e], 1.0),
                           GRB_GREATER_EQUAL,
                           GRBLinExpr(std::ceil(x[e])));
        model_ge.update();
        // Check if the soft limit for the number of LPs is hit.
        if (problems.size() < lp_soft_limit) {
          // Add new problems at the end of the queue.
          // This corresponds to breadth-first search.
          problems.push_back(model_le);
          problems.push_back(model_ge);
        }
        else {
          // Add new problems at the beginning of the queue.
          // This corresponds to depth-first search.
          problems.push_front(model_le);
          problems.push_front(model_ge);
        }
        // Print information about the queue.
        std::cout << "Branching; there are now "
                  << problems.size()
                  << " models in the queue."
                  << std::endl << std::endl;
        // Stop after creating one branch!
        break;
      }
    }
    // Update optimal cost and optimal solution if integral solution was found.
    if (integral_sol) {
      cost_opt = cost;
      x_opt = x;
    }
  }
}
int main(int argc, char **argv)
{
  // Read problem from stdin.
  std::string line;
  // Read problem size.
  std::getline(std::cin, line);
  // Remove leading spaces.
  while (line[0] == ' ') line = line.substr(1);
  const int num_nodes = std::stoi(line.substr(0, line.find(" ")));
  const int num_edges = std::stoi(line.substr(line.find(" ") + 1));
  // Read graph.
  std::vector<Edge> edges;
  edges.resize(num_edges);
  for (int e = 0; e < num_edges; e++) {
    std::getline(std::cin, line);
    // Remove leading spaces.
    while (line[0] == ' ') line = line.substr(1);
    edges[e].end1 = std::stoi(line.substr(0, line.find(" ")));
    line = line.substr(line.find(" ") + 1);
    edges[e].end2 = std::stoi(line.substr(0, line.find(" ")));
    line = line.substr(line.find(" ") + 1);
    edges[e].weight = std::stod(line);
  }
  std::cout << "Loaded TSP with " << num_nodes << " nodes and "
                                  << num_edges << " edges.\n";
  // Check if the edge weights are integral.
  // If so, we can optimize a bit more aggressively in some places.
  bool integral_weights = true;
  for (int e = 0; e < num_edges; e++) {
    if (!is_integral(edges[e].weight)) {
      integral_weights = false;
      break;
    }
  }
  std::cout << "Computation begins.\n";
  // Start timer.
  const auto t_start = std::chrono::high_resolution_clock::now();
  // Solve TSP using Gurobi (for the LPs).
  std::vector<double> x_opt;
  int lp_solves;
  int subtour_constraints;
  try {
    solve_TSP(num_nodes,
              num_edges,
              edges,
              integral_weights,
              x_opt,
              lp_solves,
              subtour_constraints);
  }
  catch (const GRBException &e) {
    std::cerr << "Gurobi exception: " << e.getMessage() << std::endl;
    std::exit(1);
  }
  // Stop timer.
  const auto t_end = std::chrono::high_resolution_clock::now();
  const std::chrono::duration<double> dtime = t_end - t_start;
  std::cout << "Computation finished (" 
            << std::fixed << std::setprecision(3)
            << dtime.count() << "s).\n";
  // Print additional information.
  std::cout << "Solved a total of " << lp_solves << " LPs." << std::endl;
  std::cout << "Added a total of " << subtour_constraints
            << " subtour elimination constraints." << std::endl;
  // Print optimal solution.
  std::cout << "The best tour is:\n";
  double c_optimal = 0.0;
  // Set output format.
  if (integral_weights) {
    std::cout << std::setprecision(0);
  }
  else {
    std::cout << std::setprecision(1);
  }
  for (int e = 0; e < num_edges; e++) {
    // See if the edge is used.
    if (x_opt[e] > 0.0) {
      std::cout << edges[e].end1 << " "
                << edges[e].end2 << " "
                << edges[e].weight << std::endl;
      c_optimal += x_opt[e]*edges[e].weight;
    }
  }
  std::cout << "The cost of the best tour is: " << c_optimal << std::endl;
  return 0;
}

reference


  • [1]  (https://en.wikipedia.org/wiki/Branch_and_cut)
  • [2] (https://optimization.mccormick.northwestern.edu/index.php/Branch_and_cut)


相关文章
|
6天前
|
机器学习/深度学习 存储 算法
近端策略优化(PPO)算法的理论基础与PyTorch代码详解
近端策略优化(PPO)是深度强化学习中高效的策略优化方法,广泛应用于大语言模型的RLHF训练。PPO通过引入策略更新约束机制,平衡了更新幅度,提升了训练稳定性。其核心思想是在优势演员-评论家方法的基础上,采用裁剪和非裁剪项组成的替代目标函数,限制策略比率在[1-ϵ, 1+ϵ]区间内,防止过大的策略更新。本文详细探讨了PPO的基本原理、损失函数设计及PyTorch实现流程,提供了完整的代码示例。
114 10
近端策略优化(PPO)算法的理论基础与PyTorch代码详解
|
1月前
|
机器学习/深度学习 算法 PyTorch
深度强化学习中SAC算法:数学原理、网络架构及其PyTorch实现
软演员-评论家算法(Soft Actor-Critic, SAC)是深度强化学习领域的重要进展,基于最大熵框架优化策略,在探索与利用之间实现动态平衡。SAC通过双Q网络设计和自适应温度参数,提升了训练稳定性和样本效率。本文详细解析了SAC的数学原理、网络架构及PyTorch实现,涵盖演员网络的动作采样与对数概率计算、评论家网络的Q值估计及其损失函数,并介绍了完整的SAC智能体实现流程。SAC在连续动作空间中表现出色,具有高样本效率和稳定的训练过程,适合实际应用场景。
216 7
深度强化学习中SAC算法:数学原理、网络架构及其PyTorch实现
|
2月前
|
算法 Java 数据库
理解CAS算法原理
CAS(Compare and Swap,比较并交换)是一种无锁算法,用于实现多线程环境下的原子操作。它通过比较内存中的值与预期值是否相同来决定是否进行更新。JDK 5引入了基于CAS的乐观锁机制,替代了传统的synchronized独占锁,提升了并发性能。然而,CAS存在ABA问题、循环时间长开销大和只能保证单个共享变量原子性等缺点。为解决这些问题,可以使用版本号机制、合并多个变量或引入pause指令优化CPU执行效率。CAS广泛应用于JDK的原子类中,如AtomicInteger.incrementAndGet(),利用底层Unsafe库实现高效的无锁自增操作。
理解CAS算法原理
|
1月前
|
算法 决策智能
基于SA模拟退火优化算法的TSP问题求解matlab仿真,并对比ACO蚁群优化算法
本项目基于MATLAB2022A,使用模拟退火(SA)和蚁群优化(ACO)算法求解旅行商问题(TSP),对比两者的仿真时间、收敛曲线及最短路径长度。SA源于金属退火过程,允许暂时接受较差解以跳出局部最优;ACO模仿蚂蚁信息素机制,通过正反馈发现最优路径。结果显示SA全局探索能力强,ACO在路径优化类问题中表现优异。
|
2月前
|
存储 算法 程序员
C 语言递归算法:以简洁代码驾驭复杂逻辑
C语言递归算法简介:通过简洁的代码实现复杂的逻辑处理,递归函数自我调用解决分层问题,高效而优雅。适用于树形结构遍历、数学计算等领域。
|
2月前
|
存储 人工智能 缓存
【AI系统】布局转换原理与算法
数据布局转换技术通过优化内存中数据的排布,提升程序执行效率,特别是对于缓存性能的影响显著。本文介绍了数据在内存中的排布方式,包括内存对齐、大小端存储等概念,并详细探讨了张量数据在内存中的排布,如行优先与列优先排布,以及在深度学习中常见的NCHW与NHWC两种数据布局方式。这些布局方式的选择直接影响到程序的性能,尤其是在GPU和CPU上的表现。此外,还讨论了连续与非连续张量的概念及其对性能的影响。
99 3
|
3月前
|
并行计算 算法 测试技术
C语言因高效灵活被广泛应用于软件开发。本文探讨了优化C语言程序性能的策略,涵盖算法优化、代码结构优化、内存管理优化、编译器优化、数据结构优化、并行计算优化及性能测试与分析七个方面
C语言因高效灵活被广泛应用于软件开发。本文探讨了优化C语言程序性能的策略,涵盖算法优化、代码结构优化、内存管理优化、编译器优化、数据结构优化、并行计算优化及性能测试与分析七个方面,旨在通过综合策略提升程序性能,满足实际需求。
99 1
|
2月前
|
算法 决策智能
基于遗传优化算法的TSP问题求解matlab仿真
本项目使用遗传算法解决旅行商问题(TSP),目标是在四个城市间找到最短路径。算法通过编码、选择、交叉、变异等步骤,在MATLAB2022A上实现路径优化,最终输出最优路径及距离。
|
4天前
|
算法 数据安全/隐私保护 计算机视觉
基于FPGA的图像双线性插值算法verilog实现,包括tb测试文件和MATLAB辅助验证
本项目展示了256×256图像通过双线性插值放大至512×512的效果,无水印展示。使用Matlab 2022a和Vivado 2019.2开发,提供完整代码及详细中文注释、操作视频。核心程序实现图像缩放,并在Matlab中验证效果。双线性插值算法通过FPGA高效实现图像缩放,确保质量。
|
1月前
|
算法 数据安全/隐私保护 计算机视觉
基于Retinex算法的图像去雾matlab仿真
本项目展示了基于Retinex算法的图像去雾技术。完整程序运行效果无水印,使用Matlab2022a开发。核心代码包含详细中文注释和操作步骤视频。Retinex理论由Edwin Land提出,旨在分离图像的光照和反射分量,增强图像对比度、颜色和细节,尤其在雾天条件下表现优异,有效解决图像去雾问题。