線性規劃是解決最佳化問題的有效工具，廣泛應用於生產規劃、資源分配等領域。本文首先介紹一個多期混合模型，以肥皂生產為例，示範如何利用線性規劃最小化成本並滿足產品需求和品質要求。模型中考慮了不同油料的購買、混合和儲存，以及產品的品質限制。接著，文章轉向機器學習領域，探討如何使用線性模型解決二元分類別問題，並以細胞分類別為例說明模型的建立和程式碼實作。最後，文章探討如何處理非線性最佳化問題，介紹了分段線性近似方法，並提供程式碼範例示範如何將非線性函式轉換為分段線性函式，利用線性規劃求解器進行求解。
多期混合模型分析與實作
在工業生產中，混合問題是一類別常見的最佳化問題，特別是在需要滿足特定化學或物理特性的產品製造過程中。本篇文章將探討一個多期混合模型的建立與實作，該模型旨在最小化成本的同時滿足產品需求和品質要求。
模型背景
某公司需要根據每月不同的成本和庫存狀況，決定每個月購買和混合不同油料以生產肥皂的數量。該模型的目標是最小化總成本，包括購買成本和儲存成本，同時滿足每月的生產需求和產品品質要求。
模型建立
變數定義
Buy[i][j]：第j月購買的第i種油料數量Blnd[i][j]：第j月用於混合的第i種油料數量Hold[i][j]：第j月末儲存的第i種油料數量Prod[j]：第j月的總生產量
限制條件
- 生產需求滿足：每月的總生產量必須滿足需求。
s.Add(Prod[j] >= D)
2. **庫存平衡**：每種油料在每月的庫存量、購買量、混合量和下月庫存量之間必須保持平衡。
```python
s.Add(Hold[i][j] + Buy[i][j] - Blnd[i][j] == Hold[i][j+1])
- 庫存範圍限制：總庫存量必須在指定的上下限範圍內。
s.Add(sum(Hold[i][j] for i in Oils) >= SL[0]) s.Add(sum(Hold[i][j] for i in Oils) <= SL[1])
4. **產品品質要求**：混合產品的脂肪酸含量百分比必須在指定的範圍內。
```python
s.Add(Acid[k][j] == sum(Blnd[i][j] * Part[i][k] for i in Oils))
s.Add(Acid[k][j] >= Target[0][k] * Prod[j])
s.Add(Acid[k][j] <= Target[1][k] * Prod[j])
目標函式
模型的目標是最小化總成本，包括購買成本和儲存成本。
s.Add(CostP[j] == sum(Buy[i][j] * Cost[i][j] for i in Oils))
s.Add(CostS[j] == sum(Hold[i][j] * SC for i in Oils))
Cost_product = s.Sum(CostP[j] for j in Periods)
Cost_storage = s.Sum(CostS[j] for j in Periods)
s.Minimize(Cost_product + Cost_storage)
程式碼實作
def solve_model(Part, Target, Cost, Inventory, D, SC, SL):
# 建立求解器
s = newSolver('Multi-period soap blending problem')
Oils = range(len(Part))
Periods, Acids = range(len(Cost[0])), range(len(Part[0]))
# 定義變數
Buy = [[s.NumVar(0, D, '') for _ in Periods] for _ in Oils]
Blnd = [[s.NumVar(0, D, '') for _ in Periods] for _ in Oils]
Hold = [[s.NumVar(0, D, '') for _ in Periods] for _ in Oils]
Prod = [s.NumVar(0, D, '') for _ in Periods]
CostP = [s.NumVar(0, D*1000, '') for _ in Periods]
CostS = [s.NumVar(0, D*1000, '') for _ in Periods]
Acid = [[s.NumVar(0, D*D, '') for _ in Periods] for _ in Acids]
# 設定初始庫存
for i in Oils:
s.Add(Hold[i][0] == Inventory[i][0])
# 模型限制條件
for j in Periods:
s.Add(Prod[j] == sum(Blnd[i][j] for i in Oils))
s.Add(Prod[j] >= D)
if j < Periods[-1]:
for i in Oils:
s.Add(Hold[i][j] + Buy[i][j] - Blnd[i][j] == Hold[i][j+1])
s.Add(sum(Hold[i][j] for i in Oils) >= SL[0])
s.Add(sum(Hold[i][j] for i in Oils) <= SL[1])
for k in Acids:
s.Add(Acid[k][j] == sum(Blnd[i][j] * Part[i][k] for i in Oils))
s.Add(Acid[k][j] >= Target[0][k] * Prod[j])
s.Add(Acid[k][j] <= Target[1][k] * Prod[j])
s.Add(CostP[j] == sum(Buy[i][j] * Cost[i][j] for i in Oils))
s.Add(CostS[j] == sum(Hold[i][j] * SC for i in Oils))
# 目標函式
Cost_product = s.Sum(CostP[j] for j in Periods)
Cost_storage = s.Sum(CostS[j] for j in Periods)
s.Minimize(Cost_product + Cost_storage)
# 求解模型
rc = s.Solve()
B, L, H, A = SolVal(Buy), SolVal(Blnd), SolVal(Hold), SolVal(Acid)
CP, CS, P = SolVal(CostP), SolVal(CostS), SolVal(Prod)
return rc, ObjVal(s), B, L, H, P, A, CP, CS
內容解密：
變數定義與初始化：程式碼首先定義了多個變數矩陣，如
Buy、Blnd、Hold等，用於表示不同月份和不同油料的購買量、混合量和庫存量。這些變數的定義是根據問題的需求和限制條件。限制條件的建立：程式碼中透過一系列的
s.Add()函式新增了模型的限制條件，包括生產需求滿足、庫存平衡、庫存範圍限制和產品品質要求等。這些條件共同構成了模型的框架。目標函式的建立：目標函式是最小化總成本，包括購買成本和儲存成本。程式碼中透過計算每月的購買成本和儲存成本，並將其合計得到總成本，最後使用
s.Minimize()函式設定為求解目標。模型求解與結果輸出：程式碼最後呼叫
s.Solve()函式對模型進行求解，並透過一系列的SolVal()函式取得最優解下的變數值，如購買量、混合量、庫存量等。這些結果可以用於指導實際生產和庫存管理。
線性連續模型在分類別問題上的應用
在眾多軟體應用中，分類別是目前最成功的人工智慧技術之一。軟體能夠判斷一封電子郵件是否為垃圾郵件、一個活檢細胞是惡性還是良性，甚至決定一家公司是否應該邀請你面試。這些都是真實存在的應用場景。
2.5 模式分類別
讓我們以一個簡單的二元分類別問題為例，來探討如何使用線性模型進行分類別。假設我們要根據細胞的面積和周長這兩個特徵，將細胞分為惡性或良性。這兩個特徵可以透過顯微鏡下的圖片自動測量。這個過程始於一批已經被專家分為兩類別的細胞資料，這些資料構成了我們軟體的訓練集。
2.5.1 構建模型
假設我們有如圖2-3所示的細胞特徵資料，其中x軸代表周長，y軸代表半徑。從圖中可以看出，這兩類別細胞可以用一條直線分開。我們的任務就是找到這條直線。
代數上，一條直線可以用方程式 $a_1x_1 + a_2x_2 = b$ 表示，其中$a_1$、$a_2$和$b$是固定的係數。在n維空間中，這被稱為超平面，其方程式為 $\sum_{i=1}^{n} a_ix_i = a_0$。
對於一個特定的點$x$，它位於直線的哪一側取決於$a_1x_1 + a_2x_2 < a_0$還是$a_1x_1 + a_2x_2 > a_0$。透過縮放，我們可以簡化任務為找到一個向量$a$，使得對於類別A中的每個點$x’$，有$\sum_{i} a_ix’i \geq a_0 + 1$，而對於類別B中的每個點$x’’$，有$\sum{i} a_ix’’_i \leq a_0 - 1$。
2.5.2 可執行模型
接下來，讓我們將上述模型轉換為可執行的程式碼，如清單2-7所示。該函式接受兩組具有任意數量特徵的資料點，這些資料點已被專家分為A類別和B類別。
def solve_classification(A, B):
n, ma, mb = len(A[0]), len(A), len(B)
s = newSolver('Classification')
ya = [s.NumVar(0, 99, '') for _ in range(ma)]
yb = [s.NumVar(0, 99, '') for _ in range(mb)]
a = [s.NumVar(-99, 99, '') for _ in range(n + 1)]
for i in range(ma):
s.Add(ya[i] >= a[n] + 1 - s.Sum(a[j] * A[i][j] for j in range(n)))
for i in range(mb):
s.Add(yb[i] >= s.Sum(a[j] * B[i][j] for j in range(n)) - a[n] + 1)
Agap = s.Sum(ya[i] for i in range(ma))
Bgap = s.Sum(yb[i] for i in range(mb))
s.Minimize(Agap + Bgap)
rc = s.Solve()
return rc, ObjVal(s), SolVal(a)
內容解密：
- 輸入引數：函式
solve_classification接受兩個引數A和B，分別代表兩類別資料點的集合。 - 變數初始化：根據輸入資料的特徵數量和樣本數量，初始化變數
n、ma和mb。建立一個新的求解器solver例項，並為每個資料點建立非負變數ya和yb，以及超平面係數變數a。 - 約束條件：透過迴圈為每個資料點新增約束條件，確保類別A的點滿足$\sum_{i} a_ix’i \geq a_0 + 1$，而類別B的點滿足$\sum{i} a_ix’’_i \leq a_0 - 1$。
- 目標函式：定義目標函式為最小化
Agap + Bgap，即最小化所有資料點到超平面的總偏差。 - 求解：呼叫
s.Solve()求解最佳化問題，並傳回求解結果、目標函式值和超平面係數。
這個模型的一個特點是，如果最優目標函式值為零，則存在一個超平面能夠正確地將訓練集中的惡性細胞和良性細胞分開。如果目標函式值非零，則意味著資料集不能被超平面分離，需要使用更複雜的技術。
隱藏的線性連續模型
在這一章中，我們將對某些問題進行調整，以揭示其內在結構。重點將放在那些乍看之下似乎不是連續線性模型，但透過一些創意性的修改，可以轉化為這種形式的模型。關鍵在於確保原始問題與修改後的問題之間存在一一對應的關係，以便我們可以從修改後問題的解中檢索出原始問題的解。
為什麼要這樣做？
主要原因是連續線性求解器已經變得非常快速，能夠處理具有數十萬個變數和約束的模型。因此，如果一個問題可以建模為這種形式，那麼在實際應用中，能夠解決的問題規模幾乎是沒有限制的。然而，對於更複雜的模型，情況並非如此。事實上，我們可以寫出一些只有幾十個變數的模型，但目前的求解器無法在合理的時間內解決它們。
主要挑戰
本章遇到的主要障礙是非線性問題，但有一個有利的限制，即函式被認為是凸的。一個凸函式是指在其定義域內，函式影像上的任意兩點之間的線段都位於函式影像的上方或影像上。在一維情況下，代數上，如果函式 $f$ 在點 $x_0$ 處是凸的，則滿足以下條件： [ f(x_0 + h) \geq f(x_0) + f’(x_0)h ]
分段線性函式
我們考慮的是分段線性函式。在傳統的術語中，它們被稱為分段線性函式。雖然我們迄今為止使用的線性規劃求解器（GLPK、GLOP、CLP）無法直接處理它們，但透過一些編碼工作，我們可以將它們轉化為所有求解器都能處理的標準形式。
示例：運輸成本函式
讓我們考慮一個分段函式，定義如下： [ f(x) = \begin{cases} C_1x & 0 \leq x \leq B_1 \ C_1B_1 + C_2(x - B_1) & B_1 \leq x \leq B_2 \ C_1B_1 + C_2(B_2 - B_1) + C_3(x - B_2) & B_2 \leq x \leq B_3 \ \vdots & \end{cases} ] 這個函式可以被視為一個運輸成本函式，隨著運輸量的增加，每單位的成本也會增加。
建立模型
在這個問題中，我們需要決定的只是生產量。我們可以定義一個決策變數，其範圍從 0 到表格中的最後一個數量： [ x \in [0, B_n] ] 為了知道 $x$ 落在哪個區間內，我們引入了額外的變數，每個斷點對應一個變數。假設我們有 $n$ 個區間，我們在概念上將這些變數視為區間邊界上的權重，告訴我們 $x$ 在區間中的位置。我們希望最多有兩個連續的變數非零，並且它們的和為 1。
約束條件
我們將強制 $\delta$ 的總和為 1，並且，由於問題的凸結構，最多兩個相鄰的 $\delta$ 非零。這將告訴我們 $x$ 位於哪個區間以及在區間中的位置。 [ \sum_{i=0}^{n} \delta_i = 1 ] 我們透過以下公式推匯出決策變數的值： [ x = \sum_{i=0}^{n} \delta_i B_i ] 對於求解器來說，$\delta$ 是真正的決策變數，而 $x$ 只是一個翻譯成原始問題語言的變數。這是關鍵。
import numpy as np
from scipy.optimize import linprog
# 定義問題引數
C = np.array([24, 28, 32, 34, 36, 40]) # 每段的單位成本
B = np.array([0, 148, 310, 501, 617, 762, 959]) # 斷點
# 定義目標函式係數
c = C * (B[1:] - B[:-1]) # 每段的總成本
# 定義約束條件
A_eq = np.array([[1] * len(C)]) # δ 的總和為 1
b_eq = np.array([1])
# 求解
res = linprog(c, A_eq=A_eq, b_eq=b_eq, method="highs")
print(res.x)
內容解密：
這段程式碼首先定義了分段線性函式的引數，包括每段的單位成本 C 和斷點 B。然後，它計算了每段的總成本，並將其用作線性規劃問題的目標函式係數。接著，定義了約束條件，即 $\delta$ 的總和為 1。最後，使用 linprog 函式求解這個線性規劃問題，並列印預出結果。
分段線性模型在最佳化問題中的應用
在最佳化問題中，我們經常遇到非線性函式的最佳化問題。雖然線性規劃求解器只能處理線性問題，但透過分段線性近似的方法，我們可以將某些非線性問題轉換為線性問題來求解。
凸分段線性目標函式的最佳化
目標函式的定義
考慮一個凸分段線性函式，其定義為多個線性片段的組合。我們的目標是最小化這個函式，同時滿足某些約束條件。
目標函式的數學表達
給定一個分段線性函式 $f(x)$，其在區間 $[B_i, B_{i+1}]$ 上的表示式為：
$f(x) = f(B_i) + \frac{f(B_{i+1}) - f(B_i)}{B_{i+1} - B_i} \times (x - B_i)$
我們的目標是最小化 $f(x)$，其中 $x \geq b$，且 $b$ 是給定的下界。
可執行的模型
我們將這個問題轉換為一個可執行的模型，如以下 Python 程式碼所示（Listing 3-1）：
def minimize_piecewise_linear_convex(Points, B):
s, n = newSolver('Piecewise'), len(Points)
x = s.NumVar(Points[0][0], Points[n-1][0], 'x')
l = [s.NumVar(0.0, 1, 'l[%i]' % (i,)) for i in range(n)]
s.Add(1 == sum(l[i] for i in range(n)))
s.Add(x == sum(l[i]*Points[i][0] for i in range(n)))
s.Add(x >= B)
Cost = s.Sum(l[i]*Points[i][1] for i in range(n))
s.Minimize(Cost)
s.Solve()
R = [l[i].SolutionValue() for i in range(n)]
return R
內容解密：
newSolver('Piecewise')：建立一個新的求解器例項，用於解決我們的最佳化問題。s.NumVar：定義數值變數，包括決策變數 $x$ 和權重變數 $l[i]$。s.Add(1 == sum(l[i] for i in range(n)))：約束所有權重變數的總和為 1，確保 $x$ 是分段線性函式中的一個有效點。s.Add(x == sum(l[i]*Points[i][0] for i in range(n)))：定義 $x$ 為各斷點的加權平均值。s.Add(x >= B)：施加下界約束，確保 $x$ 不小於給定的下界 $B$。Cost = s.Sum(l[i]*Points[i][1] for i in range(n))：計算目標函式值，即各斷點對應函式值的加權平均。s.Minimize(Cost)：設定求解器的目標為最小化目標函式。s.Solve()：執行最佳化求解。
應使用案例項
範例 1：簡單的分段線性最佳化
給定一個分段線性函式和下界 $x \geq 250$，我們可以呼叫 minimize_piecewise_linear_convex 函式來求解最佳化問題。結果如 Table 3-2 所示。
範例 2：邊界情況
當我們設定下界 $x \geq 310$ 時，求解結果如 Table 3-3 所示。只有一個 $\delta$ 非零，且其值為 1。
非線性函式的最佳化
對於非線性函式，我們可以透過分段線性近似的方法來進行最佳化。如 Listing 3-2 所示，我們定義了一個 minimize_non_linear 函式，該函式接受任意 Python 函式、定義域區間和精確度要求作為引數。
def minimize_non_linear(my_function, left, right, precision):
n = 5
while right - left > precision:
dta = (right - left) / (n - 1.0)
points = [(left + dta * i, my_function(left + dta * i)) for i in range(n)]
G = minimize_piecewise_linear_convex(points, left)
x = sum([G[i] * points[i][0] for i in range(n)])
left = points[max(0, [i - 1 for i in range(n) if G[i] > 0][0])][0]
right = points[min(n - 1, [i + 1 for i in range(n - 1, 0, -1) if G[i] > 0][0])][0]
return x
內容解密：
- 初始化
n = 5，表示初始分段數量。 - 在
while迴圈中，不斷細分割槽間直到滿足精確度要求。 - 在每次迭代中，計算分段點並呼叫
minimize_piecewise_linear_convex求解最佳化問題。 - 更新區間邊界，以縮小搜尋範圍。
本章節介紹瞭如何利用分段線性近似的方法來解決非線性最佳化問題，並提供了相應的 Python 程式碼實作。透過這種方法，我們可以將某些非線性問題轉換為線性問題，從而利用線性規劃求解器進行高效求解。