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