浅谈Python数学建模之整数规划-巨人网络通讯

浅谈Python数学建模之整数规划

一、从线性规划到整数规划

1.1、为什么会有整数规划？

线性规划问题的最优解可能是分数或小数。整数规划是指变量的取值只能是整数的规划。

这在实际问题中很常见，例如车间人数、设备台数、行驶次数，这些变量显然必须取整数解。

根据对变量的不同情况，整数规划又可以分为：

完全整数规划，全部变量都要求是整数；
混合整数规划，部分变量要求是整数；
0-1整数规划，变量的取值只能是 0 或 1；
混合0-1规划，部分变量的取值只能是 0 或 1。

0-1整数规划是非常重要也非常特殊的整数规划，需要在另外的文章进行讨论。

1.2、四舍五入就能得到整数解吗？

整数规划问题与线性规划问题的区别只是增加了整数约束。这看上去好像只要把线性规划得到的非整数解舍入化整，就可以得到整数解，并不是多么复杂的问题。

但是问题并没有这么简单。化整后的解不仅不一定是最优解，甚至不一定是可行解的——线性规划的最优解，取整后可能就不满足约束条件了。

那么，不要按四舍五入取整，而是向满足约束条件的方向取整，是不是就可以呢？这是很好的想法，通常这样可以获得可行解，但却不一定是最优解了。

因此，整数规划问题比线性规划复杂的多，以至于至今还没有通用的多项式解法，也就是说算法复杂度与问题规模成指数关系（NP问题）。还没有意识到与问题规模指数关系意味着什么吗？就是那个在象棋棋盘上放麦子，每格比前一格加倍的故事。

问题区别一点点，难度却相差千万里。小白与学霸，差距其实并不大。

二、整数规划的求解方法

2.1、分支定界法（Branch and bound）

分支定界法的基本思想是把原问题（整数规划问题）转换为一个个线性规划问题来处理，并在求解这些线性规划问题的过程中不断追踪原问题的上界（最优可行解）和下界（最优线性松弛解）。

分支定界法把全部可行解空间反复地分割为越来越小的子集，称为分枝；并且对每个子集内的解集计算一个目标上界，称为定界。每次分枝后，对于超出已知可行解集目标值的那些子集不再进一步分枝，就可以删减很多子集，这称为剪枝。

数学课代表的说法是：设有最大化的整数规划问题 A，先解与之相应的线性规划问题 B，若 B 的最优解不符合 A 的整数条件，则 B 的最优目标函数必是 A 的最优目标函数 z 的上界，记为 z2，而 A 的任意可行解的目标函数值将是 z 的一个下界 z1。分支定界法就是将 B 的可行域分成子区域（分支）的方法，逐步减小 z2 和增大 z1，最终求到 z*。

分支定界法是一个迭代算法，随着迭代过程不断更新上界和下界，直到上界和下界非常接近时结束。通常设置 Gap 0.1%，就可把当前的最优可行解近似为问题的全局最优解了。因此，分支定界法的“收敛” 不是分析意义上的而是算法意义上的，优化结果是近似解而不是精确解。

分支定界法不用区分完全整数规划与混合整数规划，算法便于实现，但计算量比较大。

2.2、割平面法（Cutting plane）

割平面法的基本思路是先求解普通线性规划问题的最优解，再对非整数解添加约束条件使可行域缩小，如此反复求解添加了约束条件的普通线性规划问题，直到得到整数解。

也就是说，先不考虑整数约束条件，直接求解松弛问题的最优解，如果满足整数条件就结束了，如果不满足整数条件，就在此非整数解的基础上增加新的约束条件重新求解。这个新增加的约束条件称为割平面，对松弛问题的可行域割一刀，割去松弛问题的部分非整数解。经过有限次的反复切割，必定可在缩小的可行域的一个整数极点上达到整数规划问题的最优解。

割平面法的计算量比较小，但对问题的结构及求解的要求较高，算法比较复杂。

2.3、整数规划的编程方案

在各种算法的介绍和评价中，有时会说“算法比较简单，编程比较容易”。对此小白千万不要当真。不论分支定界法还是割平面法，小白不要说自己按照算法步骤一步步编程实现，就是给你现成的程序估计你也看不懂的。这很正常，就算大神也没几个人能看懂哪怕是自己写出来的算法。

