运筹系列3：整数规划分支定界法python代码

1. 模型

整数规划的模型与线性规划基本相同，只是额外的添加了部分变量为整数的约束。
我们假设问题为：
Min $cx$
s.t. $Ax\le b$
$A_{eq}x=b_{eq}$
$x\in Z$

2. 求解步骤

整数规划求解的基本框架是分支定界法（Branch and bound，BnB）。首先去除整数约束得到“松弛模型”，使用线性规划的方法求解。若有某个变量不是整数，在松弛模型上分别添加约束：

x ≤ floor(A)

和

x ≥ ceil(A)

然后再分别求解，这个过程叫做分支。当节点求解结果中所有变量都是整数时，停止分支。这样不断迭代，形成了一棵树。
所谓的定界，指的是叶子节点产生后，相当于给max问题定了一个下界（当然也可以提前用启发式算法给出下界）。之后在求解过程中一旦某个节点的目标函数值小于这个下界，那就直接pass，不用再进行分支了，这个步骤称为fathoming；每次新产生叶子节点，则更新下界。
fathoming的原则有三:

1. 分枝无解；
2. 分枝最优解<=lb；
3. 分枝最优解>lb，并且同时是整数解。此时更新lb=分枝最优解（bound）。
   第四种情况就是branch了。

伪代码如下：

如果我们要求的是近似解，比如Z大于等于Z*-K即可，那么我们只需要修改一下fathoming的原则：
2. 分枝最优解<=lb-K；并且在3之后执行

3. 算法实现

3.1 python代码

import math
from scipy.optimize import linprog
import sys

def integerPro(c, A, b, Aeq, beq,t=1.0E-12):
	# 求解松弛问题
	res = linprog(c, A_ub=A, b_ub=b, A_eq=Aeq, b_eq=beq)
    bestVal = sys.maxsize
    bestX = res.x
    if not(type(res.x) is float or res.status != 0): 
        bestVal = sum([x*y for x,y in zip(c, bestX)])
    # 停止条件 & bound
    if all(((x-math.floor(x))<t or (math.ceil(x)-x)<t) for x in bestX):
        return (bestVal,bestX)
    else:
    	# 进行branch，这里简单选择第一个非整数变量
        ind = [i for i, x in enumerate(bestX) if (x-math.floor(x))>t and (math.ceil(x)-x)>t][0]
        # branch出两个子问题
        newCon1 = [0]*len(A[0])
        newCon2 = [0]*len(A[0])
        newCon1[ind] = -1
        newCon2[ind] = 1
        newA1 = A.copy()
        newA2 = A.copy()
        newA1.append(newCon1)
        newA2.append(newCon2)
        newB1 = b.copy()
        newB2 = b.copy()
        newB1.append(-math.ceil(bestX[ind]))
        newB2.append(math.floor(bestX[ind]))
        r1 = integerPro(c, newA1, newB1, Aeq, beq)
        r2 = integerPro(c, newA2, newB2, Aeq, beq)
        # tree search，这里使用width first
        if r1[0] < r2[0]:
            return r1
        else:
            return r2

例子：输入

c = [3,4,1]
A = [[-1,-6,-2],[-2,0,0]]
b = [-5,-3]
Aeq = [[0,0,0]]
beq = [0]
print(integerPro(c, A, b, Aeq, beq))

输出

(8.0, array([2., 0., 2.]))

其中8是目标函数值，2，0，2是3个整数变量的值。
这里推荐使用对偶单纯形法或者内、外点法，不用每次都求解一个新问题。

3.2 Julia代码

using JuMP,GLPK,DataStructures
struct Node
    C::Vector{Int64} 
    A::Vector{Vector{Int64} }
    b::Vector{Int64} 
    x::Vector{Float64} 
    lb::Float64
end

function linPro(c, A, b)
    m = Model(GLPK.Optimizer) 
    @variable(m, mx[1:length(c)]>=0)
    @objective(m, Min, sum(c.*mx))
    @constraint(m, [i = 1:length(b)], sum(A[i].*mx) <= b[i])
    optimize!(m)
    if Int(termination_status(m))==1
        r_x = JuMP.value.(mx)
        r_val = objective_value(m)
        return r_x,r_val
    end
    return nothing,nothing
end

function integerPro(c, A, b, t=1.0E-12)
    PQ = PriorityQueue{Node,Float64}()
    x,lb = linPro(c,A,b)
    v = Node(c,A,b,x,lb)
    enqueue!(PQ,v,v.lb) 
    while length(PQ)>0
        v = dequeue!(PQ)
        ind = findfirst(xi->abs(round(xi)-xi)>t,v.x)
        if ind == nothing
            return v.x,v.lb
        else
            An=zeros(Int,length(c));An[ind]=1
            A1=[v.A;[-An]]
            A2=[v.A;[An]]
            b1=[v.b;-ceil(v.x[ind])]
            b2=[v.b;floor(v.x[ind])]
            x2,lb2 = linPro(c,A2,b2)
            x1,lb1 = linPro(c,A1,b1)
            if x1!=nothing;enqueue!(PQ,Node(c,A1,b1,x1,lb1),lb1);end
            if x2!=nothing;enqueue!(PQ,Node(c,A2,b2,x2,lb2),lb2);end
        end
    end
    return nothing,nothing
end

c = [3,4,1]
A = [[-1,-6,-2],[-2,0,0]]
b = [-5,-3]
integerPro(c, A, b)
## 输出([2.0, 0.0, 2.0], 8.0)

4. 官方文档示例图解

求解问题为：

4.1 求解松弛问题

4.2 对小数变量进行分支

4.3 迭代，剪枝

4.4 继续迭代

4.5 继续迭代，完成搜索