运筹系列4：整数规划割平面法python代码

1. 分支割平面（branch and cut）

割平面简单来说，就是添加约束条件。
Cuts are constraints added to a model to restrict (cut away) noninteger solutions that would otherwise be solutions of the continuous relaxation. The addition of cuts usually reduces the number of branches needed to solve a MIP.

在使用分支定界法时，我们一般是首先尝试添加各种割平面后看能不能求出整数解，如果不行再分支。这种方法叫做Branch and Cut。
首先介绍一些基本的cut方法：

1. rounding：
   比如整数变量$x\le 1.5$可以转化为$x\le 1$；
   还有GCD reduction：比如$3x+6y+9z\le 11$，同时除以3后有$x+2y+3z\le 3$；
   Gomory rouding cut：比如$3x+3y+5z\le 8$，同时除以3有$x+y+5/3z\le 8/3$，两边rounding得到$x+y+z\le 2$
2. lifting：
   比如$4x+y\ge 2,x$ binary，x从0提升到1，左边slack为2，x前面的系数可以减去2变为$2x+y\ge 2$
3. disjunction：
   比如$x+y\ge 3.5,x,y\ge 0,y$ integer，可以把y拆分为$y\ge 4$和$y\le 3$两部分。

2. cut介绍

剪枝方法分为对约束形式有要求的特殊剪枝以及通用的剪枝：

• Generic Cuts (valid for any MILP)
  Gomory Mixed Integer
  Mixed Integer Rounding
  Disjunctive cut
• Special Structures (valid for certain relaxations of MILPs)
  Knapsack / Gub Cover, Pack (many applications)
  Flow Cover / Path (fixed charge network flow, lot-sizing, …)
  Cliques / Odd-Hole (set partitioning, covering, packing)
  Implied Bound (logical implications between binary variables)

2.1 Disjunctive cut

在分支定界算法中，添加的x≤floor[x_s]和x≥ceil[x_s]便是两个用来割平面的约束条件，floor[x]和ceil[x]之间的整个可行域在对x进行分支的过程中被切割掉了，称为disjunctive cuts。

2.2 Gomory cut

Gomory的思想是将等式两边系数的小数部分拿出来做成新的约束，主要适用于纯整数规划。
假设$ax=b,x\in Z$，则一定有$\lfloor a\rfloor x\le \lfloor b\rfloor$，然后两者相减，得到：$f_{a}x\ge f_b$；反过来有$y+\lceil a\rceil x\ge \lceil b\rceil$，然后两者相减，得到：$(1-f_a)x\ge 1-f_b$。下面是个例子：

Gomory可以与其他方法组合使用。假设整数规划的线性松弛问题求解结果中有一个基变量$y=b$不是整数，对应约束：$y+\Sigma a_j{x}_j+\Sigma a_k{x}_k=b$
令：$y+\Sigma \lfloor a_j\rfloor x_j+\Sigma \lceil a_k\rceil x_k=t$（两边取整）
相减得到：$\Sigma f_{a_j}{x}_j-\Sigma (1-f_{a_k})x_k=b-t$
disjunction，有
$b-t\ge 0=>\Sigma f_{a_j}{x}_j\ge f_b$，
$b-t<0=>\Sigma (1-f_{a_k})x_k\ge 1-f_b$
合并得到 Σ f a j x j / f b + Σ ( 1 − f a k ) x k / ( 1 − f b ) ≥ max ⁡ { Σ f a j x j / f b , Σ ( 1 − f a k ) x k / ( 1 − f b ) } ≥ 1 \Sigma f_{a_j}x_j / f_b+\Sigma (1-f_{a_k})x_k/(1-f_b)\ge \max\{\Sigma f_{a_j}x_j / f_b,\Sigma (1-f_{a_k})x_k/(1-f_b)\}\ge1 Σfaj​​xj​/fb​+Σ(1−fak​​)xk​/(1−fb​)≥max{Σfaj​​xj​/fb​,Σ(1−fak​​)xk​/(1−fb​)}≥1
我们一般选取$f_{a_j}\le f_b,f_{a_k}>f_b$，这样系数都小于1，约束比较紧。

