張佳琛，郭慶來，王志偉，孫 勇，李寶聚，尹冠雄，孫宏斌

（1.清華大學(xué) 電機(jī)工程與應(yīng)用電子技術(shù)系，北京 100084；2.國網(wǎng)吉林省電力有限公司，吉林 長春 130012）

0 引言

傳統(tǒng)的電力系統(tǒng)、供熱系統(tǒng)、天然氣系統(tǒng)通常獨(dú)立規(guī)劃運(yùn)行，彼此間缺乏協(xié)調(diào)，存在能源利用效率低、故障跨系統(tǒng)蔓延等問題［1］。隨著現(xiàn)代工業(yè)生產(chǎn)中各類能源轉(zhuǎn)換設(shè)備的快速發(fā)展與部署，不同能源子系統(tǒng)間的耦合進(jìn)一步增強(qiáng)［2］。為了更高效地實現(xiàn)對能源系統(tǒng)的調(diào)度優(yōu)化，學(xué)者們提出了電、氣、熱、冷集成的綜合能源系統(tǒng)（integrated energy system，IES），能有力推動能源系統(tǒng)低碳化轉(zhuǎn)型［3］。

為了保證IES 的安全、可靠、優(yōu)質(zhì)、經(jīng)濟(jì)運(yùn)行，需要建立對多能流統(tǒng)一管理和科學(xué)調(diào)度的IES 能量管理系統(tǒng)。IES 能量管理系統(tǒng)的運(yùn)行依賴于對IES 全面、實時的精確感知，因此需要引入IES 狀態(tài)估計。IES狀態(tài)估計對原始量測數(shù)據(jù)進(jìn)行濾波，剔除量測數(shù)據(jù)中的壞數(shù)據(jù)，為IES 能量管理系統(tǒng)提供數(shù)據(jù)基礎(chǔ)。在IES 狀態(tài)估計中，熱網(wǎng)狀態(tài)轉(zhuǎn)移過程由偏微分方程（partial differential equation，PDE）描述，其時間常數(shù)遠(yuǎn)大于電力系統(tǒng)。因此，熱網(wǎng)動態(tài)狀態(tài)估計的主要難點(diǎn)在于PDE 實時求解。根據(jù)簡化方法的不同，衍生出了不同的模型驅(qū)動方法。文獻(xiàn)［4］使用節(jié)點(diǎn)法，將管道中的水流劃分為質(zhì)塊，以質(zhì)塊的流動來模擬管道的熱動態(tài)過程，而無須進(jìn)行空間上的差分。文獻(xiàn)［5-7］使用差分法，將管道內(nèi)的動態(tài)過程進(jìn)行不同格式的離散，將無窮維的函數(shù)空間的優(yōu)化轉(zhuǎn)變?yōu)橛邢蘧W(wǎng)格點(diǎn)中的優(yōu)化。文獻(xiàn)［8-9］使用能路法，將時域中的非線性關(guān)系轉(zhuǎn)換為頻域中的線性關(guān)系，顯著加快了計算速度?；谝陨虾喕Ｐ停墨I(xiàn)［10-11］類比電力系統(tǒng)傳統(tǒng)狀態(tài)估計方法，使用加權(quán)最小二乘法構(gòu)建了多能系統(tǒng)的優(yōu)化問題。文獻(xiàn)［12-13］將管道動態(tài)方程差分后得到一組線性方程，進(jìn)一步運(yùn)用卡爾曼濾波器對管道內(nèi)部溫度／流量進(jìn)行動態(tài)估計。已有的研究對熱網(wǎng)狀態(tài)估計做出了巨大的貢獻(xiàn)，但是仍然存在著以下的不足。

1）考慮熱網(wǎng)動態(tài)特性時必須引入PDE。依賴差分法的動態(tài)估計模型，計算精度取決于PDE 的差分格式和差分步長的選取。為了保證實時計算效率，差分步長往往較大，這將產(chǎn)生較大截斷誤差［5］。

2）在已有研究中，往往需要熱網(wǎng)提供全量測配置［14］，同時要求所有的量測完全同步［15］，否則需要進(jìn)行插值處理，這使得已有方法在實際的應(yīng)用中受到了限制。

3）模型驅(qū)動方法對參數(shù)依賴性較高。當(dāng)前關(guān)于熱網(wǎng)參數(shù)的不確定性和校準(zhǔn)方法的研究都非常有限［16］，所以當(dāng)熱網(wǎng)提供的模型參數(shù)不準(zhǔn)確或緩慢時變時，將顯著影響估計精度。

4）已有熱網(wǎng)狀態(tài)估計方法的壞數(shù)據(jù)辨識主要借鑒電力系統(tǒng)中的最大正則化殘差法［17］，并未考慮熱網(wǎng)慢動態(tài)特性下部分狀態(tài)變量的可預(yù)測性對壞數(shù)據(jù)辨識的作用。

神經(jīng)網(wǎng)絡(luò)等數(shù)據(jù)驅(qū)動方法具備從數(shù)據(jù)中發(fā)現(xiàn)信息的能力。文獻(xiàn)［18-19］利用多任務(wù)學(xué)習(xí)和長短時記憶網(wǎng)絡(luò)，提升了IES 負(fù)荷預(yù)測精度，且通過與已有模型結(jié)合，能得到更加有效的結(jié)果。文獻(xiàn)［20］提出了一種物理信息神經(jīng)網(wǎng)絡(luò)（physics-informed neural networks，PINNs），巧妙地將已有的PDE 模型融合到神經(jīng)網(wǎng)絡(luò)的訓(xùn)練之中，實現(xiàn)對PDE 的求解。實際物理系統(tǒng)因為數(shù)據(jù)獲取成本較大，能獲得的數(shù)據(jù)集往往較小，而PINNs 能夠用較小的數(shù)據(jù)集進(jìn)行高效訓(xùn)練，在實際應(yīng)用中具有優(yōu)勢。電力系統(tǒng)中發(fā)電機(jī)的運(yùn)行也通過微分方程來描述。文獻(xiàn)［21］首次將PINNs 應(yīng)用于電力系統(tǒng)，通過該方法成功預(yù)測了發(fā)電機(jī)暫態(tài)過程中的轉(zhuǎn)子角，并實現(xiàn)了對系統(tǒng)慣量和阻尼系數(shù)的識別，取得了良好的研究成果。隨后，文獻(xiàn)［22］在此基礎(chǔ)上進(jìn)一步構(gòu)建了基于PINNs 框架的機(jī)器學(xué)習(xí)方法，將其應(yīng)用于電力系統(tǒng)穩(wěn)定性分析。