但是如果给你程序也不会使用，那就是问题了。不幸的是，这是数学建模学习和参赛中经常遇到的问题：有了调试好的程序，例程运行结果也正常，但换个问题仍然不会使用。

这并不是你的错。程序有漏洞，接口不标准，文档对不上，教程说不清，这就是你所拿到的例程。你的错误，是选择了这样的例程，或者说选择了这样的编程方案。

这也是本系列教程希望解决的问题。就拿线性规划、整数规划来说，算法还不是很复杂，第三方软件包也很丰富。但是，Scipy 只能求解线性规划，不能求解整数规划，如果选择 Scipy 做线性规划，那在学整数规划时就要再学另一种工具包，二者的模型描述、函数定义、参数设置肯定也是不同的。接下来遇到非线性规划问题再学一种软件包，最后别说熟练掌握算法函数，连什么时候该用哪个工具包都搞晕了。

闲话少说，我们还是用上节求解线性规划问题的 PuLP 工具包。

三、PuLP 求解整数规划问题

我们不仅继续用 PuLP 工具包，而且解题过程和编程步骤也与求解线性规划问题完全一致。

下面我们以一个简单的数学模型练习，来讲解整个解题过程，而不仅给出例程。

3.1、案例问题描述

例题 1：

某厂生产甲乙两种饮料，每百箱甲饮料需用原料 6千克、工人 10名，获利 10万元；每百箱乙饮料需用原料 5千克、工人 20名，获利 9万元。

今工厂共有原料 60千克、工人 150名，又由于其他条件所限甲饮料产量不超过8百箱。

问题 1：问如何安排生产计划，即两种饮料各生产多少使获利最大？

问题 2：若投资0.8万元可增加原料1千克，是否应作这项投资？投资多少合理？

问题 3：若不允许散箱（按整百箱生产），如何安排生产计划，即两种饮料各生产多少使获利最大？

问题 4：若不允许散箱（按整百箱生产），若投资0.8万元可增加原料1千克，是否应作这项投资？投资多少合理？

3.2、建模过程分析

线性规划和整数规划类的问题的建模和求解，通常可以按问题定义、模型构建、模型求解的步骤进行。

3.2.1、问题定义

问题定义，确定决策变量、目标函数和约束条件。

1.决策变量是问题中可以在一定范围内进行变化而获得不同结果的变量。

对于问题 1，问题描述中说的很明确，希望通过改变甲、乙两种饮料的产量使总利润最大，甲、乙两种饮料的产量就是决策变量。

对于问题 2 则要注意，如果只看前一句，就是比较问题 1 与问题 2 的利润，还是把甲、乙两种饮料的产量作为决策变量。但要回答后一句“投资多少合理”，这就出现了一个新的变量“投资额”，因此对问题 2 要建立 3个决策变量：甲产量、乙产量和投资额。

2.目标函数是决策变量的函数，我们希望通过改变决策变量的值而获得目标函数的最大值或最小值，通常是总成本（最小）、总利润（最大）、总时间（最短）。

对于本案例，每个问题都是希望获得最大利润，目标函数都是总利润，问题是求目标函数即总利润的最大值。

3.约束条件是决策变量所要满足的限制条件。

约束条件 3 种情况：

一是不等式约束，例如题目指出共有原料 60千克、工人 150名，因此生产计划所用的原料、工人的需求不能大于题目中数值。

二是等式约束，本题没有等式约束条件。

三是决策变量取值范围的约束。

通常，题目隐含着决策变量大于等于 0 的条件，例如工人人数、原料数量都要大于等于 0。

另外，如果能通过分析前面的等式约束或不等式约束，得出决策变量的上限，将会极大的提高问题求解的速度和性能。后文将对此举例说明。

3.2.2、模型构建

模型构建，由问题描述建立数学方程，并转化为标准形式的数学模型。

对于问题 1，目标函数是生产甲、乙两种饮料的总利润，约束条件是原料总量、工人总数的约束，而且原料、工人都要大于等于 0。

进一步分析决策变量取值范围的约束条件，由原料数量、工人数量的不等式约束可以推出：

对于问题 2，可以通过增加投资来获得更多的原料，投资额是一个新的变量。要注意的是，此时目标函数虽然也是生产两种饮料的总利润，但总利润不等于总收入，而是总收入减去总成本，在本例中就是要减去购买原料的投资。

