化CCG與Benders代碼實戰(zhàn)：原理、實現(xiàn)與避坑)
簡介這份代碼包圍繞兩階段魯棒優(yōu)化中的CCG列與約束生成及Benders分解算法面向使用MATLAB與YALMIP的研究者和工程師既適合快速復現(xiàn)原論文基線也便于按實際決策場景擴展改編。壓縮包共5個文件、體積1.44MB含3篇PDF文獻與2個MATLAB腳本文獻涵蓋CCG與Benders對比、兩階段魯棒優(yōu)化求解方法以及Benders分解算法綜述腳本則提供可運行的案例代碼實現(xiàn)主問題構建、切割生成與迭代求解的完整鏈路。資源已有2287人學習下載在魯棒優(yōu)化與分解算法主題上具備較好的實踐參考熱度。通過學習包內(nèi)代碼與配套文獻讀者可掌握在YALMIP中搭建兩階段主從問題框架、利用CCG/Benders生成最優(yōu)割并識別最壞場景的關鍵思路繼而將原論文方法遷移至新能源出力、供應鏈需求等含不確定參數(shù)的優(yōu)化問題中為后續(xù)算法對比和工程落地提供可復現(xiàn)基礎。1. 兩階段魯棒優(yōu)化的CCG與Benders代碼包拿到就能跑但先搞懂割長什么樣做過兩階段魯棒優(yōu)化的人應該都有體會模型寫出來不難難的是把CCG或Benders的迭代循環(huán)跑通并讓結果能復現(xiàn)原論文的數(shù)值。這套帶MATLABYALMIP源碼和對比文獻的壓縮包解決的就是“從論文公式到可運行代碼”這段最磨人的路。它不是一個教學Demo而是把兩階段魯棒優(yōu)化里兩個主流求解框架——列與約束生成CCG和Benders分解——做成了可以直接改參數(shù)、換模型、對結果的可擴展腳本。適合正在寫相關論文的研究生也適合要把確定性優(yōu)化改成魯棒模型的工程師。一個反直覺的結論先放在這里CCG和Benders看起來是兩條完全不同的技術路線但實際上它們共用同一個“主問題-子問題”迭代骨架差別只在“切割”的構造方式理解這一點改造代碼能少走一半彎路。2. CCG列與約束生成主問題-子問題迭代與YALMIP代碼骨架2.1 CCG為什么快列生成與約束生成的雙重含義兩階段魯棒優(yōu)化的一般形式是 min c?x max_{u∈U} min_{y∈F(x,u)} d?y第一階段變量x要現(xiàn)在拍板第二階段變量y要等不確定參數(shù)u實現(xiàn)后再調(diào)整。CCG的思路非常直接不要試圖遍歷連續(xù)的不確定集U而是維護一個“已發(fā)現(xiàn)的極端場景”集合每輪找到當前最壞的那個u*把它對應的第二階段變量和約束整組塞進主問題。這里有個關鍵點CCG同時在做列生成和約束生成。主問題每加入一個新場景u*就會新增一組第二階段變量y^k這是“列”同時新增這組變量對應的約束需求滿足、容量限制、與θ的連線這是“約束”。因為不確定集的極值一定發(fā)生在頂點上線性目標函數(shù)在凸多面體上取極值必在頂點所以每輪找出的u*都是一個個頂點主問題不斷增強對“最壞情況”的描述。相比傳統(tǒng)Benders只往主問題丟一條對偶割CCG犧牲了主問題的規(guī)模換來了更緊的下界通常迭代次數(shù)明顯更少。2.2 CCG主問題與最壞場景搜索可直接修改的YALMIP實現(xiàn)下面用一個小規(guī)模的設施選址加不確定需求運輸問題演示CCG主循環(huán)。問題設定是I個候選設施建設成本c?容量cap?J個客戶需求u?在基準值上下波動第二階段決策是從設施i到客戶j的運輸量y??單位成本d??。第二階段目標是最小化運費整體目標是在最壞需求下最小化“建設成本運費”。% demo_ccg_2stage.m % 兩階段魯棒優(yōu)化CCG 主循環(huán)YALMIP MILP 求解器 % 模型設施選址 不確定需求運輸 % 第一階段0/1 決策是否建設施 % 第二階段滿足需求的最小運費 nI 2; % 設施數(shù) nJ 3; % 客戶數(shù) c [50; 60]; % 建設成本 cap [10; 10]; % 設施容量 d [1 2 3; 2 1 2]; % 單位運費d(i,j) U0 [4; 5; 6]; % 需求基準值 Delta [2; 2; 2]; % 需求波動范圍 % 盒式不確定集頂點枚舉2^3 8 個節(jié)點 vertices []; for u1 [U0(1)-Delta(1), U0(1)Delta(1)] for u2 [U0(2)-Delta(2), U0(2)Delta(2)] for u3 [U0(3)-Delta(3), U0(3)Delta(3)] vertices [vertices, [u1; u2; u3]]; end end end U_scen U0; % 初始場景用基準值 LB -inf; UB inf; tol 1e-3; Kmax 20; for k 1:Kmax % 主問題 MP x binvar(nI, 1); % 是否建設施 theta sdpvar(1, 1); % 第二階段成本上界 cons [sum(x) 1]; % 至少建一個設施 obj c*x theta; for s 1:size(U_scen, 2) % 每個已發(fā)現(xiàn)場景加一組變量 y sdpvar(nJ, nI, full); cons [cons, sum(y, 2) U_scen(:, s)]; % 需求等式 cons [cons, sum(y, 1) cap .* x]; % 容量上限 cons [cons, y 0]; cons [cons, theta sum(sum(d .* y))]; % θ吃掉運費 end optimize(cons, obj, sdpsettings(solver, gurobi, verbose, 0)); x_val value(x); theta_val value(theta); LB c*x_val theta_val; % 下界主問題最優(yōu)值 % 子問題 SP固定 x找最壞需求 best_sp -inf; u_star U0; for vidx 1:size(vertices, 2) u vertices(:, vidx); y sdpvar(nJ, nI, full); cons_sp [sum(y, 2) u; sum(y, 1) cap .* x_val; y 0]; optimize(cons_sp, sum(sum(d .* y)), ... sdpsettings(solver, gurobi, verbose, 0)); sp_val value(sum(sum(d .* y))); if sp_val best_sp best_sp sp_val; u_star u; end end UB min(UB, c*x_val best_sp); % 上界不斷收緊 fprintf(iter %2d, LB%.4f, UB%.4f, gap%.6f\n, ... k, LB, UB, (UB-LB)/UB); if (UB - LB) / abs(UB) tol break; end if ismember(U_scen, u_star, rows) break; % 場景重復已經(jīng)收斂 end U_scen [U_scen, u_star]; end這段代碼里最需要理解的是主問題里的θ約束。θ被設計成“所有已發(fā)現(xiàn)場景下運費的最大值”所以θ≥d?y^k對每個場景k都要成立MP的目標c?xθ就是在最小化“建設成本最壞運費”。子問題則是固定x_val后在盒式不確定集頂點上枚舉找出讓最小運費最大的那個需求場景u*這個u*下一輪會被追加進主問題。參數(shù)里tol控制收斂精度Kmax是最大迭代次數(shù)實際跑的時候建議先看gap變化曲線再決定要不要加大Kmax。一個很容易忽略的細節(jié)solver我寫的是gurobi如果本機沒有裝YALMIP會自動找可用求解器但binvar變量必須有MILP求解器才能處理。MATLAB自帶的intlinprog也能跑把sdpsettings里的solver改成intlinprog即可代價是求解速度慢一些模型調(diào)通后建議換商用求解器。3. Benders分解落地對偶割的推導與迭代實現(xiàn)3.1 從max-min到對偶問題Benders割是怎么來的Benders分解不像CCG那樣往主問題里塞整個場景的第二階段變量而是把第二階段問題取對偶形成一條只含第一階段變量x和對偶變量的線性約束。這條約束叫Benders最優(yōu)性割它的作用是在主問題里給θ一個不斷收緊的下界。核心推導是這樣的對固定的x和u第二階段內(nèi)層min問題是一個線性規(guī)劃取對偶后目標函數(shù)變成λ?u - μ?(cap ⊙ x)其中λ是需求等式的對偶變量μ是容量不等式≤約束的對偶變量。因為對偶變量只依賴對偶約束而不依賴x所以對任何新的x(λ*, μ*)仍然可行于是 λ?u- μ?(cap ⊙ x) 始終是第二階段最優(yōu)值的一個下界。把這個下界做成主問題里θ的約束就得到了Benders割。每一輪從子問題里撈出一組(λ, μ*, u*)就能給主問題加一條比上一輪更緊的割。3.2 Benders主問題與割累積一套比CCG更省的MPBenders的主問題里只有x和θ外加一堆累積的割不等式變量規(guī)模遠小于CCG的主問題這是它處理大規(guī)模問題的本錢。下面是同一選址問題的Benders循環(huán)SP部分仍然枚舉盒式頂點但重點是從約束句柄里取出對偶變量。% demo_benders_2stage.m % Benders 分解主循環(huán)同一套設施選址問題 % 延續(xù)上一節(jié)的 nI/nJ/c/cap/d/U0/Delta/vertices 定義 x binvar(nI, 1); theta sdpvar(1, 1); cons0 [sum(x) 1, theta 0]; % θ 先給一個平凡下界 cuts []; % 累積割的 cell for k 1:Kmax cons cons0; for t 1:length(cuts) lam cuts(t).lam; u cuts(t).u; mu cuts(t).mu; % Benders 最優(yōu)性割θ ≥ λ*u - μ*(cap ⊙ x) cons [cons, theta lam*u - mu*(cap .* x)]; end optimize(cons, c*x theta, ... sdpsettings(solver, gurobi, verbose, 0)); x_val value(x); theta_val value(theta); LB c*x_val theta_val; % SP固定 x 枚舉頂點找最壞場景并取對偶 best_sp -inf; for vidx 1:size(vertices, 2) u vertices(:, vidx); y sdpvar(nJ, nI, full); cons_sp [sum(y, 2) u; sum(y, 1) cap .* x_val; y 0]; optimize(cons_sp, sum(sum(d .* y)), ... sdpsettings(solver, gurobi, verbose, 0)); sp_val value(sum(sum(d .* y))); if sp_val best_sp best_sp sp_val; u_star u; % 取出需求等式和容量約束的對偶變量 lam dual(cons_sp(1)); mu dual(cons_sp(2)); end end UB min(UB, c*x_val best_sp); cuts(end1) struct(lam, lam, u, u_star, mu, mu); fprintf(iter %2d, LB%.4f, UB%.4f, gap%.6f\n, ... k, LB, UB, (UB-LB)/abs(UB)); if (UB - LB) / abs(UB) tol break; end end注意dual的提取順序非常嚴格它只認你在cons_sp里寫約束時的排列順序cons_sp(1)是需求等式對應λcons_sp(2)是容量不等式對應μ。建模順序一變對偶變量就全錯了這是Benders實現(xiàn)里最容易翻車的點。另外Benders的主問題在不含任何割的時候θ沒有上界約束目標函數(shù)會把θ往負無窮推所以我在cons0里加了theta 0作為平凡的初始下界否則第一輪MP就可能異常。實際對比起來Benders的主問題變量少每輪只加一條約束計算很快但它的割相對“松散”經(jīng)常要迭代幾十輪才收斂而CCG往往十幾輪就停了。這也是為什么現(xiàn)在主流的兩階段魯棒文獻更偏愛CCG。4. CCG與Benders如何選收斂速度、割類型與模型規(guī)模的取舍4.1 兩種割的本質(zhì)差異幾何視角與文獻結論從幾何上看CCG和Benders都在做同一件事用一個逐步精細化的多面體去逼近第二階段值函數(shù)。Benders的割是從對偶空間投影下來的支撐超平面每輪加一個主問題始終是個小規(guī)模LP/MILPCCG則是直接把不確定場景對應的原始變量和約束整個投射進主問題相當于把值函數(shù)的關鍵點連同局部結構一起復制進來。文獻里的共識是CCG的收斂速度在迭代輪數(shù)上幾乎總是優(yōu)于Benders尤其是在第二階段子問題約束較多的模型中CCG往往用不到Benders一半的迭代次數(shù)就達到同樣gap。但代價是CCG主問題的規(guī)模隨場景數(shù)線性膨脹當?shù)诙A段變量本身很多時主問題會變得非常龐大。Benders恰好相反每輪割很廉價主問題長期保持精簡適合割數(shù)量能控制在可接受范圍的大規(guī)模問題。實踐里還發(fā)現(xiàn)一個規(guī)律當不確定集維數(shù)低、第二階段LP規(guī)模適中時Benders的累積求解時間可能反而占優(yōu)而當不確定集結構復雜、需要少輪次快速收斂時直接選CCG。4.2 適用場景對照表與選型建議對比項CCG列與約束生成Benders 分解割的內(nèi)容完整場景變量約束一條對偶線性不等式主問題規(guī)模增長快每輪加一組y和約束慢每輪只加一條割典型迭代輪數(shù)少通常520輪多常需數(shù)十輪子問題處理只需解內(nèi)層min找最壞場景必須取對偶變量生成割第二階段變量多主問題膨脹求解變慢割仍然是一維影響小實現(xiàn)與調(diào)試難度低直觀中對偶提取易出錯適合場景論文復現(xiàn)、中小規(guī)模、快速收斂大規(guī)模LP結構、MP必須精簡我自己的選型習慣是復現(xiàn)論文階段先寫CCG因為它的子問題只用求最優(yōu)值不需要對偶變量代碼出錯率低等模型定型、確需要推到很大規(guī)模時再改成Benders。反過來一開始就寫B(tài)enders容易在dual提取和割正負號上浪費大量時間。5. 復現(xiàn)避坑指南五個讓兩階段模型翻車的常見錯誤5.1 MP有界但gap長期不降現(xiàn)象主問題每輪都在正常求解LB和UB各自變化但gap卡在某個值附近震蕩達到設定迭代上限還沒收斂。原因通常是兩個一是子問題找到的u其實已經(jīng)出現(xiàn)過但主問題重復添加同一場景變量重復但θ并不變緊二是MP和SP的求解器容差設置不一致MP用默認1e-7SP用1e-6兩邊數(shù)值精度互相拉扯。解決方法是加場景去重判斷代碼里用ismember(U_scen, u_star, rows)檢查u是否已存在同時把tol設為1e-3到1e-4之間并在MP和SP的sdpsettings里顯式設置一樣的求解器容差。5.2 求解器報“No suitable solver”現(xiàn)象明明裝了MATLAB模型也建好了一運行就提示找不到可用的求解器或者提示MILP問題需要MILP求解器。原因YALMIP默認的求解器搜索列表里沒有能處理binvar的求解器很多人裝了gurobi卻沒在YALMIP里配置路徑或者本機只有l(wèi)inprog這類連續(xù)LP求解器。解決臨時改用MATLAB自帶intlinprogsdpsettings(solver,intlinprog)長期則是正確安裝并配置gurobi或cplex的YALMIP接口。注意intlinprog能解MILP但大規(guī)模時性能差距明顯調(diào)通后盡快換。5.3 SP無界或者MP給出荒謬的x現(xiàn)象SP求解器返回“problem is unbounded”或者主問題求出的x明顯不合理比如一個設施都不建導致第二階段無可行解。原因子問題里某些需求頂點在當前x下根本滿足不了運輸問題不可行對偶問題自然無界更深一層是主問題的容量約束沒有和第一階段x耦合好比如cap ⊙ x在某些維度上為0而對應需求又必須被滿足。解決在SP里加人工松弛變量并賦予大懲罰系數(shù)優(yōu)先保證SP永遠有界同時在MP里加總容量約束比如sum(cap .* x) max(sum(vertices))從根上杜絕不可行。5.4 dual讀出來是空的現(xiàn)象Benders代碼運行到dual(cons_sp(1))時報錯提示dual extraction不可用或者返回的變量里全是NaN。原因YALMIP不是對所有求解器都支持dual提取有些開源求解器或intlinprog不會回傳對偶信息另一個常見原因是約束句柄順序和建模順序不一致取到了錯誤的約束。解決先用check(cons_sp)確認對偶可行性更換支持對偶提取的求解器比如gurobi、cplex、mosek同時嚴格按建模順序記錄約束句柄最好在注釋里標明第幾個約束對應哪個對偶變量。5.5 把代碼套進自己的模型后維度對不上現(xiàn)象按照壓縮包里的模板改數(shù)據(jù)結果主問題報維度不匹配或者SP里sum(y, 2)和u的長度對不齊。原因資源包的代碼是面向原論文的特定模型結構客戶數(shù)、設施數(shù)、約束排列都寫死在腳本里直接替換數(shù)據(jù)不會自動適配新結構。解決把nI、nJ、cap、U0、Delta這些全部抽成參數(shù)結構體凡是尺寸相關的循環(huán)一律用size()動態(tài)獲取特別要注意需求等式sum(y,2)是按客戶維度匯總容量約束sum(y,1)是按設施維度匯總兩個方向的維度含義千萬別搞反。6. 參數(shù)調(diào)優(yōu)與收斂驗證把CCG腳本改造成你自己的模型拿到資源包之后直接改數(shù)據(jù)是最低效的用法。我一般會把主問題和子問題分別封裝成函數(shù)輸入一個data結構體輸出x和u*這樣換模型時只需要動數(shù)據(jù)定義和約束行。data.c、data.cap、data.d、data.U0、data.Delta放進去函數(shù)內(nèi)部用length(data.U0)獲取維度所有約束用循環(huán)生成而不是手寫固定下標。改造的第一個驗證點是讓子問題目標值和你手算的某個確定場景一致把不確定集固定成基準值U0跑一遍SP如果運費不等于手算結果說明約束方向或者sum維度有問題。收斂驗證方面我習慣在迭代循環(huán)里保存每一輪的LB、UB和x_val退出循環(huán)后畫三條曲線LB曲線、UB曲線、gap曲線。三張圖能直接看出算法是卡在割質(zhì)量上還是卡在數(shù)值容差上。如果LB曲線一直不動多半是主問題里θ的約束不夠緊檢查新場景是否真的加了進去如果UB曲線波動很大則是SP每輪找到的最壞場景不穩(wěn)定通常是頂點枚舉時漏了某些頂點或者不確定集定義有誤。另外有一個容易被忽視的參數(shù)初始場景的選擇。資源包代碼用U0做初始場景這在很多模型里足夠但如果不確定集的波動范圍很大從U0出發(fā)會浪費好幾輪迭代。我通常會同時嘗試全上界、全下界、基準值三個初始場景對比總迭代輪數(shù)選效果最好的作為默認值。這個操作就幾行代碼但對收斂速度的影響往往很顯著。最后說一個我自己踩過的坑CCG代碼調(diào)通后直接拿去跑大規(guī)模數(shù)據(jù)J增到15以上時頂點枚舉會爆炸2^15個頂點根本枚舉不完。遇到這種情況把盒式集改成預算不確定集SP不必枚舉頂點而是解一個帶排序邏輯的線性規(guī)劃或者直接對偶后按λ符號定界。從那以后我每次拿到新的兩階段魯棒問題都先確認不確定集結構再決定SP用枚舉還是解析法這個習慣幫我省掉了大量無效迭代。希望幫到你。本文還有配套的精品資源點擊獲取 兩階段魯棒優(yōu)化CCG與Benders代碼實戰(zhàn)：原理、實現(xiàn)與避坑
簡介這份代碼包圍繞兩階段魯棒優(yōu)化中的CCG列與約束生成及Benders分解算法面向使用MATLAB與YALMIP的研究者和工程師既適合快速復現(xiàn)原論文基線也便于按實際決策場景擴展改編。壓縮包共5個文件、體積1.44MB含3篇PDF文獻與2個MATLAB腳本文獻涵蓋CCG與Benders對比、兩階段魯棒優(yōu)化求解方法以及Benders分解算法綜述腳本則提供可運行的案例代碼實現(xiàn)主問題構建、切割生成與迭代求解的完整鏈路。資源已有2287人學習下載在魯棒優(yōu)化與分解算法主題上具備較好的實踐參考熱度。通過學習包內(nèi)代碼與配套文獻讀者可掌握在YALMIP中搭建兩階段主從問題框架、利用CCG/Benders生成最優(yōu)割并識別最壞場景的關鍵思路繼而將原論文方法遷移至新能源出力、供應鏈需求等含不確定參數(shù)的優(yōu)化問題中為后續(xù)算法對比和工程落地提供可復現(xiàn)基礎。1. 兩階段魯棒優(yōu)化的CCG與Benders代碼包拿到就能跑但先搞懂割長什么樣做過兩階段魯棒優(yōu)化的人應該都有體會模型寫出來不難難的是把CCG或Benders的迭代循環(huán)跑通并讓結果能復現(xiàn)原論文的數(shù)值。這套帶MATLABYALMIP源碼和對比文獻的壓縮包解決的就是“從論文公式到可運行代碼”這段最磨人的路。它不是一個教學Demo而是把兩階段魯棒優(yōu)化里兩個主流求解框架——列與約束生成CCG和Benders分解——做成了可以直接改參數(shù)、換模型、對結果的可擴展腳本。適合正在寫相關論文的研究生也適合要把確定性優(yōu)化改成魯棒模型的工程師。一個反直覺的結論先放在這里CCG和Benders看起來是兩條完全不同的技術路線但實際上它們共用同一個“主問題-子問題”迭代骨架差別只在“切割”的構造方式理解這一點改造代碼能少走一半彎路。2. CCG列與約束生成主問題-子問題迭代與YALMIP代碼骨架2.1 CCG為什么快列生成與約束生成的雙重含義兩階段魯棒優(yōu)化的一般形式是 min c?x max_{u∈U} min_{y∈F(x,u)} d?y第一階段變量x要現(xiàn)在拍板第二階段變量y要等不確定參數(shù)u實現(xiàn)后再調(diào)整。CCG的思路非常直接不要試圖遍歷連續(xù)的不確定集U而是維護一個“已發(fā)現(xiàn)的極端場景”集合每輪找到當前最壞的那個u*把它對應的第二階段變量和約束整組塞進主問題。這里有個關鍵點CCG同時在做列生成和約束生成。主問題每加入一個新場景u*就會新增一組第二階段變量y^k這是“列”同時新增這組變量對應的約束需求滿足、容量限制、與θ的連線這是“約束”。因為不確定集的極值一定發(fā)生在頂點上線性目標函數(shù)在凸多面體上取極值必在頂點所以每輪找出的u*都是一個個頂點主問題不斷增強對“最壞情況”的描述。相比傳統(tǒng)Benders只往主問題丟一條對偶割CCG犧牲了主問題的規(guī)模換來了更緊的下界通常迭代次數(shù)明顯更少。2.2 CCG主問題與最壞場景搜索可直接修改的YALMIP實現(xiàn)下面用一個小規(guī)模的設施選址加不確定需求運輸問題演示CCG主循環(huán)。問題設定是I個候選設施建設成本c?容量cap?J個客戶需求u?在基準值上下波動第二階段決策是從設施i到客戶j的運輸量y??單位成本d??。第二階段目標是最小化運費整體目標是在最壞需求下最小化“建設成本運費”。% demo_ccg_2stage.m % 兩階段魯棒優(yōu)化CCG 主循環(huán)YALMIP MILP 求解器 % 模型設施選址 不確定需求運輸 % 第一階段0/1 決策是否建設施 % 第二階段滿足需求的最小運費 nI 2; % 設施數(shù) nJ 3; % 客戶數(shù) c [50; 60]; % 建設成本 cap [10; 10]; % 設施容量 d [1 2 3; 2 1 2]; % 單位運費d(i,j) U0 [4; 5; 6]; % 需求基準值 Delta [2; 2; 2]; % 需求波動范圍 % 盒式不確定集頂點枚舉2^3 8 個節(jié)點 vertices []; for u1 [U0(1)-Delta(1), U0(1)Delta(1)] for u2 [U0(2)-Delta(2), U0(2)Delta(2)] for u3 [U0(3)-Delta(3), U0(3)Delta(3)] vertices [vertices, [u1; u2; u3]]; end end end U_scen U0; % 初始場景用基準值 LB -inf; UB inf; tol 1e-3; Kmax 20; for k 1:Kmax % 主問題 MP x binvar(nI, 1); % 是否建設施 theta sdpvar(1, 1); % 第二階段成本上界 cons [sum(x) 1]; % 至少建一個設施 obj c*x theta; for s 1:size(U_scen, 2) % 每個已發(fā)現(xiàn)場景加一組變量 y sdpvar(nJ, nI, full); cons [cons, sum(y, 2) U_scen(:, s)]; % 需求等式 cons [cons, sum(y, 1) cap .* x]; % 容量上限 cons [cons, y 0]; cons [cons, theta sum(sum(d .* y))]; % θ吃掉運費 end optimize(cons, obj, sdpsettings(solver, gurobi, verbose, 0)); x_val value(x); theta_val value(theta); LB c*x_val theta_val; % 下界主問題最優(yōu)值 % 子問題 SP固定 x找最壞需求 best_sp -inf; u_star U0; for vidx 1:size(vertices, 2) u vertices(:, vidx); y sdpvar(nJ, nI, full); cons_sp [sum(y, 2) u; sum(y, 1) cap .* x_val; y 0]; optimize(cons_sp, sum(sum(d .* y)), ... sdpsettings(solver, gurobi, verbose, 0)); sp_val value(sum(sum(d .* y))); if sp_val best_sp best_sp sp_val; u_star u; end end UB min(UB, c*x_val best_sp); % 上界不斷收緊 fprintf(iter %2d, LB%.4f, UB%.4f, gap%.6f\n, ... k, LB, UB, (UB-LB)/UB); if (UB - LB) / abs(UB) tol break; end if ismember(U_scen, u_star, rows) break; % 場景重復已經(jīng)收斂 end U_scen [U_scen, u_star]; end這段代碼里最需要理解的是主問題里的θ約束。θ被設計成“所有已發(fā)現(xiàn)場景下運費的最大值”所以θ≥d?y^k對每個場景k都要成立MP的目標c?xθ就是在最小化“建設成本最壞運費”。子問題則是固定x_val后在盒式不確定集頂點上枚舉找出讓最小運費最大的那個需求場景u*這個u*下一輪會被追加進主問題。參數(shù)里tol控制收斂精度Kmax是最大迭代次數(shù)實際跑的時候建議先看gap變化曲線再決定要不要加大Kmax。一個很容易忽略的細節(jié)solver我寫的是gurobi如果本機沒有裝YALMIP會自動找可用求解器但binvar變量必須有MILP求解器才能處理。MATLAB自帶的intlinprog也能跑把sdpsettings里的solver改成intlinprog即可代價是求解速度慢一些模型調(diào)通后建議換商用求解器。3. Benders分解落地對偶割的推導與迭代實現(xiàn)3.1 從max-min到對偶問題Benders割是怎么來的Benders分解不像CCG那樣往主問題里塞整個場景的第二階段變量而是把第二階段問題取對偶形成一條只含第一階段變量x和對偶變量的線性約束。這條約束叫Benders最優(yōu)性割它的作用是在主問題里給θ一個不斷收緊的下界。核心推導是這樣的對固定的x和u第二階段內(nèi)層min問題是一個線性規(guī)劃取對偶后目標函數(shù)變成λ?u - μ?(cap ⊙ x)其中λ是需求等式的對偶變量μ是容量不等式≤約束的對偶變量。因為對偶變量只依賴對偶約束而不依賴x所以對任何新的x(λ*, μ*)仍然可行于是 λ?u- μ?(cap ⊙ x) 始終是第二階段最優(yōu)值的一個下界。把這個下界做成主問題里θ的約束就得到了Benders割。每一輪從子問題里撈出一組(λ, μ*, u*)就能給主問題加一條比上一輪更緊的割。3.2 Benders主問題與割累積一套比CCG更省的MPBenders的主問題里只有x和θ外加一堆累積的割不等式變量規(guī)模遠小于CCG的主問題這是它處理大規(guī)模問題的本錢。下面是同一選址問題的Benders循環(huán)SP部分仍然枚舉盒式頂點但重點是從約束句柄里取出對偶變量。% demo_benders_2stage.m % Benders 分解主循環(huán)同一套設施選址問題 % 延續(xù)上一節(jié)的 nI/nJ/c/cap/d/U0/Delta/vertices 定義 x binvar(nI, 1); theta sdpvar(1, 1); cons0 [sum(x) 1, theta 0]; % θ 先給一個平凡下界 cuts []; % 累積割的 cell for k 1:Kmax cons cons0; for t 1:length(cuts) lam cuts(t).lam; u cuts(t).u; mu cuts(t).mu; % Benders 最優(yōu)性割θ ≥ λ*u - μ*(cap ⊙ x) cons [cons, theta lam*u - mu*(cap .* x)]; end optimize(cons, c*x theta, ... sdpsettings(solver, gurobi, verbose, 0)); x_val value(x); theta_val value(theta); LB c*x_val theta_val; % SP固定 x 枚舉頂點找最壞場景并取對偶 best_sp -inf; for vidx 1:size(vertices, 2) u vertices(:, vidx); y sdpvar(nJ, nI, full); cons_sp [sum(y, 2) u; sum(y, 1) cap .* x_val; y 0]; optimize(cons_sp, sum(sum(d .* y)), ... sdpsettings(solver, gurobi, verbose, 0)); sp_val value(sum(sum(d .* y))); if sp_val best_sp best_sp sp_val; u_star u; % 取出需求等式和容量約束的對偶變量 lam dual(cons_sp(1)); mu dual(cons_sp(2)); end end UB min(UB, c*x_val best_sp); cuts(end1) struct(lam, lam, u, u_star, mu, mu); fprintf(iter %2d, LB%.4f, UB%.4f, gap%.6f\n, ... k, LB, UB, (UB-LB)/abs(UB)); if (UB - LB) / abs(UB) tol break; end end注意dual的提取順序非常嚴格它只認你在cons_sp里寫約束時的排列順序cons_sp(1)是需求等式對應λcons_sp(2)是容量不等式對應μ。建模順序一變對偶變量就全錯了這是Benders實現(xiàn)里最容易翻車的點。另外Benders的主問題在不含任何割的時候θ沒有上界約束目標函數(shù)會把θ往負無窮推所以我在cons0里加了theta 0作為平凡的初始下界否則第一輪MP就可能異常。實際對比起來Benders的主問題變量少每輪只加一條約束計算很快但它的割相對“松散”經(jīng)常要迭代幾十輪才收斂而CCG往往十幾輪就停了。這也是為什么現(xiàn)在主流的兩階段魯棒文獻更偏愛CCG。4. CCG與Benders如何選收斂速度、割類型與模型規(guī)模的取舍4.1 兩種割的本質(zhì)差異幾何視角與文獻結論從幾何上看CCG和Benders都在做同一件事用一個逐步精細化的多面體去逼近第二階段值函數(shù)。Benders的割是從對偶空間投影下來的支撐超平面每輪加一個主問題始終是個小規(guī)模LP/MILPCCG則是直接把不確定場景對應的原始變量和約束整個投射進主問題相當于把值函數(shù)的關鍵點連同局部結構一起復制進來。文獻里的共識是CCG的收斂速度在迭代輪數(shù)上幾乎總是優(yōu)于Benders尤其是在第二階段子問題約束較多的模型中CCG往往用不到Benders一半的迭代次數(shù)就達到同樣gap。但代價是CCG主問題的規(guī)模隨場景數(shù)線性膨脹當?shù)诙A段變量本身很多時主問題會變得非常龐大。Benders恰好相反每輪割很廉價主問題長期保持精簡適合割數(shù)量能控制在可接受范圍的大規(guī)模問題。實踐里還發(fā)現(xiàn)一個規(guī)律當不確定集維數(shù)低、第二階段LP規(guī)模適中時Benders的累積求解時間可能反而占優(yōu)而當不確定集結構復雜、需要少輪次快速收斂時直接選CCG。4.2 適用場景對照表與選型建議對比項CCG列與約束生成Benders 分解割的內(nèi)容完整場景變量約束一條對偶線性不等式主問題規(guī)模增長快每輪加一組y和約束慢每輪只加一條割典型迭代輪數(shù)少通常520輪多常需數(shù)十輪子問題處理只需解內(nèi)層min找最壞場景必須取對偶變量生成割第二階段變量多主問題膨脹求解變慢割仍然是一維影響小實現(xiàn)與調(diào)試難度低直觀中對偶提取易出錯適合場景論文復現(xiàn)、中小規(guī)模、快速收斂大規(guī)模LP結構、MP必須精簡我自己的選型習慣是復現(xiàn)論文階段先寫CCG因為它的子問題只用求最優(yōu)值不需要對偶變量代碼出錯率低等模型定型、確需要推到很大規(guī)模時再改成Benders。反過來一開始就寫B(tài)enders容易在dual提取和割正負號上浪費大量時間。5. 復現(xiàn)避坑指南五個讓兩階段模型翻車的常見錯誤5.1 MP有界但gap長期不降現(xiàn)象主問題每輪都在正常求解LB和UB各自變化但gap卡在某個值附近震蕩達到設定迭代上限還沒收斂。原因通常是兩個一是子問題找到的u其實已經(jīng)出現(xiàn)過但主問題重復添加同一場景變量重復但θ并不變緊二是MP和SP的求解器容差設置不一致MP用默認1e-7SP用1e-6兩邊數(shù)值精度互相拉扯。解決方法是加場景去重判斷代碼里用ismember(U_scen, u_star, rows)檢查u是否已存在同時把tol設為1e-3到1e-4之間并在MP和SP的sdpsettings里顯式設置一樣的求解器容差。5.2 求解器報“No suitable solver”現(xiàn)象明明裝了MATLAB模型也建好了一運行就提示找不到可用的求解器或者提示MILP問題需要MILP求解器。原因YALMIP默認的求解器搜索列表里沒有能處理binvar的求解器很多人裝了gurobi卻沒在YALMIP里配置路徑或者本機只有l(wèi)inprog這類連續(xù)LP求解器。解決臨時改用MATLAB自帶intlinprogsdpsettings(solver,intlinprog)長期則是正確安裝并配置gurobi或cplex的YALMIP接口。注意intlinprog能解MILP但大規(guī)模時性能差距明顯調(diào)通后盡快換。5.3 SP無界或者MP給出荒謬的x現(xiàn)象SP求解器返回“problem is unbounded”或者主問題求出的x明顯不合理比如一個設施都不建導致第二階段無可行解。原因子問題里某些需求頂點在當前x下根本滿足不了運輸問題不可行對偶問題自然無界更深一層是主問題的容量約束沒有和第一階段x耦合好比如cap ⊙ x在某些維度上為0而對應需求又必須被滿足。解決在SP里加人工松弛變量并賦予大懲罰系數(shù)優(yōu)先保證SP永遠有界同時在MP里加總容量約束比如sum(cap .* x) max(sum(vertices))從根上杜絕不可行。5.4 dual讀出來是空的現(xiàn)象Benders代碼運行到dual(cons_sp(1))時報錯提示dual extraction不可用或者返回的變量里全是NaN。原因YALMIP不是對所有求解器都支持dual提取有些開源求解器或intlinprog不會回傳對偶信息另一個常見原因是約束句柄順序和建模順序不一致取到了錯誤的約束。解決先用check(cons_sp)確認對偶可行性更換支持對偶提取的求解器比如gurobi、cplex、mosek同時嚴格按建模順序記錄約束句柄最好在注釋里標明第幾個約束對應哪個對偶變量。5.5 把代碼套進自己的模型后維度對不上現(xiàn)象按照壓縮包里的模板改數(shù)據(jù)結果主問題報維度不匹配或者SP里sum(y, 2)和u的長度對不齊。原因資源包的代碼是面向原論文的特定模型結構客戶數(shù)、設施數(shù)、約束排列都寫死在腳本里直接替換數(shù)據(jù)不會自動適配新結構。解決把nI、nJ、cap、U0、Delta這些全部抽成參數(shù)結構體凡是尺寸相關的循環(huán)一律用size()動態(tài)獲取特別要注意需求等式sum(y,2)是按客戶維度匯總容量約束sum(y,1)是按設施維度匯總兩個方向的維度含義千萬別搞反。6. 參數(shù)調(diào)優(yōu)與收斂驗證把CCG腳本改造成你自己的模型拿到資源包之后直接改數(shù)據(jù)是最低效的用法。我一般會把主問題和子問題分別封裝成函數(shù)輸入一個data結構體輸出x和u*這樣換模型時只需要動數(shù)據(jù)定義和約束行。data.c、data.cap、data.d、data.U0、data.Delta放進去函數(shù)內(nèi)部用length(data.U0)獲取維度所有約束用循環(huán)生成而不是手寫固定下標。改造的第一個驗證點是讓子問題目標值和你手算的某個確定場景一致把不確定集固定成基準值U0跑一遍SP如果運費不等于手算結果說明約束方向或者sum維度有問題。收斂驗證方面我習慣在迭代循環(huán)里保存每一輪的LB、UB和x_val退出循環(huán)后畫三條曲線LB曲線、UB曲線、gap曲線。三張圖能直接看出算法是卡在割質(zhì)量上還是卡在數(shù)值容差上。如果LB曲線一直不動多半是主問題里θ的約束不夠緊檢查新場景是否真的加了進去如果UB曲線波動很大則是SP每輪找到的最壞場景不穩(wěn)定通常是頂點枚舉時漏了某些頂點或者不確定集定義有誤。另外有一個容易被忽視的參數(shù)初始場景的選擇。資源包代碼用U0做初始場景這在很多模型里足夠但如果不確定集的波動范圍很大從U0出發(fā)會浪費好幾輪迭代。我通常會同時嘗試全上界、全下界、基準值三個初始場景對比總迭代輪數(shù)選效果最好的作為默認值。這個操作就幾行代碼但對收斂速度的影響往往很顯著。最后說一個我自己踩過的坑CCG代碼調(diào)通后直接拿去跑大規(guī)模數(shù)據(jù)J增到15以上時頂點枚舉會爆炸2^15個頂點根本枚舉不完。遇到這種情況把盒式集改成預算不確定集SP不必枚舉頂點而是解一個帶排序邏輯的線性規(guī)劃或者直接對偶后按λ符號定界。從那以后我每次拿到新的兩階段魯棒問題都先確認不確定集結構再決定SP用枚舉還是解析法這個習慣幫我省掉了大量無效迭代。希望幫到你。本文還有配套的精品資源點擊獲取