本文結(jié)合已有熱網(wǎng)模型，提出了基于PINNs 的熱網(wǎng)動態(tài)狀態(tài)估計方法，本文的創(chuàng)新點(diǎn)歸納如下。

1）構(gòu)建了基于PINNs 的熱網(wǎng)狀態(tài)估計模型，將帶PDE 約束的優(yōu)化問題轉(zhuǎn)化為神經(jīng)網(wǎng)絡(luò)的訓(xùn)練問題。后者避免了對PDE 的簡化，能夠得到更精確的結(jié)果。同時PINNs 的損失函數(shù)也是待辨識參數(shù)的函數(shù)，將熱網(wǎng)參數(shù)辨識過程和損失函數(shù)優(yōu)化過程相融合，實現(xiàn)了對參數(shù)的實時辨識校正，降低了對模型參數(shù)的依賴性。

2）提出了一種滾動時域熱網(wǎng)動態(tài)狀態(tài)估計方法，實現(xiàn)相鄰時間窗的內(nèi)部變量初值和神經(jīng)網(wǎng)絡(luò)參數(shù)的信息傳遞，使基于PINNs 的熱網(wǎng)狀態(tài)估計模型能夠應(yīng)用于在線實時估計，并適應(yīng)不同量測配置與量測頻率，動態(tài)追蹤熱網(wǎng)狀態(tài)量的變化。

3）提出了一種基于PINNs 預(yù)測能力的高效壞數(shù)據(jù)辨識方法。通過擴(kuò)展PINNs的PDE 配置點(diǎn)的時間范圍，使得通過訓(xùn)練得到的模型具有對未來一段時間狀態(tài)的預(yù)測能力。通過對比預(yù)測值和量測值完成壞數(shù)據(jù)的辨識，提高了壞數(shù)據(jù)辨識的效率。

1 PINNs簡介

PINNs 是一種求解監(jiān)督學(xué)習(xí)任務(wù)的神經(jīng)網(wǎng)絡(luò)，可以作為非線性PDE 求解器。一般使用的神經(jīng)網(wǎng)絡(luò)通過大量的輸入和輸出數(shù)據(jù)來訓(xùn)練得到直接的映射，但是并不會用到某些領(lǐng)域中的先驗知識。PINNs 的核心思想在于將已有的PDE 模型編碼到學(xué)習(xí)算法中，這將增強(qiáng)神經(jīng)網(wǎng)絡(luò)從數(shù)據(jù)中獲取信息的能力，使其能夠更快逼近真實解并在訓(xùn)練數(shù)據(jù)較少的情況下具有良好的泛化性。

對于一般形式的PDE，可以用式（1）表示。

式中：t為時間變量；x為空間變量；u(t，x)為隱式的解；N[u(t，x)]表示對u(t，x)進(jìn)行的各種非線性運(yùn)算；Ω為空間位置的邊界；[0，Tw]為時間范圍。大多數(shù)的PDE 無法獲得解析的u(t，x)，可以利用神經(jīng)網(wǎng)絡(luò)U(t，x)這一通用函數(shù)近似器來逼近u(t，x)，U(t，x)的表達(dá)式如下：

式中：U(t，x)共有D層（不含輸入層）；“°”為復(fù)合函數(shù)算子；Lk(zk-1)（k=1，2，…，D）為仿射變換，將k-1層的輸出zk-1使用權(quán)重wk和偏置bk仿射變換為k層非線性激活函數(shù)σ的輸入，通過式（4）計算得到k層的輸出zk，而z0是輸入的(t，x)。激活函數(shù)的非線性賦予了神經(jīng)網(wǎng)絡(luò)強(qiáng)大的函數(shù)表達(dá)能力。

自然地，U(t，x)在近似u(t，x)的同時，也應(yīng)滿足對應(yīng)的PDE。由于U(t，x)具有顯式的函數(shù)表達(dá)式，可以獲得解析的對t和x的偏導(dǎo)數(shù)。在主流的機(jī)器學(xué)習(xí)庫如torch 庫中，神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)以計算圖的形式存儲，自動微分技術(shù)使得求取神經(jīng)網(wǎng)絡(luò)對輸入的偏導(dǎo)數(shù)變得非常便捷。因此，用U(t，x)近似u(t，x)時在PDE 中的失配量f( )t，x如式（5）所示，具體而言，f(t，x)可以通過torch庫中的autograd函數(shù)計算得到。

大多數(shù)情況下，u(t，x)不具有解析的表達(dá)式，所以只能通過差分得到數(shù)值導(dǎo)數(shù)。熱網(wǎng)動態(tài)狀態(tài)估計是帶PDE 約束的優(yōu)化問題，因此，常規(guī)做法需要使用差分法處理PDE 約束，得到大量需要求解的狀態(tài)量和線性方程約束，最后通過求解優(yōu)化問題得到狀態(tài)量的估計值。而PINNs 的巧妙之處在于使用U(t，x)近似u(t，x)后獲得了解析的導(dǎo)數(shù)，在不依賴差分的情況下天然內(nèi)嵌了PDE 約束，從而將帶PDE約束的優(yōu)化問題轉(zhuǎn)化為神經(jīng)網(wǎng)絡(luò)的訓(xùn)練問題，后者依賴于機(jī)器學(xué)習(xí)的優(yōu)化器，求解更為容易。

