做電力系統(tǒng)優(yōu)化的人應(yīng)該都有這種體感儲(chǔ)能電站的論文這兩年多到看標(biāo)題就想劃走但真正能直接拿去改改參數(shù)、換換數(shù)據(jù)就能跑的代碼反而少見。這篇要聊的模型標(biāo)題很長(zhǎng)——《考慮碳交易與電網(wǎng)交互波動(dòng)懲罰的共享儲(chǔ)能電站優(yōu)化配置與調(diào)度模型研究Matlab代碼實(shí)現(xiàn)》但落點(diǎn)其實(shí)很實(shí)在當(dāng)儲(chǔ)能電站不是給單一風(fēng)電場(chǎng)或光伏電站服務(wù)而是要同時(shí)應(yīng)對(duì)多個(gè)用戶、還要把碳交易收益和電網(wǎng)交互功率波動(dòng)一起塞進(jìn)優(yōu)化目標(biāo)時(shí)容量怎么配、功率怎么調(diào)、收益怎么算這三個(gè)問題會(huì)互相糾纏。這篇文章我會(huì)用MatlabYalmip的視角把模型從問題拆解到代碼實(shí)現(xiàn)講清楚包括每一層的決策變量、約束條件怎么寫、波動(dòng)懲罰項(xiàng)怎么處理才不會(huì)把求解器搞崩以及跑算例時(shí)最容易踩的幾個(gè)坑。這篇內(nèi)容的適用對(duì)象很明確正在做儲(chǔ)能優(yōu)化配置方向畢業(yè)設(shè)計(jì)或論文復(fù)現(xiàn)的研究生以及剛接觸共享儲(chǔ)能商業(yè)模式、想把碳交易機(jī)制納入規(guī)劃的工程師。如果你只是想知道共享儲(chǔ)能能賺多少錢這篇可能偏硬但如果你需要的是一個(gè)能落地、能復(fù)現(xiàn)、能改參數(shù)出圖表的建模思路那這篇應(yīng)該能幫你省下不少翻文獻(xiàn)的時(shí)間。1. 這個(gè)模型到底在解決什么問題共享儲(chǔ)能背后的三重矛盾1.1 共享儲(chǔ)能為什么比自建儲(chǔ)能更值得研究先說(shuō)共享儲(chǔ)能電站和傳統(tǒng)儲(chǔ)能電站的本質(zhì)區(qū)別。傳統(tǒng)儲(chǔ)能通常是新能源場(chǎng)站自建目的很單純棄風(fēng)棄光多了就充電峰谷電價(jià)差大了就放電收益來(lái)源單一容量配置只需要服務(wù)自家場(chǎng)站。而共享儲(chǔ)能電站是獨(dú)立投資、獨(dú)立運(yùn)營(yíng)的第三方主體它同時(shí)服務(wù)多個(gè)新能源場(chǎng)站或多個(gè)用戶側(cè)主體通過容量租賃、調(diào)峰輔助服務(wù)、峰谷套利等渠道回收成本。這種模式的好處是提高了儲(chǔ)能利用率。自建儲(chǔ)能的利用率受限于單個(gè)場(chǎng)站的發(fā)電曲線而共享儲(chǔ)能可以把多個(gè)場(chǎng)站的負(fù)荷曲線錯(cuò)峰疊加電池的充放電次數(shù)和能量吞吐量都能上去。但問題也隨之而來(lái)服務(wù)對(duì)象多了優(yōu)化配置就不能再拍腦袋定容量必須考慮不同用戶的用電特性、消納需求、以及儲(chǔ)能調(diào)度策略如何影響長(zhǎng)期收益。所以共享儲(chǔ)能的優(yōu)化配置天然是個(gè)配置調(diào)度耦合的問題——配置決定調(diào)度的物理邊界調(diào)度決定配置的經(jīng)濟(jì)回報(bào)。1.2 碳交易怎么切入儲(chǔ)能的經(jīng)濟(jì)性評(píng)估碳交易是這兩年電力系統(tǒng)優(yōu)化模型里最常被加進(jìn)去的機(jī)制。它的邏輯鏈條是這樣新能源發(fā)電替代火電出力減少了碳排放這部分減排量可以在碳交易市場(chǎng)轉(zhuǎn)化為收益反過來(lái)如果系統(tǒng)還得依賴火電平衡那就要購(gòu)買碳排放配額。共享儲(chǔ)能在這個(gè)過程中扮演的角色是減排催化劑——儲(chǔ)能吸收多余新能源電量、減少棄電相當(dāng)于間接替代了火電出力所以在目標(biāo)函數(shù)里可以給儲(chǔ)能一個(gè)碳收益項(xiàng)。但注意儲(chǔ)能的碳收益不是直接少排了多少噸碳而是圍繞系統(tǒng)的碳排放配額與實(shí)際排放量差值來(lái)計(jì)算的。這個(gè)差值怎么定義直接決定模型的目標(biāo)函數(shù)形式。有的論文用基準(zhǔn)線法有的用配額法我這里用的是比較常見的配額法系統(tǒng)根據(jù)發(fā)電量或用電量獲得初始配額實(shí)際碳排放量超過配額就要在碳市場(chǎng)購(gòu)買低于配額就可以出售獲利。1.3 電網(wǎng)交互波動(dòng)懲罰被很多論文忽略的隱性成本第三個(gè)矛盾是儲(chǔ)能電站與外部電網(wǎng)的交互功率波動(dòng)。很多基礎(chǔ)配置模型只關(guān)心儲(chǔ)能自身的充放電功率、SOC約束、以及收益最大化卻忽略了一個(gè)工程現(xiàn)實(shí)儲(chǔ)能電站從電網(wǎng)取電或向電網(wǎng)送電的功率曲線如果波動(dòng)過于劇烈會(huì)對(duì)電網(wǎng)造成沖擊尤其在接入點(diǎn)比較薄弱的區(qū)域這種沖擊會(huì)被調(diào)度中心視作不合格的并網(wǎng)行為。電網(wǎng)交互波動(dòng)懲罰就是對(duì)這種工程約束的經(jīng)濟(jì)化表達(dá)。模型里加一個(gè)懲罰項(xiàng)當(dāng)儲(chǔ)能與電網(wǎng)交互功率的變化率超過某個(gè)閾值時(shí)懲罰成本上升從而在優(yōu)化過程中逼著儲(chǔ)能調(diào)度曲線變得更平滑。這也是這個(gè)標(biāo)題里最有區(qū)分度的一點(diǎn)——大部分基礎(chǔ)儲(chǔ)能配置模型不會(huì)考慮這個(gè)約束而實(shí)際工程中它往往比碳交易更影響調(diào)度結(jié)果。2. 雙層模型拆解上層定容量下層定曲線2.1 上層決策變量容量配置不是拍腦袋共享儲(chǔ)能電站的優(yōu)化配置通常采用雙層模型結(jié)構(gòu)英文里叫Bi-level Optimization。上層是規(guī)劃層決策變量可以簡(jiǎn)化為一組儲(chǔ)能額定功率 (P_{ess}^{rated})單位MW儲(chǔ)能額定容量 (E_{ess}^{rated})單位MWh各用戶/新能源場(chǎng)站的容量租賃比例或接入容量上層的目標(biāo)函數(shù)一般是全壽命周期內(nèi)的凈收益最大化或綜合成本最小化。成本項(xiàng)包括初始投資成本功率成本容量成本、運(yùn)行維護(hù)成本、充電購(gòu)電成本收益項(xiàng)包括向用戶收取的容量租賃費(fèi)、峰谷套利收益、輔助服務(wù)收益、碳交易收益。這里要注意投資成本是一次性成本但收益是逐時(shí)段的所以需要用等年值法或折現(xiàn)率把全壽命周期的現(xiàn)金流折算到同一年上去。常見做法是使用等年值系數(shù)把初始投資成本折算成年值再與年運(yùn)行收益相加。2.2 下層調(diào)度約束充放電功率、SOC與生命周期下層是運(yùn)行層決策變量是每個(gè)調(diào)度時(shí)段 (t) 的充放電功率、是否充電的0-1狀態(tài)變量、以及各時(shí)段與電網(wǎng)交互的購(gòu)售電功率。約束條件比上層復(fù)雜得多主要包括儲(chǔ)能SOC遞推約束[ SOC_{t1} SOC_t (P_{ch,t} \cdot \eta_{ch} - \frac{P_{dis,t}}{\eta_{dis}}) \cdot \Delta t / E_{ess}^{rated} ]充放電功率上下限[ 0 \le P_{ch,t} \le P_{ess}^{rated} \cdot u_{ch,t}, \quad 0 \le P_{dis,t} \le P_{ess}^{rated} \cdot u_{dis,t} ] 其中 (u_{ch,t} u_{dis,t} \le 1)保證同一時(shí)段不能同時(shí)充放電。SOC上下限[ SOC_{min} \le SOC_t \le SOC_{max} ]調(diào)度周期始末SOC相等保證儲(chǔ)能可持續(xù)運(yùn)行。各新能源場(chǎng)站/用戶的功率平衡約束。這里有個(gè)工程經(jīng)驗(yàn)要提SOC的上下限不建議直接取0和1鋰離子電池的實(shí)際運(yùn)行區(qū)間一般取0.1~0.9。這個(gè)區(qū)間設(shè)置直接影響電池循環(huán)壽命也讓調(diào)度模型更接近真實(shí)。稍微復(fù)雜一點(diǎn)的模型還會(huì)把SOC分段線性化用來(lái)估算電池健康狀態(tài)衰減但如果只是做配置層面的優(yōu)化固定效率固定SOC區(qū)間的模型已經(jīng)夠用。2.3 雙層耦合關(guān)系投資決策如何影響運(yùn)行邊界雙層模型的核心難點(diǎn)在于上下層之間的耦合關(guān)系。上層選定的 (P_{ess}^{rated}) 和 (E_{ess}^{rated}) 會(huì)直接進(jìn)入下層約束作為充放電功率和SOC遞推公式的邊界參數(shù)而下層的調(diào)度結(jié)果——各時(shí)段的充放電量、碳交易量、與電網(wǎng)交互功率曲線——又會(huì)反過來(lái)決定上層目標(biāo)函數(shù)里的運(yùn)行收益項(xiàng)。這種雙向耦合導(dǎo)致沒法把兩層各自單獨(dú)求解必須通過迭代或數(shù)學(xué)轉(zhuǎn)化手段統(tǒng)一求解。在Matlab實(shí)現(xiàn)中處理雙層模型有三種常見策略第一種是用KKT條件把下層模型轉(zhuǎn)化為約束后并入上層第二種是采用粒子群、遺傳算法等啟發(fā)式算法做外層尋優(yōu)內(nèi)層調(diào)用線性規(guī)劃或二次規(guī)劃求解調(diào)度第三種是如果模型規(guī)模不大直接枚舉容量候選集用配置-調(diào)度-評(píng)估的循環(huán)暴力搜索。實(shí)際代碼實(shí)現(xiàn)中第二種最多因?yàn)镵KT轉(zhuǎn)化對(duì)非凸問題很容易出問題而暴力搜索在網(wǎng)格粒度稍細(xì)一點(diǎn)時(shí)就慢到懷疑人生。后面我會(huì)專門講為什么我最終選了啟發(fā)式規(guī)劃器嵌套的方案。3. 碳交易機(jī)制的數(shù)學(xué)建模配額、碳價(jià)與收益3.1 碳交易的核算口徑碳交易機(jī)制要在優(yōu)化模型里落地首先要確定核算口徑。共享儲(chǔ)能場(chǎng)景下我建議把整個(gè)園區(qū)或微網(wǎng)視為一個(gè)碳排放主體系統(tǒng)的碳排放主要來(lái)自從電網(wǎng)購(gòu)電所對(duì)應(yīng)的火電發(fā)電排放。設(shè)電網(wǎng)的綜合碳排放因子為 (EF_{grid})單位是 kgCO2/kWh那么t時(shí)段系統(tǒng)從電網(wǎng)購(gòu)電 (P_{buy,t}) 對(duì)應(yīng)的碳排放量為[ E_{emissions,t} EF_{grid} \cdot P_{buy,t} \cdot \Delta t ]系統(tǒng)獲得的初始配額 (E_{allowance}) 可以按年度總用電量的一定比例設(shè)定也可以與新能源發(fā)電量掛鉤[ E_{allowance} \gamma \cdot (E_{load}^{total} E_{ess}^{charges}) ]其中 (E_{load}^{total}) 是用戶總負(fù)荷電量(E_{ess}^{charges}) 是儲(chǔ)能充電電量(\gamma) 是配額系數(shù)。實(shí)際碳排放量低于配額時(shí)剩余配額可以在碳市場(chǎng)出售高于配額時(shí)不足部分需要購(gòu)買。這樣碳成本/收益項(xiàng)就是[ C_{carbon} P_{carbon} \cdot (E_{emissions,total} - E_{allowance}) ]當(dāng) (E_{emissions,total} E_{allowance}) 時(shí) (C_{carbon} 0)也就是收益。3.2 碳價(jià)與參數(shù)靈敏度碳價(jià)的設(shè)置對(duì)配置結(jié)果影響很大?，F(xiàn)階段國(guó)內(nèi)碳市場(chǎng)成交價(jià)波動(dòng)區(qū)間比較大從每噸幾十元到上百元都有。文獻(xiàn)里常用的是50~100元/噸的區(qū)間。如果你做敏感性分析可以從30元/噸掃到150元/噸看儲(chǔ)能配置容量和總收益的變化趨勢(shì)。以下是一組算例參數(shù)的參考取值我自己跑下來(lái)結(jié)果比較穩(wěn)定參數(shù)取值說(shuō)明電網(wǎng)碳排放因子 (EF_{grid})0.581 kg/kWh區(qū)域電網(wǎng)平均排放因子初始配額系數(shù) (\gamma)0.85按總用電量的85%免費(fèi)發(fā)放配額碳價(jià) (P_{carbon})80 元/噸可按敏感性分析調(diào)整儲(chǔ)能效率 (\eta_{ch}/\eta_{dis})0.95 / 0.95磷酸鐵鋰典型值儲(chǔ)能單位功率成本800 元/kW含PCS等設(shè)備儲(chǔ)能單位容量成本1200 元/kWh含電池組與BMS這里的配額系數(shù)0.85是個(gè)關(guān)鍵假設(shè)。如果取1.0那就是完全免費(fèi)配額儲(chǔ)能減少購(gòu)電反而不會(huì)產(chǎn)生額外碳收益碳交易機(jī)制對(duì)優(yōu)化結(jié)果的影響就弱很多。取0.8以下懲罰意味太強(qiáng)可能導(dǎo)致模型過度激勵(lì)儲(chǔ)能擴(kuò)容。建議在這個(gè)參數(shù)上做敏感性分析而不是拍腦袋固定一個(gè)數(shù)。3.3 碳交易項(xiàng)在目標(biāo)函數(shù)里的工程處理很多第一次建碳模型的同學(xué)會(huì)把碳排放約束寫成硬約束例如系統(tǒng)碳排放不得超過配額X噸。實(shí)際求解時(shí)這種硬約束經(jīng)常導(dǎo)致模型無(wú)解或可行域極小尤其當(dāng)新能源出力差、必須從電網(wǎng)大量購(gòu)電時(shí)。更好的做法是把碳交易項(xiàng)寫成目標(biāo)函數(shù)中的軟性成本項(xiàng)讓模型自動(dòng)權(quán)衡購(gòu)買配額和減少購(gòu)電/增加儲(chǔ)能哪個(gè)更經(jīng)濟(jì)。這也是商業(yè)上真實(shí)的決策邏輯——碳成本只是眾多成本中的一項(xiàng)沒有哪個(gè)企業(yè)會(huì)為了碳達(dá)標(biāo)直接把生產(chǎn)線停了。4. 電網(wǎng)交互波動(dòng)懲罰的刻畫方式與參數(shù)討論4.1 波動(dòng)懲罰的三種常見數(shù)學(xué)形式電網(wǎng)交互波動(dòng)懲罰是標(biāo)題里比較有技術(shù)含量的一部分。先定義變量(P_{grid,t}) 表示t時(shí)段儲(chǔ)能電站與電網(wǎng)的交互功率正值表示從電網(wǎng)購(gòu)電負(fù)值表示向電網(wǎng)送電。波動(dòng)懲罰要抑制的是相鄰時(shí)段交互功率的劇烈變化也就是[ \Delta P_t |P_{grid,t} - P_{grid,t-1}| ]這個(gè) (\Delta P_t) 的懲罰項(xiàng)有三種常見的寫法我列個(gè)對(duì)比形式數(shù)學(xué)表達(dá)式優(yōu)點(diǎn)缺點(diǎn)線性懲罰(C_{fluc} \lambda \cdot \sum |P_{grid,t} - P_{grid,t-1}|)線性模型好求解Yalmip直接用懲罰力度偏軟可能出現(xiàn)小幅高頻波動(dòng)二次懲罰(C_{fluc} \lambda \cdot \sum (P_{grid,t} - P_{grid,t-1})^2)對(duì)大波動(dòng)懲罰更重曲線更平滑引入二次項(xiàng)MILP變成MIQP閾值懲罰(C_{fluc} \lambda \cdot \sum \max(0, |P_{grid,t} - P_{grid,t-1}| - \Delta P_{max}))與實(shí)際并網(wǎng)考核規(guī)則一致需要引入輔助變量模型復(fù)雜度增加我實(shí)際項(xiàng)目中最常用的是線性懲罰原因很簡(jiǎn)單能夠在保持MILP結(jié)構(gòu)的前提下把懲罰項(xiàng)通過引入輔助變量 (u_t^, u_t^-) 拆成兩個(gè)正偏差變量[ P_{grid,t} - P_{grid,t-1} u_t^ - u_t^-, \quad u_t^, u_t^- \ge 0 ][ |P_{grid,t} - P_{grid,t-1}| u_t^ u_t^- ]這樣波動(dòng)懲罰項(xiàng)就變成了線性約束Gurobi或Cplex求解時(shí)不會(huì)破壞MILP結(jié)構(gòu)。如果只是單純?cè)谀繕?biāo)函數(shù)里加絕對(duì)值不引入偏差變量Yalmip里也能通過abs()函數(shù)直接寫但個(gè)人建議還是顯式寫出輔助變量和約束方便調(diào)試也方便看拉格朗日乘子和對(duì)偶變量的值。4.2 懲罰因子選多大才合理懲罰因子 (\lambda_{fluc}) 的取值沒有統(tǒng)一標(biāo)準(zhǔn)它本質(zhì)上是一個(gè)權(quán)重用來(lái)權(quán)衡平滑曲線和犧牲套利收益兩者之間的關(guān)系。我用過的經(jīng)驗(yàn)范圍是(\lambda_{fluc}) 取值為峰谷電價(jià)差單價(jià)的10%~50%。舉個(gè)例子如果峰谷套利的電價(jià)差是0.7元/kWh那么 (\lambda_{fluc}) 設(shè)在0.07~0.35元/kWh之間比較合理。設(shè)太小懲罰項(xiàng)可有可無(wú)儲(chǔ)能調(diào)度曲線基本跟不設(shè)懲罰一樣可能在午間光伏大發(fā)時(shí)段出現(xiàn)功率大起大落設(shè)太大模型會(huì)過度平滑導(dǎo)致儲(chǔ)能幾乎不在電價(jià)尖峰時(shí)段放電配置容量也跟著縮水。最終目標(biāo)應(yīng)該是保證調(diào)度曲線平滑度明顯提升同時(shí)總收益下降不超過5%~10%在這個(gè)區(qū)間內(nèi)取一個(gè)平衡點(diǎn)。這里有一個(gè)快速檢驗(yàn)方法先跑一遍不帶波動(dòng)懲罰的模型統(tǒng)計(jì)相鄰時(shí)段交互功率差值的平均絕對(duì)值再跑帶懲罰的模型如果平均波動(dòng)幅度下降了30%以上且收益損失在可接受范圍內(nèi)說(shuō)明懲罰系數(shù)取得合適。4.3 波動(dòng)懲罰項(xiàng)對(duì)配置結(jié)果的影響邏輯波動(dòng)懲罰不僅影響調(diào)度曲線還會(huì)反過來(lái)影響容量配置。這個(gè)傳導(dǎo)邏輯值得展開說(shuō)一下。上層在做容量?jī)?yōu)化時(shí)它要評(píng)估多配一組儲(chǔ)能電池能不能帶來(lái)額外收益。沒有波動(dòng)懲罰時(shí)新增容量可以在電價(jià)低谷大量充電、高峰大量放電邊際收益明顯但有波動(dòng)懲罰后儲(chǔ)能在高峰放電后不能立刻在低谷充電充放電之間要有爬坡過渡這就導(dǎo)致儲(chǔ)能利用率下降等效循環(huán)次數(shù)減少邊際收益也跟著下降。所以在高波動(dòng)懲罰系數(shù)下上層模型最終選出的額定容量往往會(huì)比無(wú)懲罰時(shí)低10%~20%同時(shí)最優(yōu)功率容量比即儲(chǔ)能時(shí)長(zhǎng)會(huì)略微拉長(zhǎng)。這個(gè)趨勢(shì)如果沒在論文里討論到審稿人大概率會(huì)問。實(shí)際復(fù)現(xiàn)時(shí)你應(yīng)該畫一張不同懲罰因子下最優(yōu)配置結(jié)果對(duì)比表非常直觀也很有說(shuō)服力。5. Matlab代碼實(shí)現(xiàn)路線從目標(biāo)函數(shù)到求解器調(diào)用5.1 建模前的數(shù)據(jù)準(zhǔn)備Matlab實(shí)現(xiàn)的第一步不是寫代碼而是把數(shù)據(jù)整理成能被Yalmip識(shí)別的格式。我們需要的典型數(shù)據(jù)文件包括負(fù)荷數(shù)據(jù)24h或8760h的用戶側(cè)負(fù)荷曲線單位kW新能源出力數(shù)據(jù)風(fēng)電/光伏的歸一化出力序列分時(shí)電價(jià)表峰、平、谷時(shí)段的購(gòu)電價(jià)和售電價(jià)電網(wǎng)交互功率的初始上下限值碳交易參數(shù)配額系數(shù)、碳價(jià)、排放因子儲(chǔ)能技術(shù)參數(shù)充放電效率、SOC范圍、成本系數(shù)這里最花時(shí)間的往往是新能源出力數(shù)據(jù)的清洗。實(shí)際數(shù)據(jù)里經(jīng)常有缺失值、負(fù)值異常和零值密集段。我在第一次跑這個(gè)模型時(shí)直接用原始光伏出力數(shù)據(jù)結(jié)果調(diào)度模型把儲(chǔ)能安排得亂七八糟——因?yàn)榱璩繒r(shí)段偶爾冒出幾個(gè)尖峰毛刺模型以為是可充電源。后來(lái)做了移動(dòng)平均平滑和限幅處理結(jié)果才正常。建議所有時(shí)序數(shù)據(jù)在進(jìn)入模型前先畫圖看一眼確認(rèn)趨勢(shì)合理再做歸一化。5.2 用Yalmip寫上層配置-下層調(diào)度循環(huán)的骨架下面這段是我常用的共享儲(chǔ)能優(yōu)化配置代碼框架YalmipGurobi組合。這里只展示核心結(jié)構(gòu)重點(diǎn)看雙層迭代和懲罰項(xiàng)是怎么落地的。%% 參數(shù)初始化 T 24; % 調(diào)度時(shí)段數(shù) N_pop 30; % 粒子群種群數(shù)量 max_iter 50; % 外層迭代次數(shù) % 儲(chǔ)能成本參數(shù) c_power 800; % 單位功率成本 元/kW c_energy 1200; % 單位容量成本 元/kWh life_year 10; % 壽命年限 r_discount 0.08; % 折現(xiàn)率 % 電網(wǎng)交互波動(dòng)懲罰系數(shù) lambda_fluc 0.15; % 元/kWh % 分時(shí)電價(jià)峰平谷三段 price_buy [0.35*ones(T/3,1); 0.8*ones(T/3,1); 1.2*ones(T/3,1)]; price_sell price_buy * 0.85; %% 外層粒子群算法優(yōu)化上層配置 for iter 1:max_iter for k 1:N_pop % 從種群中取出候選配置 P_rated pop_power(k); % 候選額定功率 kW E_rated pop_energy(k); % 候選額定容量 kWh % 內(nèi)層求解當(dāng)前配置下的最優(yōu)調(diào)度MILP [revenue, P_grid] solve_dispatch(P_rated, E_rated, ... load_data, pv_data, price_buy, price_sell, ... lambda_fluc, carbon_params); % 計(jì)算上層凈收益 invest_cost c_power * P_rated c_energy * E_rated; annual_invest invest_cost * (r_discount*(1r_discount)^life_year) ... / ((1r_discount)^life_year - 1); fitness(k) revenue - annual_invest; end % 更新粒子群位置、速度此處省略標(biāo)準(zhǔn)PSO更新代碼 end內(nèi)層調(diào)度函數(shù)solve_dispatch用Yalmip建模關(guān)鍵片段如下function [revenue, P_grid] solve_dispatch(P_rated, E_rated, ...) %% 定義決策變量 P_ch sdpvar(T, 1); % 充電功率 P_dis sdpvar(T, 1); % 放電功率 u_ch binvar(T, 1); % 充電狀態(tài) u_dis binvar(T, 1); % 放電狀態(tài) SOC sdpvar(T1, 1); % SOC狀態(tài)變量 P_buy sdpvar(T, 1); % 從電網(wǎng)購(gòu)電功率 P_sell sdpvar(T, 1); % 向電網(wǎng)售電功率 delta_up sdpvar(T-1, 1); % 交互功率正向波動(dòng) delta_down sdpvar(T-1, 1); % 交互功率負(fù)向波動(dòng) %% 約束條件 C []; C [C, 0 P_ch P_rated * u_ch]; C [C, 0 P_dis P_rated * u_dis]; C [C, u_ch u_dis 1]; SOC(1) 0.2 * E_rated; % 初始SOC C [C, 0.1*E_rated SOC 0.9*E_rated]; for t 1:T C [C, SOC(t1) SOC(t) (P_ch(t)*0.95 - P_dis(t)/0.95)]; end C [C, SOC(T1) SOC(1)]; % 周期始末SOC一致 %% 電網(wǎng)交互功率與波動(dòng)懲罰約束 P_grid P_buy - P_sell; C [C, 0 P_buy 5000, 0 P_sell 5000]; for t 2:T % 相鄰時(shí)段交互功率差值拆成正負(fù)偏差 C [C, P_grid(t) - P_grid(t-1) delta_up(t-1) - delta_down(t-1)]; C [C, delta_up(t-1) 0, delta_down(t-1) 0]; end %% 目標(biāo)函數(shù) revenue_elec sum(P_dis .* price_sell) - sum(P_ch .* price_buy); carbon_cost carbon_params.carbon_price * ... (sum(P_buy) * carbon_params.grid_ef - carbon_params.allowance); fluctuation_cost lambda_fluc * sum(delta_up delta_down); objective revenue_elec - carbon_cost - fluctuation_cost; %% 求解 ops sdpsettings(solver, gurobi, verbose, 0); optimize(C, -objective, ops); end這段代碼里我要特別說(shuō)明兩個(gè)設(shè)計(jì)。第一SOC遞推約束我直接用E_rated做了歸一化這樣SOC變量保持在0~1區(qū)間數(shù)值穩(wěn)定性好第二電網(wǎng)交互功率我這里用P_buy和P_sell兩個(gè)非負(fù)變量相減來(lái)表示這比用一個(gè)有正有負(fù)的變量更規(guī)范因?yàn)楹罄m(xù)要區(qū)分購(gòu)電和售電價(jià)格。5.3 為什么不直接做KKT轉(zhuǎn)化很多人問為什么不用KKT條件把雙層模型轉(zhuǎn)化成單層MILP一起求解這個(gè)問題我實(shí)踐之后的體會(huì)是KKT轉(zhuǎn)化在理論上完美但工程實(shí)現(xiàn)時(shí)非常痛苦。下層如果包含二元變量KKT轉(zhuǎn)化的強(qiáng)對(duì)偶條件基本沒法直接寫就算下層全是連續(xù)變量轉(zhuǎn)化后的互補(bǔ)松弛約束會(huì)引入大量非線性項(xiàng)還得做大M法線性化而大M的取值又直接影響求解穩(wěn)定性。相比之下粒子群/遺傳算法外層尋優(yōu)K層內(nèi)MILP精確求解的組合更貼近工程實(shí)踐。內(nèi)層MILP本身是凸問題Gurobi可以在幾十毫秒到幾秒內(nèi)求解完成外層種群設(shè)30~50個(gè)個(gè)體迭代50次也就是求解1500~2500次MILP總耗時(shí)在幾分鐘到一小時(shí)內(nèi)。對(duì)于配置規(guī)劃這類離線計(jì)算場(chǎng)景這個(gè)速度完全可接受而且代碼可控性高出了問題容易定位。6. 跑通算例之后我踩過的幾個(gè)坑6.1 非線性項(xiàng)處理的第一反應(yīng)第一次寫目標(biāo)函數(shù)時(shí)我順手把 (\Delta P_t) 的平方項(xiàng)直接寫進(jìn)Yalmip想著Gurobi能處理MIQP。結(jié)果在SOC區(qū)間約束和二進(jìn)制狀態(tài)變量疊加下求解時(shí)間從幾秒飆到幾百秒甚至頻繁報(bào)數(shù)值問題。后來(lái)?yè)Q成線性懲罰輔助變量拆分求解時(shí)間降回原來(lái)的水平結(jié)果曲線也很平滑。我的建議是除非二次型對(duì)結(jié)果有不可替代的意義否則優(yōu)先選線性懲罰尤其在做大規(guī)模長(zhǎng)時(shí)間序列8760h場(chǎng)景時(shí)線性化是必須的。6.2 雙層迭代不收斂問題出在離散噪聲外層用PSO時(shí)經(jīng)常遇到迭代曲線鋸齒狀劇烈波動(dòng)的現(xiàn)象。剛開始我以為算法參數(shù)沒調(diào)好后來(lái)分析發(fā)現(xiàn)問題不在PSO而在內(nèi)層調(diào)度問題對(duì)配置參數(shù)不連續(xù)。當(dāng)額定功率從1200kW變成1210kW時(shí)內(nèi)層MILP的最優(yōu)調(diào)度可能跳變到完全不同的運(yùn)行模式比如從每天兩充兩放跳到深度一充一放導(dǎo)致目標(biāo)函數(shù)值突變。這種離散噪聲會(huì)讓標(biāo)準(zhǔn)PSO收斂很慢。我用的處理辦法是在適應(yīng)度評(píng)估里加一個(gè)小懲罰如果候選配置相比上一代全局最優(yōu)的容量變化幅度超過30%在適應(yīng)度函數(shù)里加一個(gè)探索懲罰項(xiàng)抑制粒子盲目跳躍。另外把粒子群粒子的速度上限設(shè)小一檔比如速度限幅設(shè)為搜索空間寬度的10%收斂穩(wěn)定性明顯改善。6.3 結(jié)果合理性檢驗(yàn)三個(gè)必須看的指標(biāo)模型跑完不要直接拿結(jié)果去寫論文先做三件事第一繪制儲(chǔ)能的充放電功率曲線檢查是否頻繁在相鄰時(shí)段出現(xiàn)充電-放電-充電的鋸齒模式。如果出現(xiàn)說(shuō)明或者電價(jià)數(shù)據(jù)有問題或者輔助服務(wù)收益項(xiàng)界定不清真實(shí)系統(tǒng)里不會(huì)有人這樣操作電池。第二檢查儲(chǔ)能年循環(huán)次數(shù)。按我的經(jīng)驗(yàn)日循環(huán)次數(shù)一般在0.8~1.5次之間折合成年循環(huán)次數(shù)為300~500次。如果你算出來(lái)的年循環(huán)次數(shù)超過800次說(shuō)明模型在過度使用電池全壽命周期成本核算一定有問題實(shí)際電池早衰減了。第三做懲罰系數(shù)靈敏度測(cè)試。固定其他參數(shù)把 (\lambda_{fluc}) 從0逐步增加到目標(biāo)值看最優(yōu)配置容量是否單調(diào)不增。如果出現(xiàn)容量忽大忽小的非單調(diào)現(xiàn)象多半是雙層嵌套求解時(shí)陷入了局部最優(yōu)需要增大PSO種群規(guī)?；蛘{(diào)整參數(shù)范圍重新搜索。6.4 關(guān)于Matlab版本和求解器的幾點(diǎn)備注YalmipGurobi的組合在Matlab環(huán)境里是跑這類模型的主力方案。我自己在用的Matlab版本是R2023a之后的版本Yalmip選的是GitHub上持續(xù)維護(hù)的版本Gurobi用的是學(xué)術(shù)授權(quán)。需要提醒的是幾個(gè)環(huán)境兼容問題Yalmip對(duì)較新的Matlab版本偶爾會(huì)有適配滯后如果遇到求解器識(shí)別不了先檢查Yalmip是否最新。Gurobi的Matlab接口在Windows和Linux下的安裝路徑不一樣Linux下建議配置環(huán)境變量GUROBI_HOME后在startup.m里手動(dòng)添加路徑。如果機(jī)器內(nèi)存有限8760h全年場(chǎng)景建議把SOC約束寫成矩陣方式一次性添加不要用for循環(huán)逐條[C, C [C, ...]]拼接否則構(gòu)建約束的時(shí)間比求解還長(zhǎng)。這些坑大多不會(huì)出現(xiàn)在論文的模型描述里但復(fù)現(xiàn)過的人應(yīng)該都懂?；氐轿易约旱膽?yīng)用感受。共享儲(chǔ)能配置與調(diào)度模型最迷人的地方不是單點(diǎn)技術(shù)多高深而是碳交易、電網(wǎng)波動(dòng)懲罰、容量規(guī)劃和時(shí)序調(diào)度這幾件原本分散在不同體系里的事情能統(tǒng)一收攏到一個(gè)目標(biāo)函數(shù)下互相權(quán)衡。把Matlab代碼跑通只是第一步真正有價(jià)值的是在跑算例的過程中理解每個(gè)參數(shù)背后的工程含義——配額的松緊如何影響儲(chǔ)能投資意愿波動(dòng)懲罰的權(quán)重如何改變調(diào)度行為這些結(jié)論在同質(zhì)化嚴(yán)重的儲(chǔ)能論文里反而是最能讓工作被記住的部分。希望這篇文章能幫你少走一些彎路。