技術(shù)特征：

1.一種天然氣系統(tǒng)能流計(jì)算的改進(jìn)方法，其特征在于，包括以下步驟:

(1)輸入基礎(chǔ)數(shù)據(jù)及初始化

1.1輸入基礎(chǔ)數(shù)據(jù)

獲取天然氣系統(tǒng)網(wǎng)絡(luò)結(jié)構(gòu)數(shù)據(jù)：氣源參數(shù)、輸氣管道參數(shù)、壓縮機(jī)參數(shù)、氣負(fù)荷、平衡節(jié)點(diǎn)壓力以及待求節(jié)點(diǎn)壓力范圍。

1.2參數(shù)初始化

設(shè)置遺傳算法的交叉概率P_C，變異概率P_m，變量精度esp，種群規(guī)模Pop，表示種群中全部個(gè)體所在集合；牛頓法收斂精度ε，牛頓法的最大迭代次數(shù)iterNWmax。

設(shè)置天然氣系統(tǒng)中除平衡節(jié)點(diǎn)外的節(jié)點(diǎn)個(gè)數(shù)為N；表示網(wǎng)絡(luò)N個(gè)節(jié)點(diǎn)所在的集合。

(2)基于遺傳算法的初值優(yōu)選

2.1編碼并產(chǎn)生初始種群

第(1)步完成之后，以一個(gè)非平衡節(jié)點(diǎn)上的節(jié)點(diǎn)壓力為一個(gè)變量，則一共有N個(gè)變量，根據(jù)變量范圍和變量精度，按照二進(jìn)制編碼方法，確定第s(s＝1,2,···,N)個(gè)變量對應(yīng)的二進(jìn)制串長度為n_s，其數(shù)值滿足下列關(guān)系式：

$2^{n_s-1}\≤\frac{{\mathrm{\Π}}_{\max ,s}-{\mathrm{\Π}}_{\min ,s}}{esp}\≤2^{n_s},(s=1,2,...,N)---\left(1\right)$

式中esp的值在步驟1.2中給出；Π_max,s、Π_min,s分別表示第s個(gè)節(jié)點(diǎn)壓力變量的上限和下限；

確定每個(gè)變量對應(yīng)的二進(jìn)制串長度后，每個(gè)個(gè)體包含N個(gè)變量的二進(jìn)制串，則確定每個(gè)個(gè)體的染色體長度為

$L=\underset{s=1}{\overset{N}{\Σ}}{n}_s---\left(2\right)$

運(yùn)用Matlab產(chǎn)生基于狀態(tài)變量的初始種群TP，種群規(guī)模均為Pop。具體步驟為：運(yùn)用Matlab命令rand(Pop,L)產(chǎn)生一個(gè)Pop×L的矩陣T，其元素值為區(qū)間(0,1)上的隨機(jī)數(shù)；然后判斷矩陣T中的每個(gè)元素，若值小于0.5，則該元素置為0，否則置為1；從而得到二進(jìn)制編碼的矩陣TP。其中，TP的每一行表示攜帶N個(gè)節(jié)點(diǎn)壓力變量信息的一個(gè)個(gè)體，是位數(shù)為L的二進(jìn)制碼。

2.2譯碼并得變量的實(shí)際值

將各個(gè)二進(jìn)制編碼轉(zhuǎn)換為對應(yīng)的變量的實(shí)際值Π_{m_s}，轉(zhuǎn)換公式如下：

${\mathrm{\Π}}_{m\_s}=\left({\mathrm{\Π}}_{\min ,s}+{\left({\mathrm{TP}}_{m\_s}\right)}_D\frac{{\mathrm{\Π}}_{\max ,s}-{\mathrm{\Π}}_{\min ,s}}{2^{n_s}-1}\right),m\∈\left\{1,2,\mathrm{...},200\right\}---\left(3\right)$

式中，TP_{m_s}表示第m個(gè)個(gè)體上第s個(gè)節(jié)點(diǎn)壓力變量的二進(jìn)制編碼，即二進(jìn)制編碼矩陣TP的第m行第至第列二進(jìn)制編碼，Π_{m_s}表示TP_{m_s}代表的變量實(shí)際值，(TP_{m_s})_D表示TP_{m_s}對應(yīng)的變量的十進(jìn)制數(shù)值，，求取公式為：

${\left({\mathrm{TP}}_{m\_s}\right)}_D=A_{n_s}\·2^{n_s-1}+A_{n_s-1}\·2^{n_s-2}+...+A_1\·2^0---\left(4\right)$