傳統(tǒng)差分法將偏導(dǎo)數(shù)轉(zhuǎn)化為數(shù)值導(dǎo)數(shù)。在本文應(yīng)用中，時間步長取值即為量測間隔，考慮熱網(wǎng)的實際量測配置，這一取值往往較大，同時為了保證差分穩(wěn)定性，空間步長取值也不會太小，這增加了數(shù)值方法計算的截斷誤差，影響整體計算精度。而本文方法使用了神經(jīng)網(wǎng)絡(luò)的解析導(dǎo)數(shù)，從而避免了這一截斷誤差，能獲得相對更高的精度。

PINNs 示意圖如圖1 所示。圖1 的左半部分是全連接神經(jīng)網(wǎng)絡(luò)，輸入為歸一化后的時間和空間坐標(biāo)，經(jīng)過由權(quán)重參數(shù)、偏置參數(shù)和非線性激活函數(shù)組成的多層隱藏層后，輸出得到神經(jīng)網(wǎng)絡(luò)的計算結(jié)果U(t，x)。

圖1 PINNs示意圖Fig.1 Schematic diagram of PINNs

而PINNs 與一般神經(jīng)網(wǎng)絡(luò)的主要差異在于損失函數(shù)的構(gòu)造。由于適定的PDE 需要給定邊界條件，所以神經(jīng)網(wǎng)絡(luò)U(t，x)的訓(xùn)練損失來源于兩部分：邊界采樣數(shù)據(jù)和PDE 約束。其參數(shù)可以通過最小化式（6）—（8）所示的均方誤差學(xué)習(xí)得到。

式中：L為總損失函數(shù)；δ為權(quán)重系數(shù)；為U(t，x)的邊界采樣數(shù)據(jù)，其中i=1，2，…，Nu，Nu為邊

界采樣總數(shù)；Lu為邊界采樣數(shù)據(jù)的訓(xùn)練誤差；為f(t，x)的PDE約束配置點(diǎn)，其中i=1，2，…，Nf，Nf為配置點(diǎn)總數(shù)；Lf為在有限個配置點(diǎn)處根據(jù)PDE 計算得到的失配量。文獻(xiàn)［23］指出PDE 中存在不合理的系數(shù)時，方程呈現(xiàn)病態(tài)。訓(xùn)練過程中，損失函數(shù)對不同的損失項的梯度差異較大，導(dǎo)致其對初值和邊界條件滿足情況較差。所以在式（6）中引入權(quán)重系數(shù)δ，將不同損失項對參數(shù)的梯度平衡到接近同一數(shù)量級，提升訓(xùn)練效果。

在將2種類型的損失項加權(quán)求和后，通過圖1下半部分所示的梯度下降方法完成對參數(shù)θ的訓(xùn)練更新。在torch 庫中內(nèi)置了各種參數(shù)更新的優(yōu)化器，本文使用其中的L-BFGS優(yōu)化器，學(xué)習(xí)率設(shè)置為1。

另外，PINNs 對噪聲數(shù)據(jù)輸入具有較強(qiáng)的魯棒性。一方面，由于訓(xùn)練神經(jīng)網(wǎng)絡(luò)的損失函數(shù)中包含邊值條件和PDE 殘差損失項，防止了網(wǎng)絡(luò)對噪聲的過擬合；另一方面，使用對神經(jīng)網(wǎng)絡(luò)的自動微分計算偏導(dǎo)數(shù)，而非使用含噪聲數(shù)據(jù)的數(shù)值導(dǎo)數(shù)，避免了傳統(tǒng)方法中低精度差分格式帶來的不良影響。

2 基于PINNs的熱網(wǎng)動態(tài)狀態(tài)估計

2.1 基于PINNs的管道等效建模

一個典型的5節(jié)點(diǎn)熱網(wǎng)拓?fù)浼捌鋮?shù)見附錄A，包括熱源（熱電聯(lián)產(chǎn)機(jī)組）、熱負(fù)荷（換熱站）、供水和回水管道，供水網(wǎng)絡(luò)和回水網(wǎng)絡(luò)對稱設(shè)計運(yùn)行。

熱網(wǎng)由水力系統(tǒng)和熱力系統(tǒng)構(gòu)成。目前廣泛采用的調(diào)度方法為質(zhì)調(diào)節(jié)，即在熱源處只改變網(wǎng)路的供水溫度，而網(wǎng)路的循環(huán)流量維持設(shè)計流量不變。同時，區(qū)域供熱網(wǎng)中壓力變化比溫度變化的傳播速度快了1 000 倍，在熱動態(tài)過程中，水力工況可近似認(rèn)為是穩(wěn)態(tài)。綜上考慮，水力系統(tǒng)和熱力系統(tǒng)相對解耦，可以分別進(jìn)行計算。因此本文的狀態(tài)估計僅針對熱力系統(tǒng)，并假設(shè)水力工況已知，在本文中考慮的量測配置對應(yīng)為節(jié)點(diǎn)的溫度量測和熱源及熱負(fù)荷的熱功率量測。

1）普通管道建模。

盡管神經(jīng)網(wǎng)絡(luò)有著強(qiáng)大的模型表達(dá)能力，但是如果用一個神經(jīng)網(wǎng)絡(luò)構(gòu)建整個熱網(wǎng)量測量和狀態(tài)量間的映射關(guān)系，則模型參數(shù)和約束都會過多，進(jìn)一步引發(fā)收斂性和實時性的問題，所以為熱網(wǎng)中每條管道構(gòu)建各自對應(yīng)的PINNs模型。

在忽略鄰近流體間熱傳導(dǎo)項的情況下，一維動態(tài)熱網(wǎng)管道的溫度時空分布物理特性可以用如下PDE表示［24］：

式中：c為水的比熱容；ρ為水的密度；A為管道橫截面積；m為管道流量；T(t，x)為t時刻距管道首端x處的凈溫度，表示水溫與環(huán)境溫度差值；μ為管道散熱系數(shù)。

