燃?xì)夤芫W(wǎng)的混合整數(shù)規(guī)劃建模及其等價(jià)連續(xù)化

張 晴

(北京郵電大學(xué) 理學(xué)院，北京 100876)

隨著天然氣管道系統(tǒng)的規(guī)模越來(lái)越大，越來(lái)越復(fù)雜，導(dǎo)致這些設(shè)施的優(yōu)化運(yùn)行和規(guī)劃的重要性也越來(lái)越大.此外，管網(wǎng)的投資成本和運(yùn)行費(fèi)用非常大，即使是系統(tǒng)利用率的小幅提高，也會(huì)涉及大量資金.

天然氣行業(yè)的服務(wù)包括生產(chǎn)、輸送和銷售天然氣.本研究的主要興趣集中在通過(guò)管網(wǎng)輸送天然氣方面.在輸送天然氣的過(guò)程中，氣體流經(jīng)各種設(shè)備，包括管道及壓縮機(jī)等管網(wǎng)元素，氣體在管道中流動(dòng)時(shí)，由于與管壁的摩擦，導(dǎo)致氣體壓力降低，不能穩(wěn)定傳輸，因此需在管網(wǎng)的適當(dāng)位置建立壓氣站給氣體增壓，使之穩(wěn)定傳輸.天然氣管網(wǎng)系統(tǒng)的總運(yùn)行成本很大程度上取決于管網(wǎng)中壓氣站的運(yùn)行成本.而壓氣站的運(yùn)行成本，一般是以壓氣站消耗的燃料來(lái)衡量的.根據(jù)Luongo等人[1]的研究，壓氣站的運(yùn)行成本占公司總運(yùn)行預(yù)算的25%～50%.因此，傳輸網(wǎng)絡(luò)的目標(biāo)是在滿足規(guī)定的輸送流量和輸送終端最低壓力和溫度要求的前提下，盡量減少壓氣站的總?cè)剂舷腫2].

根據(jù)氣體流量隨時(shí)間的變化情況，管網(wǎng)系統(tǒng)可分為穩(wěn)態(tài)和瞬態(tài).當(dāng)系統(tǒng)中表征氣體流動(dòng)的數(shù)值與時(shí)間無(wú)關(guān)時(shí)，系統(tǒng)被稱為穩(wěn)態(tài).在此情況下，系統(tǒng)的約束條件，特別是描述氣體流經(jīng)管道的約束條件，可用代數(shù)非線性方程來(lái)描述.本文主要討論穩(wěn)態(tài)輸氣網(wǎng)絡(luò)優(yōu)化問(wèn)題.

最近20年的工作旨在解決所有網(wǎng)絡(luò)元素的非線性和離散方面.然而，由于目前的優(yōu)化算法仍然無(wú)法解決真實(shí)網(wǎng)絡(luò)的完整混合整數(shù)非線性規(guī)劃(MINLP)模型，因此主要采用線性近似的方法，以獲得混合整數(shù)線性規(guī)劃(MILP)模型(參見(jiàn)Hacklnder[3]，2002年；Martin和M?ller[4]，2005年；Martin等[5]，2006年；GeiBler 和Martin2015[6]，2015年；黃亞魁等[7]，2017年；Robert Burlacu等[8]，2019年).

1 管網(wǎng)模型

本節(jié)將介紹天然氣運(yùn)輸網(wǎng)絡(luò)優(yōu)化中的一個(gè)基本問(wèn)題.考慮一個(gè)天然氣運(yùn)輸網(wǎng)絡(luò)和預(yù)先已知的需要輸送的天然氣數(shù)量，即在每個(gè)入口，天然氣可以進(jìn)入網(wǎng)絡(luò)；而在每個(gè)出口，天然氣可以從網(wǎng)絡(luò)中輸出.此外，由于長(zhǎng)距離輸送管道上的摩擦引起的壓降，必須考慮可增加氣體壓力的壓縮機(jī).上述內(nèi)容導(dǎo)致了在找到滿足所有物理和技術(shù)約束的可控元件的可接受配置的同時(shí)，最小化壓縮機(jī)燃料消耗量的目標(biāo).由于氣體物理學(xué)是通過(guò)非線性函數(shù)來(lái)描述的，考慮到可控元素的離散特性，該問(wèn)題可被建模為MINLP.

