多時間尺度電熱綜合能源系統狀態估計研究

慈文斌，顧海飛，朱勁松

（1.國網山東省電力公司，山東 濟南 250000；2.東南大學電氣工程學院，江蘇 南京 210096；3.國網江蘇省電力有限公司，江蘇 南京 210096）

日益增長的能源危機和環境問題引發了構建高效、清潔能源結構的新需求[1-2]。作為能源變革的一種主要形式，電熱綜合能源系統通過能源梯級利用和協調優化將電力系統和熱力系統聯系在一起，從而促進可再生能源消納和運行靈活性[3]，并在工程實踐中得到廣泛應用[4-5]。因此，電熱綜合能源系統的統一建模和精確分析，對促進綜合能源系統發展十分必要。

狀態估計作為一種評估系統運行及預測系統發展的數值方法，為系統在線優化和控制提供理論指導。文獻[6]首先利用最小二乘法構建了熱力系統狀態估計模型，但由于僅考慮用戶側數據，其估計精度較低。文獻[7]通過構建基于全量測配置的電熱綜合能源系統狀態模型，有效提高了估計精度，但該方法抗擾動能力較弱。文獻[8]引入雙線性抗差方法用于解決電熱綜合能源系統中的非線性問題。文獻[9]考慮了多主體特性，實現了綜合能源系統的分布式狀態估計。文獻[10-11]將電熱聯合狀態估計與獨立估計對比，并指出聯合估計有利于提高量測冗余從而改善估計精度。

目前對于電熱綜合能源系統的技術應用已經開展了一定研究，如協調優化[12]、新能源消納[13]等，但這些研究均忽略了系統的動態特性。考慮到不同的運行時間尺度，在一定控制指令下，電力系統能迅速達到穩定狀態，而熱力系統仍處于動態過程[14]。若采用穩態模型對熱力系統進行建模，則忽略了熱力狀態分布的時空演變特性，將產生較大的估計誤差[15]。多時間尺度特性使得傳統的針對單一時間斷面的狀態估計方法不再適用于電熱綜合能源系統，由于熱力系統的狀態量在時序上相互影響，計及多時間尺度的狀態估計問題需研究一段時間內的狀態分布，在滿足各狀態量在不同控制時刻下的交互特性的前提下，最小化一段時間內量測量與狀態量之間的誤差，從而提高估計精度。

為此，本文構建了一種考慮多時間尺度的電熱綜合能源系統狀態估計方法，通過刻畫熱力系統狀態量的時序分布特性，對電熱綜合能源系統的狀態量進行準確估計，并利用拉格朗日乘子進行非線性優化問題的求解。最后，通過算例分析驗證了本文所提模型的有效性，并與穩態模型所得結果進行了對比。

1 電熱綜合能源系統動態建模

1.1 電力系統模型

電力系統主要關注節點和支路處的有功和無功分布，本文采用交流模型對其進行建模。節點注入功率平衡和支路傳輸功率平衡分別為[14]：

式中：Pi和Qi分別為節點i的注入有功和無功功率，Ui和θi分別為節點i的電壓幅值和相角，Gij和Bij分別為節點i與節點j之間的電導和電納，θij為節點i與節點j之間的相角差，Pl和Ql分別為線路傳輸的有功功率和無功功率。

1.2 熱力系統模型

熱力系統中，熱源產生的熱能通過工質流動在供/回水網中傳遞，因此，熱力系統模型可分為水力部分和熱力部分。本文主要關注運行在質調節的熱網，在此場景下，熱負荷需求通過調節在熱源處供水溫度得到滿足，并維持網絡中管道流量不變。由于水力狀態穩定且操作簡單，該調節方式在我國得到大量應用[9,15]。

熱力部分主要描述系統內的供/回水溫度和功率分布。根據換熱器模型，節點注入功率[13]為

式中：cp為水的比熱容，φ為節點注入功率，q為節點注入流量，Ts和Tr分別為節點供水和回水溫度。

混合后的節點溫度可表示為

式中：Tout為節點處的混合溫度，Tin為節點注入管道的末端溫度，mout為節點流出管道的流量，min為節點流入管道的流量。

