礦井風網雅可比矩陣對稱特性及并行求解模型

吳奉亮，高佳南，常心坦，李龍清

(1．西安科技大學安全學院，陜西西安　710054;2．西安科技大學教育部西部礦井開采與災害防治重點實驗室，陜西西安　710054)

?

礦井風網雅可比矩陣對稱特性及并行求解模型

吳奉亮1，2，高佳南1，2，常心坦1，2，李龍清1，2

(1．西安科技大學安全學院，陜西西安710054;2．西安科技大學教育部西部礦井開采與災害防治重點實驗室，陜西西安710054)

摘要:為了提高礦井通風網絡解算軟件的可靠性和求解大型風網時的性能，研究了通風網絡雅可比矩陣的對稱特性，引入并行計算方法求雅可比矩陣與回路風量修正值。分析了牛頓法求解礦井通風網絡的原理，發現并證明了通風網絡雅可比矩陣的對稱特性，提出用LDLT分解法求解回路風量修正值以有效減少每次迭代的時間。根據多CPU計算機的特點，研究了通風網絡雅可比矩陣以及回路風量修正值方程組的并行求解模型。采用VC語言的多線程開發技術實現了通風網絡的并行求解，通過兩個算例驗證了本模型的正確性與高效性。

關鍵詞:礦井通風網絡解算;并行計算;牛頓法;雅可比矩陣;對稱矩陣

吳奉亮，高佳南，常心坦，等．礦井風網雅可比矩陣對稱特性及并行求解模型［J］．煤炭學報，2016，41(6):1454－1459．doi:10．13225/j．cnki．jccs．2015．1302

Wu Fengliang，Gao Jia’nan，Chang Xintan，et al．Symmetry property of Jacobian matrix of mine ventilation network and its parallel calculation model［J］．Journal of China Coal Society，2016，41(6):1454－1459．doi:10．13225/j．cnki．jccs．2015．1302

通風網絡解算在礦井通風設計、管理中具有重要作用，因此具備良好性能和人機接口的礦井通風網絡解算軟件一直倍受學者關注［1－4］。網絡解算方法主要有斯拷特－恒斯雷法、牛頓法等。斯拷特－恒斯雷法來源于手工計算礦井通風網絡的需要，它在牛頓法的原理上進行了不嚴密的假設，存在迭代求解回路風量修正值不收斂的可能，為此許多學者對此進行了改進研究［5－10］，有效降低了其不收斂的可能性，但理論上仍存在無解的情況。牛頓法在理論上是嚴格的［11］，不存在無解的情況，但在求解時必須計算雅可比矩陣的逆陣，對于小型通風網絡計算時間較少，當網絡規模較大時牛頓法雅可比矩陣及其逆陣的計算量倍增，算法性能明顯下降。除了用以模擬風流在井巷中流動的規律外［12－14］，網絡解算也被用于模擬采空區流場［15－16］，對于前者風網分支數一般不超過1 000，而對于后者可形成大規模風網，其分支數可超過5 000，其對網絡解算的性能提出了更高的要求。

目前基于多CPU計算機的并行計算方興未艾［17－18］，而礦井通風網絡解算理論仍是建立在串行計算流程上的，基于此開發的軟件只需要一個CPU的支持，其他CPU都處于閑置。為此本文對礦井通風網絡的雅可比矩陣特性及其并行求解技術進行研究，以解決牛頓法在大型通風網絡解算方面的不足。

1　通風網絡雅可比矩陣的對稱特性

1.1牛頓法通風網絡解算原理

對于有 n條邊、m個節點的礦井風網，設b= n－m+1，其節點風量平衡與回路風壓平衡方程見式(1)，(2)。

其中，i為整數;cij=1表示分支j與回路i同向，cij=－1表示與回路i反向，cij=0表示分支j不在回路i中; bij=0表示節點i與分支j不相連，bij=1表示分支j風流流入i節點，bij=－1表示j分支風流流出i節點; Rj是分支 j的風阻，N·s2/m8;qj為分支 j的流量，m3/s;hfj為j分支風機風壓(是風量的函數)，Pa。根據通風網絡理論，風網分支可分為余樹支與樹支，由結點風量平衡方程(1)可以推得樹支風量都可由余樹支表示，將風網余樹枝風量向量記為Qy=(qy1，qy2，…，qyb)T，則任一分支j的風量qj可表示為