在求解器中，一般只应用于根节点，挑选非整数最严重的一些变量（比如100个）添加gomory割平面到松弛问题上，然后重复两遍。

2.3 MIR cut

MIR cuts are generated by applying integer rounding on the coefficients of integer variables and the righthand side of a constraint.
核心思想是对整数变量与非整数变量>=0的交点进行distructive cut，然后连接cut的两个交点形成一个新的约束条件。

Mix integer rounding针对的是如下问题：
若$y\le b+x,y\in Z$，我们可以添加cut：
$y\le \lfloor b\rfloor +x/(1-f_b)$
证明：
若$f_x+f_b<1$，则原约束满足$y\le \lfloor b\rfloor +\lfloor x\rfloor \le \lfloor b\rfloor +x/(1-f_b)$
若$f_x+f_b\ge 1$，则原约束满足$y\le \lfloor b\rfloor +1+\lfloor x\rfloor \le \lfloor b\rfloor +(f_x+\lfloor x\rfloor )/f_x\le \lfloor b\rfloor +x/(1-f_b)$

反过来，我们也有：
若$y+x\ge b,y\in Z$，我们可以添加cut：
$y+x/f_b\ge \lceil b\rceil$
如下图：

MIR的一个扩展是：
对于问题$\Sigma a_j{y}_j+x-z\le \lfloor b\rfloor$，我们可以添加cut：
$\Sigma (\lfloor a_j\rfloor +(f_{a_j}-f_b{)}^{+}/(1-f_b))y_j\le \lfloor b\rfloor +z/(1-f_b)$
证明：
若$f_{a_j}\ge f_b$，则$(f_{a_j}-f_b)/(1-f_b)\le (f_{a_j}-f_b{f}_{a_j})/(1-f_b)=f_{a_j}$
若$f_{a_j}<f_b$，使用rounding，必有$\lfloor a_j\rfloor y_j\le \lfloor b\rfloor$

2.4 Cover cut

cover cut有多种，这里介绍knapsack cover cut，针对如下约束：
$\Sigma ax\le b,x$ binary
若集合$C$满足${\Sigma }_{C}a>b$，则把$C$称为一个cover，cover cut为：
${\sigma }_C{x}_j\le |C|-1$

cover cut也可以和其他方法进行结合，如下：

2.5 clique cut

若一系列binary变量两两互斥，则可以生成clique cut，即：
若$x+y\le 1,z+y\le 1,x+z\le 1$
则$x+y+z\le 1$

3.从User cut到Lasy cut

无论是普通B&B中的对变量进行切分，还是B&C用约束条件的小数部分形成切分，切分条件对于原问题都是符合的，称为User cut。如果原问题还有一些隐含的额外约束可以作为cut，这些cut称为Lasy cut。
另外有种情况会使用到lasy cut，那就是有很多约束在最开始的时候是冗余的，这些约束会被放进一个pool中，时不时拿出来检查一下，如果被违反了，就加入约束中；如果有一段时间不违反了，再把它放回pool中。这样可以减少每次迭代的计算量

4.python代码

基本框架还是用分支定界法，每次求解完之后添加割平面的约束条件：

def add_new_restriction(matrix):
    new_column = np.zeros(matrix.shape[0]+1)
    new_line = np.zeros(matrix.shape[1])
    new_column[-1] = -1 
    #这里简单使用第一行约束条件为基础生成新约束条件。
    new_line = matrix[1, :] 
    for index in range(0, len(new_line)):
        number = np.array(new_line[index], dtype=float)
        if number.tolist().is_integer() == False:
            new_line[index] = math.floor(new_line[index])
    matrix = np.insert(matrix, matrix.shape[0], new_line, axis=0)
    matrix = np.insert(matrix, -1, new_column, axis=1)
    return matrix