簡介這份資源是面向算法學(xué)習(xí)者與C開發(fā)者的單純形法完整實(shí)現(xiàn)程序用于求解線性規(guī)劃問題的最優(yōu)解。單純形法由丹齊格于1947年提出通過迭代在基本可行解之間轉(zhuǎn)換逐步逼近目標(biāo)函數(shù)最優(yōu)值是運(yùn)籌學(xué)與數(shù)值計(jì)算中的經(jīng)典算法。資源包共3個(gè)文件包含2個(gè)cpp源文件與1個(gè)h頭文件壓縮包約2KB代碼圍繞單純形表、基變量與非基變量替換、迭代終止條件等核心環(huán)節(jié)展開結(jié)構(gòu)緊湊便于直接閱讀與二次開發(fā)。目前已有742人學(xué)習(xí)下載適合作為課程實(shí)驗(yàn)、算法練習(xí)或工程入門的參考實(shí)現(xiàn)。讀者可從中理解線性規(guī)劃標(biāo)準(zhǔn)形式的轉(zhuǎn)換思路、單純形表的構(gòu)造方式以及C中矩陣與向量的組織技巧并借鑒異常處理與退出條件檢查等工程細(xì)節(jié)為后續(xù)性能優(yōu)化和測試調(diào)試打下基礎(chǔ)。1. 用 C 手搓單純形法從建模到求解的完整計(jì)算程序很多做運(yùn)籌學(xué)作業(yè)或線性規(guī)劃求解的開發(fā)者第一次接觸單純形法都是在紙上推演表格一旦變量數(shù)超過五個(gè)手工迭代基本就崩了。這份 C 實(shí)現(xiàn)的單純形算法計(jì)算程序解決的正是這個(gè)痛點(diǎn)把標(biāo)準(zhǔn)型線性規(guī)劃問題丟進(jìn)去程序自動完成基變換、入基出基判斷、迭代收斂最后輸出最優(yōu)解和目標(biāo)函數(shù)值。它適合兩類人一類是正在學(xué)運(yùn)籌學(xué)、需要驗(yàn)證手工計(jì)算結(jié)果的學(xué)生另一類是想把線性規(guī)劃求解嵌進(jìn)自己工程里的 C 開發(fā)者拿過去改改輸入輸出就能用。整個(gè)項(xiàng)目不依賴第三方數(shù)學(xué)庫純標(biāo)準(zhǔn) C 實(shí)現(xiàn)編譯即跑這也是它最實(shí)在的地方。2. 項(xiàng)目結(jié)構(gòu)與技術(shù)選型為什么用純 C 而不是調(diào)庫2.1 目錄組織與文件職責(zé)拿到這份源碼第一件事是看目錄。這類單純形法程序的結(jié)構(gòu)通常很緊湊常見做法是拆成三個(gè)部分主入口、算法核心、以及輸入輸出輔助。我一般會按下面的方式組織讀者對照自己手上的源碼可以快速定位simplex/ ├── main.cpp // 程序入口負(fù)責(zé)讀取問題、調(diào)用求解、打印結(jié)果 ├── simplex.h // 單純形算法類聲明定義數(shù)據(jù)結(jié)構(gòu)與接口 ├── simplex.cpp // 算法核心實(shí)現(xiàn)迭代邏輯都在這里 └── README.md // 編譯說明與輸入格式約定main.cpp只做三件事構(gòu)造問題、調(diào)用SimplexSolver::solve()、輸出結(jié)果。真正的計(jì)算全部壓在simplex.cpp里。這種拆分的好處是你想把求解器嵌到別的項(xiàng)目只需要拿走simplex.h和simplex.cppmain.cpp直接丟掉換自己的調(diào)用代碼。2.2 數(shù)據(jù)結(jié)構(gòu)設(shè)計(jì)表格法還是修正單純形單純形法有兩種主流實(shí)現(xiàn)路線表格法Tableau和修正單純形法Revised Simplex。這份程序用的是表格法原因是它直觀、易調(diào)試適合教學(xué)和中小規(guī)模問題。核心數(shù)據(jù)結(jié)構(gòu)是一張二維表常見做法是用std::vectorstd::vectordouble存整張單純形表行代表約束列代表變量加松弛變量最后一列存右端項(xiàng) b。// simplex.h class SimplexSolver { public: // 構(gòu)造時(shí)傳入目標(biāo)函數(shù)系數(shù)、約束系數(shù)矩陣、右端項(xiàng) SimplexSolver(const std::vectordouble c, const std::vectorstd::vectordouble A, const std::vectordouble b); bool solve(); // 執(zhí)行迭代返回是否找到最優(yōu)解 double getOptimalValue() const; // 取最優(yōu)目標(biāo)函數(shù)值 std::vectordouble getSolution() const; // 取最優(yōu)解向量 private: int m, n; // 約束數(shù)、原始變量數(shù) std::vectorstd::vectordouble tableau; // 單純形表 std::vectorint basis; // 基變量索引 void buildTableau(); // 初始化表格加入松弛變量 int findEnteringVariable(); // 入基變量檢驗(yàn)數(shù)為正且最大 int findLeavingVariable(int enterCol); // 出基變量最小比值檢驗(yàn) void pivot(int row, int col); // 高斯消元做基變換 };這里幾個(gè)參數(shù)要說明白m是約束個(gè)數(shù)n是決策變量個(gè)數(shù)。tableau的列數(shù)是n m 1多出來的m列是松弛變量最后 1 列是右端項(xiàng)。basis記錄每一行當(dāng)前哪個(gè)變量在基里初始時(shí)全是松弛變量。findEnteringVariable返回檢驗(yàn)數(shù)最大的列最大化問題findLeavingVariable用最小比值法則避免可行解跑出可行域。2.3 為什么不用第三方庫有開發(fā)者會問為什么不直接上某優(yōu)化庫原因很實(shí)際這份程序的目標(biāo)是教學(xué)和輕量嵌入。第三方庫動輒幾十兆編譯鏈復(fù)雜而單純形法本身邏輯清晰純 C 實(shí)現(xiàn)不到三百行。對于變量數(shù)在幾十以內(nèi)的線性規(guī)劃問題表格法的性能完全夠用。當(dāng)然如果你的問題規(guī)模上千那確實(shí)該換修正單純形法或?qū)I(yè)求解器這是后話。3. 核心算法拆解一次迭代到底發(fā)生了什么3.1 初始化把標(biāo)準(zhǔn)型變成單純形表單純形法的第一步是把問題轉(zhuǎn)成標(biāo)準(zhǔn)型。假設(shè)原始問題是最大化c^T x約束是Ax ≤ bx ≥ 0。程序通過加入松弛變量把不等式變成等式構(gòu)造初始基可行解。這一步在buildTableau()里完成。void SimplexSolver::buildTableau() { // 表格列數(shù) 決策變量 n 松弛變量 m 右端項(xiàng) 1 tableau.assign(m, std::vectordouble(n m 1, 0.0)); for (int i 0; i m; i) { for (int j 0; j n; j) { tableau[i][j] A[i][j]; // 填入約束系數(shù) } tableau[i][n i] 1.0; // 第 i 行松弛變量系數(shù)為 1 tableau[i][n m] b[i]; // 右端項(xiàng) } // 目標(biāo)函數(shù)行檢驗(yàn)數(shù)行放在最后一行初始為 -c std::vectordouble objRow(n m 1, 0.0); for (int j 0; j n; j) { objRow[j] -c[j]; // 最大化問題取負(fù)便于統(tǒng)一判斷 } tableau.push_back(objRow); // 初始基全是松弛變量 basis.resize(m); for (int i 0; i m; i) basis[i] n i; }邏輯說明tableau比約束多一行最后一行是檢驗(yàn)數(shù)行。把-c放進(jìn)去是因?yàn)樽畲蠡瘑栴}里檢驗(yàn)數(shù)為負(fù)表示還有改進(jìn)空間迭代到全非負(fù)就最優(yōu)。basis初始化為松弛變量索引對應(yīng)單位矩陣保證初始解可行。參數(shù)上A[i][j]是第 i 個(gè)約束里第 j 個(gè)變量的系數(shù)b[i]必須非負(fù)否則初始基不可行需要兩階段法這份程序默認(rèn)輸入滿足b ≥ 0。3.2 入基與出基最小比值法則的代碼落地迭代的核心是選誰進(jìn)基、誰出基。入基選檢驗(yàn)數(shù)最負(fù)的那一列對應(yīng)最大化改進(jìn)最快出基用最小比值法則保證右端項(xiàng)在消元后仍然非負(fù)。int SimplexSolver::findEnteringVariable() { int enterCol -1; double minVal -1e-9; // 容差避免浮點(diǎn)誤差誤判 int lastRow m; // 檢驗(yàn)數(shù)行索引 for (int j 0; j n m; j) { if (tableau[lastRow][j] minVal) { minVal tableau[lastRow][j]; enterCol j; } } return enterCol; // 返回 -1 表示已最優(yōu) } int SimplexSolver::findLeavingVariable(int enterCol) { int leaveRow -1; double minRatio 1e18; for (int i 0; i m; i) { double coef tableau[i][enterCol]; if (coef 1e-9) { // 只考慮正系數(shù)否則比值無意義 double ratio tableau[i][n m] / coef; if (ratio minRatio) { minRatio ratio; leaveRow i; } } } return leaveRow; // 返回 -1 表示問題無界 }參數(shù)說明1e-9是浮點(diǎn)比較容差單純形法迭代多了容易積累誤差沒有容差會出現(xiàn)死循環(huán)或誤判。findEnteringVariable返回 -1 時(shí)說明所有檢驗(yàn)數(shù)非負(fù)當(dāng)前解就是最優(yōu)解。findLeavingVariable返回 -1 時(shí)說明入基列所有系數(shù)非正問題無界目標(biāo)函數(shù)可以無限增大。這兩個(gè)返回值在solve()里分別對應(yīng)「最優(yōu)」和「無界」兩種終止?fàn)顟B(tài)。3.3 基變換pivot 消元的實(shí)現(xiàn)細(xì)節(jié)選定入基列和出基行后做高斯消元把入基列變成單位向量。這一步是數(shù)值計(jì)算最容易翻車的地方主元太小會導(dǎo)致誤差爆炸。void SimplexSolver::pivot(int row, int col) { double pivotVal tableau[row][col]; // 歸一化主元行 for (int j 0; j n m; j) { tableau[row][j] / pivotVal; } // 其他行消元包括檢驗(yàn)數(shù)行 for (int i 0; i m; i) { if (i row) continue; double factor tableau[i][col]; if (std::fabs(factor) 1e-12) continue; for (int j 0; j n m; j) { tableau[i][j] - factor * tableau[row][j]; } } basis[row] col; // 更新基變量索引 }邏輯說明先把主元行整體除以主元讓主元變成 1然后對其他每一行減去factor倍的主元行把該列其他元素清零。注意循環(huán)上界是i m因?yàn)闄z驗(yàn)數(shù)行也要消元。basis[row] col更新基索引保證后續(xù)能正確還原解向量。主元判斷用1e-12是防止除以極小數(shù)實(shí)際工程里如果主元絕對值小于閾值通常要觸發(fā)重新選主元或報(bào)數(shù)值不穩(wěn)定。3.4 主循環(huán)與終止條件把上面幾個(gè)函數(shù)串起來solve()就是一個(gè) while 循環(huán)找入基、找出基、pivot直到最優(yōu)或無界。bool SimplexSolver::solve() { buildTableau(); const int MAX_ITER 10000; // 防止死循環(huán)的后悔藥 int iter 0; while (iter MAX_ITER) { int enterCol findEnteringVariable(); if (enterCol -1) return true; // 最優(yōu) int leaveRow findLeavingVariable(enterCol); if (leaveRow -1) return false; // 無界 pivot(leaveRow, enterCol); } return false; // 迭代超限視為失敗 }MAX_ITER是血淚經(jīng)驗(yàn)浮點(diǎn)誤差可能讓檢驗(yàn)數(shù)在 0 附近反復(fù)橫跳沒有迭代上限程序可能永遠(yuǎn)不退出。設(shè)成 10000 對中小規(guī)模問題足夠真超了說明輸入有問題或需要換數(shù)值更穩(wěn)的算法。4. 編譯運(yùn)行與輸入格式從源碼到可執(zhí)行文件4.1 編譯命令與依賴這份程序沒有外部依賴只要有 C11 以上的編譯器就能編。Linux 或 macOS 下g -stdc11 -O2 -o simplex main.cpp simplex.cpp ./simplexWindows 下用 MinGW 或 Visual Studio 命令行同理把源文件一起編譯即可。-O2開優(yōu)化對單純形法這種循環(huán)密集的代碼有實(shí)際提升。如果編譯報(bào)fabs未定義檢查有沒有#include cmath這是新手最常踩的坑。4.2 輸入格式約定程序讀取問題的方式通常有兩種硬編碼在main.cpp里或者從標(biāo)準(zhǔn)輸入讀。常見做法是前者方便調(diào)試。一個(gè)典型輸入長這樣// main.cpp 中的問題定義 std::vectordouble c {3, 5}; // 目標(biāo)函數(shù)系數(shù) std::vectorstd::vectordouble A { {1, 0}, // 約束1: x1 4 {0, 2}, // 約束2: 2*x2 12 {3, 2} // 約束3: 3*x1 2*x2 18 }; std::vectordouble b {4, 12, 18}; // 右端項(xiàng) SimplexSolver solver(c, A, b); if (solver.solve()) { std::cout 最優(yōu)值: solver.getOptimalValue() std::endl; auto sol solver.getSolution(); for (size_t i 0; i sol.size(); i) std::cout x i 1 sol[i] std::endl; }參數(shù)說明c是目標(biāo)函數(shù)系數(shù)長度等于決策變量數(shù)A是約束系數(shù)矩陣行數(shù)等于約束數(shù)列數(shù)等于變量數(shù)b是右端項(xiàng)必須非負(fù)。這個(gè)例子是經(jīng)典的線性規(guī)劃入門題最優(yōu)解應(yīng)該是x1 2, x2 6目標(biāo)值 36。跑出來對不上先檢查b有沒有負(fù)數(shù)再檢查約束方向是不是都統(tǒng)一成了≤。4.3 結(jié)果還原從表格里取出解向量迭代結(jié)束后解藏在tableau的右端項(xiàng)列里但只有基變量有值非基變量全是 0。getSolution()要做的是按basis索引把值填回原變量位置。std::vectordouble SimplexSolver::getSolution() const { std::vectordouble x(n, 0.0); // 決策變量初始為 0 for (int i 0; i m; i) { int var basis[i]; if (var n) { // 只關(guān)心原始決策變量 x[var] tableau[i][n m]; } } return x; }邏輯說明遍歷每一行看基變量索引var如果它落在原始變量范圍[0, n)內(nèi)就把該行右端項(xiàng)賦給對應(yīng)變量。松弛變量的值不需要輸出。getOptimalValue()更簡單直接把檢驗(yàn)數(shù)行最后一列取出來注意符號——因?yàn)槌跏蓟瘯r(shí)目標(biāo)函數(shù)取了負(fù)這里要取反。5. 避坑與常見問題迭代不收斂、結(jié)果對不上怎么辦5.1 現(xiàn)象程序跑出無界但問題明明有最優(yōu)解原因通常是約束方向沒統(tǒng)一。單純形法標(biāo)準(zhǔn)型要求所有約束是≤且b ≥ 0。如果原始問題里有≥或等式約束直接塞進(jìn)去會導(dǎo)致初始基不可行程序誤判無界。解決辦法是先把問題轉(zhuǎn)成標(biāo)準(zhǔn)型≥約束兩邊乘 -1 變≤但這樣b可能變負(fù)需要引入人工變量走兩階段法。這份程序默認(rèn)不做兩階段所以輸入前必須自己保證b非負(fù)。5.2 現(xiàn)象迭代次數(shù)很多但檢驗(yàn)數(shù)一直在 0 附近震蕩這是浮點(diǎn)誤差導(dǎo)致的退化循環(huán)。原因在于主元太小消元后誤差被放大。解決方式有兩個(gè)一是加容差把1e-9調(diào)大一點(diǎn)比如1e-7二是加 Bland 規(guī)則入基時(shí)選下標(biāo)最小的負(fù)檢驗(yàn)數(shù)列而不是選最負(fù)的能理論上避免循環(huán)。我一般會在findEnteringVariable里加個(gè)開關(guān)退化嚴(yán)重時(shí)切 Bland 規(guī)則。5.3 現(xiàn)象最優(yōu)值和手工算的對不上差一點(diǎn)點(diǎn)先別懷疑算法九成是浮點(diǎn)精度問題。檢查輸出時(shí)有沒有做四舍五入0.9999999和1.0在程序里是兩個(gè)數(shù)。常見做法是在輸出前統(tǒng)一處理小于1e-6的當(dāng)成 0接近整數(shù)的保留整數(shù)。另外檢查目標(biāo)函數(shù)值有沒有取反初始化時(shí)-c進(jìn)去取結(jié)果時(shí)忘了取反符號就反了。5.4 現(xiàn)象變量數(shù)一多程序變慢甚至卡死表格法的復(fù)雜度是 O(m*n) 每輪迭代迭代輪數(shù)又和變量數(shù)相關(guān)規(guī)模上去確實(shí)會慢。如果變量數(shù)超過幾百建議換修正單純形法只維護(hù)基矩陣的逆內(nèi)存和計(jì)算量都小一個(gè)量級。這份程序定位是教學(xué)和中小規(guī)模硬扛大規(guī)模問題不是它的設(shè)計(jì)目標(biāo)。5.5 現(xiàn)象編譯通過但運(yùn)行崩潰提示 vector 越界多半是A的行列數(shù)和c、b對不上。A的行數(shù)必須等于b的長度列數(shù)必須等于c的長度。構(gòu)造SimplexSolver時(shí)最好加個(gè)斷言檢查assert(A.size() b.size())和assert(A[0].size() c.size())把錯(cuò)誤擋在構(gòu)造階段比運(yùn)行到一半崩潰好排查得多。6. 二次開發(fā)改造從教學(xué)程序到工程可用的求解器這份程序拿來直接跑沒問題但要嵌進(jìn)真實(shí)工程還得做幾處改造。第一處是輸入輸出解耦。現(xiàn)在問題硬編碼在main.cpp改成從文件或標(biāo)準(zhǔn)輸入讀格式用簡單的文本協(xié)議比如第一行放m n接著c再接著A和b。這樣你的上層業(yè)務(wù)代碼可以用任意語言生成問題文件C 程序只負(fù)責(zé)求解。第二處是錯(cuò)誤處理?，F(xiàn)在的solve()只返回 bool信息量太少。我一般會改成返回枚舉OPTIMAL、UNBOUNDED、INFEASIBLE、MAX_ITER_REACHED再配一個(gè)getErrorMessage()。這樣調(diào)用方能區(qū)分「問題本身無解」和「程序沒算出來」排查方向完全不同。第三處是數(shù)值穩(wěn)定性。把主元閾值、檢驗(yàn)數(shù)容差做成可配置參數(shù)退化時(shí)自動切 Bland 規(guī)則。下面是一個(gè)改造后的接口示意enum class SolveStatus { OPTIMAL, UNBOUNDED, INFEASIBLE, MAX_ITER }; struct SolverConfig { double pivotTol 1e-9; // 主元最小閾值 double objTol 1e-7; // 檢驗(yàn)數(shù)容差 int maxIter 10000; // 迭代上限 bool useBland false; // 是否強(qiáng)制 Bland 規(guī)則 }; class SimplexSolver { public: SimplexSolver(const std::vectordouble c, const std::vectorstd::vectordouble A, const std::vectordouble b, const SolverConfig cfg {}); SolveStatus solve(); // ... 其余接口不變 };改造完怎么驗(yàn)證我的習(xí)慣是準(zhǔn)備一組已知答案的測試用例覆蓋最優(yōu)、無界、退化三種情況每次改完跑一遍。特別是退化用例構(gòu)造一個(gè)檢驗(yàn)數(shù)反復(fù)為 0 的問題看程序能不能在有限步內(nèi)退出。這一步別省單純形法的數(shù)值坑基本都藏在退化場景里。從那以后我每次改求解器核心都強(qiáng)制走一遍「已知答案用例 退化用例 大規(guī)模隨機(jī)用例」三件套少一件都不敢提交。希望這份拆解幫到你拿到源碼先跑通那個(gè)x12, x26的例子再動手改輸入輸出路子就順了。本文還有配套的精品資源點(diǎn)擊獲取