由式(3)可得風網任意分支風量qj對某余樹枝l風量qyl的偏導數為

將式(3)代入式(2)得式(5)。采用牛頓法迭代求解式(5)時，記第k次近似解為，…，，將式(5)按泰勒級數展開，并將二次及以上的高階部分忽略不計，第 k+1次線性化近似式為式(6)。

記各回路的風壓不平衡值向量F(Qy)=(f1，f2，…，fb)，回路風量修正值向量ΔQy=(Δqy1，Δqy2，…，Δqyb)T，則式(6)可整理為式(7)。

將矩陣?F/?Qy記為A，A即為雅可比矩陣。現有基于牛頓法的網絡解算都是采用式(8)對A求逆來求。第k+1次的迭代風量，此即為牛頓法解算通風網絡的迭代計算模型。當最大元素的絕對值小于給定的計算精度ε時，可被視為風網風量的準確解。

1.2雅可比矩陣的稀疏對稱特性

雅可比矩陣展開見式(9)。其矩陣的元素是式(10)在Qy=處取值。由于式(10)中i與l的取值均為(1，2，…，b)，且i和l互換后式(10)保持不變，因此可得式(11)。

式(11)表明雅可比矩陣為對稱陣。易知雅可比矩陣元素個數N為b2、非對角線元素的個數為b2－b，則考慮對稱特性時需計算的元素個數N1為

故對稱特性下，雅可比矩陣的計算量對原計算量的占比S為

從式(12)可知，當b趨于無窮時，S=50%。如取b=100時，S=50.5%，可見當風網具備一定規模時，雅可比矩陣的對稱特性可將其計算量降低約50%。

從式(10)可知，當i回路與l回路無公共分支時，可得式(10)等于0，因為若cij不等于0(j分支在i回路中)，則clj必為0(j分支必定不在l回路中)。根據通風網絡獨立回路的特性，在大型風網中兩回路無公共分支的可能性是較大的，因此大型風網的雅可比矩陣也是稀疏陣。

2　基于LDLT矩陣分解的回路風量修正值求解

文獻［19］表明當線性方程組的系數矩陣為對稱陣時，宜用LDLT分解法求解。對式(7)以為未知數的方程組AΔ=－F，使用LDLT分解法時，將對稱矩陣A分解為一個下三角矩陣L、一個對角線矩陣 D和一個上三角矩陣 LT的乘積，即A=LDLT，L，D見式(13)，(14)，各元素由式(15)～(17)確定。當L和D確定后，令DLT=z，首先由回代過程求解方程組Lz=－F而得到z，再由方程組DLT=z解出，其計算公式見式(18)～(21)。

3　通風網絡解算的并行計算模型

3.1牛頓法礦井通風網絡解算流程

根據以上結果可得牛頓法礦井通風網絡解算的流程如圖1所示，圖中劃分的各流程之間均是依賴關系，因此并行計算只能存在于各流程內部。從圖1可知，雅可比矩陣A的計算與LDLT分解法存在于回路風量修正值的迭代求解循環中(流程6～10)，其計算量占據網絡解算的主體，為此網絡解算的并行求解模型將基于雅可比矩陣與LDLT分解法的并行計算來建立。

圖1　牛頓法礦井通風網絡解算流程Fig.1　Flow chart of ventilation network calculation based on Newtown method

3.2網絡解算并行計算模型

根據式(10)可知，雅可比矩陣各元素沒有相互依賴關系，可以將所有非0元素平均分配到多個線程中并行計算。由于L矩陣對角線元素為常數1，可將D矩陣對角線元素保存于L矩陣中，如圖2所示，圖中箭頭表示由式(15)～(17)中形成的元素依賴關系。由圖2可知，D，L矩陣之間以及各矩陣元素內部都存在依賴關系，但L矩陣的同一列元素(圖2中用矩形框出)沒有依賴關系，可以并行計算。根據以上分析將圖1中的第6～8步細化為圖3所示流程。圖3中第6.1～6.3步完成雅可比矩陣的并行計算:主線程首先獲取CPU個數C，并將雅可比矩陣非0

圖2　D，L矩陣元素依賴關系Fig.2　Dependence of the elements of the D and L matrixes

圖3　雅可比矩陣與L矩陣的并行求解流程Fig.3　Flow chart of parallel calculation of the Jacobian matrix and L matrix