管道溫度分布的動態過程一般通過常系數偏微分方程描述，即管道溫度分布與流速v、管道熱阻r以及環境溫度Ta相關[16]：

應用微分方程進行求解極大地增加了分析的復雜性，本文采用有限差分法對描述管道溫度分布動態特性的偏微分方程進行處理。首先，在區間0≤x≤L，0≤t≤th內將式(5)進行離散化，所得離散時空步長τ和h可分別表示為

式中：L和Γ分別表示管道長度和研究時間區間，N和M分別表示差分的時空步數。

在(xi,tk)處將式(5)的偏微分方程差分化，式(5)中各項偏微分可表示為：

式中O(h)和O(τ)分別為關于空間步長和時間步長的高階余項。

將式(7)和式(8)代入式(5)，則式(5)表示的管道溫度分布方程可轉換為：

1.3 耦合環節模型

電力系統和熱力系統通過耦合設備實現電能和熱能的雙向傳遞。常見設備包括熱電聯產(CHP)機組和電鍋爐(electric boiler，EB)等，電熱綜合能源系統耦合關系如圖1所示。

圖1 電熱綜合能源系統耦合關系示意Fig.1 Schematic diagram of the coupling relationship in the integrated electricity and heat system

由圖1 可見，電鍋爐通過消耗電能為熱力系統供熱，可將其視為熱源和等效電負荷；熱電聯產機組同時產生電能和熱能，可將其視為等效熱源和電源。其能量轉換關系為[15]

式中：?EB、?CHP分別為EB 或CHP 機組的輸出熱功率，PEB和PCHP分別為EB 和CHP 機組的輸出電功率，ηEB、ηCHP分別為EB 和CHP 機組的效率系數。

從耦合特性上看，電力和熱力系統主要通過CHP 機組和EB 耦合，耦合機組的電力和熱力輸出分別根據電力系統和熱力系統的運行時間尺度進行運行調整。因此，熱力系統的動態特性會通過耦合機組最終傳遞至電力系統，從而使得電熱綜合能源系統在整體上呈現動態特性[9-10,14-15,17]。

2 多時間尺度狀態估計

2.1 狀態量與量測量

電力系統中的量測量ze,t一般包括節點電壓幅值Ut,me、節點注入功率Pt,me和Qt,me以及支路傳輸功率P1,t,me和Q1,t,me；狀態量xe,t包括節點電壓幅值Ut和相角θt。熱力系統中量測量一般布置于熱源、用戶以及網絡中的部分管道，在質調節下，系統的水力狀態是確定的[7]，即可認為管道流量mme和節點壓力pme全可量測。在此基礎上，熱力系統中的量測量zh,t還包括供水溫度Ts,t,me和回水溫度Tr,t,me；狀態量xh,t包括節點供水溫度Ts,t和節點回水溫度Tr,t以及節點熱功率φt。電熱綜合能源系統中量測量和狀態量分布如式(12)所示，式中下標me 代表量測量。

2.2 多時間尺度特性

考慮系統運行的動態特性，在一定的控制指令下，電力系統能迅速達到穩定狀態而熱力系統仍然處于動態過程。因此，電力系統的調度周期通常較短，為分鐘級；而熱力系統調度周期相對較長，通常為幾十分鐘甚至數小時[9-10]。圖2 為電熱綜合能源系統耦合關系動態特性，圖2 中各網格為電力和熱力系統的調度周期。由圖2 可見，計及電力系統和熱力系統的運行時間尺度特性差異，在相同的調度周期內，電力系統調度周期較短、分布較密集；而熱力系統調度周期較長，分布較為稀疏。

圖2 電熱綜合能源系統耦合關系動態特性Fig.2 Dynamic characteristics of the coupling relationship in the integrated electricity and heat system

2.3 狀態估計模型

狀態估計模型的本質是：在一定約束條件下，最小化狀態量估計值和量測值的誤差。其目標函數可通過最小二乘法構建，表示為[7]：