其中A_l表示TP_{m_s}的第l位二進(jìn)制數(shù)；其余變量的含義同公式(1)-(2)。

2.3計(jì)算個(gè)體對應(yīng)的目標(biāo)函數(shù)值

2.3.1.計(jì)算管道流量

根據(jù)步驟個(gè)體中各個(gè)節(jié)點(diǎn)的壓力變量實(shí)際值，得管道i-j流量f_p,ij的計(jì)算公式如下：

$f_{p,ij}=k_{ij}{s}_{p,ij}\sqrt{s_{p,ij}({{\mathrm{\Π}}_i}^2-{{\mathrm{\Π}}_j}^2)}---\left(5\right)$

式中，下標(biāo)p表示管道，i表示系統(tǒng)中的某一節(jié)點(diǎn)i，j表示與節(jié)點(diǎn)i直接相連的另一個(gè)節(jié)點(diǎn)，j∈φ_i，φ_i為與節(jié)點(diǎn)i直接相聯(lián)但不包括節(jié)點(diǎn)i的節(jié)點(diǎn)集合；k_ij為管道的傳輸參數(shù),Π_i和Π_j分別表示節(jié)點(diǎn)i和節(jié)點(diǎn)j的節(jié)點(diǎn)壓力；s_p,ij表示的是節(jié)點(diǎn)i，j間管道中的氣流方向，其取值如下所示：

$s_{p,ij}=\left\{\begin{array}{cc}+1,& {\mathrm{\Π}}_i-{\mathrm{\Π}}_j\≥0\\ -1,& {\mathrm{\Π}}_i-{\mathrm{\Π}}_j\<0\end{array}\right.---\left(6\right)$

2.3.2計(jì)算燃?xì)廨啓C(jī)消耗的氣流量

由管道方程可知，天然氣在傳輸過程中會(huì)有壓力損失，一般需要在天然氣系統(tǒng)中配置壓縮機(jī)來提升壓力。選擇燃?xì)廨啓C(jī)作為壓縮機(jī)，選擇作用于節(jié)點(diǎn)i和節(jié)點(diǎn)j間的壓縮機(jī)，記為壓縮機(jī)i-j，則其消耗的氣流量公式如下：

τ_c,ij＝α_c+β_cH_c,ij+γ_cH_c,ij² (7)

式中，下標(biāo)c代表壓縮機(jī)編號(hào)；α_c、β_c和γ_c為能量轉(zhuǎn)化效率常數(shù)；

H_c,ij為壓縮機(jī)消耗的功率，其求取公式如下：

$H_{c,ij}=B_{c,ij}{f}_{c,ij}\[{\left(\frac{{\mathrm{\Π}}_j}{{\mathrm{\Π}}_i}\right)}^{Z_{c,ij}}-1\]---\left(8\right)$

式中，B_c,ij和Z_c,ij為壓縮機(jī)參數(shù)；f_c,ij為流過壓縮機(jī)i-j的氣流量，其余變量的含義同公式(5)。

2.3.3計(jì)算天然氣系統(tǒng)的節(jié)點(diǎn)能流平衡方程

天然氣系統(tǒng)中，節(jié)點(diǎn)i的能流平衡方程如下所示：

$f_i=\underset{j\∈{\mathrm{\φ}}_i}{\mathrm{\Σ}}{f}_{p,ij}+\underset{j\∈{\mathrm{\φ}}_i}{\mathrm{\Σ}}{s}_{c,ij}{f}_{c,ij}+\underset{j\∈{\mathrm{\φ}}_i}{\mathrm{\Σ}}{\mathrm{\τ}}_{c,ij}---\left(9\right)$

式中，f_i為節(jié)點(diǎn)i的注入氣流量，s_c,ij是用來表示壓縮機(jī)的方向的符號(hào)變量，其公式如下：

式(9)其余變量含義與式(5)到式(8)一樣。

2.3.4.計(jì)算目標(biāo)函數(shù)值

將天然氣系統(tǒng)中能流方程中各節(jié)點(diǎn)不平衡量的總和作為目標(biāo)函數(shù)值。假設(shè)對于天然氣系統(tǒng)，其能流方程的一般形式為：

