混合建模方法在蒸汽供熱管網(wǎng)建模中的應(yīng)用

胡佳怡，鄭夢蓮

(浙江大學(xué) 熱工與動力系統(tǒng)研究所，浙江 杭州 310027)

0 引 言

隨著我國能源轉(zhuǎn)型進程的進行，能源利用轉(zhuǎn)向追求更加高效、清潔和智能化的用能。區(qū)域供熱系統(tǒng)是重要的能源消費者之一，目前大多數(shù)區(qū)域供熱系統(tǒng)仍采用燃燒化石燃料供暖，化石燃料的過多燃燒會增加碳排放從而影響環(huán)境，因此，提高我國區(qū)域供熱系統(tǒng)的能源利用效率十分重要，對供熱系統(tǒng)的準確建模有利于準確估計系統(tǒng)的能源使用情況，從而進行運行分析和調(diào)度控制[1-2]。

對供熱系統(tǒng)進行建模的方法包括水力計算、熱力計算和水力熱力耦合計算。水力計算通常采用圖論法[3]，建立水力瞬態(tài)模型可以為管網(wǎng)的有效運行提供基礎(chǔ)[4]。而阻力系數(shù)是水力建模過程中的一個重要參數(shù)，對其進行辨識有助于提高建模準確性[5]。熱力計算則通過熱力公式和模型對管道進行數(shù)值求解從而求得溫度分布[6-7]。在實際供熱管網(wǎng)中，水力和熱力參數(shù)相互影響，因此還需要進行水力熱力耦合建模[2,8]。

由于管道結(jié)構(gòu)及材料老化等影響，采用上述機理方法仍可能存在建模誤差。數(shù)據(jù)驅(qū)動方法是一種黑箱模型，根據(jù)歷史數(shù)據(jù)進行建模，避免了復(fù)雜的機理建模，已被用于評估供熱管網(wǎng)的熱損失和熱需求[9]。但是數(shù)據(jù)驅(qū)動方法由于物理模型缺失會導(dǎo)致可解釋性弱和外推性差等問題，因此本文提出通過數(shù)據(jù)—物理混合建模來對供熱管網(wǎng)進行建模。混合建模方法在建筑節(jié)能[10]、電池健康度預(yù)測[11]和電力系統(tǒng)[12]等方面都有應(yīng)用，根據(jù)其融合方法的不同可以分為以下三種類型：串聯(lián)模型將一種模型的結(jié)果作為另一種模型的輸入[13]；并聯(lián)模型將數(shù)據(jù)驅(qū)動模型和機理模型的結(jié)果通過加權(quán)等方式進行疊加得到最后結(jié)果[14]；引導(dǎo)模式以機理方法為主，通過數(shù)據(jù)驅(qū)動方法預(yù)測其中的某些未知參數(shù)或環(huán)節(jié)[15]。現(xiàn)有的文獻更多地討論了混合模型與數(shù)據(jù)驅(qū)動模型或物理模型之間的差異[16]，而對不同混合建模方法的比較較少。因此本文著重探討不同混合模型的建模效果比較。

目前工業(yè)園區(qū)的蒸汽供熱管網(wǎng)建模主要以機理建模為主，本文考慮了機理建模與數(shù)據(jù)驅(qū)動混合的建模方法進行建模，從而提高模型精度，并比較了三種混合建模方法的區(qū)別。

1 方法

1.1 水力熱力模型

水力模型基于質(zhì)量守恒定律和能量守恒定律，離散的管段模型不利于計算，圖論方法可以將其抽象為矩陣從而方便水力計算。以下為圖論計算方法，首先定義節(jié)點—管段關(guān)聯(lián)矩陣A 和基環(huán)-管段關(guān)聯(lián)矩陣B：

其中，m為管段數(shù)；n為節(jié)點數(shù)；k為基環(huán)數(shù)。

質(zhì)量守恒定律可以表示為：

其中，G為流量向量，kg/s；g為凈流量，無泄漏和疏水情況下為0，kg/s。

