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