1.1 穩(wěn)態(tài)管網(wǎng)的MINLP模型

為了模型的簡(jiǎn)化，本文討論的管網(wǎng)僅包含節(jié)點(diǎn)、管道和壓縮機(jī).

使用一個(gè)有向有限圖G=(V,A)來(lái)模擬氣體網(wǎng)絡(luò).節(jié)點(diǎn)集V由入口節(jié)點(diǎn)集Vs、出口節(jié)點(diǎn)集Vd和內(nèi)部節(jié)點(diǎn)集Vi組成.弧的集合A由管道的集合Ap和壓縮機(jī)的集合Ac組成.此外，弧a=(u,v)∈A上的質(zhì)量流量用qa表示，節(jié)點(diǎn)v∈V處的壓力和溫度分別用pv和Tv來(lái)表示.考慮到一個(gè)規(guī)范的網(wǎng)絡(luò)，對(duì)于每個(gè)節(jié)點(diǎn)v∈V，需給定壓力變量pv和溫度變量Tv的上下界約束：

(1)

(2)

此外，每個(gè)節(jié)點(diǎn)v∈V處需滿足質(zhì)量流量守恒：

(3)

其中：dv為v處的天然氣需求量，滿足dv≥0，若v∈Vd；dv≤0，若v∈Vs；dv=0，若v∈Vi.

1.1.1 管道

我國(guó)燃?xì)夤芫W(wǎng)中最常見(jiàn)的部件是管道.每一根管道a=(u,v)∈AP由其長(zhǎng)度La、直徑Da和粗糙度ka決定.管道內(nèi)的氣體流動(dòng)可以用質(zhì)量流量q、壓力p、溫度T、密度ρ和歐拉方程體系與具體狀態(tài)方程來(lái)描述.在本文中，狀態(tài)方程由真實(shí)氣體的熱力學(xué)標(biāo)準(zhǔn)方程給出：

ρRszT=p，

其中:Rs為通用氣體常數(shù)，z為壓縮系數(shù).后者可用下式所示：

上式中的假臨界溫度pc和假臨界溫度Tc均為常數(shù).通過(guò)一定的簡(jiǎn)化，管道上的壓降可以用著名的Weymouth方程[9]以解析形式給出：

(4)

(5)

此外，管道上的流量上下界由以下約束給出：

(6)

在考慮焦耳-湯姆遜效應(yīng)的影響下，管道上的溫降可以用著名的Sukhov公式加上焦耳-湯姆遜效應(yīng)的影響項(xiàng)[10]以解析形式給出：

?a=(u,v)∈Ap，

(7)

1.1.2 壓縮機(jī)

為了補(bǔ)償氣網(wǎng)中的壓降，管網(wǎng)中安裝了將入口壓力提高到較高出口壓力的壓縮機(jī).每臺(tái)壓縮機(jī)a=(u,v)∈Ac由其絕熱效率ηad,a和能耗成本系數(shù)ca指定.本文僅考慮渦輪壓縮機(jī)，其可行工作范圍由其特性圖來(lái)描述，見(jiàn)圖1中的灰色區(qū)域[11].

圖1 渦輪壓縮機(jī)的特性圖Figure 1 Characteristic diagram of turbo compressor

其中：Qa代表體積流量，與qa有如下關(guān)系：

(8)

其中：ρ0為標(biāo)況下的氣體密度.此外，壓縮機(jī)的絕熱過(guò)程會(huì)導(dǎo)致絕熱焓Had,a的變化：

(9)

其中：κ是絕熱指數(shù)，假定為常數(shù).此時(shí)壓縮系數(shù)za可近似為[10]：

(10)

壓縮機(jī)將進(jìn)口壓力pu提高到較高的出口壓力pv所消耗的功率Na為：

(11)

