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