對于非熱源或換熱站等效的普通管道，用深度神經(jīng)網(wǎng)絡(luò)u(t，x)近似T(t，x)，即式（9）的解。同時，定義f(t，x)為給定方程式（9）等號左側(cè)的式子：

2）等效管道建模。

對于連接供水網(wǎng)絡(luò)及回水網(wǎng)絡(luò)的熱源和換熱站，為保證拓?fù)浞治龅囊恢滦?，將其等效為管道支路處理。在熱源和換熱站中存在熱功率量測，與溫度的對應(yīng)關(guān)系為：

式中：Φp(t)表示t時刻等效為管道p的熱源或換熱站的熱功率，其值大于0 表示消耗熱功率，反之則表示產(chǎn)生熱功率；Tp(t，0)和Tp(t，1)分別為t時刻管道p在首端和末端的凈溫度。下文中，上標(biāo)p表示管道p對應(yīng)的相關(guān)變量。需要注意的是，在熱功率計算中使用的溫度為熱源或換熱站等效管道的首末端溫度，而非節(jié)點(diǎn)溫度。

熱源和換熱站等效管道并沒有PDE 的約束，只有式（11）所示的線性等式約束，所以可以用較小的神經(jīng)網(wǎng)絡(luò)u(t，x)，結(jié)合溫度量測及熱功率量測，近似時間窗內(nèi)換熱站首末端的溫度變化。u(t，x)中t的取值同普通管道，但是由于等效管道不具有實際物理長度，x有意義的取值設(shè)定為0 和1，u( )t，x在這2 個位置的取值分別對應(yīng)進(jìn)入、離開熱源或換熱站后的溫度。

2.2 PINNs損失函數(shù)構(gòu)建

在2.1 節(jié)中建立的模型僅針對單個熱網(wǎng)管道或等效管道，而整個供熱系統(tǒng)還存在其他用靜態(tài)方程描述的約束，包括節(jié)點(diǎn)溫度方程式（12）和匯流混合方程式（13）。類比式（6）—（8）構(gòu)建的損失函數(shù)，這些靜態(tài)方程也可以以損失函數(shù)的形式加入PINNs 的訓(xùn)練過程中，從而保證整個網(wǎng)絡(luò)估計結(jié)果的合理性。

式中：N為熱網(wǎng)節(jié)點(diǎn)集合；為從節(jié)點(diǎn)n流出的管道集合；Tn(t)為t時刻節(jié)點(diǎn)n溫度，當(dāng)節(jié)點(diǎn)存在量測時，溫度也取該值為從節(jié)點(diǎn)n流入的管道集合；xe為管道長度。式（12）表示某節(jié)點(diǎn)流出管道首端溫度等于該節(jié)點(diǎn)溫度；式（13）表示當(dāng)存在匯流節(jié)點(diǎn)時，節(jié)點(diǎn)溫度等于匯流管道末端水流混合后的溫度。式（13）等號左側(cè)涉及多根管道的PINNs 模型，為了充分利用這部分量測數(shù)據(jù)，這些網(wǎng)絡(luò)將共享損失函數(shù)。因此根據(jù)網(wǎng)絡(luò)的拓?fù)浣Y(jié)構(gòu)，訓(xùn)練模式也分為對單個神經(jīng)網(wǎng)絡(luò)的訓(xùn)練和對多個神經(jīng)網(wǎng)絡(luò)的聯(lián)合訓(xùn)練。

在拓?fù)渥R別中，首先識別有多根管道輸入的匯流節(jié)點(diǎn)，若匯流節(jié)點(diǎn)存在量測，則所有輸入該節(jié)點(diǎn)的管道的神經(jīng)網(wǎng)絡(luò)將進(jìn)行聯(lián)合訓(xùn)練，其余管道的神經(jīng)網(wǎng)絡(luò)將進(jìn)行單獨(dú)訓(xùn)練。

1）單獨(dú)訓(xùn)練。

單獨(dú)訓(xùn)練的損失函數(shù)由邊值條件損失項（Lu）、PDE損失項（Lf）和溫度差量測損失項（Lh）組成。

2）聯(lián)合訓(xùn)練。

聯(lián)合訓(xùn)練的損失函數(shù)由邊值條件損失項（Lu）、PDE損失項（Lf）、溫度量測損失項（Lh）和匯流方程損失項（Lk）組成。

在熱網(wǎng)匯流節(jié)點(diǎn)處需要引入額外的約束。在匯流節(jié)點(diǎn)測得的溫度不等于管道末端的溫度，而是多根管道水流匯合后得到的溫度。使用該點(diǎn)測得的數(shù)據(jù)來對流入管道的溫度分布進(jìn)行估計，應(yīng)該利用式（13）這一代數(shù)方程約束對多個網(wǎng)絡(luò)同時進(jìn)行訓(xùn)練。對于匯流節(jié)點(diǎn)n所涉及的所有流入管道，它們除了自身的均方誤差之外，還將共享利用匯流混合約束得到的均方誤差，即Lk。

式中：tik為匯流節(jié)點(diǎn)溫度量測時刻數(shù)據(jù)，i=1，2，…，Nk，Nk為溫度量測時刻總數(shù)；μk為Lk對應(yīng)的權(quán)重系數(shù)。

傳統(tǒng)狀態(tài)估計方法要求所有數(shù)據(jù)都有統(tǒng)一的時間坐標(biāo)，對于時間不一致的量測數(shù)據(jù)，需要通過插值獲得統(tǒng)一時間坐標(biāo)的結(jié)果。而從Lu的計算式式（14）可以看出，對使用的樣本數(shù)據(jù)ui并沒有時間坐標(biāo)統(tǒng)一的要求。由于PINNs 訓(xùn)練過程的本質(zhì)是擬合最符合量測數(shù)據(jù)的連續(xù)函數(shù)，所以不需要提供同步的量測值，異步的量測數(shù)據(jù)對模型都有著修正的作用。這對于實際中量測時間或間隔不一致、量測延遲等問題有著很好的適用性。