式中：Je和Jh分別為研究區間內的電力系統和熱力系統的狀態估計誤差；Re和Rh分別為電力系統和熱力系統量測誤差的協方差矩陣，由歷史數據集確定；re和rh分別為電力系統和熱力系統中的各項狀態估計誤差；ξP,t、ξQ,t和ξU,t分別為t時刻節點有功功率、無功功率和電壓幅值的估計誤差；ξPl,t和ξQl,t分別為t時刻支路有功功率和無功功率的估計誤差；ξcp,t和ξ?,t分別為t時刻耦合機組功率和節點熱功率的估計誤差；ξT,t和ξTout,t分別為t時刻管道溫度方程和節點溫度混合方程的估計誤差；ξTs/r,t為t時刻節點供/回水溫度的估計誤差；Δte和Δth分別為電力和熱力系統的調控時間間隔；Γe和Γh分別為電力和熱力系統的調控時間標志，通常根據實際系統的運行特征作為常數給定[14-17]，其選擇不影響模型本質特性。

倘若所研究時間區間Γ為10 min，電力和熱力系統調控間隔Δte和Δth分別為1 min 和5 min，則根據式(13)，Γe和Γh分別取10、2 min。在給定的調控時刻集下，各狀態量的估計誤差可表示為：

式(14)和式(15)分別表示電力系統中由偽量測量構建的節點、支路功率估計誤差，其中Pi、Qi、Pl,i和Ql,i分別根據式(1)和式(2)構建。式(16)表示直接量測的節點電壓幅值的估計誤差[7]。式(17)反映了電力系統和熱力系統在不同時間尺度下，量測量和狀態量的耦合，由耦合設備運行特性決定。式(18)反映了熱力系統中的管道溫度的動態特性，即相鄰時間點管道溫度的相互關系，管道溫度分布的動態特性影響著耦合機組的輸出熱功率，進而通過式(17)依次影響著機組的輸出電功率與電力系統中的狀態量分布，其中，TI,k+1根據式(9)構建。式(19)和式(20)表示熱力系統中由偽量測量構建的節點熱功率、混合溫度的估計誤差。式(21)表示直接量測的節點供/回水溫度估計誤差[7-10]。

2.4 求解算法

由于式(14)和式(15)的非線性特征，電熱綜合能源系統狀態估計模型本質上是一個難以直接求解的非線性優化問題。因此，引入拉格朗日乘子對非線性約束進行處理。處理后的目標函數可表示為

式中，F(xe,xh)對應于式(14)—式(21)中的各個等式約束，ω為對應于各個約束的乘子系數向量。

處理后的模型為一種典型的非線性優化問題，而本文所提模型僅在優化模型中考慮了線性的多時間尺度耦合約束和本地量測約束，并不影響狀態估計模型和求解算法的自身特性。因此，可通過在MATLAB 下調用求解器IPOPT 進行直接求解，其可解性已得到充分驗證[7,9-10]。計算機配置環境為Intel Core i7、8GB RAM。

3 算例分析

本文通過2 個算例對所提模型進行分析：1）熱力系統動態模型精度驗證算例，用于分析所構建差分化熱力系統動態模型的計算精度，并且與實測數據、不計及動態特性的熱網模型進行對比驗證；2）電熱綜合能源系統狀態估計驗證算例用于分析所提多時間尺度狀態估計模型的準確性，并在不同的量測誤差下與不計及動態特性的熱網模型進行對比驗證，估計精度通過平均相對誤差衡量。

3.1 算例1

算例1 采用石家莊某供熱工程進行分析[13]，熱源為聯供機組，管道全長9.25 km，管道直徑為1.4 m，室外溫度為–10 ℃，熱阻系數為0.67 W/(m·℃)。熱源通過對供熱溫度進行調整來滿足熱負荷需求。此外，為驗證模型適用性，供水網中的管道流量是時變的，即構建差分式動態模型對于不同水力狀態分布的熱力系統同樣適用。

動態模型的差分時空步長影響建模精度，本文通過常用的仿真試探法進行選擇，因而在模型精度驗證和狀態估計問題中作為常數給定[14-15,17]。在熱源處添加階躍擾動，研究不同差分時空步長下的管道末端溫度的仿真結果，不同差分步長下的仿真精度如圖3所示。

圖3 不同差分步長下的仿真精度Fig.3 The simulation accuracy with various differential sizes

