球磨機磨礦系統(tǒng)多變量辨識建模的研究

蔡改貧，許 琴，熊 洋，楊麗榮

（江西理工大學(xué) 機電工程學(xué)院，江西 贛州 341000）

球磨機磨礦系統(tǒng)多變量辨識建模的研究

蔡改貧，許 琴，熊 洋，楊麗榮

（江西理工大學(xué) 機電工程學(xué)院，江西 贛州 341000）

針對球磨機磨礦系統(tǒng)難以理論推導(dǎo)出磨礦系統(tǒng)的數(shù)學(xué)模型問題，應(yīng)用最小二乘法估計的多變量ARX模型辨識的算法，給出了球磨機系統(tǒng)模型的離散形式。同時結(jié)合試驗中球磨機的內(nèi)部負(fù)荷參數(shù)（介質(zhì)填充率、料球比、磨礦濃度）與功率消耗建立出磨礦系統(tǒng)的動態(tài)數(shù)學(xué)模型，詳細(xì)介紹了利用MATLAB辨識工具箱建模方法及其辨識的實現(xiàn)過程。通過模型輸出曲線與實際輸出曲線的對比，表明該模型具有較高的準(zhǔn)確性，為球磨機實現(xiàn)高效低能過程控制提供理論依據(jù)。

最小二乘法；多變量ARX模型；辨識建模；球磨機負(fù)荷；功率消耗；系統(tǒng)辨識

0 引言

鎢在全球范圍內(nèi)被列為重要的戰(zhàn)略金屬，隨著近幾十年對鎢的不斷開采，我國鎢礦資源消耗過快，貧礦比例大、難選礦石多，礦石的綜合利用水平低，各企業(yè)相繼提出提高鎢礦資源的回收及加強礦山企業(yè)節(jié)能降耗主題[1-2]。球磨機的磨礦作業(yè)在鎢選礦工藝流程中是非常關(guān)鍵的環(huán)節(jié)，合理的磨礦細(xì)度是選別作業(yè)中獲得高品位，高回收率鎢精礦的重要因素，磨礦效率和能耗直接影響選礦廠選別技術(shù)指標(biāo)[3-4]。磨礦作業(yè)是整個鎢礦選礦廠動力消耗和金屬材料消耗最大的作業(yè)，維護和修理費用也比較高，據(jù)相關(guān)數(shù)據(jù)統(tǒng)計，其生產(chǎn)費用約占全廠總費用的30%以上[5]，因此降低磨礦成本已顯得尤為迫切。在鎢礦選礦過程中，球磨機磨礦作業(yè)的好壞，對于提高生產(chǎn)能力、產(chǎn)品質(zhì)量、回收率和降低生產(chǎn)成本都起著重要的作用。球磨機負(fù)荷（包括鋼球負(fù)荷、物料負(fù)荷以及水量）是磨礦過程的一個重要參數(shù)[6]，而礦石屬性、磨機結(jié)構(gòu)參數(shù)以及磨機轉(zhuǎn)速率在磨礦生產(chǎn)過程中一般比較穩(wěn)定，因此，合理選擇和優(yōu)化球磨機內(nèi)部負(fù)荷參數(shù)，是提高球磨機磨礦效率和降低能耗的關(guān)鍵。