對損失函數(shù)做進(jìn)一步的討論，可以發(fā)現(xiàn)用PINNs 近似狀態(tài)變量時，Lu相當(dāng)于原狀態(tài)估計問題的優(yōu)化目標(biāo)，Lf相當(dāng)于PDE 約束，Lh和Lk相當(dāng)于其他偽量測約束。PINNs 的訓(xùn)練過程實際上是對這一優(yōu)化問題的另一種求解思路。

2.3 基于PINNs的熱網(wǎng)參數(shù)辨識方法

當(dāng)PDE 中含有未知參數(shù)時，通過觀測數(shù)據(jù)進(jìn)行參數(shù)反演被稱為PDE 的逆問題。逆問題比正問題的求解往往更為困難，因為它要從有限的數(shù)據(jù)中直接推斷方程的信息。微分方程參數(shù)反演的基本數(shù)值方法為構(gòu)建待反演參數(shù)的最小二乘泛函，結(jié)合正則化技術(shù)，通過極小化來逼近原問題的解。但是傳統(tǒng)方法存在的主要問題包括：正則化參數(shù)對擬合效果的影響極大，選擇合適的正則化參數(shù)較為困難；求解該問題需要迭代優(yōu)化，計算精度要求較高時產(chǎn)生的大量網(wǎng)格導(dǎo)致大量的計算開銷；僅對特定逆問題的求解有效，并不是普適方法。

PINNs 為逆問題的求解提供了新的思路，它更自然地將數(shù)據(jù)融合到求解過程中，對PDE 正問題和逆問題的求解在形式上并無太大差別，適用范圍更廣，這正是其巨大優(yōu)勢之一。PINNs 在含噪聲觀測數(shù)據(jù)中依然能夠發(fā)現(xiàn)最能擬合觀測數(shù)據(jù)的PDE 參數(shù)，同時獲得利用神經(jīng)網(wǎng)絡(luò)得到的近似的解。含有未知參數(shù)的PDE可以表示為：

式中：λ為f的未知參數(shù)。類似神經(jīng)網(wǎng)絡(luò)中的權(quán)重和偏置參數(shù)的訓(xùn)練方式，同樣可以通過自動微分獲得損失函數(shù)對λ的梯度，進(jìn)而能通過梯度下降等方法訓(xùn)練得到估計結(jié)果。文獻(xiàn)［20］指出采用PINNs 進(jìn)行參數(shù)辨識時，其對噪聲干擾有很好的魯棒性，相較于高斯過程回歸和稀疏回歸有更好的效果。

3 基于PINNs的滾動時域動態(tài)狀態(tài)估計

文獻(xiàn)［20］提出的PINNs 只是作為定邊界條件下的PDE 求解器來使用。熱網(wǎng)運(yùn)行過程中節(jié)點(diǎn)溫度都在實時發(fā)生變化，PINNs 訓(xùn)練的邊界條件也在不斷發(fā)生變化，模型需要實時訓(xùn)練。為了使PINNs 能夠追蹤熱網(wǎng)狀態(tài)量的變化，下面提出基于PINNs 的滾動時域動態(tài)狀態(tài)估計方法。

定義時間窗為當(dāng)前動態(tài)狀態(tài)估計計算的時間范圍。當(dāng)出現(xiàn)新的量測時，時間窗需要不斷向前滾動。定義時間窗窗寬為式（1）中的最長時間Tw，時間窗間隔為ΔTw。

在時間窗滾動的過程中，本次估計在中間某時刻的估計值將作為初值提供給下一次的估計，從而保證了估計的連續(xù)性。

式中：Tn(0，x)為第n個時間窗0時刻距管道首端x處的邊界采樣數(shù)據(jù)；un-1(ΔTw，x)為第n-1 個時間窗ΔTw時刻距管道首端x處的狀態(tài)估計結(jié)果。

在熱網(wǎng)量測有限的情況下，動態(tài)狀態(tài)估計的冗余度體現(xiàn)在時間窗窗寬Tw的選取上。隨著Tw的增加，歷史量測斷面逐漸增多，通過PDE 產(chǎn)生的時空耦合約束也隨之增多，相應(yīng)的測量冗余度也會提高，計算效果更佳。而時間窗窗寬Tw的選取涉及權(quán)衡取舍的問題。一方面，時間窗太短會導(dǎo)致時間窗中訓(xùn)練數(shù)據(jù)不足，難以反映PDE 所描述的熱動態(tài)特性，影響計算精度；另一方面，時間窗太長會導(dǎo)致訓(xùn)練代價增加，而且在時間窗內(nèi)可能出現(xiàn)的參數(shù)變化會導(dǎo)致參數(shù)辨識精度降低。所以需要在保證估計精度的情況下，根據(jù)熱網(wǎng)運(yùn)行周期選擇盡可能短的時間窗。

借鑒遷移學(xué)習(xí)中的經(jīng)驗，PINNs 的訓(xùn)練同樣也不需要在時間窗滾動后隨機(jī)初始化全部參數(shù)重新訓(xùn)練。由于熱網(wǎng)溫度變化較慢，所以相鄰時間窗的邊值條件相差較小。當(dāng)前時間窗網(wǎng)絡(luò)參數(shù)初始化過程中，直接使用上一個時間窗訓(xùn)練得到的網(wǎng)絡(luò)參數(shù)，將顯著加速PINNs 的實時訓(xùn)練過程。滾動時間窗PINNs狀態(tài)估計示意圖如圖2所示。

圖2 滾動時間窗PINNs狀態(tài)估計示意圖Fig.2 Schematic diagram of rolling time window PINNs state estimation

4 基于PINNs的狀態(tài)預(yù)測及壞數(shù)據(jù)辨識

4.1 基于PINNs的狀態(tài)預(yù)測