環(huán)能量方程可以表示為：

其中，ΔP 為各管段的壓降，MPa。

蒸汽供熱管道的壓降計算公式如下：

其中，λ是沿程阻力系數(shù)；L為管段長度，m；d0為管道內(nèi)徑，m；ρ為工質(zhì)密度，kg/m3；v為管內(nèi)蒸汽流速，m/s；α為局部阻力與沿程阻力的比值。

管段摩擦阻力公式可以表示為：

其中，K為管壁絕對粗糙度，m。

管道的阻力系數(shù)定義如式（11）所示，已知阻力系數(shù)及流量則可以方便地計算管段壓降。

其中，S為管段阻力系數(shù)。

熱力模型通過計算內(nèi)外溫差及保溫?zé)嶙鑱碛嬎愎芏紊崃浚瑢τ陔p層保溫管道，其單位管長散熱量可以表示為：

其中，q為單位管長的散熱量，W/m；η為管道散熱修正系數(shù)，與管道結(jié)構(gòu)、保溫老化和環(huán)境熱阻等因素相關(guān)，tavg為管內(nèi)流體平均溫度，℃；ta為環(huán)境溫度，℃；λ1為內(nèi)保溫層導(dǎo)熱系數(shù)，W/(m·℃)；λ2為外保溫層導(dǎo)熱系數(shù)，W/(m·℃)；αo為環(huán)境散熱系數(shù)，W/(m2·℃)；d1為管道外徑，m；d2為內(nèi)保溫層外徑，m；d3為外保溫層外徑，m。

通過下式可以計算蒸汽管道的沿途溫降：

其中，h1為管道起點焓值，kJ/kg；h2為管道終點焓值，kJ/kg。

對于實際供熱管道，水力參數(shù)與熱力參數(shù)存在著相互影響的關(guān)系，單一的水力模型或熱力模型可能導(dǎo)致結(jié)果不準確，因此，本研究經(jīng)過熱力水力參數(shù)迭代后得到最終的管網(wǎng)機理模型。

1.2 混合建模方法

本文主要采用串聯(lián)式、并聯(lián)式和引導(dǎo)式三種混合建模方法對管段的壓力和溫度進行預(yù)測，其建模流程如圖1 所示。并聯(lián)建模方法先采用機理方法進行水力熱力計算，然后采用數(shù)據(jù)驅(qū)動方法對實際值與機理值之間的誤差E 進行預(yù)測，最后將機理結(jié)果和誤差E 疊加得到最后的預(yù)測結(jié)果。串聯(lián)方法將機理方法得到的計算結(jié)果作為數(shù)據(jù)驅(qū)動方法的輸入，再通過數(shù)據(jù)驅(qū)動方法對壓力和溫度進行預(yù)測。引導(dǎo)方法先通過數(shù)據(jù)驅(qū)動的方式，對水力熱力模型中的阻力系數(shù)和散熱修正系數(shù)進行辨識，再將辨識結(jié)果代入機理模型中進行計算。

圖1 混合建模方法示意圖

1.3 最小二乘辨識

最小二乘方法具有簡單和運算速度快的特點，傳統(tǒng)的非線性最小二乘求解方法高斯牛頓法容易出現(xiàn)雅可比矩陣的問題，而Levenberg-Marquard（LM）方法則可以避免相關(guān)問題的發(fā)生[17]，因此，本文采用LM 方法對阻力系數(shù)和散熱修正系數(shù)進行辨識。

辨識阻力系數(shù)的目的是使得辨識后的壓力觀測值與預(yù)測值之間的誤差最小，其目標(biāo)函數(shù)可以表示為：

其中，pp是壓力預(yù)測值，MPa；pr是壓力觀測值，MPa；是權(quán)重系數(shù)，這里取1；M是觀測點數(shù)量；NLC是總數(shù)據(jù)量。

考慮到阻抗的實際值與計算值的偏差應(yīng)該在一定范圍內(nèi)，設(shè)定阻抗的約束條件可以表示為：