目前，許多研究學(xué)者針對球磨機內(nèi)部負(fù)荷檢測及控制已開展了大量研究，如利用神經(jīng)網(wǎng)絡(luò)建立球磨機內(nèi)部負(fù)荷與外部響應(yīng)之間的關(guān)系模型以實時地預(yù)測球磨機的內(nèi)部負(fù)荷參數(shù)[7]，在球磨機軟測量技術(shù)上提出支持向量機結(jié)合模糊控制實現(xiàn)球磨機優(yōu)化運行[8]，以及綜合多種智能算法將專家控制、模糊控制與自尋優(yōu)算法相結(jié)合對球磨機負(fù)荷進行控制[9]，以上研究工作都取得較好預(yù)測與控制效果，但是它們并沒有給出內(nèi)部負(fù)荷與外部狀態(tài)之間的數(shù)學(xué)模型。近幾年，隨著選礦控制技術(shù)的發(fā)展，有關(guān)科研部門將系統(tǒng)辨識技術(shù)應(yīng)用到一些先進控制中，如國外學(xué)者提出了系統(tǒng)辨識軟測量技術(shù)應(yīng)用到選礦工業(yè)控制中[10]。系統(tǒng)辨識作為以試驗數(shù)據(jù)建立系統(tǒng)模型的方法，依據(jù)系統(tǒng)辨識的方法能對磨礦過程進行描述，可以獲得比較理想的磨礦系統(tǒng)數(shù)學(xué)模型，有望用于磨礦過程的控制。

在磨礦過程中控制好能耗，將有利于我們合理利用資源，降低能耗，為此，利用系統(tǒng)辨識的方法，探索球磨機內(nèi)部負(fù)荷參數(shù)與能量消耗之間的內(nèi)在聯(lián)系，指導(dǎo)磨礦設(shè)備實現(xiàn)節(jié)能降耗，對降低能耗、節(jié)約能源具有非常重要的現(xiàn)實意義。

1 球磨機磨礦過程

某選礦廠球磨機現(xiàn)場工況如圖1所示。球磨機主要由進出料裝置、主軸承部件、傳動部件、轉(zhuǎn)動部件、減速部件和電動機等組成，球磨機主體是一個長形圓筒，筒體一端是給料入口，另一端是磨礦后的礦物和水混合物的出口，筒內(nèi)裝有直徑大小不一的鋼球，不同型號球磨機內(nèi)壁襯上裝有不同形狀的錳鋼襯板。球磨機在運行過程中，在離心力和摩擦力作用下，礦石和鋼球由襯板提升到一定高度，然后在自身重力作用下沿拋物線下落。礦石主要受鋼球沖擊破碎，同時受鋼球與鋼球間、鋼球與襯板間的擠壓與研磨。

圖1 球磨機現(xiàn)場工況Fig.1 Field working condition diagram of ball mill

球磨機系統(tǒng)是一個典型的多變量、強耦合、動態(tài)特性復(fù)雜的系統(tǒng)[11]。球磨機是否控制在最優(yōu)工作狀態(tài)直接影響磨礦質(zhì)量和球磨機工作效率。由前期大量磨礦試驗分析可知，影響球磨機磨礦生產(chǎn)過程的因素比較多，如礦石屬性、磨機結(jié)構(gòu)參數(shù)、磨礦操作因素等，難以建立或推導(dǎo)出精確的數(shù)學(xué)模型。由于磨礦過程中磨機結(jié)構(gòu)一旦確定，參數(shù)將不會改變，生產(chǎn)過程中磨機轉(zhuǎn)速率也一般保持不變，而介質(zhì)填充率、料球比、磨礦濃度是經(jīng)常變化的，由此，可將介質(zhì)填充率、料球比和磨礦濃度作為多輸入?yún)?shù)，功率消耗作為輸出參數(shù)。系統(tǒng)有較強的耦合性，當(dāng)其中三個輸入?yún)?shù)中任何一個變化都將引起輸出參數(shù)變化，通過對上述參數(shù)的檢測和最優(yōu)調(diào)整就能實現(xiàn)球磨機磨礦系統(tǒng)工作在最佳工況。為提高效率，降低能耗，使球磨機運行的最佳工作狀態(tài)，本文通過人工加入不同比例的礦料、鋼球和水量，并檢測球磨機的電機功率，應(yīng)用多變量辨識方法建立介質(zhì)充填率、料球比和磨礦濃度與對應(yīng)功率消耗之間的球磨機數(shù)學(xué)模型。

2 球磨機系統(tǒng)辨識控制模型

2.1 系統(tǒng)辨識理論