神 經(jīng) 網(wǎng) 絡(luò)u(t，x) 的 默 認(rèn) 輸 入 范 圍 是x∈Ω，t∈[0，Tw]。但在實際過程中，由于熱慣性的存在，當(dāng)管道首端溫度發(fā)生變化時，并不會即時影響末端的溫度變化。換言之，通過添加Tw＜tf＜ξTw(ξ＞1)的PDE 約束配置點(diǎn)訓(xùn)練得到的模型，當(dāng)u(t，x)中的輸入為Tw＜t＜ξTw時，仍然能得到準(zhǔn)確的結(jié)果。

溫度變化的傳遞伴隨著熱水質(zhì)塊的流動，所以PINNs 能預(yù)測的時間范圍可以由流速大致估計得到。定義管道的特征時間函數(shù)，即熱水從管道首端流到距首端x處的時間t(x)為：

為了使模型具有預(yù)測能力，PDE 約束點(diǎn)取值范圍變?yōu)閤f∈Ω，tf∈[0，Tw+γt(xe)]，其中γ＞1，為保證計算精度而留有一定安全裕度。

4.2 基于PINNs的壞數(shù)據(jù)辨識

狀態(tài)估計的另一個重要作用是辨識壞數(shù)據(jù)?？紤]實際過程中容易出現(xiàn)的2 類壞數(shù)據(jù)：由于儀表故障、通信干擾等導(dǎo)致的較大的量測誤差；由于數(shù)據(jù)丟失導(dǎo)致量測值為0或丟失。

在4.1 節(jié)中，PINNs 方法通過對PDE 約束點(diǎn)的配置，在當(dāng)前時間窗實現(xiàn)了對下一時間窗狀態(tài)量的預(yù)測。當(dāng)下一時間窗的量測與基于這一預(yù)測步的結(jié)果之間存在較大的差異時，可據(jù)此完成壞數(shù)據(jù)的辨識。該壞數(shù)據(jù)辨識算法基本示意圖及流程圖見附錄B。

在傳統(tǒng)電力系統(tǒng)狀態(tài)估計方法中，先進(jìn)行狀態(tài)估計，再利用正則化殘差法剔除壞數(shù)據(jù)。循環(huán)以上過程，直到?jīng)]有新的壞數(shù)據(jù)產(chǎn)生或達(dá)到迭代次數(shù)上限。由于本文方法對新的時間窗內(nèi)量測的變化已經(jīng)有了一定的估計，可以設(shè)定一個閾值，將預(yù)測和量測之間的絕對值誤差超過閾值的量測判定為壞數(shù)據(jù)，在訓(xùn)練前就進(jìn)行剔除，這樣避免了多次狀態(tài)估計的重復(fù)計算，同時對多種壞數(shù)據(jù)的形式都有較好的濾除作用。閾值的選取可以根據(jù)概率分布的參數(shù)選擇。假設(shè)噪聲服從高斯分布時，可以根據(jù)歷史經(jīng)驗確定噪聲大致分布的方差，進(jìn)一步根據(jù)期望置信度得到的分位數(shù)確定門檻值。

但這種預(yù)測能力僅限于流入管道不含熱源／換熱站等效支路的節(jié)點(diǎn)，這類節(jié)點(diǎn)量測定義為Ⅰ類節(jié)點(diǎn)量測，其他節(jié)點(diǎn)量測定義為Ⅱ類節(jié)點(diǎn)量測。Ⅰ類節(jié)點(diǎn)在新的時間窗內(nèi)的預(yù)測值可以用上一個時間窗內(nèi)得到的PINNs 模型直接計算得到；Ⅱ類節(jié)點(diǎn)預(yù)測通過等效支路相連的Ⅰ類節(jié)點(diǎn)預(yù)測疊加熱源／換熱站的熱功率量測得到。由于熱功率量測中可能含有噪聲，所以Ⅱ類節(jié)點(diǎn)對于壞數(shù)據(jù)的判定范圍也會相應(yīng)增大。

5 算例分析

5.1 算例設(shè)置及計算指標(biāo)

下面首先使用5 節(jié)點(diǎn)算例驗證本文方法在動態(tài)狀態(tài)估計、參數(shù)辨識及壞數(shù)據(jù)辨識方面的有效性，5節(jié)點(diǎn)算例拓?fù)浼皡?shù)見附錄A。算例數(shù)據(jù)和神經(jīng)網(wǎng)絡(luò)的參數(shù)設(shè)定見附錄C。同時，使用量測誤差統(tǒng)計值和估計誤差統(tǒng)計值來表征狀態(tài)估計的性能。

量測誤差統(tǒng)計值可以表示為：

估計誤差統(tǒng)計值可以表示為：

式中：N為蒙特卡羅仿真次數(shù)；l為時間窗內(nèi)量測的數(shù)量；分別為第n*次仿真中第i個量測值、估計值和真值；σi為第i個量測值量測誤差標(biāo)準(zhǔn)差。

由于假定量測噪聲服從高斯分布，所以SM的值接近于1。因此，可以用SH/SM來評價狀態(tài)估計性能，該值小于1 表示進(jìn)行狀態(tài)估計后效果優(yōu)于不進(jìn)行狀態(tài)估計，并且該值越小，狀態(tài)估計效果越好。

本文仿真程序均基于Python 編寫，處理器為Intel?CoreTMi5-8279U@2.4 GHz，8 GB內(nèi)存。

5.2 狀態(tài)估計方法有效性驗證

考慮非全量測配置情況，在節(jié)點(diǎn)0、1、3 — 5、8、9引入溫度量測；在所有熱源／換熱站等效支路設(shè)置熱功率（溫度差）量測。通過拓?fù)浞治龅玫降墓烙嫹绞綖椋汗艿?、3、4、6 — 9 單獨(dú)訓(xùn)練；管道1 和2、5 和11、10 和14、12 和13 聯(lián)合訓(xùn)練。同時考慮量測并不同步，節(jié)點(diǎn)1 量測間隔設(shè)為80 s，節(jié)點(diǎn)4 量測間隔設(shè)為90 s，其余節(jié)點(diǎn)量測間隔設(shè)為60 s。溫度量測中加入了標(biāo)準(zhǔn)差為0.35 ℃ 的高斯噪聲，溫度差量測中加入了標(biāo)準(zhǔn)差為0.14 ℃ 的高斯噪聲。