當(dāng)壓縮機(jī)關(guān)閉時(shí)，沒(méi)有氣體流經(jīng)壓縮機(jī)，其入口和出口的壓力和溫度是解耦的，顯然也不消耗功率.當(dāng)壓縮機(jī)運(yùn)行時(shí)，其在消耗功率的同時(shí)也增加了沿弧線方向流動(dòng)的氣體的壓力和溫度，即qa≥0，pv≥pu，Tv≥Tu，逆壓縮機(jī)弧線方向的流動(dòng)是不可能的.引入二元變量sa來(lái)表示壓縮機(jī)的不同狀態(tài)：sa=1表示壓縮機(jī)a是開(kāi)的，反之亦然.除了能耗約束(12)外，模型中還加入了以下約束條件：

(12)

(13)

(14)

(15)

puybmin≤pv≤puybmax，?a=(u,v)∈Ac，

(16)

sa∈{0,1}，a∈Ac.

(17)

1.2 MINLP問(wèn)題的總結(jié)

管網(wǎng)的運(yùn)行成本主要為壓縮機(jī)的燃料消耗，因此，本文假設(shè)每臺(tái)壓縮機(jī)a∈Ac的能耗成本系數(shù)ca≥0.考慮到前兩小節(jié)中描述的燃?xì)饩W(wǎng)絡(luò)中的不同元素，可得到下述MINLP問(wèn)題：

s.t.(1)～(17).

2 網(wǎng)絡(luò)約簡(jiǎn)、線性化及算法框架

本節(jié)將介紹一個(gè)用于構(gòu)建MINLP模型的MILP近似的方法及算法框架，以解決第二節(jié)中的問(wèn)題.本節(jié)討論的內(nèi)容是文獻(xiàn)[7]中的方法，為了內(nèi)容的完整性，在此節(jié)進(jìn)行簡(jiǎn)要回顧.

2.1 網(wǎng)絡(luò)約簡(jiǎn)

用虛擬壓縮機(jī)代替順序連接的壓縮機(jī)，該虛擬壓縮機(jī)的行為表現(xiàn)為其組件以最佳方式運(yùn)行[12].此處使用的約簡(jiǎn)方法為只替換順序相連的壓縮機(jī)站，如圖2轉(zhuǎn)變?yōu)閳D3所示.

圖2 循環(huán)網(wǎng)絡(luò)Figure 2 Cyclic network

圖3 約簡(jiǎn)網(wǎng)絡(luò)Figure 3 Reduced network

2.2 線性化-凸組合法

對(duì)原始網(wǎng)絡(luò)進(jìn)行約簡(jiǎn)后的簡(jiǎn)化網(wǎng)絡(luò)的規(guī)模雖然會(huì)減少，但由于目標(biāo)函數(shù)和某些約束的非凸非線性，導(dǎo)致直接調(diào)用MINLP求解器得到的解的質(zhì)量不會(huì)很高.其中一種處理該問(wèn)題的方法是基于建立混合整數(shù)線性模型的思想.用到的線性化方法是凸組合法[13]，其具體形式為：

λj≥0,j=1,…,m;zi∈{0,1},i=1,…,n.

其中S:={S1,…,Sn}為非線性函數(shù)y(x)的可行域剖分成的單純形的集合，m是S中的頂點(diǎn)個(gè)數(shù).此外，引入了輔助變量λ和z及其相應(yīng)的線性約束，其中z是0-1變量.

2.3 算法框架

求解穩(wěn)態(tài)燃?xì)饩W(wǎng)絡(luò)的基于線性化的算法如下所述：

算法3.1 優(yōu)化燃?xì)饩W(wǎng)絡(luò)的基于線性化的算法

步驟1 網(wǎng)絡(luò)約簡(jiǎn)

1)將網(wǎng)絡(luò)分解為僅包含順序壓縮器的小型子網(wǎng)站；

2)用虛擬壓縮器替換每個(gè)子網(wǎng)并獲得簡(jiǎn)化的網(wǎng)絡(luò)；

步驟2 建立管網(wǎng)的MINLP模型

步驟3 線性化

1)交替固定流量和其他變量，并選擇離散點(diǎn)；