球磨機磨礦系統(tǒng)是一個復(fù)雜的多變量系統(tǒng)，筆者應(yīng)用多變量系統(tǒng)辨識方法建立用于控制描述球磨機磨礦動態(tài)特性的數(shù)學(xué)模型。多變量系統(tǒng)[12]可以看作是單輸入單輸出系統(tǒng)的直接推廣和延伸。理想的多輸入多輸出過程不包含噪聲，在現(xiàn)實的多輸入多輸出系統(tǒng)中，由于外部環(huán)境和檢測精度等因素使系統(tǒng)的輸出有噪聲存在。用u（k）、r（k）表示系統(tǒng)的輸入向量與理想輸出向量，ω（k）表示系統(tǒng)噪聲，y（k）表示系統(tǒng)的實際輸出向量，圖2可以來描述實際的包含噪聲的多變量過程。

圖2 包含噪聲的多變量過程Fig.2 Noise-containing multivariable process

多變量系統(tǒng)理想的輸出響應(yīng)可描述成[13]：

其中G（z）為傳遞函數(shù)矩陣，也即多變量系統(tǒng)模型。因此，包含噪聲的多變量過程可表示為：

傳遞函數(shù)矩陣中，每個元素Gij（z）（i=1，2，…，r）都是z的有理真分式函數(shù)，并且Gij（z）的分母多項式是系統(tǒng)的特征方程，可用A*（z）來表示，若Gij（z）的分子用B*ij（z）（i，j=1，2，…，r）來表示，則Gij（z）可表示成：

其中：

若將式（2）右邊各元素分子分母同除以zn，式（2）可演變?yōu)椋?/p>

其中：式（3）即系統(tǒng)動態(tài)過程的多輸入多輸出誤差模型。

2.2 球磨機磨礦系統(tǒng)數(shù)學(xué)模型辨識

球磨機磨礦系統(tǒng)可以描述為如圖3所示的多變量系統(tǒng)，其中介質(zhì)充填率、料球比和磨礦濃度為輸入?yún)?shù)，功率消耗為輸出參數(shù)，由于系統(tǒng)同時受多個輸入信號作用，為了使系統(tǒng)具有良好的可辨識性，各輸入信號之間需滿足一定的關(guān)系，即各輸入信號需要相互獨立[14]。

圖3 球磨機多變量控制模型Fig.3 Multivariable control model of ball mill

其中：式中的輸入項u（k）、輸出項y（k）即可代表實際測量數(shù)據(jù)，引入時間平移算子q，則q-1y（k）=y（k-1），另外，為了準(zhǔn)確的描述對象的純滯后特性，可在式（4）中加入時間滯后算子q-d，d為滯后的拍數(shù)[16]，記q-du（k）=u（k-d），由式（4）可得球磨機磨礦系統(tǒng)的數(shù)學(xué)模型表示為：系統(tǒng)辨識工具箱函數(shù)ARX（）可以用于多變量系統(tǒng)的辨識[15]，本文選擇基于最小二乘法估計的

ARX參數(shù)辨識方法，以ARX模型對球磨機磨礦系統(tǒng)的數(shù)學(xué)模型進行描述如下：

離散后處理可得該模型的離散形式為：

其中na、njb、njk分別表示輸出量、輸入量和時間延時階次，式中將各種各樣的誤差均歸并到噪聲項e（k）中去。

3 球磨機磨礦系統(tǒng)辨識過程

系統(tǒng)辨識過程步驟為：辨識試驗設(shè)計、數(shù)據(jù)預(yù)處理、選擇模型結(jié)構(gòu)及參數(shù)估計、模型檢驗。

3.1 辨識試驗設(shè)計