其中，Si是第i條管線的阻抗值；ci和di是第i條管線的阻力系數(shù)上下界搜索系數(shù)。

同理，辨識散熱修正系數(shù)時目標(biāo)函數(shù)可以表示為：

其中，tp是溫度預(yù)測值，℃；tr是溫度觀測值，℃。散熱修正系數(shù)的約束條件可以表示為：

其中，ηi是第i條管線的阻抗值；s和t是管線的散熱修正系數(shù)上下界。

1.4 多輸出高斯過程回歸

供熱系統(tǒng)建模過程中需對熱用戶側(cè)壓力溫度進行預(yù)測，各個熱用戶處的溫度和壓力可能存在潛在的相關(guān)性。高斯過程回歸是一種基于高斯過程先驗知識得到預(yù)測結(jié)果的方法，單輸出高斯過程回歸輸出為一個值，多輸出高斯過程回歸輸出為一組向量，其核函數(shù)可以捕獲輸出之間相關(guān)性，因此，與多輸出高斯過程回歸相比，單輸出高斯過程回歸在解決多輸出問題時可以得到更好的結(jié)果[18]，本文采用多元高斯過程回歸方法建模。

2 算例

2.1 數(shù)據(jù)集

本文選取某工業(yè)園區(qū)的一部分支狀供熱管網(wǎng)作為研究對象，采集了該管網(wǎng)2020 年12 月某三天的流量、壓力和溫度數(shù)據(jù)，共4175 個數(shù)據(jù)點作為數(shù)據(jù)集。管網(wǎng)簡化圖如圖2 所示，包括14 個溫度壓力測點u1-u14。計算所采用的管道和保溫參數(shù)如表1 所示，參考城鎮(zhèn)供熱管網(wǎng)設(shè)計規(guī)范選擇局部與沿程阻力比值[19]。

表1 管道設(shè)計參數(shù)

圖2 支狀供熱管網(wǎng)示意圖

2.2 水力熱力模型參數(shù)辨識

計算阻力系數(shù)設(shè)計值時，考慮造成阻力系數(shù)辨識值與設(shè)計值偏差的原因是絕對當(dāng)量粗糙度和局部沿程阻力比值的取值偏差，因此假設(shè)了絕對當(dāng)量粗糙度和局部沿程阻力比值的上下限，作為辨識約束條件的系數(shù)ci和系數(shù)di，選取合適的上下限范圍有助于辨識到更接近真實情況的阻力系數(shù)值，具體取值如表2 所示。

表2 管道當(dāng)量粗糙度及沿程局部阻力比值

阻力系數(shù)的設(shè)計值和辨識值如表3 所示，部分阻力系數(shù)辨識值與設(shè)計值相差較大，可能是由于壓力測點的儀表測量誤差或未考慮到的局部阻力系數(shù)導(dǎo)致的。

表3 阻力系數(shù)計算值

設(shè)置散熱修正系數(shù)的上下邊界s和t為1 和15。散熱修正系數(shù)辨識值如表4 所示，部分管道散熱修正系數(shù)偏大，原因是管道可能存在保溫老化或局部破壞等問題，散熱修正系數(shù)較小的管道保溫情況相對良好。

表4 散熱修正系數(shù)辨識值

2.3 混合建模結(jié)果

