李卓然，南君培，王 超，孫 可，周 浩

（1.浙江大學(xué) 電氣工程學(xué)院，浙江 杭州310027；2.國網(wǎng)浙江省電力有限公司，浙江 杭州310007）

0 引 言

溫室氣體[1]排放的最大貢獻(xiàn)來自能源供應(yīng)中的煤炭、天然氣和石油的燃燒發(fā)電和供熱。為了實(shí)現(xiàn)減少溫室氣體排放量，必須大幅增加可再生能源的利用效率。隨著可再生能源[2-3]的區(qū)域電熱綜合能源系統(tǒng)中以燃?xì)廨啓C(jī)、分布式能源站為代表的城市天然氣應(yīng)用工程加速天然氣網(wǎng)絡(luò)的擴(kuò)展，燃?xì)廨啓C(jī)發(fā)電站點(diǎn)的建設(shè)以及天然氣網(wǎng)絡(luò)的擴(kuò)張，電-氣互聯(lián)綜合能源系統(tǒng)將成為未來重要的供能系統(tǒng)。同時(shí)，電轉(zhuǎn)氣裝置（power to gas，P2G）的出現(xiàn)進(jìn)一步促進(jìn)了電網(wǎng)和天然氣網(wǎng)絡(luò)的融合與耦合。但由于當(dāng)前電力網(wǎng)絡(luò)及天然氣網(wǎng)絡(luò)的能流管理及調(diào)度分屬不同主體，其聯(lián)合調(diào)度受到了物理因素的制約。

對區(qū)域電熱互聯(lián)的綜合能源系統(tǒng)合理建模與求解，從而獲取電熱互聯(lián)的區(qū)域綜合能源系統(tǒng)的穩(wěn)態(tài)運(yùn)行方式，對于電熱系統(tǒng)規(guī)劃、優(yōu)化運(yùn)行、事故檢修以及靜態(tài)安全分析[4]都有著重要指導(dǎo)作用。

本文提出了一種電-熱互聯(lián)綜合能源系統(tǒng)穩(wěn)態(tài)模型，能夠較為準(zhǔn)確地獲取考慮電熱耦合的綜合能源系統(tǒng)的電熱潮流[5]數(shù)據(jù)。以巴厘島電熱綜合能源系統(tǒng)算例驗(yàn)證了模型的準(zhǔn)確性，并在求解過程中比較了P-Q分解法和Fsolve兩種求解算法的優(yōu)劣。

1 區(qū)域電熱綜合能源系統(tǒng)穩(wěn)態(tài)潮流模型

1.1 熱力系統(tǒng)水力模型

水力模型約束方程描述為：

式中：A為網(wǎng)絡(luò)矩陣；˙m為管道流量；˙mq為節(jié)點(diǎn)負(fù)荷流量；B為回路矩陣；hf為管道壓降；K為管道阻力系數(shù)。

1.2 熱力系統(tǒng)熱力模型

熱力模型約束方程描述為：

式中：Φ為熱負(fù)荷；Ts為供水溫度；To為回水溫度；Tstart為管道首端溫度；Tend為管道末端溫度；λ為傳熱系數(shù)；min為流入節(jié)點(diǎn)管道流量；mout為流出節(jié)點(diǎn)管道流量；Tin為管道末端的溫度；Tout為混合溫度。

1.3 電熱耦合元件模型

電熱耦合元件的代表是熱電聯(lián)產(chǎn)機(jī)組[6]，由燃?xì)廨啓C(jī)及發(fā)電機(jī)與余熱鍋爐、蒸汽輪機(jī)或供熱式蒸汽輪機(jī)（抽汽式或背壓式）共同組成的循環(huán)系統(tǒng)，它將燃?xì)廨啓C(jī)做功后的高溫乏煙氣熱量通過余熱鍋爐回收轉(zhuǎn)換為蒸汽熱量，蒸汽注入蒸汽輪機(jī)發(fā)電，或?qū)⒉糠职l(fā)電作功后的乏汽用于供熱。

根據(jù)熱電聯(lián)產(chǎn)機(jī)組的熱電比是否隨著出力改變，可以分為定熱電比和變熱電比機(jī)組兩類。