根據(jù)圖3球磨機多變量控制模型中多輸入單輸出參數(shù)，需要分別辨識出介質(zhì)充填率u1、料球比u2和磨礦濃度u3作為系統(tǒng)輸入量，功率y作為輸出量，建立球磨機系統(tǒng)數(shù)學(xué)模型。前期通過對球磨機磨礦現(xiàn)場考察獲得經(jīng)驗知識，根據(jù)試驗條件，輸入信號采用階躍信號，試驗過程的介質(zhì)填充率、料球比均設(shè)置5個階躍水平，磨礦濃度設(shè)置4個階躍水平，如表1所示。采樣周期設(shè)置為600 s，共采集了100組試驗輸入輸出數(shù)據(jù)。

表1 輸入量階躍水平分布 %Tab.1 distribution table of Input step level

3.2 試驗數(shù)據(jù)處理

系統(tǒng)辨識實際測量直接得到的數(shù)據(jù)是按時間序排列，須進行數(shù)據(jù)預(yù)處理，主要包括重采樣、消除數(shù)據(jù)趨勢項和對測量數(shù)據(jù)進行濾波[17]。本次試驗數(shù)據(jù)因事先按周期采集，可省去數(shù)據(jù)重采樣。數(shù)據(jù)預(yù)處理前，先將試驗獲得的輸入、輸出數(shù)據(jù)導(dǎo)入到MATLAB工作空間中，然后將其轉(zhuǎn)換成為iddata對象形式。數(shù)據(jù)導(dǎo)入程序如下：

數(shù)據(jù)導(dǎo)入后，輸入輸出序列的均值不為零且隨時間變化，此處采用MATLAB系統(tǒng)辨識工具箱中的dtrend函數(shù)消除數(shù)據(jù)的趨勢項，把測量的數(shù)據(jù)變成一系列在零值上下波動的數(shù)據(jù)序列。采用idfilt函數(shù)利用Butterworth濾波器對數(shù)據(jù)進行濾波，若數(shù)據(jù)頻率緩慢變化則沒有高頻成分。運行后的結(jié)果如圖4所示。

3.3 辨識結(jié)果與模型檢驗

使用ARX()函數(shù)可以直接辨識出系統(tǒng)模型參數(shù)，該模型辨識的階次為：na=[5]，nb=[5 5 5]，nk=[10 10 3]，辨識結(jié)果如式（7）。

圖4 輸入輸出數(shù)據(jù)曲線及其消除趨勢項后曲線Fig.4 Input/output data curve and the graph after eliminating trends item（a）—功率消耗消除趨勢項曲線；（b）—介質(zhì)填充率消除趨勢項曲線；（c）—料球比消除趨勢項曲線；（d）—磨礦濃度消除趨勢項曲線

獲得辨識模型后，可以通過tf（）函數(shù)來提取系統(tǒng)的傳遞函數(shù)矩陣。最后，將球磨機磨礦試驗中負(fù)荷參數(shù)量作為辨識模型的輸入量，由辨識模型計算出的輸出量與磨礦試驗實際獲得的輸出數(shù)據(jù)進行對比，比較結(jié)果如圖5所示。

觀察圖5可知，將試驗輸入數(shù)據(jù)導(dǎo)入辨識模型得到的模型輸出數(shù)據(jù)曲線與實際輸出數(shù)據(jù)曲線很接近，基本重現(xiàn)實際輸出變化規(guī)律，其中出現(xiàn)的誤差很可能是受其他不可測干擾對其輸出功率影響。辨識結(jié)果表明，該球磨機系統(tǒng)辨識模型是合適的，滿足要求。利用辨識方法所得的數(shù)學(xué)模型能較好地反映球磨機磨礦系統(tǒng)的實際動態(tài)特性。另外，由于試驗過程操作或測量誤差的存在，也可能存在部分模型輸出與實際輸出不相符的情況，通過多次試驗辨識出模型可使模型輸出盡可能接近實際輸出。

圖5 試驗輸出功率與辨識模型輸出功率曲線Fig.5 Curve graphs of experimental output power and identification model output power

4 結(jié)語

