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