基于靈敏度分析與增強捕食-食餌優化的重介質選煤過程動態模型

王光輝，彭 勇，代 偉，董 良，馬小平

(1.中國礦業大學 化工學院，江蘇 徐州 221116; 2.中國礦業大學 地下空間智能控制教育部工程研究中心，江蘇 徐州 221116; 3.中國礦業大學 信息與控制工程學院，江蘇 徐州 221116)

重介質選煤是常用的煤炭洗選工藝之一，其是將密度介于精煤和矸石之間的重介質懸浮液和原煤混合后送入重介質旋流器，在離心力和重力的作用下實現精煤和矸石的分離。重介質旋流器因分選精度高、分選密度調節范圍寬、適應能力強等特點，在選煤工藝中廣泛使用[1-2]。

先進控制是提升重介質選煤效率和穩定性的重要手段之一，其實現往往需要依托過程動態數學模型，因此重介質選煤過程模型的研究一直是礦物分選領域的重要研究內容之一。現有的研究成果主要研究集中在靜態模型、流體力學模型、數據驅動模型、質量守恒模型等。如文獻[3]通過多項方程式擬合實驗數據法的方法，建立了一種將配分因子作為相對密度和礦石粒度函數的重介質選煤靜態模型。文獻[4]采用同樣的建模與策略構建了一種能夠表示分選效率的靜態模型。然而靜態模型僅能用于工藝設計與優化，難以基于此設計閉環控制控制系統。文獻[5]通過使用Fluent等軟件求解Navier-Stokes方程建立了重介質旋流器的計算機流體動力學(Computational Fluid Dynamics,CFD)模型。文獻[6-8]將CFD法與離散元法相結合，建立了煤粒在重介質旋流器中的動態模型，但CFD模型過于復雜而難以用于控制器設計[9]。當前，數據驅動建模方法由于較少依賴先驗知識，難以機理建?；蛉狈C理知識的復雜工業過程中被廣泛使用。文獻[10]采用機器學習模型，建立了基于變塊寬增量隨機配置網絡的精煤灰分模型。文獻[11]針對重介質選煤過程動態時變與精煤灰分數據稀缺感知的問題，提出了基于在線自適應半監督學習的精煤灰分模型。但上述基于機器學習方法的數據驅動模型是在模型輸入特征空間映射到復雜高維非線性特征空間基礎上，再進行回歸計算，模型結構相對復雜，無法直接用于控制器設計。

質量平衡模型是一種描述原料投入、產量與損耗之間關系的機理模型，文獻[12]基于礦漿中各物料的質量平衡關系，建立了重介質選煤動態模型，被用于設計了重介質選煤過程模型預測控制器與自適應控制器[13-14]。但所建立的質量平衡模型對不確定的模型參數，主要是依靠經驗或實驗試湊來獲得，難以精確描述實際動態過程，給實際控制系統的設計帶來難題。因此，重介質選煤過程控制系統性能的提高迫切需要進一步開展模型參數的優化方法研究，以提高模型質量。

受這一問題驅動，筆者致力于采用捕食-食餌優化(Prey-Predator Optimization,PPO)算法[15]對重介質選煤過程動態模型參數進行優化。PPO是一種新型元啟發式優化算法，相較于傳統算法，具有收斂速度快、對最優值敏感等優點。在PPO算法框架中，搜索步長是決定對最優值搜索程度的重要因素，傳統方法將Best Prey以恒定的搜索步長對最優值進行搜索。此時，如果搜索步長過小將導致收斂過慢，過大則可能導致在臨近最優解時被跳過，從而陷入局部最優。文獻[16]根據判斷Prey與最優值之間的距離來選擇步長，可快速調整2者之間距離，達到快速收斂的目的，但依然需要依托在人為設定步長的基礎上;文獻[17]利用遞減函數來實現搜索步長隨迭代的自適應變小，但搜索步長只能機械性變小，無法自適應調整算法的搜索能力。

