编程手养成计划 Day17-19-CSDN博客

本文链接：https://blog.csdn.net/weixin_51939900/article/details/131749201

说明：为了防止我写的笔记误导新入门其他人，所以后续笔记就不会公开了。

本文参考自知乎王源博主的文章，推荐大家去看，真的很好【整数规划(一)】整数规划问题综述

一、什么是整数规划

（一）整数规划问题定义

在现实生活当中，整数规划的使用率很高，他主要解决的是在数学规划问题中一部分变量或者全部变量为整数变量的情况。例如：要表示人数，产品数等无法用小数。
线性规划的问题可以表示为：
$\min_xc^Tx \\ s.t. \left\{ \begin{aligned} Ax&\leq b\\ x&\geq0\\ x & \in Z^n \end{aligned} \right.$

（二）常见整数规划问题举例

1. 背包问题

设有一个背包，其最大承重为 $b$ ，考虑 $n$ 件物品，其中第 $j$ 件重量为 $a_j$ ，其价值为 $c_j$ 。问如何选取物品装入背包中，使得背包内物品的总价值最大？

2. 广义指派问题

设有 $m$ 台机器， $n$ 个工件，第 $i$ 台机器的可用工时数为 $b_i$ ，第 $i$ 台机器完成第 $j$ 件工件需要的工时数为 $a_ij$ ，费用为 $c_ij$ 。问如何最优指派机器生产。

3. 集合覆盖问题

设某地区划分为若干个区域，需要建立若干个应急服务中心（如消防站，急救中心等），每个中心的建立都需要一笔建站费用，设候选中心的位置已知，每个中心可以服务的区域预先知道，问如何选取中心使得应急服务能覆盖整个地区且使得建站费用最小。

二、解决方案

以下内容参考：【整数规划(二)】整数规划中的"简单"问题之全单模矩阵(totally unimodular matrix)
详细证明请看原文

（一）全单模矩阵

1. 整数规划问题的连续松弛和原整数规划问题的关系

假设有如下整数规划问题S： $\min\{c^Tx:Ax\leq b, x\in Z^n_+\}$
其连续松弛问题S1： $\min\{c^Tx:Ax\leq b, x\in R^n_+\}$
可以理解为将整数规划拓展到线性规划中

不难发现，一定有S1的最优解 $\geq$ S的最优解。
那么如果两式的最优解等价，就可以转换模型，极大的降低了求解难度。

2. 等号成立的条件

如果矩阵A为全单模，那么整数规划等于线性规划。

其中全单模矩阵为：设矩阵 $A_{m\times n}$ ，若矩阵 $A$ 的任意子方阵的行列式为0,-1,+1，则称矩阵 $A$ 为全单模矩阵。（子方阵是指把一个矩阵的一部分行和一部分列的交点所构成的方阵）

但是，由于全单模矩阵的判定，如果使用定义，由于需要计算每一个子方阵，算法复杂度呈指数级别，无法直接使用。于是我们想到了用性质来判别：

性质1： 若矩阵 $A$ 是全单模矩阵，则矩阵中的元素只能为0，-1或者+1。

其逆否命题为：若矩阵中有任意一个元素不等于0，-1或者+1，则矩阵 $A$ 不是全单模矩阵。好了这个逆否命题是非常实用的，一下子就把很多矩阵排除出了全单模矩阵的范畴。

性质2： 设整数矩阵 $A$ 是全单模矩阵，对 $A$ 进行一下运算不改变其全单模性质：

(1) 对矩阵 $A$ 进行转置
(2) 矩阵 $(A, I)$ 是全单模的
(3) 去掉 $A$ 的一行或者一列
(4) 将 $A$ 的一行或者一列乘以 $- 1$
(5) 互换 $A$ 的两行或者两列
(6) 对 $A$ 进转轴运算

推论1： 设矩阵 $A$ 的任意元素都是0，-1或者+1，并且每列至多有2个非0元素，则矩阵 $A$ 是全单模矩阵当且仅当存在 $A$ 的行分割 $Q_1,Q_2$ 使得同一列中的两个非0元素满足以下条件：

(1) 若符号相同，则一个元素位于 $Q_1$ ，另一个元素位于 $Q_2$
(2) 若符号相反，则这两个元素同时属于 $Q_1$ 或者同时属于 $Q_2$

推论2： 设矩阵 $A$ 的任意元素都是0，-1或者+1，若 $A$ 满足以下两个条件则矩阵 $A$ 是全单模的；

(1) $A$ 的每一列至多含有2个非0元素
(2) 若某列含有2个非0元素，则两个元素之和为0

最后强调一点就是推论1和2都是充分条件，也就是说满足推论1和2的那肯定是全单模矩阵，但是不满足的也可能是全单模矩阵也可能不是全单模矩阵，这一点大家在使用的时候需要注意。

另外，还可以使用如下方法辅助判定：
给定一个矩阵 $A$ ，然后随机生成向量 $c^T$ ，然后求解如下的线性规划问题： $\min\{c^Tx:Ax\leq b,x\in R^n_+\}$ 。重复以上步骤多次，如果每次得到的线性规划的最优解都是整数解的话，那么大概率 $A$ 是全单模的。

注意，本质上它是一种蒙特卡洛采样的方法，如果采样的 $N$ 越大我们就越有信心相信该方法的正确性。当然要说的是本方法只是一个技巧，不能严谨的去保证全单模性。

3. 常见的全单模矩阵整数规划问题实例

（1）二部图

$G = (V, E)$ 表示无向图， $M$ 表示图的 $V\times E$ 关联矩阵，则矩阵 $M$ 是全单模当且仅当图 $G$ 是二部图。

二部图：通俗来讲就是，把一个图的顶点划分为两个不相交子集，使得每一条边都分别连接两个集合中的顶点。如果存在这样的划分，则此图为一个二分图。

关于二部图有一个重要的定理：G为二部图的充要条件是G中的每一个圈的长度都是偶数。

下面举一些二部图的例子：
在这里插入图片描述

（2）指派问题

将 $n$ 项任务分配给 $n$ 个工人，每个人有且只能分配到一个任务。工人 $i$ 完成任务 $j$ 的成本为 $c_{ij}$ ，如何使得分配任务的总成本最小。
$\underset{x}{\min}\sum_{i=1}^n{\sum_{j=1}^n{c_{ij}x_{ij}}}\\ s.t. \;\; \sum_{j=1}^n{x_{ij}}=1,\ \forall i=1,...,n\\ \sum_{i=1}^n{x_{ij}}=1,\ \forall j=1,...,n\\ x_{ij}\in \left\{ 0,1 \right\} ,\ \forall i,j=1,...,n$
本质上它就是前面二部图问题的一种特殊情况，工人和任务分别就是二部图的两部分，然后我们实际上就是对这两部分做一个匹配。所以指派问题是全单模整数规划问题，因此可以通过求解线性规划来得到整数最优解。

（3）最小费用最大流问题

这篇文章对最小费用最大流问题做了很详细的解释，我会在后面针对网络流问题问题写一篇笔记，这里就不再重复了。
【网络流优化(一)】网络流问题综述

（二）分支定界法

1. 算法概述

上面我们在全单模矩阵中计算了原整数规划问题的连续松弛问题，如果恰好是整数，那么就可以转化为线性规划。
下面我们假设一种情况，我们解出来的最优解 $x_0 = 3.5$ 。因为我们希望结果是整数，那么此时我们可以基于 $x_0=3.5$ 将原可行域划分为两部分，一个是 $x_0\leq3$ ，一个是 $x_0\geq4$ ，并以这两个条件与原规划问题形成两个新的规划问题。可以表示为：

在这里插入图片描述

我们将可行域 $S$ ，转化为 $S_1$ 和 $S_2$ 。对于多变量，我们还可以继续向下迭代，然后遍历所有节点，那么我们就能找到最优解。但是很明显对于这种指数增长的模型，遍历所有节点的复杂度极高，那么有没有什么办法能够改进呢？

以上图为例，以 $S$ 为可行域得到的最优解一定会优于或等于以 $S_1$ 或以 $S_2$ 为可行域的解。因为 $S>S_1\cup S_2$ 恒成立，就像本来我们要从一堆产品中找出最优的，现在是要在原来的基础上剔除掉一些再找一个最优的，那么不难想象，我们此时最理想的状态就是原先的最优解，不会更好。

那么，我们不难得出，越向下迭代，最优解只会越来越差，因为迭代便意味着要新增加限制条件，可行域范围缩小。我们可以想到每计算一次便更新一次最优解的范围，那么如果我后面算到的最优解比我先前得到的要差，便可以舍弃那一条支路。那么最后经过迭代，我们就可以得到最优解。

这里只做一个简单的算法原理解释，更详细的解释和证明请看：
【整数规划(三)】分支定界法及其代码实现

2. 代码实现

然后这里以一道例题来放出代码：
$\max z=40x_1 + 90x_2\\ s.t.\left\{ \begin{aligned} 9x_1+7x_2&\leq 56\\ 7x_1+20x_2&\leq 70\\ x_1,x_2 &\geq 0,且为整数 \end{aligned} \right.$

from scipy.optimize import linprog
import numpy as np
import math
import sys
from queue import Queue


class ILP():
    def __init__(self, c, A_ub, b_ub, A_eq, b_eq, bounds):
        # 全局参数
        self.LOWER_BOUND = -sys.maxsize
        self.UPPER_BOUND = sys.maxsize
        self.opt_val = None
        self.opt_x = None
        self.Q = Queue()

        # 这些参数在每轮计算中都不会改变
        self.c = -c
        self.A_eq = A_eq
        self.b_eq = b_eq
        self.bounds = bounds

        # 首先计算一下初始问题
        r = linprog(-c, A_ub, b_ub, A_eq, b_eq, bounds)

        # 若最初问题线性不可解
        if not r.success:
            raise ValueError('Not a feasible problem!')

        # 将解和约束参数放入队列
        self.Q.put((r, A_ub, b_ub))

    def solve(self):
        while not self.Q.empty():
            # 取出当前问题
            res, A_ub, b_ub = self.Q.get(block=False)

            # 当前最优值小于总下界，则排除此区域
            if -res.fun < self.LOWER_BOUND:
                continue

            # 若结果 x 中全为整数，则尝试更新全局下界、全局最优值和最优解
            if all(list(map(lambda f: f.is_integer(), res.x))):
                if self.LOWER_BOUND < -res.fun:
                    self.LOWER_BOUND = -res.fun

                if self.opt_val is None or self.opt_val < -res.fun:
                    self.opt_val = -res.fun
                    self.opt_x = res.x

                continue

            # 进行分枝
            else:
                # 寻找 x 中第一个不是整数的，取其下标 idx
                idx = 0
                for i, x in enumerate(res.x):
                    if not x.is_integer():
                        break
                    idx += 1

                # 构建新的约束条件（分割
                new_con1 = np.zeros(A_ub.shape[1])
                new_con1[idx] = -1
                new_con2 = np.zeros(A_ub.shape[1])
                new_con2[idx] = 1
                new_A_ub1 = np.insert(A_ub, A_ub.shape[0], new_con1, axis=0)
                new_A_ub2 = np.insert(A_ub, A_ub.shape[0], new_con2, axis=0)
                new_b_ub1 = np.insert(
                    b_ub, b_ub.shape[0], -math.ceil(res.x[idx]), axis=0)
                new_b_ub2 = np.insert(
                    b_ub, b_ub.shape[0], math.floor(res.x[idx]), axis=0)

                # 将新约束条件加入队列，先加最优值大的那一支
                r1 = linprog(self.c, new_A_ub1, new_b_ub1, self.A_eq,
                             self.b_eq, self.bounds)
                r2 = linprog(self.c, new_A_ub2, new_b_ub2, self.A_eq,
                             self.b_eq, self.bounds)
                if not r1.success and r2.success:
                    self.Q.put((r2, new_A_ub2, new_b_ub2))
                elif not r2.success and r1.success:
                    self.Q.put((r1, new_A_ub1, new_b_ub1))
                elif r1.success and r2.success:
                    if -r1.fun > -r2.fun:
                        self.Q.put((r1, new_A_ub1, new_b_ub1))
                        self.Q.put((r2, new_A_ub2, new_b_ub2))
                    else:
                        self.Q.put((r2, new_A_ub2, new_b_ub2))
                        self.Q.put((r1, new_A_ub1, new_b_ub1))


def test():
    """ 此测试的真实最优解为 [4, 2] """
    c = np.array([40, 90])
    A = np.array([[9, 7], [7, 20]])
    b = np.array([56, 70])
    Aeq = None
    beq = None
    bounds = [(0, None), (0, None)]

    solver = ILP(c, A, b, Aeq, beq, bounds)
    solver.solve()

    print("Test 's result:", solver.opt_val, solver.opt_x)
    print("Test 's true optimal x: [4, 2]\n")


if __name__ == '__main__':
    test()

代码摘自：整数线性规划实现（lingo，python分枝界定法）

（三）割平面法

1. 算法概述

我们以一道题为例来简单阐述一下割平面法的整体思路。
$\max Z=x_1+x_2\\ s.t.\left\{ \begin{aligned} -x_1+x_2&\leq1\\ 3x_1+x_2&\leq4\\ x_1,x_2 &\geq 0,且为整数 \end{aligned} \right.$
为了方便理解，我们直接放图：

在这里插入图片描述

左侧图中蓝色部分即为原整数规划问题的松弛问题的可行域，由图解法，我们得到红色线即为最优解。

随后，我们将右图中蓝色阴影区域切割，即添加约束条件： $x < 1$ ，因为 $x_1$ ， $x_2$ 都为整数，所以我们去掉的这部分并没有减少原整数规划可行解的数量，最优解自然也不会变化。

接下来我们继续对新的式子做线性规划求解，如果解仍然不满足条件，那么我们可以继续进行切割，直到算出最优解。

上面的可能不够严谨，仅代表个人理解，详细的割平面法还请大家移步：【整数规划(六)】割平面法

2. 代码实现

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

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

（四）动态规划

动态规划问题我会再单独发出来一篇文章，这里就省略了。

三、Matlab求解

官方文档：intlinprog

完整函数：

[x,fval,exitflag]= intlinprog(f,intcon,A,b,Aeq,beq,lb,ub)

整体来讲，intlinprog 与 linprog 使用方法保持一致，只是多了参数 intcon。这个参数表示有那些变量是整数，如：intcon = [1,2,7] 表示 $x_1$ 、 $x_2$ 和 $x_7$ 仅取整数值。

示例代码:

f_13=[-1,-1];
ic_13=[1,2];
A_13=[-4,2;4,2;0,-2];
b_13=[-1;11;-1];
lb_13=zeros(2,1);
[x_13,fval_13,flag_13]=intlinprog(f_13,ic_13,A_13,b_13,[],[],lb_13,[])