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

胡佳怡，鄭夢(mèng)蓮

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

0 引 言

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

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

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

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

1 方法

1.1 水力熱力模型

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

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

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

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

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

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

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

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

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

其中，K為管壁絕對(duì)粗糙度，m。

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

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

熱力模型通過(guò)計(jì)算內(nèi)外溫差及保溫?zé)嶙鑱?lái)計(jì)算管段散熱量，對(duì)于雙層保溫管道，其單位管長(zhǎng)散熱量可以表示為：

其中，q為單位管長(zhǎng)的散熱量，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。

通過(guò)下式可以計(jì)算蒸汽管道的沿途溫降：

其中，h1為管道起點(diǎn)焓值，kJ/kg；h2為管道終點(diǎn)焓值，kJ/kg。

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

1.2 混合建模方法

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

圖1 混合建模方法示意圖

1.3 最小二乘辨識(shí)

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

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

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

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

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

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

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

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

1.4 多輸出高斯過(guò)程回歸

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

2 算例

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

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

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

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

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

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

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

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

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

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

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

2.3 混合建模結(jié)果

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

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

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

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

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

3 結(jié)論

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