簡介這份資源面向學習最優(yōu)化方法、線性規(guī)劃課程的高校學生與工程技術人員聚焦單純形法求解中的大M法與兩階段法編程實現(xiàn)幫助讀者理解算法原理并完成課程作業(yè)或實際建模求解。壓縮包共3個文件均為MATLAB的m腳本整體約4KB體量輕便包含主程序入口、兩階段法函數(shù)與大M法單純形函數(shù)代碼注釋詳盡、邏輯通俗易懂。主程序負責輸入約束方程與目標函數(shù)自動調用兩階段法或大M法完成規(guī)劃求解其中兩階段法的單純形迭代通過調用大M法函數(shù)實現(xiàn)便于對照兩種方法的異同。目前已有2547人學習下載讀者可借此掌握單純形表構建、人工變量處理、迭代終止判斷等關鍵環(huán)節(jié)并在此基礎上修改參數(shù)求解自己的線性規(guī)劃問題是入門與復習單純形法的實用參考。1. 單純形法為什么需要人造變量從標準型到初始基可行解的鴻溝線性規(guī)劃的標準型看起來人畜無害目標函數(shù)求極大約束全是等式變量非負。教科書上寫得很輕松——找到初始基可行解然后迭代換基直到檢驗數(shù)全部非正。但真正動手寫代碼的人會撞上一堵墻約束里那些大于等于和等式根本湊不出一組現(xiàn)成的單位基向量。沒有單位基單純形表連第一行都填不出來。這就是大M法和兩階段法存在的理由。它們不是單純形法的高級變體而是專門用來跨越初始基可行解不存在這道鴻溝的工程手段。很多人學最優(yōu)化方法時把這兩個方法當成兩套獨立算法背下來其實它們解決的是同一個問題只是代價函數(shù)的設計哲學不同。搞懂這一點才算真正理解了單純形法的完整鏈路。這份資源把大M法和兩階段法的程序實現(xiàn)放在一起恰好覆蓋了從數(shù)學模型到可運行代碼的關鍵環(huán)節(jié)。下面從為什么需要人造變量開始一層層拆到代碼里的具體處理。2. 人造變量的本質為什么不能直接對不等式做消元2.1 松弛變量與剩余變量的分工先把約束類型理清楚。標準型要求所有約束是等式所以不等式要轉化小于等于約束左邊加一個松弛變量系數(shù) 1直接充當基變量大于等于約束左邊減一個剩余變量系數(shù) -1但 -1 不能當基變量基變量要求對應列是單位向量等式約束沒有任何多余變量可以借用問題就出在后兩種。剩余變量系數(shù)是 -1等式約束連個多余的變量都沒有。單純形法的第一步是構造初始單純形表要求基矩陣是單位陣。如果約束里湊不出單位陣迭代就無法啟動。常見的實現(xiàn)思路是對每個缺基變量的約束人為添加一個變量讓它在該約束中系數(shù)為 1在其他約束中系數(shù)為 0。這樣所有添加的變量就構成了一組單位基。這些變量就是人造變量。但人造變量是憑空造出來的它們并不屬于原問題。如果最終解里人造變量不為零說明原問題無可行解。所以人造變量必須被趕出基或者被迫取零值。大M法和兩階段法的區(qū)別就在于用什么手段逼人造變量歸零。2.2 大M法的懲罰邏輯大M法的思路很直接在目標函數(shù)里給人造變量加一個巨大的懲罰系數(shù) M。如果是求極大值人造變量的系數(shù)就是 -M求極小值則相反。M 是一個理論上無窮大的正數(shù)。這樣做的效果是只要人造變量還留在基里且取正值目標函數(shù)就會被扣掉一個巨大的量單純形法在迭代中會優(yōu)先把它換出去。代碼實現(xiàn)時M 不能真的取無窮大通常取一個遠大于所有原始系數(shù)的數(shù)比如 10^6 或 10^8。# 大M法構造初始單純形表的核心邏輯 def build_big_m_table(c, A, b, constraint_types, M1e6): c: 原始目標函數(shù)系數(shù)求極大 A: 約束系數(shù)矩陣 b: 右端項 constraint_types: 每個約束的類型 、、 M: 懲罰系數(shù) m, n A.shape total_vars n artificial_indices [] for i, ctype in enumerate(constraint_types): if ctype : # 加松弛變量系數(shù) 1可作基變量 slack np.zeros((m, 1)) slack[i] 1 A np.hstack([A, slack]) c np.append(c, 0) total_vars 1 elif ctype : # 加剩余變量系數(shù) -1不能作基變量 surplus np.zeros((m, 1)) surplus[i] -1 A np.hstack([A, surplus]) c np.append(c, 0) total_vars 1 # 再加人造變量系數(shù) 1 art np.zeros((m, 1)) art[i] 1 A np.hstack([A, art]) c np.append(c, -M) # 求極大懲罰為 -M artificial_indices.append(total_vars) total_vars 1 else: # art np.zeros((m, 1)) art[i] 1 A np.hstack([A, art]) c np.append(c, -M) artificial_indices.append(total_vars) total_vars 1 return A, b, c, artificial_indices這段代碼的關鍵在于松弛變量系數(shù) 1 可以直接當基變量剩余變量系數(shù) -1 不行所以大于等于約束需要同時加剩余變量和人造變量。等式約束只需要加人造變量。人造變量在目標函數(shù)里的系數(shù)是 -M這個 -M 會隨著迭代進入檢驗數(shù)計算驅動人造變量出基。M 的取值是個工程問題。取太小懲罰力度不夠人造變量可能賴在基里不走取太大浮點運算時容易和其他系數(shù)產生數(shù)量級差異導致數(shù)值不穩(wěn)定。常見的做法是取原始系數(shù)最大絕對值的 1000 到 10000 倍。2.3 兩階段法的分步策略兩階段法不依賴 M 這個魔法數(shù)字而是把問題拆成兩個階段第一階段構造一個輔助問題目標函數(shù)是極小化所有人造變量之和。如果輔助問題的最優(yōu)值為零說明人造變量可以全部取零原問題有可行解如果最優(yōu)值大于零原問題無可行解。第二階段在第一階段得到的基可行解基礎上去掉人造變量列恢復原始目標函數(shù)繼續(xù)單純形迭代。# 兩階段法第一階段極小化人造變量之和 def phase_one(A, b, artificial_indices): 構造輔助目標函數(shù)min sum(artificial_vars) 等價于 max -sum(artificial_vars) n A.shape[1] c_phase1 np.zeros(n) for idx in artificial_indices: c_phase1[idx] -1 # 求極大化 -sum等價于極小化 sum # 用單純形法求解輔助問題 result simplex(A, b, c_phase1, basisartificial_indices) if abs(result[objective]) 1e-8: return None # 原問題無可行解 # 檢查人造變量是否還在基里若在且值為零嘗試換出 return result兩階段法的優(yōu)勢在于不需要選 M。M 選得不好會導致數(shù)值問題而兩階段法把可行性判斷和最優(yōu)性優(yōu)化徹底分開邏輯更清晰。代價是要跑兩遍單純形計算量略大。從設計取舍看大M法代碼更緊湊適合手寫實現(xiàn)和教學演示兩階段法數(shù)值更穩(wěn)健適合放進通用求解器。這份資源把兩者都實現(xiàn)了對比著看能理解為什么工業(yè)級求解器更傾向兩階段法。3. 單純形表的數(shù)據(jù)結構基變量、檢驗數(shù)與非基變量怎么組織3.1 單純形表的矩陣表示單純形法的每次迭代本質上是一次矩陣變換。把約束寫成 $Ax b$把 $A$ 按基變量和非基變量分塊$A [B | N]$其中 $B$ 是基矩陣$N$ 是非基矩陣。對應的變量也分成 $x_B$ 和 $x_N$。約束變?yōu)?$B x_B N x_N b$所以 $x_B B^{-1}b - B^{-1}N x_N$。代入目標函數(shù) $z c_B^T x_B c_N^T x_N$得到$$z c_B^T B^{-1}b (c_N^T - c_B^T B^{-1}N) x_N$$令非基變量 $x_N 0$得到基可行解 $x_B B^{-1}b$目標值 $z c_B^T B^{-1}b$。檢驗數(shù)就是 $\sigma_N c_N^T - c_B^T B^{-1}N$。代碼里不需要每次都求 $B^{-1}$而是用高斯-若爾當消元維護整張單純形表。表的結構是基變量$x_1$$x_2$...$x_n$右端項$x_{B1}$$a_{11}$$a_{12}$...$a_{1n}$$b_1$$x_{B2}$$a_{21}$$a_{22}$...$a_{2n}$$b_2$..................檢驗數(shù)$\sigma_1$$\sigma_2$...$\sigma_n$$z$每次迭代選一個檢驗數(shù)為正的非基變量入基求極大時再按最小比值法則選一個基變量出基然后做一次初等行變換把入基變量對應的列變成單位向量。3.2 入基與出基的判定細節(jié)入基變量的選擇直接影響收斂速度。最簡單的規(guī)則是選檢驗數(shù)最大的那個叫 Dantzig 規(guī)則。但 Dantzig 規(guī)則容易在退化問題上循環(huán)。常見的改進是 Bland 規(guī)則選下標最小的正檢驗數(shù)變量入基同時出基也選下標最小的保證不循環(huán)但收斂慢。出基判定用最小比值法則對每個系數(shù) $a_{ik} 0$ 的行計算 $b_i / a_{ik}$選最小的那個對應的基變量出基。如果所有 $a_{ik} \leq 0$說明問題無界。def simplex_iteration(table, basis, c, maximizeTrue): table: 單純形表最后一行為檢驗數(shù) basis: 基變量下標列表 c: 目標函數(shù)系數(shù) m, n table.shape n - 1 # 最后一列是右端項 # 1. 選入基變量找最大正檢驗數(shù)求極大 reduced_costs table[-1, :n] if maximize: entering np.argmax(reduced_costs) if reduced_costs[entering] 1e-8: return None # 已最優(yōu) else: entering np.argmin(reduced_costs) if reduced_costs[entering] -1e-8: return None # 2. 最小比值法則選出基變量 ratios [] for i in range(m - 1): if table[i, entering] 1e-8: ratios.append(table[i, -1] / table[i, entering]) else: ratios.append(np.inf) if min(ratios) np.inf: return unbounded # 無界解 leaving np.argmin(ratios) # 3. 初等行變換把 entering 列變成單位向量 pivot table[leaving, entering] table[leaving, :] / pivot for i in range(m): if i ! leaving: table[i, :] - table[i, entering] * table[leaving, :] basis[leaving] entering return table, basis這段代碼里有個容易翻車的點浮點數(shù)比較不能用 0必須用容差。檢驗數(shù)判斷最優(yōu)時用 1e-8比值判斷時用 1e-8。容差取太小會誤判取太大會丟精度。常見的做法是取 1e-9 到 1e-6 之間具體看問題規(guī)模。3.3 退化與循環(huán)的處理退化是指基可行解里出現(xiàn)零分量導致最小比值出現(xiàn)多個相同的最小值。退化本身不致命但可能引發(fā)循環(huán)——基變量組合重復出現(xiàn)迭代永遠不終止。處理循環(huán)的常見手段有三種Bland 規(guī)則、字典序法、擾動法。Bland 規(guī)則實現(xiàn)最簡單但收斂慢字典序法在出基判定時引入額外的比較維度實現(xiàn)復雜但效率好擾動法給右端項加微小擾動破壞退化條件實現(xiàn)也不復雜。代碼里如果只做教學演示用 Bland 規(guī)則就夠了。如果要做通用求解器建議用字典序法。這份資源里的實現(xiàn)如果沒處理循環(huán)遇到退化問題可能會卡住這是使用時需要留意的。4. 從數(shù)學公式到可運行代碼大M法與兩階段法的實現(xiàn)差異4.1 大M法的完整迭代流程大M法的代碼結構可以拆成四步構造初始表、迭代換基、判斷終止、提取解。構造初始表時需要根據(jù)約束類型決定加松弛變量、剩余變量還是人造變量。這一步的代碼在 2.2 節(jié)已經展示過。迭代換基就是反復調用 3.2 節(jié)的simplex_iteration。終止條件有兩個檢驗數(shù)全部非正最優(yōu)或者入基列所有系數(shù)非正無界。提取解的時候要注意如果最優(yōu)解里人造變量不為零說明原問題無可行解。這個判斷不能省否則會把一個無可行解的問題當成有解輸出。def solve_big_m(c, A, b, constraint_types, M1e6): A, b, c_ext, art_indices build_big_m_table(c, A, b, constraint_types, M) m, n A.shape # 初始基變量松弛變量和人造變量 basis [] for i, ctype in enumerate(constraint_types): if ctype : basis.append(original_n i) # 松弛變量下標 else: basis.append(art_indices.pop(0)) # 構造單純形表 table np.hstack([A, b.reshape(-1, 1)]) # 計算初始檢驗數(shù)sigma c - c_B * B^{-1} * A c_B c_ext[basis] table[-1, :] c_ext - c_B table[:-1, :] # 迭代 while True: result simplex_iteration(table, basis, c_ext) if result is None: break if result unbounded: return {status: unbounded} # 檢查人造變量是否為零 solution np.zeros(n) for i, b_idx in enumerate(basis): solution[b_idx] table[i, -1] for idx in art_indices: if abs(solution[idx]) 1e-8: return {status: infeasible} return {status: optimal, solution: solution[:original_n], objective: table[-1, -1]}這段代碼里有個細節(jié)初始檢驗數(shù)的計算不能直接抄目標函數(shù)系數(shù)因為基變量的檢驗數(shù)必須是零。需要用 $c_B$ 做一次修正。很多人手寫單純形法時在這里翻車迭代半天發(fā)現(xiàn)檢驗數(shù)不對其實是初始表沒構造對。4.2 兩階段法的階段銜接兩階段法的代碼比大M法多一個階段切換的邏輯。第一階段結束后如果輔助問題最優(yōu)值為零需要做兩件事去掉人造變量列恢復原始目標函數(shù)重新計算檢驗數(shù)。去掉人造變量列時要注意如果某個人造變量還在基里但值為零不能直接刪列否則基矩陣會缺一個單位向量。常見的處理是找一個非基的非人造變量通過行變換把它換進基再刪掉人造變量列。如果找不到這樣的變量說明該約束是冗余的可以直接刪掉整行。def solve_two_phase(c, A, b, constraint_types): # 第一階段 A1, b1, c1, art_indices build_phase_one_table(A, b, constraint_types) result1 simplex_solve(A1, b1, c1, basisart_indices) if abs(result1[objective]) 1e-8: return {status: infeasible} # 階段切換去掉人造變量列恢復原始目標 keep_cols [i for i in range(A1.shape[1]) if i not in art_indices] A2 A1[:, keep_cols] c2 np.zeros(A2.shape[1]) c2[:len(c)] c # 原始變量系數(shù) # 用第一階段的基作為第二階段的初始基 basis2 [keep_cols.index(b) for b in result1[basis] if b in keep_cols] # 第二階段迭代 result2 simplex_solve(A2, b1, c2, basisbasis2) return result2階段切換時最容易出問題的地方是基變量的映射。第一階段的基變量下標是原始變量、松弛變量、剩余變量、人造變量的混合下標去掉人造變量列后下標需要重新映射。如果映射錯了第二階段的初始基就不是可行基迭代會從錯誤的地方開始。4.3 兩種方法的數(shù)值表現(xiàn)對比對比維度大M法兩階段法代碼復雜度低一次構造高需要階段切換數(shù)值穩(wěn)定性依賴 M 的取值較穩(wěn)定迭代次數(shù)通常較少通常較多無可行解判斷檢查人造變量是否為零檢查第一階段最優(yōu)值適用場景教學演示、小規(guī)模問題通用求解器、大規(guī)模問題從工程角度看大M法的 M 是個魔法數(shù)字取大了數(shù)值不穩(wěn)定取小了懲罰不夠。兩階段法沒有這個問題但代碼量更大。如果只是手寫一個求解器驗證算法大M法更快如果要處理實際數(shù)據(jù)兩階段法更可靠。5. 避坑與排查單純形法實現(xiàn)中最容易翻車的五個地方5.1 檢驗數(shù)符號搞反導致迭代方向錯誤現(xiàn)象程序迭代幾次后目標函數(shù)值越來越小或者直接報無界但問題明明有有限最優(yōu)解。原因求極大和求極小的檢驗數(shù)判定方向相反。求極大時選正檢驗數(shù)入基求極小時選負檢驗數(shù)入基。如果目標函數(shù)是求極小但代碼里用了求極大的判定邏輯迭代方向就反了。解決在代碼里顯式區(qū)分 maximize 和 minimize或者在構造目標函數(shù)時統(tǒng)一轉成求極大。轉法很簡單min z 等價于 max -z把目標函數(shù)系數(shù)全部取反即可。5.2 最小比值法則漏判無界現(xiàn)象程序在某次迭代中找不到出基變量但沒有報無界而是繼續(xù)迭代或者崩潰。原因最小比值法則要求只對入基列系數(shù)為正的行計算比值。如果所有系數(shù)都非正說明入基變量可以無限增大而不違反約束問題無界。代碼里如果沒判斷這個條件argmin會返回一個無效下標。解決在計算比值前先檢查入基列是否有正系數(shù)。如果沒有直接返回無界。代碼里用if min(ratios) np.inf或者if all(a 1e-8 for a in column)來判斷。5.3 浮點誤差累積導致檢驗數(shù)誤判現(xiàn)象程序在接近最優(yōu)時反復迭代檢驗數(shù)在零附近震蕩永遠不滿足終止條件。原因浮點運算有誤差理論上應該為零的檢驗數(shù)可能變成 1e-15 或者 -1e-15。如果終止條件用sigma 0嚴格判斷就會因為一個極小的正數(shù)繼續(xù)迭代。解決所有浮點比較都加容差。檢驗數(shù)判斷用sigma 1e-8比值判斷用a 1e-8。容差不能太大否則會丟失精度也不能太小否則起不到作用。1e-8 是個常用的經驗值。5.4 人造變量未完全出基導致誤判無可行解現(xiàn)象大M法求解一個明明有可行解的問題程序卻報無可行解。原因人造變量在最優(yōu)解里可能還留在基里但值為零。如果代碼只檢查人造變量是否在基里而不檢查其值是否為零就會誤判。或者 M 取得太小人造變量沒有被完全趕出基值不為零。解決判斷無可行解的條件是人造變量的值大于容差而不是人造變量是否在基里。同時檢查 M 的取值是否足夠大通常取原始系數(shù)最大絕對值的 1000 倍以上。5.5 階段切換時基變量映射錯誤現(xiàn)象兩階段法第一階段正常結束第二階段一開始就報無界或者迭代到錯誤的最優(yōu)解。原因第一階段結束后基變量列表里可能包含人造變量。去掉人造變量列時如果直接刪列而不調整基變量下標基矩陣就不再是單位陣第二階段的初始基不可行。解決階段切換時先檢查基里是否有人造變量。如果有且值為零找一個非基的非人造變量換進去如果找不到說明該約束冗余刪掉整行。然后再刪人造變量列重新映射基變量下標。6. 設計取舍為什么工業(yè)求解器很少用大M法大M法和兩階段法在教科書里是并列的但翻開源碼級的求解器實現(xiàn)會發(fā)現(xiàn)兩階段法或者它的變體才是主流。這不是偶然而是數(shù)值計算領域幾十年踩坑的結果。M 的選取本質上是一個沒有正確答案的問題。取 10^6遇到系數(shù)本身就是 10^5 量級的問題懲罰力度不夠取 10^12浮點數(shù)的有效位數(shù)只有 15 到 16 位M 和其他系數(shù)的加減運算會直接丟精度。更麻煩的是M 的影響會滲透到整個單純形表的每一次行變換中誤差不斷累積最終可能導致基變量判定錯誤。我一般會這么理解大M法是用一個人造的大數(shù)去模擬無窮大但計算機里的浮點數(shù)沒有無窮大只有有限精度這個模擬從一開始就是有損的。兩階段法把可行性判斷和最優(yōu)性優(yōu)化拆開第一階段只關心人造變量能不能歸零不涉及原始目標函數(shù)。這樣做的代價是多跑一遍單純形但換來的是數(shù)值穩(wěn)定性。對于大規(guī)模問題多跑一遍的代價遠小于數(shù)值不穩(wěn)定導致的錯誤結果。工業(yè)求解器還會在第一階段用專門的可行性松弛技術進一步降低數(shù)值風險。從代碼實現(xiàn)的角度看兩階段法的階段切換邏輯確實更復雜需要處理基變量映射、冗余約束刪除、人造變量換出等邊界情況。但這些復雜度是一次性的——寫對了就對了。大M法的復雜度是參數(shù)性的——M 取多少永遠需要根據(jù)問題調整換個數(shù)據(jù)集可能就要重新調。對于需要反復求解不同問題的場景兩階段法的總成本更低。這份資源把兩種方法都實現(xiàn)了價值在于能對比著看。先跑大M法理解懲罰機制再跑兩階段法理解可行性判斷的獨立性最后對比兩者在同一組測試數(shù)據(jù)上的迭代次數(shù)和解的質量。如果發(fā)現(xiàn)大M法在某些數(shù)據(jù)上結果異常大概率是 M 的取值問題換成兩階段法就能驗證。從那以后我每次實現(xiàn)單純形法都會先用兩階段法跑一遍確認可行性和最優(yōu)解再用大M法做交叉驗證。如果兩者結果不一致優(yōu)先信兩階段法然后回頭查大M法的 M 是不是取小了。這個習慣幫我省了很多調試時間。希望幫到你。本文還有配套的精品資源點擊獲取