球磨機磨礦效果好壞直接影響鎢礦選礦效率，從影響球磨機磨礦效率及功率消耗因素出發(fā)，針對球磨機多變量磨礦系統(tǒng)的控制運用系統(tǒng)辨識方法進行辨識和研究。

根據(jù)前期工作知球磨機負(fù)荷參數(shù)是影響球磨機磨礦效率和功率消耗主要因素之一，通過將球磨機負(fù)荷參數(shù)（介質(zhì)填充率、料球比、磨礦濃度）與功率消耗分別作為辨識的輸入量與輸出量，設(shè)計了辨識試驗方案，建立出球磨機磨礦系統(tǒng)內(nèi)部負(fù)荷與功率消耗之間的數(shù)學(xué)模型，最后將辨識模型輸出功率曲線

與實測功率曲線對比，證明了將系統(tǒng)辨識方法應(yīng)用到球磨機磨礦系統(tǒng)數(shù)學(xué)模型預(yù)測中，具有較好的建模效果，有利于控制球磨機負(fù)荷參數(shù)工作在最佳工況條件下，指導(dǎo)鎢礦回收工藝過程中使磨礦效率最高，功率消耗最低，從而推進鎢礦山強化資源的高效利用。

[1] 盛繼福.我國鎢礦資源形勢[J].中國鎢業(yè)，1992，（9）：5-9. SHENG Ji-fu.The situation of tungsten resources in China[J]. China Tungsten Industry，1992，7（9）：5-9.

[2] 余良暉，馬茁卉，周海東.我國鎢礦資源開發(fā)利用現(xiàn)狀與發(fā)展建議[J].中國鎢業(yè)，2013，28（4）：6-9. YU Liang-hui，MA Zhuo-hui，ZHOU Hai-dong.Prensent situation and development suggestion of exploiting and utilizing China's tungsten resource[J].China Tungsten Industry，2013，（4）：6-9.

[3] 常富強，段德華，宋 龑.生產(chǎn)中提高球磨機磨礦效率的方法[J].現(xiàn)代礦業(yè)，2011，（3）：81-84. CHANGFu-qiang，DUANde-hua，SONGYan.Methodforimproving grinding efficiency of ball mill in production[J].Modern Mining，2011，（3）：81-84.

[4] 李英霞，王國生.某白鎢礦回收工藝的試驗研究[J].中國鎢業(yè)，2012，27（4）：23-25. LI Ying-xia，WANG Guo-sheng.Recovery process experiment of a scheelite[J].China Tungsten Industry，2012，27（4）：23-25.

[5] 河南豫暉礦山.球磨機的鎢礦磨礦作業(yè)的好壞有關(guān)的三個因素[EB/OL].[2014-06-25]http：//www.yzw.cc/space/blog/283532.html.

[6] 王澤紅，陳炳辰.球磨機負(fù)荷檢測的現(xiàn)狀與發(fā)展趨勢[J].中國粉體技術(shù)，2001，7（1）：19-23. WANG Ze-hong，CHEN Bing-chen.The current situation and trend of development of load testing[J].China Powder Science and Technology，2001，7（1）：19-23.

[7] 楊麗榮，蔡改貧.基于神經(jīng)網(wǎng)絡(luò)與SOPC的球磨機負(fù)荷檢測系統(tǒng)[J].礦業(yè)研究與開發(fā)，2007，27（3）：46-48. YANG Li-rong，CAI Gai-pin.Design of ball mill load mesuring system based on the neural network and SOPC[J].Mining Research and Development，2007，27（3）：46-48.

[8] 王 恒，賈民平，陳左亮.基于軟測量的球磨機優(yōu)化控制方法[J].控制工程，2011，18（5）：762-766. WANG Heng，JIA Min-ping，CHEN Zuo-liang.Method of optimize control for ball mill based on soft sensor[J].Control Engineering of China，2011，18（5）：762-766.