混合建模方法中，并聯(lián)和串聯(lián)建模方法采用多輸出高斯過程回歸對壓力和溫度進行預(yù)測，由于核函數(shù)的選取會影響建模精度，本文采用組合核函數(shù)來提高建模精度，選取的核函數(shù)由徑向基核函數(shù)、Matern32 核函數(shù)和線性核函數(shù)組成。此外，數(shù)據(jù)集中訓(xùn)練集和測試集的比例為85%和15%，訓(xùn)練集用于模型訓(xùn)練，測試集用于最終的模型評估，模型評價指標(biāo)包括均方根誤差（Root Mean Square Error, RMSE）、平均絕對誤差（Mean Absolute Error, MAE）和平均絕對百分比誤差（Mean Absolute Percentage Error, MAPE）。圖3和圖4 分別是不同建模方法下溫度和壓力的預(yù)測結(jié)果，可以看到不同方法下預(yù)測結(jié)果與觀測結(jié)果的變化趨勢基本相同。幾種建模方法中機理計算誤差最大。圖中大多數(shù)測點的溫度機理計算值大于觀測值，是因為機理計算低估了管道由于老化等原因造成的散熱增加量，而大多數(shù)測點的壓力機理計算值小于觀測值，是由計算得到的阻力系數(shù)值偏大引起的。采用混合建模方法預(yù)測溫度時，引導(dǎo)建模方法和其他兩種方法預(yù)測值較為接近，采用混合建模方法預(yù)測壓力時，引導(dǎo)建模方法結(jié)果與機理計算值更為接近，是因為引導(dǎo)方法以機理計算模型為計算基礎(chǔ)，使得采用引導(dǎo)方法預(yù)測壓力時更接近機理值，因此設(shè)定的機理模型的準確度會影響引導(dǎo)建模結(jié)果。

圖3 不同建模方法下的溫度預(yù)測結(jié)果比較

為了在統(tǒng)計意義上獲得更準確的結(jié)果，采用隨機劃分的方法對測試集和訓(xùn)練集進行30 次不同的劃分，溫度和壓力的預(yù)測結(jié)果誤差如圖5 和圖6 所示。其中串聯(lián)方法的預(yù)測誤差較低，溫度和壓力測點的平均RMSE 分別為0.36 和0.0055，引導(dǎo)方法的預(yù)測誤差相對較高，溫度和壓力測點的平均RMSE 為1.47 和0.043。由于串聯(lián)模型結(jié)構(gòu)的最后一層是數(shù)據(jù)驅(qū)動模型，本文數(shù)據(jù)集質(zhì)量較好且機理建模結(jié)果可以較好地輔助數(shù)據(jù)驅(qū)動模型預(yù)測，因此得到了較好的預(yù)測結(jié)果。并聯(lián)模型是機理結(jié)果與數(shù)據(jù)驅(qū)動結(jié)果的疊加，其精度取決于數(shù)據(jù)集是否能準確預(yù)測機理結(jié)果與實際值之間的誤差。引導(dǎo)模型預(yù)測結(jié)果受到機理模型影響，因此機理模型誤差導(dǎo)致了引導(dǎo)模型預(yù)測精度低，但是往往引導(dǎo)模型更能反映真實物理趨勢且具有一定的外推性。

圖5 混合建模方法的溫度預(yù)測誤差

圖6 混合建模方法的壓力預(yù)測誤差

3 結(jié)論

本文采用不同建模方法對工業(yè)園區(qū)的供熱管網(wǎng)進行熱力水力計算，由于機理—數(shù)據(jù)混合建模可以避免單一機理模型和數(shù)據(jù)驅(qū)動模型可能存在的問題，本文提出了三種不同的混合建模方法，比較了不同建模方法之間的模型結(jié)構(gòu)差異和建模結(jié)果差異。結(jié)果顯示提出的三種數(shù)據(jù)—機理混合建模方法相比單一的機理計算方法可以得到更好的計算結(jié)果。采用本文數(shù)據(jù)集進行建模時，串聯(lián)和并聯(lián)建模方法預(yù)測精度較高，引導(dǎo)建模方法預(yù)測精度較低，由于串聯(lián)模型和并聯(lián)模型更加依賴數(shù)據(jù)驅(qū)動模型，因此在數(shù)據(jù)集質(zhì)量較高的情況下可以采用這兩種建模方法，而并聯(lián)模型相比串聯(lián)模型需要考慮目前數(shù)據(jù)集中的特征是否可以較好地用于預(yù)測誤差，對于引導(dǎo)建模則需要考慮機理模型的準確性。