為此，筆者應用增強學習思想，提出增強PPO(Reinforcement PPO,RPPO)算法，主要思想是，搜索個體利用神經網絡將自身的狀態信息映射到動作集合，通過動作來控制搜索步長的增減以及網絡權重的更新，從而使優化模型中搜索個體通過不斷學習自身歷史信息來選擇模型參數，建立數據驅動的優化模型參數自適應調整策略，最終提高優化性能。此外，充分考慮機理模型參數眾多，難以對模型參數直接進行優化的難題，采用了低差異的Sobol’序列與基于方差的Sobol’參數靈敏分析相結合的方法[18]，分析了在建模過程中產生的不確定性參數對精煤中灰分含量的影響，確定各參數對模型輸出的重要性，進而對這些“重要的”模型采用RPPO算法進行尋優。通過與其他傳統啟發式優化算法進行對比，分析了所提算法的性能;且采用實際數據驗證了所提模型的準確性。

1 重介質選煤過程與動態模型

1.1 重介質選煤過程描述

典型的重介質選煤過程如圖1所示[19]，經過篩選加工的原煤和由合格介質桶傳輸出的密度適中的重介質懸浮液共同加入混料桶中。經充分混合后，礦漿由礦漿泵打入到重介質旋流器，在離心力和重力的作用下，密度小于重介質懸浮液的煤礦在重介質旋流器上方聚集，相反，密度大于重介質懸浮液的礦物雜質由重介質旋流器底流口流出。重介質旋流器的溢流和底流經脫水脫介處理后，分別形成精煤和尾礦，送下一道工序;脫水脫介篩將殘余的介質溶液送入磁選機中進行回收;從磁選機回收的介質與從高濃介質桶中的高濃度介質、稀釋水以一定比例共同加入合格介質桶中，經過混合后進行再次工藝循環。

圖1 重介質選煤過程工藝流程Fig.1 Flow diagram of classical dense medium coal separation process

1.2 重介質選煤動態模型

重介質選煤過程主要分為礦漿混合、重介質分選以及重介質回收3個動態過程。

1.2.1混合過程模型

原煤經過破碎等工藝后，與來自合格介質桶傳輸的重介質懸浮液在混料桶中進行充分混合。由質量平衡定理可得混合過程模型為

(1)

(2)

(3)

Wi=Qmbρmb

(4)

1.2.2重介質分選過程模型

從混料桶中流出的混合礦漿送入重介質旋流器，在離心力和重力的作用下，煤與礦物雜質分別在底流和溢流礦漿中被排出，根據質量平衡定理可建立如下的重介質分選過程模型:

ρm)xi.ash]-xu,c[Wi-Qoρo-KoVo(ρo-ρm)xi,C]-

Ku,cVuρu(ρc-ρm)(xi,c-xu,c)}

(5)

ρm)xi.ash]-xu,m[Wi-Qoρo-KoVo(ρo-ρm)xi,C]-

Ku,mVuρu(ρo,m-ρm)(xi,m-xu,m)}

(6)

ρm)xi.ash]-xu,c[Wi-Qoρo-KoVo(ρo-ρm)xi,C]-

Ko,cVoρo(ρm-ρc)(xi,c-xo,c)}

(7)

ρm)xi.ash]-xu,m[Wi-Qoρo-KoVo(ρo-ρm)xi,C]-

Ko,mVuρu(ρm-ρu,m)(xi,m-xo,m)}

(8)

(9)

(10)

(11)

(12)

式中，Vc為重介質旋流器中礦漿體積，假設其不變。

1.2.3重介質回收過程模型

重介質選煤過程需要通過對旋流器產物進行脫介、脫水來回收磁性介質，回收的介質與高濃介質在合格介質桶中經稀釋水按照期望的介質濃度進行調節。為保證模型的精確性，需考慮介質回收中的損耗。因此假設介質從重介質旋流器到磁選機的回收率為β，介質從磁選機到合格介質桶的回收率為γ，磁選機出口的重介質密度維持在ρrm，從而可建立的重介質回收過程模型為

磁選機可回收到的介質質量mm:

(13)

磁選機可回收到的介質體積流量Qrm:

(14)

合格介質桶內部動態模型:

