多階段魯棒調(diào)度模型MATLAB復(fù)現(xiàn)與CCG算法實(shí)踐)
做微電網(wǎng)調(diào)度的復(fù)現(xiàn)我印象最深的一課是在MATLAB里把確定性經(jīng)濟(jì)調(diào)度跑通之后信心滿滿地把風(fēng)光預(yù)測誤差從±5%改到±30%結(jié)果同一套計(jì)劃直接讓儲能越限、功率平衡失配。從那之后我就明白只要系統(tǒng)里裝了可再生能源和儲能調(diào)度模型里不考慮不確定性的最優(yōu)大概率只是某個理想場景下的最優(yōu)。這篇內(nèi)容圍繞的就是一套含可再生能源和儲能的區(qū)域微電網(wǎng)多階段魯棒調(diào)度模型基于MATLAB完整復(fù)現(xiàn)核心思路是建立日前計(jì)劃日內(nèi)調(diào)整的兩階段魯棒框架用區(qū)間集合刻畫風(fēng)光和負(fù)荷的預(yù)測偏差再通過列和約束生成算法CCG迭代求解最壞場景下的最優(yōu)運(yùn)行方案。文章會從模型為什么這樣設(shè)計(jì)、不確定性集合怎么構(gòu)建、CCG算法如何落地到MATLABYALMIPGurobi的代碼骨架和復(fù)現(xiàn)時(shí)容易踩的坑完整過一遍適合正在做微電網(wǎng)優(yōu)化調(diào)度、儲能配置或魯棒優(yōu)化應(yīng)用研究的朋友參考。1. 可再生滲透率攀升后傳統(tǒng)確定性調(diào)度為什么越來越不好用1.1 確定性調(diào)度的假設(shè)與失效場景傳統(tǒng)的微電網(wǎng)經(jīng)濟(jì)調(diào)度模型無論是目標(biāo)函數(shù)還是約束條件都默認(rèn)風(fēng)電、光伏和負(fù)荷在某個時(shí)段的出力是一個已知的確定值。模型看起來非常干凈目標(biāo)函數(shù)是最小化運(yùn)行成本包括向上級電網(wǎng)購電費(fèi)用、微型燃?xì)廨啓C(jī)的燃料費(fèi)用、儲能損耗成本約束包括功率平衡、儲能SOC上下限、充放電功率限制、機(jī)組出力限制、聯(lián)絡(luò)線功率限制。這種模型在可再生能源滲透率比較低、系統(tǒng)備用充足的時(shí)候還能湊合用。但一旦風(fēng)機(jī)和光伏在微電網(wǎng)里的裝機(jī)占比上來了問題就暴露了。風(fēng)電預(yù)測誤差通常在15%到25%之間波動光伏在天氣驟變時(shí)誤差可以超過30%負(fù)荷預(yù)測雖然相對準(zhǔn)但在極端天氣下也會出現(xiàn)明顯偏離。確定性調(diào)度把所有預(yù)測值當(dāng)成真實(shí)值來處理等于把全部風(fēng)險(xiǎn)都壓在了預(yù)測完全準(zhǔn)確這個脆弱的假設(shè)上。我做過一個對比測試在一個24節(jié)點(diǎn)微電網(wǎng)測試系統(tǒng)里用確定性模型求出一組日前計(jì)劃然后把實(shí)際風(fēng)光出力設(shè)定為預(yù)測值的70%和130%兩種偏差場景重新跑功率平衡校驗(yàn)結(jié)果有多組約束直接違反儲能SOC甚至越過了下限。這就是典型的計(jì)劃趕不上變化。1.2 隨機(jī)優(yōu)化、機(jī)會約束與魯棒優(yōu)化三條技術(shù)路線的實(shí)際取舍面對不確定性學(xué)術(shù)界和工程界常用的解決方案有三類很多初學(xué)者一上來容易混在一起方法基本思路優(yōu)點(diǎn)主要問題隨機(jī)規(guī)劃用大量場景近似不確定性對場景求期望最優(yōu)結(jié)果貼合概率分布經(jīng)濟(jì)性較好需要精確的概率分布場景數(shù)量大時(shí)計(jì)算量爆炸分布估計(jì)不準(zhǔn)時(shí)結(jié)果不可靠機(jī)會約束規(guī)劃允許約束以較高概率滿足如95%概率滿足功率平衡可以在安全性和經(jīng)濟(jì)性之間折中概率參數(shù)如何標(biāo)定比較主觀分布信息不足時(shí)難以求解魯棒優(yōu)化假設(shè)不確定性落在某個集合內(nèi)保證最壞情況下約束也滿足不需要精確概率分布結(jié)果對極端場景穩(wěn)健計(jì)算效率相對可控結(jié)果偏保守成本一般高于隨機(jī)規(guī)劃在微電網(wǎng)這種規(guī)模不大、但安全約束要求高的場景下魯棒優(yōu)化是比較務(wù)實(shí)的選擇。原因也很直接你通常拿不到足夠多的歷史數(shù)據(jù)來精確擬合風(fēng)電和光伏誤差的分布函數(shù)就算擬合出來分布參數(shù)本身也有估計(jì)誤差。與其在概率上賭一個精度不如直接說我知道誤差不會超過某個范圍在這個范圍內(nèi)怎么折騰我的調(diào)度方案都兜得住。1.3 標(biāo)題里魯棒性和不確定性到底指什么標(biāo)題里的魯棒性英文對應(yīng)Robustness放在調(diào)度模型里指的就是當(dāng)實(shí)際的風(fēng)電、光伏出力偏離預(yù)測值時(shí)系統(tǒng)依然能通過儲能充放電調(diào)整、機(jī)組出力再調(diào)度、與主網(wǎng)功率交換等手段維持功率平衡并滿足所有安全約束。衡量魯棒性的標(biāo)準(zhǔn)不是平均工況下的成本而是最壞工況下的約束滿足能力。不確定性在這個問題里主要來自三個方向源端不確定性風(fēng)機(jī)出力和光伏出力受氣象條件影響預(yù)測模型再準(zhǔn)也有偏差荷端不確定性區(qū)域負(fù)荷受生產(chǎn)生活規(guī)律、天氣溫度等因素影響同樣存在預(yù)測誤差設(shè)備狀態(tài)不確定性比如某臺機(jī)組臨時(shí)故障、某條線路檢修但在調(diào)度主問題里通常把這類視為N-1預(yù)想故障單獨(dú)處理不會混進(jìn)源荷不確定性里。在工程復(fù)現(xiàn)時(shí)通常把源荷不確定性統(tǒng)一建模為一個不確定集合U集合里的每一個元素代表一組可能的源荷出力偏差向量。魯棒優(yōu)化要求對U里的任意一個偏差向量系統(tǒng)都能通過第二階段調(diào)整找到可行解。這就是后面整個模型的核心出發(fā)點(diǎn)。2. 多階段魯棒調(diào)度模型的數(shù)學(xué)框架與不確定性集合設(shè)計(jì)2.1 兩階段模型的前后配合日前決策與日內(nèi)調(diào)整標(biāo)題里寫的是多階段魯棒調(diào)度模型實(shí)際復(fù)現(xiàn)時(shí)最常見的落地形式是兩階段魯棒優(yōu)化模型也就是經(jīng)典的三層結(jié)構(gòu)min-max-min。之所以叫多階段是因?yàn)檎{(diào)度決策本身是分層分時(shí)段的第一階段日前決策階段在獲知次日風(fēng)光和負(fù)荷預(yù)測值之后提前決定機(jī)組的啟停狀態(tài)、儲能充放電基準(zhǔn)計(jì)劃、向上級電網(wǎng)購電的日前計(jì)劃。這些決策的特點(diǎn)是必須在看到實(shí)際偏差之前做出所以也叫here-and-now決策。第二階段日內(nèi)調(diào)整階段當(dāng)實(shí)際風(fēng)光出力、負(fù)荷與預(yù)測出現(xiàn)偏差后在日前計(jì)劃的基礎(chǔ)上通過調(diào)節(jié)儲能出力、機(jī)組出力微調(diào)、可削減負(fù)荷等手段消除功率不平衡。這個階段的決策是在偏差揭曉之后才能做的叫wait-and-see決策。兩階段之間通過功率平衡約束和設(shè)備運(yùn)行約束耦合在一起。模型的目標(biāo)函數(shù)也相應(yīng)地拆成兩塊min (第一階段運(yùn)行成本) max_{u∈U} min_{y∈Ω(x,u)} (第二階段調(diào)整成本)外層min決定第一階段決策x內(nèi)層max尋找使第二階段成本最高的最壞偏差場景u中間層min則是在給定x和u之后尋找最優(yōu)的調(diào)整策略y。整個模型求的是在所有可能出現(xiàn)在的偏差場景下系統(tǒng)總成本的最壞情況最優(yōu)值。2.2 盒式、預(yù)算約束與橢球式不確定集合的適用場景不確定集合怎么建直接決定了模型的保守程度和求解難度。常見的有三種形式盒式集合Box直接把每個時(shí)段的不確定量限制在一個區(qū)間內(nèi)如u(t)∈[u_min(t), u_max(t)]。這種集合最簡單但要求所有時(shí)段同時(shí)達(dá)到最大偏差結(jié)果是最壞中的最壞過于保守。預(yù)算約束集合Budget在盒式集合基礎(chǔ)上增加一個總偏差比例限制形如Σ|u(t)-u_hat(t)|/δ(t) ≤ Γ。這里的Γ就是預(yù)算參數(shù)表示在調(diào)度周期內(nèi)最多允許多少個時(shí)段的偏差同時(shí)達(dá)到上限。Γ越小魯棒約束越寬松成本越低Γ越大方案越保守。這個設(shè)置非常貼合實(shí)際因?yàn)檎鎸?shí)的風(fēng)光出力偏差不太可能24個小時(shí)每個小時(shí)都頂在最大值上。橢球式集合Ellipsoidal用2-范數(shù)約束替代1-范數(shù)額外限制形如‖W(u-u_hat)‖? ≤ Ω。橢球集合考慮了變量之間的相關(guān)性數(shù)學(xué)上更精細(xì)但會引入二次約束計(jì)算復(fù)雜度明顯上升。在微電網(wǎng)調(diào)度里我最推薦的是盒式預(yù)算約束的組合。原因有三個物理意義清晰Γ可以理解為最壞情況下同時(shí)偏離預(yù)測的時(shí)間段數(shù)方便向工程人員解釋模型保持線性結(jié)構(gòu)經(jīng)過對偶變換后仍然是LP或MILP可以用YALMIPGurobi穩(wěn)定求解通過調(diào)節(jié)Γ的取值可以畫出從激進(jìn)到保守的完整成本曲線這是論文和報(bào)告里很有說服力的一張圖。2.3 目標(biāo)函數(shù)與約束的完整數(shù)學(xué)表達(dá)以含風(fēng)電、光伏、儲能和微型燃?xì)廨啓C(jī)的區(qū)域微電網(wǎng)為例完整的魯棒調(diào)度模型可以寫成決策變量第一階段燃?xì)廨啓C(jī)開機(jī)狀態(tài)z(t)、儲能充電功率P_ess_c(t)、儲能放電功率P_ess_d(t)、日前購電功率P_grid(t)、儲能SOC(t)。第二階段針對偏差場景u的儲能調(diào)整量ΔP_ess(t)、機(jī)組調(diào)整量ΔP_g(t)、可削減負(fù)荷量ΔP_load(t)。目標(biāo)函數(shù)min Σ_t [ C_grid * P_grid(t) C_g * P_g(t) C_ess * (P_ess_c(t)P_ess_d(t)) ] max_{u∈U} min_Δ Σ_t [ C_penalty * (ΔP_load(t) 棄風(fēng)棄光懲罰項(xiàng)) ]約束條件功率平衡約束 P_grid(t) P_wind(t) P_pv(t) P_g(t) P_ess_d(t) P_load(t) P_ess_c(t) ΔP_load(t)其中P_wind(t)和P_pv(t)在第二階段表示為P_wind_hat(t) u_w(t)、P_pv_hat(t) u_pv(t)。儲能約束 SOC(t1) SOC(t) η_c * P_ess_c(t) * Δt - P_ess_d(t) * Δt / η_d SOC_min ≤ SOC(t) ≤ SOC_max 0 ≤ P_ess_c(t) ≤ P_ess_c_max * z_ess(t) 0 ≤ P_ess_d(t) ≤ P_ess_d_max * (1 - z_ess(t))這里需要注意的是充放電狀態(tài)變量z_ess(t)的引入它保證儲能不會同時(shí)充電和放電這在后面代碼實(shí)現(xiàn)時(shí)是個關(guān)鍵細(xì)節(jié)。機(jī)組約束 z(t) * P_g_min ≤ P_g(t) ≤ z(t) * P_g_max 爬坡約束|P_g(t) - P_g(t-1)| ≤ Ramp_g與主網(wǎng)交互約束 P_grid_min ≤ P_grid(t) ≤ P_grid_max可削減負(fù)荷約束 0 ≤ ΔP_load(t) ≤ α * P_load(t)2.4 為什么多階段不等于多周期很多初學(xué)者會把多階段魯棒理解成把24小時(shí)的調(diào)度周期分成很多段每段做一個魯棒優(yōu)化這是一個明顯的誤區(qū)。多周期優(yōu)化只是把時(shí)間維度拉長決策變量在整個周期內(nèi)一次性規(guī)劃本質(zhì)上仍是開環(huán)調(diào)度。多階段魯棒優(yōu)化的關(guān)鍵不在于時(shí)間段的多少而在于決策的可調(diào)整性——第一階段決策做出后可以看到不確定性的實(shí)現(xiàn)然后第二階段再做出調(diào)整決策。這個特征對應(yīng)的是閉環(huán)調(diào)度也就是日前計(jì)劃日內(nèi)實(shí)時(shí)調(diào)整的真實(shí)運(yùn)行流程。用數(shù)學(xué)語言說真正的多階段魯棒模型必須具有min-max-min的嵌套結(jié)構(gòu)而不是簡單的min-max或者包含很多時(shí)段的一次性min。在復(fù)現(xiàn)時(shí)如果你的模型只是每個時(shí)段都預(yù)留備用容量然后用魯棒約束包住所有時(shí)段這雖然也算魯棒但不是多階段魯棒。要體現(xiàn)多階段就必須顯式區(qū)分可調(diào)整變量如儲能的實(shí)時(shí)充放電調(diào)整量和不可調(diào)整變量如機(jī)組啟停、日前購電計(jì)劃并在目標(biāo)函數(shù)中體現(xiàn)第二階段對最壞場景的應(yīng)對成本。3. 列和約束生成CCG算法迭代框架的收斂邏輯3.1 主問題的最小費(fèi)用加最壞場景切割線兩階段魯棒模型直接求解極其困難因?yàn)閮?nèi)層嵌套了max-min優(yōu)化。工程上最常用的解法是列和約束生成算法核心思路是把原問題拆成一個主問題MP和一個子問題SP通過迭代逐步逼近最優(yōu)解。主問題的形式如下min c * x η s.t. 第一階段約束 η ≥ d * y_l, 對于已知的場景u_l x、y_l之間的耦合約束這里的關(guān)鍵是引入了輔助變量η它表示第二階段成本的上界估計(jì)。每輪迭代中子問題會找到一個最壞場景u_l主問題就把這個場景對應(yīng)的第二階段變量y_l和約束作為切割加入模型讓η不斷收緊。主問題因?yàn)橹豢紤]了已發(fā)現(xiàn)的最壞場景所以實(shí)際上是原問題的松弛版本它的最優(yōu)值給出原問題的下界LB。3.2 子問題的max-min對偶轉(zhuǎn)化子問題是在給定第一階段解x_star之后尋找使第二階段成本最大的場景u以及對應(yīng)的最優(yōu)調(diào)整成本形式是f(x_star) max_{u∈U} min_{y∈Ω(x_star, u)} d * y這個max-min結(jié)構(gòu)不能直接丟給求解器。標(biāo)準(zhǔn)做法是對內(nèi)層的min問題取強(qiáng)對偶。由于內(nèi)層是關(guān)于y的線性規(guī)劃只要可行域非空且有界強(qiáng)對偶成立就可以把min問題等價(jià)轉(zhuǎn)化為其對偶max問題然后和外層的max合并變成一個單層max問題f(x_star) max_{u∈U, λ∈Λ} λ * (b - A * x_star - C * u)其中λ是對偶變量它的符號由原約束方向決定。如果原約束是≤型λ非負(fù)如果是等式約束λ自由變量。合并之后子問題變成帶二次項(xiàng)的非線性問題因?yàn)棣撕蛈相乘。在實(shí)際求解中這個雙線性項(xiàng)通?？梢杂么驧法或直接交給支持MIQP的非線性求解器處理Gurobi可以直接求解這類帶雙線性項(xiàng)的小規(guī)模問題。需要強(qiáng)調(diào)一點(diǎn)內(nèi)層min問題必須滿足強(qiáng)對偶條件。因此第二階段調(diào)整問題必須建模為線性規(guī)劃不能隨意加入整數(shù)變量。如果第二階段需要決策是否切除某條線路這類0-1變量標(biāo)準(zhǔn)CCG就不適用了需要改成嵌套CCG或者引入其他分解方法。3.3 收斂判據(jù)與最大迭代次數(shù)CCG算法的主循環(huán)流程如下初始化UBinf, LB-inf, k0設(shè)置初始場景u_0通常取預(yù)測均值 while k ≤ K_max 且 UB - LB tol: 1. 求解主問題得到x_k*和η_k*更新LB c*x_k* η_k* 2. 將x_k*代入子問題求解最壞場景u_k*和最壞場景成本f(x_k*) 3. 更新UB min(UB, c*x_k* f(x_k*)) 4. if UB - LB ≤ tol: 收斂退出 5. 否則將u_k*作為新場景加入主問題添加對應(yīng)的第二階段變量y_k和割約束 k k 1收斂判據(jù)一般取相對間隙或絕對間隙。我在復(fù)現(xiàn)時(shí)常用的組合是絕對誤差取1e-4最大迭代次數(shù)取20輪。在24時(shí)段微電網(wǎng)模型里通常3到8輪就能收斂到很緊的間隙很少會真的跑到20輪。3.4 CCG與Benders分解的選擇建議很多文獻(xiàn)會同時(shí)提到Benders分解和CCG兩者都是求解兩階段問題和分解主-子問題的經(jīng)典算法但適用場景有明顯差異對比維度Benders分解CCG切割信息通過對偶函數(shù)添加割平面直接添加最壞場景對應(yīng)的完整約束迭代收斂對小規(guī)模問題收斂較快通常需要更少迭代次數(shù)尤其有連續(xù)變量時(shí)對整數(shù)變量的適應(yīng)性一般需要處理整數(shù)對偶gap主問題中可以自然處理整數(shù)變量實(shí)現(xiàn)難度較低中等在微電網(wǎng)魯棒調(diào)度里第一階段通常包含機(jī)組啟停這類整數(shù)變量而CCG的切割是把具體場景u_l下的約束完整塞進(jìn)主問題對整數(shù)變量天然友好所以實(shí)際復(fù)現(xiàn)時(shí)我基本都用CCG。只有當(dāng)子問題規(guī)模特別大、每次求解SP耗時(shí)很高時(shí)才會回頭考慮Benders或者近似方法。4. MATLABYALMIPGurobi的關(guān)鍵代碼實(shí)現(xiàn)4.1 環(huán)境準(zhǔn)備與求解器適配復(fù)現(xiàn)前先確認(rèn)環(huán)境MATLAB版本R2020b及以上即可建議R2022b以上對sdpvar和optimizer的兼容性更好YALMIP直接從GitHub拉取最新release不要用十年前的老版本求解器Gurobi 10.x或11.x學(xué)術(shù)許可免費(fèi)申請Mosek 10.x也可以但Gurobi在MILP上的速度通常更優(yōu)。在MATLAB里先跑一句驗(yàn)證yalmiptest能列出Gurobi的LP、QP、MILP、MIQP支持狀態(tài)就說明環(huán)境沒問題。如果用了舊版MATLAB配新版Gurobi經(jīng)常會出現(xiàn)java相關(guān)報(bào)錯那是MATLAB版本太老、Java版本不匹配導(dǎo)致的屬于環(huán)境問題先升級MATLAB再排查。4.2 主問題的YALMIP建模代碼下面是主問題的核心建模骨架采用24時(shí)段、1小時(shí)為步長的設(shè)定。nT 24; % 第一階段變量 z_g binvar(1, nT); % 燃?xì)廨啓C(jī)啟停 P_g sdpvar(1, nT); % 燃?xì)廨啓C(jī)出力 P_ess_c sdpvar(1, nT); % 儲能充電功率 P_ess_d sdpvar(1, nT); % 儲能放電功率 z_ess binvar(1, nT); % 儲能充放電狀態(tài)1充電0放電 SOC sdpvar(1, nT1); % 荷電狀態(tài) P_grid sdpvar(1, nT); % 購電功率 eta sdpvar(1); % 最壞場景成本上界 % 目標(biāo)函數(shù) objective_mp sum(grid_price .* P_grid) ... sum(g_fuel_cost * P_g) ... sum(ess_cost * (P_ess_c P_ess_d)) ... eta; % 約束 constraints_mp []; constraints_mp [constraints_mp, SOC(1) SOC_init]; for t 1:nT constraints_mp [constraints_mp, ... SOC(t1) SOC(t) eta_charge * P_ess_c(t) - P_ess_d(t) / eta_discharge, ... SOC_min SOC(t1) SOC_max, ... 0 P_ess_c(t) P_ess_c_max * z_ess(t), ... 0 P_ess_d(t) P_ess_d_max * (1 - z_ess(t)), ... z_g(t) * P_g_min P_g(t) z_g(t) * P_g_max, ... P_grid_min P_grid(t) P_grid_max, ... P_wind_forecast(t) P_pv_forecast(t) P_g(t) P_grid(t) P_ess_d(t) ... P_load_forecast(t) P_ess_c(t)]; end ops sdpsettings(solver, gurobi, verbose, 0, gurobi.NumericFocus, 3); optimize(constraints_mp, objective_mp, ops);這里有個細(xì)節(jié)給Gurobi設(shè)置NumericFocus為3可以明顯降低大M約束帶來的數(shù)值病態(tài)問題。如果模型里用到了大M這個設(shè)置幾乎是必加的。4.3 子問題的對偶化與最壞場景求解代碼第二階段的最壞場景求解是核心代碼分三步走第一步固定第一階段變量第二步構(gòu)造第二階段LP并取對偶第三步求解單層max問題。% 固定第一階段變量 x_star value([P_g; P_grid; P_ess_c; P_ess_d; SOC]); % 不確定變量 u_w sdpvar(1, nT); % 風(fēng)電偏差 u_pv sdpvar(1, nT); % 光伏偏差 u_l sdpvar(1, nT); % 負(fù)荷偏差 % 第二階段調(diào)整變量 delta_ess sdpvar(1, nT); % 儲能功率調(diào)整 delta_g sdpvar(1, nT); % 機(jī)組出力調(diào)整 delta_load sdpvar(1, nT); % 可削減負(fù)荷 % 不確定集合盒式預(yù)算 U_con [-u_w_max u_w u_w_max, ... -u_pv_max u_pv u_pv_max, ... -u_l_max u_l u_l_max, ... sum(abs(u_w) ./ u_w_max) sum(abs(u_pv) ./ u_pv_max) Gamma]; % 第二階段約束僅對應(yīng)給定x_star SP_con [...]; SP_obj sum(penalty_load * delta_load penalty_curtail * (u_w u_pv) ... C_delta_g * abs(delta_g) C_delta_ess * abs(delta_ess)); % 通過YALMIP取對偶或直接用求解器求解最大最小問題 % 更通用的做法是手動構(gòu)造對偶變量把max-min轉(zhuǎn)為單層max dual_lambda sdpvar(size(A_sp, 1), 1); dual_mu sdpvar(size(A_eq, 1), 1); % 對偶目標(biāo)函數(shù)注意線性化雙線性項(xiàng) dual_obj dual_lambda * (b_sp - A_sp * x_star_vec - B_u * [u_w; u_pv; u_l]) ... dual_mu * (b_eq - A_eq * x_star_vec); % 合并后的最壞場景問題 SP_merged_con [dual_lambda 0, dual_mu free, U_con]; result_sp optimize(SP_merged_con, -dual_obj, ops_sp);實(shí)際上當(dāng)雙線性項(xiàng)存在時(shí)Gurobi會以MIQP方式求解變量數(shù)不大時(shí)速度可以接受。如果模型規(guī)模較大建議用big-M法把雙線性項(xiàng)線性化把問題轉(zhuǎn)成MILP。4.4 迭代主循環(huán)與結(jié)果保存主循環(huán)代碼框架如下UB inf; LB -inf; k 0; tol 1e-4; Kmax 20; cut_results {}; while k Kmax (UB - LB) tol k k 1; % 1. 求解主問題 optimize(constraints_mp, objective_mp, ops); LB value(objective_mp); x_cur value([P_g; P_grid; P_ess_c; P_ess_d; SOC]); % 2. 求解子問題 optimize(SP_merged_con, -dual_obj, ops_sp); sp_cost value(dual_obj); u_cur value([u_w; u_pv; u_l]); day_ahead_cost value(sum(grid_price .* P_grid) ... sum(g_fuel_cost * P_g) ... sum(ess_cost * (P_ess_c P_ess_d))); UB min(UB, day_ahead_cost sp_cost); fprintf(迭代%d: LB%.4f, UB%.4f, gap%.6f\n, k, LB, UB, UB-LB); % 3. 添加新場景切割到主問題 y_new sdpvar(size_y, 1); constraints_mp [constraints_mp, ... y_new對應(yīng)的第二階段約束把u_cur作為常數(shù)代入, ... eta 第二階段成本表達(dá)式(y_new)]; cut_results{k}.u u_cur; cut_results{k}.cost sp_cost; end有幾個實(shí)現(xiàn)上的細(xì)節(jié)值得說不要把主問題的sdpvar變量每輪重新創(chuàng)建。變量結(jié)構(gòu)要保持不變只是追加約束。否則MATLAB每輪都要重新解析整個符號表達(dá)式速度直接掉一個數(shù)量級。主問題的optimize結(jié)果如果返回problemsolver 信息先打印optimize的第二個輸出看看是primal infeasible還是numerical issues再決定調(diào)參數(shù)還是調(diào)模型。保存每輪的最壞場景u_cur最后畫圖時(shí)可以直接把最壞場景下的風(fēng)電光伏曲線和基準(zhǔn)預(yù)測曲線對比很直觀。4.5 代碼性能優(yōu)化的幾個細(xì)節(jié)復(fù)現(xiàn)完跑第一版10輪迭代可能要跑幾分鐘如果覺得慢按下面這幾個方向優(yōu)化通常能把時(shí)間壓到原來的三分之一用矩陣約束代替for循環(huán)。上面為了可讀性用了for循環(huán)實(shí)際寫代碼時(shí)盡量把所有時(shí)段的約束拼接成大矩陣一次性傳給constraints。YALMIP在拼接大規(guī)模符號矩陣時(shí)效率一般減少循環(huán)能明顯提升解析速度。子問題求解時(shí)不要用默認(rèn)的optimize改用optimizer對象把x_star作為參數(shù)傳入省去每輪重建約束的開銷。子問題的對偶變量數(shù)量通常和約束數(shù)量一致如果第二階段約束很多先做約束化簡把明顯冗余的約束刪掉。Gurobi的線程數(shù)設(shè)為4到8ops.gurobi.Threads默認(rèn)是0自動在共享計(jì)算環(huán)境里可以顯式設(shè)置。模型里大量使用abs()時(shí)YALMIP會引入輔助變量和不等式變量規(guī)模膨脹很快盡量在建模階段就用正負(fù)分解來避開abs。5. 完全復(fù)現(xiàn)中的參數(shù)設(shè)定、結(jié)果校核與常見坑5.1 推薦參考輸入?yún)?shù)參數(shù)設(shè)定直接影響收斂性和結(jié)果合理性下面這組參數(shù)是我在多個復(fù)現(xiàn)案例里校準(zhǔn)過的參考值可以直接作為起點(diǎn)參數(shù)推薦值說明調(diào)度周期24時(shí)段1h/時(shí)段多階段魯棒最基本的設(shè)定風(fēng)電裝機(jī)500 kW預(yù)測誤差取±15%光伏裝機(jī)300 kW預(yù)測誤差取±10%負(fù)荷峰值600 kW預(yù)測誤差取±5%儲能容量400 kWhSOC范圍20%~95%儲能最大充放電功率100 kW充放電效率均取0.95燃?xì)廨啓C(jī)容量200 kW最小出力50 kW爬坡率50kW/h購電價(jià)格上限1.2 元/kWh分時(shí)電價(jià)峰谷差越大越能體現(xiàn)儲能價(jià)值預(yù)算參數(shù)Γ3~5表示最多3~5個時(shí)段同時(shí)達(dá)到最大偏差收斂閾值1e-4絕對間隙最大迭代次數(shù)20實(shí)際通常8輪內(nèi)收斂5.2 怎么判斷你復(fù)現(xiàn)的模型確實(shí)魯棒跑通模型只是第一步關(guān)鍵是要會驗(yàn)證。我復(fù)現(xiàn)完之后一般做三組校驗(yàn)任何一組不過都說明模型有問題第一組最壞場景重放。把CCG最后一輪得到的最壞場景u_cur代入第二階段模型在固定第一階段決策的情況下重新求解第二階段調(diào)整問題。如果找不到可行解說明子問題建模或者切割添加有問題。第二組蒙特卡洛模擬。從風(fēng)電、光伏、負(fù)荷的預(yù)測誤差分布中隨機(jī)抽樣生成1000個場景對每個場景調(diào)用第二階段模型檢驗(yàn)是否可行。如果可行率低于98%說明不確定性集合設(shè)計(jì)偏激進(jìn)需要調(diào)大?；驍U(kuò)大盒式范圍。第三組Γ對比分析。分別取Γ0、1、3、6、12觀察總成本和最壞場景成本的變化。Γ0時(shí)模型退化為確定性問題總成本最低隨著Γ增大總成本單調(diào)上升這是符合預(yù)期的。如果曲線出現(xiàn)非單調(diào)或者劇烈跳變大概率是不確定集合的約束寫錯了。5.3 我實(shí)際踩過的幾個大坑第一個坑是儲能的SOC約束沒有處理好充放電同時(shí)性。模型里P_ess_c和P_ess_d是兩個獨(dú)立變量如果不加z_ess互斥約束求解器會出現(xiàn)一個時(shí)段既充電又放電的情況利用效率空轉(zhuǎn)降低成本看起來成本很漂亮實(shí)際上物理不可行。解決辦法就是代碼里用z_ess做互斥或者加入一個很小的同時(shí)充放電懲罰項(xiàng)。第二個坑是對偶變量符號搞錯。子問題在對偶化時(shí)如果原約束是Ax Bu ≤ b對偶變量必須是非負(fù)的如果是等式約束對偶變量自由。符號一錯對偶目標(biāo)的值會偏離真實(shí)最壞場景成本導(dǎo)致UB和LB無法收斂甚至發(fā)散。第三個坑是無腦擴(kuò)大M。在線性化雙線性項(xiàng)時(shí)M取值太小會剪掉有效解M取值太大會導(dǎo)致數(shù)值病態(tài)。解決方法是盡量用不確定變量的區(qū)間上下界來推導(dǎo)M的理論值而不是隨便取一個1e6。Gurobi的NumericFocus開啟后有一定緩解但不解決根本問題。第四個坑是YALMIP的符號展開。模型復(fù)雜時(shí)YALMIP默認(rèn)會把所有約束展開成完整系數(shù)矩陣變量多的時(shí)候內(nèi)存占用會爆炸。設(shè)置sdpsettings(expand, 0)可以減輕這個問題但代價(jià)是求解器內(nèi)部的presolve工作量增加需要實(shí)測權(quán)衡。5.4 復(fù)現(xiàn)后如何擴(kuò)展這個模型跑通基礎(chǔ)版之后有幾個很自然的擴(kuò)展方向都能在不推翻框架的前提下增加價(jià)值分布魯棒優(yōu)化把固定盒式集合換成Wasserstein球利用歷史數(shù)據(jù)構(gòu)建經(jīng)驗(yàn)分布再以球內(nèi)最壞分布為目標(biāo)。這個擴(kuò)展只需要改不確定集合的建模和對偶約束CCG框架可以復(fù)用。多能互補(bǔ)在電儲能的基礎(chǔ)上加入蓄熱罐、燃?xì)忮仩t把區(qū)域微電網(wǎng)擴(kuò)展成綜合能源系統(tǒng)。多能耦合會引入熱功率平衡約束但兩階段魯棒的結(jié)構(gòu)完全一致。碳約束在目標(biāo)函數(shù)中加入碳排放懲罰或碳配額約束分析碳交易價(jià)格對魯棒調(diào)度策略的影響。這類擴(kuò)展對接雙碳主題在實(shí)際匯報(bào)中更有說服力。儲能容量優(yōu)化把儲能的額定容量和功率作為第一階段決策變量目標(biāo)函數(shù)中加入儲能投資年化成本模型就變成了容量規(guī)劃魯棒調(diào)度聯(lián)合優(yōu)化問題CCG依然適用。我個人在多次復(fù)現(xiàn)這個模型后的體會是兩階段魯棒框架本身并不算特別難難的是把每個環(huán)節(jié)的物理意義和數(shù)學(xué)細(xì)節(jié)對齊。比如Γ怎么解釋給非優(yōu)化背景的同事聽比如儲能SOC方程里為什么充電效率在乘法位置、放電效率在除法位置這些細(xì)節(jié)如果不理清楚代碼復(fù)制過來改了參數(shù)結(jié)果不對勁都不知道去哪查。建議第一次復(fù)現(xiàn)的朋友先用小系統(tǒng)比如6節(jié)點(diǎn)、4時(shí)段把CCG迭代曲線畫出來確認(rèn)LB和UB的收斂趨勢沒問題了再放大到24時(shí)段完整系統(tǒng)。小系統(tǒng)里能一眼看出的邏輯錯誤放大以后往往會被數(shù)值噪音掩蓋排查難度翻倍。這個模型跑通之后你再去讀其他魯棒調(diào)度相關(guān)的論文思路會清晰很多。 微電網(wǎng)多階段魯棒調(diào)度模型MATLAB復(fù)現(xiàn)與CCG算法實(shí)踐
做微電網(wǎng)調(diào)度的復(fù)現(xiàn)我印象最深的一課是在MATLAB里把確定性經(jīng)濟(jì)調(diào)度跑通之后信心滿滿地把風(fēng)光預(yù)測誤差從±5%改到±30%結(jié)果同一套計(jì)劃直接讓儲能越限、功率平衡失配。從那之后我就明白只要系統(tǒng)里裝了可再生能源和儲能調(diào)度模型里不考慮不確定性的最優(yōu)大概率只是某個理想場景下的最優(yōu)。這篇內(nèi)容圍繞的就是一套含可再生能源和儲能的區(qū)域微電網(wǎng)多階段魯棒調(diào)度模型基于MATLAB完整復(fù)現(xiàn)核心思路是建立日前計(jì)劃日內(nèi)調(diào)整的兩階段魯棒框架用區(qū)間集合刻畫風(fēng)光和負(fù)荷的預(yù)測偏差再通過列和約束生成算法CCG迭代求解最壞場景下的最優(yōu)運(yùn)行方案。文章會從模型為什么這樣設(shè)計(jì)、不確定性集合怎么構(gòu)建、CCG算法如何落地到MATLABYALMIPGurobi的代碼骨架和復(fù)現(xiàn)時(shí)容易踩的坑完整過一遍適合正在做微電網(wǎng)優(yōu)化調(diào)度、儲能配置或魯棒優(yōu)化應(yīng)用研究的朋友參考。1. 可再生滲透率攀升后傳統(tǒng)確定性調(diào)度為什么越來越不好用1.1 確定性調(diào)度的假設(shè)與失效場景傳統(tǒng)的微電網(wǎng)經(jīng)濟(jì)調(diào)度模型無論是目標(biāo)函數(shù)還是約束條件都默認(rèn)風(fēng)電、光伏和負(fù)荷在某個時(shí)段的出力是一個已知的確定值。模型看起來非常干凈目標(biāo)函數(shù)是最小化運(yùn)行成本包括向上級電網(wǎng)購電費(fèi)用、微型燃?xì)廨啓C(jī)的燃料費(fèi)用、儲能損耗成本約束包括功率平衡、儲能SOC上下限、充放電功率限制、機(jī)組出力限制、聯(lián)絡(luò)線功率限制。這種模型在可再生能源滲透率比較低、系統(tǒng)備用充足的時(shí)候還能湊合用。但一旦風(fēng)機(jī)和光伏在微電網(wǎng)里的裝機(jī)占比上來了問題就暴露了。風(fēng)電預(yù)測誤差通常在15%到25%之間波動光伏在天氣驟變時(shí)誤差可以超過30%負(fù)荷預(yù)測雖然相對準(zhǔn)但在極端天氣下也會出現(xiàn)明顯偏離。確定性調(diào)度把所有預(yù)測值當(dāng)成真實(shí)值來處理等于把全部風(fēng)險(xiǎn)都壓在了預(yù)測完全準(zhǔn)確這個脆弱的假設(shè)上。我做過一個對比測試在一個24節(jié)點(diǎn)微電網(wǎng)測試系統(tǒng)里用確定性模型求出一組日前計(jì)劃然后把實(shí)際風(fēng)光出力設(shè)定為預(yù)測值的70%和130%兩種偏差場景重新跑功率平衡校驗(yàn)結(jié)果有多組約束直接違反儲能SOC甚至越過了下限。這就是典型的計(jì)劃趕不上變化。1.2 隨機(jī)優(yōu)化、機(jī)會約束與魯棒優(yōu)化三條技術(shù)路線的實(shí)際取舍面對不確定性學(xué)術(shù)界和工程界常用的解決方案有三類很多初學(xué)者一上來容易混在一起方法基本思路優(yōu)點(diǎn)主要問題隨機(jī)規(guī)劃用大量場景近似不確定性對場景求期望最優(yōu)結(jié)果貼合概率分布經(jīng)濟(jì)性較好需要精確的概率分布場景數(shù)量大時(shí)計(jì)算量爆炸分布估計(jì)不準(zhǔn)時(shí)結(jié)果不可靠機(jī)會約束規(guī)劃允許約束以較高概率滿足如95%概率滿足功率平衡可以在安全性和經(jīng)濟(jì)性之間折中概率參數(shù)如何標(biāo)定比較主觀分布信息不足時(shí)難以求解魯棒優(yōu)化假設(shè)不確定性落在某個集合內(nèi)保證最壞情況下約束也滿足不需要精確概率分布結(jié)果對極端場景穩(wěn)健計(jì)算效率相對可控結(jié)果偏保守成本一般高于隨機(jī)規(guī)劃在微電網(wǎng)這種規(guī)模不大、但安全約束要求高的場景下魯棒優(yōu)化是比較務(wù)實(shí)的選擇。原因也很直接你通常拿不到足夠多的歷史數(shù)據(jù)來精確擬合風(fēng)電和光伏誤差的分布函數(shù)就算擬合出來分布參數(shù)本身也有估計(jì)誤差。與其在概率上賭一個精度不如直接說我知道誤差不會超過某個范圍在這個范圍內(nèi)怎么折騰我的調(diào)度方案都兜得住。1.3 標(biāo)題里魯棒性和不確定性到底指什么標(biāo)題里的魯棒性英文對應(yīng)Robustness放在調(diào)度模型里指的就是當(dāng)實(shí)際的風(fēng)電、光伏出力偏離預(yù)測值時(shí)系統(tǒng)依然能通過儲能充放電調(diào)整、機(jī)組出力再調(diào)度、與主網(wǎng)功率交換等手段維持功率平衡并滿足所有安全約束。衡量魯棒性的標(biāo)準(zhǔn)不是平均工況下的成本而是最壞工況下的約束滿足能力。不確定性在這個問題里主要來自三個方向源端不確定性風(fēng)機(jī)出力和光伏出力受氣象條件影響預(yù)測模型再準(zhǔn)也有偏差荷端不確定性區(qū)域負(fù)荷受生產(chǎn)生活規(guī)律、天氣溫度等因素影響同樣存在預(yù)測誤差設(shè)備狀態(tài)不確定性比如某臺機(jī)組臨時(shí)故障、某條線路檢修但在調(diào)度主問題里通常把這類視為N-1預(yù)想故障單獨(dú)處理不會混進(jìn)源荷不確定性里。在工程復(fù)現(xiàn)時(shí)通常把源荷不確定性統(tǒng)一建模為一個不確定集合U集合里的每一個元素代表一組可能的源荷出力偏差向量。魯棒優(yōu)化要求對U里的任意一個偏差向量系統(tǒng)都能通過第二階段調(diào)整找到可行解。這就是后面整個模型的核心出發(fā)點(diǎn)。2. 多階段魯棒調(diào)度模型的數(shù)學(xué)框架與不確定性集合設(shè)計(jì)2.1 兩階段模型的前后配合日前決策與日內(nèi)調(diào)整標(biāo)題里寫的是多階段魯棒調(diào)度模型實(shí)際復(fù)現(xiàn)時(shí)最常見的落地形式是兩階段魯棒優(yōu)化模型也就是經(jīng)典的三層結(jié)構(gòu)min-max-min。之所以叫多階段是因?yàn)檎{(diào)度決策本身是分層分時(shí)段的第一階段日前決策階段在獲知次日風(fēng)光和負(fù)荷預(yù)測值之后提前決定機(jī)組的啟停狀態(tài)、儲能充放電基準(zhǔn)計(jì)劃、向上級電網(wǎng)購電的日前計(jì)劃。這些決策的特點(diǎn)是必須在看到實(shí)際偏差之前做出所以也叫here-and-now決策。第二階段日內(nèi)調(diào)整階段當(dāng)實(shí)際風(fēng)光出力、負(fù)荷與預(yù)測出現(xiàn)偏差后在日前計(jì)劃的基礎(chǔ)上通過調(diào)節(jié)儲能出力、機(jī)組出力微調(diào)、可削減負(fù)荷等手段消除功率不平衡。這個階段的決策是在偏差揭曉之后才能做的叫wait-and-see決策。兩階段之間通過功率平衡約束和設(shè)備運(yùn)行約束耦合在一起。模型的目標(biāo)函數(shù)也相應(yīng)地拆成兩塊min (第一階段運(yùn)行成本) max_{u∈U} min_{y∈Ω(x,u)} (第二階段調(diào)整成本)外層min決定第一階段決策x內(nèi)層max尋找使第二階段成本最高的最壞偏差場景u中間層min則是在給定x和u之后尋找最優(yōu)的調(diào)整策略y。整個模型求的是在所有可能出現(xiàn)在的偏差場景下系統(tǒng)總成本的最壞情況最優(yōu)值。2.2 盒式、預(yù)算約束與橢球式不確定集合的適用場景不確定集合怎么建直接決定了模型的保守程度和求解難度。常見的有三種形式盒式集合Box直接把每個時(shí)段的不確定量限制在一個區(qū)間內(nèi)如u(t)∈[u_min(t), u_max(t)]。這種集合最簡單但要求所有時(shí)段同時(shí)達(dá)到最大偏差結(jié)果是最壞中的最壞過于保守。預(yù)算約束集合Budget在盒式集合基礎(chǔ)上增加一個總偏差比例限制形如Σ|u(t)-u_hat(t)|/δ(t) ≤ Γ。這里的Γ就是預(yù)算參數(shù)表示在調(diào)度周期內(nèi)最多允許多少個時(shí)段的偏差同時(shí)達(dá)到上限。Γ越小魯棒約束越寬松成本越低Γ越大方案越保守。這個設(shè)置非常貼合實(shí)際因?yàn)檎鎸?shí)的風(fēng)光出力偏差不太可能24個小時(shí)每個小時(shí)都頂在最大值上。橢球式集合Ellipsoidal用2-范數(shù)約束替代1-范數(shù)額外限制形如‖W(u-u_hat)‖? ≤ Ω。橢球集合考慮了變量之間的相關(guān)性數(shù)學(xué)上更精細(xì)但會引入二次約束計(jì)算復(fù)雜度明顯上升。在微電網(wǎng)調(diào)度里我最推薦的是盒式預(yù)算約束的組合。原因有三個物理意義清晰Γ可以理解為最壞情況下同時(shí)偏離預(yù)測的時(shí)間段數(shù)方便向工程人員解釋模型保持線性結(jié)構(gòu)經(jīng)過對偶變換后仍然是LP或MILP可以用YALMIPGurobi穩(wěn)定求解通過調(diào)節(jié)Γ的取值可以畫出從激進(jìn)到保守的完整成本曲線這是論文和報(bào)告里很有說服力的一張圖。2.3 目標(biāo)函數(shù)與約束的完整數(shù)學(xué)表達(dá)以含風(fēng)電、光伏、儲能和微型燃?xì)廨啓C(jī)的區(qū)域微電網(wǎng)為例完整的魯棒調(diào)度模型可以寫成決策變量第一階段燃?xì)廨啓C(jī)開機(jī)狀態(tài)z(t)、儲能充電功率P_ess_c(t)、儲能放電功率P_ess_d(t)、日前購電功率P_grid(t)、儲能SOC(t)。第二階段針對偏差場景u的儲能調(diào)整量ΔP_ess(t)、機(jī)組調(diào)整量ΔP_g(t)、可削減負(fù)荷量ΔP_load(t)。目標(biāo)函數(shù)min Σ_t [ C_grid * P_grid(t) C_g * P_g(t) C_ess * (P_ess_c(t)P_ess_d(t)) ] max_{u∈U} min_Δ Σ_t [ C_penalty * (ΔP_load(t) 棄風(fēng)棄光懲罰項(xiàng)) ]約束條件功率平衡約束 P_grid(t) P_wind(t) P_pv(t) P_g(t) P_ess_d(t) P_load(t) P_ess_c(t) ΔP_load(t)其中P_wind(t)和P_pv(t)在第二階段表示為P_wind_hat(t) u_w(t)、P_pv_hat(t) u_pv(t)。儲能約束 SOC(t1) SOC(t) η_c * P_ess_c(t) * Δt - P_ess_d(t) * Δt / η_d SOC_min ≤ SOC(t) ≤ SOC_max 0 ≤ P_ess_c(t) ≤ P_ess_c_max * z_ess(t) 0 ≤ P_ess_d(t) ≤ P_ess_d_max * (1 - z_ess(t))這里需要注意的是充放電狀態(tài)變量z_ess(t)的引入它保證儲能不會同時(shí)充電和放電這在后面代碼實(shí)現(xiàn)時(shí)是個關(guān)鍵細(xì)節(jié)。機(jī)組約束 z(t) * P_g_min ≤ P_g(t) ≤ z(t) * P_g_max 爬坡約束|P_g(t) - P_g(t-1)| ≤ Ramp_g與主網(wǎng)交互約束 P_grid_min ≤ P_grid(t) ≤ P_grid_max可削減負(fù)荷約束 0 ≤ ΔP_load(t) ≤ α * P_load(t)2.4 為什么多階段不等于多周期很多初學(xué)者會把多階段魯棒理解成把24小時(shí)的調(diào)度周期分成很多段每段做一個魯棒優(yōu)化這是一個明顯的誤區(qū)。多周期優(yōu)化只是把時(shí)間維度拉長決策變量在整個周期內(nèi)一次性規(guī)劃本質(zhì)上仍是開環(huán)調(diào)度。多階段魯棒優(yōu)化的關(guān)鍵不在于時(shí)間段的多少而在于決策的可調(diào)整性——第一階段決策做出后可以看到不確定性的實(shí)現(xiàn)然后第二階段再做出調(diào)整決策。這個特征對應(yīng)的是閉環(huán)調(diào)度也就是日前計(jì)劃日內(nèi)實(shí)時(shí)調(diào)整的真實(shí)運(yùn)行流程。用數(shù)學(xué)語言說真正的多階段魯棒模型必須具有min-max-min的嵌套結(jié)構(gòu)而不是簡單的min-max或者包含很多時(shí)段的一次性min。在復(fù)現(xiàn)時(shí)如果你的模型只是每個時(shí)段都預(yù)留備用容量然后用魯棒約束包住所有時(shí)段這雖然也算魯棒但不是多階段魯棒。要體現(xiàn)多階段就必須顯式區(qū)分可調(diào)整變量如儲能的實(shí)時(shí)充放電調(diào)整量和不可調(diào)整變量如機(jī)組啟停、日前購電計(jì)劃并在目標(biāo)函數(shù)中體現(xiàn)第二階段對最壞場景的應(yīng)對成本。3. 列和約束生成CCG算法迭代框架的收斂邏輯3.1 主問題的最小費(fèi)用加最壞場景切割線兩階段魯棒模型直接求解極其困難因?yàn)閮?nèi)層嵌套了max-min優(yōu)化。工程上最常用的解法是列和約束生成算法核心思路是把原問題拆成一個主問題MP和一個子問題SP通過迭代逐步逼近最優(yōu)解。主問題的形式如下min c * x η s.t. 第一階段約束 η ≥ d * y_l, 對于已知的場景u_l x、y_l之間的耦合約束這里的關(guān)鍵是引入了輔助變量η它表示第二階段成本的上界估計(jì)。每輪迭代中子問題會找到一個最壞場景u_l主問題就把這個場景對應(yīng)的第二階段變量y_l和約束作為切割加入模型讓η不斷收緊。主問題因?yàn)橹豢紤]了已發(fā)現(xiàn)的最壞場景所以實(shí)際上是原問題的松弛版本它的最優(yōu)值給出原問題的下界LB。3.2 子問題的max-min對偶轉(zhuǎn)化子問題是在給定第一階段解x_star之后尋找使第二階段成本最大的場景u以及對應(yīng)的最優(yōu)調(diào)整成本形式是f(x_star) max_{u∈U} min_{y∈Ω(x_star, u)} d * y這個max-min結(jié)構(gòu)不能直接丟給求解器。標(biāo)準(zhǔn)做法是對內(nèi)層的min問題取強(qiáng)對偶。由于內(nèi)層是關(guān)于y的線性規(guī)劃只要可行域非空且有界強(qiáng)對偶成立就可以把min問題等價(jià)轉(zhuǎn)化為其對偶max問題然后和外層的max合并變成一個單層max問題f(x_star) max_{u∈U, λ∈Λ} λ * (b - A * x_star - C * u)其中λ是對偶變量它的符號由原約束方向決定。如果原約束是≤型λ非負(fù)如果是等式約束λ自由變量。合并之后子問題變成帶二次項(xiàng)的非線性問題因?yàn)棣撕蛈相乘。在實(shí)際求解中這個雙線性項(xiàng)通?？梢杂么驧法或直接交給支持MIQP的非線性求解器處理Gurobi可以直接求解這類帶雙線性項(xiàng)的小規(guī)模問題。需要強(qiáng)調(diào)一點(diǎn)內(nèi)層min問題必須滿足強(qiáng)對偶條件。因此第二階段調(diào)整問題必須建模為線性規(guī)劃不能隨意加入整數(shù)變量。如果第二階段需要決策是否切除某條線路這類0-1變量標(biāo)準(zhǔn)CCG就不適用了需要改成嵌套CCG或者引入其他分解方法。3.3 收斂判據(jù)與最大迭代次數(shù)CCG算法的主循環(huán)流程如下初始化UBinf, LB-inf, k0設(shè)置初始場景u_0通常取預(yù)測均值 while k ≤ K_max 且 UB - LB tol: 1. 求解主問題得到x_k*和η_k*更新LB c*x_k* η_k* 2. 將x_k*代入子問題求解最壞場景u_k*和最壞場景成本f(x_k*) 3. 更新UB min(UB, c*x_k* f(x_k*)) 4. if UB - LB ≤ tol: 收斂退出 5. 否則將u_k*作為新場景加入主問題添加對應(yīng)的第二階段變量y_k和割約束 k k 1收斂判據(jù)一般取相對間隙或絕對間隙。我在復(fù)現(xiàn)時(shí)常用的組合是絕對誤差取1e-4最大迭代次數(shù)取20輪。在24時(shí)段微電網(wǎng)模型里通常3到8輪就能收斂到很緊的間隙很少會真的跑到20輪。3.4 CCG與Benders分解的選擇建議很多文獻(xiàn)會同時(shí)提到Benders分解和CCG兩者都是求解兩階段問題和分解主-子問題的經(jīng)典算法但適用場景有明顯差異對比維度Benders分解CCG切割信息通過對偶函數(shù)添加割平面直接添加最壞場景對應(yīng)的完整約束迭代收斂對小規(guī)模問題收斂較快通常需要更少迭代次數(shù)尤其有連續(xù)變量時(shí)對整數(shù)變量的適應(yīng)性一般需要處理整數(shù)對偶gap主問題中可以自然處理整數(shù)變量實(shí)現(xiàn)難度較低中等在微電網(wǎng)魯棒調(diào)度里第一階段通常包含機(jī)組啟停這類整數(shù)變量而CCG的切割是把具體場景u_l下的約束完整塞進(jìn)主問題對整數(shù)變量天然友好所以實(shí)際復(fù)現(xiàn)時(shí)我基本都用CCG。只有當(dāng)子問題規(guī)模特別大、每次求解SP耗時(shí)很高時(shí)才會回頭考慮Benders或者近似方法。4. MATLABYALMIPGurobi的關(guān)鍵代碼實(shí)現(xiàn)4.1 環(huán)境準(zhǔn)備與求解器適配復(fù)現(xiàn)前先確認(rèn)環(huán)境MATLAB版本R2020b及以上即可建議R2022b以上對sdpvar和optimizer的兼容性更好YALMIP直接從GitHub拉取最新release不要用十年前的老版本求解器Gurobi 10.x或11.x學(xué)術(shù)許可免費(fèi)申請Mosek 10.x也可以但Gurobi在MILP上的速度通常更優(yōu)。在MATLAB里先跑一句驗(yàn)證yalmiptest能列出Gurobi的LP、QP、MILP、MIQP支持狀態(tài)就說明環(huán)境沒問題。如果用了舊版MATLAB配新版Gurobi經(jīng)常會出現(xiàn)java相關(guān)報(bào)錯那是MATLAB版本太老、Java版本不匹配導(dǎo)致的屬于環(huán)境問題先升級MATLAB再排查。4.2 主問題的YALMIP建模代碼下面是主問題的核心建模骨架采用24時(shí)段、1小時(shí)為步長的設(shè)定。nT 24; % 第一階段變量 z_g binvar(1, nT); % 燃?xì)廨啓C(jī)啟停 P_g sdpvar(1, nT); % 燃?xì)廨啓C(jī)出力 P_ess_c sdpvar(1, nT); % 儲能充電功率 P_ess_d sdpvar(1, nT); % 儲能放電功率 z_ess binvar(1, nT); % 儲能充放電狀態(tài)1充電0放電 SOC sdpvar(1, nT1); % 荷電狀態(tài) P_grid sdpvar(1, nT); % 購電功率 eta sdpvar(1); % 最壞場景成本上界 % 目標(biāo)函數(shù) objective_mp sum(grid_price .* P_grid) ... sum(g_fuel_cost * P_g) ... sum(ess_cost * (P_ess_c P_ess_d)) ... eta; % 約束 constraints_mp []; constraints_mp [constraints_mp, SOC(1) SOC_init]; for t 1:nT constraints_mp [constraints_mp, ... SOC(t1) SOC(t) eta_charge * P_ess_c(t) - P_ess_d(t) / eta_discharge, ... SOC_min SOC(t1) SOC_max, ... 0 P_ess_c(t) P_ess_c_max * z_ess(t), ... 0 P_ess_d(t) P_ess_d_max * (1 - z_ess(t)), ... z_g(t) * P_g_min P_g(t) z_g(t) * P_g_max, ... P_grid_min P_grid(t) P_grid_max, ... P_wind_forecast(t) P_pv_forecast(t) P_g(t) P_grid(t) P_ess_d(t) ... P_load_forecast(t) P_ess_c(t)]; end ops sdpsettings(solver, gurobi, verbose, 0, gurobi.NumericFocus, 3); optimize(constraints_mp, objective_mp, ops);這里有個細(xì)節(jié)給Gurobi設(shè)置NumericFocus為3可以明顯降低大M約束帶來的數(shù)值病態(tài)問題。如果模型里用到了大M這個設(shè)置幾乎是必加的。4.3 子問題的對偶化與最壞場景求解代碼第二階段的最壞場景求解是核心代碼分三步走第一步固定第一階段變量第二步構(gòu)造第二階段LP并取對偶第三步求解單層max問題。% 固定第一階段變量 x_star value([P_g; P_grid; P_ess_c; P_ess_d; SOC]); % 不確定變量 u_w sdpvar(1, nT); % 風(fēng)電偏差 u_pv sdpvar(1, nT); % 光伏偏差 u_l sdpvar(1, nT); % 負(fù)荷偏差 % 第二階段調(diào)整變量 delta_ess sdpvar(1, nT); % 儲能功率調(diào)整 delta_g sdpvar(1, nT); % 機(jī)組出力調(diào)整 delta_load sdpvar(1, nT); % 可削減負(fù)荷 % 不確定集合盒式預(yù)算 U_con [-u_w_max u_w u_w_max, ... -u_pv_max u_pv u_pv_max, ... -u_l_max u_l u_l_max, ... sum(abs(u_w) ./ u_w_max) sum(abs(u_pv) ./ u_pv_max) Gamma]; % 第二階段約束僅對應(yīng)給定x_star SP_con [...]; SP_obj sum(penalty_load * delta_load penalty_curtail * (u_w u_pv) ... C_delta_g * abs(delta_g) C_delta_ess * abs(delta_ess)); % 通過YALMIP取對偶或直接用求解器求解最大最小問題 % 更通用的做法是手動構(gòu)造對偶變量把max-min轉(zhuǎn)為單層max dual_lambda sdpvar(size(A_sp, 1), 1); dual_mu sdpvar(size(A_eq, 1), 1); % 對偶目標(biāo)函數(shù)注意線性化雙線性項(xiàng) dual_obj dual_lambda * (b_sp - A_sp * x_star_vec - B_u * [u_w; u_pv; u_l]) ... dual_mu * (b_eq - A_eq * x_star_vec); % 合并后的最壞場景問題 SP_merged_con [dual_lambda 0, dual_mu free, U_con]; result_sp optimize(SP_merged_con, -dual_obj, ops_sp);實(shí)際上當(dāng)雙線性項(xiàng)存在時(shí)Gurobi會以MIQP方式求解變量數(shù)不大時(shí)速度可以接受。如果模型規(guī)模較大建議用big-M法把雙線性項(xiàng)線性化把問題轉(zhuǎn)成MILP。4.4 迭代主循環(huán)與結(jié)果保存主循環(huán)代碼框架如下UB inf; LB -inf; k 0; tol 1e-4; Kmax 20; cut_results {}; while k Kmax (UB - LB) tol k k 1; % 1. 求解主問題 optimize(constraints_mp, objective_mp, ops); LB value(objective_mp); x_cur value([P_g; P_grid; P_ess_c; P_ess_d; SOC]); % 2. 求解子問題 optimize(SP_merged_con, -dual_obj, ops_sp); sp_cost value(dual_obj); u_cur value([u_w; u_pv; u_l]); day_ahead_cost value(sum(grid_price .* P_grid) ... sum(g_fuel_cost * P_g) ... sum(ess_cost * (P_ess_c P_ess_d))); UB min(UB, day_ahead_cost sp_cost); fprintf(迭代%d: LB%.4f, UB%.4f, gap%.6f\n, k, LB, UB, UB-LB); % 3. 添加新場景切割到主問題 y_new sdpvar(size_y, 1); constraints_mp [constraints_mp, ... y_new對應(yīng)的第二階段約束把u_cur作為常數(shù)代入, ... eta 第二階段成本表達(dá)式(y_new)]; cut_results{k}.u u_cur; cut_results{k}.cost sp_cost; end有幾個實(shí)現(xiàn)上的細(xì)節(jié)值得說不要把主問題的sdpvar變量每輪重新創(chuàng)建。變量結(jié)構(gòu)要保持不變只是追加約束。否則MATLAB每輪都要重新解析整個符號表達(dá)式速度直接掉一個數(shù)量級。主問題的optimize結(jié)果如果返回problemsolver 信息先打印optimize的第二個輸出看看是primal infeasible還是numerical issues再決定調(diào)參數(shù)還是調(diào)模型。保存每輪的最壞場景u_cur最后畫圖時(shí)可以直接把最壞場景下的風(fēng)電光伏曲線和基準(zhǔn)預(yù)測曲線對比很直觀。4.5 代碼性能優(yōu)化的幾個細(xì)節(jié)復(fù)現(xiàn)完跑第一版10輪迭代可能要跑幾分鐘如果覺得慢按下面這幾個方向優(yōu)化通常能把時(shí)間壓到原來的三分之一用矩陣約束代替for循環(huán)。上面為了可讀性用了for循環(huán)實(shí)際寫代碼時(shí)盡量把所有時(shí)段的約束拼接成大矩陣一次性傳給constraints。YALMIP在拼接大規(guī)模符號矩陣時(shí)效率一般減少循環(huán)能明顯提升解析速度。子問題求解時(shí)不要用默認(rèn)的optimize改用optimizer對象把x_star作為參數(shù)傳入省去每輪重建約束的開銷。子問題的對偶變量數(shù)量通常和約束數(shù)量一致如果第二階段約束很多先做約束化簡把明顯冗余的約束刪掉。Gurobi的線程數(shù)設(shè)為4到8ops.gurobi.Threads默認(rèn)是0自動在共享計(jì)算環(huán)境里可以顯式設(shè)置。模型里大量使用abs()時(shí)YALMIP會引入輔助變量和不等式變量規(guī)模膨脹很快盡量在建模階段就用正負(fù)分解來避開abs。5. 完全復(fù)現(xiàn)中的參數(shù)設(shè)定、結(jié)果校核與常見坑5.1 推薦參考輸入?yún)?shù)參數(shù)設(shè)定直接影響收斂性和結(jié)果合理性下面這組參數(shù)是我在多個復(fù)現(xiàn)案例里校準(zhǔn)過的參考值可以直接作為起點(diǎn)參數(shù)推薦值說明調(diào)度周期24時(shí)段1h/時(shí)段多階段魯棒最基本的設(shè)定風(fēng)電裝機(jī)500 kW預(yù)測誤差取±15%光伏裝機(jī)300 kW預(yù)測誤差取±10%負(fù)荷峰值600 kW預(yù)測誤差取±5%儲能容量400 kWhSOC范圍20%~95%儲能最大充放電功率100 kW充放電效率均取0.95燃?xì)廨啓C(jī)容量200 kW最小出力50 kW爬坡率50kW/h購電價(jià)格上限1.2 元/kWh分時(shí)電價(jià)峰谷差越大越能體現(xiàn)儲能價(jià)值預(yù)算參數(shù)Γ3~5表示最多3~5個時(shí)段同時(shí)達(dá)到最大偏差收斂閾值1e-4絕對間隙最大迭代次數(shù)20實(shí)際通常8輪內(nèi)收斂5.2 怎么判斷你復(fù)現(xiàn)的模型確實(shí)魯棒跑通模型只是第一步關(guān)鍵是要會驗(yàn)證。我復(fù)現(xiàn)完之后一般做三組校驗(yàn)任何一組不過都說明模型有問題第一組最壞場景重放。把CCG最后一輪得到的最壞場景u_cur代入第二階段模型在固定第一階段決策的情況下重新求解第二階段調(diào)整問題。如果找不到可行解說明子問題建模或者切割添加有問題。第二組蒙特卡洛模擬。從風(fēng)電、光伏、負(fù)荷的預(yù)測誤差分布中隨機(jī)抽樣生成1000個場景對每個場景調(diào)用第二階段模型檢驗(yàn)是否可行。如果可行率低于98%說明不確定性集合設(shè)計(jì)偏激進(jìn)需要調(diào)大?；驍U(kuò)大盒式范圍。第三組Γ對比分析。分別取Γ0、1、3、6、12觀察總成本和最壞場景成本的變化。Γ0時(shí)模型退化為確定性問題總成本最低隨著Γ增大總成本單調(diào)上升這是符合預(yù)期的。如果曲線出現(xiàn)非單調(diào)或者劇烈跳變大概率是不確定集合的約束寫錯了。5.3 我實(shí)際踩過的幾個大坑第一個坑是儲能的SOC約束沒有處理好充放電同時(shí)性。模型里P_ess_c和P_ess_d是兩個獨(dú)立變量如果不加z_ess互斥約束求解器會出現(xiàn)一個時(shí)段既充電又放電的情況利用效率空轉(zhuǎn)降低成本看起來成本很漂亮實(shí)際上物理不可行。解決辦法就是代碼里用z_ess做互斥或者加入一個很小的同時(shí)充放電懲罰項(xiàng)。第二個坑是對偶變量符號搞錯。子問題在對偶化時(shí)如果原約束是Ax Bu ≤ b對偶變量必須是非負(fù)的如果是等式約束對偶變量自由。符號一錯對偶目標(biāo)的值會偏離真實(shí)最壞場景成本導(dǎo)致UB和LB無法收斂甚至發(fā)散。第三個坑是無腦擴(kuò)大M。在線性化雙線性項(xiàng)時(shí)M取值太小會剪掉有效解M取值太大會導(dǎo)致數(shù)值病態(tài)。解決方法是盡量用不確定變量的區(qū)間上下界來推導(dǎo)M的理論值而不是隨便取一個1e6。Gurobi的NumericFocus開啟后有一定緩解但不解決根本問題。第四個坑是YALMIP的符號展開。模型復(fù)雜時(shí)YALMIP默認(rèn)會把所有約束展開成完整系數(shù)矩陣變量多的時(shí)候內(nèi)存占用會爆炸。設(shè)置sdpsettings(expand, 0)可以減輕這個問題但代價(jià)是求解器內(nèi)部的presolve工作量增加需要實(shí)測權(quán)衡。5.4 復(fù)現(xiàn)后如何擴(kuò)展這個模型跑通基礎(chǔ)版之后有幾個很自然的擴(kuò)展方向都能在不推翻框架的前提下增加價(jià)值分布魯棒優(yōu)化把固定盒式集合換成Wasserstein球利用歷史數(shù)據(jù)構(gòu)建經(jīng)驗(yàn)分布再以球內(nèi)最壞分布為目標(biāo)。這個擴(kuò)展只需要改不確定集合的建模和對偶約束CCG框架可以復(fù)用。多能互補(bǔ)在電儲能的基礎(chǔ)上加入蓄熱罐、燃?xì)忮仩t把區(qū)域微電網(wǎng)擴(kuò)展成綜合能源系統(tǒng)。多能耦合會引入熱功率平衡約束但兩階段魯棒的結(jié)構(gòu)完全一致。碳約束在目標(biāo)函數(shù)中加入碳排放懲罰或碳配額約束分析碳交易價(jià)格對魯棒調(diào)度策略的影響。這類擴(kuò)展對接雙碳主題在實(shí)際匯報(bào)中更有說服力。儲能容量優(yōu)化把儲能的額定容量和功率作為第一階段決策變量目標(biāo)函數(shù)中加入儲能投資年化成本模型就變成了容量規(guī)劃魯棒調(diào)度聯(lián)合優(yōu)化問題CCG依然適用。我個人在多次復(fù)現(xiàn)這個模型后的體會是兩階段魯棒框架本身并不算特別難難的是把每個環(huán)節(jié)的物理意義和數(shù)學(xué)細(xì)節(jié)對齊。比如Γ怎么解釋給非優(yōu)化背景的同事聽比如儲能SOC方程里為什么充電效率在乘法位置、放電效率在除法位置這些細(xì)節(jié)如果不理清楚代碼復(fù)制過來改了參數(shù)結(jié)果不對勁都不知道去哪查。建議第一次復(fù)現(xiàn)的朋友先用小系統(tǒng)比如6節(jié)點(diǎn)、4時(shí)段把CCG迭代曲線畫出來確認(rèn)LB和UB的收斂趨勢沒問題了再放大到24時(shí)段完整系統(tǒng)。小系統(tǒng)里能一眼看出的邏輯錯誤放大以后往往會被數(shù)值噪音掩蓋排查難度翻倍。這個模型跑通之后你再去讀其他魯棒調(diào)度相關(guān)的論文思路會清晰很多。