2)優(yōu)化每個(gè)子網(wǎng)以獲得每個(gè)離散點(diǎn)的能耗和出口溫度；

3)使用可行網(wǎng)格點(diǎn)劃分可行區(qū)域；

4)建立MILP模型；

步驟4 解決MILP模型

步驟5 如果通過(guò)步驟4獲得的解滿足一定條件，則停止算法；若不滿足條件，轉(zhuǎn)到步驟3.

3 等價(jià)替換—連續(xù)化方法

模型中的整數(shù)變量的存在加大了問(wèn)題的求解難度，而已知的0-1規(guī)劃的求解方法可以分成離散方法和連續(xù)化方法兩類.

離散方法是直接解整數(shù)規(guī)劃的解法，連續(xù)化方法就是把問(wèn)題中的整數(shù)變量變?yōu)檫B續(xù)變量之后進(jìn)行求解，最簡(jiǎn)單的就是把0-1規(guī)劃中的0-1整數(shù)變量松弛到[0,1]區(qū)間上[14].本節(jié)的內(nèi)容就是通過(guò)在模型中添加一個(gè)關(guān)于0-1變量的非線性等式約束，從而把原問(wèn)題中的一部分0-1變量變?yōu)檫B續(xù)變量，并證明了進(jìn)行替換后的問(wèn)題與原問(wèn)題的等價(jià)性，與此同時(shí)利用罰函數(shù)將約束吸收到目標(biāo)函數(shù).該替換可以降低問(wèn)題求解難度，并提高問(wèn)題求解速度.

把1.1.2節(jié)中出現(xiàn)的變量sa通過(guò)如下所述的連續(xù)化方法變?yōu)閇0,1]區(qū)間上的連續(xù)變量，對(duì)?j均滿足下式：

證明：首先證明左式推右式：

由sa,j=0/1?sa,j(1-sa,j)=0?∑jαjsa,j(1-sa,j)=0，且滿足0≤sa,j≤1.

現(xiàn)在證明右式推左式：

∵0≤sa,j≤1，故0≤sa,j(1-sa,j)≤0.25，

又∵∑jαjsa,j(1-sa,j)=0?對(duì)?j，都有αj>0，故對(duì)?j,sa,j(1-sa,j)=0均成立，?sa,j=0/1，證畢.

因此，可以把模型中的整數(shù)變量連續(xù)化為[0,1]區(qū)間上的連續(xù)變量，進(jìn)而減小問(wèn)題的求解難度，并加快問(wèn)題求解速度.

4 數(shù)值結(jié)果

將前兩節(jié)中的算法應(yīng)用于我國(guó)西部天然氣管網(wǎng)系統(tǒng)，報(bào)告數(shù)值結(jié)果.

如圖4所示，該測(cè)試管網(wǎng)由27個(gè)壓縮機(jī)站，32條管道，64個(gè)節(jié)點(diǎn)和2個(gè)循環(huán)組成[7].輸氣量為1.4×106m3/d，壓縮機(jī)出入口的壓力范圍均是[8,12](MPa)，溫度取值范圍為[0,50](℃).在圖4中，豎著的黑色梯形塊符號(hào)代表的是壓縮機(jī)，箭頭表示的是氣體的流動(dòng)方向.圖5是相應(yīng)的約簡(jiǎn)管網(wǎng).

圖4 測(cè)試管網(wǎng)Figure 4 Test network

圖5 約簡(jiǎn)管網(wǎng)Figure 5 Reduced network

建立約簡(jiǎn)管網(wǎng)的MINLP模型，采用算法3.1求解該模型，并將此結(jié)果與使用等價(jià)連續(xù)化的結(jié)果進(jìn)行比較.為了表達(dá)的便利，算法3.1稱為方法1，對(duì)整數(shù)變量進(jìn)行等價(jià)連續(xù)化的算法3.1稱為方法2.