圖3 展示了經(jīng)過100 次蒙特卡羅仿真后的不同時間窗的SM、SH計算結(jié)果。全時間段SH/SM的平均值為0.124 2，同時相較于差分法的狀態(tài)估計結(jié)果，本文方法有著更低的SH，驗證了本文方法具有很好的狀態(tài)估計精度。另一方面，本文方法在不同時間窗計算得到的SH一直在較低的值附近波動，沒有出現(xiàn)較大的偏移，驗證了本文方法在不同樣本數(shù)據(jù)下能夠得到相似的狀態(tài)估計精度。

圖3 蒙特卡羅仿真誤差統(tǒng)計值計算結(jié)果Fig.3 Calculation results of Monte Carlo simulation error statistics

基于本文方法得到的部分估計結(jié)果與含噪聲量測和真值之間的對比如圖4 所示?？梢姳疚姆椒ǖ玫降墓烙嫿Y(jié)果與真值基本一致。而且在非全量測配置的情況下，節(jié)點(diǎn)2、6、7 的溫度不能直接通過量測得到，本文方法可以根據(jù)歷史數(shù)據(jù)并結(jié)合熱網(wǎng)模型，準(zhǔn)確推斷得到非量測配置點(diǎn)的溫度動態(tài)。

圖4 部分節(jié)點(diǎn)狀態(tài)估計結(jié)果Fig.4 State estimation results of certain nodes

本文方法同時提供了管道內(nèi)部任意時空的溫度估計。圖5 展示了管道0 內(nèi)部全時空分布的溫度估計結(jié)果與真實分布的對比，估計絕對誤差最大值為0.208 22 ℃，平均值為-0.000 56 ℃，可以發(fā)現(xiàn)本文方法對內(nèi)部不可觀變量有著非常準(zhǔn)確的估計結(jié)果。

圖5 管道0內(nèi)部溫度估計結(jié)果Fig.5 State estimation results of internal temperature in Pipeline 0

在計算效率方面，由于滑動時間窗之間參數(shù)的傳遞，相比于全隨機(jī)初始化效率更高。平均每次時間窗估計時間約為10 s，遠(yuǎn)小于時間窗之間的間隔，能適用于在線應(yīng)用的需要。并且表征不同管道的神經(jīng)網(wǎng)絡(luò)在沒有偽量測生成的依賴時，能夠并行訓(xùn)練，計算的時間將進(jìn)一步降低。

5.3 熱網(wǎng)不精確參數(shù)辨識及狀態(tài)估計結(jié)果

1）恒定參數(shù)辨識。

仍然使用5.2 節(jié)中的5 節(jié)點(diǎn)算例，溫度量測中加入了標(biāo)準(zhǔn)差為0.10 ℃ 的高斯噪聲，溫度差量測中加入了標(biāo)準(zhǔn)差為0.10 ℃ 的高斯噪聲。首先指定辨識參數(shù)為各管道中的流量。由于水力工況已知，所以此處流量被視作熱網(wǎng)動態(tài)狀態(tài)估計的模型參數(shù)而非狀態(tài)量。

在參數(shù)估計中一方面需要更多的數(shù)據(jù)對參數(shù)進(jìn)行擬合，所以將估計的時間窗從原來的2 h 延長為4 h；另一方面需要神經(jīng)網(wǎng)絡(luò)有更強(qiáng)的函數(shù)逼近能力，所以在神經(jīng)網(wǎng)絡(luò)層數(shù)不變的情況下適當(dāng)增加了每層神經(jīng)元的個數(shù)，設(shè)置為每層有15 個神經(jīng)元。網(wǎng)絡(luò)模型使用的流量初值均指定為100 kg／s。

下面對比在模型參數(shù)不精確的情況下，本文方法與基于差分法的動態(tài)狀態(tài)估計結(jié)果。以供水網(wǎng)絡(luò)中的管道為例，差分法計算中給定的流量參數(shù)均為真值的70 %。

部分節(jié)點(diǎn)的估計結(jié)果如圖6 所示?？梢姳疚姆椒ㄔ谔峁﹨?shù)不準(zhǔn)確的情況下較好地追蹤了熱網(wǎng)狀態(tài)量的變化。而差分法則由于對參數(shù)的強(qiáng)依賴性，給定流量參數(shù)比真值偏小時，整體估計結(jié)果向右偏移，計算結(jié)果出現(xiàn)較為明顯的偏差。

圖6 不精確參數(shù)下狀態(tài)估計結(jié)果對比Fig.6 Comparison of state estimation results under imprecise parameters

恒定參數(shù)辨識結(jié)果如圖7 所示?？梢钥闯黾词怪付ú痪_參數(shù)初值離真值較遠(yuǎn)，本文方法依然較好地追蹤到了網(wǎng)絡(luò)參數(shù)。這說明本文方法能夠自適應(yīng)修正網(wǎng)絡(luò)參數(shù)，相較差分法可以得到更好的估計結(jié)果。

圖7 恒定參數(shù)辨識結(jié)果Fig.7 Identification results of constant parameter

2）時變參數(shù)辨識。

本文方法也可適用于辨識緩慢時變的參數(shù)。辨識結(jié)果為時變參數(shù)在時間窗內(nèi)變化的平均結(jié)果。以單根管道為例，取5 節(jié)點(diǎn)熱網(wǎng)中的管道0 估計，管道流量緩慢時變。

時變參數(shù)辨識結(jié)果如圖8 所示?？梢姳疚姆椒ㄒ廊惠^為準(zhǔn)確地追蹤了參數(shù)的變化，最大相對誤差為7.06 %，平均相對誤差為0.49 %。從信息的角度出發(fā)，時變參數(shù)在分立時刻提供了不同值對應(yīng)的信息，要從觀測數(shù)據(jù)反演每個時刻時變參數(shù)的準(zhǔn)確值幾乎是不可能的，所以本文方法提供了相對折中的方案。