由圖3 可見：當差分的時間步長為6 min 時，不論差分的空間步長如何選擇，仿真結果均出現了不同程度的抖動；而當差分的時間步長為3 min 時，仿真結果的穩定性較高。因此，本算例中最終選擇動態模型的差分時間步長和空間步長分別為3 min和370 m。

未計及動態特性所得結果和本文所提動態模型的換熱站供水溫度計算結果對比如圖4所示。由圖4可見，動態模型準確地反映了供水溫度在一段時間內的變化趨勢，其平均相對誤差為0.27%。當不考慮動態特性時，系統中的狀態分布主要取決于實時的功率平衡且忽略了延時特性。因此，所得的溫度分布與動態模型相比提前了若干時刻。此外，不計動態特性時，模型所得結果完全由熱源處的供水溫度決定，其計算結果誤差較大，平均相對誤差為1.03%。

圖4 未計及動態特性模型和動態模型的供水溫度結果對比Fig.4 The supply water temperatures of the models with and without dynamic features

3.2 算例2

算例2 采用文獻[18]中電熱綜合能源系統進行分析，并對其中數據進行修正。算例2 系統結構如圖5所示。

圖5 算例2 系統結構Fig.5 The structure of IEHS in case 2

圖5 中，電力系統節點1 和節點2 為發電機組，節點6 為電熱耦合點，通過CHP 機組進行耦合。電力系統的調控間隔為1 min，熱力系統的調控間隔為5 min。與算例1 類似，動態模型中的差分時空步長通過提前測試以保證熱力系統動態模型具有較高的建模精度，本算例采用時間步長為1 min，每根管道分為25 段可獲得較高的仿真精度，并作為常數在狀態估計中給定。

算例中，以能流計算結果為真實值，在此基礎上疊加50 dB 的高斯白噪聲作為量測值進行分析，研究不同量測噪聲下的狀態估計結果從而驗證方法的通用性。研究時間區間為10:00—11:00，共1 h。

利用式(14)模型對研究區間進行狀態估計，所得電力系統和熱力系統估計結果分別如圖6—圖9所示。圖6 和圖8 分別為全時段內電力系統節點5電壓幅值和熱力系統節點3 的供水溫度時序分布；圖7 和圖9 分別為10:30 時刻電力系統和熱力系統中整體的節點電壓幅值和供水溫度估計結果分布。

圖6 和圖8 以及圖7 和圖9 分別反映了電力系統中的電壓幅值和熱力系統中供水溫度在時間和空間上的分布規律，即刻畫了電熱綜合能源系統內狀態量分布的時空相關性。

圖6 各時刻電力系統節點5 電壓幅值估計值Fig.6 The estimated U5 during the whole period

圖7 10:30 時刻電力系統電壓幅值估計值Fig.7 The estimated voltage amplitude at 10:30

圖8 各時刻熱力系統節點3 供水溫度估計值Fig.8 The estimated Ts3 during the whole period

圖9 10:30 時刻熱力系統供水溫度估計值Fig.9 The estimated supply water temperature at 10:30

由圖6—圖9 可見，本文所提方法進行狀態估計所得結果具有較高精度，且從整體上看，估計值的相對誤差小于量測值的相對誤差，從而驗證了在給定誤差的量測數據下，本文多時間尺度狀態估計模型的有效性。

此外，通過電力系統和熱力系統的結果對比可以看出：電力系統中所得到的狀態估計誤差明顯小于量測誤差；而在熱力系統中，在部分時刻下的估計誤差大于量測誤差。這是因為在多時間尺度狀態估計模型中，對于電力系統的分析本質上是對多時間斷面的研究，且各個時間斷面相互獨立；而熱力系統中的狀態分布相互影響，在較早時刻產生的估計誤差會通過管道溫度動態方程逐步向后續時刻傳遞，即形成了誤差傳遞鏈，從而影響了部分時刻的狀態估計精度。這是由熱力系統的物理特性決定的，但從整體上看，本文所提方法仍取得了較好的估計效果。

