據(jù)驅(qū)動分布魯棒優(yōu)化在電熱綜合能源系統(tǒng)調(diào)度中的Matlab實(shí)現(xiàn))
項(xiàng)目概述電熱綜合能源系統(tǒng)優(yōu)化本質(zhì)上是在一個(gè)同時(shí)包含電力網(wǎng)絡(luò)和熱力網(wǎng)絡(luò)的復(fù)雜系統(tǒng)里去解決“怎么調(diào)度設(shè)備、分配能量才能又省錢又可靠”的問題。這類系統(tǒng)的典型特征就是設(shè)備類型多——燃?xì)廨啓C(jī)、電鍋爐、儲熱罐、熱泵、余熱回收裝置等等而且電和熱之間還存在強(qiáng)耦合關(guān)系。最讓人頭疼的是系統(tǒng)運(yùn)行環(huán)境里的不確定性太多了風(fēng)電出力波動、光伏預(yù)測誤差、負(fù)荷變化這些都給調(diào)度決策帶來了很大的麻煩。傳統(tǒng)的做法無非是兩條路一條是隨機(jī)規(guī)劃假設(shè)不確定性參數(shù)服從某個(gè)已知概率分布然后算期望收益另一條是魯棒優(yōu)化干脆不考慮分布只守住不確定集合的最壞情況。前者的痛點(diǎn)在于真實(shí)場景下的概率分布很難準(zhǔn)確得知尤其是小樣本數(shù)據(jù)下估計(jì)出來的分布和真實(shí)分布差距可能非常大后者的痛點(diǎn)則在于過于保守為了覆蓋極端情況往往把系統(tǒng)運(yùn)行成本抬得很高實(shí)際經(jīng)濟(jì)性很差。分布魯棒優(yōu)化Distributionally Robust Optimization, DRO就是在這種背景下被越來越多研究者盯上的一個(gè)折中方案它不需要精確知道概率分布而是在一個(gè)“可能的分布集合”里找最壞情況下的最優(yōu)決策。如果再進(jìn)一步用數(shù)據(jù)驅(qū)動的方式去構(gòu)造這個(gè)“分布集合”——也就是模糊集Ambiguity Set就形成了標(biāo)題里提到的“數(shù)據(jù)驅(qū)動多離散場景分布魯棒”的技術(shù)路線。這篇博文我來把整個(gè)方案的思路拆解清楚從模型構(gòu)建到Matlab代碼實(shí)現(xiàn)從算法原理到實(shí)際跑代碼時(shí)踩過的坑事無巨細(xì)地分享出來。無論你是正在做綜合能源系統(tǒng)方向的研究生還是已經(jīng)入行做能源調(diào)度的工程師這篇文章都能幫你少走不少彎路。1. 先把問題說清楚電熱綜合能源系統(tǒng)優(yōu)化到底在優(yōu)化什么1.1 電熱耦合系統(tǒng)的核心設(shè)備與能量流在動手寫代碼之前第一步一定是把物理模型理清楚。電熱綜合能源系統(tǒng)不像單純的電力系統(tǒng)那樣只管有功無功它多了一張熱力網(wǎng)絡(luò)而且兩張網(wǎng)絡(luò)之間通過熱電聯(lián)產(chǎn)機(jī)組CHP、電鍋爐、熱泵這些耦合設(shè)備緊密聯(lián)系在一起。以我常用來做仿真的一個(gè)典型系統(tǒng)為例結(jié)構(gòu)大致是這樣的電源側(cè)外部電網(wǎng)可以買電、風(fēng)電機(jī)組出力不確定、燃?xì)廨啓C(jī)可控?zé)嵩磦?cè)CHP機(jī)組產(chǎn)電同時(shí)產(chǎn)熱、燃?xì)忮仩t純產(chǎn)熱、電鍋爐用電產(chǎn)熱儲能側(cè)電儲能電池、熱儲能儲熱罐負(fù)荷側(cè)電負(fù)荷、熱負(fù)荷。能量流的方向就是電網(wǎng)買電風(fēng)電CHP發(fā)電電池放電 → 供給電負(fù)荷CHP余熱燃?xì)忮仩t電鍋爐儲熱罐放熱 → 供給熱負(fù)荷。這里有個(gè)很有意思的耦合點(diǎn)電鍋爐和CHP把電和熱兩個(gè)系統(tǒng)聯(lián)系起來了你可以用電去產(chǎn)熱也可以讓CHP多發(fā)電順便產(chǎn)熱在調(diào)度上形成了很強(qiáng)的靈活性。建模的時(shí)候設(shè)備模型并不復(fù)雜。比如CHP機(jī)組通常用一個(gè)熱電比來約束它的電出力和熱出力之間的關(guān)系電出力范圍滿足上下限約束熱出力不大于熱電比乘以電出力 且熱出力本身也有上下限。再比如儲熱罐就是典型的狀態(tài)轉(zhuǎn)移方程儲熱量(下一時(shí)刻) 儲熱量(當(dāng)前時(shí)刻) × 散熱損失系數(shù) 充熱功率×效率 - 放熱功率/放熱效率這里有一個(gè)經(jīng)驗(yàn)性的提示很多初學(xué)的人會把熱網(wǎng)管道建模得很復(fù)雜加一堆溫度動態(tài)方程如果只是做日前調(diào)度層面的優(yōu)化其實(shí)沒必要把熱網(wǎng)簡化為節(jié)點(diǎn)熱功率平衡就夠了。過度精細(xì)化只會讓問題大得根本解不出來。1.2 不確定性從哪來為什么處理方式?jīng)Q定了方案質(zhì)量整個(gè)模型里最麻煩的是風(fēng)電出力和負(fù)荷預(yù)測誤差這些不確定量。你不能假設(shè)它們乖乖地等于預(yù)測值否則實(shí)際運(yùn)行的時(shí)候風(fēng)電突然少了電負(fù)荷就得切這在真實(shí)場景里是絕對不允許的。不確定性參數(shù)在模型里一般體現(xiàn)在機(jī)組出力約束、功率平衡約束里——說白了就是某些約束里帶的參數(shù)不是一個(gè)確定的數(shù)而是一個(gè)隨機(jī)量。這時(shí)候如何描述它就直接決定了你的優(yōu)化模型長什么樣也決定了求解難度和解的質(zhì)量。比如風(fēng)電出力為隨機(jī)變量那么功率平衡約束就得寫成電網(wǎng)購電 風(fēng)電出力 CHP電出力 電池放電 電負(fù)荷 電鍋爐耗電這個(gè)方程左側(cè)帶了一個(gè)無法精確預(yù)知的量。你要是用期望值替代等于告訴系統(tǒng)“風(fēng)電永遠(yuǎn)等于預(yù)測值”這在不確定性大的場景下是危險(xiǎn)的。你要是把所有可能值都考慮一遍又會導(dǎo)致決策過于保守。所以要在這里引入分布魯棒——它不假設(shè)風(fēng)電的分布是某個(gè)精確已知的函數(shù)而是給定一個(gè)包含真實(shí)分布的候選分布集合然后在這組候選分布里找最壞情況下的最優(yōu)決策。這個(gè)思路邏輯上確實(shí)比隨機(jī)規(guī)劃和傳統(tǒng)魯棒都要穩(wěn)。2. 為什么選“數(shù)據(jù)驅(qū)動分布魯棒”三種方法論的對比2.1 隨機(jī)規(guī)劃理想但不現(xiàn)實(shí)隨機(jī)規(guī)劃的思路是給每個(gè)不確定參數(shù)指定一個(gè)概率分布然后優(yōu)化目標(biāo)函數(shù)關(guān)于這個(gè)分布的期望值。比如風(fēng)電出力假設(shè)服從正態(tài)分布那么可以采樣生成大量場景每個(gè)場景帶一個(gè)概率權(quán)重構(gòu)建一個(gè)大規(guī)模的場景樹模型去求解。這個(gè)方法邏輯上沒問題前提是你對分布足夠了解。但問題恰恰出在這里——實(shí)際工程里風(fēng)電出力的分布形態(tài)往往是多峰、偏態(tài)的哪是簡單一個(gè)正態(tài)分布就能描述的你辛辛苦苦用歷史數(shù)據(jù)去擬合分布參數(shù)結(jié)果在小樣本情況下估計(jì)出來的分布和真實(shí)分布偏差很大算出來的方案自然也就不可靠。這也正是所謂“小樣本場景下數(shù)據(jù)驅(qū)動模型易過擬合”的典型體現(xiàn)——你把訓(xùn)練數(shù)據(jù)里的概率分布當(dāng)作真實(shí)分布去用了但數(shù)據(jù)少的時(shí)候這二者之間的差可能非常大。2.2 傳統(tǒng)魯棒優(yōu)化過于保守魯棒優(yōu)化的思路就簡單粗暴了不關(guān)心分布怎么樣只關(guān)心不確定量落在什么范圍內(nèi)。你給出一個(gè)不確定集合比如風(fēng)電出力在預(yù)測值上下20%浮動算法就在最壞的情況下做決策。好處是魯棒性強(qiáng)、計(jì)算簡單、不需要任何概率信息。壞處是它把集合里每個(gè)點(diǎn)都當(dāng)成等可能發(fā)生的事件來對待實(shí)際上有些極端情況發(fā)生的概率極低你為了這些極小概率事件讓方案變得非常保守——成本高到離譜甚至可能沒有可行解。這就好比為了防百年一遇的洪水把房子建在山頂上但代價(jià)是每天上下班都極其不方便。2.3 分布魯棒優(yōu)化站在兩者中間的平衡點(diǎn)分布魯棒優(yōu)化的思路是我不需要你告訴我精確分布但我會從歷史數(shù)據(jù)中構(gòu)造一組候選分布然后在這組分布里尋找使得系統(tǒng)運(yùn)行成本期望值最大的那個(gè)分布并針對它做出最優(yōu)決策。數(shù)學(xué)上可以寫成目標(biāo)函數(shù) min(第一階段的投資/調(diào)度成本 max_{分布∈模糊集} E[第二階段的運(yùn)行成本])看到這個(gè)兩層結(jié)構(gòu)沒有內(nèi)層是一個(gè)最大化問題在模糊集中找最壞分布外層是最小化問題在所有可能的分布下找一個(gè)綜合成本最低的調(diào)度方案。這個(gè)min-max結(jié)構(gòu)完美地結(jié)合了隨機(jī)規(guī)劃對分布信息的利用和魯棒優(yōu)化對不確定性的保守防護(hù)。而“數(shù)據(jù)驅(qū)動”在這里的角色是用歷史場景數(shù)據(jù)來構(gòu)造那個(gè)模糊集。說白了就是保證真實(shí)分布以較高的置信度落在這個(gè)集合里面。樣本越多集合越小方案越精確樣本越少集合越大方案越保守——但無論樣本多少都不至于讓你的方案因?yàn)榉植脊烙?jì)錯(cuò)誤而徹底失效。這個(gè)方法論的優(yōu)勢在風(fēng)電出力這類不確定性強(qiáng)的場景下體現(xiàn)得非常明顯。數(shù)據(jù)量足夠時(shí)方案幾乎可以和隨機(jī)規(guī)劃媲美數(shù)據(jù)量不足時(shí)也不至于像傳統(tǒng)魯棒那樣保守到?jīng)]法用。3. 核心機(jī)制拆解模糊集、場景生成與min-max求解策略3.1 數(shù)據(jù)驅(qū)動模糊集的構(gòu)造邏輯模糊集是整個(gè)分布魯棒優(yōu)化模型的心臟。它的作用是界定“哪些分布是可接受的”。最常用的一種構(gòu)造方式是基于矩的模糊集和基于Wasserstein距離的模糊集這篇博文重點(diǎn)講后者因?yàn)樵诙嚯x散場景的框架下Wasserstein距離的構(gòu)造更加自然而且有很好的理論性質(zhì)。基于Wasserstein距離的模糊集定義如下模糊集 { Q : Wasserstein距離(Q, 經(jīng)驗(yàn)分布) ≤ ε }什么意思呢就是說我們有一個(gè)由歷史數(shù)據(jù)得到的經(jīng)驗(yàn)分布所有和這個(gè)經(jīng)驗(yàn)分布的Wasserstein距離不超過半徑 ε 的分布都算在候選集合里。Wasserstein距離可以通俗地理解成“把一個(gè)概率分布搬運(yùn)成另一個(gè)概率分布的最小代價(jià)”它比KL散度之類的指標(biāo)更合適因?yàn)榫退銉蓚€(gè)分布的支撐集沒有重疊這個(gè)距離仍然是有限且有意義的。這里有個(gè)關(guān)鍵參數(shù) ε也就是模糊集半徑。它決定了你有多保守ε0時(shí)模糊集里只有經(jīng)驗(yàn)分布本身模型退化成了普通隨機(jī)規(guī)劃ε無窮大時(shí)模型退化成傳統(tǒng)魯棒優(yōu)化。實(shí)際中怎么選一般根據(jù)樣本數(shù)量、置信水平要求來定。樣本量越大ε可以取得越小。常見的一種做法是取經(jīng)驗(yàn)分布和真實(shí)分布之間的距離置信界也可以通過交叉驗(yàn)證來調(diào)參。我在實(shí)際代碼實(shí)現(xiàn)中通常會用這樣的公式來確定εε C / sqrt(N)其中N是場景數(shù)量C是一個(gè)和置信水平相關(guān)的常數(shù)。具體推導(dǎo)基于一個(gè)統(tǒng)計(jì)結(jié)論——經(jīng)驗(yàn)分布和真實(shí)分布的Wasserstein距離在概率意義下可以被上下界控制符合大數(shù)定律的收斂速率。這樣做的好處是你的模糊集大小不會拍腦袋拍出來而是有統(tǒng)計(jì)依據(jù)的。3.2 數(shù)據(jù)驅(qū)動多離散場景的生成與約減標(biāo)題里提到的“多離散場景”實(shí)際上就是把連續(xù)的不確定參數(shù)空間離散化為一系列帶有概率權(quán)重的典型場景。這步在工程實(shí)踐里是必須的因?yàn)橛?jì)算機(jī)沒法直接處理連續(xù)分布下的優(yōu)化問題但可以很輕松地處理“場景序號”這種離散變量。場景生成的流程我建議按以下步驟走收集原始數(shù)據(jù)風(fēng)電出力的歷史數(shù)據(jù)一般取過去1-2年的逐小時(shí)數(shù)據(jù)或者根據(jù)預(yù)測誤差的歷史統(tǒng)計(jì)來生成場景采樣如果已經(jīng)有了預(yù)測誤差的概率分布信息可以用蒙特卡洛采樣生成大量原始場景。采樣數(shù)量建議在1000-5000個(gè)左右先保證覆蓋面足夠廣場景約減用K-means聚類或者同步回代消除法Scenario Reduction把大量場景約減到幾十個(gè)有代表性的場景概率重分配每個(gè)聚類中心作為典型場景它包含的原始場景數(shù)量占總數(shù)的比例就是它的概率權(quán)重。我實(shí)際測試下來K-means聚類在大多數(shù)情況下都能用速度快、效果好。同步回代消除法在場景之間有很強(qiáng)相關(guān)性的情況下更合適但計(jì)算復(fù)雜度略高。做完這個(gè)步驟你得到的是一組場景集合形式大約是這樣場景1 (概率0.15): [風(fēng)電1, 風(fēng)電2, ..., 風(fēng)電24] 的24小時(shí)出力序列 場景2 (概率0.08): [風(fēng)電1, 風(fēng)電2, ..., 風(fēng)電24] 的24小時(shí)出力序列 ... 場景K (概率0.03): ...這些場景直接喂給分布魯棒模型做下一步的min-max優(yōu)化。這里有一個(gè)經(jīng)驗(yàn)值場景數(shù)量通常在10-30個(gè)之間就能在計(jì)算復(fù)雜度和解的精度之間取得不錯(cuò)的平衡。我曾經(jīng)試過用5個(gè)場景和30個(gè)場景分別做結(jié)果最優(yōu)成本只差了3%左右但計(jì)算時(shí)間卻差了一個(gè)數(shù)量級。3.3 兩階段分布魯棒模型的數(shù)學(xué)表達(dá)與求解策略現(xiàn)在把整個(gè)優(yōu)化模型完整地寫出來。兩階段分布魯棒優(yōu)化的標(biāo)準(zhǔn)形式是這樣的第一階段這里對應(yīng)日前調(diào)度決策: min ∑ { 第一階段成本(x) } max_{Q∈模糊集} E_Q[ 第二階段成本(y, ξ) ] 約束條件: 第一階段決策變量的運(yùn)行約束比如機(jī)組開停機(jī)、儲能初始狀態(tài)等 第二階段對應(yīng)實(shí)時(shí)調(diào)整決策依賴不確定參數(shù) ξ 的實(shí)現(xiàn): 給定 x 和 ξ 的實(shí)現(xiàn)值求解: min 第二階段成本(y) 約束條件: 功率平衡約束、設(shè)備出力上下限約束、儲能動態(tài)約束等取決于具體場景求解這個(gè)min-max問題主流的方法有兩大類一是對偶轉(zhuǎn)化把內(nèi)層最大化問題轉(zhuǎn)化為易處理的形式二是基于Benders分解或列與約束生成法CCG的迭代求解。在實(shí)際Matlab代碼中我最推薦CCG方法它比Benders分解收斂快得多。核心思路是主問題求解一個(gè)包含當(dāng)前已有場景的最優(yōu)調(diào)度問題給出決策x和最優(yōu)值下界子問題在給定x的情況下在所有場景里找最壞的那個(gè)場景及其對應(yīng)的成本并把結(jié)果反饋給主問題作為新的約束條件加進(jìn)去反復(fù)迭代直到上界和下界之間的gap小于設(shè)定閾值比如0.01%。這個(gè)流程一開始聽起來有點(diǎn)繞但寫代碼的時(shí)候很清晰。主問題是一個(gè)混合整數(shù)線性規(guī)劃因?yàn)槔锩嬗?-1變量比如機(jī)組啟停子問題是一個(gè)線性規(guī)劃。用Matlab調(diào)YalmipCplex/Gurobi幾分鐘就能搭出框架。4. Matlab代碼實(shí)現(xiàn)從建模到求解的完整實(shí)操流程4.1 求解前的環(huán)境配置與準(zhǔn)備工作Matlab環(huán)境下做這類問題最舒服的組合是Yalmip做建模層Cplex或Gurobi做底層求解器。Yalmip不是求解器它是一個(gè)建模工具箱能讓你用面向?qū)ο蟮姆绞綄懢€性規(guī)劃、混合整數(shù)規(guī)劃然后自動翻譯成求解器能吃的標(biāo)準(zhǔn)格式。關(guān)于工具箱如果你沒有Cplex用Gurobi也行兩者都支持MATLAB接口。如果連商業(yè)求解器都沒有先用免費(fèi)的Cbc或者GLPK頂著也行但求解混合整數(shù)規(guī)劃的速度會明顯慢很多大規(guī)模場景下不建議。安裝這里不詳細(xì)展開但有一條重要提示Yalmip和求解器版本的兼容性經(jīng)常出問題。我遇到過很多次代碼沒問題但結(jié)果不對最后發(fā)現(xiàn)是Cplex版本和Matlab版本不兼容導(dǎo)致的。建議使用Matlab R2021a以上的版本搭配Cplex 12.10或Gurobi 9.5以上這個(gè)組合比較穩(wěn)。4.2 場景數(shù)據(jù)生成模塊的代碼實(shí)現(xiàn)先給出場景生成部分的Matlab核心代碼框架。% 風(fēng)電出力場景生成和約減 % hist_data: 歷史風(fēng)電出力數(shù)據(jù), 維度為 N_history x T % N_scene: 需要保留的典型場景數(shù) rng(2025); % 固定隨機(jī)種子保證可復(fù)現(xiàn) N_history size(hist_data, 1); T size(hist_data, 2); % 時(shí)段數(shù)一般取24 % 步驟1: 蒙特卡洛采樣生成大量原始場景 % 以某時(shí)刻歷史數(shù)據(jù)的均值噪聲為例 mu mean(hist_data, 1); sigma std(hist_data, 1); N_sample 2000; scenarios_raw zeros(N_sample, T); for t 1:T % 用截?cái)嗾龖B(tài)分布防止出現(xiàn)負(fù)的風(fēng)電出力 pd makedist(Normal, mu, mu(t), sigma, sigma(t)); pd truncate(pd, 0, 1); scenarios_raw(:, t) random(pd, N_sample, 1); end % 步驟2: K-means聚類約減 [idx, centers] kmeans(scenarios_raw, N_scene); prob zeros(N_scene, 1); for k 1:N_scene prob(k) sum(idx k) / N_sample; end % 輸出: centers是典型場景矩庫(N_scene x T)prob是對應(yīng)概率 save(scenario_data.mat, centers, prob);這段代碼的邏輯很簡單生成大量樣本用K-means聚類找出幾個(gè)代表性中心點(diǎn)用每個(gè)簇的樣本比例作為概率權(quán)值。這里我特意用了截?cái)嗾龖B(tài)分布來防止負(fù)的風(fēng)電出力這是實(shí)際項(xiàng)目中很容易被忽略的細(xì)節(jié)——如果你不對隨機(jī)變量做截?cái)嗌沙鰜淼膱鼍翱赡芡耆环衔锢韺?shí)際。4.3 主問題與子問題的Yalmip建模實(shí)現(xiàn)接下來是核心部分兩階段分布魯棒模型的Matlab實(shí)現(xiàn)。由于完整代碼太長這里給出最關(guān)鍵的主問題和子問題結(jié)構(gòu)框架。先看主問題部分% 主問題: 調(diào)度決策 一個(gè)臨時(shí)變量eta表示最壞情況下的運(yùn)行成本 % x 是第一階段決策變量(機(jī)組出力、儲能充放電、購電等) % 需要定義u_cchp, p_chp, h_chp, u_gb, h_gb, p_eb, soc_es, soc_hs 等 ops sdpsettings(solver, cplex, verbose, 0); Constraints []; % 第一階段約束: 設(shè)備出力上下限、儲能動態(tài)、功率平衡期望場景下 Constraints [Constraints, 0 p_chp P_CHP_MAX]; Constraints [Constraints, 0 h_chp H_CHP_MAX]; Constraints [Constraints, h_chp R_CHP * p_chp]; % 熱電比約束 % ... 其他設(shè)備約束省略 % 目標(biāo)函數(shù): 第一階段成本 eta Objective sum(C_gas * (p_chp h_chp / R_CHP)) ... sum(C_buy .* p_grid) eta; % 迭代過程中不斷添加的CCG最優(yōu)割約束 for k 1:num_cuts Constraints [Constraints, eta sum(C_oper .* y_k) sum(Lagrange_mul_k .* (xi_k - x_expected))]; end optimize(Constraints, Objective, ops);這里面的關(guān)鍵是CCG思想的體現(xiàn)每迭代一次就會增加一個(gè)關(guān)于eta的割約束這個(gè)割約束里包含了來自子問題的最壞場景和拉格朗日乘子信息。隨著迭代進(jìn)行這些割約束逐漸逼近真實(shí)的最壞情況成本。再看子問題部分% 子問題: 給定主問題的決策 x_fixed在每個(gè)場景下求最優(yōu)運(yùn)行成本 % 然后選擇成本最高的場景作為最壞場景返回給主問題 costs zeros(N_scene, 1); for k 1:N_scene % 取當(dāng)前場景的風(fēng)電出力 wind centers(k, :); % 定義第二階段決策變量 y sdpvar(T, 1); % 棄風(fēng)量 shed sdpvar(T, 1); % 切負(fù)荷量 % 功率平衡約束 Constraints2 [p_grid - y - shed load_elec p_eb - wind - p_chp]; % ... 其他第二階段約束 % 目標(biāo)函數(shù): 棄風(fēng)懲罰 切負(fù)荷懲罰 Objective2 sum(C_curtail * y C_shed * shed); optimize(Constraints2, Objective2, ops); costs(k) value(Objective2); end % 找到最壞場景 [worst_cost, worst_idx] max(costs);子問題的本質(zhì)就是在每個(gè)離散場景下算一遍最優(yōu)運(yùn)行成本找到最大的那個(gè)——這就是“max”部分。注意這里的子問題我用了“棄風(fēng)和切負(fù)荷”這種松弛手段目的是一方面讓問題在極端場景下依然有可行解另一方面通過懲罰系數(shù)反映系統(tǒng)對不確定性的承受成本。切負(fù)荷懲罰系數(shù)通常設(shè)得很高比如1000元/MWh棄風(fēng)懲罰可以稍微低一點(diǎn)比如100元/MWh這兩個(gè)系數(shù)的設(shè)定直接影響調(diào)度策略的傾向性需要仔細(xì)權(quán)衡。4.4 CCG迭代求解的完整循環(huán)把主問題和子問題串起來就是完整的CCG迭代求解循環(huán)% 初始化 LB -inf; UB inf; gap 1; max_iter 100; tol 1e-4; while gap tol iter max_iter % 1. 求解主問題得到當(dāng)前最優(yōu)決策x和最優(yōu)值下界LB optimize(Constraints, Objective, ops); LB value(Objective); x_current value(x); % 2. 求解子問題得到最壞場景下運(yùn)行成本和上界UB [worst_cost, worst_idx] solve_subproblem(x_current); UB min(UB, first_stage_cost worst_cost); % 3. 將最壞場景生成的最優(yōu)割約束加入主問題 add_cut_to_master_problem(worst_idx, x_current); % 4. 更新迭代信息 gap abs((UB - LB) / UB); iter iter 1; end整個(gè)流程中還有一個(gè)實(shí)現(xiàn)細(xì)節(jié)值得注意主問題中如果也有二進(jìn)制變量比如機(jī)組啟停狀態(tài)那么主問題本身就是一個(gè)MILP問題子問題在求解時(shí)給定二進(jìn)制變量的值是已知的因此退化為一個(gè)LP問題。這種結(jié)構(gòu)下CCG方法依然能保證收斂而且收斂速度通常不錯(cuò)。我在測試中一般用24個(gè)時(shí)段、30個(gè)場景、再加4臺機(jī)組和2個(gè)儲能設(shè)備CCG迭代大約在10-25次內(nèi)就能收斂到0.01%的精度整個(gè)流程跑完以分鐘計(jì)。這個(gè)效率在論文復(fù)現(xiàn)和工程預(yù)算是完全夠用的。5. 避坑指南與常見問題排查5.1 求解極端緩慢收斂不了的典型案例我在跑代碼的時(shí)候如果說只遇到一個(gè)問題那就迭代特別慢、甚至震蕩。后來排查發(fā)現(xiàn)原因并不難找但往往藏得很隱蔽。第一個(gè)常見原因是場景數(shù)量太多。一開始我把場景數(shù)設(shè)成200個(gè)結(jié)果子問題每輪要解200次LP光這一步就非常慢。后來把場景約減到20個(gè)計(jì)算量直接降了一個(gè)數(shù)量級結(jié)果精度只損失了不超過2%。我的建議是先用10個(gè)場景跑通流程再逐步增加場景來觀察解的敏感性不要一開始就貪多。第二個(gè)常見原因是主問題是一個(gè)病態(tài)的MILP。比如機(jī)組啟停變量的Big-M約束中的M值取得太大會導(dǎo)致求解器數(shù)值穩(wěn)定性下降迭代效率嚴(yán)重降低。解決辦法是盡可能用小的合理M值或者用具有明確物理含義的約束來替換Big-M約束。第三個(gè)常見原因是子問題在某個(gè)特定場景下不可行。如果子問題無可行解那么整個(gè)CCG循環(huán)就會報(bào)錯(cuò)或進(jìn)入死循環(huán)。解決方式就是在子問題里加入松弛變量和懲罰項(xiàng)保證任何場景下都至少有一個(gè)可行解。這也正是我在4.3節(jié)里特意加入棄風(fēng)和切負(fù)荷松弛的另一個(gè)原因——它不僅僅是一個(gè)經(jīng)濟(jì)懲罰更是數(shù)學(xué)上保證算法穩(wěn)定性的保險(xiǎn)絲。5.2 結(jié)果不合常理怎么快速定位問題有時(shí)候算完了結(jié)果讓人一頭霧水。比如該買電的時(shí)候不買反而高價(jià)用氣發(fā)電或者儲能設(shè)備的行為完全反直覺。這些問題排查起來是有套路可循的。我先會去檢查約束條件是不是寫錯(cuò)了尤其是等式約束里的符號方向。Yalmip不報(bào)錯(cuò)不代表模型正確很多時(shí)候模型本身有問題但語法無誤照樣能給出一個(gè)“優(yōu)化結(jié)果”。第二個(gè)要排查的是參數(shù)的量綱一致性問題。電功率單位是MW熱功率單位可能誤寫成了kW比例系數(shù)差了1000倍結(jié)果必然千奇百怪。建議在代碼開頭集中定義所有參數(shù)并且統(tǒng)一單位比如全網(wǎng)都用MW和MWh絕不混用。第三我會把某個(gè)典型場景下各個(gè)設(shè)備的出力曲線全部畫出來疊加在同一個(gè)圖上。如果某個(gè)設(shè)備的出力長期頂在邊界上大概率是它的約束有問題如果某個(gè)設(shè)備完全沒有出力看看是不是啟動成本的懲罰系數(shù)太大導(dǎo)致模型寧愿不用它。5.3 參數(shù)敏感性模糊集半徑怎么確定才靠譜關(guān)于模糊集半徑ε的選取這是分布魯棒優(yōu)化里幾乎必被問的一個(gè)問題。我自己的經(jīng)驗(yàn)是分三步走第一步用理論公式計(jì)算初始值即前文提到的ε C/sqrt(N)第二步在這個(gè)初始值附近改變ε的值比如取0.1倍、0.5倍、1倍、2倍、5倍分別求解模型觀察最優(yōu)成本的變化曲線第三步選擇成本變化由陡變緩的拐點(diǎn)處對應(yīng)的ε作為最終取值。這個(gè)方法背后的邏輯是如果ε很小系統(tǒng)把不確定性看得過于樂觀成本低但風(fēng)險(xiǎn)大如果ε很大系統(tǒng)過于保守成本高但風(fēng)險(xiǎn)小。實(shí)際工程中你總能在中間找到一個(gè)合理的折中。關(guān)于模糊集半徑的設(shè)定有一個(gè)經(jīng)常被忽略的細(xì)節(jié)ε和場景數(shù)量N的匹配關(guān)系。理論上講N越多ε應(yīng)該越小。如果樣本量大卻選擇了一個(gè)很大的ε相當(dāng)于浪費(fèi)了數(shù)據(jù)信息如果樣本量小還選擇很小的ε模型就會變得過度自信失去分布魯棒的意義。這就是為什么很多論文中都會畫一張“ε vs 最優(yōu)成本”的敏感性分析圖目的就是為了驗(yàn)證參數(shù)選取得是否合理。模糊集半徑ε計(jì)算結(jié)果特征適用場景ε0退化為隨機(jī)規(guī)劃成本最低但忽視分布誤差歷史數(shù)據(jù)量極大且分布穩(wěn)定ε較小基于統(tǒng)計(jì)置信界成本適中兼顧穩(wěn)健性和經(jīng)濟(jì)性樣本數(shù)較多如500推薦ε中等人工調(diào)參結(jié)果成本偏高魯棒性較強(qiáng)樣本數(shù)一般如50-200常見選擇ε較大接近魯棒優(yōu)化成本高極端保守?cái)?shù)據(jù)極少或極端風(fēng)險(xiǎn)厭惡場景6. 從代碼到論文/項(xiàng)目落地結(jié)果驗(yàn)證與擴(kuò)展思考6.1 你的結(jié)果需要對比才更有說服力如果你是在做學(xué)術(shù)研究或者需要向團(tuán)隊(duì)證明這個(gè)方案的優(yōu)越性光有一個(gè)分布魯棒優(yōu)化的結(jié)果是不夠的你必須設(shè)置對照實(shí)驗(yàn)形成對比曲線和表格。正常情況下至少需要跑以下三組模型標(biāo)準(zhǔn)隨機(jī)規(guī)劃模型即ε0的退化情形傳統(tǒng)魯棒優(yōu)化模型不確定集合取各時(shí)刻風(fēng)電預(yù)測的上下界數(shù)據(jù)驅(qū)動分布魯棒模型即你實(shí)現(xiàn)的這個(gè)方案用不同場景數(shù)和模糊集半徑做多次實(shí)驗(yàn)。對比的指標(biāo)除了總運(yùn)行成本還要關(guān)注棄風(fēng)率、切負(fù)荷風(fēng)險(xiǎn)、以及不同分布偏差下的性能表現(xiàn)。一個(gè)常見做法是構(gòu)造一個(gè)“真實(shí)但未知”的分布讓三種方案的決策都在這個(gè)真實(shí)分布下做蒙特卡洛模擬測試看看誰的綜合表現(xiàn)最好。我在測試中典型的結(jié)果是分布魯棒優(yōu)化的總成本比隨機(jī)規(guī)劃只高3%-8%但切負(fù)荷風(fēng)險(xiǎn)大幅降低相比傳統(tǒng)魯棒優(yōu)化總成本能降低10%-20%且切負(fù)荷水平保持一致。這種結(jié)果圖一出來方案的價(jià)值一目了然。6.2 擴(kuò)展方向這份代碼還能怎么改如果做好了基礎(chǔ)版本還可以往幾個(gè)方向做擴(kuò)展一是把熱網(wǎng)動態(tài)特性加回來?；A(chǔ)版本里熱網(wǎng)被簡化為平衡約束如果加入管道傳輸延遲和熱損失模型決策會更精確但問題規(guī)模會顯著增大。二是引入置信區(qū)間自適應(yīng)調(diào)整。即模糊集半徑不是固定不變而是根據(jù)日前預(yù)測誤差的大小動態(tài)調(diào)整。比如天氣穩(wěn)定的日子半徑取小一點(diǎn)極端天氣時(shí)調(diào)大讓模型在不同場景下有不同保守程度。三是考慮多階段決策。日前調(diào)度是主問題日內(nèi)再通過模型預(yù)測控制MPC滾動修正。也就是說第一步的調(diào)度方案并不是一成不變執(zhí)行24小時(shí)而是每過一個(gè)小時(shí)就重新抬出最新的預(yù)測數(shù)據(jù)和場景修正方案。這種兩階段加滾動修正的組合方案在實(shí)際工程中應(yīng)用最廣也是我認(rèn)為這個(gè)方向最有落地前景的方向。6.3 常見問題速查表最后把我在整個(gè)開發(fā)過程中遇到的最典型、最有代表性的問題整理成一張速查表希望能幫你少走彎路。問題表現(xiàn)可能原因排查與解決思路模型求不出來提示不可行約束條件過強(qiáng)或互相矛盾加松弛變量和懲罰項(xiàng)檢查約束符號和量綱迭代收斂慢gap一直震蕩場景數(shù)過多或模糊集半徑過大適當(dāng)減少場景數(shù)檢查CCG割約束是否寫對結(jié)果異常設(shè)備出力全頂在上限出力范圍約束沒寫全或M值過小檢查設(shè)備上下限約束確認(rèn)Big-M值合理運(yùn)行時(shí)間太長主問題MILP規(guī)模太大嘗試固定啟停變量做松弛或減少整數(shù)變量維度不同隨機(jī)種子下結(jié)果差異大場景生成過程未固定種子設(shè)置rng固定種子并適當(dāng)增加場景采樣數(shù)Yalmip報(bào)錯(cuò)缺少求解器未安裝或未正確配置求解器運(yùn)行yalmiptest檢查求解器路徑和工作狀態(tài)風(fēng)電出力出現(xiàn)負(fù)值隨機(jī)采樣時(shí)未做截?cái)嚯S機(jī)數(shù)生成后用max(0, x)或截?cái)嗾龖B(tài)分布處理在代碼開發(fā)過程中我個(gè)人的習(xí)慣是每完成一個(gè)模塊就做一次短暫保存和結(jié)果輸出測試而不是等到代碼全寫完才整體調(diào)試。分布式魯棒優(yōu)化的代碼牽扯到主問題、子問題、場景生成和迭代邏輯四個(gè)大模塊相互之間耦合度很高一旦出現(xiàn)bug在全鏈路中排查會很痛苦。分段驗(yàn)證雖然多花了一點(diǎn)時(shí)間但發(fā)現(xiàn)問題時(shí)定位極快這個(gè)方法我用了很多年非常推薦。分布魯棒優(yōu)化的價(jià)值不只是發(fā)論文或者做一個(gè)好看的仿真結(jié)果。在電力市場改革深入、“雙碳”目標(biāo)持續(xù)推進(jìn)的背景下電熱綜合能源系統(tǒng)在需求側(cè)響應(yīng)、新能源消納這些實(shí)際工程場景里越來越常見。把小樣本數(shù)據(jù)下分布估計(jì)的不確定性考慮進(jìn)調(diào)度模型里讓方案既不盲目樂觀也不過分離譜這是從理論走向工程落地必須邁過的一關(guān)。這套Matlab代碼框架搭好之后后續(xù)不管是接入真實(shí)的風(fēng)電歷史數(shù)據(jù)還是擴(kuò)展成多能源品種的聯(lián)合調(diào)度都會快得多。希望這篇分享能幫你真正把算法運(yùn)行起來踩過的坑你都順利避開。 數(shù)據(jù)驅(qū)動分布魯棒優(yōu)化在電熱綜合能源系統(tǒng)調(diào)度中的Matlab實(shí)現(xiàn)
項(xiàng)目概述電熱綜合能源系統(tǒng)優(yōu)化本質(zhì)上是在一個(gè)同時(shí)包含電力網(wǎng)絡(luò)和熱力網(wǎng)絡(luò)的復(fù)雜系統(tǒng)里去解決“怎么調(diào)度設(shè)備、分配能量才能又省錢又可靠”的問題。這類系統(tǒng)的典型特征就是設(shè)備類型多——燃?xì)廨啓C(jī)、電鍋爐、儲熱罐、熱泵、余熱回收裝置等等而且電和熱之間還存在強(qiáng)耦合關(guān)系。最讓人頭疼的是系統(tǒng)運(yùn)行環(huán)境里的不確定性太多了風(fēng)電出力波動、光伏預(yù)測誤差、負(fù)荷變化這些都給調(diào)度決策帶來了很大的麻煩。傳統(tǒng)的做法無非是兩條路一條是隨機(jī)規(guī)劃假設(shè)不確定性參數(shù)服從某個(gè)已知概率分布然后算期望收益另一條是魯棒優(yōu)化干脆不考慮分布只守住不確定集合的最壞情況。前者的痛點(diǎn)在于真實(shí)場景下的概率分布很難準(zhǔn)確得知尤其是小樣本數(shù)據(jù)下估計(jì)出來的分布和真實(shí)分布差距可能非常大后者的痛點(diǎn)則在于過于保守為了覆蓋極端情況往往把系統(tǒng)運(yùn)行成本抬得很高實(shí)際經(jīng)濟(jì)性很差。分布魯棒優(yōu)化Distributionally Robust Optimization, DRO就是在這種背景下被越來越多研究者盯上的一個(gè)折中方案它不需要精確知道概率分布而是在一個(gè)“可能的分布集合”里找最壞情況下的最優(yōu)決策。如果再進(jìn)一步用數(shù)據(jù)驅(qū)動的方式去構(gòu)造這個(gè)“分布集合”——也就是模糊集Ambiguity Set就形成了標(biāo)題里提到的“數(shù)據(jù)驅(qū)動多離散場景分布魯棒”的技術(shù)路線。這篇博文我來把整個(gè)方案的思路拆解清楚從模型構(gòu)建到Matlab代碼實(shí)現(xiàn)從算法原理到實(shí)際跑代碼時(shí)踩過的坑事無巨細(xì)地分享出來。無論你是正在做綜合能源系統(tǒng)方向的研究生還是已經(jīng)入行做能源調(diào)度的工程師這篇文章都能幫你少走不少彎路。1. 先把問題說清楚電熱綜合能源系統(tǒng)優(yōu)化到底在優(yōu)化什么1.1 電熱耦合系統(tǒng)的核心設(shè)備與能量流在動手寫代碼之前第一步一定是把物理模型理清楚。電熱綜合能源系統(tǒng)不像單純的電力系統(tǒng)那樣只管有功無功它多了一張熱力網(wǎng)絡(luò)而且兩張網(wǎng)絡(luò)之間通過熱電聯(lián)產(chǎn)機(jī)組CHP、電鍋爐、熱泵這些耦合設(shè)備緊密聯(lián)系在一起。以我常用來做仿真的一個(gè)典型系統(tǒng)為例結(jié)構(gòu)大致是這樣的電源側(cè)外部電網(wǎng)可以買電、風(fēng)電機(jī)組出力不確定、燃?xì)廨啓C(jī)可控?zé)嵩磦?cè)CHP機(jī)組產(chǎn)電同時(shí)產(chǎn)熱、燃?xì)忮仩t純產(chǎn)熱、電鍋爐用電產(chǎn)熱儲能側(cè)電儲能電池、熱儲能儲熱罐負(fù)荷側(cè)電負(fù)荷、熱負(fù)荷。能量流的方向就是電網(wǎng)買電風(fēng)電CHP發(fā)電電池放電 → 供給電負(fù)荷CHP余熱燃?xì)忮仩t電鍋爐儲熱罐放熱 → 供給熱負(fù)荷。這里有個(gè)很有意思的耦合點(diǎn)電鍋爐和CHP把電和熱兩個(gè)系統(tǒng)聯(lián)系起來了你可以用電去產(chǎn)熱也可以讓CHP多發(fā)電順便產(chǎn)熱在調(diào)度上形成了很強(qiáng)的靈活性。建模的時(shí)候設(shè)備模型并不復(fù)雜。比如CHP機(jī)組通常用一個(gè)熱電比來約束它的電出力和熱出力之間的關(guān)系電出力范圍滿足上下限約束熱出力不大于熱電比乘以電出力 且熱出力本身也有上下限。再比如儲熱罐就是典型的狀態(tài)轉(zhuǎn)移方程儲熱量(下一時(shí)刻) 儲熱量(當(dāng)前時(shí)刻) × 散熱損失系數(shù) 充熱功率×效率 - 放熱功率/放熱效率這里有一個(gè)經(jīng)驗(yàn)性的提示很多初學(xué)的人會把熱網(wǎng)管道建模得很復(fù)雜加一堆溫度動態(tài)方程如果只是做日前調(diào)度層面的優(yōu)化其實(shí)沒必要把熱網(wǎng)簡化為節(jié)點(diǎn)熱功率平衡就夠了。過度精細(xì)化只會讓問題大得根本解不出來。1.2 不確定性從哪來為什么處理方式?jīng)Q定了方案質(zhì)量整個(gè)模型里最麻煩的是風(fēng)電出力和負(fù)荷預(yù)測誤差這些不確定量。你不能假設(shè)它們乖乖地等于預(yù)測值否則實(shí)際運(yùn)行的時(shí)候風(fēng)電突然少了電負(fù)荷就得切這在真實(shí)場景里是絕對不允許的。不確定性參數(shù)在模型里一般體現(xiàn)在機(jī)組出力約束、功率平衡約束里——說白了就是某些約束里帶的參數(shù)不是一個(gè)確定的數(shù)而是一個(gè)隨機(jī)量。這時(shí)候如何描述它就直接決定了你的優(yōu)化模型長什么樣也決定了求解難度和解的質(zhì)量。比如風(fēng)電出力為隨機(jī)變量那么功率平衡約束就得寫成電網(wǎng)購電 風(fēng)電出力 CHP電出力 電池放電 電負(fù)荷 電鍋爐耗電這個(gè)方程左側(cè)帶了一個(gè)無法精確預(yù)知的量。你要是用期望值替代等于告訴系統(tǒng)“風(fēng)電永遠(yuǎn)等于預(yù)測值”這在不確定性大的場景下是危險(xiǎn)的。你要是把所有可能值都考慮一遍又會導(dǎo)致決策過于保守。所以要在這里引入分布魯棒——它不假設(shè)風(fēng)電的分布是某個(gè)精確已知的函數(shù)而是給定一個(gè)包含真實(shí)分布的候選分布集合然后在這組候選分布里找最壞情況下的最優(yōu)決策。這個(gè)思路邏輯上確實(shí)比隨機(jī)規(guī)劃和傳統(tǒng)魯棒都要穩(wěn)。2. 為什么選“數(shù)據(jù)驅(qū)動分布魯棒”三種方法論的對比2.1 隨機(jī)規(guī)劃理想但不現(xiàn)實(shí)隨機(jī)規(guī)劃的思路是給每個(gè)不確定參數(shù)指定一個(gè)概率分布然后優(yōu)化目標(biāo)函數(shù)關(guān)于這個(gè)分布的期望值。比如風(fēng)電出力假設(shè)服從正態(tài)分布那么可以采樣生成大量場景每個(gè)場景帶一個(gè)概率權(quán)重構(gòu)建一個(gè)大規(guī)模的場景樹模型去求解。這個(gè)方法邏輯上沒問題前提是你對分布足夠了解。但問題恰恰出在這里——實(shí)際工程里風(fēng)電出力的分布形態(tài)往往是多峰、偏態(tài)的哪是簡單一個(gè)正態(tài)分布就能描述的你辛辛苦苦用歷史數(shù)據(jù)去擬合分布參數(shù)結(jié)果在小樣本情況下估計(jì)出來的分布和真實(shí)分布偏差很大算出來的方案自然也就不可靠。這也正是所謂“小樣本場景下數(shù)據(jù)驅(qū)動模型易過擬合”的典型體現(xiàn)——你把訓(xùn)練數(shù)據(jù)里的概率分布當(dāng)作真實(shí)分布去用了但數(shù)據(jù)少的時(shí)候這二者之間的差可能非常大。2.2 傳統(tǒng)魯棒優(yōu)化過于保守魯棒優(yōu)化的思路就簡單粗暴了不關(guān)心分布怎么樣只關(guān)心不確定量落在什么范圍內(nèi)。你給出一個(gè)不確定集合比如風(fēng)電出力在預(yù)測值上下20%浮動算法就在最壞的情況下做決策。好處是魯棒性強(qiáng)、計(jì)算簡單、不需要任何概率信息。壞處是它把集合里每個(gè)點(diǎn)都當(dāng)成等可能發(fā)生的事件來對待實(shí)際上有些極端情況發(fā)生的概率極低你為了這些極小概率事件讓方案變得非常保守——成本高到離譜甚至可能沒有可行解。這就好比為了防百年一遇的洪水把房子建在山頂上但代價(jià)是每天上下班都極其不方便。2.3 分布魯棒優(yōu)化站在兩者中間的平衡點(diǎn)分布魯棒優(yōu)化的思路是我不需要你告訴我精確分布但我會從歷史數(shù)據(jù)中構(gòu)造一組候選分布然后在這組分布里尋找使得系統(tǒng)運(yùn)行成本期望值最大的那個(gè)分布并針對它做出最優(yōu)決策。數(shù)學(xué)上可以寫成目標(biāo)函數(shù) min(第一階段的投資/調(diào)度成本 max_{分布∈模糊集} E[第二階段的運(yùn)行成本])看到這個(gè)兩層結(jié)構(gòu)沒有內(nèi)層是一個(gè)最大化問題在模糊集中找最壞分布外層是最小化問題在所有可能的分布下找一個(gè)綜合成本最低的調(diào)度方案。這個(gè)min-max結(jié)構(gòu)完美地結(jié)合了隨機(jī)規(guī)劃對分布信息的利用和魯棒優(yōu)化對不確定性的保守防護(hù)。而“數(shù)據(jù)驅(qū)動”在這里的角色是用歷史場景數(shù)據(jù)來構(gòu)造那個(gè)模糊集。說白了就是保證真實(shí)分布以較高的置信度落在這個(gè)集合里面。樣本越多集合越小方案越精確樣本越少集合越大方案越保守——但無論樣本多少都不至于讓你的方案因?yàn)榉植脊烙?jì)錯(cuò)誤而徹底失效。這個(gè)方法論的優(yōu)勢在風(fēng)電出力這類不確定性強(qiáng)的場景下體現(xiàn)得非常明顯。數(shù)據(jù)量足夠時(shí)方案幾乎可以和隨機(jī)規(guī)劃媲美數(shù)據(jù)量不足時(shí)也不至于像傳統(tǒng)魯棒那樣保守到?jīng)]法用。3. 核心機(jī)制拆解模糊集、場景生成與min-max求解策略3.1 數(shù)據(jù)驅(qū)動模糊集的構(gòu)造邏輯模糊集是整個(gè)分布魯棒優(yōu)化模型的心臟。它的作用是界定“哪些分布是可接受的”。最常用的一種構(gòu)造方式是基于矩的模糊集和基于Wasserstein距離的模糊集這篇博文重點(diǎn)講后者因?yàn)樵诙嚯x散場景的框架下Wasserstein距離的構(gòu)造更加自然而且有很好的理論性質(zhì)。基于Wasserstein距離的模糊集定義如下模糊集 { Q : Wasserstein距離(Q, 經(jīng)驗(yàn)分布) ≤ ε }什么意思呢就是說我們有一個(gè)由歷史數(shù)據(jù)得到的經(jīng)驗(yàn)分布所有和這個(gè)經(jīng)驗(yàn)分布的Wasserstein距離不超過半徑 ε 的分布都算在候選集合里。Wasserstein距離可以通俗地理解成“把一個(gè)概率分布搬運(yùn)成另一個(gè)概率分布的最小代價(jià)”它比KL散度之類的指標(biāo)更合適因?yàn)榫退銉蓚€(gè)分布的支撐集沒有重疊這個(gè)距離仍然是有限且有意義的。這里有個(gè)關(guān)鍵參數(shù) ε也就是模糊集半徑。它決定了你有多保守ε0時(shí)模糊集里只有經(jīng)驗(yàn)分布本身模型退化成了普通隨機(jī)規(guī)劃ε無窮大時(shí)模型退化成傳統(tǒng)魯棒優(yōu)化。實(shí)際中怎么選一般根據(jù)樣本數(shù)量、置信水平要求來定。樣本量越大ε可以取得越小。常見的一種做法是取經(jīng)驗(yàn)分布和真實(shí)分布之間的距離置信界也可以通過交叉驗(yàn)證來調(diào)參。我在實(shí)際代碼實(shí)現(xiàn)中通常會用這樣的公式來確定εε C / sqrt(N)其中N是場景數(shù)量C是一個(gè)和置信水平相關(guān)的常數(shù)。具體推導(dǎo)基于一個(gè)統(tǒng)計(jì)結(jié)論——經(jīng)驗(yàn)分布和真實(shí)分布的Wasserstein距離在概率意義下可以被上下界控制符合大數(shù)定律的收斂速率。這樣做的好處是你的模糊集大小不會拍腦袋拍出來而是有統(tǒng)計(jì)依據(jù)的。3.2 數(shù)據(jù)驅(qū)動多離散場景的生成與約減標(biāo)題里提到的“多離散場景”實(shí)際上就是把連續(xù)的不確定參數(shù)空間離散化為一系列帶有概率權(quán)重的典型場景。這步在工程實(shí)踐里是必須的因?yàn)橛?jì)算機(jī)沒法直接處理連續(xù)分布下的優(yōu)化問題但可以很輕松地處理“場景序號”這種離散變量。場景生成的流程我建議按以下步驟走收集原始數(shù)據(jù)風(fēng)電出力的歷史數(shù)據(jù)一般取過去1-2年的逐小時(shí)數(shù)據(jù)或者根據(jù)預(yù)測誤差的歷史統(tǒng)計(jì)來生成場景采樣如果已經(jīng)有了預(yù)測誤差的概率分布信息可以用蒙特卡洛采樣生成大量原始場景。采樣數(shù)量建議在1000-5000個(gè)左右先保證覆蓋面足夠廣場景約減用K-means聚類或者同步回代消除法Scenario Reduction把大量場景約減到幾十個(gè)有代表性的場景概率重分配每個(gè)聚類中心作為典型場景它包含的原始場景數(shù)量占總數(shù)的比例就是它的概率權(quán)重。我實(shí)際測試下來K-means聚類在大多數(shù)情況下都能用速度快、效果好。同步回代消除法在場景之間有很強(qiáng)相關(guān)性的情況下更合適但計(jì)算復(fù)雜度略高。做完這個(gè)步驟你得到的是一組場景集合形式大約是這樣場景1 (概率0.15): [風(fēng)電1, 風(fēng)電2, ..., 風(fēng)電24] 的24小時(shí)出力序列 場景2 (概率0.08): [風(fēng)電1, 風(fēng)電2, ..., 風(fēng)電24] 的24小時(shí)出力序列 ... 場景K (概率0.03): ...這些場景直接喂給分布魯棒模型做下一步的min-max優(yōu)化。這里有一個(gè)經(jīng)驗(yàn)值場景數(shù)量通常在10-30個(gè)之間就能在計(jì)算復(fù)雜度和解的精度之間取得不錯(cuò)的平衡。我曾經(jīng)試過用5個(gè)場景和30個(gè)場景分別做結(jié)果最優(yōu)成本只差了3%左右但計(jì)算時(shí)間卻差了一個(gè)數(shù)量級。3.3 兩階段分布魯棒模型的數(shù)學(xué)表達(dá)與求解策略現(xiàn)在把整個(gè)優(yōu)化模型完整地寫出來。兩階段分布魯棒優(yōu)化的標(biāo)準(zhǔn)形式是這樣的第一階段這里對應(yīng)日前調(diào)度決策: min ∑ { 第一階段成本(x) } max_{Q∈模糊集} E_Q[ 第二階段成本(y, ξ) ] 約束條件: 第一階段決策變量的運(yùn)行約束比如機(jī)組開停機(jī)、儲能初始狀態(tài)等 第二階段對應(yīng)實(shí)時(shí)調(diào)整決策依賴不確定參數(shù) ξ 的實(shí)現(xiàn): 給定 x 和 ξ 的實(shí)現(xiàn)值求解: min 第二階段成本(y) 約束條件: 功率平衡約束、設(shè)備出力上下限約束、儲能動態(tài)約束等取決于具體場景求解這個(gè)min-max問題主流的方法有兩大類一是對偶轉(zhuǎn)化把內(nèi)層最大化問題轉(zhuǎn)化為易處理的形式二是基于Benders分解或列與約束生成法CCG的迭代求解。在實(shí)際Matlab代碼中我最推薦CCG方法它比Benders分解收斂快得多。核心思路是主問題求解一個(gè)包含當(dāng)前已有場景的最優(yōu)調(diào)度問題給出決策x和最優(yōu)值下界子問題在給定x的情況下在所有場景里找最壞的那個(gè)場景及其對應(yīng)的成本并把結(jié)果反饋給主問題作為新的約束條件加進(jìn)去反復(fù)迭代直到上界和下界之間的gap小于設(shè)定閾值比如0.01%。這個(gè)流程一開始聽起來有點(diǎn)繞但寫代碼的時(shí)候很清晰。主問題是一個(gè)混合整數(shù)線性規(guī)劃因?yàn)槔锩嬗?-1變量比如機(jī)組啟停子問題是一個(gè)線性規(guī)劃。用Matlab調(diào)YalmipCplex/Gurobi幾分鐘就能搭出框架。4. Matlab代碼實(shí)現(xiàn)從建模到求解的完整實(shí)操流程4.1 求解前的環(huán)境配置與準(zhǔn)備工作Matlab環(huán)境下做這類問題最舒服的組合是Yalmip做建模層Cplex或Gurobi做底層求解器。Yalmip不是求解器它是一個(gè)建模工具箱能讓你用面向?qū)ο蟮姆绞綄懢€性規(guī)劃、混合整數(shù)規(guī)劃然后自動翻譯成求解器能吃的標(biāo)準(zhǔn)格式。關(guān)于工具箱如果你沒有Cplex用Gurobi也行兩者都支持MATLAB接口。如果連商業(yè)求解器都沒有先用免費(fèi)的Cbc或者GLPK頂著也行但求解混合整數(shù)規(guī)劃的速度會明顯慢很多大規(guī)模場景下不建議。安裝這里不詳細(xì)展開但有一條重要提示Yalmip和求解器版本的兼容性經(jīng)常出問題。我遇到過很多次代碼沒問題但結(jié)果不對最后發(fā)現(xiàn)是Cplex版本和Matlab版本不兼容導(dǎo)致的。建議使用Matlab R2021a以上的版本搭配Cplex 12.10或Gurobi 9.5以上這個(gè)組合比較穩(wěn)。4.2 場景數(shù)據(jù)生成模塊的代碼實(shí)現(xiàn)先給出場景生成部分的Matlab核心代碼框架。% 風(fēng)電出力場景生成和約減 % hist_data: 歷史風(fēng)電出力數(shù)據(jù), 維度為 N_history x T % N_scene: 需要保留的典型場景數(shù) rng(2025); % 固定隨機(jī)種子保證可復(fù)現(xiàn) N_history size(hist_data, 1); T size(hist_data, 2); % 時(shí)段數(shù)一般取24 % 步驟1: 蒙特卡洛采樣生成大量原始場景 % 以某時(shí)刻歷史數(shù)據(jù)的均值噪聲為例 mu mean(hist_data, 1); sigma std(hist_data, 1); N_sample 2000; scenarios_raw zeros(N_sample, T); for t 1:T % 用截?cái)嗾龖B(tài)分布防止出現(xiàn)負(fù)的風(fēng)電出力 pd makedist(Normal, mu, mu(t), sigma, sigma(t)); pd truncate(pd, 0, 1); scenarios_raw(:, t) random(pd, N_sample, 1); end % 步驟2: K-means聚類約減 [idx, centers] kmeans(scenarios_raw, N_scene); prob zeros(N_scene, 1); for k 1:N_scene prob(k) sum(idx k) / N_sample; end % 輸出: centers是典型場景矩庫(N_scene x T)prob是對應(yīng)概率 save(scenario_data.mat, centers, prob);這段代碼的邏輯很簡單生成大量樣本用K-means聚類找出幾個(gè)代表性中心點(diǎn)用每個(gè)簇的樣本比例作為概率權(quán)值。這里我特意用了截?cái)嗾龖B(tài)分布來防止負(fù)的風(fēng)電出力這是實(shí)際項(xiàng)目中很容易被忽略的細(xì)節(jié)——如果你不對隨機(jī)變量做截?cái)嗌沙鰜淼膱鼍翱赡芡耆环衔锢韺?shí)際。4.3 主問題與子問題的Yalmip建模實(shí)現(xiàn)接下來是核心部分兩階段分布魯棒模型的Matlab實(shí)現(xiàn)。由于完整代碼太長這里給出最關(guān)鍵的主問題和子問題結(jié)構(gòu)框架。先看主問題部分% 主問題: 調(diào)度決策 一個(gè)臨時(shí)變量eta表示最壞情況下的運(yùn)行成本 % x 是第一階段決策變量(機(jī)組出力、儲能充放電、購電等) % 需要定義u_cchp, p_chp, h_chp, u_gb, h_gb, p_eb, soc_es, soc_hs 等 ops sdpsettings(solver, cplex, verbose, 0); Constraints []; % 第一階段約束: 設(shè)備出力上下限、儲能動態(tài)、功率平衡期望場景下 Constraints [Constraints, 0 p_chp P_CHP_MAX]; Constraints [Constraints, 0 h_chp H_CHP_MAX]; Constraints [Constraints, h_chp R_CHP * p_chp]; % 熱電比約束 % ... 其他設(shè)備約束省略 % 目標(biāo)函數(shù): 第一階段成本 eta Objective sum(C_gas * (p_chp h_chp / R_CHP)) ... sum(C_buy .* p_grid) eta; % 迭代過程中不斷添加的CCG最優(yōu)割約束 for k 1:num_cuts Constraints [Constraints, eta sum(C_oper .* y_k) sum(Lagrange_mul_k .* (xi_k - x_expected))]; end optimize(Constraints, Objective, ops);這里面的關(guān)鍵是CCG思想的體現(xiàn)每迭代一次就會增加一個(gè)關(guān)于eta的割約束這個(gè)割約束里包含了來自子問題的最壞場景和拉格朗日乘子信息。隨著迭代進(jìn)行這些割約束逐漸逼近真實(shí)的最壞情況成本。再看子問題部分% 子問題: 給定主問題的決策 x_fixed在每個(gè)場景下求最優(yōu)運(yùn)行成本 % 然后選擇成本最高的場景作為最壞場景返回給主問題 costs zeros(N_scene, 1); for k 1:N_scene % 取當(dāng)前場景的風(fēng)電出力 wind centers(k, :); % 定義第二階段決策變量 y sdpvar(T, 1); % 棄風(fēng)量 shed sdpvar(T, 1); % 切負(fù)荷量 % 功率平衡約束 Constraints2 [p_grid - y - shed load_elec p_eb - wind - p_chp]; % ... 其他第二階段約束 % 目標(biāo)函數(shù): 棄風(fēng)懲罰 切負(fù)荷懲罰 Objective2 sum(C_curtail * y C_shed * shed); optimize(Constraints2, Objective2, ops); costs(k) value(Objective2); end % 找到最壞場景 [worst_cost, worst_idx] max(costs);子問題的本質(zhì)就是在每個(gè)離散場景下算一遍最優(yōu)運(yùn)行成本找到最大的那個(gè)——這就是“max”部分。注意這里的子問題我用了“棄風(fēng)和切負(fù)荷”這種松弛手段目的是一方面讓問題在極端場景下依然有可行解另一方面通過懲罰系數(shù)反映系統(tǒng)對不確定性的承受成本。切負(fù)荷懲罰系數(shù)通常設(shè)得很高比如1000元/MWh棄風(fēng)懲罰可以稍微低一點(diǎn)比如100元/MWh這兩個(gè)系數(shù)的設(shè)定直接影響調(diào)度策略的傾向性需要仔細(xì)權(quán)衡。4.4 CCG迭代求解的完整循環(huán)把主問題和子問題串起來就是完整的CCG迭代求解循環(huán)% 初始化 LB -inf; UB inf; gap 1; max_iter 100; tol 1e-4; while gap tol iter max_iter % 1. 求解主問題得到當(dāng)前最優(yōu)決策x和最優(yōu)值下界LB optimize(Constraints, Objective, ops); LB value(Objective); x_current value(x); % 2. 求解子問題得到最壞場景下運(yùn)行成本和上界UB [worst_cost, worst_idx] solve_subproblem(x_current); UB min(UB, first_stage_cost worst_cost); % 3. 將最壞場景生成的最優(yōu)割約束加入主問題 add_cut_to_master_problem(worst_idx, x_current); % 4. 更新迭代信息 gap abs((UB - LB) / UB); iter iter 1; end整個(gè)流程中還有一個(gè)實(shí)現(xiàn)細(xì)節(jié)值得注意主問題中如果也有二進(jìn)制變量比如機(jī)組啟停狀態(tài)那么主問題本身就是一個(gè)MILP問題子問題在求解時(shí)給定二進(jìn)制變量的值是已知的因此退化為一個(gè)LP問題。這種結(jié)構(gòu)下CCG方法依然能保證收斂而且收斂速度通常不錯(cuò)。我在測試中一般用24個(gè)時(shí)段、30個(gè)場景、再加4臺機(jī)組和2個(gè)儲能設(shè)備CCG迭代大約在10-25次內(nèi)就能收斂到0.01%的精度整個(gè)流程跑完以分鐘計(jì)。這個(gè)效率在論文復(fù)現(xiàn)和工程預(yù)算是完全夠用的。5. 避坑指南與常見問題排查5.1 求解極端緩慢收斂不了的典型案例我在跑代碼的時(shí)候如果說只遇到一個(gè)問題那就迭代特別慢、甚至震蕩。后來排查發(fā)現(xiàn)原因并不難找但往往藏得很隱蔽。第一個(gè)常見原因是場景數(shù)量太多。一開始我把場景數(shù)設(shè)成200個(gè)結(jié)果子問題每輪要解200次LP光這一步就非常慢。后來把場景約減到20個(gè)計(jì)算量直接降了一個(gè)數(shù)量級結(jié)果精度只損失了不超過2%。我的建議是先用10個(gè)場景跑通流程再逐步增加場景來觀察解的敏感性不要一開始就貪多。第二個(gè)常見原因是主問題是一個(gè)病態(tài)的MILP。比如機(jī)組啟停變量的Big-M約束中的M值取得太大會導(dǎo)致求解器數(shù)值穩(wěn)定性下降迭代效率嚴(yán)重降低。解決辦法是盡可能用小的合理M值或者用具有明確物理含義的約束來替換Big-M約束。第三個(gè)常見原因是子問題在某個(gè)特定場景下不可行。如果子問題無可行解那么整個(gè)CCG循環(huán)就會報(bào)錯(cuò)或進(jìn)入死循環(huán)。解決方式就是在子問題里加入松弛變量和懲罰項(xiàng)保證任何場景下都至少有一個(gè)可行解。這也正是我在4.3節(jié)里特意加入棄風(fēng)和切負(fù)荷松弛的另一個(gè)原因——它不僅僅是一個(gè)經(jīng)濟(jì)懲罰更是數(shù)學(xué)上保證算法穩(wěn)定性的保險(xiǎn)絲。5.2 結(jié)果不合常理怎么快速定位問題有時(shí)候算完了結(jié)果讓人一頭霧水。比如該買電的時(shí)候不買反而高價(jià)用氣發(fā)電或者儲能設(shè)備的行為完全反直覺。這些問題排查起來是有套路可循的。我先會去檢查約束條件是不是寫錯(cuò)了尤其是等式約束里的符號方向。Yalmip不報(bào)錯(cuò)不代表模型正確很多時(shí)候模型本身有問題但語法無誤照樣能給出一個(gè)“優(yōu)化結(jié)果”。第二個(gè)要排查的是參數(shù)的量綱一致性問題。電功率單位是MW熱功率單位可能誤寫成了kW比例系數(shù)差了1000倍結(jié)果必然千奇百怪。建議在代碼開頭集中定義所有參數(shù)并且統(tǒng)一單位比如全網(wǎng)都用MW和MWh絕不混用。第三我會把某個(gè)典型場景下各個(gè)設(shè)備的出力曲線全部畫出來疊加在同一個(gè)圖上。如果某個(gè)設(shè)備的出力長期頂在邊界上大概率是它的約束有問題如果某個(gè)設(shè)備完全沒有出力看看是不是啟動成本的懲罰系數(shù)太大導(dǎo)致模型寧愿不用它。5.3 參數(shù)敏感性模糊集半徑怎么確定才靠譜關(guān)于模糊集半徑ε的選取這是分布魯棒優(yōu)化里幾乎必被問的一個(gè)問題。我自己的經(jīng)驗(yàn)是分三步走第一步用理論公式計(jì)算初始值即前文提到的ε C/sqrt(N)第二步在這個(gè)初始值附近改變ε的值比如取0.1倍、0.5倍、1倍、2倍、5倍分別求解模型觀察最優(yōu)成本的變化曲線第三步選擇成本變化由陡變緩的拐點(diǎn)處對應(yīng)的ε作為最終取值。這個(gè)方法背后的邏輯是如果ε很小系統(tǒng)把不確定性看得過于樂觀成本低但風(fēng)險(xiǎn)大如果ε很大系統(tǒng)過于保守成本高但風(fēng)險(xiǎn)小。實(shí)際工程中你總能在中間找到一個(gè)合理的折中。關(guān)于模糊集半徑的設(shè)定有一個(gè)經(jīng)常被忽略的細(xì)節(jié)ε和場景數(shù)量N的匹配關(guān)系。理論上講N越多ε應(yīng)該越小。如果樣本量大卻選擇了一個(gè)很大的ε相當(dāng)于浪費(fèi)了數(shù)據(jù)信息如果樣本量小還選擇很小的ε模型就會變得過度自信失去分布魯棒的意義。這就是為什么很多論文中都會畫一張“ε vs 最優(yōu)成本”的敏感性分析圖目的就是為了驗(yàn)證參數(shù)選取得是否合理。模糊集半徑ε計(jì)算結(jié)果特征適用場景ε0退化為隨機(jī)規(guī)劃成本最低但忽視分布誤差歷史數(shù)據(jù)量極大且分布穩(wěn)定ε較小基于統(tǒng)計(jì)置信界成本適中兼顧穩(wěn)健性和經(jīng)濟(jì)性樣本數(shù)較多如500推薦ε中等人工調(diào)參結(jié)果成本偏高魯棒性較強(qiáng)樣本數(shù)一般如50-200常見選擇ε較大接近魯棒優(yōu)化成本高極端保守?cái)?shù)據(jù)極少或極端風(fēng)險(xiǎn)厭惡場景6. 從代碼到論文/項(xiàng)目落地結(jié)果驗(yàn)證與擴(kuò)展思考6.1 你的結(jié)果需要對比才更有說服力如果你是在做學(xué)術(shù)研究或者需要向團(tuán)隊(duì)證明這個(gè)方案的優(yōu)越性光有一個(gè)分布魯棒優(yōu)化的結(jié)果是不夠的你必須設(shè)置對照實(shí)驗(yàn)形成對比曲線和表格。正常情況下至少需要跑以下三組模型標(biāo)準(zhǔn)隨機(jī)規(guī)劃模型即ε0的退化情形傳統(tǒng)魯棒優(yōu)化模型不確定集合取各時(shí)刻風(fēng)電預(yù)測的上下界數(shù)據(jù)驅(qū)動分布魯棒模型即你實(shí)現(xiàn)的這個(gè)方案用不同場景數(shù)和模糊集半徑做多次實(shí)驗(yàn)。對比的指標(biāo)除了總運(yùn)行成本還要關(guān)注棄風(fēng)率、切負(fù)荷風(fēng)險(xiǎn)、以及不同分布偏差下的性能表現(xiàn)。一個(gè)常見做法是構(gòu)造一個(gè)“真實(shí)但未知”的分布讓三種方案的決策都在這個(gè)真實(shí)分布下做蒙特卡洛模擬測試看看誰的綜合表現(xiàn)最好。我在測試中典型的結(jié)果是分布魯棒優(yōu)化的總成本比隨機(jī)規(guī)劃只高3%-8%但切負(fù)荷風(fēng)險(xiǎn)大幅降低相比傳統(tǒng)魯棒優(yōu)化總成本能降低10%-20%且切負(fù)荷水平保持一致。這種結(jié)果圖一出來方案的價(jià)值一目了然。6.2 擴(kuò)展方向這份代碼還能怎么改如果做好了基礎(chǔ)版本還可以往幾個(gè)方向做擴(kuò)展一是把熱網(wǎng)動態(tài)特性加回來?；A(chǔ)版本里熱網(wǎng)被簡化為平衡約束如果加入管道傳輸延遲和熱損失模型決策會更精確但問題規(guī)模會顯著增大。二是引入置信區(qū)間自適應(yīng)調(diào)整。即模糊集半徑不是固定不變而是根據(jù)日前預(yù)測誤差的大小動態(tài)調(diào)整。比如天氣穩(wěn)定的日子半徑取小一點(diǎn)極端天氣時(shí)調(diào)大讓模型在不同場景下有不同保守程度。三是考慮多階段決策。日前調(diào)度是主問題日內(nèi)再通過模型預(yù)測控制MPC滾動修正。也就是說第一步的調(diào)度方案并不是一成不變執(zhí)行24小時(shí)而是每過一個(gè)小時(shí)就重新抬出最新的預(yù)測數(shù)據(jù)和場景修正方案。這種兩階段加滾動修正的組合方案在實(shí)際工程中應(yīng)用最廣也是我認(rèn)為這個(gè)方向最有落地前景的方向。6.3 常見問題速查表最后把我在整個(gè)開發(fā)過程中遇到的最典型、最有代表性的問題整理成一張速查表希望能幫你少走彎路。問題表現(xiàn)可能原因排查與解決思路模型求不出來提示不可行約束條件過強(qiáng)或互相矛盾加松弛變量和懲罰項(xiàng)檢查約束符號和量綱迭代收斂慢gap一直震蕩場景數(shù)過多或模糊集半徑過大適當(dāng)減少場景數(shù)檢查CCG割約束是否寫對結(jié)果異常設(shè)備出力全頂在上限出力范圍約束沒寫全或M值過小檢查設(shè)備上下限約束確認(rèn)Big-M值合理運(yùn)行時(shí)間太長主問題MILP規(guī)模太大嘗試固定啟停變量做松弛或減少整數(shù)變量維度不同隨機(jī)種子下結(jié)果差異大場景生成過程未固定種子設(shè)置rng固定種子并適當(dāng)增加場景采樣數(shù)Yalmip報(bào)錯(cuò)缺少求解器未安裝或未正確配置求解器運(yùn)行yalmiptest檢查求解器路徑和工作狀態(tài)風(fēng)電出力出現(xiàn)負(fù)值隨機(jī)采樣時(shí)未做截?cái)嚯S機(jī)數(shù)生成后用max(0, x)或截?cái)嗾龖B(tài)分布處理在代碼開發(fā)過程中我個(gè)人的習(xí)慣是每完成一個(gè)模塊就做一次短暫保存和結(jié)果輸出測試而不是等到代碼全寫完才整體調(diào)試。分布式魯棒優(yōu)化的代碼牽扯到主問題、子問題、場景生成和迭代邏輯四個(gè)大模塊相互之間耦合度很高一旦出現(xiàn)bug在全鏈路中排查會很痛苦。分段驗(yàn)證雖然多花了一點(diǎn)時(shí)間但發(fā)現(xiàn)問題時(shí)定位極快這個(gè)方法我用了很多年非常推薦。分布魯棒優(yōu)化的價(jià)值不只是發(fā)論文或者做一個(gè)好看的仿真結(jié)果。在電力市場改革深入、“雙碳”目標(biāo)持續(xù)推進(jìn)的背景下電熱綜合能源系統(tǒng)在需求側(cè)響應(yīng)、新能源消納這些實(shí)際工程場景里越來越常見。把小樣本數(shù)據(jù)下分布估計(jì)的不確定性考慮進(jìn)調(diào)度模型里讓方案既不盲目樂觀也不過分離譜這是從理論走向工程落地必須邁過的一關(guān)。這套Matlab代碼框架搭好之后后續(xù)不管是接入真實(shí)的風(fēng)電歷史數(shù)據(jù)還是擴(kuò)展成多能源品種的聯(lián)合調(diào)度都會快得多。希望這篇分享能幫你真正把算法運(yùn)行起來踩過的坑你都順利避開。