[9] 王建民，賀曉巧，曹艷忙，等.基于專家控制的磨機負(fù)荷控制算法[J].河北聯(lián)合大學(xué)學(xué)報：自然科學(xué)版，2015，（1）：77-82. WANG Jian-min，HE Xiao-qiao，CAO Yan-mang，et al.Control algorithm for the mill load based on expert control[J].Journal of Hebei UniteUniversity：NaturalScienceEdition，2015，（1）：77-82.

[10]SBARBARO D，VILLAR R D.Advanced control and supervision of mineral processing plants[M].London：Springer-Verlag，2010.

[11]楊志剛，張 杰，李艷姣.磨礦過程綜合自動化系統(tǒng)[J].河北聯(lián)合大學(xué)學(xué)報：自然科學(xué)版，2014，36（1）：66-70. YANG Zhi-gang，ZHANG Jie，LI Yan-jiao.Integrated automation system of grinding process [J].Journal of Hebei Polytechnic University：Social Science Edition，2014，36（1）：66-70.

[12]丁 鋒.系統(tǒng)辨識（2）：系統(tǒng)描述的基本模型[J].南京信息工程大學(xué)學(xué)報：自然科學(xué)版，2011，3（2）：97-117. DING Feng.System identification.Part B：Basic model for system description[J].Journal of Nanjing University of Information Science and Technology：Natural Science Edition，2011，3（2）：97-117.

[13]蕭德云.系統(tǒng)辨識理論及應(yīng)用[M].北京：清華大學(xué)出版社，2014.

[14]吳 昊，余岳峰，徐星星.多輸入多輸出熱工系統(tǒng)的辨識與建模研究[J].動力工程學(xué)報，2010，30（3）：196-200. WU Hao，YU Yue-feng，XU Xing-xing.Study on identification and modeling of multi-input and multi-output thermal system[J]. Journal of Chinese Society of Power Engineering，2010，30（3）：196-200.

[15]薛定宇.控制系統(tǒng)計算機輔導(dǎo)設(shè)計—MATLAB語言與應(yīng)用[M].北京：清華大學(xué)出版社，2012.

[16] 張偉偉.工業(yè)鍋爐燃燒系統(tǒng)辨識與建模研究[D].上海：上海交通大學(xué)，2007. ZHANG Wei-wei.Study on the identification and modeling of combustion system of industrial boiler[D].Shanghai：Shanghai JiaoTong University，2007.

[17]蔡改貧，熊 洋，姜志宏，等.基于預(yù)磨機入料排料粒級分布建模和參數(shù)辨識研究[J].礦山機械，2014，（9）：71-74. CAI Gai-ping，XIONG Yang，JIANG Zhi-hong，et al.Research on modeling and parameter identification based on pre-mill feedings and discharge particle size distribution[J].Mining and Processing Equipment，2014，（9）：71-74.

Multivariable System Identification Modeling of Ball Mill Grinding Systems

CAI Gai-pin,XU Qin,XIONG Yang,YANG Li-rong
(School of Mechanical and Electrical Engineering,Jiangxi University of Science and Technology,Ganzhou 341000,Jiangxi,China)

The discrete form of the ball mill system model was derived by applying the multivariate ARX model identification method based on the least squares estimation algorithm.This paper established the mill dynamic mathematical model of the mill grinding system by applying internal load parameters of ball mill,including media fill rates,ratio of material and ball and grinding concentration)and the power consumption in the grinding process.The use of MATLAB system identification toolbox modeling method and implementation of the identification process in detail were described as well.The comparison of the model output curve and the actual output curve demonstrated that the high accuracy model provided favorable theoretical basis for the realization of high-efficiency and lowenergy process control.

lease squares;multivariate ARX model;identification modeling;mill load;power consumption;system identification

TD453；TP273+.5

A

10.3969/j.issn.1009－0622.2015.05.013

2015-08-10

江西省自然科學(xué)基金項目（20151BAB206033）

蔡改貧（1964-），男，江西贛州人，博士，教授，主要從事近凈成形新技術(shù)、物料高效破碎技術(shù)研究與裝備開發(fā)。