(15)

(16)

(17)

其中，Qmm，ρmm分別為高濃介質的體積流量及密度;Vcor為重介質溶液在合格介質桶中的體積;Qw為稀釋水的體積流量;ΔP，Cw和lw分別為閥門系數、閥門開度和閥門壓差。

在上述動態模型中，待確定的模型參數有α，Ku，Ko，Ku,ash，Ko,ash等13個，實際過程中難以直接測量。由模型可以看出，對上述參數進行優化，是一個復雜的非線性優化問題，啟發式優化算法是解決這一問題的有效方法。然而，過多的模型參數會大大增加優化算法的復雜性和難度。實際上，模型輸出可能僅對一部分參數敏感，其它參數對模型品質不起決定性作用。因此，可通過參數靈敏度分析來確定模型參數的“重要性”，然后只對“重要的”模型參數進行尋優，大大簡化計算負荷，并有利于最優值的求解[20]。

2 基于參數靈敏度與增強捕食-食餌的模型參數優化

2.1 參數靈敏度分析

參數靈敏度分析研究了機理模型輸出的變化如何歸因于其輸入參數的變化。從參數變化范圍的角度可將靈敏度分析方法分成局部靈敏度分析和全局靈敏度分析方法。局部靈敏度只能在參數范圍內進行小范圍靈敏度分析，難以應用在非線性模型中;而全局靈敏度分析則考慮了參數在整個取值范圍內的變化。同時，考慮到重介質選煤過程動態模型具有強非線性、不確定性參數多等特點，一般的局部靈敏度分析方法難以用于對此模型的分析，而Sobol’靈敏分析方法作為一種單次可處理多輸入變量的全局靈敏度分析方法，與重介質選煤過程動態模型特點十分契合。因此，文本采用基于方差的Sobol’參數靈敏度方法對上述重介質選煤過程動態模型參數進行靈敏度分析。

首先在各個模型參數允許的取值范圍內進行采樣，并求得模型輸出，將模型輸出的方差歸因于某些參數或某些參數的集合，對模型輸出的影響大小定義為模型靈敏度。

Sobol’方法定義只含各參數或各階參數集合的模型輸出偏方差與模型輸出總方差的比值為各參數以及參數之間的靈敏度:

(18)

其中，1≤i1<…

(19)

其中，S下標的個數為靈敏度的階數，如Si為一階靈敏度，用于估計單獨參數i對模型輸出的影響;Sij為二階靈敏度，表示參數i，j之間的組合對模型輸出的影響;S12…d為d階靈敏度，用于計算不同參數組合對模型輸出影響。另定義參數i在f(x)中的全階靈敏度為STi=1-V-i/V，用于計算含有參數i的靈敏度總和，其中，V-i為不包含參數i的其余參數偏方差。

一般情況下,上述方差通過積分求得[21]，但由于重介質選煤動態模型復雜難以直接解析求得積分，因此本文基于蒙特卡羅法估計總方差和偏方差，具體步驟為:

(20)

(21)

(22)

式中，f(A)j下標j為采樣矩陣的第j行。

式(22)用于計算某參數的全階靈敏度，以此作為模型參數的靈敏度的最終判斷依據，可確定各參數對模型的“重要性”。基于此，將相對“不重要”的模型參數固定在經驗值，繼而采用元啟發式優化算法對“重要的”模型參數進行尋優。

2.2 增強捕食-食餌優化算法

本文提出的增強捕食-食餌優化(Reinforcement Prey-Predator Optimization,RPPO)算法是一種基于捕食者(Predator)與食餌(Prey)之間相互作用的元啟發式算法。其首先需要定義m個可行解，{?1,?2,…,?i,…,?m}，并將每個可行解賦予生存值(Sv)，Sv可通過計算目標函數得到，且與目標函數性能成正比例關系。通過Sv大小做出以下分類:最小可行解稱作Predator，最大可行解稱作Best Prey，其余稱作Prey。在完成各角色分配后，Predator會追捕Prey，而Prey自然會遠離Predator，而且會傾向于跟隨生存值更高且位置較近的Prey運動，如果不滿足跟隨概率，則只會在保證遠離Predator的情況下做隨機方向運動。不同于Prey，Best Prey則不需要擔憂Predator的追捕，只考慮尋找更好的躲避處，即專注于對最優解的搜索(Exploitation)即可。采用增強學習方法不斷尋找最優的搜索步長，從而保證了算法的收斂速度和準確性。