由圖8 可知，所建立的動態模型能有效反映熱力系統內溫度變化的動態過程，即各節點的供水溫度隨熱源處的溫度調節而變化，但由于前序時刻狀態差異，所體現的變化趨勢仍有部分偏差。將本文所提多時間尺度狀態估計模型與傳統的基于穩態模型所獲取的狀態估計結果進行對比，所得估計結果見表1。

表1 多時間尺度/穩態模型狀態估計平均相對誤差 %Tab.1 The mean relative error of state estimation results of the multi-timescale and dynamic model

由表1 可見：考慮電熱綜合能源系統的多時間尺度特性后，狀態估計的精度得到顯著提高。其改善主要在熱力系統部分，電力系統估計精度的改善相對較小。這是因為電力系統和熱力系統根據熱力系統調控尺度進行交互，熱力系統模型的變化僅在部分時刻對電力系統產生影響。此外，在算例2 中，熱力系統對電力系統的影響完全由聯供機組反映，而聯供機組在電力系統中可視為等效PV 節點，其功率變化對電力系統運行狀態影響較小，因而熱力系統模型精度的提高對于電力系統部分狀態估計精度的影響相對較小。對于熱力系統，由于動態模型利用了熱力系統不同時間、位置上量測數據的相關性，使得狀態估計所得結果精度得到明顯提高。而在質調節熱網中，由于節點注入流量是固定的，節點功率與溫度呈線性關系。因此，節點功率的估計誤差與溫度的估計誤差整體上近似。

對不同量測誤差對動態狀態估計結果的計算精度和計算時間的影響進行進一步研究。取量測誤差的噪聲分別為40、50、60 dB，不同量測誤差下，基于穩態模型狀態估計方法和動態狀態估計方法所得的計算誤差和計算時間對比見表2。

表2 不同量測誤差下的狀態估計誤差和計算時間結果對比Tab.2 The mean relative error and computation time with different measurement errors

由表2 可見，不同的量測誤差對于本文所提狀態估計方法計算時間差異的影響較小。這是因為狀態估計問題的計算時間主要由計算規模決定，而優化問題中所涉及的優化變量數目和約束條件數不隨量測誤差的增減而變化；而多時間尺度狀態估計所需的計算時間明顯多于穩態估計方法，其原因在于差分化的熱力系統動態模型中引入了較多的分段變量，使得優化問題規模顯著增大。在本案例中，多時間尺度狀態估計方法所需的計算時間約為穩態狀態估計方法的3～4 倍。

對比不同量測噪聲下的狀態估計結果，多時間尺度狀態估計方法和穩態狀態估計方法所得結果均隨著量測誤差的增大而增大，且穩態狀態估計結果精度受量測誤差的增大下降更為明顯。這是因為所研究系統規模較大，系統的延時特性十分明顯，而穩態狀態估計方法忽略了各個時間尺度下的耦合關系，使得時序的狀態估計過程中各個時刻點的估計結果較差；在多時間尺度方法中，除了網絡運行約束外額外在優化問題中添加了不同估計時刻之間的數值聯系，該處理方法等效在優化模型中添加了約束條件，使得多時間尺度狀態估計方法所得結果精度相對較高。

另一方面，由于量測項在優化目標中所對應的權重隨著量測誤差的增大而減小，因此，量測項在狀態估計問題中所產生的影響隨著量測誤差的增大而減小，使得狀態估計結果的誤差增長趨勢相對于量測誤差的變化趨勢較為緩慢。

4 結語

本文提出了一種多時間尺度電熱綜合能源系統狀態估計方法，利用電熱系統多時間尺度特性，分析系統狀態量的時空分布特征，進而研究了電熱綜合能源系統系統狀態演變的時序關聯性與交互機制。結合有限的量測數據，準確地感知了處于動態過程中的電熱綜合能源系統狀態演變趨勢，為實時控制和優化提供理論指導。

算例分析驗證了本文所提方法可準確描述熱力系統動態過程，與實測數據相比，其相對誤差為0.27%。通過與穩態狀態估計模型的對比，表明考慮系統多時間尺度特性可有效提高狀態估計精度，指出考慮系統動態特征進行建模分析的必要性，并從整體上驗證了本文所提方法在不同運行場景下的通用性。