元素均分成C組(當無法均分時，每組之間的元素個數相差不大于1)，然后創建C個子線程并行計算各組元素，主線程等待各子線程結束后，進入第7步。從7.1～7.9完成L矩陣的并行計算:從L矩陣的第1列(i=1)開始，主線程首先計算D矩陣第i個對角線元素dii，然后判斷L矩陣第i列需要計算元素個數b－i是否小于C，如果是則取C=b－i，之后將第i列待求元素分成C組，再創建C個子線程并行計算各組元素，主線程等待各子線程均結束后，再開始計算下一列，如此循環，直到所有列均完成后，進入第8步。綜上，A，L矩陣并行計算時子線程之間無信息交互，并發性高，其并行計算是網絡解算并行計算模型的核心。

4　算例分析

在已有的網絡解算軟件基礎之上，采用VC++實現了以上模型，現通過一簡單風網與一大型風網分別對模型的正確性與性能進行驗證。

4.1簡單風網算例

圖4風網具有12條邊，7個節點，分支風阻均取0.1 N·s2/m8，編號12是風機分支，風機性能函數hf=3 204.2－4.91q12+0.025q122－0.000 2q123，有效風量范圍為［50，220］。分支7～12是6條余樹枝，故Qy=(qy1，qy2，…，qy6)T=(q7，q8，q9，q10，q11，q12)T，各回路編號及所含分支分別為1(7，－6，－5，－4，2，3)，2(8，－5，－4，2)，3(9，－3，－2，4)，4(10，－4，2)，5 (11，－3，－2，4，5)，6(12，1，4，5，6)，分支號為負表示其方向與余樹支方向相反。分支12的風量初值取135 m3/s，其他余樹支風量初值取1 m3/s，則分支1～12風量初值依次為(135.0，－1.0，1.0，134.0，134.0，134.0，1.0，1.0，1.0，1.0，1.0，135.0)。以第1次迭代計算雅可比矩陣a23，a32元素為例，按式(3)寫出2，3回路中樹支風量的計算式分別為:q5=－q8－q7+q11+q12，q4=－q8－q7+q9－q10+q11+q12，q3=q7－q9－q11，q2=q7+q8－ q9+q10－q11;回路風壓代數和函數分別為:f2=R8－。易知均為0，故

可見此對稱位置兩元素相等，經驗證所有對稱位置的元素均相等。采用運行于4核CPU中的程序計算，需計算雅可比矩陣元素總數為21，軟件創建的第1個線程計算6個元素，其他3個線程各計算5個元素，得到整個雅可比矩陣見式(22)，對應圖2中的D，L矩陣元素見式(23)。

圖4 簡單風網Fig.4　A simple ventilation network sample

經過7次迭代后求得各分支風量向量為(129.4，64.7，64.7，32.3，32.3，32.3，32.3，0，0，64.7，64.7，129.4)。

4.2大型風網實例

引用文獻［15］中使用網絡解算求采空區流場的實例，如圖5所示，采煤工作面長度為120 m，采空區冒落帶走向長度為200 m，采空區被劃分為許多縱橫交錯的濾流分支，濾流分支特有的阻力特性方程［15－16］將采空區多孔介質轉換成了網絡系統，可用網絡解算方法計算其流場。根據其原理，濾流分支長度越小模擬精度越高。圖5是將濾流橫向長度取為2 m、縱向取為4 m形成的風網，其分支數達到6 135，節點數為 3 133，獨立回路數 3 003。采用 4核CPU(頻率2.5 GHz)、在Windows 7操作系統下解算以上風網，得到以下結果:

圖5　網絡化采空區后的大型風網Fig.5　A big ventilation network in a gob area divided by branches

(1)雅可比矩陣共有9 012 004個元素，其中非0元素880 402個，僅占總數的9.8%。

(2)考慮雅可比矩陣對稱性，但不考慮稀疏性時，單線程計算時間是4 min，多線程并行計算時間減至94 s。考慮對稱性與稀疏性(只計算非0元素)、采用并行計算，雅可比矩陣的計算時間減至8 s。

(3)單線程LDLT分解法求一次回路風量修正值計算時間為13 min，基于多線程并行計算時間為4 min。