式中：φCHP是CHP機(jī)組的熱出力；PCHP是CHP機(jī)組的電出力；η為CHP機(jī)組冷凝效率；Fin為燃料輸入速率；Cm為一恒定量；CZ是隨著燃料輸入速率變化的量。

2 算例分析

本文以巴厘島電熱模型[7]為例。如圖1所示，巴厘島電熱模型為32節(jié)點(diǎn)熱網(wǎng)，9節(jié)點(diǎn)電網(wǎng)通過3臺熱電聯(lián)產(chǎn)機(jī)組耦合組成。

圖1 巴厘島電熱模型系統(tǒng)結(jié)構(gòu)圖

其中，CHP1和CHP3為定電熱比熱電機(jī)組，熱電比參數(shù)如表1所示。

表1 機(jī)組熱電比參數(shù)

電網(wǎng)基準(zhǔn)功率為1 MW，基準(zhǔn)電壓為11 kV；導(dǎo)線每千米阻抗為z=（0.164+j0.08）Ω；主變壓器采用33/11.5 kV，額定容量15 MVA，額定電壓百分?jǐn)?shù)為18%；母線I、III、IV、V和VI的有功負(fù)荷分別為0.2 MW、0.5 MW、0.5 MW、0.2 MW、0.2 MW；各用電負(fù)荷功率因數(shù)為1；三個(gè)電源電壓幅值為分別為1.02pu，1.05pu，1.05pu；電源1的電壓相角為0°；電力管道長度如圖1所示；電力節(jié)點(diǎn)IX為平衡節(jié)點(diǎn)。

熱網(wǎng)總負(fù)荷為2.164 MW，節(jié)點(diǎn)1、31和節(jié)點(diǎn)32分別對應(yīng)三個(gè)電源CHP2、CHP1和CHP3；CHP3的電出力P3=0.3 MW，根據(jù)機(jī)組3的熱電比可以求出機(jī)組3的熱出力φCHP3=0.3797 MW；熱力管網(wǎng)供熱溫度為70℃，熱力管網(wǎng)回水溫度為30℃；管道所在環(huán)境溫度為10℃；水密度ρ=958.4 kg/m3；水的運(yùn)動粘度μ=0.294×10-6m2/s；熱力管網(wǎng)參數(shù)如圖1所示；熱網(wǎng)節(jié)點(diǎn)1為平衡節(jié)點(diǎn)。

3 區(qū)域電熱綜合能源系統(tǒng)穩(wěn)態(tài)潮流模型求解方法比較

3.1 運(yùn)用牛頓拉斐遜法對電熱約束方程來回迭代求解的算法

電熱互聯(lián)的牛頓拉斐遜迭代法[8]的求解過程是先根據(jù)熱網(wǎng)約束中的水力方程計(jì)算出一組管道流量m0作為迭代初值；然后之后的循環(huán)過程是先由熱網(wǎng)約束方程計(jì)算出對應(yīng)的雅各比矩陣J1并帶入初值m0得到m1、T、TM；根據(jù)管道1的質(zhì)量流量計(jì)算出機(jī)組CHP2的熱出力，并由機(jī)組CHP2的熱電比參數(shù)得到相應(yīng)機(jī)組CHP2的電出力；然后將機(jī)組CHP2的電出力作為已知量，計(jì)算電網(wǎng)約束方程的雅各比矩陣并求出電力潮流參數(shù)；然后根據(jù)求解出的母線2的電出力算出機(jī)組CHP1的電出力，并由CHP1的熱電比參數(shù)機(jī)組CHP1的熱出力；接著判斷所求管道流量等參數(shù)是否收斂并滿足求解精度；然后將機(jī)組CHP1的熱出力作為已知量帶入水力熱力計(jì)算流程開始下一次的迭代過程。牛頓拉斐遜法迭代法的計(jì)算流程如圖2所示。

圖2 牛頓拉斐遜法迭代法求解穩(wěn)態(tài)潮流的流程示意圖

3.2 運(yùn)用Fsolve函數(shù)統(tǒng)一處理電熱非線性方程組的算法