对于问题 3 和问题 4，区别只是不允许散箱，明确提出了决策变量 x1、x2 的取值要取整数值，所以是整数规划问题。
需要注意的是，问题 4 中对增加的投资额即购买的原料数量并没有整数限制，因此 x1、x2 的取值范围是正整数，但 x3 的取值范围是正数，这是一个混合整数规划问题。
还要说明的是，对于问题 1 和问题 2，虽然题目中没有明确要求生产甲、乙饮料的工人人数为整数，但是人数也不可能是小数的，那么这是不是也是整数规划问题呢？
如果你能提出这个问题，那么恭喜你，你已经从小白升级为菜鸟了。
我的理解是，这个问题怎么说都可以。如果要简化问题，使用线性规划模型，最好在问题假设中说一句，假设甲乙饮料在同一车间先后生产，只要允许甲乙饮料散箱生产，即使根据产量所求出的工人数是小数，也可以解释的通。如果你掌握了整数规划问题的求解，那就先按线性规划建模，再补充讨论工人人数也必须是整数的条件，按整数规划建模求解，这就是妥妥的获奖论文了。

3.2.3、模型求解

模型求解，用标准模型的优化算法对模型求解，得到优化结果。

在线性规划问题中已经讲过使用 PuLP 的求解步骤：

（0）导入 PuLP库函数

import pulp

（1）定义一个规划问题

ProbLP1 = pulp.LpProblem("ProbLP1", sense=pulp.LpMaximize)    # 定义问题 1，求最大值

pulp.LpProblem 用来定义问题的构造函数。"ProbLP1"是用户定义的问题名。
参数 sense 指定问题求目标函数的最小值/最大值。本例求最大值，选择 “pulp.LpMaximize” 。

（2）定义决策变量

对于问题 1：

x1 = pulp.LpVariable('x1', lowBound=0, upBound=15, cat='Continuous')  # 定义 x1
x2 = pulp.LpVariable('x2', lowBound=0, upBound=7.5, cat='Continuous')  # 定义 x2

pulp.LpVariable 用来定义决策变量的函数。'x1'、'x2' 是用户定义的变量名。
参数 lowBound、upBound 用来设定决策变量的下界、上界；可以不定义下界/上界，默认的下界/上界是负无穷/正无穷。本例中 x1、x2 的取值区间分别为 [0,15]、[0,7.5]。
参数 cat 用来设定变量类型，可选参数值：'Continuous' 表示连续变量（默认值）、' Integer ' 表示离散变量（用于整数规划问题）、' Binary ' 表示0/1变量（用于0/1规划问题）。

对于问题 3，甲乙饮料产量 x1、x2 必须取整数，是整数规划问题，因此要设置变量类型为离散变量（整数变量）：

x1 = pulp.LpVariable('x1', lowBound=0, upBound=15, cat='Integer')  # 定义 x1，变量类型：整数
x2 = pulp.LpVariable('x2', lowBound=0, upBound=7.5, cat='Integer')  # 定义 x2，变量类型：整数

（3）添加目标函数

ProbLP1 += (10*x1 + 9*x2)  # 设置目标函数 f(x)

添加目标函数使用 "问题名 += 目标函数式" 格式。

（4）添加约束条件

ProbLP1 += (6*x1 + 5*x2 = 60)  # 不等式约束
ProbLP1 += (10*x1 + 20*x2 = 150)  # 不等式约束

添加约束条件使用 "问题名 += 约束条件表达式" 格式。
约束条件可以是等式约束或不等式约束，不等式约束可以是小于等于或大于等于，分别使用关键字">="、"="和"=="。

（5）求解

ProbLP1.solve()
print(ProbLP1.name)  # 输出求解状态
print("Status:", pulp.LpStatus[ProbLP1.status])  # 输出求解状态
for v in ProbLP1.variables():
    print(v.name, "=", v.varValue)  # 输出每个变量的最优值
print("F1(x) =", pulp.value(ProbLP1.objective))  # 输出最优解的目标函数值

solve() 是求解函数，可以对求解器、求解精度进行设置。
PuLP默认采用 CBC 求解器来求解优化问题，也可以调用其它的优化器来求解，但需要另外安装。　