圖6是程序運行時4個CPU的使用率。每個CPU的曲線均分為3段:第1段是軟件單線程運行圖1中第1～5步，CPU的平均使用率不到30%;第2段是計算雅可比矩陣，時間為8 s左右，CPU使用率均為100%;第3段是運行圖3中7.2～7.9步，曲線的振蕩是循環中不斷的創建線程與退出線程形成的，但CPU的平均使用率均超過90%。整個解算過程迭代了25次收斂，共耗時約110 min。

圖6　基于并行計算的網絡解算CPU使用率Fig.6　The usage rate of 4 CPUs when parallel calculating a ventilation network

按25次迭代、每次迭代耗時38 min計算，普通牛頓法將耗時15.8 h。可見本文研究對于大型風網的解算具有明顯優勢。

5　結　　論

(1)礦井通風網絡雅可比矩陣是稀疏對稱陣，采用LDLT矩陣分解求回路風量修正值具有更高的性能。通過分析基于牛頓法的礦井通風網絡解算原理，發現了雅可比矩陣的對稱特性。用LDLT矩陣分解法求以對稱雅可比矩陣為系數矩陣的回路風量修正值方程組，提高了牛頓法解算礦井通風網絡的性能。通過一簡單算例對雅可比矩陣的對稱性及LDLT分解法應用于網絡解算的正確性進行了驗證。

(2)礦井通風網絡解算的回路風量修正值求解具有顯著的并行特征。網絡解算主要耗時于回路風量修正值的迭代計算，每一次迭代包括求雅可比矩陣與解以回路風量修正值為未知數的線性方程組。雅可比矩陣、L矩陣并行計算過程中，子線程之間無信息交互，并發性高。通過實現雅可比矩陣與L矩陣的并行計算，構建了礦井通風網絡解算的并行求解模型。

(3)與普通網絡解算軟件相比，基于牛頓法并行計算的網絡解算在當前多CPU個人計算機上性能大幅提升。通過對一含有6 135條分支的風網進行驗證，軟件可以在2 h左右給出常規網絡解算難以求解的結果。

參考文獻:

［1］王德明，周福寶．基于WINDOWS的礦井通風網絡解算軟件的研制［J］．中國礦業大學學報，2000，29(1):41－44．Wang Deming，Zhou Fubao．Study on Windows-Based software for mine ventilation network solution［J］．Journal of China University of Mining＆Technology，2009，29(1):41－44．

［2］劉劍．流體網絡理論［M］．北京:煤炭工業出版社，2002．

［3］王德明，王俊，周福寶．基于面向對象技術開發的礦井通風圖形系統［J］．煤炭學報，2000，25(5):510－513．Wang Deming，Wang Jun，Zhou Fubao．Developed mine ventilation graphics system based on the object-oriented technologies［J］．Journal of China Coal Society，2000，25(5):510－513．

［4］吳奉亮，周澎，李暉，等．基于智能對象的通風CAD模型研究［J］．煤炭科學技術，2009，37(5):54－57．Wu Fengliang，Zhou Peng，Li Hui，et al．Research on the model of mine ventilation CAD based on the intelligent object［J］．Coal Science and Technology，2009，37(5):54－57．

［5］周心權，吳兵．定向圈劃回路法探討［J］．煤礦安全，1993(4): 10－12．Zhou Xinquan，Wu Bing．The exploration of directional circle loop method［J］．Safety in Coal Mines，1993(4):10－12．

［6］魏連江，周福寶，朱華新．通風網絡拓撲理論及通路算法研究［J］．煤炭學報，2008，33(8):926－930．Wei Lianjiang，Zhou Fubao，Zhu Huaxin．Topology theory of ventilation network and path algorithm［J］．Journal of China Coal Society，2008，33(8):926－930．

［7］孫臣良，題正義，趙鐵文．基于改進Cross算法的礦井復雜風網可視化解算系統［J］．遼寧工程技術大學學報(自然科學版)， 2008，28(S1):25－27．Sun Chenliang，Ti Zhengyi，Zhao Tiewen．Visual calculation system of complicated ventilation network of coal mine based on improvedcrossalgorithm［J］．JournalofLiaoningTechnical University(Natural Science)，2008，28(S1):25－27．

［8］鐘德云，王李管，畢林，等．基于回路風量法的復雜礦井通風網絡解算算法［J］．煤炭學報，2015，40(2):365－370．Zhong Deyun，Wang Liguan，Bi Lin，et al．Algorithm of complex ventilation network solution based on circuit air quantity method［J］．Journal of China Coal Society，2015，40(2):365－370．

