趙自豪 任玉輝 陳世江

（1.內(nèi)蒙古科技大學(xué)礦業(yè)工程學(xué)院，內(nèi)蒙古自治區(qū)包頭市，014010；2.中國礦業(yè)大學(xué)（北京）資源與安全工程學(xué)院，北京市海淀區(qū)，100083）

利用同倫法解算自然分風(fēng)網(wǎng)絡(luò)

趙自豪1，2任玉輝1，2陳世江1

（1.內(nèi)蒙古科技大學(xué)礦業(yè)工程學(xué)院，內(nèi)蒙古自治區(qū)包頭市，014010；2.中國礦業(yè)大學(xué)（北京）資源與安全工程學(xué)院，北京市海淀區(qū)，100083）

為了解決牛頓法解算自然分風(fēng)網(wǎng)絡(luò)時(shí)對初值的依賴性，并解決利用諸如Matlab和VB等編程語言進(jìn)行風(fēng)網(wǎng)解算時(shí)所付出的版權(quán)費(fèi)用，利用自由軟件Octave對解非線性方程組的常用算法——同倫法進(jìn)行了討論，論證了其在解算自然分風(fēng)網(wǎng)絡(luò)時(shí)的有效性。

Octave軟件 通風(fēng)網(wǎng)絡(luò) 同倫法

目前針對地下礦通風(fēng)網(wǎng)絡(luò)解算的軟件有很多，但是大多是由科研院所或?qū)I(yè)的軟件開發(fā)公司開發(fā)的受版權(quán)保護(hù)的商業(yè)軟件，如Mvent，Avent等，這些軟件使用時(shí)要收費(fèi)，且其解算模塊是固定的，用戶不能根據(jù)自己的需要自行調(diào)整、刪減。因此，筆者尋找一種能夠?qū)崿F(xiàn)快速開發(fā)、快速解算，且是免費(fèi)的通風(fēng)網(wǎng)絡(luò)解算方案。本文以同倫法為例介紹了利用Octave軟件開發(fā)通風(fēng)解算軟件的可行性及方便性。

1 Octave軟件簡介

Octave是一款用于數(shù)值計(jì)算和繪圖的開源軟件，除一般的數(shù)學(xué)運(yùn)算外，尤其精于矩陣運(yùn)算、求解聯(lián)立方程組、計(jì)算矩陣特征值和特征向量等，而在許多的工程實(shí)際問題中，數(shù)據(jù)都可以用矩陣或向量表示出來從而轉(zhuǎn)化為對這類矩陣的求解。Octave能夠通過多種形式實(shí)現(xiàn)數(shù)據(jù)可視化，Octave本身也是一門編程語言，所以其擴(kuò)展性非常強(qiáng)。由于Octave與科研項(xiàng)目中普遍使用的Matlab基本兼容，因此其易用性得到了廣泛的認(rèn)可。

2 同倫法解算原理

眾所周知，通風(fēng)網(wǎng)絡(luò)的解算實(shí)際上就是解算通風(fēng)網(wǎng)絡(luò)回路能量平衡方程和節(jié)點(diǎn)風(fēng)量平衡方程，見式（1）、式（2），以上兩式均為矩陣形式。

式中：C——通風(fēng)網(wǎng)絡(luò)獨(dú)立回路矩陣；

CT——通風(fēng)網(wǎng)絡(luò)獨(dú)立回路矩陣的轉(zhuǎn)置；

Rdiag——各分支風(fēng)阻向量的對角陣形式；

Q——各分支的風(fēng)量列向量；

Qy——各余支的風(fēng)量列向量，為待求量；

Hf——風(fēng)機(jī)風(fēng)壓列向量；

HN——自然風(fēng)壓列向量。

式（1）是一個(gè)多元二次方程組，傳統(tǒng)上，一般多采用牛頓（Newton）迭代法和它的變體進(jìn)行求解，但該方法的收斂半徑很小，需要很好的初值x0。而在20世紀(jì)70年代發(fā)展出了一種同倫算法，它有大范圍的收斂性，能夠?qū)⒊踔档娜菰S范圍顯著擴(kuò)大。

同倫算法的思想是人為地引進(jìn)一個(gè)參數(shù)t，構(gòu)造一個(gè)函數(shù)族H（x，t），使得：

