Python小白的數學建模課-04.整數規劃

youcans發表於2021-06-03
Python

整數規劃與線性規劃的差別只是變數的整數約束。

問題區別一點點，難度相差千萬裡。

選擇簡單通用的程式設計方案，讓求解器去處理吧。

『Python小白的數學建模課 @ Youcans』帶你從數模小白成為國賽達人。


1. 從線性規劃到整數規劃

1.1 為什麼會有整數規劃？

線性規劃問題的最優解可能是分數或小數。整數規劃是指變數的取值只能是整數的規劃。

這在實際問題中很常見，例如車間人數、裝置臺數、行駛次數，這些變數顯然必須取整數解。

整數規劃並不一定是線性規劃問題的變數取整限制，對於二次規劃、非線性規劃問題也有變數取整限制而引出的整數規劃。但在數學建模問題中所說的整數規劃，通常是指整數線性規劃。

根據對變數的不同情況，整數規劃又可以分為：

  • 完全整數規劃，全部變數都要求是整數；
  • 混合整數規劃，部分變數要求是整數；
  • 0-1整數規劃，變數的取值只能是 0 或 1；
  • 混合0-1規劃，部分變數的取值只能是 0 或 1。

0-1整數規劃 是非常重要也非常特殊的整數規劃，需要在另外的文章進行討論。


1.2 四捨五入就能得到整數解嗎？

整數規劃問題與線性規劃問題的區別只是增加了整數約束。這看上去好像只要把線性規劃得到的非整數解舍入化整，就可以得到整數解，並不是多麼複雜的問題。

但是問題並沒有這麼簡單。化整後的解不僅不一定是最優解，甚至不一定是可行解的——線性規劃的最優解，取整後可能就不滿足約束條件了。

那麼，不要按四捨五入取整，而是向滿足約束條件的方向取整，是不是就可以呢？這是很好的想法，通常這樣可以獲得可行解，但卻不一定是最優解了。

Python小白的數學建模課-04.整數規劃

因此，整數規劃問題比線性規劃複雜的多，以至於至今還沒有通用的多項式解法，也就是說演算法複雜度與問題規模成指數關係（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，區別只是不允許散箱，明確提出了決策變數 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 例程

# 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小白的數學建模課-A1.國賽賽題型別分析
Python小白的數學建模課-A2.2021年數維杯C題探討
Python數模筆記-PuLP庫
Python數模筆記-StatsModels統計迴歸
Python數模筆記-Sklearn
Python數模筆記-NetworkX
Python數模筆記-模擬退火演算法


相關文章