Python小白的数学建模课-04.整数规划

时间  2021-07-23
标签 html python 算法 编程 函数 工具 性能 学习 优化 栏目 Python 繁體版
原文   原文链接

整数规划与线性规划的差异只是变量的整数约束。html

问题区别一点点，难度相差千万里。python

选择简单通用的编程方案，让求解器去处理吧。算法

『Python小白的数学建模课 @ Youcans』带你从数模小白成为国赛达人。编程


1. 从线性规划到整数规划

1.1 为何会有整数规划？

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

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

整数规划并不必定是线性规划问题的变量取整限制，对于二次规划、非线性规划问题也有变量取整限制而引出的整数规划。但在数学建模问题中所说的整数规划，一般是指整数线性规划。性能

根据对变量的不一样状况，整数规划又能够分为：学习

  • 彻底整数规划，所有变量都要求是整数；
  • 混合整数规划，部分变量要求是整数；
  • 0-1整数规划，变量的取值只能是 0 或 1；
  • 混合0-1规划，部分变量的取值只能是 0 或 1。

0-1整数规划 是很是重要也很是特殊的整数规划，须要在另外的文章进行讨论。优化


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

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

可是问题并无这么简单。化整后的解不只不必定是最优解，甚至不必定是可行解的——线性规划的最优解，取整后可能就不知足约束条件了。

那么，不要按四舍五入取整，而是向知足约束条件的方向取整，是否是就能够呢？这是很好的想法，一般这样能够得到可行解，但却不必定是最优解了。

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

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

欢迎关注 『Python小白的数学建模课 @ Youcans』，每周更新数模笔记
Python小白的数学建模课-01.新手必读
Python小白的数学建模课-02.数据导入
Python小白的数学建模课-03.线性规划
Python小白的数学建模课-04.整数规划
Python数模笔记-PuLP库
Python数模笔记-StatsModels统计回归
Python数模笔记-Sklearn
Python数模笔记-NetworkX
Python数模笔记-模拟退火算法


2. 整数规划的求解方法

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 工具包。


3. 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。

\[max\;f(x) = 10*x_1 + 9*x_2\\ s.t.:\begin{cases} 6*x_1 + 5*x_2 \leq 60\\ 10*x_1 + 20*x_2 \leq 150\\ 0 \leq x_1 \leq 8\\ x_2 \geq 0 \end{cases} \]

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

\[x_1 \leq 15\\ x_2 \leq 7.5 \]

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

\[max\;f(x) = 10*x_1 + 9*x_2 - x_3\\ s.t.:\begin{cases} 6*x_1 + 5*x_2 \leq 60 + x_3/0.8\\ 10*x_1 + 20*x_2 \leq 150\\ 0 \leq x_1 \leq 15\\ 0 \leq x_2 \leq 7.5\\ x_3 \geq 0 \end{cases} \]

对于问题 3 和问题 4，区别只是不容许散箱，明确提出了决策变量 x一、x2 的取值要取整数值，因此是整数规划问题。
须要注意的是，问题 4 中对增长的投资额即购买的原料数量并无整数限制，所以 x一、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 用来设定决策变量的下界、上界；能够不定义下界/上界，默认的下界/上界是负无穷/正无穷。本例中 x一、x2 的取值区间分别为 [0,15]、[0,7.5]。
参数 cat 用来设定变量类型，可选参数值：'Continuous' 表示连续变量（默认值）、' Integer ' 表示离散变量（用于整数规划问题）、' Binary ' 表示0/1变量（用于0/1规划问题）。

对于问题 3， 甲乙饮料产量 x一、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 例程

# mathmodel05_v1.py
# Demo05 of mathematical modeling algorithm
# Solving integer programming with PuLP.
# Copyright 2021 Youcans, XUPT
# Crated：2021-05-31
# Python小白的数学建模课 @ Youcans

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()  # Python小白的数学建模课 @ Youcans

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

【本节完】


版权说明：

欢迎关注『Python小白的数学建模课 @ Youcans』 原创做品

原创做品，转载必须标注原文连接:(https://www.cnblogs.com/youcans/p/14844841.html)。

Copyright 2021 Youcans, XUPT

Crated：2021-05-31


欢迎关注 『Python小白的数学建模课 @ Youcans』，每周更新数模笔记
Python小白的数学建模课-01.新手必读
Python小白的数学建模课-02.数据导入
Python小白的数学建模课-03.线性规划
Python小白的数学建模课-04.整数规划
Python小白的数学建模课-05.0-1规划
Python小白的数学建模课-06.固定费用问题
Python小白的数学建模课-07.选址问题
Python小白的数学建模课-09.微分方程模型
Python数模笔记-PuLP库
Python数模笔记-StatsModels统计回归
Python数模笔记-Sklearn
Python数模笔记-NetworkX
Python数模笔记-模拟退火算法****

相关文章
相关标签/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公众号
   欢迎关注本站公众号,获取更多信息
联系我们 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程