2.2.1更新可行解位置

(1)Best Prey

(23)

其中，ξ1∈[0,1]分布;dl為一個隨機向量;‖·‖為歐氏距離;Best Prey沿著dl方向在自身周圍采用較小的搜索步長(λmin)對最優值進行搜索，防止跳過最優值。

(2)Prey

(24)

其中，ξ2,ξ3,ξ4∈[0,1]且服從均勻分布;Pf為Prey向其他Prey靠攏的概率，如果Pf≤ξ4，Prey只會沿隨機方向dr運動來躲避Predator;反之，Prey還會以di方向跟隨其他Prey;λmax為Prey的探索步長;根據文獻[15]得

(25)

在式(25)中，di由其他Prey的方向及其生存值共同決定。其中，u為比例參數，對Sv(?j)在di中的貢獻進行縮放。

在式(24)中，dr為Prey躲避Predator的方向，Prey應沿著離Predator較遠的方向移動，因此，需先做以下判斷:

(26)

dr1=‖?p.d-(?i+dr)‖

(27)

dr2=‖?p.d-(?i-dr)‖

(28)

其中，dr1，dr2為Prey分別沿著±dr方向移動后與Predator的歐式距離。

(3)Predator

(29)

算法在運行過程中會涉及到探索、搜索2種步長，探索步長主要是在迭代前期對可行解空間的探索，擴大種群個體的可搜索空間;而搜索步長是在迭代后期針對可行解現有空間的搜索，搜索步長過大不利于對現有空間充分搜索，過小則需要冗余的迭代次數才能尋得最優解，因此適當的搜索步長對算法最終尋得的最優解起著至關重要的作用。為避免因搜索步長固定不變而導致算法尋優能力不足的問題，本文采用增強學習的思想，用種群中的個體代表智能體，通過讓其對自身歷史經驗的學習來選擇動作，依據個體采取的動作來判斷搜索步長的變化，整個算法流程如圖2所示。

2.2.2基于增強學習的自適應搜索步長

選擇Prey作為增強學習個體，通過設定個體的動作、狀態及獎懲機制，建立增強學習模型，實現搜索步長的自適應調整。增強學習的思想是個體接收當前狀態，并通過神經網絡將狀態映射到相應的動作，環境根據所選動作更新狀態，與此同時，通過獎懲機制計算所采取的動作相應的獎懲，并據此更新神經網絡，逐漸獲取與環境相一致的狀態-動作映射。本文所定義的動作、狀態及獎懲機制如下:

(1)動作。個體設定3種動作:步長加速、步長減速、步長不變，即

(30)

(2)狀態。個體的T個歷史狀態組成狀態集合Statet， Statet會隨著選擇動作而產生相應的變化，并由式(31)來更新。

(31)

(32)

(3)獎懲機制。每次迭代中，當個體與最優解的距離變小時，應當給予獎勵;反之，應當收到懲罰:

(33)

(4)狀態-動作映射。選用單隱含層神經網絡來映射狀態與動作之間的關系。以式(31)得到的狀態集合作為神經網絡的輸入，即輸入節點有T個;將動作集合中的3種動作作為3個輸出節點，并通過式(34)選擇將要執行的動作:

y=min|1-yi| 1≤i≤3

(34)

其中，y為與1距離最小的輸出節點;yi為神經網絡第i個輸出節點，當i取1,2,3時，分別對應步長加速、減速與不變3個動作，即

(35)

本文將隱含層和輸出層的激活函數設為Sigmoid函數，隱含層到輸出層的權重更新公式:

(36)

