Python求解線性規劃——PuLP使用教程

Only(AR)發表於2022-04-26

簡潔是智慧的靈魂，冗長是膚淺的藻飾。——莎士比亞《哈姆雷特》

1 PuLP 庫的安裝

如果您使用的是 Anaconda^[1] 的話（事實上我也更推薦這樣做），需要先啟用你想要安裝的虛擬環境，之後在 Prompt 輸入

pip install pulp

不出意外的話等一會就安裝完畢。

2 線性規劃簡介

想必大家能點開這篇文章一定都知道線性規劃是什麼意思吧……那麼我用兩個例子再簡單說一下。

2.1 線性規劃

2.1.1 題目描述^[2]

若變數 \(x, y\) 滿足約束條件：

\[\left\{ \begin{aligned} & 2x + 3y - 6\geq 0\\ & x + y - 3 \leq 0\\ & y - 2 \leq 0 \end{aligned} \right. \]

求 \(z = 3x + y\) 的最大值。

2.1.2 基本概念

首先，我們要認清在這道題中，\(x\) 和 \(y\) 是可以變的，所以把它們叫做決策變數。三個不等式叫做約束條件，即 \(x\) 和 \(y\) 必須同時滿足這三個不等式。我們若畫出圖來：

其中不滿足約束條件的區域被我標上了顏色，所以 \(x, y\) 可以取得值只能在純白區域內，這一片區域稱作可行域。

再看最後的我們的目標：求 \(z = x + 3y\) 的最大值。

於是 \(z=x+3y\) 就被稱作目標函式，我們的工作就是求這個目標函式的最大值。

整個問題描述為：

\[\begin{eqnarray*} &\max &z = x+3y \tag{1}\\ &\mathrm{s.t.} & \quad 2x + 3y - 6 \geq0 \tag{2}\\ & & \quad x + 3y - 3 \leq 0 \tag{3}\\ & & \quad y - 2 \leq 0 \tag{4} \end{eqnarray*} \]

然後怎麼算？別急我們再看一個例子。

2.2 整數規劃

2.2.1 題目描述^[3]

汽車廠生產小、中、大三種型別的汽車，已知各型別每輛車對鋼材、勞動時間的需求以及利潤如下表所示。要求每月的鋼材消耗不超過 600 t，總勞動時間不超過 60 000 h。試指定生產計劃使得工廠每月的利潤最大。

	小型車	中型車	大型車
鋼材 / t	1.5	3	5
勞動時間 / h	280	250	400
利潤 / 萬元	2	3	4

2.2.2 解題思路

首先，設三個決策變數，用 \(x_1, x_2, x_3\) 分別表示生產小型車、中型車、大型車的數量，但是注意要滿足：

車的數量只能是整數；
車的數量大於等於 0。

其他約束條件看題直接列：

\[\left\{\begin{aligned} & 1.5 x_1 + 3 x_2 + 5 x_3 \leq 600\\ & 280 x_1 + 250 x_2 + 400 x_2 \leq 60000 \end{aligned}\right. \]

最後寫出目標函式：

\[z = 2x_1 + 3x_2 + 4x_3 \]

綜合起來整個問題描述為：

\[\begin{eqnarray*} &\max & z = 2x_1 + 3x_2 + 4x_3 \tag{1}\\ &\mathrm{s.t.} & 1.5 x_1 + 3 x_2 + 5 x_3 \leq 600\tag{2}\\ & & 280 x_1 + 250 x_2 + 400 x_2 \leq 60000\tag{3}\\ & & x_1, x_2, x_3 \geq 0\tag{4}\\ & & x_1, x_2, x_3 均為整數\tag{5} \end{eqnarray*} \]

另外可以看出這個題由於涉及到三個決策變數，可行域是相當抽象的，這裡就不畫了 hhh~

3 求解過程

首先在最前面引入所需的pulp工具庫：

import pulp as pl

這句話是引入 pulp 庫並簡寫為 pl，一個 python 庫只有在開始 import 了之後才能在後面使用。這樣後面凡是用到 pulp 的功能都要寫成 pl.xxx。

接下來是以下幾個步驟：

定義模型
定義決策變數
新增約束條件
新增目標函式
模型求解
列印結果

3.1 定義模型

# Define the model
model = pl.LpProblem(name="My-Model", sense=pl.LpMaximize)

這個操作是使用 pl.LpProblem 建立了一個模型並賦值給變數 model，接收兩個引數：

name：模型的名字，隨便起一個；
sense：模型的型別，pl.LpMinimize是求目標函式的最小值，pl.LpMaximize 是求最大值

3.2 定義決策變數

# Define the decision variables
x = pl.LpVariable(name='x')
y = pl.LpVariable(name='y')

如果你的變數比較少的話可以簡單這麼寫。這個意思是定義了兩個浮點數變數，取值範圍是整個實數域。注意等號左邊的變數才是你在之後的計算式中使用的符號，而引數 name 只有在最後列印結果的時候才會被列印出來。另外如果你對變數有其他要求的話可以新增以下引數：

lowBound：變數的最小取值（不寫的話預設負無窮）；
upBound：變數的最大取值（預設正無窮）；
cat：變數的型別，有 pl.Binary 邏輯變數、pl.Integer 整數、pl.Continuous 實數（預設值）；

如果你的變數比較多而不得不用 1, 2, 3…… 來編號，可以採用類似這樣的寫法：

# Define the decision variables
x = {i: pl.LpVariable(name=f"x{i}", lowBound=0, cat=pl.LpInteger) for i in range(1, 9)}

這是一次定義 8 個變數並儲存在一個類似陣列的結構中，變數都是正整數，分別用 x[1], x[2], ..., x[8] 表示，依次命名為 x1, x2,..., x8。

注意 range(left, right) 表示的區間是左閉右開。

3.3 新增約束條件

# Add constraints
model += (2 * x + 3 * y - 6 >= 0, "constrain_1")
model += (x + 3 * y - 3 == 0, "constrain_2")

沒錯！如你所見就是這麼簡單，括號裡第一個變數就是你的約束不等式或等式，第二個變數是你的自定義的約束名（可以起一個有意義的名字，當然也可以省略）。

由於一些比較數學的原因，約束條件裡是不能使用大於號“>”或小於號“<”的。

如果你像前面一樣把變數定義在了陣列中，那麼可以直接用方括號呼叫：

model += (2 * x[1] + 3 * x[2] - 6 >= 0)

3.4 新增目標函式

# Set the objective
model += x + 3 * y

與前面新增約束條件不同，新增目標函式這一步不用加最外層的括號。

3.5 模型求解

# Solve the optimization problem
status = model.solve()

就寫這一句話，呼叫 model 的 solve() 方法，並把結果儲存在 status 中。

3.4 列印結果

# Get the results
print(f"status: {model.status}, {pl.LpStatus[model.status]}")
print(f"objective: {model.objective.value()}")

for var in model.variables():
    print(f"{var.name}: {var.value()}")

for name, constraint in model.constraints.items():
    print(f"{name}: {constraint.value()}")

然後你就能看到模型求解的結果了。

4 示例程式碼

4.1 高考題程式碼

首先解決一下 3.1 的高考題：

import pulp as pl

# 定義一個模型，命名為 "Model_3.1"，求最大值
model = pl.LpProblem(name="Model_3.1", sense=pl.LpMaximize)

# 定義兩個決策變數，取值為整個實數域
x = pl.LpVariable(name='x')
y = pl.LpVariable(name='y')

# 新增三個約束條件
model += (2 * x + 3 * y - 6 >= 0)
model += (x + y - 3 <= 0)
model += (y - 2 <= 0)

# 目標函式
model += x + 3 * y

# 求解
status = model.solve()

# 列印結果
print(f"status: {model.status}, {pl.LpStatus[model.status]}")
print(f"objective: {model.objective.value()}")

for var in model.variables():
    print(f"{var.name}: {var.value()}")

for name, constraint in model.constraints.items():
    print(f"{name}: {constraint.value()}")

檢視結果的最後幾行：

status: 1, Optimal
objective: 7.0
x: 1.0
y: 2.0
_C1: 2.0
_C2: 0.0
_C3: 0.0

最大值是 \(7.0\)，在 \(x=1.0, y=2.0\) 時取到。

4.2 汽車廠程式碼

import pulp as pl

# 定義一個模型，命名為 "Model_3.2"，求最大值
model = pl.LpProblem(name="Model_3.2", sense=pl.LpMaximize)

# 定義三個決策變數，取值正整數
x = {i: pl.LpVariable(name=f"x{i}", lowBound=0, cat=pl.LpInteger) for i in range(1, 4)}

# 新增約束條件
model += (1.5 * x[1] + 3 * x[2] + 5 * x[3] <= 600)
model += (280 * x[1] + 250 * x[2] + 400 * x[3] <= 60000)

# 目標函式
model += 2 * x[1] + 3 * x[2] + 4 * x[3]

# 求解
status = model.solve()

# 列印結果
print(f"status: {model.status}, {pl.LpStatus[model.status]}")
print(f"objective: {model.objective.value()}")

for var in model.variables():
    print(f"{var.name}: {var.value()}")

for name, constraint in model.constraints.items():
    print(f"{name}: {constraint.value()}")