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

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

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 总体描述

03 举个例子

branch and cut（左）和branch and bound（右）求解过程如下：

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

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

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 算法过程

// 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
}

// 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);
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);
for (int n = 0; n < num_nodes; n++) {
GRBLinExpr lhs;
for (int e = 0; e < num_edges; e++) {
if (edges[e].end1 == n || edges[e].end2 == n) {
lhs += GRBLinExpr(vars[e], 1.0);
}
}
}
for (int e = 0; e < num_edges; e++) {
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]]++;
}
for (int c = 0; c < num_components; 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);
GRB_LESS_EQUAL,
GRBLinExpr(std::floor(x[e])));
model_le.update();
// Add model with >= constraint for fractional variable.
GRBModel model_ge(model);
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)
{
std::string line;
std::getline(std::cin, line);
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));
std::vector<Edge> edges;
edges.resize(num_edges);
for (int e = 0; e < num_edges; e++) {
std::getline(std::cin, line);
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";
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)

|
7天前
|

【8月更文挑战第6天】在机器学习领域，支持向量机（SVM）犹如璀璨明珠。它是一种强大的监督学习算法，在分类、回归及异常检测中表现出色。SVM通过在高维空间寻找最大间隔超平面来分隔不同类别的数据，提升模型泛化能力。为处理非线性问题，引入了核函数将数据映射到高维空间。SVM在文本分类、图像识别等多个领域有广泛应用，展现出高度灵活性和适应性。
59 2
|
2天前
|

python与朴素贝叶斯算法（附示例和代码）

5 0
|
9天前
|

【8月更文挑战第4天】字符串最小周期问题旨在找出字符串中最短重复子串的长度。KPM（实为KMP，Knuth-Morris-Pratt）算法，虽主要用于字符串匹配，但其生成的前缀函数（next数组）也可用于求解最小周期。核心思想是构建LPS数组，记录模式串中每个位置的最长相等前后缀长度。对于长度为n的字符串S，其最小周期T可通过公式ans = n - LPS[n-1]求得。通过分析周期字符串的特性，可证明该方法的有效性。提供的C++示例代码展示了如何计算给定字符串的最小周期，体现了KPM算法在解决此类问题上的高效性。
19 0
|
23天前
|

18 0
|
27天前
|

MCKP-MMF算法是一种启发式流量估计方法，用于寻找无线传感器网络的局部最优解。它从最小配置开始，逐步优化部分解，调整访问点的状态。算法处理访问点的动态影响半径，根据带宽需求调整，以避免拥塞。在MATLAB 2022a中进行了仿真，显示了访问点半径请求变化和代价函数随时间的演变。算法分两阶段：慢启动阶段识别瓶颈并重设半径，随后进入周期性调整阶段，追求最大最小公平性。
32 1
|
11天前
|

30 4
|
18天前
|

36 12
|
9天前
|

23 2
|
19天前
|

1. **功能描述**: 提出了一个去抖动缓冲区感知调度器，结合用户终端的缓冲状态减少服务中断。该算法通过动态调整数据包发送速率以优化网络延迟和吞吐量。 2. **测试结果**: 使用MATLAB 2022a进行了仿真测试，结果显示De-Jitter Buffer算法在网络拥塞时比RR调度算法更能有效利用资源，减少延迟，并能根据网络状态动态调整发送速率。 3. **核心程序**: MATLAB代码实现了调度逻辑，包括排序、流量更新、超时和中断处理等功能。 仿真结果和算法原理验证了De-Jitter Buffer算法在无线网络调度中的优势。
36 4
|
23天前
|

41 4