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