圖8 時變參數(shù)辨識結(jié)果Fig.8 Identification results of time-varying parameters

5.4 壞數(shù)據(jù)辨識結(jié)果

壞數(shù)據(jù)辨識仍然使用5.2 節(jié)中的5 節(jié)點(diǎn)算例。在每個節(jié)點(diǎn)全時間段（00:00 — 24:00）共1 440 個量測數(shù)據(jù)中，隨機(jī)選擇20 個數(shù)據(jù)更改為壞數(shù)據(jù)。壞數(shù)據(jù)中有1/3 的概率為0 值量測（表示丟失量測）；有1/3 的概率額外引入了1.75～2.8 ℃ 的較大噪聲；有1/3 的概率額外引入了-2.8～-1.75 ℃ 的較大噪聲。壞數(shù)據(jù)辨識閾值設(shè)定為1.4 ℃，部分溫度量測壞數(shù)據(jù)辨識結(jié)果如圖9所示。

圖9 部分溫度量測壞數(shù)據(jù)辨識結(jié)果Fig.9 Bad data identification results of certain temperature measurements

進(jìn)一步，定義壞數(shù)據(jù)檢測指標(biāo)如下：

式中：P為查準(zhǔn)率；NTP為所辨識壞數(shù)據(jù)是預(yù)設(shè)壞數(shù)據(jù)的數(shù)量；NFP為所辨識壞數(shù)據(jù)不是預(yù)設(shè)壞數(shù)據(jù)的數(shù)量；R為查全率；NFN為未能辨識的預(yù)設(shè)壞數(shù)據(jù)數(shù)量。查準(zhǔn)率表示判定正確結(jié)果占判定結(jié)果的比例，而查全率表示判定正確結(jié)果占預(yù)設(shè)結(jié)果的比例。

本文方法在所有等效節(jié)點(diǎn)辨識結(jié)果如表1 所示?？梢姳疚姆椒ㄡ槍臄?shù)據(jù)辨識的查準(zhǔn)率和查全率都有很好的結(jié)果。事實上，額外判定的壞數(shù)據(jù)也往往與真值相差較大，可以認(rèn)為是量測質(zhì)量較低的數(shù)據(jù)，在壞數(shù)據(jù)檢測部分剔除反而可以提升估計的精度。

表1 壞數(shù)據(jù)辨識對應(yīng)指標(biāo)計算結(jié)果Table 1 Calculation results of indicators of bad data identification

為了進(jìn)一步驗證本文方法對壞數(shù)據(jù)的辨識能力，更改上述算例中的壞數(shù)據(jù)占比，并進(jìn)行100 次蒙特卡羅仿真，辨識結(jié)果如表2 所示。可以看出，隨著壞數(shù)據(jù)占比的增加，查準(zhǔn)率變化幅度較小，保證了壞數(shù)據(jù)辨識的正確性。另一方面，隨著壞數(shù)據(jù)占比的增加，查全率有所下降，這是由于壞數(shù)據(jù)占比的增加減少了用于估計的合理數(shù)據(jù)的數(shù)量，導(dǎo)致估計效果下降，壞數(shù)據(jù)辨識范圍更加模糊。整體而言，SM隨著壞數(shù)據(jù)占比增加顯著增大，而SH依然保持在較低的水平，這說明在壞數(shù)據(jù)占比較高的情況下，本文方法依然有效過濾了量測噪聲，保持著狀態(tài)估計的精度。

表2 不同壞數(shù)據(jù)占比下的壞數(shù)據(jù)辨識結(jié)果Table 2 Bad data identification results under different bad data proportion

27節(jié)點(diǎn)熱網(wǎng)測試算例結(jié)果見附錄D。

6 結(jié)論

熱網(wǎng)作為IES 的重要組成部分，其狀態(tài)估計結(jié)果的準(zhǔn)確性將影響整個系統(tǒng)后續(xù)依賴熟數(shù)據(jù)的高級應(yīng)用部分。本文構(gòu)建了一種模型-數(shù)據(jù)聯(lián)合驅(qū)動的熱網(wǎng)狀態(tài)估計方法，利用PINNs 和熱網(wǎng)管道PDE 模型，建立了熱網(wǎng)的滾動時域狀態(tài)估計模型，準(zhǔn)確感知管道內(nèi)部溫度變化情況。在此基礎(chǔ)上，將其應(yīng)用于熱網(wǎng)的狀態(tài)估計中，并且得到了較為準(zhǔn)確的結(jié)果。同時驗證了本文方法能夠利用歷史數(shù)據(jù)完成對熱網(wǎng)管道參數(shù)的估計，能在運(yùn)行過程中不斷對模型內(nèi)部參數(shù)進(jìn)行修正，提升模型估計的精度。最后利用PINNs 的預(yù)測能力，實現(xiàn)了溫度量測中對壞數(shù)據(jù)的更高效辨識。

實際上，在IES 中，同樣極為重要的天然氣系統(tǒng)管道運(yùn)行中的壓力和流量之間的動態(tài)關(guān)系也可以通過PDE 描述。因此，可以進(jìn)一步擴(kuò)展并適當(dāng)修改這一方法，將其應(yīng)用于天然氣系統(tǒng)的動態(tài)狀態(tài)估計中，從而為不同類型的能源系統(tǒng)提供更準(zhǔn)確的狀態(tài)估計，有助于優(yōu)化能源系統(tǒng)的運(yùn)行和管理。

未來的研究工作中，筆者將進(jìn)一步在實際更大規(guī)模算例中對本文方法進(jìn)行驗證，以探究系統(tǒng)規(guī)模對計算效率的影響。并研究新的物理融合的模式，使本文方法能適應(yīng)更為復(fù)雜的實際情況。

附錄見本刊網(wǎng)絡(luò)版（http：//www.epae.cn）。