其中，η為神經網絡的學習率;g為隱含層的輸出;di為輸出節點i的期望輸出，設定當個體受到獎賞，即Feedbackt=+1時，di=1，否則為0;yi為輸出節點i的實際輸出。同樣，可得到輸入層到隱含層權重公式:

(37)

RPPO算法每次迭代都會先確定搜索步長，然后通過所描述方程進行更新Best Prey，普通Prey及Predator的位置。具體描述如下:

輸入:可行解(x1,x2,…,xm)，最大迭代次數Niter，狀態向量T，步長變化因子σ，神經網絡學習率η，目標函數Fit，PPO算法參數(λmax，λmin，Pf)

(1)將可行解作為Fit的輸入，得到x的生存值Sv;

(2)初始化神經網絡;

(3)fort=1:Niter

按Sv從小到大對可行解排序;

利用式(29)更新Predator位置;

fori=2:m-1

根據式(24)更新Prey位置;

i=i+1;

End

根據式(23)更新Best Prey位置;

當前所受獎懲Feedbackt由式(33)計算;

利用式(36),(37)更新神經網絡權重;

根據式(31)更新狀態集合;

下一代步長可由式(30)得到;

t=t+1;

End

(4)選擇Sv的可行解作為最終結果;

輸出:Best Prey的位置及生存值

3 實 驗

精煤中灰分是決定煤炭的質量主要因素，因此，本文主要研究對象為分選動態過程模型中精煤灰分模型，并以其為目標函數。在給煤量不斷變化情況下，本文將探究通過所提方法對模型參數進行優化后，精煤灰分模型是否能夠較好跟隨實測精煤灰分。

3.1 靈敏度分析

原煤給煤量Wore的變化情況如圖3所示;同時，根據各參數的實際物理含義[12]，并通過工藝知識或實驗確定其取值范圍，見表1;而選煤模型運行過程中所需設備參數值見表2。

表1 模型參數取值范圍Table 1 Range of each model parameters

表2 設備參數取值Table 2 Values of equipment parameters

圖3 原煤的給煤量變化曲線Fig.3 Curve of raw coal feed rate

在此范圍內，利用Sobol’序列進行采樣，以采樣所得參數集為輸入，運行重介質分選過程模型。通過Sobol’方法對多參數變化下的各模型進行分析，得到各模型參數的一階靈敏度及全階靈敏度，如圖4，5所示。

圖4，5利用箱型圖表示各參數在完整工藝過程中的靈敏度變化范圍，其中，菱形為所有時刻靈敏度的均值。從圖4，5可知，雖各參數的一階靈敏度和全階靈敏度有所差別，但不影響參數“重要性”的判斷，因此文中將不考慮各參數之間的高階靈敏度。同時，本文認定靈敏度均值大于0.05且變化范圍較大的參數為“重要”參數，需要將其進行進一步的參數優化。由圖5可知重介質分選過程模型中只有α，Ku,ash，Ko,ash符合上述要求，說明在設計重介質回收過程模型參數時，應考慮對這些參數進行優化。而對于輸出結果影響較小的參數，將通過工程經驗確定其近似值。

圖4 模型參數的一階靈敏度Fig.4 First-order indices of each model parameters

圖5 模型參數的全階靈敏度Fig.5 Total-effect indices of each model parameters

3.2 參數優化

為使所建立模型能更好的描述選煤過程，本文以每一時刻模型仿真輸出與實際過程24 h測量數據之間的差距為目標函數，整個工業過程中每一時刻的差距之和越小說明優化所得的參數越合適，即所有差距和的倒數越大越好:

(38)

其中，Fit為目標函數;t為模型運行過程中的某時刻;fs為所建立的模型輸出;fm(t)為實際測量數據;τ為模型運行時間。

為驗證RPPO在重介質選煤模型參數優化過程中的有效性，本文將其與PPO、遺傳算法(Genetic Algorithm,GA)[23]以及粒子群算法(Particle Swarm Optimization,PSO)[24]在相同初始條件下運行，各算法對目標函數的尋優過程如圖6所示。本文根據實驗經驗得出各參數取值，見表3。

表3 各算法參數設置Table 3 Parameters of each algorithm

