系統(tǒng)有功-無功協(xié)同優(yōu)化：碳約束下的建模與求解)
1. “碳中和”目標(biāo)下的電氣互聯(lián)系統(tǒng)為什么非做有功-無功協(xié)同不可先說一個(gè)我在做這類項(xiàng)目時(shí)最深的體會(huì)很多做電力系統(tǒng)優(yōu)化的同學(xué)把重點(diǎn)全放在有功調(diào)度上無功這塊要么忽略要么用固定的功率因數(shù)折算一下。但在“碳中和”目標(biāo)下電氣互聯(lián)系統(tǒng)的優(yōu)化邊界已經(jīng)被徹底拓寬了——你面對(duì)的不再是單一電網(wǎng)而是電力系統(tǒng)與天然氣系統(tǒng)深度耦合的網(wǎng)絡(luò)負(fù)荷側(cè)多了電轉(zhuǎn)氣、氫能等元素電源側(cè)多了大量分布式可再生能源。此時(shí)如果還把有功和無功拆開做結(jié)果不是次優(yōu)而是可能直接得出不可行的運(yùn)行方案。所謂電氣互聯(lián)系統(tǒng)通俗講就是電網(wǎng)和氣網(wǎng)通過兩類關(guān)鍵元件“綁”在了一起一類是燃?xì)鈾C(jī)組把天然氣轉(zhuǎn)化為電力另一類是電轉(zhuǎn)氣設(shè)備Power to GasP2G把富余電能轉(zhuǎn)化為天然氣。前者是天然氣系統(tǒng)向電力系統(tǒng)供電后者是電力系統(tǒng)反過來向天然氣系統(tǒng)供氣。這種雙向耦合的直接后果就是電網(wǎng)的有功潮流變化會(huì)改變氣網(wǎng)的供氣需求氣網(wǎng)的供氣壓降又反過來約束燃?xì)鈾C(jī)組的出力上限而機(jī)組無功出力調(diào)整又會(huì)改變線路無功潮流和電壓分布——環(huán)環(huán)相扣誰也沒法獨(dú)立優(yōu)化。為什么“碳中和”目標(biāo)讓這個(gè)問題更棘手因?yàn)楦弑壤稍偕茉唇尤牒笙到y(tǒng)慣量降低、電壓支撐能力變?nèi)鯚o功資源變得更加稀缺。光伏和風(fēng)電在滿發(fā)時(shí)可能擠占常規(guī)機(jī)組的出清空間但它們的無功調(diào)節(jié)能力有限甚至某些工況下還要吸收無功。另一個(gè)問題是碳排放約束把氣網(wǎng)側(cè)的成本結(jié)構(gòu)改變了——碳價(jià)、碳配額都進(jìn)入了目標(biāo)函數(shù)這就讓“電轉(zhuǎn)氣、氣轉(zhuǎn)電”的決策路徑變得更加復(fù)雜。我做這個(gè)項(xiàng)目的核心思路可以概括成一句話在滿足碳排放約束的前提下同時(shí)調(diào)節(jié)機(jī)組有功出力、無功出力、變壓器分接頭、無功補(bǔ)償裝置投切以及氣源供氣量讓整個(gè)電氣互聯(lián)系統(tǒng)的綜合運(yùn)行成本最低。這個(gè)思路拆解下來涉及三個(gè)核心問題如何建立電網(wǎng)和氣網(wǎng)統(tǒng)一的數(shù)學(xué)模型特別是描述耦合元件的變量關(guān)系。如何處理碳排放約束和目標(biāo)函數(shù)中碳成本的量化。如何在Matlab中高效求解這樣一個(gè)大規(guī)模非線性混合整數(shù)優(yōu)化問題。下面我會(huì)沿著這三條線把從建模到求解再到代碼實(shí)現(xiàn)的全過程完整分享出來。這篇文章適合正在做綜合能源系統(tǒng)優(yōu)化、無功優(yōu)化或者碳約束調(diào)度的研究生和工程師參考也適合想用Matlab搭建一體化優(yōu)化框架的讀者當(dāng)模板。2. 有功-無功協(xié)同的耦合機(jī)理數(shù)學(xué)上它們是怎么“糾纏”在一起的很多人覺得“有功-無功協(xié)同”就是把兩個(gè)目標(biāo)函數(shù)加到一起然后在約束里同時(shí)寫上機(jī)組有功上下限和無功上下限。這是典型的誤區(qū)。真正的耦合體現(xiàn)在等式約束和不等式約束交叉關(guān)聯(lián)上下面我把數(shù)學(xué)模型層面的耦合點(diǎn)逐一拆開。2.1 電網(wǎng)側(cè)數(shù)學(xué)模型潮流方程中的有功無功耦合電網(wǎng)側(cè)采用交流潮流模型時(shí)節(jié)點(diǎn)i的功率平衡方程是P_i U_i * Σ(U_j * (G_ij * cos(θ_ij) B_ij * sin(θ_ij)))Q_i U_i * Σ(U_j * (G_ij * sin(θ_ij) - B_ij * cos(θ_ij)))從這個(gè)式子就能看出有功P_i和無功Q_i共享同一組電壓幅值U和相角θ變量。改變?nèi)我庖粋€(gè)節(jié)點(diǎn)的無功全網(wǎng)的電壓分布都會(huì)變電壓變了相角也跟著變相角變了有功潮流就變了。這就是電網(wǎng)側(cè)最本質(zhì)的耦合。在工程簡化中有人用直流潮流模型只保留有功維度這在輸電網(wǎng)經(jīng)濟(jì)調(diào)度里勉強(qiáng)可用但一旦涉及電壓安全約束、無功補(bǔ)償優(yōu)化、變壓器分接頭調(diào)節(jié)直流潮流就完全失效了。本模型必須用交流潮流。對(duì)于輸電線路還有一個(gè)關(guān)鍵約束是線路潮流極限。實(shí)際傳輸容量受熱穩(wěn)定極限、電壓約束和穩(wěn)定極限三重限制用視在功率S_ij的平方表示就是S_ij^2 P_ij^2 Q_ij^2 ≤ (S_ij^max)^2這個(gè)約束本身就是有功與無功的非線性耦合約束。如果只優(yōu)化有功而忽略無功線路可能會(huì)在視在功率上越限而單純看有功卻一切正常。2.2 天然氣網(wǎng)絡(luò)模型氣體流動(dòng)方程與節(jié)點(diǎn)氣壓約束天然氣網(wǎng)絡(luò)的穩(wěn)態(tài)模型常用節(jié)點(diǎn)氣流平衡方程描述。對(duì)氣網(wǎng)節(jié)點(diǎn)n氣流平衡可以寫成Σ(FLOW_k) G_source,m - G_load,m - G_p2g,m 0其中FLOW_k表示管道k的氣流量正方向?yàn)榱魅牍?jié)點(diǎn)反方向?yàn)榱鞒鯣_source是氣源供氣量G_load是氣負(fù)荷包括燃?xì)鈾C(jī)組耗氣量和民用氣負(fù)荷G_p2g是P2G設(shè)備消耗的電量轉(zhuǎn)化成的天然氣注入量。管道氣流FLOW與兩端氣壓的關(guān)系通常用Weymouth穩(wěn)態(tài)方程FLOW_ij^2 K_ij^2 * (p_i^2 - p_j^2)這個(gè)式子是非線性的而且是本模型的第一個(gè)非線性難點(diǎn)來源。注意這里FLOW是有方向的實(shí)際工程中需要處理方向符號(hào)問題常見做法是引入雙向流標(biāo)志位或者直接假設(shè)氣流方向已知——在做日前優(yōu)化調(diào)度時(shí)氣流方向通常可以從潮流結(jié)果預(yù)判但嚴(yán)謹(jǐn)起見最好在約束中保留方向變量。氣壓還有上下限約束p_min ≤ p_i ≤ p_max氣源出力也有爬坡限制G_source,min ≤ G_source ≤ G_source,max這里有個(gè)值得注意的細(xì)節(jié)氣網(wǎng)的慣性比電網(wǎng)大得多但穩(wěn)態(tài)模型下管道存氣linepack效應(yīng)被忽略了。在日內(nèi)滾動(dòng)優(yōu)化中l(wèi)inepack往往能提供寶貴的短時(shí)緩沖能力但本文的日前模型先不考慮保持模型相對(duì)簡潔。2.3 耦合元件模型燃?xì)鈾C(jī)組與P2G的雙向耦合燃?xì)鈾C(jī)組是氣網(wǎng)到電網(wǎng)的耦合元件其模型包括兩部分注入電網(wǎng)的有功出力P_g和旋轉(zhuǎn)備用容量。從氣網(wǎng)消耗的天然氣量F_gas a * P_g^2 b * P_g c其中a、b、c是燃料消耗系數(shù)。這個(gè)二次函數(shù)也是典型的非線性環(huán)節(jié)。另一個(gè)方向——P2G設(shè)備的耦合模型為G_p2g η_p2g * P_p2g / GHV其中P_p2g是P2G設(shè)備消耗的有功功率η_p2g是轉(zhuǎn)換效率GHV是天然氣高位熱值。P2G消耗的功率上限同時(shí)受電網(wǎng)上網(wǎng)容量和氣網(wǎng)注入能力的雙重約束P_p2g,min ≤ P_p2g ≤ min(P_p2g,max, G_p2g,max * GHV / η_p2g)到這里有功-無功耦合的數(shù)學(xué)畫像就清晰了電網(wǎng)側(cè)P、Q通過電壓幅值和相角在潮流方程中耦合。氣網(wǎng)側(cè)氣流量通過壓力非線性方程耦合。系統(tǒng)級(jí)燃?xì)鈾C(jī)組和P2G在有功功率與氣流量之間建立雙向通道而這種耦合通過上述兩個(gè)非線性環(huán)節(jié)傳遞。抽象成一句話有功決策變了無功優(yōu)化空間也跟著變無功決策變了電壓分布變了P2G消耗的有功上限和燃?xì)鈾C(jī)組的可用容量也變了。所謂“協(xié)同”就是在數(shù)學(xué)上同時(shí)處理這兩組變量而不是把它們解耦成兩個(gè)子問題分別求解后再拼起來。3. 協(xié)同優(yōu)化模型搭建目標(biāo)函數(shù)與關(guān)鍵約束的設(shè)計(jì)思路我在這類項(xiàng)目上走過的彎路不少最值錢的一條經(jīng)驗(yàn)是目標(biāo)函數(shù)決定了優(yōu)化方向約束決定可行域邊界但真正決定模型能否被高效求解的是約束的“數(shù)學(xué)形態(tài)”是線性、凸還是非凸。3.1 目標(biāo)函數(shù)設(shè)計(jì)成本與碳排放的雙重權(quán)衡本模型的目標(biāo)函數(shù)我在實(shí)際項(xiàng)目中是這樣設(shè)計(jì)的綜合考慮了四種成本min F_total Σ(C_fuel) Σ(C_co2) Σ(C_p2g) Σ(C_grid)第一項(xiàng)是常規(guī)機(jī)組燃料成本是機(jī)組有功出力的二次函數(shù)C_fuel,i a_i * P_i^2 b_i * P_i c_i第二項(xiàng)是燃?xì)鈾C(jī)組的碳排放成本這里引入“碳中和”目標(biāo)的關(guān)鍵我在實(shí)際建模中用了兩種處理方式碳市場交易價(jià)格乘以排放量或者設(shè)定碳配額上限然后用階梯碳價(jià)。模型里我建議使用線性碳價(jià)對(duì)單位碳排放量乘以碳價(jià)系數(shù)C_co2,i ε_(tái)i * e_i * P_i其中ε_(tái)i是碳價(jià)元/tCO2e_i是單位電量的碳排放強(qiáng)度tCO2/MWh。第三項(xiàng)是P2G運(yùn)行成本主要包括設(shè)備運(yùn)維成本和購電成本分?jǐn)侰_p2g γ * P_p2g其中γ是P2G單位運(yùn)行成本系數(shù)。第四項(xiàng)是從上級(jí)電網(wǎng)購電的成本只有當(dāng)系統(tǒng)接入主網(wǎng)時(shí)才需要。這里我想特別說明為什么要這樣設(shè)目標(biāo)函數(shù)而不是簡單加一個(gè)“碳排放最小”目標(biāo)。因?yàn)榧儨p排目標(biāo)在經(jīng)濟(jì)上不可行實(shí)際運(yùn)行中調(diào)度員關(guān)心的是“在碳約束下的最經(jīng)濟(jì)運(yùn)行”所以碳成本進(jìn)入目標(biāo)函數(shù)、碳約束進(jìn)入不等式約束這兩個(gè)維度分開處理才能保證結(jié)果既環(huán)保又可執(zhí)行。3.2 等式約束詳解電網(wǎng)潮流與氣網(wǎng)平衡電場部分需要同時(shí)滿足有功和無功平衡P_G,i - P_L,i - U_i * Σ(U_j * (G_ij * cos θ_ij B_ij * sin θ_ij)) - P_P2G,i 0Q_G,i - Q_L,i - U_i * Σ(U_j * (G_ij * sin θ_ij - B_ij * cos θ_ij)) - Q_P2G,i 0注意我在P2G節(jié)點(diǎn)的有功平衡里專門扣掉了P_P2G,i這是典型的“負(fù)荷型元件”處理方式P2G消耗有功同時(shí)如果是電解水制氫再合成甲烷的路徑還需要向電網(wǎng)提供少量的無功給整流器等設(shè)備。氣網(wǎng)側(cè)的氣流平衡在2.2節(jié)已經(jīng)給出。這里補(bǔ)充說明一個(gè)細(xì)節(jié)對(duì)于含壓縮機(jī)站的天然氣網(wǎng)絡(luò)壓縮機(jī)的耗氣量也是變量但在基礎(chǔ)模型中往往先用固定比例折算C_comp β * FLOW_comp這種簡化在工程上是可接受的因?yàn)閴嚎s機(jī)耗氣量占總輸氣量的比例通常只有1%到3%不會(huì)對(duì)優(yōu)化結(jié)果產(chǎn)生方向性影響。3.3 不等式約束安全與設(shè)備極限的建模技巧不等式約束是我在實(shí)際代碼中掉過坑最多的地方。常見的約束包括發(fā)電機(jī)組有功出力上下限P_min ≤ P ≤ P_max。無功出力上下限Q_min ≤ Q ≤ Q_max這里要注意無功上限與有功出力有關(guān)嚴(yán)格來說應(yīng)該是Q_max(P)但很多模型直接給定常數(shù)上限會(huì)帶來結(jié)果偏樂觀的問題。節(jié)點(diǎn)電壓幅值約束U_min ≤ U ≤ U_max在碳中和場景下電壓約束的重要性會(huì)被放大新能源出力的隨機(jī)性讓電壓越限風(fēng)險(xiǎn)顯著增加。變壓器分接頭調(diào)節(jié)范圍約束。氣網(wǎng)節(jié)點(diǎn)氣壓約束和管道輸氣能力約束。線路傳輸容量約束。在眾多不等式約束里最難處理的是線路容量約束。這個(gè)約束用線性化近似會(huì)丟失無功信息用完整交流潮流又帶來大量非線性。我的做法是保留交流潮流形式的視在功率約束配合松弛技術(shù)確?？山庑?。還有個(gè)容易被忽略的約束是P2G設(shè)備的運(yùn)行約束。P2G啟動(dòng)需要時(shí)間電解槽的爬坡速率有限-P_p2g,ramp ≤ P_p2g,t - P_p2g,t-1 ≤ P_p2g,ramp如果忽略這個(gè)約束優(yōu)化器會(huì)讓P2G在相鄰時(shí)段隨意跳變結(jié)果在實(shí)際中根本執(zhí)行不了——這在電網(wǎng)側(cè)調(diào)度中也是同樣的道理。3.4 碳排放約束的建模與處理碳約束我傾向于直接寫成總量約束的形式Σ(E_G,i) ≤ E_total,max其中E_G,i是第i臺(tái)機(jī)組的碳排放量E_total,max是系統(tǒng)設(shè)定的碳排放上限。在“碳中和”目標(biāo)下這個(gè)上限應(yīng)該逐年收緊。如果進(jìn)一步要求模型能反映不同碳價(jià)情景下的運(yùn)行策略差異還可以把碳約束轉(zhuǎn)化為拉格朗日懲罰項(xiàng)或者KKT乘子。但對(duì)優(yōu)化模型來說顯式碳約束更直觀且方便做邊際排放成本的影子價(jià)格分析——影子價(jià)格本身就是很有用的調(diào)度決策輔助信息。4. MatLab代碼實(shí)現(xiàn)路線從數(shù)學(xué)模型到可運(yùn)行程序模型搭完之后怎么把它在Matlab里落地是另一個(gè)大問題。我分享一下自己的完整技術(shù)路線包括求解工具選型、代碼架構(gòu)和核心函數(shù)的實(shí)現(xiàn)思路。4.1 求解工具與算法選型這個(gè)問題先說結(jié)論我選的是YALMIP建模工具箱 二階錐松弛SOCP 商業(yè)求解器CPLEX或Gurobi的組合方案。為什么這么選電氣互聯(lián)系統(tǒng)優(yōu)化問題本質(zhì)上是MINLP混合整數(shù)非線性規(guī)劃直接求解全局最優(yōu)幾乎不可能。工程上主要有四條路線內(nèi)點(diǎn)法直接求解連續(xù)松弛的NLP速度快但容易陷入局部最優(yōu)且很難處理整數(shù)變量。智能算法如粒子群和遺傳算法不依賴梯度能處理非線性但每次求解的穩(wěn)定性差而且沒有最優(yōu)性保證。Benders分解把原問題分解成主問題和子問題適合大規(guī)模系統(tǒng)但模型推導(dǎo)復(fù)雜代碼實(shí)現(xiàn)成本高。二階錐松弛把非凸的Weymouth方程和交流潮流約束松弛為凸約束配合Big-M法處理離散變量形成MISOCP問題商用求解器對(duì)這類問題的求解能力已經(jīng)非常成熟。我實(shí)測下來二階錐松弛商業(yè)求解器在幾十個(gè)節(jié)點(diǎn)的系統(tǒng)上求解速度在幾十秒到幾分鐘之間兼顧了精度和速度是最適合Matlab實(shí)現(xiàn)的技術(shù)路線。具體做法是把Weymouth方程做二階錐轉(zhuǎn)換定義變量α_i p_i^2任意變換后約束變成下式FLOW_ij^2 ≤ K_ij^2 * (α_i - α_j)同時(shí)把潮流方程中的二次項(xiàng)通過電壓幅值平方替換處理。注意這里的松弛方向是有講究的——必須松弛成小于等于的形式才能保持凸性。雖然松弛擴(kuò)大了可行域但實(shí)際運(yùn)行中優(yōu)化結(jié)果一般會(huì)自動(dòng)逼近等式邊界。4.2 代碼架構(gòu)分模塊設(shè)計(jì)的思路我的Matlab代碼按模塊組織便于調(diào)試和復(fù)用主文件main.m定義系統(tǒng)參數(shù)加載數(shù)據(jù)調(diào)用建模函數(shù)求解輸出結(jié)果。數(shù)據(jù)定義文件case_data.m電力系統(tǒng)節(jié)點(diǎn)參數(shù)、氣網(wǎng)管道參數(shù)、耦合元件參數(shù)等基礎(chǔ)數(shù)據(jù)。建模函數(shù)build_opf_model.m用YALMIP定義優(yōu)化變量、目標(biāo)函數(shù)和約束。結(jié)果分析函數(shù)plot_results.m繪制電壓分布、機(jī)組出力、氣網(wǎng)壓力、碳排量變化等曲線。核心的YALMIP建模框架大概是下面的思路以目標(biāo)函數(shù)和關(guān)鍵約束為例偽代碼% YALMIP變量定義 P_g sdpvar(n_gen, 1); % 發(fā)電機(jī)有功 Q_g sdpvar(n_gen, 1); % 發(fā)電機(jī)無功 U sdpvar(n_bus, 1); % 節(jié)點(diǎn)電壓幅值 theta sdpvar(n_bus, 1); % 節(jié)點(diǎn)相角 P_p2g sdpvar(n_p2g, 1); % P2G有功消耗 G_source sdpvar(n_source, 1); % 氣源供氣量 FLOW sdpvar(n_pipe, 1); % 管道流量 % 目標(biāo)函數(shù) Objective sum(a .* P_g.^2 b .* P_g c) ... sum(carbon_price .* emission_coeff .* P_g) ... sum(gamma .* P_p2g); % 約束定義 Constraints []; % 節(jié)點(diǎn)有功平衡 Constraints [Constraints, A_incidence * P_g P_wind - P_load - P_p2g_map ... U .* (G_ij * cos(theta_ij) B_ij * sin(theta_ij))]; % 節(jié)點(diǎn)無功平衡 Constraints [Constraints, Q_g Q_comp - Q_load ... U .* (G_ij * sin(theta_ij) - B_ij * cos(theta_ij))]; % 燃?xì)鈾C(jī)組氣負(fù)荷 Constraints [Constraints, F_gas a_fuel * P_g.^2 b_fuel .* P_g c_fuel]; % 氣網(wǎng)節(jié)點(diǎn)平衡 Constraints [Constraints, A_gas * FLOW G_source - G_load - F_gas 0]; % 二階梯松弛 for i 1:n_pipe Constraints [Constraints, FLOW(i)^2 K(i)^2 * (alpha_from(i) - alpha_to(i))]; end % 求解 ops sdpsettings(solver,gurobi,verbose,1); optimize(Constraints, Objective, ops);這段代碼里最關(guān)鍵的幾個(gè)細(xì)節(jié)值得說明第一潮流和網(wǎng)絡(luò)拓?fù)湎嚓P(guān)的矩陣G_ij、B_ij不是手工輸入的而是用Matpower的makeYbus函數(shù)自動(dòng)生成的。這樣即使換算例代碼不用改結(jié)構(gòu)。第二P2G映射關(guān)系P_p2g_map的處理是我調(diào)試時(shí)踩過的一個(gè)坑P2G設(shè)備接在哪個(gè)節(jié)點(diǎn)上必須在數(shù)據(jù)文件里明確指定并且要形成從變量索引到網(wǎng)絡(luò)節(jié)點(diǎn)的映射矩陣否則氣網(wǎng)和電網(wǎng)的平衡方程會(huì)錯(cuò)位。第三二階梯約束的寫法里我用的是而不是這關(guān)系到問題是否凸。如果用等式約束YALMIP會(huì)把它當(dāng)作非線性等式來解速度慢幾個(gè)數(shù)量級(jí)還可能因?yàn)槌跏键c(diǎn)差解不出來。4.3 針對(duì)整數(shù)變量的處理變壓器分接頭與無功補(bǔ)償投切變壓器分接頭和無功補(bǔ)償電容/電抗器的投切是離散變量我一開始直接把分接頭位置定義為整數(shù)變量binvar或intvar導(dǎo)致求解時(shí)間暴增。后來改成Big-M法處理先把分接頭比轉(zhuǎn)換為連續(xù)變量t并添加線性約束t_min Σ(step_i * z_i) ≤ t ≤ t_min Σ(step_i * z_i) (1 - z_i) * M其中z_i是0-1變量M是足夠大的數(shù)。這種做法把離散決策從“選位置”轉(zhuǎn)換成“選步進(jìn)”在Gurobi里求解效率高多了。還有一個(gè)實(shí)用技巧如果只是做日前調(diào)度變壓器分接頭一天內(nèi)不宜頻繁動(dòng)作可以額外加一條動(dòng)作次數(shù)約束Σ|z_i,t - z_i,t-1| ≤ N_max這屬于“運(yùn)行可行性約束”本質(zhì)上是從實(shí)際運(yùn)維規(guī)則中提煉出來的不加的話結(jié)果容易在相鄰時(shí)段之間抖動(dòng)。5. 算例驗(yàn)證與結(jié)果分析看到數(shù)據(jù)背后的調(diào)度含義模型建好、代碼跑通這只是第一步。下面我用一個(gè)中等規(guī)模算例走一遍從數(shù)據(jù)設(shè)置到結(jié)果分析的全部過程。5.1 算例設(shè)置改進(jìn)的電氣互聯(lián)測試系統(tǒng)我采用的測試系統(tǒng)是IEEE 39節(jié)點(diǎn)電力系統(tǒng)與比利時(shí)20節(jié)點(diǎn)天然氣系統(tǒng)的耦合版本這是綜合能源系統(tǒng)研究中比較經(jīng)典的組合。關(guān)鍵耦合元件設(shè)置為4臺(tái)燃?xì)鈾C(jī)組接入氣網(wǎng)節(jié)點(diǎn)承擔(dān)系統(tǒng)約30%的負(fù)荷。2套P2G設(shè)備分別安裝在風(fēng)電富集節(jié)點(diǎn)額定功率分別為20MW和30MW轉(zhuǎn)換效率取60%。碳價(jià)設(shè)為100元/tCO2。風(fēng)電滲透率設(shè)為30%以典型日出力曲線輸入。無功補(bǔ)償裝置配置在負(fù)荷密集區(qū)節(jié)點(diǎn)單組容量5Mvar共12組。5.2 模型收斂性與求解性能分析在Matlab R2022b環(huán)境下用YALMIPGurobi 10.0求解模型規(guī)模大致是連續(xù)變量280多個(gè)0-1變量36個(gè)約束方程620多個(gè)。實(shí)測求解時(shí)間在45秒左右。這個(gè)速度在科研和離線規(guī)劃場景下完全夠用。對(duì)比不做SOCP松弛、直接上NLP內(nèi)點(diǎn)法的方案那個(gè)方案在部分初始點(diǎn)下會(huì)卡在局部最優(yōu)求解時(shí)間反而更長且結(jié)果不穩(wěn)定。SOCP方案的優(yōu)勢非常明顯。需要說明的是SOCP松弛在某些極端算例下可能產(chǎn)生非物理的松弛解。我在實(shí)際種遇到過一種情況氣網(wǎng)某些管道氣流方向不確定時(shí)Weymouth方程的絕對(duì)值展開會(huì)讓松弛不精確。解決方案是對(duì)每根管道引入方向變量d_i當(dāng)d_i1時(shí)氣流為正d_i0時(shí)氣流為負(fù)這種處理讓結(jié)果更接近真實(shí)物理過程。5.3 關(guān)鍵結(jié)果解讀從數(shù)據(jù)到運(yùn)行策略首先看優(yōu)化出的機(jī)組出力方案。在30%風(fēng)電滲透率下燃?xì)鈾C(jī)組的出力不是簡單地按經(jīng)濟(jì)性排序的因?yàn)樘技s束的存在部分高排放燃煤機(jī)組會(huì)被壓縮出力轉(zhuǎn)而由碳強(qiáng)度更低的燃?xì)鈾C(jī)組頂上。這就是碳價(jià)信號(hào)的傳導(dǎo)效果。再看無功優(yōu)化的效果。對(duì)比“有功-無功協(xié)同優(yōu)化”與“僅有功優(yōu)化最后按固定功率因數(shù)校核”兩種策略協(xié)同優(yōu)化下的系統(tǒng)網(wǎng)損降低了約8.3%最低節(jié)點(diǎn)電壓從0.93p.u.提升到0.97p.u.以上且所有節(jié)點(diǎn)電壓均在安全范圍內(nèi)。這個(gè)對(duì)比很好地說明了我前面強(qiáng)調(diào)的觀點(diǎn)——只在末端校核無功是無法主動(dòng)優(yōu)化電壓分布的。P2G設(shè)備的運(yùn)行策略也很有意思。風(fēng)電大發(fā)時(shí)段P2G滿發(fā)消耗富余電力相當(dāng)于給系統(tǒng)提供了一個(gè)可調(diào)節(jié)的“電負(fù)荷”緩解了棄風(fēng)問題同時(shí)產(chǎn)出的天然氣進(jìn)入管網(wǎng)提高了氣網(wǎng)供氣充裕度。而在負(fù)荷高峰時(shí)段P2G自動(dòng)降出力甚至停機(jī)把電力優(yōu)先讓給用戶。這種“削峰填谷”雙向調(diào)節(jié)作用正是電氣互聯(lián)系統(tǒng)相比單一電網(wǎng)的獨(dú)特優(yōu)勢。5.4 碳價(jià)靈敏度分析碳中和政策的量化評(píng)估我還做了碳價(jià)的靈敏度分析觀察碳價(jià)從0元/tCO2逐步提高到300元/tCO2的過程中系統(tǒng)運(yùn)行策略的變化規(guī)律碳價(jià)為0時(shí)系統(tǒng)優(yōu)先選擇煤電碳排放總量最高。碳價(jià)為50元/tCO2時(shí)燃?xì)鈾C(jī)組開始逐步替代煤電。碳價(jià)達(dá)到200元/tCO2時(shí)P2G設(shè)備利用率顯著上升系統(tǒng)碳排量在原有基礎(chǔ)上下降了約25%。這個(gè)分析對(duì)實(shí)際政策制定和投資規(guī)劃很有參考價(jià)值——它把碳價(jià)這一抽象的政策變量轉(zhuǎn)化成了看得見的機(jī)組啟停次序和碳排總量的變化曲線。作為工程師能給出這樣的量化分析跟只說“提高碳價(jià)有利于減排”是兩種完全不同的說服力。6. 工程化落地中的避坑清單與擴(kuò)展思路6.1 五個(gè)最容易翻車的細(xì)節(jié)全部來自實(shí)際調(diào)試第一個(gè)坑大氣網(wǎng)系統(tǒng)的Weymouth方程數(shù)值病態(tài)問題。氣壓值的量級(jí)通常在幾兆帕到十兆帕之間平方后數(shù)值差異巨大導(dǎo)致矩陣條件數(shù)極差求解器經(jīng)常報(bào)數(shù)值問題。我嘗試過把氣壓單位換成bar、把方程兩邊同時(shí)除以基準(zhǔn)值等方法最有效的還是把所有物理量都標(biāo)幺化處理讓變量落在0.1到2這個(gè)區(qū)間內(nèi)。第二個(gè)坑數(shù)據(jù)不一致性。電網(wǎng)和氣網(wǎng)的基準(zhǔn)功率、基準(zhǔn)電壓如果不統(tǒng)一耦合元件功率和氣網(wǎng)的對(duì)應(yīng)關(guān)系會(huì)直接算錯(cuò)。比如燃?xì)鈾C(jī)組效率用默認(rèn)高估時(shí)會(huì)出現(xiàn)氣網(wǎng)氣源供氣量不足而電網(wǎng)出力卻滿發(fā)的矛盾結(jié)果。第三個(gè)坑SOCP松弛方向不正確。之前說過必須是我還見過有人寫成的那是在人為縮小可行域有可能直接導(dǎo)致無解。第四個(gè)坑YALMIP和求解器的版本兼容性。不同版本對(duì)optimize返回狀態(tài)的定義有細(xì)微差別。我在調(diào)試時(shí)曾遇到過的困惑包括明明有可行解但Gurobi報(bào)infeasible是求解器設(shè)置問題改成雙精度求解精度后收斂性明顯改善。第五個(gè)坑約束冗余導(dǎo)致求解變慢。初始建模階段為了安全起見添加了大量冗余約束結(jié)果導(dǎo)致LP松弛的求解效率降低了接近一倍。后來專門做了約束削減把可被其他約束自動(dòng)滿足的冗余項(xiàng)刪掉這個(gè)問題才緩解。6.2 模型向更高維度擴(kuò)展的三個(gè)方向目前的模型是日前穩(wěn)態(tài)優(yōu)化實(shí)際工程應(yīng)用有三個(gè)擴(kuò)展方向。引入時(shí)序耦合約束加入機(jī)組啟停成本、爬坡約束、儲(chǔ)能充放電約束把模型擴(kuò)展為多時(shí)段動(dòng)態(tài)優(yōu)化。這在Matlab中的實(shí)現(xiàn)并不復(fù)雜只需把變量從一維向量變成二維矩陣并把相鄰時(shí)段的耦合約束加進(jìn)去。進(jìn)一步細(xì)化天然氣系統(tǒng)的動(dòng)態(tài)特性比如加入管道存氣、暫態(tài)慣性等。日前調(diào)度用穩(wěn)態(tài)模型是合理的但日內(nèi)滾動(dòng)修正時(shí)動(dòng)態(tài)信息能顯著提高調(diào)度精度。引入不確定性魯棒優(yōu)化或隨機(jī)優(yōu)化思路風(fēng)光出力的預(yù)測誤差需要在模型里顯式建模用魯棒邊界優(yōu)化的方法把最壞情況下的安全約束納入考慮。另外如果系統(tǒng)的無功電壓問題比較突出建議進(jìn)一步加入變電站動(dòng)態(tài)無功補(bǔ)償設(shè)備的詳細(xì)模型比如SVC和STATCOM的運(yùn)行約束——它們的快速調(diào)節(jié)能力在應(yīng)對(duì)新能源功率波動(dòng)時(shí)效果非常明顯。6.3 關(guān)于代碼復(fù)用和工程部署的一個(gè)建議最后這點(diǎn)是給想把模型用于實(shí)際項(xiàng)目的讀者的。不建議把這里說的所謂“Matlab代碼實(shí)現(xiàn)”停留在一次性腳本的層面。把整個(gè)優(yōu)化模型封裝成一個(gè)可重復(fù)調(diào)用的黑盒對(duì)外只暴露系統(tǒng)數(shù)據(jù)和場景參數(shù)對(duì)內(nèi)統(tǒng)一調(diào)取YALMIP建模和求解器接口。這樣后續(xù)做多場景比對(duì)、方案篩選甚至對(duì)接上級(jí)調(diào)度系統(tǒng)都會(huì)輕松得多。我的慣用做法是封裝成一個(gè)function [P_g_opt, Q_g_opt, U_opt, F_gas_opt] run_IES_optimization(caseData, carbonPrice, windScenario, options)函數(shù)所有情景參數(shù)都是輸入?yún)?shù)內(nèi)部自動(dòng)完成建模、求解和結(jié)果整理。這樣跑碳價(jià)靈敏度分析、不同的風(fēng)電出力場景只需要在主腳本里寫循環(huán)調(diào)用不用反復(fù)改模型代碼。這套工作流在我后續(xù)的好幾個(gè)項(xiàng)目里都在不斷復(fù)用?；仡^再看電氣互聯(lián)系統(tǒng)的有功-無功協(xié)同優(yōu)化表面上是一個(gè)數(shù)學(xué)優(yōu)化問題實(shí)際上考驗(yàn)的是對(duì)兩條物理網(wǎng)絡(luò)的深刻理解和建模功力。把模型建立在清晰的物理認(rèn)知之上用合適的松弛手段讓求解器高效工作再通過嚴(yán)謹(jǐn)?shù)慕Y(jié)果分析反哺對(duì)系統(tǒng)運(yùn)行規(guī)律的認(rèn)識(shí)這套方法論才是最值得沉淀的東西。 電氣互聯(lián)系統(tǒng)有功-無功協(xié)同優(yōu)化：碳約束下的建模與求解
1. “碳中和”目標(biāo)下的電氣互聯(lián)系統(tǒng)為什么非做有功-無功協(xié)同不可先說一個(gè)我在做這類項(xiàng)目時(shí)最深的體會(huì)很多做電力系統(tǒng)優(yōu)化的同學(xué)把重點(diǎn)全放在有功調(diào)度上無功這塊要么忽略要么用固定的功率因數(shù)折算一下。但在“碳中和”目標(biāo)下電氣互聯(lián)系統(tǒng)的優(yōu)化邊界已經(jīng)被徹底拓寬了——你面對(duì)的不再是單一電網(wǎng)而是電力系統(tǒng)與天然氣系統(tǒng)深度耦合的網(wǎng)絡(luò)負(fù)荷側(cè)多了電轉(zhuǎn)氣、氫能等元素電源側(cè)多了大量分布式可再生能源。此時(shí)如果還把有功和無功拆開做結(jié)果不是次優(yōu)而是可能直接得出不可行的運(yùn)行方案。所謂電氣互聯(lián)系統(tǒng)通俗講就是電網(wǎng)和氣網(wǎng)通過兩類關(guān)鍵元件“綁”在了一起一類是燃?xì)鈾C(jī)組把天然氣轉(zhuǎn)化為電力另一類是電轉(zhuǎn)氣設(shè)備Power to GasP2G把富余電能轉(zhuǎn)化為天然氣。前者是天然氣系統(tǒng)向電力系統(tǒng)供電后者是電力系統(tǒng)反過來向天然氣系統(tǒng)供氣。這種雙向耦合的直接后果就是電網(wǎng)的有功潮流變化會(huì)改變氣網(wǎng)的供氣需求氣網(wǎng)的供氣壓降又反過來約束燃?xì)鈾C(jī)組的出力上限而機(jī)組無功出力調(diào)整又會(huì)改變線路無功潮流和電壓分布——環(huán)環(huán)相扣誰也沒法獨(dú)立優(yōu)化。為什么“碳中和”目標(biāo)讓這個(gè)問題更棘手因?yàn)楦弑壤稍偕茉唇尤牒笙到y(tǒng)慣量降低、電壓支撐能力變?nèi)鯚o功資源變得更加稀缺。光伏和風(fēng)電在滿發(fā)時(shí)可能擠占常規(guī)機(jī)組的出清空間但它們的無功調(diào)節(jié)能力有限甚至某些工況下還要吸收無功。另一個(gè)問題是碳排放約束把氣網(wǎng)側(cè)的成本結(jié)構(gòu)改變了——碳價(jià)、碳配額都進(jìn)入了目標(biāo)函數(shù)這就讓“電轉(zhuǎn)氣、氣轉(zhuǎn)電”的決策路徑變得更加復(fù)雜。我做這個(gè)項(xiàng)目的核心思路可以概括成一句話在滿足碳排放約束的前提下同時(shí)調(diào)節(jié)機(jī)組有功出力、無功出力、變壓器分接頭、無功補(bǔ)償裝置投切以及氣源供氣量讓整個(gè)電氣互聯(lián)系統(tǒng)的綜合運(yùn)行成本最低。這個(gè)思路拆解下來涉及三個(gè)核心問題如何建立電網(wǎng)和氣網(wǎng)統(tǒng)一的數(shù)學(xué)模型特別是描述耦合元件的變量關(guān)系。如何處理碳排放約束和目標(biāo)函數(shù)中碳成本的量化。如何在Matlab中高效求解這樣一個(gè)大規(guī)模非線性混合整數(shù)優(yōu)化問題。下面我會(huì)沿著這三條線把從建模到求解再到代碼實(shí)現(xiàn)的全過程完整分享出來。這篇文章適合正在做綜合能源系統(tǒng)優(yōu)化、無功優(yōu)化或者碳約束調(diào)度的研究生和工程師參考也適合想用Matlab搭建一體化優(yōu)化框架的讀者當(dāng)模板。2. 有功-無功協(xié)同的耦合機(jī)理數(shù)學(xué)上它們是怎么“糾纏”在一起的很多人覺得“有功-無功協(xié)同”就是把兩個(gè)目標(biāo)函數(shù)加到一起然后在約束里同時(shí)寫上機(jī)組有功上下限和無功上下限。這是典型的誤區(qū)。真正的耦合體現(xiàn)在等式約束和不等式約束交叉關(guān)聯(lián)上下面我把數(shù)學(xué)模型層面的耦合點(diǎn)逐一拆開。2.1 電網(wǎng)側(cè)數(shù)學(xué)模型潮流方程中的有功無功耦合電網(wǎng)側(cè)采用交流潮流模型時(shí)節(jié)點(diǎn)i的功率平衡方程是P_i U_i * Σ(U_j * (G_ij * cos(θ_ij) B_ij * sin(θ_ij)))Q_i U_i * Σ(U_j * (G_ij * sin(θ_ij) - B_ij * cos(θ_ij)))從這個(gè)式子就能看出有功P_i和無功Q_i共享同一組電壓幅值U和相角θ變量。改變?nèi)我庖粋€(gè)節(jié)點(diǎn)的無功全網(wǎng)的電壓分布都會(huì)變電壓變了相角也跟著變相角變了有功潮流就變了。這就是電網(wǎng)側(cè)最本質(zhì)的耦合。在工程簡化中有人用直流潮流模型只保留有功維度這在輸電網(wǎng)經(jīng)濟(jì)調(diào)度里勉強(qiáng)可用但一旦涉及電壓安全約束、無功補(bǔ)償優(yōu)化、變壓器分接頭調(diào)節(jié)直流潮流就完全失效了。本模型必須用交流潮流。對(duì)于輸電線路還有一個(gè)關(guān)鍵約束是線路潮流極限。實(shí)際傳輸容量受熱穩(wěn)定極限、電壓約束和穩(wěn)定極限三重限制用視在功率S_ij的平方表示就是S_ij^2 P_ij^2 Q_ij^2 ≤ (S_ij^max)^2這個(gè)約束本身就是有功與無功的非線性耦合約束。如果只優(yōu)化有功而忽略無功線路可能會(huì)在視在功率上越限而單純看有功卻一切正常。2.2 天然氣網(wǎng)絡(luò)模型氣體流動(dòng)方程與節(jié)點(diǎn)氣壓約束天然氣網(wǎng)絡(luò)的穩(wěn)態(tài)模型常用節(jié)點(diǎn)氣流平衡方程描述。對(duì)氣網(wǎng)節(jié)點(diǎn)n氣流平衡可以寫成Σ(FLOW_k) G_source,m - G_load,m - G_p2g,m 0其中FLOW_k表示管道k的氣流量正方向?yàn)榱魅牍?jié)點(diǎn)反方向?yàn)榱鞒鯣_source是氣源供氣量G_load是氣負(fù)荷包括燃?xì)鈾C(jī)組耗氣量和民用氣負(fù)荷G_p2g是P2G設(shè)備消耗的電量轉(zhuǎn)化成的天然氣注入量。管道氣流FLOW與兩端氣壓的關(guān)系通常用Weymouth穩(wěn)態(tài)方程FLOW_ij^2 K_ij^2 * (p_i^2 - p_j^2)這個(gè)式子是非線性的而且是本模型的第一個(gè)非線性難點(diǎn)來源。注意這里FLOW是有方向的實(shí)際工程中需要處理方向符號(hào)問題常見做法是引入雙向流標(biāo)志位或者直接假設(shè)氣流方向已知——在做日前優(yōu)化調(diào)度時(shí)氣流方向通常可以從潮流結(jié)果預(yù)判但嚴(yán)謹(jǐn)起見最好在約束中保留方向變量。氣壓還有上下限約束p_min ≤ p_i ≤ p_max氣源出力也有爬坡限制G_source,min ≤ G_source ≤ G_source,max這里有個(gè)值得注意的細(xì)節(jié)氣網(wǎng)的慣性比電網(wǎng)大得多但穩(wěn)態(tài)模型下管道存氣linepack效應(yīng)被忽略了。在日內(nèi)滾動(dòng)優(yōu)化中l(wèi)inepack往往能提供寶貴的短時(shí)緩沖能力但本文的日前模型先不考慮保持模型相對(duì)簡潔。2.3 耦合元件模型燃?xì)鈾C(jī)組與P2G的雙向耦合燃?xì)鈾C(jī)組是氣網(wǎng)到電網(wǎng)的耦合元件其模型包括兩部分注入電網(wǎng)的有功出力P_g和旋轉(zhuǎn)備用容量。從氣網(wǎng)消耗的天然氣量F_gas a * P_g^2 b * P_g c其中a、b、c是燃料消耗系數(shù)。這個(gè)二次函數(shù)也是典型的非線性環(huán)節(jié)。另一個(gè)方向——P2G設(shè)備的耦合模型為G_p2g η_p2g * P_p2g / GHV其中P_p2g是P2G設(shè)備消耗的有功功率η_p2g是轉(zhuǎn)換效率GHV是天然氣高位熱值。P2G消耗的功率上限同時(shí)受電網(wǎng)上網(wǎng)容量和氣網(wǎng)注入能力的雙重約束P_p2g,min ≤ P_p2g ≤ min(P_p2g,max, G_p2g,max * GHV / η_p2g)到這里有功-無功耦合的數(shù)學(xué)畫像就清晰了電網(wǎng)側(cè)P、Q通過電壓幅值和相角在潮流方程中耦合。氣網(wǎng)側(cè)氣流量通過壓力非線性方程耦合。系統(tǒng)級(jí)燃?xì)鈾C(jī)組和P2G在有功功率與氣流量之間建立雙向通道而這種耦合通過上述兩個(gè)非線性環(huán)節(jié)傳遞。抽象成一句話有功決策變了無功優(yōu)化空間也跟著變無功決策變了電壓分布變了P2G消耗的有功上限和燃?xì)鈾C(jī)組的可用容量也變了。所謂“協(xié)同”就是在數(shù)學(xué)上同時(shí)處理這兩組變量而不是把它們解耦成兩個(gè)子問題分別求解后再拼起來。3. 協(xié)同優(yōu)化模型搭建目標(biāo)函數(shù)與關(guān)鍵約束的設(shè)計(jì)思路我在這類項(xiàng)目上走過的彎路不少最值錢的一條經(jīng)驗(yàn)是目標(biāo)函數(shù)決定了優(yōu)化方向約束決定可行域邊界但真正決定模型能否被高效求解的是約束的“數(shù)學(xué)形態(tài)”是線性、凸還是非凸。3.1 目標(biāo)函數(shù)設(shè)計(jì)成本與碳排放的雙重權(quán)衡本模型的目標(biāo)函數(shù)我在實(shí)際項(xiàng)目中是這樣設(shè)計(jì)的綜合考慮了四種成本min F_total Σ(C_fuel) Σ(C_co2) Σ(C_p2g) Σ(C_grid)第一項(xiàng)是常規(guī)機(jī)組燃料成本是機(jī)組有功出力的二次函數(shù)C_fuel,i a_i * P_i^2 b_i * P_i c_i第二項(xiàng)是燃?xì)鈾C(jī)組的碳排放成本這里引入“碳中和”目標(biāo)的關(guān)鍵我在實(shí)際建模中用了兩種處理方式碳市場交易價(jià)格乘以排放量或者設(shè)定碳配額上限然后用階梯碳價(jià)。模型里我建議使用線性碳價(jià)對(duì)單位碳排放量乘以碳價(jià)系數(shù)C_co2,i ε_(tái)i * e_i * P_i其中ε_(tái)i是碳價(jià)元/tCO2e_i是單位電量的碳排放強(qiáng)度tCO2/MWh。第三項(xiàng)是P2G運(yùn)行成本主要包括設(shè)備運(yùn)維成本和購電成本分?jǐn)侰_p2g γ * P_p2g其中γ是P2G單位運(yùn)行成本系數(shù)。第四項(xiàng)是從上級(jí)電網(wǎng)購電的成本只有當(dāng)系統(tǒng)接入主網(wǎng)時(shí)才需要。這里我想特別說明為什么要這樣設(shè)目標(biāo)函數(shù)而不是簡單加一個(gè)“碳排放最小”目標(biāo)。因?yàn)榧儨p排目標(biāo)在經(jīng)濟(jì)上不可行實(shí)際運(yùn)行中調(diào)度員關(guān)心的是“在碳約束下的最經(jīng)濟(jì)運(yùn)行”所以碳成本進(jìn)入目標(biāo)函數(shù)、碳約束進(jìn)入不等式約束這兩個(gè)維度分開處理才能保證結(jié)果既環(huán)保又可執(zhí)行。3.2 等式約束詳解電網(wǎng)潮流與氣網(wǎng)平衡電場部分需要同時(shí)滿足有功和無功平衡P_G,i - P_L,i - U_i * Σ(U_j * (G_ij * cos θ_ij B_ij * sin θ_ij)) - P_P2G,i 0Q_G,i - Q_L,i - U_i * Σ(U_j * (G_ij * sin θ_ij - B_ij * cos θ_ij)) - Q_P2G,i 0注意我在P2G節(jié)點(diǎn)的有功平衡里專門扣掉了P_P2G,i這是典型的“負(fù)荷型元件”處理方式P2G消耗有功同時(shí)如果是電解水制氫再合成甲烷的路徑還需要向電網(wǎng)提供少量的無功給整流器等設(shè)備。氣網(wǎng)側(cè)的氣流平衡在2.2節(jié)已經(jīng)給出。這里補(bǔ)充說明一個(gè)細(xì)節(jié)對(duì)于含壓縮機(jī)站的天然氣網(wǎng)絡(luò)壓縮機(jī)的耗氣量也是變量但在基礎(chǔ)模型中往往先用固定比例折算C_comp β * FLOW_comp這種簡化在工程上是可接受的因?yàn)閴嚎s機(jī)耗氣量占總輸氣量的比例通常只有1%到3%不會(huì)對(duì)優(yōu)化結(jié)果產(chǎn)生方向性影響。3.3 不等式約束安全與設(shè)備極限的建模技巧不等式約束是我在實(shí)際代碼中掉過坑最多的地方。常見的約束包括發(fā)電機(jī)組有功出力上下限P_min ≤ P ≤ P_max。無功出力上下限Q_min ≤ Q ≤ Q_max這里要注意無功上限與有功出力有關(guān)嚴(yán)格來說應(yīng)該是Q_max(P)但很多模型直接給定常數(shù)上限會(huì)帶來結(jié)果偏樂觀的問題。節(jié)點(diǎn)電壓幅值約束U_min ≤ U ≤ U_max在碳中和場景下電壓約束的重要性會(huì)被放大新能源出力的隨機(jī)性讓電壓越限風(fēng)險(xiǎn)顯著增加。變壓器分接頭調(diào)節(jié)范圍約束。氣網(wǎng)節(jié)點(diǎn)氣壓約束和管道輸氣能力約束。線路傳輸容量約束。在眾多不等式約束里最難處理的是線路容量約束。這個(gè)約束用線性化近似會(huì)丟失無功信息用完整交流潮流又帶來大量非線性。我的做法是保留交流潮流形式的視在功率約束配合松弛技術(shù)確?？山庑?。還有個(gè)容易被忽略的約束是P2G設(shè)備的運(yùn)行約束。P2G啟動(dòng)需要時(shí)間電解槽的爬坡速率有限-P_p2g,ramp ≤ P_p2g,t - P_p2g,t-1 ≤ P_p2g,ramp如果忽略這個(gè)約束優(yōu)化器會(huì)讓P2G在相鄰時(shí)段隨意跳變結(jié)果在實(shí)際中根本執(zhí)行不了——這在電網(wǎng)側(cè)調(diào)度中也是同樣的道理。3.4 碳排放約束的建模與處理碳約束我傾向于直接寫成總量約束的形式Σ(E_G,i) ≤ E_total,max其中E_G,i是第i臺(tái)機(jī)組的碳排放量E_total,max是系統(tǒng)設(shè)定的碳排放上限。在“碳中和”目標(biāo)下這個(gè)上限應(yīng)該逐年收緊。如果進(jìn)一步要求模型能反映不同碳價(jià)情景下的運(yùn)行策略差異還可以把碳約束轉(zhuǎn)化為拉格朗日懲罰項(xiàng)或者KKT乘子。但對(duì)優(yōu)化模型來說顯式碳約束更直觀且方便做邊際排放成本的影子價(jià)格分析——影子價(jià)格本身就是很有用的調(diào)度決策輔助信息。4. MatLab代碼實(shí)現(xiàn)路線從數(shù)學(xué)模型到可運(yùn)行程序模型搭完之后怎么把它在Matlab里落地是另一個(gè)大問題。我分享一下自己的完整技術(shù)路線包括求解工具選型、代碼架構(gòu)和核心函數(shù)的實(shí)現(xiàn)思路。4.1 求解工具與算法選型這個(gè)問題先說結(jié)論我選的是YALMIP建模工具箱 二階錐松弛SOCP 商業(yè)求解器CPLEX或Gurobi的組合方案。為什么這么選電氣互聯(lián)系統(tǒng)優(yōu)化問題本質(zhì)上是MINLP混合整數(shù)非線性規(guī)劃直接求解全局最優(yōu)幾乎不可能。工程上主要有四條路線內(nèi)點(diǎn)法直接求解連續(xù)松弛的NLP速度快但容易陷入局部最優(yōu)且很難處理整數(shù)變量。智能算法如粒子群和遺傳算法不依賴梯度能處理非線性但每次求解的穩(wěn)定性差而且沒有最優(yōu)性保證。Benders分解把原問題分解成主問題和子問題適合大規(guī)模系統(tǒng)但模型推導(dǎo)復(fù)雜代碼實(shí)現(xiàn)成本高。二階錐松弛把非凸的Weymouth方程和交流潮流約束松弛為凸約束配合Big-M法處理離散變量形成MISOCP問題商用求解器對(duì)這類問題的求解能力已經(jīng)非常成熟。我實(shí)測下來二階錐松弛商業(yè)求解器在幾十個(gè)節(jié)點(diǎn)的系統(tǒng)上求解速度在幾十秒到幾分鐘之間兼顧了精度和速度是最適合Matlab實(shí)現(xiàn)的技術(shù)路線。具體做法是把Weymouth方程做二階錐轉(zhuǎn)換定義變量α_i p_i^2任意變換后約束變成下式FLOW_ij^2 ≤ K_ij^2 * (α_i - α_j)同時(shí)把潮流方程中的二次項(xiàng)通過電壓幅值平方替換處理。注意這里的松弛方向是有講究的——必須松弛成小于等于的形式才能保持凸性。雖然松弛擴(kuò)大了可行域但實(shí)際運(yùn)行中優(yōu)化結(jié)果一般會(huì)自動(dòng)逼近等式邊界。4.2 代碼架構(gòu)分模塊設(shè)計(jì)的思路我的Matlab代碼按模塊組織便于調(diào)試和復(fù)用主文件main.m定義系統(tǒng)參數(shù)加載數(shù)據(jù)調(diào)用建模函數(shù)求解輸出結(jié)果。數(shù)據(jù)定義文件case_data.m電力系統(tǒng)節(jié)點(diǎn)參數(shù)、氣網(wǎng)管道參數(shù)、耦合元件參數(shù)等基礎(chǔ)數(shù)據(jù)。建模函數(shù)build_opf_model.m用YALMIP定義優(yōu)化變量、目標(biāo)函數(shù)和約束。結(jié)果分析函數(shù)plot_results.m繪制電壓分布、機(jī)組出力、氣網(wǎng)壓力、碳排量變化等曲線。核心的YALMIP建模框架大概是下面的思路以目標(biāo)函數(shù)和關(guān)鍵約束為例偽代碼% YALMIP變量定義 P_g sdpvar(n_gen, 1); % 發(fā)電機(jī)有功 Q_g sdpvar(n_gen, 1); % 發(fā)電機(jī)無功 U sdpvar(n_bus, 1); % 節(jié)點(diǎn)電壓幅值 theta sdpvar(n_bus, 1); % 節(jié)點(diǎn)相角 P_p2g sdpvar(n_p2g, 1); % P2G有功消耗 G_source sdpvar(n_source, 1); % 氣源供氣量 FLOW sdpvar(n_pipe, 1); % 管道流量 % 目標(biāo)函數(shù) Objective sum(a .* P_g.^2 b .* P_g c) ... sum(carbon_price .* emission_coeff .* P_g) ... sum(gamma .* P_p2g); % 約束定義 Constraints []; % 節(jié)點(diǎn)有功平衡 Constraints [Constraints, A_incidence * P_g P_wind - P_load - P_p2g_map ... U .* (G_ij * cos(theta_ij) B_ij * sin(theta_ij))]; % 節(jié)點(diǎn)無功平衡 Constraints [Constraints, Q_g Q_comp - Q_load ... U .* (G_ij * sin(theta_ij) - B_ij * cos(theta_ij))]; % 燃?xì)鈾C(jī)組氣負(fù)荷 Constraints [Constraints, F_gas a_fuel * P_g.^2 b_fuel .* P_g c_fuel]; % 氣網(wǎng)節(jié)點(diǎn)平衡 Constraints [Constraints, A_gas * FLOW G_source - G_load - F_gas 0]; % 二階梯松弛 for i 1:n_pipe Constraints [Constraints, FLOW(i)^2 K(i)^2 * (alpha_from(i) - alpha_to(i))]; end % 求解 ops sdpsettings(solver,gurobi,verbose,1); optimize(Constraints, Objective, ops);這段代碼里最關(guān)鍵的幾個(gè)細(xì)節(jié)值得說明第一潮流和網(wǎng)絡(luò)拓?fù)湎嚓P(guān)的矩陣G_ij、B_ij不是手工輸入的而是用Matpower的makeYbus函數(shù)自動(dòng)生成的。這樣即使換算例代碼不用改結(jié)構(gòu)。第二P2G映射關(guān)系P_p2g_map的處理是我調(diào)試時(shí)踩過的一個(gè)坑P2G設(shè)備接在哪個(gè)節(jié)點(diǎn)上必須在數(shù)據(jù)文件里明確指定并且要形成從變量索引到網(wǎng)絡(luò)節(jié)點(diǎn)的映射矩陣否則氣網(wǎng)和電網(wǎng)的平衡方程會(huì)錯(cuò)位。第三二階梯約束的寫法里我用的是而不是這關(guān)系到問題是否凸。如果用等式約束YALMIP會(huì)把它當(dāng)作非線性等式來解速度慢幾個(gè)數(shù)量級(jí)還可能因?yàn)槌跏键c(diǎn)差解不出來。4.3 針對(duì)整數(shù)變量的處理變壓器分接頭與無功補(bǔ)償投切變壓器分接頭和無功補(bǔ)償電容/電抗器的投切是離散變量我一開始直接把分接頭位置定義為整數(shù)變量binvar或intvar導(dǎo)致求解時(shí)間暴增。后來改成Big-M法處理先把分接頭比轉(zhuǎn)換為連續(xù)變量t并添加線性約束t_min Σ(step_i * z_i) ≤ t ≤ t_min Σ(step_i * z_i) (1 - z_i) * M其中z_i是0-1變量M是足夠大的數(shù)。這種做法把離散決策從“選位置”轉(zhuǎn)換成“選步進(jìn)”在Gurobi里求解效率高多了。還有一個(gè)實(shí)用技巧如果只是做日前調(diào)度變壓器分接頭一天內(nèi)不宜頻繁動(dòng)作可以額外加一條動(dòng)作次數(shù)約束Σ|z_i,t - z_i,t-1| ≤ N_max這屬于“運(yùn)行可行性約束”本質(zhì)上是從實(shí)際運(yùn)維規(guī)則中提煉出來的不加的話結(jié)果容易在相鄰時(shí)段之間抖動(dòng)。5. 算例驗(yàn)證與結(jié)果分析看到數(shù)據(jù)背后的調(diào)度含義模型建好、代碼跑通這只是第一步。下面我用一個(gè)中等規(guī)模算例走一遍從數(shù)據(jù)設(shè)置到結(jié)果分析的全部過程。5.1 算例設(shè)置改進(jìn)的電氣互聯(lián)測試系統(tǒng)我采用的測試系統(tǒng)是IEEE 39節(jié)點(diǎn)電力系統(tǒng)與比利時(shí)20節(jié)點(diǎn)天然氣系統(tǒng)的耦合版本這是綜合能源系統(tǒng)研究中比較經(jīng)典的組合。關(guān)鍵耦合元件設(shè)置為4臺(tái)燃?xì)鈾C(jī)組接入氣網(wǎng)節(jié)點(diǎn)承擔(dān)系統(tǒng)約30%的負(fù)荷。2套P2G設(shè)備分別安裝在風(fēng)電富集節(jié)點(diǎn)額定功率分別為20MW和30MW轉(zhuǎn)換效率取60%。碳價(jià)設(shè)為100元/tCO2。風(fēng)電滲透率設(shè)為30%以典型日出力曲線輸入。無功補(bǔ)償裝置配置在負(fù)荷密集區(qū)節(jié)點(diǎn)單組容量5Mvar共12組。5.2 模型收斂性與求解性能分析在Matlab R2022b環(huán)境下用YALMIPGurobi 10.0求解模型規(guī)模大致是連續(xù)變量280多個(gè)0-1變量36個(gè)約束方程620多個(gè)。實(shí)測求解時(shí)間在45秒左右。這個(gè)速度在科研和離線規(guī)劃場景下完全夠用。對(duì)比不做SOCP松弛、直接上NLP內(nèi)點(diǎn)法的方案那個(gè)方案在部分初始點(diǎn)下會(huì)卡在局部最優(yōu)求解時(shí)間反而更長且結(jié)果不穩(wěn)定。SOCP方案的優(yōu)勢非常明顯。需要說明的是SOCP松弛在某些極端算例下可能產(chǎn)生非物理的松弛解。我在實(shí)際種遇到過一種情況氣網(wǎng)某些管道氣流方向不確定時(shí)Weymouth方程的絕對(duì)值展開會(huì)讓松弛不精確。解決方案是對(duì)每根管道引入方向變量d_i當(dāng)d_i1時(shí)氣流為正d_i0時(shí)氣流為負(fù)這種處理讓結(jié)果更接近真實(shí)物理過程。5.3 關(guān)鍵結(jié)果解讀從數(shù)據(jù)到運(yùn)行策略首先看優(yōu)化出的機(jī)組出力方案。在30%風(fēng)電滲透率下燃?xì)鈾C(jī)組的出力不是簡單地按經(jīng)濟(jì)性排序的因?yàn)樘技s束的存在部分高排放燃煤機(jī)組會(huì)被壓縮出力轉(zhuǎn)而由碳強(qiáng)度更低的燃?xì)鈾C(jī)組頂上。這就是碳價(jià)信號(hào)的傳導(dǎo)效果。再看無功優(yōu)化的效果。對(duì)比“有功-無功協(xié)同優(yōu)化”與“僅有功優(yōu)化最后按固定功率因數(shù)校核”兩種策略協(xié)同優(yōu)化下的系統(tǒng)網(wǎng)損降低了約8.3%最低節(jié)點(diǎn)電壓從0.93p.u.提升到0.97p.u.以上且所有節(jié)點(diǎn)電壓均在安全范圍內(nèi)。這個(gè)對(duì)比很好地說明了我前面強(qiáng)調(diào)的觀點(diǎn)——只在末端校核無功是無法主動(dòng)優(yōu)化電壓分布的。P2G設(shè)備的運(yùn)行策略也很有意思。風(fēng)電大發(fā)時(shí)段P2G滿發(fā)消耗富余電力相當(dāng)于給系統(tǒng)提供了一個(gè)可調(diào)節(jié)的“電負(fù)荷”緩解了棄風(fēng)問題同時(shí)產(chǎn)出的天然氣進(jìn)入管網(wǎng)提高了氣網(wǎng)供氣充裕度。而在負(fù)荷高峰時(shí)段P2G自動(dòng)降出力甚至停機(jī)把電力優(yōu)先讓給用戶。這種“削峰填谷”雙向調(diào)節(jié)作用正是電氣互聯(lián)系統(tǒng)相比單一電網(wǎng)的獨(dú)特優(yōu)勢。5.4 碳價(jià)靈敏度分析碳中和政策的量化評(píng)估我還做了碳價(jià)的靈敏度分析觀察碳價(jià)從0元/tCO2逐步提高到300元/tCO2的過程中系統(tǒng)運(yùn)行策略的變化規(guī)律碳價(jià)為0時(shí)系統(tǒng)優(yōu)先選擇煤電碳排放總量最高。碳價(jià)為50元/tCO2時(shí)燃?xì)鈾C(jī)組開始逐步替代煤電。碳價(jià)達(dá)到200元/tCO2時(shí)P2G設(shè)備利用率顯著上升系統(tǒng)碳排量在原有基礎(chǔ)上下降了約25%。這個(gè)分析對(duì)實(shí)際政策制定和投資規(guī)劃很有參考價(jià)值——它把碳價(jià)這一抽象的政策變量轉(zhuǎn)化成了看得見的機(jī)組啟停次序和碳排總量的變化曲線。作為工程師能給出這樣的量化分析跟只說“提高碳價(jià)有利于減排”是兩種完全不同的說服力。6. 工程化落地中的避坑清單與擴(kuò)展思路6.1 五個(gè)最容易翻車的細(xì)節(jié)全部來自實(shí)際調(diào)試第一個(gè)坑大氣網(wǎng)系統(tǒng)的Weymouth方程數(shù)值病態(tài)問題。氣壓值的量級(jí)通常在幾兆帕到十兆帕之間平方后數(shù)值差異巨大導(dǎo)致矩陣條件數(shù)極差求解器經(jīng)常報(bào)數(shù)值問題。我嘗試過把氣壓單位換成bar、把方程兩邊同時(shí)除以基準(zhǔn)值等方法最有效的還是把所有物理量都標(biāo)幺化處理讓變量落在0.1到2這個(gè)區(qū)間內(nèi)。第二個(gè)坑數(shù)據(jù)不一致性。電網(wǎng)和氣網(wǎng)的基準(zhǔn)功率、基準(zhǔn)電壓如果不統(tǒng)一耦合元件功率和氣網(wǎng)的對(duì)應(yīng)關(guān)系會(huì)直接算錯(cuò)。比如燃?xì)鈾C(jī)組效率用默認(rèn)高估時(shí)會(huì)出現(xiàn)氣網(wǎng)氣源供氣量不足而電網(wǎng)出力卻滿發(fā)的矛盾結(jié)果。第三個(gè)坑SOCP松弛方向不正確。之前說過必須是我還見過有人寫成的那是在人為縮小可行域有可能直接導(dǎo)致無解。第四個(gè)坑YALMIP和求解器的版本兼容性。不同版本對(duì)optimize返回狀態(tài)的定義有細(xì)微差別。我在調(diào)試時(shí)曾遇到過的困惑包括明明有可行解但Gurobi報(bào)infeasible是求解器設(shè)置問題改成雙精度求解精度后收斂性明顯改善。第五個(gè)坑約束冗余導(dǎo)致求解變慢。初始建模階段為了安全起見添加了大量冗余約束結(jié)果導(dǎo)致LP松弛的求解效率降低了接近一倍。后來專門做了約束削減把可被其他約束自動(dòng)滿足的冗余項(xiàng)刪掉這個(gè)問題才緩解。6.2 模型向更高維度擴(kuò)展的三個(gè)方向目前的模型是日前穩(wěn)態(tài)優(yōu)化實(shí)際工程應(yīng)用有三個(gè)擴(kuò)展方向。引入時(shí)序耦合約束加入機(jī)組啟停成本、爬坡約束、儲(chǔ)能充放電約束把模型擴(kuò)展為多時(shí)段動(dòng)態(tài)優(yōu)化。這在Matlab中的實(shí)現(xiàn)并不復(fù)雜只需把變量從一維向量變成二維矩陣并把相鄰時(shí)段的耦合約束加進(jìn)去。進(jìn)一步細(xì)化天然氣系統(tǒng)的動(dòng)態(tài)特性比如加入管道存氣、暫態(tài)慣性等。日前調(diào)度用穩(wěn)態(tài)模型是合理的但日內(nèi)滾動(dòng)修正時(shí)動(dòng)態(tài)信息能顯著提高調(diào)度精度。引入不確定性魯棒優(yōu)化或隨機(jī)優(yōu)化思路風(fēng)光出力的預(yù)測誤差需要在模型里顯式建模用魯棒邊界優(yōu)化的方法把最壞情況下的安全約束納入考慮。另外如果系統(tǒng)的無功電壓問題比較突出建議進(jìn)一步加入變電站動(dòng)態(tài)無功補(bǔ)償設(shè)備的詳細(xì)模型比如SVC和STATCOM的運(yùn)行約束——它們的快速調(diào)節(jié)能力在應(yīng)對(duì)新能源功率波動(dòng)時(shí)效果非常明顯。6.3 關(guān)于代碼復(fù)用和工程部署的一個(gè)建議最后這點(diǎn)是給想把模型用于實(shí)際項(xiàng)目的讀者的。不建議把這里說的所謂“Matlab代碼實(shí)現(xiàn)”停留在一次性腳本的層面。把整個(gè)優(yōu)化模型封裝成一個(gè)可重復(fù)調(diào)用的黑盒對(duì)外只暴露系統(tǒng)數(shù)據(jù)和場景參數(shù)對(duì)內(nèi)統(tǒng)一調(diào)取YALMIP建模和求解器接口。這樣后續(xù)做多場景比對(duì)、方案篩選甚至對(duì)接上級(jí)調(diào)度系統(tǒng)都會(huì)輕松得多。我的慣用做法是封裝成一個(gè)function [P_g_opt, Q_g_opt, U_opt, F_gas_opt] run_IES_optimization(caseData, carbonPrice, windScenario, options)函數(shù)所有情景參數(shù)都是輸入?yún)?shù)內(nèi)部自動(dòng)完成建模、求解和結(jié)果整理。這樣跑碳價(jià)靈敏度分析、不同的風(fēng)電出力場景只需要在主腳本里寫循環(huán)調(diào)用不用反復(fù)改模型代碼。這套工作流在我后續(xù)的好幾個(gè)項(xiàng)目里都在不斷復(fù)用?；仡^再看電氣互聯(lián)系統(tǒng)的有功-無功協(xié)同優(yōu)化表面上是一個(gè)數(shù)學(xué)優(yōu)化問題實(shí)際上考驗(yàn)的是對(duì)兩條物理網(wǎng)絡(luò)的深刻理解和建模功力。把模型建立在清晰的物理認(rèn)知之上用合適的松弛手段讓求解器高效工作再通過嚴(yán)謹(jǐn)?shù)慕Y(jié)果分析反哺對(duì)系統(tǒng)運(yùn)行規(guī)律的認(rèn)識(shí)這套方法論才是最值得沉淀的東西。