$\left\{\begin{array}{c}{f}_1\left({\mathrm{\Π}}_1,{\mathrm{\Π}}_2,\mathrm{...},{\mathrm{\Π}}_N\right)=0\\ f_2\left({\mathrm{\Π}}_1,{\mathrm{\Π}}_2,\mathrm{...},{\mathrm{\Π}}_N\right)=0\\ \begin{array}{c}.\\ .\\ .\end{array}\\ f_N\left({\mathrm{\Π}}_1,{\mathrm{\Π}}_2,\mathrm{...},{\mathrm{\Π}}_N\right)=0\end{array}\right.---\left(11\right)$

則代入每個(gè)節(jié)點(diǎn)的壓力實(shí)際值至式(9)，得出每個(gè)節(jié)點(diǎn)對應(yīng)的注入氣流量的不平衡量如下：

Δf_s(Π)＝f_s-f_s(Π₁,Π₂,…,Π_N) (12)

算出每個(gè)節(jié)點(diǎn)的注入其流量的不平衡量Δf_s(Π)，得出能流方程中不平衡量的總和為：

$\left\{\begin{array}{c}F\left(\mathrm{\Π}\right)=\underset{s=1}{\overset{N}{\mathrm{\Σ}}}\left|{\mathrm{\Δf}}_s\left(\mathrm{\Π}\right)\right|\\ \mathrm{\Π}={\[{\mathrm{\Π}}_1,{\mathrm{\Π}}_2,\mathrm{...},{\mathrm{\Π}}_N\]}^T\end{array}\right.---\left(13\right)$

式中，F(xiàn)(Π)為目標(biāo)函數(shù)，即能流方程中不平衡量的總量。

2.4收斂判據(jù)

計(jì)算出每一個(gè)個(gè)體對應(yīng)的目標(biāo)函數(shù)F_m(Π)，比較目標(biāo)函數(shù)的大小，選出其中的最小值min F(Π)。遺傳算法的收斂條件是，若min F(Π)小于閾值，則滿足收斂條件，將min F(Π)對應(yīng)的個(gè)體中的每個(gè)變量作為初值，代入步驟(3)；若不滿足收斂條件，即繼續(xù)執(zhí)行步驟2.5，直至滿足收斂條件為止。

2.5產(chǎn)生新個(gè)體

2.5.1選擇

步驟2.3計(jì)算完初始種群TP的全部個(gè)體對應(yīng)的目標(biāo)函數(shù)后，需要用適應(yīng)度大小評判個(gè)體優(yōu)劣，從而決定其遺傳機(jī)會(huì)的多少。以不平衡量總和F(Π)最小作為優(yōu)化目標(biāo)時(shí)，個(gè)體m的適應(yīng)度函數(shù)如下：

$f_m=\frac{1}{F_m\left(\mathrm{\Π}\right)}---\left(14\right)$

式中：f_m為表示第m個(gè)個(gè)體對應(yīng)的適應(yīng)度值，F(xiàn)_m(Π)為第m個(gè)個(gè)體對應(yīng)的目標(biāo)函數(shù)值。

采用輪盤賭選取，每個(gè)個(gè)體被選中的概率為

$P\left(f_m\right)=\frac{f_m}{\underset{m=1}{\overset{Pop}{\Σ}}{f}_m}---\left(15\right)$

具體操作為：在區(qū)間[0,1]產(chǎn)生一個(gè)均勻分布的隨機(jī)數(shù)r，若滿足則表示個(gè)體w被選中。重復(fù)選擇Pop次，得到Pop個(gè)經(jīng)過選擇后的個(gè)體，組成父代種群T^f。

2.5.2交叉

對父代種群T^f進(jìn)行交叉操作。具體步驟為：將種群T^f中的個(gè)體和(v＝1，2，…，100)作為父代，對父代第1個(gè)變量對應(yīng)的二進(jìn)制碼，運(yùn)用Matlab命令rand()隨機(jī)產(chǎn)生一個(gè)(0,1)上的數(shù)r,若r<Pc(Pc為交叉概率，已在步驟1.2中給出)，n₁為第一個(gè)變量的二進(jìn)制串長度，取rn₁的整數(shù)部分c₁作為第1個(gè)變量的二進(jìn)制碼交叉點(diǎn)位置，將父代個(gè)體和中第1個(gè)變量的二進(jìn)制碼的第c₁至n₁位進(jìn)行互換，按照同樣方法，對父代和中剩下的N-1個(gè)變量進(jìn)行二進(jìn)制碼交叉操作，得到子代個(gè)體和從而可以得到新的子代種群T^son，并且該種群規(guī)模保持Pop個(gè)個(gè)體不變。

2.5.3變異

對子代種群T^son進(jìn)行變異操作，具體步驟為：對種群T^son中的個(gè)體(m＝1，2，…，Pop)，對個(gè)體上第1個(gè)變量對應(yīng)的二進(jìn)制碼，運(yùn)用Matlab命令rand()隨機(jī)產(chǎn)生一個(gè)(0,1)上的數(shù)r，若r<Pm(Pm為變異概率，已在步驟1.2中給出)，取rn₁的整數(shù)部分d₁作為第1個(gè)變量二進(jìn)制碼的變異點(diǎn)位置，則該個(gè)體第1個(gè)變量二進(jìn)制碼的第d₁位二進(jìn)制數(shù)發(fā)生變異，即原來為1則變?yōu)?，原來為0則變?yōu)?，若r>Pm，則該變量對應(yīng)的二進(jìn)制碼不變，對該個(gè)體上剩下的N-1個(gè)變量對應(yīng)的二進(jìn)制碼作同樣的操作，得到新個(gè)體從而得到新種群TP＝T^New，規(guī)模仍為Pop。轉(zhuǎn)到步驟2.2繼續(xù)計(jì)算。