本節(jié)中介紹的所有計(jì)算都是在Intel Core i7，1.80 GHz的處理器上完成的.在每一次迭代中，均需選擇離散點(diǎn)，并對(duì)每一子網(wǎng)進(jìn)行優(yōu)化.測(cè)試中使用DP優(yōu)化子網(wǎng)，使用建模語(yǔ)言YALMIP[15]為約簡(jiǎn)管網(wǎng)構(gòu)建MILP模型.作為MILP求解器，本文使用CPLEX12.8.0[16]，其中所有參數(shù)都設(shè)置為默認(rèn)值.

可看到可行的解決方案是通過(guò)三次迭代獲得的.如前所述，簡(jiǎn)化網(wǎng)絡(luò)僅包含虛擬壓縮機(jī).因此，需測(cè)量耗氣量與出口溫度的準(zhǔn)確性.當(dāng)誤差小于給定的公差ε時(shí)，則接受該解.在實(shí)驗(yàn)中，設(shè)置誤差精度ε=5%，超過(guò)兩次迭代后，精度取其一半.表1中給出兩種方法得到的解.

表1 MILP約簡(jiǎn)網(wǎng)絡(luò)解的質(zhì)量比較Table 1 Quality comparison of MILP approximate network solutions

從表1中可看出，第一次迭代時(shí)，除了VCS5子網(wǎng)的能耗誤差稍大于給定誤差外，其余子網(wǎng)的誤差均達(dá)到給定精度，再根據(jù)表中的運(yùn)行時(shí)間，表明兩種方法均可在較短時(shí)間內(nèi)得到可行解，這對(duì)于實(shí)際管網(wǎng)優(yōu)化問(wèn)題來(lái)說(shuō)十分重要，而且方法2的用時(shí)比方法1的用時(shí)短.第三次迭代后，耗氣量和出口溫度均滿足給定的誤差精度，說(shuō)明兩種方法都可以提供令人滿意的解決方案.此外，方法2獲得的耗氣量較低，為22 711.94 m3/d，比方法1的耗氣量降低了133.85 m3/d.兩者均提供了一種子網(wǎng)VCS2和VCS3關(guān)閉的情況下的運(yùn)行方案.

跟蹤約簡(jiǎn)管網(wǎng)的解可得到原始網(wǎng)絡(luò)的解，表2顯示出相應(yīng)結(jié)果，其中粗體說(shuō)明對(duì)應(yīng)的壓縮機(jī)是關(guān)閉的，可看到在滿足管網(wǎng)各物理和技術(shù)約束的條件下，兩種方法均關(guān)閉了15個(gè)壓縮機(jī)，但每個(gè)壓縮機(jī)的出入口壓力、溫度不盡相同.換句話說(shuō)，通過(guò)對(duì)管道中的氣體進(jìn)行適當(dāng)?shù)脑黾踊驕p少氣體壓力或者溫度可達(dá)到降低管網(wǎng)系統(tǒng)耗氣量的效果.由此可見(jiàn)，管網(wǎng)中各變量值的微小改變可帶來(lái)非?？捎^的經(jīng)濟(jì)效益.

表2 MILP原始網(wǎng)絡(luò)解的比較Table 2 Comparison of MILP original network solutions

綜上所述，上節(jié)提出的方法可有效減少管網(wǎng)的耗氣量，運(yùn)行時(shí)間也有所減少.

5 結(jié) 語(yǔ)

本文提出了針對(duì)模型中整數(shù)變量的等價(jià)替換形式，可進(jìn)一步減少問(wèn)題求解時(shí)間，甚至可影響CPLEX的分支過(guò)程得出更優(yōu)的解.數(shù)值結(jié)果也證明了本文方法的有效性.下一步希望能找到更多的燃?xì)夤芫W(wǎng)實(shí)例來(lái)研究該方法有效性，同時(shí)進(jìn)一步提高算法的求解效率.然而，實(shí)際生產(chǎn)中的管網(wǎng)系統(tǒng)的變量隨時(shí)間是不斷變化的(即關(guān)于時(shí)間的函數(shù))，所以對(duì)非穩(wěn)態(tài)的天然氣管網(wǎng)運(yùn)行優(yōu)化的研究值得進(jìn)一步探索.