Fsolve函數(shù)的電熱互聯(lián)系統(tǒng)的穩(wěn)態(tài)潮流求解算法流程上與牛頓拉斐遜法迭代法的求解過程類似。不同的是：（1）迭代初值需要自己設(shè)定并且初值的選取可能直接影響到結(jié)果是否收斂。選取不同的初值對求解的結(jié)果影響較大，初值選取的不合適可能會導(dǎo)致不收斂的問題。（2）省去了計(jì)算非線性方程組對應(yīng)的雅各比矩陣?yán)^而去求解每一次的迭代結(jié)果的步驟。Fsolve函數(shù)可以直接求解非線性方程組，計(jì)算速度相對較快。Fsolve函數(shù)求解電熱系統(tǒng)穩(wěn)態(tài)潮流的流程如圖3所示。

圖3 Fsolve函數(shù)求解電熱系統(tǒng)穩(wěn)態(tài)潮流流程示意圖

4 算例結(jié)果分析

4.1 穩(wěn)態(tài)潮流結(jié)果分析

運(yùn)用牛頓拉斐遜法對電熱約束方程來回迭代求解的算法，軟件為MATLAB，計(jì)算機(jī)參數(shù)為Intel core i5-8300H，2.30GHz，16GB，運(yùn)行時(shí)間為27.5s。求得的巴厘島電熱互聯(lián)系統(tǒng)的熱力參數(shù)和電力參數(shù)如表2和表3所示。其中，以文獻(xiàn)[7]中SINCAL軟件仿真值作為對照。

表2 所求熱網(wǎng)管道流量與實(shí)際管道流量數(shù)據(jù)對比

表3 電壓幅值和相角與實(shí)際值數(shù)據(jù)對比

從表2、表3所示電網(wǎng)和熱網(wǎng)數(shù)據(jù)來看，通過牛頓拉斐遜迭代法對電熱互聯(lián)的綜合能源系統(tǒng)的求解的電力熱力參數(shù)非常準(zhǔn)確，考慮電熱耦合的區(qū)域電熱互聯(lián)的綜合能源系統(tǒng)模型的建立與求解比較合理。

4.2 兩種求解算法結(jié)果比較分析

為了進(jìn)一步比較和研究兩種求解算法，預(yù)設(shè)了管道流量的三組初值用Fsolve函數(shù)的電熱互聯(lián)系統(tǒng)的穩(wěn)態(tài)潮流求解算法。管道流量的三組初值分別設(shè)定為1kg/s，-1kg/s和0kg/s，每組運(yùn)行時(shí)間約8s，結(jié)果如表4所示。其中，以文獻(xiàn)[7]中SINCAL軟件仿真準(zhǔn)確值作為對照。

表4 Fsolve函數(shù)計(jì)算的三組穩(wěn)態(tài)潮流結(jié)果對比

由表4數(shù)據(jù)可以看出，1 kg/s，-1 kg/s和0 kg/s分別作為初值用Fsolve函數(shù)求解時(shí)，0 kg/s作為初值得到了與文獻(xiàn)[7]中SINCAL軟件仿真準(zhǔn)確值相一致的仿真結(jié)果，而1 kg/s和-1 kg/s作為初值得到的結(jié)果與SINCAL軟件仿真準(zhǔn)確值不吻合。由前面計(jì)算結(jié)果，運(yùn)用牛頓拉斐遜法對電熱約束方程來回迭代求解的算法的單次運(yùn)行時(shí)間平均為27.5 s，而用Fsolve函數(shù)的電熱互聯(lián)系統(tǒng)的穩(wěn)態(tài)潮流求解算法單次運(yùn)行時(shí)間平均為3.9 s，可知Fsolve函數(shù)求解速度較快，但初值的選取對Fsolve求解結(jié)果影響較大。

5 結(jié) 論

本文對考慮電熱耦合的綜合能源系統(tǒng)進(jìn)行建模，將巴厘島電熱模型作為算例驗(yàn)證了模型的合理性，并利用牛頓拉斐遜法準(zhǔn)確求解出電熱綜合能源系統(tǒng)的穩(wěn)態(tài)潮流。同時(shí)，通過比較牛頓拉斐遜迭代法和Fsolve函數(shù)求解算法的優(yōu)劣，認(rèn)為牛頓拉斐遜迭代法不需要給定初值，但求解時(shí)間相對較長，F(xiàn)solve函數(shù)求解算法求解速度快，但對初始值的選取較為嚴(yán)格。