［9］馬心校．解算礦井通風網絡不收斂的一種處理方法［J］．煤炭工程師，1990(3):8－17．Ma Xinxiao．Method of solving the problem of not convergence in mineventilationnetworkcalculation［J］．CoalEngineer，1990(3):8－17．

［10］楊鑫祥，肖偉，武猛猛，等．礦井通風網絡解算軟件的改進研究與實現［J］．礦業研究與開發，2015，35(4):80－83．Yang Xinxiang，Xiao Wei，Wu Mengmeng，et al．Improvement and realization of the software for mine improvement and realization of the software for mine ventilation network solution［J］．Mining R．＆D．，2015，35(4):80－83．

［11］李恕和，王義章．礦井通風網絡計算的牛頓法［J］．煤炭學報，1982，7(4):52－62．Li Shuhe，Wang Yizhang．The newtow method for calculating the mine ventilation networks［J］．Journal of China Coal Society，1982，7(4):52－62．

［12］張國樞．通風安全學［M］．徐州:中國礦業大學出版社，2000．

［13］吳奉亮．對角式通風系統封閉型網絡降阻分析［J］．煤礦安全，2011，42(7):135－138．Wu Fengliang．Pressure reduction analysis of closed network in a diagonal ventilation system［J］．Safety in Coal Mines，2011，42(7):135－138．

［14］El-Nagdy K A，Shoaib A M．Alternate solutions for mine ventilation network to keep a preassigned fixed quantity in a working place［J］．International Journal of Coal Science＆ Technology，2015，2(4):269－278．

［15］劉澤功．通風安全工程計算機模擬與預測［M］．北京:煤炭工業出版社，1995:172－180．

［16］魏引尚，邵煒航，王慧．基于網絡解算的采空區漏風分布計算模式［J］．煤礦安全，2011，42(1):46－49．Wei Yingshang，Shao Weihang，Wang Hui．Calculation model of gob air leakage based on the ventilation network calculation ［J］．Safety in Coal Mines，2011，42(1):46－49．

［17］陳國良．并行算法的設計與分析［M］．北京:高等教育出版社，2009．

［18］谷照升．基于多核CPU的并行計算設計［J］．長春工程學院學報(自然科學版)，2009，10(3):92－94．Gu Zhaosheng．Design for Multi-kernel CPU based parallel computing［J］．Journal of Changchun Institute of Technology(Natural Sciences Edition)，2009，10(3):92－94．

［19］徐士良．C常用算法程序集(第二版)［M］．北京:清華大學出版社，2004．

中圖分類號:TD722

文獻標志碼:A

文章編號:0253－9993(2016)06－1454－06

收稿日期:2015－09－10修回日期:2015－12－08責任編輯:畢永華

基金項目:國家自然科學基金青年基金資助項目(51204135，51104116)

作者簡介:吳奉亮(1977—)，男，四川新都人，副教授，碩士生導師。E－mail:15038537@qq．com

Symmetry property of Jacobian matrix of mine ventilation network and its parallel calculation model

WU Feng-liang1，2，GAO Jia-nan1，2，CHANG Xin-tan1，2，LI Long-qing1，2
(1．School of Safety Engineering，Xi’an University of Science and Technology，Xi’an710054，China;2．Key Laboratory of Western Mine＆Hazard Prevention，China Ministry of Education，Xi’an University of Science and Technology，Xi’an710054，China)

Abstract:A parallel calculation and symmetry property of the Jacobian matrix deduced from a ventilation network is adopted in this study to enhance the reliability and the performance of software on calculating Coalmine’s ventilation network．The theory of solving Coalmine’s ventilation network by using Newtown method is presented in some details，and then the symmetry characteristic of the Jacobian matrix of a ventilation network is illustrated．A LDLTdecomposition method is put forward to solve the corrected value of flow rate in a circuit，which can save a lot of time during each iteration．To take the most advantages of a multi-CPU computer，a parallel model is employed to solve the Jacobian matrix and the equations of corrected values in circuits of a ventilation network．This parallel model is achieved using the technology of multi-threads in VC language，whose validity and performance are verified by running it in two ventilation network samples．

Key words:mine ventilation network calculation;parallel calculation;Newtown method;Jacobian matrix;symmetric matrix