H（x，0）＝F0（x），H（x，1）＝F（x）假設(shè)F0（x）＝0的解已知，從t＝0出發(fā)可以求解

對于t∈[0，1]，假設(shè)式（3）的解為x（t）。如果x（t）可以形成一條光滑曲線，其起點(diǎn)x（0）為F0（x）＝0的解，據(jù)假設(shè)它是已知的，曲線x（t）的終點(diǎn)x（1）就是F（x）＝0的解。H（x，t）稱為一個(gè)同倫，它的解x（t）稱為同倫曲線。

對式（3）兩邊求t的偏導(dǎo)，整理得：

式（4）即為x（t）的定解問題，很明顯，這是常微分方程的初值問題，其解法張國樞在參考文獻(xiàn)4中有清楚表述，不再贅述。

構(gòu)造H（x，t）的方法有很多，這里針對式（1）做如下構(gòu)造：

將式（5）代入式（4），得：

式（6）中Jacob為雅可比矩陣，對于通風(fēng)網(wǎng)絡(luò)回路能量平衡方程，其雅可比矩陣為

對式（6）進(jìn)行積分，積分區(qū)間為[0，t]，進(jìn)一步換算即可得Qy（1），此即通風(fēng)網(wǎng)絡(luò)回路能量平衡方程的解。

3 同倫法解算自然分風(fēng)網(wǎng)絡(luò)的Octave實(shí)現(xiàn)

為簡潔起見，本例僅演示風(fēng)機(jī)風(fēng)壓恒定，不考慮自然風(fēng)壓情況下風(fēng)網(wǎng)的同倫算法解算工具的開發(fā)情況。

其具體求解代碼及說明文字如下：

%采用同倫算法解算自然分風(fēng)網(wǎng)絡(luò)，函數(shù)各輸入?yún)?shù)意義如下：

%c為獨(dú)立回路矩陣，r為各分支風(fēng)阻，hf為通風(fēng)機(jī)風(fēng)壓，x0為獨(dú)立分支的初值，N為積分步距倒數(shù)，該數(shù)值越大，則積分步距越小，結(jié)果越精確。

%本函數(shù)輸出為各分支的自然分風(fēng)量。

%下文中的fq為回路風(fēng)壓平衡方程的矩陣形式。

%jq為回路風(fēng)壓平衡方程的雅可比行列式。

%下面為利用同倫法推導(dǎo)出的式6的積分求算形式

將以上代碼存儲到Octave安裝目錄下的Work文件夾，并命名為tl＿method.m，即可在Octave交互環(huán)境下實(shí)現(xiàn)對這一函數(shù)的調(diào)用。

4 解算實(shí)例

以圖1所示的通風(fēng)網(wǎng)絡(luò)為例。該網(wǎng)絡(luò)中風(fēng)阻分別為：r1＝0.5kg／m7，r2＝0.005kg／m7，r3＝0.004kg／m7，r4＝0.005kg／m7，r5＝0.001kg／m7，r6＝0.02kg／m7，r7＝0kg／m7，風(fēng)機(jī)風(fēng)壓為500Pa。

圖1 通風(fēng)網(wǎng)絡(luò)

在Octave交互環(huán)境下鍵入如下命令：

c＝[1－1 1 0 0 0 0；0 1 0 1 1 0 1；0 0 1 1 0 1 1]；%本風(fēng)網(wǎng)的獨(dú)立回路矩陣r＝[.5.005.004.005.01.02 0]＇；%本風(fēng)網(wǎng)各分支的風(fēng)阻列向量

hf＝[0 0 0 0 0 0 500]＇；%本風(fēng)網(wǎng)的風(fēng)機(jī)風(fēng)壓列向量

x0＝[10 30 60]＇%各余枝的風(fēng)量初值，為任意值

tl＿method（c，r，hf，x0，100）

即可求解出本風(fēng)網(wǎng)各分支的風(fēng)量。在本例中，

解算結(jié)果為[7.5201，126.1737 106.9637，233.1374，133.6938，99.4437，233.1374]＇。該結(jié)果為列向量，向量中的數(shù)據(jù)分別對應(yīng)圖1中1～7號分支的風(fēng)量，單位為m3／s。

5 解算效率及初值選擇對結(jié)果精度的影響