(3)基于牛頓法的天然氣系統(tǒng)能流計(jì)算

3.1計(jì)算不平衡量

執(zhí)行步驟(1)和步驟(2)，可以獲得最優(yōu)的個(gè)體中代表的各節(jié)點(diǎn)壓力，即較好的牛頓法初值。結(jié)合式(5)至式(10)，計(jì)算出個(gè)體中所有節(jié)點(diǎn)的不平衡量Δf^(k)，即

${\mathrm{\&Delta;f}}^{\left(k\right)}=\left[\begin{array}{c}{\mathrm{\&Delta;f}}_1^{\left(k\right)}\\ {\mathrm{\&Delta;f}}_2^{\left(k\right)}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{\&Delta;f}}_N^{\left(k\right)}\end{array}\right]---\left(16\right)$

式中，上標(biāo)k表示的是迭代次數(shù)，第一次計(jì)算時(shí)，k＝1。

3.2收斂判據(jù)

牛頓法計(jì)算天然氣系統(tǒng)能流的收斂判據(jù)是：max(|Δf^(k)|)≤ε(ε為牛頓法收斂精度，在步驟1.2已給出)。若滿足收斂條件，則轉(zhuǎn)至步驟3.5；若不滿足收斂條件，則繼續(xù)執(zhí)行步驟3.3。

3.3計(jì)算雅克比矩陣

計(jì)算雅克比矩陣J^(k),即

$J^{\left(k\right)}=\left[\begin{array}{cccc}\frac{\&part;f_1}{\&part;{\mathrm{\&Pi;}}_1}{|}_k& \frac{\&part;f_1}{\&part;{\mathrm{\&Pi;}}_2}{|}_k& ...& \frac{\&part;f_1}{\&part;{\mathrm{\&Pi;}}_N}{|}_k\\ \frac{\&part;f_2}{\&part;{\mathrm{\&Pi;}}_1}{|}_k& \frac{\&part;f_2}{\&part;{\mathrm{\&Pi;}}_2}{|}_k& ...& \frac{\&part;f_2}{\&part;{\mathrm{\&Pi;}}_N}{|}_k\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ \\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{\&part;f_N}{\&part;{\mathrm{\&Pi;}}_1}{|}_k& \frac{\&part;f_N}{\&part;{\mathrm{\&Pi;}}_2}{|}_k& ...& \frac{\&part;f_N}{\&part;{\mathrm{\&Pi;}}_N}{|}_k\end{array}\right]---\left(17\right)$

式中所有變量含義與式(9)和式(11)相同。

3.4計(jì)算修正量

根據(jù)步驟3.1和步驟3.3計(jì)算的不平衡量Δf^(k)和雅克比矩陣J^(k)，可以計(jì)算N個(gè)節(jié)點(diǎn)注入壓力的修正量，如下所示：

ΔΠ^(k)＝-[J^(k)]^-1Δf^(k) (18)

然后更新各節(jié)點(diǎn)壓力，公式如下：

Π^(k+1)＝Π^(k)+ΔΠ^(k) (19)

轉(zhuǎn)到步驟3.1繼續(xù)計(jì)算。

3.5計(jì)算天然氣系統(tǒng)能流分布

根據(jù)公式(7)至式(10)，計(jì)算出管道的流量，壓縮機(jī)流量以及壓縮機(jī)消耗的流量，以此得到天然氣系統(tǒng)的能流分布。