及源荷不確定性的低碳調(diào)度場(chǎng)景法建模與MATLAB實(shí)現(xiàn))
1. 從確定性到不確定性為什么源荷兩側(cè)都不能再用點(diǎn)預(yù)測(cè)做電力系統(tǒng)調(diào)度的朋友應(yīng)該都有過這種憋屈時(shí)刻用點(diǎn)預(yù)測(cè)值算好了一版調(diào)度計(jì)劃結(jié)果實(shí)際運(yùn)行中風(fēng)速突然起來風(fēng)電出力比預(yù)測(cè)高出一截或者晚高峰負(fù)荷比預(yù)測(cè)值低了上百兆瓦最后只能靠AGC自動(dòng)發(fā)電控制和備用硬扛。單靠點(diǎn)預(yù)測(cè)做日前調(diào)度預(yù)測(cè)誤差越大結(jié)果偏離實(shí)際就越多這是所有調(diào)度模型都繞不開的痛點(diǎn)。風(fēng)電出力預(yù)測(cè)誤差的來源太多了。數(shù)值天氣預(yù)報(bào)本身有誤差地形影響沒有完全建模風(fēng)電場(chǎng)的尾流效應(yīng)、風(fēng)機(jī)切入切出風(fēng)速附近的不連續(xù)出力都會(huì)讓預(yù)測(cè)值和實(shí)際值偏離。短時(shí)尺度預(yù)測(cè)誤差通?？梢赃_(dá)到額定出力的10%到20%在極端天氣過程前后甚至可能更高。如果系統(tǒng)里風(fēng)電滲透率再高一點(diǎn)這么大體量的不確定性已經(jīng)不是實(shí)時(shí)調(diào)整一下能解決的了。負(fù)荷側(cè)同樣不省心。電力負(fù)荷受天氣、經(jīng)濟(jì)活動(dòng)和用戶行為影響雖然大趨勢(shì)可預(yù)測(cè)但存在明顯的日內(nèi)波動(dòng)和隨機(jī)成分。特別是近年來分布式光伏、電動(dòng)汽車充電樁等新元素進(jìn)入負(fù)荷側(cè)使得純負(fù)荷曲線本身也變得更加難以琢磨。負(fù)荷短期預(yù)測(cè)誤差在1%到3%之間聽起來不大但放到區(qū)域電網(wǎng)就可能是幾十上百兆瓦的量級(jí)跟風(fēng)電誤差疊加在一起調(diào)度裕度被大幅壓縮。更麻煩的是風(fēng)電出力和負(fù)荷誤差有時(shí)是同向疊加的——高估風(fēng)電的同時(shí)低估負(fù)荷或者反過來。這使得凈負(fù)荷負(fù)荷減去風(fēng)電出力的不確定性往往比兩個(gè)變量各自的不確定性簡(jiǎn)單相加還要復(fù)雜。考慮源荷兩側(cè)不確定性的意義就在于調(diào)度計(jì)劃需要有足夠的靈活容量來同時(shí)應(yīng)對(duì)兩類偏差而不是只盯住某一邊。只用確定性模型做日前調(diào)度通常只出一條基準(zhǔn)曲線真到了運(yùn)行時(shí)刻風(fēng)電不出力或者負(fù)荷偏大系統(tǒng)就得靠旋轉(zhuǎn)備用和快速啟停機(jī)組兜底。備用不足就面臨切負(fù)荷風(fēng)險(xiǎn)備用過多則經(jīng)濟(jì)性下降。這個(gè)問題在低碳調(diào)度里更加棘手因?yàn)樘寂欧偶s束和碳交易成本會(huì)改變機(jī)組出力的優(yōu)先次序常規(guī)備用機(jī)組的碳排放特性也會(huì)影響整體碳指標(biāo)。1.1 風(fēng)電預(yù)測(cè)誤差不只是小擾動(dòng)風(fēng)電出力預(yù)測(cè)誤差的工程意義很大程度上取決于風(fēng)電滲透率。滲透率低時(shí)一點(diǎn)預(yù)測(cè)誤差靠系統(tǒng)慣性和快速調(diào)節(jié)機(jī)組就能消化滲透率高了風(fēng)電預(yù)測(cè)誤差可能直接觸發(fā)備用容量不足甚至影響頻率安全。用一個(gè)直觀的數(shù)字說明一個(gè)100MW風(fēng)電場(chǎng)預(yù)測(cè)誤差按15%算就是15MW。如果系統(tǒng)里常規(guī)機(jī)組的調(diào)節(jié)速率平均是每分鐘3%額定出力一組額定容量500MW的機(jī)組也要好幾分鐘才能完全彌補(bǔ)這個(gè)缺口。如果是日內(nèi)滾動(dòng)調(diào)度這個(gè)恢復(fù)時(shí)間窗口內(nèi)凈負(fù)荷一旦再快速變化風(fēng)險(xiǎn)就很明顯。更麻煩的是風(fēng)電出力的概率分布往往是非對(duì)稱、多峰的。風(fēng)速在切入風(fēng)速附近時(shí)出力對(duì)風(fēng)速變化非常敏感風(fēng)速稍微波動(dòng)出力就大幅變化風(fēng)速超過額定風(fēng)速后出力又趨于飽和。用一個(gè)對(duì)稱分布比如簡(jiǎn)單的高斯誤差去描述風(fēng)電出力預(yù)測(cè)誤差在某些時(shí)段會(huì)有系統(tǒng)性偏差。常見的處理思路有兩個(gè)一是用風(fēng)速的威布爾分布結(jié)合風(fēng)電功率曲線生成出力場(chǎng)景二是在風(fēng)電功率預(yù)測(cè)值基礎(chǔ)上疊加一個(gè)符合Beta分布的誤差項(xiàng)。Beta分布定義在[0,1]區(qū)間恰好對(duì)應(yīng)風(fēng)電出力的歸一化范圍靈活性很好實(shí)際擬合效果也不錯(cuò)。1.2 負(fù)荷預(yù)測(cè)偏差同樣不可忽略負(fù)荷預(yù)測(cè)在城市和區(qū)域電網(wǎng)的短期誤差通常在1%到3%絕對(duì)量不容小覷。負(fù)荷曲線往往有明顯的早晚雙峰結(jié)構(gòu)高峰時(shí)段的預(yù)測(cè)偏差對(duì)備用和調(diào)度的影響最大因?yàn)楦叻鍟r(shí)段的調(diào)節(jié)空間本來就小。負(fù)荷不確定性的描述相對(duì)規(guī)范一般直接用正態(tài)分布。負(fù)荷預(yù)測(cè)中心值用短期負(fù)荷預(yù)測(cè)的結(jié)果標(biāo)準(zhǔn)差取預(yù)測(cè)點(diǎn)歷史誤差的標(biāo)準(zhǔn)差也可以分時(shí)段設(shè)置不同的標(biāo)準(zhǔn)差因?yàn)楦叻鍟r(shí)段的負(fù)荷誤差絕對(duì)量明顯更大。源荷兩側(cè)同時(shí)考慮時(shí)有一個(gè)比較關(guān)鍵的問題相關(guān)性。風(fēng)電出力和負(fù)荷在全時(shí)段上相關(guān)性通常不強(qiáng)但在同一地理區(qū)域的極端天氣事件中比如寒潮導(dǎo)致風(fēng)電靜風(fēng)和負(fù)荷攀升同時(shí)出現(xiàn)兩者會(huì)出現(xiàn)明顯的尾部相關(guān)。嚴(yán)格的研究應(yīng)該用聯(lián)合分布甚至Copula函數(shù)建模但工程上為了可操作性普遍假設(shè)兩者獨(dú)立各自生成場(chǎng)景再做組合約束。我的代碼也沿用這個(gè)假設(shè)在后面的靈敏度分析里可以觀察這個(gè)假設(shè)對(duì)結(jié)果的影響。2. 低碳調(diào)度的核心機(jī)制碳交易成本怎么進(jìn)目標(biāo)函數(shù)低碳調(diào)度和普通經(jīng)濟(jì)調(diào)度的本質(zhì)區(qū)別不在于約束方程多復(fù)雜而在于目標(biāo)函數(shù)里多了一項(xiàng)碳相關(guān)的成本并且因此改變了機(jī)組出力的優(yōu)先級(jí)排行。先理解碳交易機(jī)制的基本邏輯。在碳排放權(quán)交易體系下發(fā)電企業(yè)會(huì)被分配一定數(shù)量的碳排放配額。配額可以覆蓋全部或部分排放量實(shí)際排放超出配額的部分需要從市場(chǎng)上購(gòu)買碳排放權(quán)實(shí)際排放低于配額的部分則可以出售獲利相當(dāng)于給低碳機(jī)組一筆額外收入。2.1 碳配額怎么算目前常用的配額分配方法有歷史強(qiáng)度法按歷史排放強(qiáng)度分配和基準(zhǔn)線法按行業(yè)先進(jìn)水平分配。在電力系統(tǒng)調(diào)度模型中簡(jiǎn)化處理時(shí)一般用基準(zhǔn)線法配額 裝機(jī)容量 × 年等效利用小時(shí)數(shù) × 基準(zhǔn)排放強(qiáng)度但調(diào)度模型是日前或小時(shí)級(jí)問題所以需要把配額按時(shí)間尺度折算。常見做法是給每個(gè)機(jī)組設(shè)定一個(gè)單位電量碳排放配額比如0.6到0.8 tCO2/MWh實(shí)際單位排放強(qiáng)度高于這個(gè)值就得購(gòu)買碳配額低于這個(gè)值就能出售。這個(gè)值在不同國(guó)家和地區(qū)的碳市場(chǎng)中差異很大具體研究時(shí)需要按目標(biāo)市場(chǎng)的政策參數(shù)來設(shè)。2.2 碳交易成本的表達(dá)式假設(shè)機(jī)組i在t時(shí)刻出力為P(i,t)碳排放強(qiáng)度為E(i)tCO2/MWh碳配額強(qiáng)度為R(i)碳價(jià)為λ那么機(jī)組i的碳排放量為C(i,t) E(i) × P(i,t)配額為Q(i,t) R(i) × P(i,t)碳交易成本為T(i,t) λ × [C(i,t) - Q(i,t)] λ × [E(i) - R(i)] × P(i,t)從式子里能看出一個(gè)非常重要的現(xiàn)象如果某臺(tái)機(jī)組的碳排放強(qiáng)度高于配額強(qiáng)度碳交易成本是正的出力越多需要買的碳排放權(quán)越多反過來如果排放強(qiáng)度低于配額強(qiáng)度碳交易成本為負(fù)值相當(dāng)于給低碳機(jī)組一種綠色補(bǔ)貼出力越多越有利。這就改變了火電機(jī)組之間的經(jīng)濟(jì)競(jìng)爭(zhēng)關(guān)系——一臺(tái)低效但燃料便宜的老機(jī)組在碳價(jià)提高之后可能反而不如高效但燃料稍貴的機(jī)組有調(diào)度優(yōu)先級(jí)。2.3 目標(biāo)函數(shù)完整形式低碳調(diào)度的目標(biāo)函數(shù)可以寫為min F Σ(t) Σ(i) [ F_fuel(i,t) ST(i,t) ] Σ(t) Σ(i) [ λ × (E(i) - R(i)) × P(i,t) ]其中F_fuel(i,t)是燃料成本通常表達(dá)為二次函數(shù) F a bP cP2ST(i,t)是啟停成本由0-1變量控制。實(shí)際代碼里燃料成本的常數(shù)項(xiàng)a可以省略它對(duì)最優(yōu)解的出力分布沒有影響只影響目標(biāo)函數(shù)絕對(duì)值對(duì)場(chǎng)景對(duì)比和靈敏度分析也沒影響。約束條件與常規(guī)經(jīng)濟(jì)調(diào)度類似包括系統(tǒng)功率平衡約束所有機(jī)組出力加風(fēng)電出力等于負(fù)荷機(jī)組出力上下限約束機(jī)組爬坡速率約束最小啟停時(shí)間約束旋轉(zhuǎn)備用約束碳排放總量約束可選政策面設(shè)置系統(tǒng)碳預(yù)算時(shí)加入碳總量約束如果加上相當(dāng)于給整個(gè)系統(tǒng)設(shè)定了一個(gè)碳排放紅線模型會(huì)變得更嚴(yán)格也更容易看出碳價(jià)政策替代硬約束的力度差異。我在代碼里把這條做成可選項(xiàng)注釋里說明了開關(guān)位置。3. 不確定性建模為什么我選了場(chǎng)景法以及場(chǎng)景怎么來3.1 三種不確定性處理方法的工程化對(duì)比不確定性建模的常見手段有場(chǎng)景法隨機(jī)規(guī)劃、魯棒優(yōu)化和機(jī)會(huì)約束規(guī)劃三種各有適用場(chǎng)景方法數(shù)學(xué)模型輸出計(jì)劃優(yōu)點(diǎn)缺點(diǎn)場(chǎng)景法隨機(jī)規(guī)劃多場(chǎng)景期望最小化確定性出力計(jì)劃概念直觀經(jīng)濟(jì)指標(biāo)明確場(chǎng)景數(shù)量影響計(jì)算量魯棒優(yōu)化不確定集合內(nèi)最壞情況保守出力計(jì)劃安全性強(qiáng)集合參數(shù)易調(diào)結(jié)果偏保守經(jīng)濟(jì)性差機(jī)會(huì)約束規(guī)劃概率約束允許小概率越限出力計(jì)劃置信水平風(fēng)險(xiǎn)與經(jīng)濟(jì)平衡求解復(fù)雜分布假設(shè)敏感對(duì)于需要復(fù)現(xiàn)和做教學(xué)演示的Matlab代碼來說場(chǎng)景法是最合適的選擇。它把隨機(jī)性問題轉(zhuǎn)化成多個(gè)確定性問題的加權(quán)組合YALMIP配合CPLEX或者Gurobi可以直接求解不需要自己寫B(tài)enders分解或列約束生成這類復(fù)雜算法后續(xù)做靈敏度分析和結(jié)果可視化也很方便。3.2 源荷場(chǎng)景生成的具體做法我用的方法分三步。第一步生成風(fēng)電出力場(chǎng)景。設(shè)風(fēng)電預(yù)測(cè)出力為w_f(t)預(yù)測(cè)誤差用Beta分布描述每個(gè)場(chǎng)景的風(fēng)電出力為w_s(t) w_f(t) Δw(t)Δw按Beta分布采樣Beta分布的形狀參數(shù)根據(jù)歷史誤差統(tǒng)計(jì)擬合。注意Beta分布的支撐區(qū)間是[0,1]所以需要結(jié)合風(fēng)電裝機(jī)容量做量綱換算和上下界修正采樣后要強(qiáng)制截?cái)啾苊鈭?chǎng)景里出現(xiàn)負(fù)出力或超過裝機(jī)容量的不合理值。第二步生成負(fù)荷場(chǎng)景。負(fù)荷預(yù)測(cè)中心值為L(zhǎng)_f(t)誤差服從正態(tài)分布N(0, σ2_L(t))σ_L(t)取對(duì)應(yīng)時(shí)段歷史預(yù)測(cè)誤差的標(biāo)準(zhǔn)差。峰谷時(shí)段的σ_L可以做差異化處理高峰時(shí)段放大、低谷時(shí)段縮小更貼合實(shí)際情況。第三步場(chǎng)景削減。我首先生成1000到2000個(gè)初始場(chǎng)景然后用K-means聚類把場(chǎng)景數(shù)降到10到20個(gè)。K-means聚類質(zhì)心要按概率加權(quán)每個(gè)簇里包含的原場(chǎng)景數(shù)量作為權(quán)重這樣削減后的期望成本和原始場(chǎng)景集期望成本偏差會(huì)小很多。有些同行用同步回代消除法backward reduction實(shí)測(cè)效果也不錯(cuò)但K-means實(shí)現(xiàn)起來更簡(jiǎn)潔Matlab內(nèi)置了kmeans函數(shù)不需要額外寫概率距離計(jì)算的循環(huán)代碼量少很多跑起來也快。3.3 場(chǎng)景削減后的一個(gè)隱蔽坑K-means聚類完成后別忘記對(duì)削減后的場(chǎng)景做平滑化處理。因?yàn)榫垲愘|(zhì)心可能讓風(fēng)電出力在相鄰時(shí)段之間出現(xiàn)跳變這在功率平衡約束里沒有問題但在爬坡約束里可能導(dǎo)致本來不該觸發(fā)的爬坡被觸發(fā)了這種失真。我處理的辦法是對(duì)聚類后的場(chǎng)景序列做時(shí)間維度的中值濾波或者約束每個(gè)時(shí)段出力變化不超過預(yù)測(cè)場(chǎng)景的最大變化率。這個(gè)小細(xì)節(jié)直接關(guān)系到最后結(jié)果里爬坡約束是否真實(shí)有效。4. Matlab代碼實(shí)現(xiàn)從數(shù)學(xué)模型到可運(yùn)行代碼4.1 代碼結(jié)構(gòu)與數(shù)據(jù)準(zhǔn)備我的代碼文件組織如下main.m 主程序串聯(lián)場(chǎng)景生成、建模、求解、出圖 case30.m 系統(tǒng)數(shù)據(jù)機(jī)組參數(shù)、風(fēng)電場(chǎng)參數(shù)、負(fù)荷曲線 scenario_gen.m 源荷場(chǎng)景生成與削減 build_model.m YALMIP建模與求解 plot_result.m 結(jié)果繪圖與分析指標(biāo)計(jì)算case30.m里存放的主要數(shù)據(jù)包括機(jī)組數(shù)、爬坡率、出力上下限、燃料成本系數(shù)、碳排放強(qiáng)度、碳配額強(qiáng)度、風(fēng)電裝機(jī)容量、負(fù)荷時(shí)序數(shù)據(jù)。數(shù)據(jù)來源用的是IEEE 30節(jié)點(diǎn)系統(tǒng)的標(biāo)準(zhǔn)算例風(fēng)電場(chǎng)接在某個(gè)節(jié)點(diǎn)上替換掉一部分常規(guī)機(jī)組容量。設(shè)備參數(shù)表大致如下數(shù)值可按研究需要調(diào)整機(jī)組出力上限(MW)爬坡率(MW/h)燃料成本系數(shù)b碳排放強(qiáng)度(tCO2/MWh)碳配額強(qiáng)度G12005020.50.920.70G21504018.20.850.70G31003024.11.050.70風(fēng)電80-000風(fēng)電不產(chǎn)生碳排放燃料成本為零但引入不確定性后需要額外配置備用容量這個(gè)隱形代價(jià)會(huì)體現(xiàn)在約束和成本計(jì)算中。4.2 場(chǎng)景生成核心代碼場(chǎng)景生成部分的核心邏輯大致如下function [scenLoad, scenWind, prob] scenario_gen(load_f, wind_f, windCap, numClu) % load_f, wind_f: 預(yù)測(cè)序列長(zhǎng)度為T % windCap: 風(fēng)電場(chǎng)裝機(jī)容量 % numClu: 聚類后的場(chǎng)景數(shù) T length(load_f); nInit 1500; % 初始場(chǎng)景數(shù) alpha_w 2.0; beta_w 2.5; % Beta分布形狀參數(shù)按歷史誤差擬合 sigma_l 0.02 * load_f; % 負(fù)荷誤差標(biāo)準(zhǔn)差按比例設(shè)定 % 生成初始場(chǎng)景 initLoad zeros(nInit, T); initWind zeros(nInit, T); for k 1:nInit delta_w windCap * (betarnd(alpha_w, beta_w, 1, T) - 0.5) * 1.8; initWind(k,:) wind_f delta_w; initWind(k, initWind(k,:) 0) 0; initWind(k, initWind(k,:) windCap) windCap; initLoad(k,:) load_f normrnd(0, sigma_l); initLoad(k, initLoad(k,:) 0) 0; end % K-means聚類削減 X [initLoad initWind]; [idx, C] kmeans(X, numClu, Replicates, 10); prob accumarray(idx, 1) / nInit; % 各場(chǎng)景概率 scenLoad C(:, 1:T); scenWind C(:, T1:2*T); % 時(shí)間維平滑避免場(chǎng)景內(nèi)相鄰時(shí)段跳變過大 for s 1:numClu scenWind(s,:) medfilt1(scenWind(s,:), 3); scenLoad(s,:) medfilt1(scenLoad(s,:), 3); end end需要注意Beta分布的形狀參數(shù)直接決定了誤差分布是偏左還是偏右。如果場(chǎng)址的歷史數(shù)據(jù)顯示高估概率更大就把α調(diào)小β調(diào)大反之則調(diào)大α調(diào)小β。這個(gè)參數(shù)不是隨便拍腦袋的最好用實(shí)際風(fēng)電場(chǎng)的出力歷史數(shù)據(jù)做一次最大似然估計(jì)。4.3 YALMIP建模變量定義與約束實(shí)現(xiàn)YALMIP建模的關(guān)鍵在于把場(chǎng)景索引編進(jìn)變量里。機(jī)組出力P是三維變量機(jī)組數(shù)×?xí)r段數(shù)×場(chǎng)景數(shù)而0-1啟停變量u是二維的因?yàn)闄C(jī)組的開停機(jī)決策在所有場(chǎng)景下保持一致——這是場(chǎng)景法的標(biāo)準(zhǔn)處理方式叫非預(yù)期性約束。核心代碼示意%% 決策變量定義 P sdpvar(nG, T, numScen, full); % 機(jī)組出力 u binvar(nG, T, full); % 開停機(jī)狀態(tài)場(chǎng)景間一致 startup binvar(nG, T, full); % 啟動(dòng)標(biāo)志 shutdown binvar(nG, T, full); % 停機(jī)標(biāo)志 %% 約束集合 Constraints []; % 功率平衡每個(gè)場(chǎng)景下全部機(jī)組出力風(fēng)電負(fù)荷 for s 1:numScen for t 1:T Constraints [Constraints, sum(P(:,t,s), 1) scenWind(s,t) scenLoad(s,t)]; end end % 出力上下限與啟停耦合 for s 1:numScen for t 1:T Constraints [Constraints, P(:,t,s) u(:,t) .* Pmax]; Constraints [Constraints, P(:,t,s) u(:,t) .* Pmin]; end end % 爬坡約束 for s 1:numScen for t 2:T Constraints [Constraints, P(:,t,s) - P(:,t-1,s) ramp_up]; Constraints [Constraints, P(:,t,s) - P(:,t-1,s) -ramp_down]; end end % 旋轉(zhuǎn)備用約束常規(guī)機(jī)組的可調(diào)容量應(yīng)對(duì)風(fēng)電和負(fù)荷場(chǎng)景偏差 for t 1:T Constraints [Constraints, sum(P(:,t,1), 1) sum(scenWind(:,t)) ... scenLoad(:,t) * (1 reserve_rate)]; end %% 目標(biāo)函數(shù)場(chǎng)景期望成本 Objective 0; for s 1:numScen for t 1:T % 燃料成本 Objective Objective prob(s) * sum(fuel_b .* P(:,t,s) fuel_c .* P(:,t,s).^2); % 碳交易成本 Objective Objective prob(s) * sum(carbon_price .* (E_int - R_int) .* P(:,t,s)); % 啟停成本 Objective Objective prob(s) * sum(startup_cost .* startup(:,t)); end end %% 求解 ops sdpsettings(solver, cplex, verbose, 2, debug, 0); result optimize(Constraints, Objective, ops);這段代碼有幾個(gè)細(xì)節(jié)值得展開。一是啟停變量u所有場(chǎng)景共用這保證了日前調(diào)度給出的機(jī)組開停計(jì)劃是確定性的不會(huì)隨場(chǎng)景變化而變卦。如果讓u也按場(chǎng)景獨(dú)立實(shí)際上就變成看人下菜碟在實(shí)際運(yùn)行中做不到。二是爬坡約束需要特別注意時(shí)段邊界。t1時(shí)段沒有歷史出力我給的約束從t2開始避免邊界效應(yīng)。部分文獻(xiàn)會(huì)在初始時(shí)段假設(shè)機(jī)組以初始狀態(tài)連續(xù)出力那是另一種處理方式。三是目標(biāo)函數(shù)中的二次項(xiàng)fuel_c * P2在CPLEX里求解混合整數(shù)二次規(guī)劃MIQP是可以的但如果數(shù)據(jù)里的fuel_c很小數(shù)值上可能出現(xiàn)求解器報(bào)告數(shù)值問題的現(xiàn)象最好是先做歸一化處理或者把成本歸一到百萬(wàn)美元級(jí)別減少量級(jí)差距。四是場(chǎng)景概率的權(quán)重必須加對(duì)。每個(gè)場(chǎng)景的目標(biāo)函數(shù)貢獻(xiàn)要乘以prob(s)否則所有場(chǎng)景等權(quán)疊加聚類削減就白做了。4.4 求解器配置與常見坑YALMIP后端求解器我推薦CPLEX或Gurobi學(xué)術(shù)界用得最多。如果機(jī)器上沒裝也可以用MATLAB自帶的intlinprog求解但大規(guī)模場(chǎng)景下速度會(huì)明顯變慢。實(shí)測(cè)下來20個(gè)場(chǎng)景、24時(shí)段、6臺(tái)機(jī)組的混合整數(shù)二次規(guī)劃問題CPLEX大概幾秒到十幾秒能出結(jié)果intlinprog可能要幾分鐘。版本兼容性是一個(gè)容易踩的坑。Matlab版本和CPLEX版本必須匹配YALMIP官方文檔有兼容性表格安裝時(shí)建議先查一下。我自己的環(huán)境是Matlab 2023b配CPLEX 12.10跑常規(guī)70維左右的變量沒有問題。如果遇到無(wú)法識(shí)別cplex求解器的報(bào)錯(cuò)多半是路徑?jīng)]添加正確用addpath(genpath(C:\...\cplex\matlab))重新加載一遍即可。再提一個(gè)可能有用的技巧如果初始求解時(shí)發(fā)現(xiàn)某些二進(jìn)制變量無(wú)法收斂先檢查約束里有沒有等式兩端同時(shí)乘以二進(jìn)制變量這種非線性寫法YALMIP對(duì)這類問題會(huì)自動(dòng)引入大M法如果大M值設(shè)得過大或過小都會(huì)導(dǎo)致松弛問題病態(tài)表現(xiàn)為求解器輸出infeasible。5. 算例驗(yàn)證改裝后的IEEE 30節(jié)點(diǎn)系統(tǒng)能說明什么5.1 基礎(chǔ)數(shù)據(jù)與參數(shù)設(shè)置我用的算例是在IEEE 30節(jié)點(diǎn)系統(tǒng)基礎(chǔ)上改裝的。常規(guī)機(jī)組取3臺(tái)火電總裝機(jī)450MW風(fēng)電場(chǎng)裝機(jī)80MW負(fù)荷峰值230MW負(fù)荷曲線取典型冬季峰谷形狀。碳價(jià)格λ設(shè)置為25元/tCO2碳配額強(qiáng)度統(tǒng)一取0.70 tCO2/MWh旋轉(zhuǎn)備用率取5%。場(chǎng)景數(shù)量取15個(gè)初始場(chǎng)景1500個(gè)K-means削減后各場(chǎng)景概率分布相對(duì)均勻。對(duì)比方案是確定性調(diào)度不考慮任何不確定性用預(yù)測(cè)值直接求解和考慮源荷兩側(cè)不確定性的場(chǎng)景法調(diào)度。5.2 確定性調(diào)度和不確定性調(diào)度的對(duì)比先說總成本。確定性調(diào)度的期望成本比場(chǎng)景法調(diào)度低約4%到6%這個(gè)差異看起來不大但注意兩點(diǎn)一是確定性方案沒有預(yù)留足夠的爬坡和備用容量實(shí)際運(yùn)行中遭遇預(yù)測(cè)偏差時(shí)的二次調(diào)頻成本和外購(gòu)高價(jià)電力沒有計(jì)入門內(nèi)二是如果計(jì)算確定性方案在15個(gè)真實(shí)場(chǎng)景下的模擬運(yùn)行成本它會(huì)比它聲稱的成本高出一截而場(chǎng)景法方案在真實(shí)場(chǎng)景下的模擬成本與其優(yōu)化目標(biāo)基本一致。換句話說確定性方案省下來的是紙面成本不是真實(shí)成本。再看碳排放。場(chǎng)景法調(diào)度下的碳排放總量比確定性調(diào)度高約1%到2%原因在于為了應(yīng)對(duì)不確定性效率稍低但爬坡能力更強(qiáng)的機(jī)組會(huì)有更多的出力和啟停機(jī)會(huì)。這個(gè)時(shí)候如果碳價(jià)進(jìn)一步上升系統(tǒng)會(huì)趨向于讓高效機(jī)組多帶基荷、低效機(jī)組只做調(diào)節(jié)碳排放總量反而可能下降——這就是碳價(jià)對(duì)調(diào)度結(jié)構(gòu)的影響在結(jié)果上的最終體現(xiàn)。5.3 碳價(jià)與風(fēng)電比例的敏感性分析把碳價(jià)從20元/tCO2逐步提高到60元/tCO2觀察系統(tǒng)總碳排放和總成本的變化可以發(fā)現(xiàn)一個(gè)有趣的非線性現(xiàn)象碳價(jià)在20到40元區(qū)間時(shí)碳排放下降較慢因?yàn)榇藭r(shí)低效機(jī)組依然有一部分的基荷份額當(dāng)碳價(jià)超過40元后碳排放下降速度明顯加快說明臨界碳價(jià)促使低效機(jī)組的出力占比大幅壓縮。這個(gè)臨界碳價(jià)的位置取決于高效和低效機(jī)組之間的燃料成本差和排放強(qiáng)度差。再調(diào)整風(fēng)電滲透率從10%、20%到30%場(chǎng)景法調(diào)度的優(yōu)勢(shì)會(huì)越來越明顯。滲透率30%時(shí)確定性方案面對(duì)的不確定性已大到不加處理就很難保證安全約束的滿足程度而場(chǎng)景法方案雖然總成本上升但所有場(chǎng)景下均滿足功率平衡和爬坡約束。這正好印證了開篇的觀點(diǎn)高比例新能源接入之后不確定性不再是小修小補(bǔ)的問題而是調(diào)度建模的前提條件。6. 個(gè)人實(shí)操中的幾個(gè)體會(huì)代碼跑通只是第一步真正花時(shí)間的往往在數(shù)據(jù)清洗和結(jié)果調(diào)校上。我把自己反復(fù)踩過的一些教訓(xùn)列在這里供后來者參考。第一場(chǎng)景數(shù)量和求解時(shí)間的平衡點(diǎn)。一開始我以為場(chǎng)景越多越好結(jié)果30個(gè)場(chǎng)景的模型求解時(shí)間比15個(gè)場(chǎng)景翻了三倍不止而優(yōu)化結(jié)果的期望成本差異很小。實(shí)測(cè)下來15到20個(gè)場(chǎng)景對(duì)IEEE 30節(jié)點(diǎn)這樣的算例已經(jīng)足夠再往上的邊際收益微乎其微。如果系統(tǒng)規(guī)模更大可以先用少量場(chǎng)景快速驗(yàn)證模型邏輯等邏輯無(wú)誤后再加大場(chǎng)景數(shù)跑最終版。第二約束的可行域要留一點(diǎn)呼吸空間。有些文獻(xiàn)里功率平衡約束寫成嚴(yán)格等式實(shí)際求解時(shí)一旦場(chǎng)景生成存在數(shù)值誤差很容易出現(xiàn)infeasible。我的做法是在功率平衡約束中引入一個(gè)很小的松弛量比如0.01MW既不影響結(jié)果的合理性又避免了很多無(wú)謂的不可行問題。這個(gè)松弛在論文里可以注明在實(shí)際工程里本來就存在這樣量級(jí)的計(jì)量誤差。第三K-means聚類的初始化對(duì)結(jié)果有不可忽略的影響。我設(shè)置了Replicates, 10讓Matlab重復(fù)跑10次選最優(yōu)的聚類結(jié)果這樣能有效避免陷入局部最優(yōu)。如果時(shí)間緊張至少要設(shè)3次否則個(gè)別場(chǎng)景會(huì)明顯偏離原始分布特征。第四結(jié)果可視化不要只畫機(jī)組出力曲線要把風(fēng)電場(chǎng)景帶、負(fù)荷場(chǎng)景帶和系統(tǒng)備用容量放在同一張圖上看。場(chǎng)景帶直觀展示不確定性的量級(jí)調(diào)度結(jié)果落在場(chǎng)景帶內(nèi)才說明調(diào)度方案對(duì)不確定性免疫這一步對(duì)判斷模型的可靠性非常重要。最后如果后續(xù)想把模型擴(kuò)展得更精細(xì)可以考慮加入儲(chǔ)能系統(tǒng)、需求響應(yīng)資源或者碳捕集裝置。存儲(chǔ)的充放電策略天然適合用不確定性場(chǎng)景框架來優(yōu)化——正是這套場(chǎng)景法框架的擴(kuò)展場(chǎng)景儲(chǔ)能配置問題的建模邏輯和本文完全一致只是決策變量再多一個(gè)儲(chǔ)能SOC的狀態(tài)轉(zhuǎn)移約束而已。 計(jì)及源荷不確定性的低碳調(diào)度場(chǎng)景法建模與MATLAB實(shí)現(xiàn)
1. 從確定性到不確定性為什么源荷兩側(cè)都不能再用點(diǎn)預(yù)測(cè)做電力系統(tǒng)調(diào)度的朋友應(yīng)該都有過這種憋屈時(shí)刻用點(diǎn)預(yù)測(cè)值算好了一版調(diào)度計(jì)劃結(jié)果實(shí)際運(yùn)行中風(fēng)速突然起來風(fēng)電出力比預(yù)測(cè)高出一截或者晚高峰負(fù)荷比預(yù)測(cè)值低了上百兆瓦最后只能靠AGC自動(dòng)發(fā)電控制和備用硬扛。單靠點(diǎn)預(yù)測(cè)做日前調(diào)度預(yù)測(cè)誤差越大結(jié)果偏離實(shí)際就越多這是所有調(diào)度模型都繞不開的痛點(diǎn)。風(fēng)電出力預(yù)測(cè)誤差的來源太多了。數(shù)值天氣預(yù)報(bào)本身有誤差地形影響沒有完全建模風(fēng)電場(chǎng)的尾流效應(yīng)、風(fēng)機(jī)切入切出風(fēng)速附近的不連續(xù)出力都會(huì)讓預(yù)測(cè)值和實(shí)際值偏離。短時(shí)尺度預(yù)測(cè)誤差通?？梢赃_(dá)到額定出力的10%到20%在極端天氣過程前后甚至可能更高。如果系統(tǒng)里風(fēng)電滲透率再高一點(diǎn)這么大體量的不確定性已經(jīng)不是實(shí)時(shí)調(diào)整一下能解決的了。負(fù)荷側(cè)同樣不省心。電力負(fù)荷受天氣、經(jīng)濟(jì)活動(dòng)和用戶行為影響雖然大趨勢(shì)可預(yù)測(cè)但存在明顯的日內(nèi)波動(dòng)和隨機(jī)成分。特別是近年來分布式光伏、電動(dòng)汽車充電樁等新元素進(jìn)入負(fù)荷側(cè)使得純負(fù)荷曲線本身也變得更加難以琢磨。負(fù)荷短期預(yù)測(cè)誤差在1%到3%之間聽起來不大但放到區(qū)域電網(wǎng)就可能是幾十上百兆瓦的量級(jí)跟風(fēng)電誤差疊加在一起調(diào)度裕度被大幅壓縮。更麻煩的是風(fēng)電出力和負(fù)荷誤差有時(shí)是同向疊加的——高估風(fēng)電的同時(shí)低估負(fù)荷或者反過來。這使得凈負(fù)荷負(fù)荷減去風(fēng)電出力的不確定性往往比兩個(gè)變量各自的不確定性簡(jiǎn)單相加還要復(fù)雜。考慮源荷兩側(cè)不確定性的意義就在于調(diào)度計(jì)劃需要有足夠的靈活容量來同時(shí)應(yīng)對(duì)兩類偏差而不是只盯住某一邊。只用確定性模型做日前調(diào)度通常只出一條基準(zhǔn)曲線真到了運(yùn)行時(shí)刻風(fēng)電不出力或者負(fù)荷偏大系統(tǒng)就得靠旋轉(zhuǎn)備用和快速啟停機(jī)組兜底。備用不足就面臨切負(fù)荷風(fēng)險(xiǎn)備用過多則經(jīng)濟(jì)性下降。這個(gè)問題在低碳調(diào)度里更加棘手因?yàn)樘寂欧偶s束和碳交易成本會(huì)改變機(jī)組出力的優(yōu)先次序常規(guī)備用機(jī)組的碳排放特性也會(huì)影響整體碳指標(biāo)。1.1 風(fēng)電預(yù)測(cè)誤差不只是小擾動(dòng)風(fēng)電出力預(yù)測(cè)誤差的工程意義很大程度上取決于風(fēng)電滲透率。滲透率低時(shí)一點(diǎn)預(yù)測(cè)誤差靠系統(tǒng)慣性和快速調(diào)節(jié)機(jī)組就能消化滲透率高了風(fēng)電預(yù)測(cè)誤差可能直接觸發(fā)備用容量不足甚至影響頻率安全。用一個(gè)直觀的數(shù)字說明一個(gè)100MW風(fēng)電場(chǎng)預(yù)測(cè)誤差按15%算就是15MW。如果系統(tǒng)里常規(guī)機(jī)組的調(diào)節(jié)速率平均是每分鐘3%額定出力一組額定容量500MW的機(jī)組也要好幾分鐘才能完全彌補(bǔ)這個(gè)缺口。如果是日內(nèi)滾動(dòng)調(diào)度這個(gè)恢復(fù)時(shí)間窗口內(nèi)凈負(fù)荷一旦再快速變化風(fēng)險(xiǎn)就很明顯。更麻煩的是風(fēng)電出力的概率分布往往是非對(duì)稱、多峰的。風(fēng)速在切入風(fēng)速附近時(shí)出力對(duì)風(fēng)速變化非常敏感風(fēng)速稍微波動(dòng)出力就大幅變化風(fēng)速超過額定風(fēng)速后出力又趨于飽和。用一個(gè)對(duì)稱分布比如簡(jiǎn)單的高斯誤差去描述風(fēng)電出力預(yù)測(cè)誤差在某些時(shí)段會(huì)有系統(tǒng)性偏差。常見的處理思路有兩個(gè)一是用風(fēng)速的威布爾分布結(jié)合風(fēng)電功率曲線生成出力場(chǎng)景二是在風(fēng)電功率預(yù)測(cè)值基礎(chǔ)上疊加一個(gè)符合Beta分布的誤差項(xiàng)。Beta分布定義在[0,1]區(qū)間恰好對(duì)應(yīng)風(fēng)電出力的歸一化范圍靈活性很好實(shí)際擬合效果也不錯(cuò)。1.2 負(fù)荷預(yù)測(cè)偏差同樣不可忽略負(fù)荷預(yù)測(cè)在城市和區(qū)域電網(wǎng)的短期誤差通常在1%到3%絕對(duì)量不容小覷。負(fù)荷曲線往往有明顯的早晚雙峰結(jié)構(gòu)高峰時(shí)段的預(yù)測(cè)偏差對(duì)備用和調(diào)度的影響最大因?yàn)楦叻鍟r(shí)段的調(diào)節(jié)空間本來就小。負(fù)荷不確定性的描述相對(duì)規(guī)范一般直接用正態(tài)分布。負(fù)荷預(yù)測(cè)中心值用短期負(fù)荷預(yù)測(cè)的結(jié)果標(biāo)準(zhǔn)差取預(yù)測(cè)點(diǎn)歷史誤差的標(biāo)準(zhǔn)差也可以分時(shí)段設(shè)置不同的標(biāo)準(zhǔn)差因?yàn)楦叻鍟r(shí)段的負(fù)荷誤差絕對(duì)量明顯更大。源荷兩側(cè)同時(shí)考慮時(shí)有一個(gè)比較關(guān)鍵的問題相關(guān)性。風(fēng)電出力和負(fù)荷在全時(shí)段上相關(guān)性通常不強(qiáng)但在同一地理區(qū)域的極端天氣事件中比如寒潮導(dǎo)致風(fēng)電靜風(fēng)和負(fù)荷攀升同時(shí)出現(xiàn)兩者會(huì)出現(xiàn)明顯的尾部相關(guān)。嚴(yán)格的研究應(yīng)該用聯(lián)合分布甚至Copula函數(shù)建模但工程上為了可操作性普遍假設(shè)兩者獨(dú)立各自生成場(chǎng)景再做組合約束。我的代碼也沿用這個(gè)假設(shè)在后面的靈敏度分析里可以觀察這個(gè)假設(shè)對(duì)結(jié)果的影響。2. 低碳調(diào)度的核心機(jī)制碳交易成本怎么進(jìn)目標(biāo)函數(shù)低碳調(diào)度和普通經(jīng)濟(jì)調(diào)度的本質(zhì)區(qū)別不在于約束方程多復(fù)雜而在于目標(biāo)函數(shù)里多了一項(xiàng)碳相關(guān)的成本并且因此改變了機(jī)組出力的優(yōu)先級(jí)排行。先理解碳交易機(jī)制的基本邏輯。在碳排放權(quán)交易體系下發(fā)電企業(yè)會(huì)被分配一定數(shù)量的碳排放配額。配額可以覆蓋全部或部分排放量實(shí)際排放超出配額的部分需要從市場(chǎng)上購(gòu)買碳排放權(quán)實(shí)際排放低于配額的部分則可以出售獲利相當(dāng)于給低碳機(jī)組一筆額外收入。2.1 碳配額怎么算目前常用的配額分配方法有歷史強(qiáng)度法按歷史排放強(qiáng)度分配和基準(zhǔn)線法按行業(yè)先進(jìn)水平分配。在電力系統(tǒng)調(diào)度模型中簡(jiǎn)化處理時(shí)一般用基準(zhǔn)線法配額 裝機(jī)容量 × 年等效利用小時(shí)數(shù) × 基準(zhǔn)排放強(qiáng)度但調(diào)度模型是日前或小時(shí)級(jí)問題所以需要把配額按時(shí)間尺度折算。常見做法是給每個(gè)機(jī)組設(shè)定一個(gè)單位電量碳排放配額比如0.6到0.8 tCO2/MWh實(shí)際單位排放強(qiáng)度高于這個(gè)值就得購(gòu)買碳配額低于這個(gè)值就能出售。這個(gè)值在不同國(guó)家和地區(qū)的碳市場(chǎng)中差異很大具體研究時(shí)需要按目標(biāo)市場(chǎng)的政策參數(shù)來設(shè)。2.2 碳交易成本的表達(dá)式假設(shè)機(jī)組i在t時(shí)刻出力為P(i,t)碳排放強(qiáng)度為E(i)tCO2/MWh碳配額強(qiáng)度為R(i)碳價(jià)為λ那么機(jī)組i的碳排放量為C(i,t) E(i) × P(i,t)配額為Q(i,t) R(i) × P(i,t)碳交易成本為T(i,t) λ × [C(i,t) - Q(i,t)] λ × [E(i) - R(i)] × P(i,t)從式子里能看出一個(gè)非常重要的現(xiàn)象如果某臺(tái)機(jī)組的碳排放強(qiáng)度高于配額強(qiáng)度碳交易成本是正的出力越多需要買的碳排放權(quán)越多反過來如果排放強(qiáng)度低于配額強(qiáng)度碳交易成本為負(fù)值相當(dāng)于給低碳機(jī)組一種綠色補(bǔ)貼出力越多越有利。這就改變了火電機(jī)組之間的經(jīng)濟(jì)競(jìng)爭(zhēng)關(guān)系——一臺(tái)低效但燃料便宜的老機(jī)組在碳價(jià)提高之后可能反而不如高效但燃料稍貴的機(jī)組有調(diào)度優(yōu)先級(jí)。2.3 目標(biāo)函數(shù)完整形式低碳調(diào)度的目標(biāo)函數(shù)可以寫為min F Σ(t) Σ(i) [ F_fuel(i,t) ST(i,t) ] Σ(t) Σ(i) [ λ × (E(i) - R(i)) × P(i,t) ]其中F_fuel(i,t)是燃料成本通常表達(dá)為二次函數(shù) F a bP cP2ST(i,t)是啟停成本由0-1變量控制。實(shí)際代碼里燃料成本的常數(shù)項(xiàng)a可以省略它對(duì)最優(yōu)解的出力分布沒有影響只影響目標(biāo)函數(shù)絕對(duì)值對(duì)場(chǎng)景對(duì)比和靈敏度分析也沒影響。約束條件與常規(guī)經(jīng)濟(jì)調(diào)度類似包括系統(tǒng)功率平衡約束所有機(jī)組出力加風(fēng)電出力等于負(fù)荷機(jī)組出力上下限約束機(jī)組爬坡速率約束最小啟停時(shí)間約束旋轉(zhuǎn)備用約束碳排放總量約束可選政策面設(shè)置系統(tǒng)碳預(yù)算時(shí)加入碳總量約束如果加上相當(dāng)于給整個(gè)系統(tǒng)設(shè)定了一個(gè)碳排放紅線模型會(huì)變得更嚴(yán)格也更容易看出碳價(jià)政策替代硬約束的力度差異。我在代碼里把這條做成可選項(xiàng)注釋里說明了開關(guān)位置。3. 不確定性建模為什么我選了場(chǎng)景法以及場(chǎng)景怎么來3.1 三種不確定性處理方法的工程化對(duì)比不確定性建模的常見手段有場(chǎng)景法隨機(jī)規(guī)劃、魯棒優(yōu)化和機(jī)會(huì)約束規(guī)劃三種各有適用場(chǎng)景方法數(shù)學(xué)模型輸出計(jì)劃優(yōu)點(diǎn)缺點(diǎn)場(chǎng)景法隨機(jī)規(guī)劃多場(chǎng)景期望最小化確定性出力計(jì)劃概念直觀經(jīng)濟(jì)指標(biāo)明確場(chǎng)景數(shù)量影響計(jì)算量魯棒優(yōu)化不確定集合內(nèi)最壞情況保守出力計(jì)劃安全性強(qiáng)集合參數(shù)易調(diào)結(jié)果偏保守經(jīng)濟(jì)性差機(jī)會(huì)約束規(guī)劃概率約束允許小概率越限出力計(jì)劃置信水平風(fēng)險(xiǎn)與經(jīng)濟(jì)平衡求解復(fù)雜分布假設(shè)敏感對(duì)于需要復(fù)現(xiàn)和做教學(xué)演示的Matlab代碼來說場(chǎng)景法是最合適的選擇。它把隨機(jī)性問題轉(zhuǎn)化成多個(gè)確定性問題的加權(quán)組合YALMIP配合CPLEX或者Gurobi可以直接求解不需要自己寫B(tài)enders分解或列約束生成這類復(fù)雜算法后續(xù)做靈敏度分析和結(jié)果可視化也很方便。3.2 源荷場(chǎng)景生成的具體做法我用的方法分三步。第一步生成風(fēng)電出力場(chǎng)景。設(shè)風(fēng)電預(yù)測(cè)出力為w_f(t)預(yù)測(cè)誤差用Beta分布描述每個(gè)場(chǎng)景的風(fēng)電出力為w_s(t) w_f(t) Δw(t)Δw按Beta分布采樣Beta分布的形狀參數(shù)根據(jù)歷史誤差統(tǒng)計(jì)擬合。注意Beta分布的支撐區(qū)間是[0,1]所以需要結(jié)合風(fēng)電裝機(jī)容量做量綱換算和上下界修正采樣后要強(qiáng)制截?cái)啾苊鈭?chǎng)景里出現(xiàn)負(fù)出力或超過裝機(jī)容量的不合理值。第二步生成負(fù)荷場(chǎng)景。負(fù)荷預(yù)測(cè)中心值為L(zhǎng)_f(t)誤差服從正態(tài)分布N(0, σ2_L(t))σ_L(t)取對(duì)應(yīng)時(shí)段歷史預(yù)測(cè)誤差的標(biāo)準(zhǔn)差。峰谷時(shí)段的σ_L可以做差異化處理高峰時(shí)段放大、低谷時(shí)段縮小更貼合實(shí)際情況。第三步場(chǎng)景削減。我首先生成1000到2000個(gè)初始場(chǎng)景然后用K-means聚類把場(chǎng)景數(shù)降到10到20個(gè)。K-means聚類質(zhì)心要按概率加權(quán)每個(gè)簇里包含的原場(chǎng)景數(shù)量作為權(quán)重這樣削減后的期望成本和原始場(chǎng)景集期望成本偏差會(huì)小很多。有些同行用同步回代消除法backward reduction實(shí)測(cè)效果也不錯(cuò)但K-means實(shí)現(xiàn)起來更簡(jiǎn)潔Matlab內(nèi)置了kmeans函數(shù)不需要額外寫概率距離計(jì)算的循環(huán)代碼量少很多跑起來也快。3.3 場(chǎng)景削減后的一個(gè)隱蔽坑K-means聚類完成后別忘記對(duì)削減后的場(chǎng)景做平滑化處理。因?yàn)榫垲愘|(zhì)心可能讓風(fēng)電出力在相鄰時(shí)段之間出現(xiàn)跳變這在功率平衡約束里沒有問題但在爬坡約束里可能導(dǎo)致本來不該觸發(fā)的爬坡被觸發(fā)了這種失真。我處理的辦法是對(duì)聚類后的場(chǎng)景序列做時(shí)間維度的中值濾波或者約束每個(gè)時(shí)段出力變化不超過預(yù)測(cè)場(chǎng)景的最大變化率。這個(gè)小細(xì)節(jié)直接關(guān)系到最后結(jié)果里爬坡約束是否真實(shí)有效。4. Matlab代碼實(shí)現(xiàn)從數(shù)學(xué)模型到可運(yùn)行代碼4.1 代碼結(jié)構(gòu)與數(shù)據(jù)準(zhǔn)備我的代碼文件組織如下main.m 主程序串聯(lián)場(chǎng)景生成、建模、求解、出圖 case30.m 系統(tǒng)數(shù)據(jù)機(jī)組參數(shù)、風(fēng)電場(chǎng)參數(shù)、負(fù)荷曲線 scenario_gen.m 源荷場(chǎng)景生成與削減 build_model.m YALMIP建模與求解 plot_result.m 結(jié)果繪圖與分析指標(biāo)計(jì)算case30.m里存放的主要數(shù)據(jù)包括機(jī)組數(shù)、爬坡率、出力上下限、燃料成本系數(shù)、碳排放強(qiáng)度、碳配額強(qiáng)度、風(fēng)電裝機(jī)容量、負(fù)荷時(shí)序數(shù)據(jù)。數(shù)據(jù)來源用的是IEEE 30節(jié)點(diǎn)系統(tǒng)的標(biāo)準(zhǔn)算例風(fēng)電場(chǎng)接在某個(gè)節(jié)點(diǎn)上替換掉一部分常規(guī)機(jī)組容量。設(shè)備參數(shù)表大致如下數(shù)值可按研究需要調(diào)整機(jī)組出力上限(MW)爬坡率(MW/h)燃料成本系數(shù)b碳排放強(qiáng)度(tCO2/MWh)碳配額強(qiáng)度G12005020.50.920.70G21504018.20.850.70G31003024.11.050.70風(fēng)電80-000風(fēng)電不產(chǎn)生碳排放燃料成本為零但引入不確定性后需要額外配置備用容量這個(gè)隱形代價(jià)會(huì)體現(xiàn)在約束和成本計(jì)算中。4.2 場(chǎng)景生成核心代碼場(chǎng)景生成部分的核心邏輯大致如下function [scenLoad, scenWind, prob] scenario_gen(load_f, wind_f, windCap, numClu) % load_f, wind_f: 預(yù)測(cè)序列長(zhǎng)度為T % windCap: 風(fēng)電場(chǎng)裝機(jī)容量 % numClu: 聚類后的場(chǎng)景數(shù) T length(load_f); nInit 1500; % 初始場(chǎng)景數(shù) alpha_w 2.0; beta_w 2.5; % Beta分布形狀參數(shù)按歷史誤差擬合 sigma_l 0.02 * load_f; % 負(fù)荷誤差標(biāo)準(zhǔn)差按比例設(shè)定 % 生成初始場(chǎng)景 initLoad zeros(nInit, T); initWind zeros(nInit, T); for k 1:nInit delta_w windCap * (betarnd(alpha_w, beta_w, 1, T) - 0.5) * 1.8; initWind(k,:) wind_f delta_w; initWind(k, initWind(k,:) 0) 0; initWind(k, initWind(k,:) windCap) windCap; initLoad(k,:) load_f normrnd(0, sigma_l); initLoad(k, initLoad(k,:) 0) 0; end % K-means聚類削減 X [initLoad initWind]; [idx, C] kmeans(X, numClu, Replicates, 10); prob accumarray(idx, 1) / nInit; % 各場(chǎng)景概率 scenLoad C(:, 1:T); scenWind C(:, T1:2*T); % 時(shí)間維平滑避免場(chǎng)景內(nèi)相鄰時(shí)段跳變過大 for s 1:numClu scenWind(s,:) medfilt1(scenWind(s,:), 3); scenLoad(s,:) medfilt1(scenLoad(s,:), 3); end end需要注意Beta分布的形狀參數(shù)直接決定了誤差分布是偏左還是偏右。如果場(chǎng)址的歷史數(shù)據(jù)顯示高估概率更大就把α調(diào)小β調(diào)大反之則調(diào)大α調(diào)小β。這個(gè)參數(shù)不是隨便拍腦袋的最好用實(shí)際風(fēng)電場(chǎng)的出力歷史數(shù)據(jù)做一次最大似然估計(jì)。4.3 YALMIP建模變量定義與約束實(shí)現(xiàn)YALMIP建模的關(guān)鍵在于把場(chǎng)景索引編進(jìn)變量里。機(jī)組出力P是三維變量機(jī)組數(shù)×?xí)r段數(shù)×場(chǎng)景數(shù)而0-1啟停變量u是二維的因?yàn)闄C(jī)組的開停機(jī)決策在所有場(chǎng)景下保持一致——這是場(chǎng)景法的標(biāo)準(zhǔn)處理方式叫非預(yù)期性約束。核心代碼示意%% 決策變量定義 P sdpvar(nG, T, numScen, full); % 機(jī)組出力 u binvar(nG, T, full); % 開停機(jī)狀態(tài)場(chǎng)景間一致 startup binvar(nG, T, full); % 啟動(dòng)標(biāo)志 shutdown binvar(nG, T, full); % 停機(jī)標(biāo)志 %% 約束集合 Constraints []; % 功率平衡每個(gè)場(chǎng)景下全部機(jī)組出力風(fēng)電負(fù)荷 for s 1:numScen for t 1:T Constraints [Constraints, sum(P(:,t,s), 1) scenWind(s,t) scenLoad(s,t)]; end end % 出力上下限與啟停耦合 for s 1:numScen for t 1:T Constraints [Constraints, P(:,t,s) u(:,t) .* Pmax]; Constraints [Constraints, P(:,t,s) u(:,t) .* Pmin]; end end % 爬坡約束 for s 1:numScen for t 2:T Constraints [Constraints, P(:,t,s) - P(:,t-1,s) ramp_up]; Constraints [Constraints, P(:,t,s) - P(:,t-1,s) -ramp_down]; end end % 旋轉(zhuǎn)備用約束常規(guī)機(jī)組的可調(diào)容量應(yīng)對(duì)風(fēng)電和負(fù)荷場(chǎng)景偏差 for t 1:T Constraints [Constraints, sum(P(:,t,1), 1) sum(scenWind(:,t)) ... scenLoad(:,t) * (1 reserve_rate)]; end %% 目標(biāo)函數(shù)場(chǎng)景期望成本 Objective 0; for s 1:numScen for t 1:T % 燃料成本 Objective Objective prob(s) * sum(fuel_b .* P(:,t,s) fuel_c .* P(:,t,s).^2); % 碳交易成本 Objective Objective prob(s) * sum(carbon_price .* (E_int - R_int) .* P(:,t,s)); % 啟停成本 Objective Objective prob(s) * sum(startup_cost .* startup(:,t)); end end %% 求解 ops sdpsettings(solver, cplex, verbose, 2, debug, 0); result optimize(Constraints, Objective, ops);這段代碼有幾個(gè)細(xì)節(jié)值得展開。一是啟停變量u所有場(chǎng)景共用這保證了日前調(diào)度給出的機(jī)組開停計(jì)劃是確定性的不會(huì)隨場(chǎng)景變化而變卦。如果讓u也按場(chǎng)景獨(dú)立實(shí)際上就變成看人下菜碟在實(shí)際運(yùn)行中做不到。二是爬坡約束需要特別注意時(shí)段邊界。t1時(shí)段沒有歷史出力我給的約束從t2開始避免邊界效應(yīng)。部分文獻(xiàn)會(huì)在初始時(shí)段假設(shè)機(jī)組以初始狀態(tài)連續(xù)出力那是另一種處理方式。三是目標(biāo)函數(shù)中的二次項(xiàng)fuel_c * P2在CPLEX里求解混合整數(shù)二次規(guī)劃MIQP是可以的但如果數(shù)據(jù)里的fuel_c很小數(shù)值上可能出現(xiàn)求解器報(bào)告數(shù)值問題的現(xiàn)象最好是先做歸一化處理或者把成本歸一到百萬(wàn)美元級(jí)別減少量級(jí)差距。四是場(chǎng)景概率的權(quán)重必須加對(duì)。每個(gè)場(chǎng)景的目標(biāo)函數(shù)貢獻(xiàn)要乘以prob(s)否則所有場(chǎng)景等權(quán)疊加聚類削減就白做了。4.4 求解器配置與常見坑YALMIP后端求解器我推薦CPLEX或Gurobi學(xué)術(shù)界用得最多。如果機(jī)器上沒裝也可以用MATLAB自帶的intlinprog求解但大規(guī)模場(chǎng)景下速度會(huì)明顯變慢。實(shí)測(cè)下來20個(gè)場(chǎng)景、24時(shí)段、6臺(tái)機(jī)組的混合整數(shù)二次規(guī)劃問題CPLEX大概幾秒到十幾秒能出結(jié)果intlinprog可能要幾分鐘。版本兼容性是一個(gè)容易踩的坑。Matlab版本和CPLEX版本必須匹配YALMIP官方文檔有兼容性表格安裝時(shí)建議先查一下。我自己的環(huán)境是Matlab 2023b配CPLEX 12.10跑常規(guī)70維左右的變量沒有問題。如果遇到無(wú)法識(shí)別cplex求解器的報(bào)錯(cuò)多半是路徑?jīng)]添加正確用addpath(genpath(C:\...\cplex\matlab))重新加載一遍即可。再提一個(gè)可能有用的技巧如果初始求解時(shí)發(fā)現(xiàn)某些二進(jìn)制變量無(wú)法收斂先檢查約束里有沒有等式兩端同時(shí)乘以二進(jìn)制變量這種非線性寫法YALMIP對(duì)這類問題會(huì)自動(dòng)引入大M法如果大M值設(shè)得過大或過小都會(huì)導(dǎo)致松弛問題病態(tài)表現(xiàn)為求解器輸出infeasible。5. 算例驗(yàn)證改裝后的IEEE 30節(jié)點(diǎn)系統(tǒng)能說明什么5.1 基礎(chǔ)數(shù)據(jù)與參數(shù)設(shè)置我用的算例是在IEEE 30節(jié)點(diǎn)系統(tǒng)基礎(chǔ)上改裝的。常規(guī)機(jī)組取3臺(tái)火電總裝機(jī)450MW風(fēng)電場(chǎng)裝機(jī)80MW負(fù)荷峰值230MW負(fù)荷曲線取典型冬季峰谷形狀。碳價(jià)格λ設(shè)置為25元/tCO2碳配額強(qiáng)度統(tǒng)一取0.70 tCO2/MWh旋轉(zhuǎn)備用率取5%。場(chǎng)景數(shù)量取15個(gè)初始場(chǎng)景1500個(gè)K-means削減后各場(chǎng)景概率分布相對(duì)均勻。對(duì)比方案是確定性調(diào)度不考慮任何不確定性用預(yù)測(cè)值直接求解和考慮源荷兩側(cè)不確定性的場(chǎng)景法調(diào)度。5.2 確定性調(diào)度和不確定性調(diào)度的對(duì)比先說總成本。確定性調(diào)度的期望成本比場(chǎng)景法調(diào)度低約4%到6%這個(gè)差異看起來不大但注意兩點(diǎn)一是確定性方案沒有預(yù)留足夠的爬坡和備用容量實(shí)際運(yùn)行中遭遇預(yù)測(cè)偏差時(shí)的二次調(diào)頻成本和外購(gòu)高價(jià)電力沒有計(jì)入門內(nèi)二是如果計(jì)算確定性方案在15個(gè)真實(shí)場(chǎng)景下的模擬運(yùn)行成本它會(huì)比它聲稱的成本高出一截而場(chǎng)景法方案在真實(shí)場(chǎng)景下的模擬成本與其優(yōu)化目標(biāo)基本一致。換句話說確定性方案省下來的是紙面成本不是真實(shí)成本。再看碳排放。場(chǎng)景法調(diào)度下的碳排放總量比確定性調(diào)度高約1%到2%原因在于為了應(yīng)對(duì)不確定性效率稍低但爬坡能力更強(qiáng)的機(jī)組會(huì)有更多的出力和啟停機(jī)會(huì)。這個(gè)時(shí)候如果碳價(jià)進(jìn)一步上升系統(tǒng)會(huì)趨向于讓高效機(jī)組多帶基荷、低效機(jī)組只做調(diào)節(jié)碳排放總量反而可能下降——這就是碳價(jià)對(duì)調(diào)度結(jié)構(gòu)的影響在結(jié)果上的最終體現(xiàn)。5.3 碳價(jià)與風(fēng)電比例的敏感性分析把碳價(jià)從20元/tCO2逐步提高到60元/tCO2觀察系統(tǒng)總碳排放和總成本的變化可以發(fā)現(xiàn)一個(gè)有趣的非線性現(xiàn)象碳價(jià)在20到40元區(qū)間時(shí)碳排放下降較慢因?yàn)榇藭r(shí)低效機(jī)組依然有一部分的基荷份額當(dāng)碳價(jià)超過40元后碳排放下降速度明顯加快說明臨界碳價(jià)促使低效機(jī)組的出力占比大幅壓縮。這個(gè)臨界碳價(jià)的位置取決于高效和低效機(jī)組之間的燃料成本差和排放強(qiáng)度差。再調(diào)整風(fēng)電滲透率從10%、20%到30%場(chǎng)景法調(diào)度的優(yōu)勢(shì)會(huì)越來越明顯。滲透率30%時(shí)確定性方案面對(duì)的不確定性已大到不加處理就很難保證安全約束的滿足程度而場(chǎng)景法方案雖然總成本上升但所有場(chǎng)景下均滿足功率平衡和爬坡約束。這正好印證了開篇的觀點(diǎn)高比例新能源接入之后不確定性不再是小修小補(bǔ)的問題而是調(diào)度建模的前提條件。6. 個(gè)人實(shí)操中的幾個(gè)體會(huì)代碼跑通只是第一步真正花時(shí)間的往往在數(shù)據(jù)清洗和結(jié)果調(diào)校上。我把自己反復(fù)踩過的一些教訓(xùn)列在這里供后來者參考。第一場(chǎng)景數(shù)量和求解時(shí)間的平衡點(diǎn)。一開始我以為場(chǎng)景越多越好結(jié)果30個(gè)場(chǎng)景的模型求解時(shí)間比15個(gè)場(chǎng)景翻了三倍不止而優(yōu)化結(jié)果的期望成本差異很小。實(shí)測(cè)下來15到20個(gè)場(chǎng)景對(duì)IEEE 30節(jié)點(diǎn)這樣的算例已經(jīng)足夠再往上的邊際收益微乎其微。如果系統(tǒng)規(guī)模更大可以先用少量場(chǎng)景快速驗(yàn)證模型邏輯等邏輯無(wú)誤后再加大場(chǎng)景數(shù)跑最終版。第二約束的可行域要留一點(diǎn)呼吸空間。有些文獻(xiàn)里功率平衡約束寫成嚴(yán)格等式實(shí)際求解時(shí)一旦場(chǎng)景生成存在數(shù)值誤差很容易出現(xiàn)infeasible。我的做法是在功率平衡約束中引入一個(gè)很小的松弛量比如0.01MW既不影響結(jié)果的合理性又避免了很多無(wú)謂的不可行問題。這個(gè)松弛在論文里可以注明在實(shí)際工程里本來就存在這樣量級(jí)的計(jì)量誤差。第三K-means聚類的初始化對(duì)結(jié)果有不可忽略的影響。我設(shè)置了Replicates, 10讓Matlab重復(fù)跑10次選最優(yōu)的聚類結(jié)果這樣能有效避免陷入局部最優(yōu)。如果時(shí)間緊張至少要設(shè)3次否則個(gè)別場(chǎng)景會(huì)明顯偏離原始分布特征。第四結(jié)果可視化不要只畫機(jī)組出力曲線要把風(fēng)電場(chǎng)景帶、負(fù)荷場(chǎng)景帶和系統(tǒng)備用容量放在同一張圖上看。場(chǎng)景帶直觀展示不確定性的量級(jí)調(diào)度結(jié)果落在場(chǎng)景帶內(nèi)才說明調(diào)度方案對(duì)不確定性免疫這一步對(duì)判斷模型的可靠性非常重要。最后如果后續(xù)想把模型擴(kuò)展得更精細(xì)可以考慮加入儲(chǔ)能系統(tǒng)、需求響應(yīng)資源或者碳捕集裝置。存儲(chǔ)的充放電策略天然適合用不確定性場(chǎng)景框架來優(yōu)化——正是這套場(chǎng)景法框架的擴(kuò)展場(chǎng)景儲(chǔ)能配置問題的建模邏輯和本文完全一致只是決策變量再多一個(gè)儲(chǔ)能SOC的狀態(tài)轉(zhuǎn)移約束而已。