為了測定不同的初值選擇對解算精度的影響，仍以圖1所示的通風(fēng)網(wǎng)絡(luò)為例，假定在初值x0＝[10，30，60]＇，迭代次數(shù)為1000時(shí)的輸出為標(biāo)準(zhǔn)結(jié)果，記為Q。取Q中對應(yīng)的余枝風(fēng)量為標(biāo)準(zhǔn)初值，記為X。計(jì)算機(jī)任意生成初值xi（初值數(shù)值在－1000～1000區(qū)間），取N＝20，代入調(diào)用同倫法解算模塊，其輸出為qi，記初值偏差sum（（X－xi）.＾2）／sum（X.＾2）為橫坐標(biāo)，結(jié)果偏差sum（（Q－qi）.＾2）／sum（Q.＾2）為縱坐標(biāo)，多次重復(fù)這一過程，生成的初值偏差－結(jié)果偏差圖像見圖2。

圖2 初值偏差－結(jié)果偏差圖

從圖2可以看出，經(jīng)過20步的積分計(jì)算，對任意初值，結(jié)果偏差遠(yuǎn)小于初值偏差，表明該計(jì)算方法能夠有效地逼近精確結(jié)果。同時(shí)也表明，初值偏差越大，經(jīng)過相同的積分步數(shù)，結(jié)果偏差也相應(yīng)越大。

分別取積分步數(shù)從10到200，依次重復(fù)上述過程，取每個(gè)積分步數(shù)下對應(yīng)的結(jié)果偏差的平均值。繪制出步數(shù)－平均結(jié)果偏差圖像，見圖3。

圖3 積分步數(shù)－平均結(jié)果偏差圖

從圖3可以看出，隨著積分步數(shù)的增多，平均結(jié)果偏差逐漸逼近0。從計(jì)算結(jié)果及現(xiàn)場實(shí)踐可知，當(dāng)積分步數(shù)為100時(shí)，平均結(jié)果偏差為0.013411，已經(jīng)達(dá)到工程實(shí)用的精度。故利用同倫法解算時(shí)，采用積分步數(shù)為100完全可以滿足精度需求。

6 結(jié)論

本文通過論證證明，同倫法可以很好地應(yīng)用于通風(fēng)網(wǎng)絡(luò)解算，且具有很好的初值容許范圍，能夠克服牛頓法的初值選擇方面的缺點(diǎn)。同時(shí)，Octave軟件是免費(fèi)的、開源的，這就意味著利用該軟件進(jìn)行通風(fēng)網(wǎng)絡(luò)解算具有非常大的靈活性。Octave軟件具有卓越的矩陣運(yùn)算能力，而通風(fēng)網(wǎng)絡(luò)的解算本質(zhì)上就是矩陣的解算，所以利用Octave能很好的完成通風(fēng)網(wǎng)絡(luò)解算的任務(wù)。

[1] Soren Hauberg.Gnu Octave Manual Version 3[M].NETWORK THEORY LTD，2008

[2] 王德明.礦井通風(fēng)與安全[M].北京：中國礦業(yè)大學(xué)出版社，2007

[3] 蔡大用.現(xiàn)代科學(xué)計(jì)算[M].北京：科學(xué)出版社，2000

[4] 張國樞.通風(fēng)安全學(xué)[M].北京：中國礦業(yè)大學(xué)出版社，2004

Calculation of natural distribution of ventilation network using homotopy method

Zhao Zihao1，2，Ren Yuhui1，2，Chen Shijiang1

（1.School of Mining Engineering，Inner Mongolia University of Science and Technology，Baotou，Inner Mongolia 014010，China；2.School of Resources and Safety Engineering，China University of Mining Technology Beijing，Haidian，Beijing 100083，China）

In order to solve the reliance of Newton method on initial value when calculating the natural ventilation network and to solve the payment of copyright fees for using Matlab，VB and other programming languages，the cost－free Octave software was used to discuss the common algorithm，namely homotopy method，used for solving non－linear equation system.The effectiveness of homotopy method on solving natural ventilation network was demonstrated too.

Octave，ventilation network，homotopy method

TD725

A

趙自豪（1977－），男，山東臨沂人，講師，現(xiàn)為中國礦業(yè)大學(xué)（北京）在讀博士生，主要研究領(lǐng)域?yàn)槊旱V瓦斯防治及礦山防滅火。

（責(zé)任編輯 梁子榮）