圖6 目標函數隨迭代次數的變化情況Fig.6 Variation of the objective function with iterations

為了避免初始條件的隨機性對實驗結果的影響，本文采取了種群同一初始位置的方法，使得4種算法初始條件相同，并基于此初始位置進行30次重復試驗。圖6為某一次對比實驗結果，其他結果類似。由圖6可看出，在開始迭代的初期，RPPO和PSO可快速將目標函數收斂在理想最優值附近，這說明RPPO和PSO在迭代初期可以實現很好的探索;結束迭代過程后，以RPPO算法的結果最為突出，說明了RPPO中Best Prey的后期搜索能力相較于其他算法有所提高。在整個尋優過程中，個體的搜索步長變化如圖7所示，個體可根據自身的歷史信息來判斷所處搜索空間的形勢，進而選擇不同的動作，保證了對現處地形的充分搜索，可避免選擇次優解。

圖7 搜索步長變化Fig.7 Variation diagram of exploitation step lengths

表4給出了整個24 h的實驗統計結果。從表4可看出，RPPO和PSO尋找最優值能力相當;在平均值、最差值及標準差指標中，RPPO和PPO表現較好且接近，表明這2種算法的穩定性都較高。說明相較于其他3種算法，RPPO兼顧尋優能力及穩定性，因此，本文選用RPPO來進行重介質選煤的參數進行優化。

表4 4種算法統計結果Table 4 Statistical results of the four algorithms

為使所建立模型貼合實際工藝過程，本文基于已建立的重介質分選過程模型，進行了各優化算法對比尋優實驗，將所得的模型參數代入重介質分選過程模型公式得到整個工藝過程的精煤灰分，如圖8所示。其中，RPPO所優化的模型參數為(1.389 7，6.270 0×10-4，7.934 4×10-4)。

圖8 模型仿真結果與實際測量數據對比Fig.8 Comparison curves between the simulation results and actual measurement data

為了更清楚的表明所提算法的性能，圖8給出了具有代表性的30 min內的對比實驗結果。從工業過程來看，由RPPO所得灰分曲線圖與實際采樣數據最為接近，PPO次之;計算得各算法(按圖例從上至下)與實測數據的均方誤差分別為0.009 3，0.004 6，0.006 6，0.004 2。圖9為不同算法優化所得模型的誤差概率密度函數(Probability Density Function，PDF)對比圖，圖9中，各PDF曲線是以誤差的均值為中心值所得，顯然，RPPO與PPO算法的PDF曲線中心點更接近0，且RPPO的PDF曲線更高，與其他方法相比標準差更小;圖10為以0為中心值的各算法PDF曲線，從圖10可以看出，RPPO算法的PDF曲線最高且窄，效果明顯優于其他算法。綜上，可得RPPO相較于傳統PPO及其他傳統算法能夠更精確的估計重介質選煤模型參數，便于基于模型的進一步工作展開。

圖9 不同算法建模誤差PDF對比(以均值為中心)Fig.9 PDF comparison curves of modeling errors based on mean values by different algorithms

4 結 論

(1)針對一個典型由礦漿混合、重介質分選和重介質回收3部分組成的重介質選煤過程，建立基于質量平衡重介質選煤過程動態數學模型。

(2)基于所建立模型，首先利用Sobol’參數靈敏度分析方法分析出“重要的”模型參數，進而采用自適應步長的RPPO算法進行模型參數優化，從而提高了重介質選煤動態模型精度，使其更接近真實的工業過程。

(3)從實驗結果可以看出，本文所提的方法相較于傳統PPO等算法可自適應控制步長，搜索能力有所增強，進而能夠避免選擇次優解且穩定性高;應用在重介質選煤模型參數優化中可充分發揮RPPO的優勢，尋找到最佳的模型參數，使得模型輸出更貼合實際工業采樣數據。

(4)本文所提的靈敏度分析與增強捕食-食餌優化相結合的模型參數優化策略，不限于本文所研究的典型重介質選煤過程模型，可推廣用于解決煤炭領域其他過程或裝備的建模問題。