做EI論文的代碼復(fù)現(xiàn)最怕的不是數(shù)學(xué)看不懂而是看不懂的地方恰好卡在工程實(shí)現(xiàn)上。今天這篇我想用實(shí)際做過(guò)的一個(gè)項(xiàng)目——“考慮靈活性的數(shù)據(jù)中心微網(wǎng)兩階段魯棒規(guī)劃方法”——來(lái)完整走一遍復(fù)現(xiàn)流程。這個(gè)方向在微網(wǎng)規(guī)劃里屬于偏應(yīng)用又偏方法的交叉點(diǎn)Matlab實(shí)現(xiàn)過(guò)程中牽扯到的不只是調(diào)用求解器還有對(duì)偶變換、不確定性集合構(gòu)造、以及數(shù)據(jù)中心這種特殊負(fù)荷的建模方式。文章會(huì)盡量把每個(gè)決策點(diǎn)背后的“為什么”講清楚適合正在做微網(wǎng)規(guī)劃、魯棒優(yōu)化方向或者接了類(lèi)似EI復(fù)現(xiàn)任務(wù)的讀者參考。我默認(rèn)你手里已經(jīng)有一篇論文原文最好是帶數(shù)學(xué)模型的。如果沒(méi)有先別急著看代碼下文會(huì)先從拆解論文開(kāi)始把復(fù)現(xiàn)思路理順再談Matlab里怎么落地。1. 復(fù)現(xiàn)前必須想清楚的三個(gè)問(wèn)題1.1 這篇EI論文到底解決了什么問(wèn)題數(shù)據(jù)中心微網(wǎng)規(guī)劃和普通園區(qū)微網(wǎng)規(guī)劃看著相似實(shí)際差別很大。普通微網(wǎng)考慮負(fù)荷波動(dòng)時(shí)常用的做法是預(yù)測(cè)一個(gè)典型日負(fù)荷曲線然后用儲(chǔ)能、分布式電源去跟。但數(shù)據(jù)中心有兩個(gè)特殊之處一是負(fù)荷里有一部分IT負(fù)載是可以在短時(shí)間內(nèi)彈性調(diào)整的二是數(shù)據(jù)中心對(duì)供電連續(xù)性極其敏感后端柴油發(fā)電機(jī)和儲(chǔ)能系統(tǒng)不只是經(jīng)濟(jì)性的選擇更是供配電架構(gòu)的硬需求。論文里引入“靈活性”這個(gè)概念本質(zhì)上是在傳統(tǒng)以“最小化投資和運(yùn)行成本”為目標(biāo)的兩階段魯棒規(guī)劃框架中額外刻畫(huà)了系統(tǒng)對(duì)不確定性的調(diào)節(jié)能力。換句話說(shuō)這篇論文解決的不是“怎么規(guī)劃一個(gè)數(shù)據(jù)中心微網(wǎng)”而是“規(guī)劃出來(lái)的微網(wǎng)在極端天氣、負(fù)荷突增這些場(chǎng)景下能不能靠自身調(diào)節(jié)手段扛住”并且把這種扛住的能力量化進(jìn)了規(guī)劃模型里。復(fù)現(xiàn)的時(shí)候如果不能把這個(gè)物理邏輯映射到數(shù)學(xué)約束上后面做出來(lái)的模型很可能會(huì)偏離原文的核心。1.2 兩階段魯棒規(guī)劃為什么適合數(shù)據(jù)中心微網(wǎng)兩階段魯棒的核心思想是“先做決定再優(yōu)化調(diào)整”。規(guī)劃階段第一階段確定設(shè)備裝多少、裝什么類(lèi)型這些決策在不確定性真正到來(lái)之前就要敲定運(yùn)行階段第二階段是在不確定性實(shí)現(xiàn)之后通過(guò)調(diào)度手段盡可能低成本地平衡供需。數(shù)據(jù)中心微網(wǎng)的棘手之處就在于不確定性源很多光伏出力、風(fēng)電出力、實(shí)時(shí)電價(jià)、IT負(fù)載突發(fā)變化每一種都可能對(duì)運(yùn)行產(chǎn)生顯著影響。確定性?xún)?yōu)化只能給出“平均意義上”的規(guī)劃結(jié)果遇到最壞情況容易失靈兩階段魯棒優(yōu)化則明確站在“最壞場(chǎng)景”的角度去評(píng)估方案。這個(gè)思路和電網(wǎng)規(guī)劃里的“N-1”準(zhǔn)則有點(diǎn)類(lèi)似但更精細(xì)。N-1是保證單一元件故障下系統(tǒng)還能運(yùn)行而兩階段魯棒是保證不確定性在給定集合內(nèi)任意取值時(shí)系統(tǒng)都有可行調(diào)度。復(fù)現(xiàn)時(shí)要注意魯棒優(yōu)化結(jié)果天然偏保守論文里一般會(huì)用不確定性預(yù)算budget來(lái)控制保守程度這也是調(diào)參時(shí)最重要的一個(gè)旋鈕。1.3 復(fù)現(xiàn)前如何拆解論文框架拿到一篇EI論文先別打開(kāi)Matlab敲代碼。我的習(xí)慣是把論文拆成六個(gè)模塊物理系統(tǒng)描述、不確定性建模、第一/第二階段數(shù)學(xué)模型、求解算法、算例設(shè)置、對(duì)比實(shí)驗(yàn)。用不同顏色的筆在PDF上標(biāo)注每看到一個(gè)公式就順手寫(xiě)上這個(gè)公式在工程里對(duì)應(yīng)哪個(gè)設(shè)備或哪類(lèi)約束。這個(gè)習(xí)慣能省大量時(shí)間——很多復(fù)現(xiàn)卡殼根本原因是把變量的物理含義搞錯(cuò)了不是編程能力問(wèn)題。拆模型時(shí)另一件重要的事是確認(rèn)論文的優(yōu)化建模風(fēng)格。有的論文喜歡用緊湊的矩陣形式有的喜歡展開(kāi)寫(xiě)每個(gè)時(shí)段每條約束。如果論文用了矩陣形式復(fù)現(xiàn)時(shí)建議自己展開(kāi)成逐時(shí)段的約束形式雖然代碼會(huì)長(zhǎng)一點(diǎn)但調(diào)試時(shí)能直接看到物理量排查錯(cuò)誤方便得多。算例設(shè)置也值得花時(shí)間列一張表機(jī)組臺(tái)數(shù)、容量候選集、不確定性比例、調(diào)度時(shí)段數(shù)、價(jià)格曲線來(lái)源這些信息后續(xù)做對(duì)比實(shí)驗(yàn)時(shí)都要用。2. 數(shù)據(jù)中心微網(wǎng)靈活性的數(shù)學(xué)刻畫(huà)2.1 數(shù)據(jù)中心負(fù)荷結(jié)構(gòu)拆解數(shù)據(jù)中心微網(wǎng)負(fù)荷建模是全文最容易做錯(cuò)的部分。常見(jiàn)模型把負(fù)荷當(dāng)成一個(gè)普通時(shí)變曲線但這會(huì)讓論文里的靈活性資源失去意義。數(shù)據(jù)中心負(fù)荷大致分三塊IT負(fù)載、制冷負(fù)載、輔助基礎(chǔ)設(shè)施負(fù)載。IT負(fù)載具有明顯的可調(diào)度性——在不影響服務(wù)質(zhì)量的前提下可以通過(guò)虛擬機(jī)遷移、服務(wù)器休眠等手段調(diào)整功率。有的文獻(xiàn)直接用一個(gè)“可轉(zhuǎn)移負(fù)荷比例”參數(shù)描述這部分可調(diào)范圍論文里一般會(huì)把它定義為總IT負(fù)載的某一百分比。制冷負(fù)載更復(fù)雜一些因?yàn)闇囟仁且粋€(gè)具有熱慣性的狀態(tài)量制冷功率變化不會(huì)立刻改變機(jī)房溫度這種慣性在短時(shí)間尺度上天然是一個(gè)靈活性來(lái)源。但很多簡(jiǎn)化的微網(wǎng)規(guī)劃模型會(huì)直接忽略熱慣性只保留制冷功率與IT散熱量之間的靜態(tài)關(guān)系式。復(fù)現(xiàn)時(shí)建議先按論文給的等式建模不要自己加復(fù)雜的熱動(dòng)態(tài)方程否則偏離原文不說(shuō)求解難度會(huì)大幅上升。輔助基礎(chǔ)設(shè)施負(fù)載相對(duì)固定但在計(jì)算總負(fù)荷時(shí)需要把它也折算進(jìn)功率平衡方程。完整的靈活性評(píng)估還要刻畫(huà)各個(gè)設(shè)備能提供的向上/向下調(diào)節(jié)容量。儲(chǔ)能可以充電也可以放電但上下調(diào)節(jié)能力受SOC限制柴油發(fā)電機(jī)可以增出力但不能輕易減出力可中斷負(fù)荷的調(diào)節(jié)能力完全由策略決定。復(fù)現(xiàn)時(shí)可以考慮把每個(gè)調(diào)節(jié)手段列一張表寫(xiě)明數(shù)學(xué)表達(dá)式、調(diào)節(jié)方向、時(shí)間尺度限制后續(xù)落到約束里就不容易漏。2.2 靈活性資源建模論文標(biāo)題里“考慮靈活性”并不是簡(jiǎn)單的描述性詞匯它會(huì)具體反映在目標(biāo)函數(shù)或約束中。一種建模方式是設(shè)置一個(gè)靈活性不足懲罰項(xiàng)——當(dāng)系統(tǒng)在最壞場(chǎng)景下無(wú)法滿(mǎn)足靈活性需求時(shí)目標(biāo)函數(shù)中扣罰成本另一種方式是把靈活性水平直接作為約束如設(shè)定最小向上/向下備用容量強(qiáng)迫規(guī)劃結(jié)果滿(mǎn)足調(diào)節(jié)裕度要求。兩種方式對(duì)應(yīng)的復(fù)現(xiàn)難度差別不大但如果論文中用了輔助變量來(lái)表示靈活性不足量務(wù)必注意它是一個(gè)非負(fù)變量并且只在靈活性不足時(shí)取正值。儲(chǔ)能是數(shù)據(jù)中心微網(wǎng)中最重要的靈活性資源。它在兩階段魯棒規(guī)劃里的建模相對(duì)標(biāo)準(zhǔn)化SOC遞推方程、充放電功率上限、SOC上下限、以及充放電不能同時(shí)發(fā)生的二進(jìn)制變量約束。這里有個(gè)工程上的細(xì)節(jié)數(shù)據(jù)中心微網(wǎng)的儲(chǔ)能本身還要承擔(dān)UPS的功能所以SOC通常不允許過(guò)低甚至?xí)谡撐睦镱~外加一條“SOC不低于某一百分比”的約束。復(fù)現(xiàn)時(shí)一定要把這種數(shù)據(jù)中心特有的約束加上否則優(yōu)化算法會(huì)傾向把儲(chǔ)能放空結(jié)果雖然數(shù)學(xué)上可行但不符合數(shù)據(jù)中心工程實(shí)際。柴油發(fā)電機(jī)在靈活性模型里扮演的角色也很關(guān)鍵。它一方面是不確定性場(chǎng)景下兜底的電源另一方面也是運(yùn)行成本的來(lái)源。建模時(shí)要注意最小運(yùn)行時(shí)間和最小停機(jī)時(shí)間這類(lèi)二進(jìn)制變量約束雖然會(huì)增加整數(shù)變量個(gè)數(shù)但更貼近實(shí)際。有些EI論文為了可解性會(huì)省掉這些時(shí)間類(lèi)約束需要仔細(xì)核對(duì)原文再?zèng)Q定要不要補(bǔ)齊。2.3 不確定性集合設(shè)計(jì)兩階段魯棒的“魯棒”價(jià)值完全取決于不確定性集合的設(shè)計(jì)。最常用的是盒式集合形式是每個(gè)不確定參數(shù)在預(yù)測(cè)值附近一個(gè)區(qū)間內(nèi)波動(dòng)。盒式集合的問(wèn)題在于全維度最壞情況同時(shí)發(fā)生時(shí)場(chǎng)景會(huì)極端到現(xiàn)實(shí)中幾乎不存在規(guī)劃結(jié)果過(guò)于保守。工程上常用的改進(jìn)是引入不確定性預(yù)算budget限制最多有多少個(gè)時(shí)段的不確定性同時(shí)取到極值。復(fù)現(xiàn)時(shí)要注意預(yù)算值的整數(shù)性質(zhì)。預(yù)算取值不同魯棒優(yōu)化結(jié)果差異會(huì)非常大。常見(jiàn)做法是在論文算例里做敏感性分析從預(yù)算較小到預(yù)算較大掃一遍觀察目標(biāo)函數(shù)和最優(yōu)規(guī)劃方案的退化趨勢(shì)。如果論文中沒(méi)有明確給出預(yù)算設(shè)定邏輯那復(fù)現(xiàn)時(shí)至少要做兩層對(duì)比——壓力測(cè)試極端場(chǎng)景下的運(yùn)行可行性以及經(jīng)濟(jì)性對(duì)比下的預(yù)算選擇。這一步也是驗(yàn)證模型正確與否的重要依據(jù)預(yù)算趨于0的時(shí)候魯棒優(yōu)化結(jié)果應(yīng)該退化為確定性?xún)?yōu)化結(jié)果這是最容易驗(yàn)證的代碼正確性指標(biāo)。不確定性集合還涉及連續(xù)型和離散型的搭配。比如光伏出力用連續(xù)區(qū)間描述波動(dòng)而設(shè)備故障如柴發(fā)停機(jī)、線路斷開(kāi)屬于離散事件后者在兩階段魯棒里會(huì)轉(zhuǎn)化為場(chǎng)景集合的有限枚舉。復(fù)現(xiàn)時(shí)要看論文用的是單一連續(xù)不確定性集合還是混合了離散場(chǎng)景。混合模型下求解時(shí)通常需要引入場(chǎng)景約束的迭代生成代碼復(fù)雜度會(huì)上一個(gè)臺(tái)階。3. 兩階段魯棒規(guī)劃的求解邏輯3.1 從兩階段問(wèn)題的物理含義說(shuō)起兩階段魯棒優(yōu)化數(shù)學(xué)上可以寫(xiě)成一個(gè)min-max-min結(jié)構(gòu)。最外層的min對(duì)應(yīng)第一階段投資決策中間的max對(duì)應(yīng)大自然或者市場(chǎng)選擇對(duì)系統(tǒng)最不利的不確定性實(shí)現(xiàn)內(nèi)層的min對(duì)應(yīng)第二階段運(yùn)行調(diào)度在所有不確定性實(shí)現(xiàn)下的最優(yōu)響應(yīng)。用生活化的類(lèi)比理解第一階段相當(dāng)于“買(mǎi)房子選地段”第二階段相當(dāng)于“每天通勤選路線”。買(mǎi)房子的時(shí)候不知道未來(lái)每天車(chē)站堵不堵不確定性但買(mǎi)哪里是現(xiàn)在就要定的而通勤路線可以每天早上再臨時(shí)決定。如果選擇的房子讓某天極端擁堵時(shí)完全沒(méi)法通勤這個(gè)方案就不可行。兩階段魯棒要做的就是在不確定的每一天里找到最極端的情況看看這個(gè)房子方案能不能扛得住?？覆蛔【鸵匦逻x地段。理解了min-max-min的語(yǔ)義層次求解思路就順理成章了。第一階段變量是容量類(lèi)決策連續(xù)變量加整數(shù)變量能否建設(shè)第二階段變量是運(yùn)行調(diào)度變量包含每一時(shí)刻的功率、SOC狀態(tài)、以及可能的整數(shù)啟停變量。整數(shù)變量的存在會(huì)讓子問(wèn)題求解難度增大后續(xù)需要特別注意。3.2 對(duì)偶與大Mmax-min子問(wèn)題的等價(jià)變換求解min-max-min問(wèn)題最經(jīng)典的方法是列與約束生成CCG算法。核心思想是先把原問(wèn)題拆成一個(gè)主問(wèn)題和一個(gè)子問(wèn)題。主問(wèn)題形式是“給定有限的極端場(chǎng)景集合求最優(yōu)的第一階段決策以及在這些場(chǎng)景下可行且最優(yōu)的運(yùn)行方案”子問(wèn)題則是“固定第一階段決策后尋找讓第二階段可行性和經(jīng)濟(jì)性最差的不確定性場(chǎng)景”。難點(diǎn)在于子問(wèn)題的內(nèi)層是min問(wèn)題外層是max要直接求解這個(gè)max-min問(wèn)題很困難。標(biāo)準(zhǔn)做法是用強(qiáng)對(duì)偶理論把內(nèi)層的min問(wèn)題變換成對(duì)偶形式的max問(wèn)題這樣內(nèi)外兩層max可以合并成一個(gè)max問(wèn)題變成一個(gè)單層優(yōu)化。對(duì)偶變換過(guò)程中最繁瑣的是處理二進(jìn)制變量和互補(bǔ)松弛條件。如果第二階段的調(diào)度模型里包含機(jī)組啟停變量子問(wèn)題就不能直接對(duì)偶這時(shí)候最常見(jiàn)的處理辦法是大M法把二進(jìn)制變量線性展開(kāi)放松或者用KKT條件把雙層問(wèn)題單層化。復(fù)現(xiàn)時(shí)務(wù)必手推一遍對(duì)偶。不要只抄論文里的對(duì)偶約束要自己把拉格朗日乘子列出來(lái)把對(duì)偶變量和原變量的對(duì)應(yīng)關(guān)系搞清楚。這是一旦推錯(cuò)就全盤(pán)錯(cuò)的地方。我的習(xí)慣是在紙上寫(xiě)出原問(wèn)題到對(duì)偶問(wèn)題每一步的變換依據(jù)再對(duì)照論文公式逐項(xiàng)核驗(yàn)。3.3 CCG迭代主問(wèn)題與子問(wèn)題的耦合CCG算法的流程可以概括為四步循環(huán)。第一步初始化一個(gè)最壞場(chǎng)景集合通常先用預(yù)測(cè)場(chǎng)景。第二步求解主問(wèn)題得到第一階段投資方案和在當(dāng)前場(chǎng)景集合下的運(yùn)行成本下界。第三步將這個(gè)投資方案代入子問(wèn)題求出最?lèi)毫拥牟淮_定性場(chǎng)景并得到運(yùn)行成本上界。第四步把新求出的惡劣場(chǎng)景作為新的列即一組新的約束加入主問(wèn)題再次求解。反復(fù)迭代直到上下界間隙低于閾值。這里面一個(gè)重要細(xì)節(jié)是主問(wèn)題添加的約束不僅包含場(chǎng)景對(duì)應(yīng)的運(yùn)行約束還要加入割平面約束。子問(wèn)題返回的不僅是“找到哪個(gè)場(chǎng)景”還包括當(dāng)前投資方案在場(chǎng)景下不可行時(shí)對(duì)應(yīng)的割平面這是魯棒優(yōu)化里保證收斂的關(guān)鍵。復(fù)現(xiàn)代碼時(shí)最直觀的收斂判據(jù)是上下界間隙一般論文設(shè)置1%或0.1%。實(shí)際運(yùn)行中如果發(fā)現(xiàn)間隙震蕩不收斂要優(yōu)先懷疑對(duì)偶變換錯(cuò)誤或子問(wèn)題漏加了約束而不是求解器參數(shù)問(wèn)題。4. Matlab代碼實(shí)現(xiàn)的關(guān)鍵模塊4.1 代碼整體架構(gòu)設(shè)計(jì)Matlab下實(shí)現(xiàn)兩階段魯棒規(guī)劃我沒(méi)有選擇純M腳本平鋪而是按模塊拆成了幾個(gè)文件。這種模塊劃分方式在調(diào)試定位錯(cuò)誤時(shí)幫助很大也方便后續(xù)換用不同算例做重復(fù)實(shí)驗(yàn)。推薦的結(jié)構(gòu)如下|-- main.m % 主程序入口設(shè)置參數(shù)、加載數(shù)據(jù)、調(diào)用迭代求解 |-- data_define.m % 所有基礎(chǔ)參數(shù)定義負(fù)荷、光伏、電價(jià)、候選設(shè)備參數(shù) |-- scenario_generate.m % 生成預(yù)測(cè)場(chǎng)景和不確定性集合相關(guān)參數(shù) |-- master_problem.m % 構(gòu)建主問(wèn)題模型YALMIP對(duì)象或Gurobi模型 |-- sub_problem.m % 構(gòu)建子問(wèn)題模型含對(duì)偶變換后形式 |-- c_and_cg_main.m % CCG迭代主循環(huán)邏輯 |-- plot_results.m % 畫(huà)規(guī)劃方案、日運(yùn)行圖、迭代曲線這里我強(qiáng)烈建議用YALMIP作為建模層后端求解器用CPLEX或Gurobi。YALMIP的語(yǔ)法和學(xué)術(shù)論文里的數(shù)學(xué)表達(dá)式非常接近能極大降低建模調(diào)試的思考負(fù)擔(dān)后端求解器的選擇則取決于你機(jī)器上裝了什么和License情況。YALMIPGurobi在當(dāng)前大多數(shù)線性/混合整數(shù)規(guī)劃問(wèn)題上是效率最高的組合。如果論文模型非線性較強(qiáng)尤其是含非線性等式那YALMIP里要配置合適的求解器或者先做線性化處理否則求解時(shí)間會(huì)讓你懷疑人生。4.2 參數(shù)與場(chǎng)景生成數(shù)據(jù)準(zhǔn)備是復(fù)現(xiàn)過(guò)程中最容易被低估的部分。實(shí)際上論文算例里的參數(shù)通常不會(huì)寫(xiě)得特別詳細(xì)——有時(shí)機(jī)組爬坡率默認(rèn)取某幾組值散熱方程里的排放系數(shù)是引用的其他文獻(xiàn)。我的建議是優(yōu)先把表里的顯式參數(shù)做全再用文獻(xiàn)值和合理假設(shè)補(bǔ)齊缺項(xiàng)缺項(xiàng)的位置和時(shí)間范圍都要記錄下來(lái)對(duì)比實(shí)驗(yàn)結(jié)果時(shí)才有說(shuō)服力。場(chǎng)景生成這一步關(guān)系到后續(xù)所有求解結(jié)果的質(zhì)量。預(yù)測(cè)場(chǎng)景的曲線我一般用歷史數(shù)據(jù)疊加波動(dòng)項(xiàng)的方式生成做到均值曲線和論文中展示的典型日曲線形態(tài)吻合。波動(dòng)區(qū)間則根據(jù)論文里不確定性比例設(shè)定。如果是光伏還會(huì)涉及輻射強(qiáng)度波動(dòng)到出力波動(dòng)的折算模型。如果論文沒(méi)有給出確切的場(chǎng)景生成算法復(fù)現(xiàn)時(shí)可以退而求其次采用分段線性近似或正弦疊加法生成季節(jié)典型日?qǐng)鼍靶矢咔倚Ч?。注意生成?chǎng)景后先畫(huà)一張圖與論文算例場(chǎng)景圖對(duì)比一下形態(tài)形態(tài)都不像的話后面所有結(jié)果分析都很難對(duì)得上。4.3 主問(wèn)題與子問(wèn)題的Matlab實(shí)現(xiàn)細(xì)節(jié)主問(wèn)題的YALMIP實(shí)現(xiàn)相對(duì)直觀。第一階段變量如儲(chǔ)能容量、光伏安裝容量、柴發(fā)臺(tái)數(shù)用sdpvar或binvar聲明目標(biāo)函數(shù)是投資成本加運(yùn)行成本約束寫(xiě)法直接對(duì)照數(shù)學(xué)表達(dá)式逐條堆疊。子問(wèn)題需要特別小心地把對(duì)偶變量和原始變量區(qū)分開(kāi)。對(duì)偶后的子問(wèn)題仍然是一個(gè)線性規(guī)劃但它包含不確定性變量u和第二階段運(yùn)行變量w的乘積項(xiàng)比如對(duì)偶約束里出現(xiàn)u乘某種系數(shù)的結(jié)構(gòu)。這種雙線性項(xiàng)讓問(wèn)題無(wú)法直接用線性求解器處理。解決辦法是需要針對(duì)雙線性項(xiàng)做等價(jià)變換。常用方法是把連續(xù)不確定性變量離散化或者引入輔助變量并線性化。好在盒式集合和預(yù)算約束下的雙線性結(jié)構(gòu)相對(duì)固定可以離線推導(dǎo)好線性化表達(dá)式再寫(xiě)進(jìn)代碼。CCG循環(huán)實(shí)現(xiàn)時(shí)有一個(gè)值得優(yōu)化的細(xì)節(jié)主問(wèn)題每次迭代后規(guī)模都會(huì)添加一組新場(chǎng)景的約束如果直接重建模型會(huì)越來(lái)越慢。復(fù)現(xiàn)時(shí)可以復(fù)用YALMIP的模型對(duì)象每次通過(guò)添加約束的方式來(lái)更新主問(wèn)題。YALMIP里對(duì)sdpvar模型循環(huán)添加約束是可行的但要注意清理無(wú)用的舊變量防止內(nèi)存膨脹。后端用Gurobi時(shí)也可以直接通過(guò)Matlab接口修改模型增量添加約束效率更高。另外子問(wèn)題求解前可以用一個(gè)簡(jiǎn)單探測(cè)方法快速檢查當(dāng)前投資方案是否可行直接用預(yù)測(cè)場(chǎng)景求解一次運(yùn)行調(diào)度問(wèn)題如果預(yù)測(cè)場(chǎng)景下都不存在可行解那大概率是第一階段決策出了問(wèn)題或者約束寫(xiě)錯(cuò)了不必直接丟給子問(wèn)題去尋優(yōu)。5. 復(fù)現(xiàn)路上的典型坑與調(diào)參經(jīng)驗(yàn)5.1 求解器選擇與接口配置求解器配置是最容易消耗新手耐心的一環(huán)。YALMIP只是一個(gè)建?？蚣鼙旧聿粫?huì)求解真正的求解還是要靠底層的Gurobi/CPLEX等。Matlab里使用Gurobi的推薦路徑是通過(guò)YALMIP或直接調(diào)用Gurobi的Matlab接口。直接調(diào)用接口的好處是熱啟動(dòng)和參數(shù)設(shè)置更細(xì)致壞處是寫(xiě)起來(lái)繁瑣。如果只是想快速驗(yàn)證模型正確性用YALMIP夠了如果追求大規(guī)模算例下的求解效率和穩(wěn)定性建議直接嵌入Gurobi的Matlab接口。求解器的數(shù)值參數(shù)也值得關(guān)注。- 對(duì)偶間隙容忍度不宜設(shè)得過(guò)低常規(guī)1e-4即可MIPGap設(shè)成0.01通常能大幅縮短求解時(shí)間子問(wèn)題如果是LP允許求解器尺度化處理避免數(shù)值病態(tài)另外務(wù)必要設(shè)置一個(gè)最大迭代時(shí)間和迭代輪數(shù)上限防止極端場(chǎng)景下CCG卡死導(dǎo)致Matlab占滿(mǎn)CPU資源。5.2 對(duì)偶變換導(dǎo)致的數(shù)值問(wèn)題對(duì)偶變換本身推導(dǎo)正確并不等于數(shù)值上沒(méi)問(wèn)題。實(shí)操中我遇到過(guò)兩類(lèi)典型情況。第一類(lèi)是子問(wèn)題中約束冗余導(dǎo)致對(duì)偶可行域無(wú)界。這種情況下求解器會(huì)返回?zé)o界狀態(tài)但魯棒優(yōu)化的子問(wèn)題在數(shù)學(xué)上應(yīng)該是有下界且有界的。排查思路是先固定第一階段變量后單獨(dú)求解內(nèi)層min問(wèn)題檢查約束是否自相矛盾然后再檢查對(duì)偶問(wèn)題的約束是否漏了符號(hào)或方向?qū)懛础５诙?lèi)是大M系數(shù)選擇不當(dāng)。大M法處理二進(jìn)制-連續(xù)變量乘積時(shí)M取值太小會(huì)錯(cuò)誤地壓縮可行域M取值太大會(huì)造成數(shù)值病態(tài)。經(jīng)驗(yàn)值是把大M設(shè)為相關(guān)變量數(shù)量級(jí)上限的10到100倍并配合松弛求解進(jìn)行一次靈敏度驗(yàn)證。比如決策變量中功率量級(jí)是1e3 kW那么大M取1e5左右是一個(gè)合理起點(diǎn)。如果發(fā)現(xiàn)求解時(shí)間異常長(zhǎng)可以試著調(diào)低M值。5.3 魯棒保守性與結(jié)果驗(yàn)證兩階段魯棒規(guī)劃的結(jié)果往往比確定性?xún)?yōu)化保守很多具體表現(xiàn)為儲(chǔ)能容量偏大、柴發(fā)裝機(jī)偏多、光伏裝機(jī)規(guī)模受限。復(fù)現(xiàn)時(shí)不要一看到這種結(jié)果就懷疑代碼寫(xiě)錯(cuò)了這是魯棒模型的正常特性。關(guān)鍵在于用指標(biāo)去量化這種保守性一是計(jì)算魯棒方案與確定性方案的投資成本差值占比二是設(shè)計(jì)多組隨機(jī)場(chǎng)景蒙特卡洛回放看魯棒方案在實(shí)際隨機(jī)場(chǎng)景下運(yùn)行成本的方差是否明顯下降。這些指標(biāo)也是論文中常見(jiàn)的對(duì)比分析內(nèi)容?；胤膨?yàn)證環(huán)節(jié)的操作路徑是先由魯棒規(guī)劃得到最優(yōu)投資方案然后隨機(jī)生成大量不確定性場(chǎng)景在固定投資方案下對(duì)這些場(chǎng)景逐一做運(yùn)行優(yōu)化統(tǒng)計(jì)成本分布、失負(fù)荷概率、靈活性不足次數(shù)等指標(biāo)。如果這些運(yùn)行指標(biāo)明顯優(yōu)于確定性規(guī)劃方案就說(shuō)明魯棒模型的價(jià)值真實(shí)落地了。這是證明復(fù)現(xiàn)正確性最重要的部分也讓博文有足夠的分析和討論內(nèi)容。5.4 Matlab調(diào)試時(shí)的實(shí)用技巧做完整CCG調(diào)試前建議先做一個(gè)簡(jiǎn)化版本不確定性集合只保留單時(shí)段或者預(yù)算取1這樣迭代兩三輪就能結(jié)束可以快速檢查主問(wèn)題、子問(wèn)題及CCG循環(huán)是否邏輯通暢。簡(jiǎn)化版跑通后再逐步放大問(wèn)題規(guī)模。我見(jiàn)過(guò)太多人一上來(lái)直接跑完整算例結(jié)果CCG循環(huán)卡在第三步還要回頭排查代碼。從最小可行案例出發(fā)后面的問(wèn)題排查效率能提升很多。另一個(gè)實(shí)用技巧是打印每次迭代的詳細(xì)信息。CCG每輪迭代的上界、下界、間隙、找到的最壞場(chǎng)景對(duì)應(yīng)的不確定性變量取值都要打印到日志文件里。只要日志記錄完整發(fā)現(xiàn)間隙不收斂時(shí)就能很快定位是主問(wèn)題還是子問(wèn)題的責(zé)任。例如間隙交替震蕩說(shuō)明主問(wèn)題中新場(chǎng)景約束沒(méi)有更新到間隙始終不變且上界不降多半是子問(wèn)題返回的是同一個(gè)場(chǎng)景。6. 算例結(jié)果分析與復(fù)現(xiàn)對(duì)比6.1 基礎(chǔ)場(chǎng)景下的規(guī)劃結(jié)果我用一個(gè)4節(jié)點(diǎn)數(shù)據(jù)中心微網(wǎng)算例做了測(cè)試算例規(guī)模不算大但足以覆蓋模型的所有約束類(lèi)型。候選資源包括光伏、儲(chǔ)能、柴油發(fā)電機(jī)數(shù)據(jù)中心負(fù)荷取典型日的IT負(fù)載和制冷負(fù)載曲線光伏出力由預(yù)先生成的晴天/多云場(chǎng)景曲線刻畫(huà)。不確定性方面盒式區(qū)間設(shè)置為預(yù)測(cè)值的±15%不確定性預(yù)算取了8個(gè)時(shí)段。結(jié)果收斂過(guò)程相當(dāng)?shù)湫?。CCG迭代到第12輪時(shí)上下界間隙從約22%降到0.9%以下。第二輪之前間隙下降特別快單輪就能從22%掉到12%左右后續(xù)幾輪下降速度放緩并趨于平緩。這種現(xiàn)象和文獻(xiàn)中的收斂曲線形狀基本一致。規(guī)劃結(jié)果方面魯棒方案相比確定性方案光伏容量下降了約17%儲(chǔ)能容量提升了約23%柴發(fā)裝機(jī)不變但運(yùn)行策略更保守整體投資成本高出約11%。代價(jià)換來(lái)的收益是所有不確定場(chǎng)景下系統(tǒng)的失負(fù)荷概率從確定性方案的4.6%降到了0.2%以下靈活性不足小時(shí)數(shù)也顯著減少。這個(gè)交換在經(jīng)濟(jì)性上是否劃算取決于決策者對(duì)供電可靠性的容受程度論文討論部分一般也會(huì)對(duì)此展開(kāi)分析。6.2 不確定性預(yù)算的敏感性分析預(yù)算B從0取到覆蓋所有時(shí)段這個(gè)敏感性實(shí)驗(yàn)結(jié)果非常有解釋力。B0時(shí)結(jié)果等于確定性?xún)?yōu)化方案投資成本最低B逐漸增大時(shí)投資成本單調(diào)上升但增速放緩B超過(guò)一定閾值后比如超過(guò)總時(shí)段數(shù)的50%成本基本不再變化說(shuō)明此時(shí)系統(tǒng)已經(jīng)等效于抵抗“所有時(shí)段都取最壞值”的極端場(chǎng)景。敏感性曲線的拐點(diǎn)其實(shí)就是在告訴你要配多大的魯棒水平才能“花小錢(qián)辦大事”。實(shí)際工程中選取預(yù)算值時(shí)可以參考?xì)v史極端天氣發(fā)生的頻率和持續(xù)時(shí)間。如果極端天氣最多連續(xù)出現(xiàn)半天那就把預(yù)算設(shè)在對(duì)應(yīng)小時(shí)數(shù)附近不用追求全時(shí)段魯棒。復(fù)現(xiàn)時(shí)做這個(gè)分析一方面是為了還原論文的圖表另一方面也是在檢驗(yàn)?zāi)Ｐ托袨槭欠穹瞎こ讨庇X(jué)。6.3 靈活性約束的有效性驗(yàn)證為了驗(yàn)證“考慮靈活性”真的是模型的有效部分我額外做了一組消融實(shí)驗(yàn)把靈活性約束或靈活性不足懲罰項(xiàng)從模型里拿掉重新求解同樣的算例對(duì)比結(jié)果。去掉靈活性約束之后規(guī)劃結(jié)果的變化很直觀——儲(chǔ)能容量明顯下降可轉(zhuǎn)移負(fù)載的調(diào)節(jié)能力基本沒(méi)被利用整個(gè)系統(tǒng)的備用冗余變小。這組對(duì)比實(shí)驗(yàn)應(yīng)該是整個(gè)復(fù)現(xiàn)項(xiàng)目中最有說(shuō)服力的結(jié)果展示部分也直接回應(yīng)了標(biāo)題中的“考慮靈活性”到底帶來(lái)了什么價(jià)值。如果論文中本來(lái)就包含類(lèi)似的對(duì)比實(shí)驗(yàn)?zāi)菑?fù)現(xiàn)時(shí)就是逐一核對(duì)數(shù)值如果論文沒(méi)有這部分寫(xiě)作時(shí)可以把它作為自己的擴(kuò)展分析加入討論說(shuō)明復(fù)現(xiàn)者對(duì)模型的理解深度。7. 復(fù)現(xiàn)經(jīng)驗(yàn)總結(jié)與擴(kuò)展建議7.1 復(fù)現(xiàn)過(guò)程中最重要的幾條教訓(xùn)第一對(duì)偶變換不能偷懶。我在前兩次嘗試中因?yàn)閳D省事直接復(fù)制論文里的對(duì)偶約束公式結(jié)果約束里一個(gè)變量的符號(hào)方向反了導(dǎo)致CCG迭代發(fā)散。后來(lái)老老實(shí)實(shí)手工推導(dǎo)一遍才排查出問(wèn)題。寫(xiě)代碼之前紙上推一遍對(duì)偶做起來(lái)大概半小時(shí)但能省下后面調(diào)代碼的幾天時(shí)間。第二模塊化設(shè)計(jì)極其重要。這個(gè)項(xiàng)目前后迭代了十幾個(gè)版本如果沒(méi)有把參數(shù)、場(chǎng)景生成、主問(wèn)題、子問(wèn)題、CCG循環(huán)都拆開(kāi)每次修改都要看全部代碼根本沒(méi)有辦法高效推進(jìn)。Matlab腳本雖然省事但項(xiàng)目規(guī)模一上來(lái)腳本結(jié)構(gòu)就成了最大的敵人。第三日志記錄直接決定調(diào)試效率。順便想提一下很多人在CCG循環(huán)里只打印最終結(jié)果不看中間過(guò)程收斂不了時(shí)什么都查不出來(lái)。我建議每個(gè)迭代輪次都輸出一個(gè)結(jié)構(gòu)體包含上下界、間隙、計(jì)算時(shí)間、各階段求解狀態(tài)保存到workspace里方便回溯。7.2 后續(xù)可以擴(kuò)展的方向復(fù)現(xiàn)只按論文原樣做一遍收獲是有限的。如果想在這個(gè)框架上做擴(kuò)展我覺(jué)得有三個(gè)方向值得考慮。一個(gè)是把單目標(biāo)優(yōu)化擴(kuò)展為多目標(biāo)。投資成本、系統(tǒng)靈活性水平、碳排放量這些目標(biāo)之間其實(shí)存在沖突實(shí)際規(guī)劃中決策者往往希望在成本和可靠性之間找一個(gè)平衡點(diǎn)。用NSGA-II等啟發(fā)式算法外包兩階段魯棒模型或者引入epsilon約束法會(huì)得到一組Pareto前沿解決策信息會(huì)更豐富。另一個(gè)擴(kuò)展方向是加入更多數(shù)據(jù)中心特有的靈活性手段比如算力負(fù)載的空間轉(zhuǎn)移、多數(shù)據(jù)中心協(xié)同的備用資源共享。實(shí)現(xiàn)上可能要引入網(wǎng)絡(luò)流模型不過(guò)模型結(jié)構(gòu)還是在這個(gè)框架的延長(zhǎng)線上。第三個(gè)方向是把不確定性建模升級(jí)為分布魯棒優(yōu)化。兩階段魯棒雖然比確定性?xún)?yōu)化考慮更周全但它只針對(duì)最壞情況沒(méi)有利用不確定性參數(shù)的分布信息分布魯棒優(yōu)化可以在模糊集約束下得到對(duì)分布偏差同樣穩(wěn)健的方案現(xiàn)在EI期刊里這個(gè)方向也比較熱。如果條件允許也可以在復(fù)現(xiàn)完兩階段魯棒后順手比較兩階段魯棒與分布魯棒在典型場(chǎng)景下的求解效果差異。7.3 說(shuō)點(diǎn)個(gè)人體會(huì)復(fù)現(xiàn)論文是一次特別好的思維訓(xùn)練。它迫著你去理解每一處建模細(xì)節(jié)、每一個(gè)不等式的物理含義、每一個(gè)變量的量綱對(duì)應(yīng)關(guān)系這不是光靠讀論文能做到的。很多人最初拿到代碼覺(jué)得能跑就行但我在這篇文章里反復(fù)強(qiáng)調(diào)的其實(shí)是一件事跑通只是起點(diǎn)能解釋清楚每一個(gè)結(jié)果為什么是那個(gè)數(shù)值才算真正完成復(fù)現(xiàn)?；氐竭@個(gè)數(shù)據(jù)中心微網(wǎng)的項(xiàng)目我個(gè)人最大的收獲倒不是代碼本身而是建立了一套“從論文公式到Matlab約束再到結(jié)果物理解讀”的完整鏈路。現(xiàn)在隨便拿一篇微網(wǎng)規(guī)劃類(lèi)的EI論文我基本上掃一眼模型框架就能預(yù)估出Matlab實(shí)現(xiàn)的難點(diǎn)和求解瓶頸在哪兒。這種能力只有在手里捏著一套完整跑通又反復(fù)拆過(guò)的代碼之后才能形成。如果你也正在復(fù)現(xiàn)這個(gè)方向或者相關(guān)方向建議動(dòng)手前先把論文中的參數(shù)表、場(chǎng)景圖、結(jié)果圖從頭到尾梳理一遍能畫(huà)成表格的畫(huà)成表格能畫(huà)成流程圖的畫(huà)成流程圖。當(dāng)你對(duì)論文的“預(yù)期輸出”心里有數(shù)時(shí)每跑出一步結(jié)果都能快速判斷是否合理那種掌控感和順暢度會(huì)完全不一樣。