利用同倫法解算自然分風網絡

趙自豪 任玉輝 陳世江

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

利用同倫法解算自然分風網絡

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

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

為了解決牛頓法解算自然分風網絡時對初值的依賴性，并解決利用諸如Matlab和VB等編程語言進行風網解算時所付出的版權費用，利用自由軟件Octave對解非線性方程組的常用算法——同倫法進行了討論，論證了其在解算自然分風網絡時的有效性。

Octave軟件 通風網絡 同倫法

目前針對地下礦通風網絡解算的軟件有很多，但是大多是由科研院所或專業(yè)的軟件開發(fā)公司開發(fā)的受版權保護的商業(yè)軟件，如Mvent，Avent等，這些軟件使用時要收費，且其解算模塊是固定的，用戶不能根據自己的需要自行調整、刪減。因此，筆者尋找一種能夠實現快速開發(fā)、快速解算，且是免費的通風網絡解算方案。本文以同倫法為例介紹了利用Octave軟件開發(fā)通風解算軟件的可行性及方便性。

1 Octave軟件簡介

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

2 同倫法解算原理

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

式中：C——通風網絡獨立回路矩陣；

CT——通風網絡獨立回路矩陣的轉置；

Rdiag——各分支風阻向量的對角陣形式；

Q——各分支的風量列向量；

Qy——各余支的風量列向量，為待求量；

Hf——風機風壓列向量；

HN——自然風壓列向量。

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

同倫算法的思想是人為地引進一個參數t，構造一個函數族H（x，t），使得：

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

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

對式（3）兩邊求t的偏導，整理得：

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

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

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

式（6）中Jacob為雅可比矩陣，對于通風網絡回路能量平衡方程，其雅可比矩陣為

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

3 同倫法解算自然分風網絡的Octave實現

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

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

%采用同倫算法解算自然分風網絡，函數各輸入參數意義如下：

%c為獨立回路矩陣，r為各分支風阻，hf為通風機風壓，x0為獨立分支的初值，N為積分步距倒數，該數值越大，則積分步距越小，結果越精確。

%本函數輸出為各分支的自然分風量。

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

%jq為回路風壓平衡方程的雅可比行列式。

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

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

4 解算實例

以圖1所示的通風網絡為例。該網絡中風阻分別為：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，風機風壓為500Pa。

圖1 通風網絡

在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]；%本風網的獨立回路矩陣r＝[.5.005.004.005.01.02 0]＇；%本風網各分支的風阻列向量

hf＝[0 0 0 0 0 0 500]＇；%本風網的風機風壓列向量

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

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

即可求解出本風網各分支的風量。在本例中，

解算結果為[7.5201，126.1737 106.9637，233.1374，133.6938，99.4437，233.1374]＇。該結果為列向量，向量中的數據分別對應圖1中1～7號分支的風量，單位為m3／s。

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

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

圖2 初值偏差－結果偏差圖

從圖2可以看出，經過20步的積分計算，對任意初值，結果偏差遠小于初值偏差，表明該計算方法能夠有效地逼近精確結果。同時也表明，初值偏差越大，經過相同的積分步數，結果偏差也相應越大。

分別取積分步數從10到200，依次重復上述過程，取每個積分步數下對應的結果偏差的平均值。繪制出步數－平均結果偏差圖像，見圖3。

圖3 積分步數－平均結果偏差圖

從圖3可以看出，隨著積分步數的增多，平均結果偏差逐漸逼近0。從計算結果及現場實踐可知，當積分步數為100時，平均結果偏差為0.013411，已經達到工程實用的精度。故利用同倫法解算時，采用積分步數為100完全可以滿足精度需求。

6 結論

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

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

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

[3] 蔡大用.現代科學計算[M].北京：科學出版社，2000

[4] 張國樞.通風安全學[M].北京：中國礦業(yè)大學出版社，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－），男，山東臨沂人，講師，現為中國礦業(yè)大學（北京）在讀博士生，主要研究領域為煤礦瓦斯防治及礦山防滅火。

（責任編輯 梁子榮）