3.3、Python 例程

import pulp      # 导入 pulp 库

# 主程序
def main():

    # 模型参数设置
    """
    问题描述：
        某厂生产甲乙两种饮料，每百箱甲饮料需用原料6千克、工人10名，获利10万元；每百箱乙饮料需用原料5千克、工人20名，获利9万元。
        今工厂共有原料60千克、工人150名，又由于其他条件所限甲饮料产量不超过8百箱。
        （1）问如何安排生产计划，即两种饮料各生产多少使获利最大？
        （2）若投资0.8万元可增加原料1千克，是否应作这项投资？投资多少合理？
        （3）若不允许散箱（按整百箱生产），如何安排生产计划，即两种饮料各生产多少使获利最大？
        （4）若不允许散箱（按整百箱生产），若投资0.8万元可增加原料1千克，是否应作这项投资？投资多少合理？
    """

    # 问题 1：
    """
    问题建模：
        决策变量：
            x1：甲饮料产量（单位：百箱）
            x2：乙饮料产量（单位：百箱）
        目标函数：
            max fx = 10*x1 + 9*x2
        约束条件：
            6*x1 + 5*x2 = 60
            10*x1 + 20*x2 = 150            
            x1, x2 >= 0，x1 = 8
    此外，由 x1,x2>=0 和 10*x1+20*x2=150 可知 0=x2=7.5
    """
    ProbLP1 = pulp.LpProblem("ProbLP1", sense=pulp.LpMaximize)    # 定义问题 1，求最大值
    x1 = pulp.LpVariable('x1', lowBound=0, upBound=8, cat='Continuous')  # 定义 x1
    x2 = pulp.LpVariable('x2', lowBound=0, upBound=7.5, cat='Continuous')  # 定义 x2
    ProbLP1 += (10*x1 + 9*x2)  # 设置目标函数 f(x)
    ProbLP1 += (6*x1 + 5*x2 = 60)  # 不等式约束
    ProbLP1 += (10*x1 + 20*x2 = 150)  # 不等式约束
    ProbLP1.solve()
    print(ProbLP1.name)  # 输出求解状态
    print("Status youcans:", pulp.LpStatus[ProbLP1.status])  # 输出求解状态
    for v in ProbLP1.variables():
        print(v.name, "=", v.varValue)  # 输出每个变量的最优值
    print("F1(x) =", pulp.value(ProbLP1.objective))  # 输出最优解的目标函数值


    # 问题 2：
    """
    问题建模：
        决策变量：
            x1：甲饮料产量（单位：百箱）
            x2：乙饮料产量（单位：百箱）
            x3：增加投资（单位：万元）
        目标函数：
            max fx = 10*x1 + 9*x2 - x3
        约束条件：
            6*x1 + 5*x2 = 60 + x3/0.8
            10*x1 + 20*x2 = 150
            x1, x2, x3 >= 0，x1 = 8
    此外，由 x1,x2>=0 和 10*x1+20*x2=150 可知 0=x2=7.5
    """
    ProbLP2 = pulp.LpProblem("ProbLP2", sense=pulp.LpMaximize)    # 定义问题 2，求最大值
    x1 = pulp.LpVariable('x1', lowBound=0, upBound=8, cat='Continuous')  # 定义 x1
    x2 = pulp.LpVariable('x2', lowBound=0, upBound=7.5, cat='Continuous')  # 定义 x2
    x3 = pulp.LpVariable('x3', lowBound=0, cat='Continuous')  # 定义 x3
    ProbLP2 += (10*x1 + 9*x2 - x3)  # 设置目标函数 f(x)
    ProbLP2 += (6*x1 + 5*x2 - 1.25*x3 = 60)  # 不等式约束
    ProbLP2 += (10*x1 + 20*x2 = 150)  # 不等式约束
    ProbLP2.solve()
    print(ProbLP2.name)  # 输出求解状态
    print("Status  youcans:", pulp.LpStatus[ProbLP2.status])  # 输出求解状态
    for v in ProbLP2.variables():
        print(v.name, "=", v.varValue)  # 输出每个变量的最优值
    print("F2(x) =", pulp.value(ProbLP2.objective))  # 输出最优解的目标函数值

    # 问题 3：整数规划问题
    """
    问题建模：
        决策变量：
            x1：甲饮料产量，正整数（单位：百箱）
            x2：乙饮料产量，正整数（单位：百箱）
        目标函数：
            max fx = 10*x1 + 9*x2
        约束条件：
            6*x1 + 5*x2 = 60
            10*x1 + 20*x2 = 150
            x1, x2 >= 0，x1 = 8，x1, x2 为整数
    此外，由 x1,x2>=0 和 10*x1+20*x2=150 可知 0=x2=7.5
    """
    ProbLP3 = pulp.LpProblem("ProbLP3", sense=pulp.LpMaximize)  # 定义问题 3，求最大值
    print(ProbLP3.name)  # 输出求解状态
    x1 = pulp.LpVariable('x1', lowBound=0, upBound=8, cat='Integer')  # 定义 x1，变量类型：整数
    x2 = pulp.LpVariable('x2', lowBound=0, upBound=7.5, cat='Integer')  # 定义 x2，变量类型：整数
    ProbLP3 += (10 * x1 + 9 * x2)  # 设置目标函数 f(x)
    ProbLP3 += (6 * x1 + 5 * x2 = 60)  # 不等式约束
    ProbLP3 += (10 * x1 + 20 * x2 = 150)  # 不等式约束
    ProbLP3.solve()
    print("Shan Status:", pulp.LpStatus[ProbLP3.status])  # 输出求解状态
    for v in ProbLP3.variables():
        print(v.name, "=", v.varValue)  # 输出每个变量的最优值
    print("F3(x) =", pulp.value(ProbLP3.objective))  # 输出最优解的目标函数值


    # 问题 4：
    """
    问题建模：
        决策变量：
            x1：甲饮料产量，正整数（单位：百箱）
            x2：乙饮料产量，正整数（单位：百箱）
            x3：增加投资（单位：万元）
        目标函数：
            max fx = 10*x1 + 9*x2 - x3
        约束条件：
            6*x1 + 5*x2 = 60 + x3/0.8
            10*x1 + 20*x2 = 150
            x1, x2, x3 >= 0，x1 = 8，x1, x2 为整数
    此外，由 x1,x2>=0 和 10*x1+20*x2=150 可知 0=x2=7.5
    """
    ProbLP4 = pulp.LpProblem("ProbLP4", sense=pulp.LpMaximize)  # 定义问题 4，求最大值
    print(ProbLP4.name)  # 输出求解状态
    x1 = pulp.LpVariable('x1', lowBound=0, upBound=8, cat='Integer')  # 定义 x1，变量类型：整数
    x2 = pulp.LpVariable('x2', lowBound=0, upBound=7, cat='Integer')  # 定义 x2，变量类型：整数
    x3 = pulp.LpVariable('x3', lowBound=0, cat='Continuous')  # 定义 x3
    ProbLP4 += (10*x1 + 9*x2 - x3)  # 设置目标函数 f(x)
    ProbLP4 += (6*x1 + 5*x2 - 1.25*x3 = 60)  # 不等式约束
    ProbLP4 += (10*x1 + 20*x2 = 150)  # 不等式约束
    ProbLP4.solve()
    print("Shan Status:", pulp.LpStatus[ProbLP4.status])  # 输出求解状态
    for v in ProbLP4.variables():
        print(v.name, "=", v.varValue)  # 输出每个变量的最优值
    print("F4(x) =", pulp.value(ProbLP4.objective))  # 输出最优解的目标函数值

    return

if __name__ == '__main__':  # Copyright 2021 YouCans, XUPT
    main()

3.4、Python 例程运行结果

Welcome to the CBC MILP Solver

Version: 2.9.0

Build Date: Feb 12 2015

ProbLP1

Status: Optimal

x1 = 6.4285714

x2 = 4.2857143

F1(x) = 102.8571427

ProbLP2

Status: Optimal

x1 = 8.0

x2 = 3.5

x3 = 4.4

F2(x) = 107.1

ProbLP3

Result - Optimal solution found

Status Shan: Optimal

Status: Optimal

x1 = 8.0

x2 = 2.0

F3(x) = 98.0

ProbLP4

Result - Optimal solution found

Status: Optimal

x1 = 8.0

x2 = 3.0

x3 = 2.4

F4(x) = 104.6