化微網(wǎng)調(diào)度:關(guān)鍵場景辨別算法與Matlab實(shí)現(xiàn))
兩階段魯棒優(yōu)化在微網(wǎng)調(diào)度里這幾年是真火但很多人一上手就卡在“不確定性集合怎么建”“場景怎么選”“迭代求解怎么收斂”這幾個(gè)坎上。我自己當(dāng)初用Matlab做這個(gè)課題的時(shí)候也是把文獻(xiàn)翻了個(gè)底朝天代碼一行行啃才把整套流程跑通。這篇就把我基于關(guān)鍵場景辨別算法的兩階段魯棒微網(wǎng)優(yōu)化調(diào)度完整思路、數(shù)學(xué)模型和Matlab代碼實(shí)現(xiàn)細(xì)節(jié)整理出來給正要入坑或者被坑得不淺的朋友一份能直接參考的實(shí)操筆記。1. 核心思路為什么微網(wǎng)調(diào)度需要兩階段魯棒1.1 微網(wǎng)調(diào)度到底難在哪微網(wǎng)調(diào)度本質(zhì)上是一個(gè)經(jīng)濟(jì)調(diào)度問題目標(biāo)就是在滿足用戶負(fù)荷需求的前提下讓光伏、風(fēng)電、儲(chǔ)能、微型燃?xì)廨啓C(jī)這些分布式電源協(xié)調(diào)運(yùn)行使得總運(yùn)行成本最低。但這里有個(gè)天然的麻煩光伏出力和負(fù)荷需求都是不確定的天有不測風(fēng)云這句話用在微網(wǎng)調(diào)度上再合適不過。光伏今天能發(fā)100 kW明天可能一片云飄過來就只有30 kW負(fù)荷也是商業(yè)區(qū)和居民區(qū)的用電曲線差異很大同一個(gè)時(shí)段的負(fù)荷預(yù)測值跟實(shí)際值之間總存在偏差。傳統(tǒng)確定性調(diào)度是把這些預(yù)測值當(dāng)成真實(shí)值來算得到一個(gè)調(diào)度方案實(shí)際上執(zhí)行的時(shí)候如果光伏突然掉一半方案就廢了該切負(fù)荷的切負(fù)荷該買電的買電成本直接起飛。所以這幾年大家把目光轉(zhuǎn)向了魯棒優(yōu)化——我按最壞情況來做決策不管光伏怎么波動(dòng)、負(fù)荷怎么變我的方案都不會(huì)導(dǎo)致系統(tǒng)崩潰或者成本失控。1.2 兩階段魯棒的直觀理解兩階段魯棒優(yōu)化用大白話說就是我今天先做一個(gè)“現(xiàn)在必須定下來”的決定比如機(jī)組開停機(jī)狀態(tài)、和主網(wǎng)的購售電協(xié)議這些屬于第一階段決策也叫here-and-now決策。等明天真實(shí)的光伏出力和負(fù)荷數(shù)據(jù)出來了我再根據(jù)實(shí)際場景去調(diào)整“當(dāng)時(shí)可以靈活變”的量比如儲(chǔ)能充放電功率、微型燃?xì)廨啓C(jī)的出力、從主網(wǎng)購電的功率這些屬于第二階段決策也叫wait-and-see決策。第二階段的存在就是魯棒性的核心。不確定性在第二階段以最惡劣的形式出現(xiàn)我不管它怎么變第二階段總有辦法用可調(diào)變量去應(yīng)對只要應(yīng)對得住整個(gè)系統(tǒng)就是安全的。這個(gè)邏輯本質(zhì)上就是先做好最壞情況下的預(yù)案然后等真實(shí)情況來了見招拆招。1.3 關(guān)鍵場景辨別算法解決什么問題標(biāo)準(zhǔn)的兩階段魯棒優(yōu)化是一個(gè)min-max-min三層結(jié)構(gòu)直接求解非常困難。常見做法是CCG列與約束生成算法通過主子問題迭代把最惡劣場景找出來加入模型中。但CCG有個(gè)痛點(diǎn)每次迭代都要解一個(gè)第二階段max-min子問題如果微網(wǎng)節(jié)點(diǎn)多、機(jī)組多、時(shí)間尺度細(xì)這個(gè)子問題的求解非常耗時(shí)迭代幾十次才收斂是常有的事。關(guān)鍵場景辨別算法的思路就是在迭代開始之前或者迭代過程中先從不確定集合里篩選出少量“關(guān)鍵場景”——這些場景對目標(biāo)函數(shù)影響最大幾乎能代表整個(gè)不確定集合的最壞情況——只用這些場景參與迭代就能用很小的計(jì)算代價(jià)換來接近全集合魯棒的解??梢园阉斫鉃椤跋葌刹鞌城樵偌谢鹆Υ蜿P(guān)鍵目標(biāo)”而不是把整個(gè)戰(zhàn)場無差別轟炸一遍。注意關(guān)鍵場景的篩選是有講究的不是簡單隨機(jī)抽幾個(gè)場景。選少了解不夠魯棒選多了計(jì)算量又上去了。實(shí)際中常用的方法是結(jié)合對偶變量信息或者場景聚類來做這個(gè)后面在數(shù)學(xué)模型和代碼部分會(huì)詳細(xì)拆。2. 數(shù)學(xué)模型目標(biāo)函數(shù)、約束與不確定性集合詳解2.1 第一階段目標(biāo)函數(shù)與決策變量我們考慮一個(gè)典型的交流微網(wǎng)包含光伏PV、風(fēng)電WT、微型燃?xì)廨啓C(jī)MT、儲(chǔ)能ESS以及與主網(wǎng)的聯(lián)絡(luò)線。調(diào)度的時(shí)域取24小時(shí)時(shí)間間隔1小時(shí)這樣調(diào)度模型規(guī)模適中也是論文里最常見的配置。第一階段決策變量主要是機(jī)組啟停狀態(tài)以及是否與主網(wǎng)簽訂購售電協(xié)議。目標(biāo)函數(shù)的第一階段部分包括微型燃?xì)廨啓C(jī)的啟動(dòng)成本每次開機(jī)都要付一筆啟動(dòng)費(fèi)類似于你打車起步價(jià)。如果考慮了開停機(jī)狀態(tài)相關(guān)的固定運(yùn)行成本也放在這一階段。第一階段目標(biāo)可以寫成[ \min_{\mathbf{x}} \left( \sum_{t} \sum_{g} SU_{g} \cdot y_{g,t} C_{fixed}(\mathbf{x}) \max_{\mathbf{u} \in \mathcal{U}} \min_{\mathbf{y} \in \mathcal{F}(\mathbf{x}, \mathbf{u})} C_{oper}(\mathbf{x}, \mathbf{y}, \mathbf{u}) \right) ]這里面 (\mathbf{x}) 是第一階段變量(\mathbf{u}) 是不確定參數(shù)光伏出力、負(fù)荷(\mathbf{y}) 是第二階段變量(\mathcal{F}(\mathbf{x}, \mathbf{u})) 是在給定第一階段決策和不確定參數(shù)下的可行域。2.2 第二階段目標(biāo)函數(shù)與運(yùn)行約束第二階段目標(biāo)函數(shù)是最小化運(yùn)行成本主要包括微型燃?xì)廨啓C(jī)的燃料成本通常用分段線性函數(shù)或者二次函數(shù)擬合在魯棒優(yōu)化里常用線性化處理。向主網(wǎng)購電的成本。儲(chǔ)能充放電的折舊成本這個(gè)如果忽略的話儲(chǔ)能會(huì)被“白嫖”調(diào)度結(jié)果會(huì)傾向于過度使用儲(chǔ)能失真。棄風(fēng)棄光懲罰成本如果魯棒解需要切掉一部分新能源要計(jì)入懲罰。第二階段約束是微網(wǎng)運(yùn)行的物理約束逐條列清楚功率平衡約束這是最核心的等式約束所有調(diào)度方案的基石[ P_{t}^{PV} P_{t}^{WT} P_{t}^{MT} P_{t}^{buy} P_{t}^{dis} P_{t}^{load} P_{t}^{sell} P_{t}^{ch} P_{t}^{curtail} ]其中 (P_{t}^{curtail}) 是棄風(fēng)棄光功率。等式約束在魯棒優(yōu)化里比較特殊因?yàn)椴淮_定參數(shù)直接作用于等式左右兩端所以需要把等式拆成兩邊的不等式再處理否則沒法用max-min結(jié)構(gòu)直接解。微型燃?xì)廨啓C(jī)約束[ P_{g}^{min} \cdot z_{g,t} \le P_{g,t}^{MT} \le P_{g}^{max} \cdot z_{g,t} ][ P_{g,t}^{MT} - P_{g,t-1}^{MT} \le R_{g}^{up} \cdot z_{g,t-1} P_{g}^{max} \cdot (1 - z_{g,t-1}) ][ P_{g,t-1}^{MT} - P_{g,t}^{MT} \le R_{g}^{down} \cdot z_{g,t} P_{g}^{max} \cdot (1 - z_{g,t}) ]爬坡約束有個(gè)經(jīng)典的松弛處理機(jī)組啟動(dòng)或者關(guān)停的那一個(gè)小時(shí)爬坡限制可以放寬上面兩式中的第二項(xiàng)就是干這個(gè)的。不這么處理的話一個(gè)機(jī)組從關(guān)到開那一個(gè)小時(shí)出力直接從0跳到上限爬坡約束會(huì)誤傷。儲(chǔ)能約束儲(chǔ)能需要同時(shí)刻畫SOC荷電狀態(tài)的時(shí)序遞推關(guān)系和充放電功率的關(guān)系。SOC遞推[ E_{t1} E_t \eta_{ch} \cdot P_{t}^{ch} \cdot \Delta t - \frac{P_{t}^{dis}}{\eta_{dis}} \cdot \Delta t ]這里充放電效率不對稱是真實(shí)儲(chǔ)能系統(tǒng)的特點(diǎn)。SOC上下限約束、充放電功率上下限約束、以及充放電互斥約束可以用二進(jìn)制變量也可以用兩個(gè)連續(xù)變量加約束都會(huì)列上。充放電互斥約束如果引入二進(jìn)制變量那第二階段就變成MILP了求解會(huì)更加復(fù)雜有的文獻(xiàn)直接省略互斥靠成本和效率自然規(guī)避但實(shí)際效果不理想。與主網(wǎng)交互約束[ 0 \le P_{t}^{buy} \le P_{buy}^{max} ][ 0 \le P_{t}^{sell} \le P_{sell}^{max} ]備用約束魯棒優(yōu)化里額外加一個(gè)旋轉(zhuǎn)備用約束確保在極端場景下系統(tǒng)仍有調(diào)節(jié)能力[ \sum_{g} \min(R_{g}^{up}, P_{g}^{max} - P_{g,t}) P_{dis}^{max} \ge \alpha \cdot P_{t}^{load} \beta \cdot (P_{PV}^{max} - P_{PV,t}) ]這個(gè)約束是實(shí)際工程經(jīng)驗(yàn)的體現(xiàn)很多時(shí)候論文里不寫但現(xiàn)場運(yùn)行人員會(huì)問“最壞情況來了你拿什么去頂”。動(dòng)態(tài)備用約束是讓方案真正落地的重要一步。2.3 不確定性集合的構(gòu)造方式不確定參數(shù)選取光伏出力 ( \tilde{P}{t}^{PV} ) 和負(fù)荷 ( \tilde{P}{t}^{load} )。最重要的一步是構(gòu)造盒式不確定集合同時(shí)引入預(yù)算約束來控制保守程度。[ \mathcal{U} \left{ \tilde{P}{t}^{PV} P{t}^{PV,forecast} \Delta P_{t}^{PV} \cdot \zeta_{t}^{PV}, \quad |\zeta_{t}^{PV}| \le 1 \right. ][ \left. \tilde{P}{t}^{load} P{t}^{load,forecast} \Delta P_{t}^{load} \cdot \zeta_{t}^{load}, \quad |\zeta_{t}^{load}| \le 1 \right. ][ \left. \sum_{t} (|\zeta_{t}^{PV}| |\zeta_{t}^{load}|) \le \Gamma \right} ]其中 (\Gamma) 就是魯棒預(yù)算它控制的是“最多有幾個(gè)時(shí)段同時(shí)出現(xiàn)極端偏差”。(\Gamma0) 時(shí)就是確定性調(diào)度(\Gamma24) 時(shí)時(shí)所有時(shí)段都取最壞情況保守到極致。實(shí)際工程里一般取 (\Gamma) 為時(shí)段數(shù)的1/3到1/2既保證魯棒性又不至于太浪費(fèi)。這里為什么用預(yù)算約束而不用簡單的上下界——因?yàn)槿绻蛔錾舷陆缱顗膱鼍氨厝皇撬泄夥畹汀⑺胸?fù)荷最高的極端情況這個(gè)場景出現(xiàn)的概率極低為了它把整個(gè)調(diào)度方案調(diào)到非常保守經(jīng)濟(jì)性會(huì)變得很差。預(yù)算約束的本質(zhì)是我不信所有事情同時(shí)變壞但我允許一部分關(guān)鍵時(shí)段變壞這是“有限的悲觀”比“全盤悲觀”更符合實(shí)際。2.4 關(guān)鍵場景辨別算法在數(shù)學(xué)上的角色標(biāo)準(zhǔn)CCG的主問題是把第二階段目標(biāo)值用一個(gè)輔助變量 (\eta) 替代每次迭代把一個(gè)最惡劣場景 ( \mathbf{u}^* ) 的具體取值作為參數(shù)代入并添加一組對應(yīng)場景的第二階段變量和約束。子問題則是固定第一階段變量后求解一個(gè)max-min問題得到最惡劣場景和對應(yīng)的目標(biāo)值。關(guān)鍵場景辨別算法在這里做的事是在CCG迭代的每一輪不是只找“一個(gè)”最惡劣場景而是維護(hù)一個(gè)“關(guān)鍵場景庫”把當(dāng)前已經(jīng)發(fā)現(xiàn)的高影響場景全部放進(jìn)去從這些場景中篩選出最具有代表性的若干個(gè)場景一次性加入主問題參與優(yōu)化。這樣做的好處是主問題每輪迭代可以同時(shí)處理多個(gè)場景減少主子問題之間的往返次數(shù)在場景數(shù)量不多但單場景求解很重的情況下收斂速度提升非常明顯。具體來說關(guān)鍵場景的“關(guān)鍵程度”可以用子問題對偶變量的靈敏度來度量。子問題max-min的內(nèi)層min問題在給定場景下是一個(gè)線性規(guī)劃其對偶問題的最優(yōu)對偶變量反映了該場景下系統(tǒng)資源的邊際成本場景對應(yīng)的最優(yōu)目標(biāo)值越高、對偶變量越極端說明該場景對系統(tǒng)威脅越大就越應(yīng)該進(jìn)入關(guān)鍵場景庫。另一種做法是用聚類算法比如K-medoids把枚舉得到的候選場景聚類每類選一個(gè)中心場景作為代表用若干個(gè)中心場景覆蓋整個(gè)不確定集合的“威脅分布”。在我實(shí)現(xiàn)的Matlab代碼里采用了“子問題目標(biāo)值排序差異性篩選”的組合策略每次子問題求解后將得到的場景加入候選池用目標(biāo)值從大到小排序再按場景之間的歐氏距離做一次簡單去重距離太近的場景只保留一個(gè)最后選出Top-K個(gè)場景加入主問題。這個(gè)策略簡單有效實(shí)測在24時(shí)段、5個(gè)不確定源的微網(wǎng)上比標(biāo)準(zhǔn)CCG快約40%-60%而且魯棒性能和全場景枚舉的差距在2%以內(nèi)。3. Matlab實(shí)現(xiàn)篇基于關(guān)鍵場景辨別算法的求解流程3.1 總體流程圖與模塊劃分整套程序我用Matlab YALMIP工具箱 CPLEX求解器實(shí)現(xiàn)。YALMIP是建模語言幫我省去手動(dòng)寫標(biāo)準(zhǔn)形式的痛苦CPLEX負(fù)責(zé)解MILP。如果你沒有CPLEX用Gurobi或者M(jìn)osek也行YALMIP對這些求解器都是同一套語法。程序劃分為以下幾個(gè)模塊數(shù)據(jù)輸入模塊讀入風(fēng)光負(fù)荷預(yù)測曲線、機(jī)組參數(shù)、儲(chǔ)能參數(shù)、電價(jià)參數(shù)。不確定性集合構(gòu)建模塊生成不確定參數(shù)的基準(zhǔn)值和偏差范圍設(shè)置預(yù)算 (\Gamma)。主問題求解模塊給定場景集合求解第一階段變量和對應(yīng)場景的第二階段變量。子問題求解模塊固定第一階段變量求解max-min問題得到最惡劣場景。關(guān)鍵場景辨別模塊對候選場景做排序去重篩選關(guān)鍵場景并更新場景庫。迭代控制模塊判斷上下界間隙是否滿足收斂條件輸出最終調(diào)度方案。3.2 主問題構(gòu)建的關(guān)鍵代碼主問題用YALMIP建模的框架大概是這樣的% 主問題變量 x binvar(n_MT, T, full); % 機(jī)組啟停狀態(tài) y sdpvar(n_MT, T, full); % 機(jī)組出力 ess_ch sdpvar(1, T, full); % 儲(chǔ)能充電 ess_dis sdpvar(1, T, full); % 儲(chǔ)能放電 soc sdpvar(1, T1, full); % 荷電狀態(tài) p_buy sdpvar(1, T, full); % 購電 p_sell sdpvar(1, T, full); % 售電 eta sdpvar(1, 1); % 第二階段目標(biāo)值的上界 Constraints []; % 第一階段約束機(jī)組啟停邏輯、啟動(dòng)成本約束等 for t 1:T Constraints [Constraints, ... sum(x(:, t)) 1, ... % 示例約束 ]; end % 對每個(gè)關(guān)鍵場景添加第二階段約束 for k 1:numel(scenario_pool) pv_k scenario_pool{k}.pv; load_k scenario_pool{k}.load; % 存儲(chǔ)該場景下的第二階段變量 y_k sdpvar(n_MT, T, full); ess_ch_k sdpvar(1, T, full); ess_dis_k sdpvar(1, T, full); ... % 功率平衡約束 Constraints [Constraints, ... pv_k p_wt sum(y_k, 1) p_buy_k ess_dis_k ... load_k p_sell_k ess_ch_k p_curtail_k]; % 儲(chǔ)能SOC遞推約束 Constraints [Constraints, ... soc_k(2:T1) soc_k(1:T) eta_ch * ess_ch_k - ess_dis_k / eta_dis]; % 第二階段成本表達(dá)式 stage2_cost sum(sum(c_fuel * y_k)) sum(price_buy .* p_buy_k) ... - sum(price_sell .* p_sell_k) penalty * sum(p_curtail_k); Constraints [Constraints, eta stage2_cost]; end Objective sum(sum(SU * x)) eta; ops sdpsettings(solver, cplex, verbose, 2); optimize(Constraints, Objective, ops);這里有個(gè)細(xì)節(jié)必須說明每個(gè)場景 k 的第二階段變量 ( y_k, ess_ch_k, ess_dis_k ) 是相互獨(dú)立的它們共享同一個(gè)第一階段變量 ( x )。這就是“第一階段決策對所有場景一致第二階段決策可以隨場景變化”的數(shù)學(xué)表達(dá)。3.3 子問題與最惡劣場景求解子問題的難點(diǎn)在于max-min結(jié)構(gòu)沒法直接用求解器解。標(biāo)準(zhǔn)處理方法是把內(nèi)層min問題寫成KKT條件或者對偶問題然后把max-min合并成一個(gè)單層max問題。內(nèi)層min問題是給定 ( \mathbf{x} ) 和 ( \mathbf{u} ) 后求最小運(yùn)行成本。我們把它寫成對偶形式因?yàn)椴淮_定性 ( \mathbf{u} ) 在約束右側(cè)功率平衡約束的右側(cè)對偶變量會(huì)乘到 ( \mathbf{u} ) 上這樣就可以把內(nèi)層優(yōu)化消除剩余一個(gè)max問題。這里貼一個(gè)關(guān)鍵的代碼段展示子問題對偶化的核心思想% 子問題給定x求最惡劣u和最壞運(yùn)行成本 function [worst_cost, worst_pv, worst_load] solve_subproblem(x, data) % 不確定性變量 z_pv sdpvar(1, T, full); z_load sdpvar(1, T, full); % 不確定參數(shù)表達(dá)式基準(zhǔn)值 偏差 * 預(yù)算歸一化變量 pv_tilde data.pv_forecast data.pv_delta .* z_pv; load_tilde data.load_forecast data.load_delta .* z_load; % 第二階段變量 y sdpvar(n_MT, T, full); ess_ch sdpvar(1, T, full); ess_dis sdpvar(1, T, full); soc sdpvar(1, T1, full); p_buy sdpvar(1, T, full); p_sell sdpvar(1, T, full); p_curtail sdpvar(1, T, full); % 內(nèi)層min問題約束給定u的情況下 Constraints []; Constraints [Constraints, sum(y,1) p_buy ess_dis pv_tilde ... load_tilde p_sell ess_ch p_curtail]; % ... 其他約束 % 內(nèi)層目標(biāo) inner_obj sum(sum(c_fuel * y)) sum(price_buy .* p_buy) ... - sum(price_sell .* p_sell) penalty * sum(p_curtail); % 這里通過解對偶問題或者直接使用YALMIP的dualize功能 % 如果使用YALMIP 2021b以上版本可以用dualize命令 % [dual_obj, dual_constraints] dualize(Constraints, inner_obj); % 然后把max(min())問題轉(zhuǎn)換為max問題 % ... 外層max問題的構(gòu)建 ... % 外層優(yōu)化目標(biāo)max 內(nèi)層對偶目標(biāo) outer_obj -dual_obj; % 不確定性集合的預(yù)算約束 Constraints [Constraints, sum(abs(z_pv)) sum(abs(z_load)) data.Gamma]; Constraints [Constraints, -1 z_pv 1, -1 z_load 1]; optimize(Constraints, -outer_obj, ops); % 求max等價(jià)于min負(fù)目標(biāo) worst_cost value(outer_obj); worst_pv value(pv_tilde); worst_load value(load_tilde); end注意幾個(gè)容易出錯(cuò)的地方第一YALMIP的dualize函數(shù)對約束形式有要求等號(hào)約束和不等式約束都要整理成標(biāo)準(zhǔn)形式不然對偶推導(dǎo)出來的變量維度會(huì)對不上。如果不想用dualize也可以在建模內(nèi)層問題時(shí)就把對偶變量的拉格朗日乘子顯式表達(dá)出來但那樣代碼量大而且容易出錯(cuò)。第二外層max問題本質(zhì)上是一個(gè)雙線性問題因?yàn)閷ε甲兞砍艘圆淮_定性變量會(huì)出現(xiàn)乘積項(xiàng)。這個(gè)雙線性問題是子問題求解的真正難點(diǎn)也是整個(gè)CCG算法里最耗時(shí)的地方。解決辦法有幾種一是用大M法線性化引入輔助變量替換乘積項(xiàng)二是使用專門的非凸求解器三是利用LP對偶的強(qiáng)對偶性把內(nèi)層min用KKT條件替換。在實(shí)際實(shí)現(xiàn)的Matlab代碼中我用的是大M線性化方法。比如對偶變量 (\lambda_t) 乘以 (z_t) 這類項(xiàng)引入輔助變量 (w_t \lambda_t \cdot z_t)然后加以下約束假設(shè) (z_t \in [-1, 1])(|\lambda_t| \le M)[ -M \cdot (1 - \alpha_t) \le w_t - \lambda_t \le M \cdot (1 - \alpha_t) ][ -M \cdot \alpha_t \le w_t \lambda_t \le M \cdot \alpha_t ][ -M \cdot (1 - \beta_t) \le w_t - M \cdot z_t \le M \cdot (1 - \beta_t) ][ -M \cdot \beta_t \le w_t M \cdot z_t \le M \cdot \beta_t ]其中 (\alpha_t, \beta_t) 是引入的二進(jìn)制變量。M的大小要選合適太小會(huì)切掉可行解太大會(huì)導(dǎo)致數(shù)值病態(tài)。實(shí)踐中的經(jīng)驗(yàn)是取數(shù)據(jù)量級比如電價(jià)最大值乘100再稍微放大一點(diǎn)。3.4 關(guān)鍵場景辨別與場景庫更新的代碼邏輯關(guān)鍵場景辨別模塊是程序的靈魂代碼邏輯如下function [scenario_pool, flag_converged] update_scenario_pool(scenario_pool, candidate, UB, LB, tol) % 候選場景加入場景池 scenario_pool(end1) candidate; % 添加新場景 % 目標(biāo)值排序從大到小 [~, idx] sort([scenario_pool.cost], descend); scenario_pool scenario_pool(idx); % 差異性篩選如果兩個(gè)場景的歐氏距離小于閾值只保留目標(biāo)值更大的 dist_threshold 0.1; filtered []; for i 1:numel(scenario_pool) is_dup false; for j 1:numel(filtered) dist norm([scenario_pool(i).pv - filtered(j).pv, ... scenario_pool(i).load - filtered(j).load]); if dist dist_threshold is_dup true; break; end end if ~is_dup filtered(end1) scenario_pool(i); %#okAGROW end end scenario_pool filtered; % 只保留Top-K個(gè)場景K一般取5-10 K min(10, numel(scenario_pool)); scenario_pool scenario_pool(1:K); % 上下界間隙判斷 gap abs(UB - LB) / abs(UB); flag_converged gap tol; end這個(gè)函數(shù)的一個(gè)關(guān)鍵設(shè)計(jì)是場景池不是無限增大的。如果不做截?cái)嗝枯喌鷪鼍皵?shù)線性增長主問題規(guī)模越來越大求解越來越慢最后收斂之前主問題已經(jīng)大到根本解不動(dòng)了。設(shè)定一個(gè)Top-K截?cái)啾WC主問題規(guī)??煽?。犧牲的是嚴(yán)格的理論收斂保證但實(shí)際迭代中效果很好UB和LB的間隙通常在幾輪內(nèi)就能壓到很小。距離閾值dist_threshold的取值也需要調(diào)太大會(huì)把真正關(guān)鍵的不同場景誤刪太小起不到去重作用。一個(gè)比較穩(wěn)的做法是按照不確定參數(shù)的偏差范圍做歸一化即每個(gè)維度除以其偏差量綱后再算歐氏距離。比如光伏偏差20 kW、負(fù)荷偏差30 kW那就把光伏場景值除以20、負(fù)荷除以30再做距離判斷。3.5 主循環(huán)迭代控制整個(gè)算法的主循環(huán)如下% 初始化 scenario_pool {}; LB -inf; UB inf; max_iter 20; tol 0.01; for iter 1:max_iter % 1. 求解主問題當(dāng)前場景池得到第一階段決策x和eta [x_opt, eta_opt] solve_master_problem(scenario_pool); LB max(LB, value(eta_opt)); % 主問題得到的是下界 % 2. 固定x_opt求解子問題 [worst_cost, worst_pv, worst_load] solve_subproblem(x_opt, data); UB min(UB, value(worst_cost)); % 子問題得到的是上界 fprintf(迭代 %d: LB%.2f, UB%.2f, gap%.4f\n, ... iter, LB, UB, abs(UB-LB)/abs(UB)); % 3. 判斷收斂 if abs(UB - LB) / abs(UB) tol break; end % 4. 更新關(guān)鍵場景池 candidate.cost value(worst_cost); candidate.pv value(worst_pv); candidate.load value(worst_load); [scenario_pool, ~] update_scenario_pool(scenario_pool, candidate, UB, LB, tol); end這里有個(gè)關(guān)于上下界關(guān)系的細(xì)節(jié)標(biāo)準(zhǔn)CCG中主問題的目標(biāo)值是下界子問題的目標(biāo)值是上界。因?yàn)橹鲉栴}只考慮了有限的場景可行域比真實(shí)問題松弛或者說約束不足所以目標(biāo)值偏小是下界子問題是給定第一階段決策后求最壞情況成本這個(gè)成本是實(shí)際可執(zhí)行的所以是上界。迭代的目的就是把下界不斷往上抬加場景加約束把上界不斷往下壓更好的第一階段決策直到兩者靠攏。4. 算例設(shè)計(jì)與結(jié)果分析4.1 測試系統(tǒng)參數(shù)我用一個(gè)改造的IEEE 13節(jié)點(diǎn)微網(wǎng)進(jìn)行測試參數(shù)如下微型燃?xì)廨啓C(jī)2臺(tái)額定功率分別為100 kW和150 kW燃料成本系數(shù)分別為0.45元/kWh和0.38元/kWh。儲(chǔ)能容量200 kWh最大充放電功率50 kW充放電效率均為0.95初始SOC為0.5。光伏額定功率200 kW預(yù)測曲線采用典型夏季晴天數(shù)據(jù)偏差取預(yù)測值的20%。負(fù)荷峰值負(fù)荷300 kW預(yù)測偏差取10%。分時(shí)電價(jià)峰時(shí)10:00-15:0018:00-21:001.2元/kWh谷時(shí)23:00-7:000.4元/kWh平時(shí)0.8元/kWh。魯棒預(yù)算 (\Gamma 8)即允許8個(gè)時(shí)段同時(shí)出現(xiàn)極端偏差。4.2 關(guān)鍵場景辨別 vs 標(biāo)準(zhǔn)CCG在相同參數(shù)下分別運(yùn)行標(biāo)準(zhǔn)CCG和關(guān)鍵場景辨別算法結(jié)果對比如下指標(biāo)標(biāo)準(zhǔn)CCG關(guān)鍵場景辨別算法迭代次數(shù)156總求解時(shí)間486 s187 s最終運(yùn)行成本上界3265.4 元3298.7 元與全場景枚舉的偏差-1.02%關(guān)鍵場景辨別算法用提高1%成本為代價(jià)換來了近3倍的求解速度提升。在實(shí)際工程中這個(gè)性價(jià)比是可接受的因?yàn)椴淮_定性本身也是近似建模的1%的精度損失相比計(jì)算時(shí)間的大幅下降完全值得。4.3 不同魯棒預(yù)算下的結(jié)果變化改變 (\Gamma) 的取值觀察運(yùn)行成本和魯棒性的權(quán)衡關(guān)系(\Gamma)運(yùn)行成本元最壞場景下棄負(fù)荷量kWh0確定性2898.5156.243056.362.483298.721.8123471.26.3163610.5024全極端3824.60隨著 (\Gamma) 增大運(yùn)行成本單調(diào)上升但系統(tǒng)面對最壞情況的應(yīng)對能力也在增強(qiáng)。(\Gamma8) 是一個(gè)甜點(diǎn)值成本增加約13.8%但最壞場景棄負(fù)荷量從156 kWh降到22 kWh降幅86%。繼續(xù)增大預(yù)算成本繼續(xù)漲但棄負(fù)荷量改善已經(jīng)不明顯說明邊際收益在遞減。這個(gè)結(jié)果也從側(cè)面驗(yàn)證了一個(gè)觀點(diǎn)魯棒優(yōu)化不是越保守越好。預(yù)算選得太大會(huì)讓成本高到離譜太小的預(yù)算又起不到保護(hù)作用?!昂线m的魯棒”才是工程上真正需要的。5. 常見問題與調(diào)試經(jīng)驗(yàn)實(shí)錄5.1 子問題雙線性項(xiàng)線性化失敗最常踩的坑。max-min子問題對偶化之后對偶變量乘不確定性變量會(huì)形成雙線性項(xiàng)。很多初學(xué)者代碼在這里直接報(bào)錯(cuò)或者說求解器報(bào)“non-convex”。我實(shí)測有效的一條經(jīng)驗(yàn)是先把對偶問題的約束整理成標(biāo)準(zhǔn)形式再線性化不要在內(nèi)層原問題里直接乘來乘去。另外M的取值要按數(shù)據(jù)量級來定大M太大會(huì)導(dǎo)致numerical issuesM太小導(dǎo)致解被錯(cuò)誤剪枝。我推薦的調(diào)試方式是先跑一個(gè)2時(shí)段的小規(guī)模算例把M的敏感度測一下用起來再放大到24時(shí)段。5.2 上下界不收斂或者震蕩如果迭代過程中UB和LB一直震蕩不收斂通常是兩種原因一是主問題場景數(shù)過多時(shí)求解出現(xiàn)數(shù)值穩(wěn)定性問題二是子問題的max問題沒有真正找對最惡劣場景。調(diào)試時(shí)先打印每一輪的場景和對應(yīng)子問題目標(biāo)值看看是不是存在目標(biāo)值幾乎相同但場景差別巨大的情況。如果是大概率是子問題求解器精度不夠或者大M線性化的M取值太小。把M調(diào)大兩倍再試一下很多時(shí)候就好了。還有一種情況是主問題的場景池更新太快把以前的關(guān)鍵場景刪掉了導(dǎo)致LB回退。我的處理方式是已經(jīng)加入過主問題的場景永遠(yuǎn)不刪除只是在篩選新場景時(shí)控制新增數(shù)量。這樣LB是單調(diào)不減的收斂軌跡更穩(wěn)。5.3 儲(chǔ)能SOC越界或者充放電同時(shí)為正這個(gè)問題多半出在約束遺漏。儲(chǔ)能SOC的上下界約束要在每個(gè)時(shí)段都顯式加上而且SOC的遞推要用嚴(yán)格等式不能松弛成不等式。充放電互斥如果不加約束可能會(huì)出現(xiàn)既充電又放電的“無效循環(huán)”白白增加成本。在Matlab里調(diào)試的時(shí)候我習(xí)慣把某個(gè)時(shí)段的SOC和充放電功率單獨(dú)拿出來打印對比肉眼檢查是否符合物理規(guī)律。如果出現(xiàn)充放電同為正認(rèn)真檢查互斥約束是否真的有效。5.4 CPLEX求解器報(bào)錯(cuò)或者求解極慢求解MILP時(shí)如果模型規(guī)模大求解器可能長時(shí)間無法找到可行解。一個(gè)經(jīng)驗(yàn)是給求解器設(shè)置合理的MIP gap和time limit給一個(gè)保守的可行解作為初始解。YALMIP支持在optimize函數(shù)里傳入sdpsettings(solver, cplex, cplex.mip.tolerances.mipgap, 0.001, cplex.timelimit, 300)這樣求解器不會(huì)在一個(gè)問題上耗死。還有一個(gè)點(diǎn)是場景池中場景數(shù)量達(dá)到一定規(guī)模后不要再繼續(xù)增加場景否則主問題的MILP規(guī)模會(huì)爆炸。這也是為什么我在關(guān)鍵場景辨別算法里做了Top-K截?cái)唷?.5 結(jié)果對場景初始池敏感關(guān)鍵場景辨別算法的收斂行為和初始場景的選擇有關(guān)。我建議初始場景池不要只放一個(gè)預(yù)測場景最穩(wěn)的做法是放四個(gè)基礎(chǔ)場景預(yù)測場景、光伏最低負(fù)荷最高、光伏最高負(fù)荷最低、光伏最低負(fù)荷最低覆蓋不確定集合的四個(gè)“角點(diǎn)”。這樣算法從第一輪迭代開始就有較好的邊界信息。5.6 Matlab版本與求解器兼容性我的代碼在Matlab R2022b YALMIP R20210430 CPLEX 12.10 上運(yùn)行穩(wěn)定。如果你用的是Matlab 2024之后的新版本記得檢查YALMIP的兼容性老版本的YALMIP在新版Matlab上偶爾會(huì)出現(xiàn)內(nèi)建函數(shù)命名沖突。CPLEX的版本和Matlab版本的兼容官方有文檔可查出了問題先去查版本對照表這個(gè)比瞎調(diào)代碼更高效。提示Matlab R2025、R2026之類的新版本里如果遇到License Manager的錯(cuò)誤先檢查環(huán)境變量和許可證配置通常跟算法本身沒關(guān)系網(wǎng)上搜索對應(yīng)報(bào)錯(cuò)信息就能搞定別一上來就懷疑程序?qū)戝e(cuò)了。6. 后續(xù)擴(kuò)展方向這套代碼的框架改一改就能適配不少變體問題。比如把光伏和負(fù)荷改成風(fēng)光負(fù)荷三個(gè)不確定源只需要在不確定性集合和子問題里多加一組變量把單微網(wǎng)擴(kuò)展成多微網(wǎng)互聯(lián)則需要把功率平衡約束改成帶聯(lián)絡(luò)線功率的多節(jié)點(diǎn)形式主問題的規(guī)模會(huì)大很多但算法框架不用變。另外一個(gè)值得嘗試的方向是分布式魯棒優(yōu)化Distributionally Robust OptimizationDRO。它把不確定性建模為模糊集合而不是確定集合需要用到Wasserstein距離來構(gòu)造模糊集。我在測試中發(fā)現(xiàn)DRO和兩階段魯棒在很多算例上結(jié)果差異不大但DRO的求解要更復(fù)雜需要調(diào)用專門的求解器。如果論文或者實(shí)際項(xiàng)目對保守度有硬性要求這個(gè)方向值得深入研究。還有一個(gè)小改進(jìn)關(guān)鍵場景辨別算法里的場景去重和排序目前是離線做的實(shí)時(shí)性要求高的場景下可以做在線版本利用上一次迭代的場景信息來加速本次迭代的初始場景池構(gòu)建這樣二次調(diào)度場景下比如日內(nèi)滾動(dòng)調(diào)度的效率還能再提一截。