本发明涉及一种基于选择性融合多通道机械信号频谱多特征子集的磨机负荷参数预测方法。
背景技术:
基于球磨机系统机械振动和振声信号的磨机负荷参数预测(mlpf)已经成为实现磨矿过程运行优化与反馈控制的技术瓶颈[1]。研究表明不同的磨机负荷参数(料球比(mbvr)、磨矿浓度(pd)和充填率(cvr))与不同通道机械信号频谱特征间存在着不同的映射关系[2]。对磨机的机械振动和振声信号直接进行快速傅里叶变换(fft)是一种常用分析方法,所获得频谱称之为单尺度频谱[2]。但是,fft在本质上并不适于分析具有非平稳和多组分特性的机械信号。经验模态分解(emd)能够将原始机械信号分解为系列具有不同时间尺度的内禀模态函数(imf)[3]。理论上,这些不同时间尺度的imf子信号具有不同的物理含义。基于imf的fft被应用于轴承故障诊断和磨机筒体振动分析[4,5],这些基于imf的频谱被称之为多尺度频谱。通过噪声辅助技术,集成emd(eemd)能够克服emd所固定的模态混合问题[6],并已被用于对磨机筒体振动和振声信号进行预处理[7]。另外的一种方法,希尔伯特振动分解(hvd)能够将原始信号分解为具有不同能量尺度分布的子信号[8]。通常,这些不同的时频方法可从多个不同的分解角度获得有价值的子信号。
通常,在较高的频率分辨率下,单尺度或多尺度频谱数据通常包含数以千计的输入特征。许多维数约简方法(包括特性选择和特征提取)可对频谱数据进行预处理[9,10]。这些方法可以有效避免预测模型的过拟合、抵制噪声影响和增强模型预测性能。基于互信息(mi)的特征选择方法比其他特征选择方法更易于理解[11,12],但同主成分分析(pca)方法类似,必须单独选择另外的学习模型对这些选择的特征进行建模[13]。潜结构映射或偏最小二乘(pls)模型能够通过提取与建模数据输入输出均相关的潜在特征(lvs)建立线性模型。针对过多的输入变量导致pls模型的过拟合问题[14],遗传算法-pls(ga-pls)用于高维数据和磨机筒体振动单尺度频谱的特征子集选择[15,16];但该方法需要较大的计算消耗才能克服ga的随机性获得优化的特征子集。基于l1正则化的最小二乘方法能够建立基于少量输入特征的稀疏模型[17],即最小绝对收缩与选择算子(lasso)算法[18]。在高维的微阵列和基因数据上的分析表明,lasso在作为理想的特征子集选择器方面具有重要潜力[19,20]。
为了有效利用单尺度筒体振动频谱中不同类特征子集的优点,核pca(kpca)、mi、频谱聚类和自适应ga(aga)用于构建基于最小二乘支撑向量机的单模型[21]。针对基于emd、eemd和hvd等时频分析算法得到的多尺度频谱,采用选择性的信息融合机制提高了mlpf模型的预测性能和简化了模型学习参数的选择过程[22]。文献[23]提出基于lasso选择多尺度频谱特征子集以构建mlpf模型,但是该建模方法未考虑如何选择合适的惩罚参数,也未考虑同一频谱中的不同特征子集的贡献率。显然,仅仅是采用多尺度频谱的单一特征子集难以构建有效的mlpf模型。最近,从选择性融合多工况样本和多源信息的视角出发,文献[24]提出了基于多层sen的mlpf模型,表明基于单尺度频谱特征子集的模型在预测性能上要高于基于多尺度频谱的模型。此外,在球磨机系统的不同位置(通道)采集的机械信号,如磨机筒体和轴承座振动,磨机筒体近表面和研磨区域下方振声等,对构建可靠mlpf模型的贡献显然不同。因此,非常有必要选择更有价值的通道信号及其频谱特征。
技术实现要素:
本发明提供一种基于选择性融合多通道机械信号频谱多特征子集的磨机负荷参数预测方法,基于最小绝对收缩与选择算子(lasso)和选择性集成(sen)学习机制的mlpf算法,首先,采用快速傅里叶变换(fft)得到多通道机械振动和振声信号的单尺度频谱。接着,针对每个通道的单尺度频谱,采用基于候选惩罚参数的集成构造策略得到面向候选单尺度特征子集构建的候选lasso子模型,结合sen学习机制构建全部通道的候选单通道sen-lasso模型。最后,再次的采用sen学习机制对候选单通道sen-lasso模型进行选择和合并,得到多通道机械信号及其频谱特征子集同时优化选择的mlpf模型。
附图说明
图1gci和基于运行专家的人工操作示意图;
图2(a)ch1的单尺度频谱;
图2(b)ch2的单尺度频谱;
图2(c)ch3的单尺度频谱;
图2(d)ch4的单尺度频谱;
图2(e)ch5的单尺度频谱;
图2(f)ch6的单尺度频谱;
图2(g)ch7的单尺度频谱;
图2(h)ch8的单尺度频谱;
图3(a).ch1通道为pd构建的lasso子模型选择的频谱特征所对应的回归系数值;
图3(b).ch2通道为pd构建的lasso子模型选择的频谱特征所对应的回归系数值;
图3(c).ch3通道为pd构建的lasso子模型选择的频谱特征所对应的回归系数值;
图3(d).ch4通道为pd构建的lasso子模型选择的频谱特征所对应的回归系数值;
图3(e).ch5通道为pd构建的lasso子模型选择的频谱特征所对应的回归系数值;
图3(f).ch6通道为pd构建的lasso子模型选择的频谱特征所对应的回归系数值;
图3(g).ch7通道为pd构建的lasso子模型选择的频谱特征所对应的回归系数值;
图3(h).ch8通道为pd构建的lasso子模型选择的频谱特征所对应的回归系数值;
图4单通道sen-lasso模型的预测误差;
图5(a).基于验证数据的ch1-4的单通道sen-lasso模型的预测曲线;
图5(b).基于验证数据的ch5-6的单通道sen-lasso模型的预测曲线;
图5(c).基于测试数据的ch1-4的单通道sen-lasso模型的预测曲线;
图5(d).基于测试数据的ch1-4的单通道sen-lasso模型的预测曲线;
图6(a)mbvr模型的预测曲线;
图6(b)pd模型的预测曲;
图6(c)cvr模型的预测曲线;
图7为本发明的流程图。
具体实施方式
磨矿过程描述和基于运行专家的磨机负荷估计
磨矿过程磨矿回路(gc)的目的是得到解离的有价值矿物,其控制策略通常是:在某一固定的磨矿粒度下尽可能的最大化磨矿生产率(gpr)[25]。国内的铁矿选矿过程通常采用两段式的闭环gc,最大化gpr通常是通过设定循环负荷的最优值得到,而循环负荷通常是由一段gc(gci)的球磨机负荷决定的。因此,球磨机应该维持在最优负荷状态以保证磨矿过程的运行优化。gci和基于运行专家的人工操作如图1所示。
如图1所示,原矿通过振动给料机给到运输皮带,然后输送至湿式预选机,进入一段磨矿回路(gci);在gci内,湿式预选机通过磁力选择有用矿石和抛尾矿,然后混合来自一段旋流器的沉砂以及周期性添加的钢球,通过给矿器进入一段球磨机;球磨机依靠筒体旋转带动钢球对矿石进行冲击破碎,形成矿浆;矿浆依靠自身的流动性排出磨机,进入一段泵池,与泵池内的新加水混合后的矿浆被泵入一段旋流器;一段旋流器将矿浆分为粒度较细的溢流和较粗的沉砂,后者进入一段球磨机再磨,构成一段球磨的闭路循环;前者进入一次磁选机进行选别,选别的溢流为尾矿,沉砂则进入二段磨矿回路(gcii)。gcii的研磨过程是与gci相同的闭路循环过程。
由图1可知,磨机的入口负荷包括:新给矿、新给水、水力旋流器的沉砂及周期性添加的钢球,磨机的出口负荷包括矿浆及磨碎的钢球。实际生产中,为保证矿浆具有合适的粘度或后续处理过程的需要,还需在磨机入口添加化学药剂。由上述分析可知,即使gci安装了所有需要的检测仪表,检测了泵池液位、磨机电流、磨机转速、泵池内矿浆的浓度、水力旋流器的流量和压力等过程变量,磨机负荷仍然难以依据物料流所蕴含的金属平衡确定。主要原因如下:(1)水力旋流器的溢流和沉砂中的球负荷,以及磨机衬板的磨损与腐蚀量难以检测;(2)用于测量泵池内矿浆的浓度、水力旋流器的流量和压力等过程变量的仪表精度难以保证;(3)球磨机的新给矿具有随机变化的特性,如原料的粒度、品位、硬度、表面的微小裂纹及特性分布等;(4)磨机内部的物料和钢球的粒度及其分布、矿浆的流变特性等也难以检测。
综上,磨矿过程的复杂机理、研磨工况的频繁波动,以及磨机内部众多复杂多变、难以检难的研磨参数等因素使我们难以采用解析方法建立磨机负荷机理模型。
在实际工业过程中,为保证磨机运行在最佳工况,领域专家需融合不同来源、不同类型的信息判断磨机负荷状态,这个过程如图1的右上方所示:(1)获取信息:包括在线仪表检测的数据(磨机振声信号、磨机前后轴承的振动信号、水力旋流器驱动电机的电流信号、泵池内的矿浆浓度和水力旋流器的溢流浓度)、离线化验数据(原矿、尾矿、精矿的硬度和粒度等)及现场技术人员听到或看到的异常工况信息(现场技术人员观测的磨机出口状态、磨矿过程关键岗位的操作和巡视人员提供的其它相关信息);(2)结合历史经验和多源信息,判断磨机是否工作在最优工况即判断磨机负荷状态;(3)调整新给矿、磨机给水、钢球、泵池给水量。由于专家经验的差别和其有限的精力,难以长期保持磨机运行在最佳工况[26],导致能耗和钢耗的增加。
实际上,基于专家经验的磨机负荷估计过程是一个选择有价值信息进行有效融合的过程。磨机内部的负荷参数(mbvr,pd和cvr)能够有效表征磨机负荷,且与磨机研磨产生的多通道机械振动/振声信号的频谱密切相关。因此,可以采用更多来源(通道)的机械信号构建mlpf模型,并通过模拟运行专家的选择性信息融合机制对多通道机械信号及其多频谱特征进行优化融合。
如图7所示,本发明提供一种基于选择性融合多通道机械信号频谱多特征子集的磨机负荷参数预测方法,包括以下步骤:
步骤1、对多通道机械振动和振声的时域信号采用快速傅里叶变换(fft)得到多通道单尺度频谱
将所采集的全部机械振动/振声信号的通道数量记为j,并进一步将全部机械信号表示为
其中,
假定在时域内针对每个通道的机械信号均有k个长度为n的建模样本变换至频域,全部j个通道的单尺度频谱可表示为:
其中p=1,…,p,p是单尺度频谱特征的数量;l=1,…,k,k是建模样本的数量;j=1,…,j,j是机械振动和振声信号通道的数量;
步骤2、针对每个通道的单尺度频谱,采用基于候选惩罚参数的集成构造策略得到面向候选单尺度特征子集构建的候选lasso子模型,结合sen学习机制构建全部通道的候选单通道sen-lasso模型
采用高分辨率的单尺度频谱数据建模的最佳方案是所构建的mlpf模型能够同时完成维数约简和模型学习。此处,采用lasso算法实现单尺度频谱的特征选择和mlpf回归模型的建立,并预先设定一组数量为ican的lasso模型的候选惩罚参数
以第jth个通道的单尺度频谱为例,基于第ith个候选惩罚参数λji构建lasso子模型。针对本发明,lasso算法通过求解如下l1范数回归问题以获得回归向量
其中,
通过求解如上公式,将基于λji和
进而,全部的候选lasso子模型可以表示为
研究表明,不同的频谱特征对磨机负荷参数的贡献不同。为提高mlpf模型的泛化性能,选择和合并部分lasso子模型以构建sen-lasso模型是必要的。
假定有isel(sen-lasso集成模型的集成尺寸)个-lasso子模型被选择,采用自适应加权融合(awf)算法合并。针对第lth个建模样本,第jth个通道的sen-lasso模型
其中
其中,σji是第ith个被选择的-lasso子模型预测输出
在确定了上述集成模型的集成尺寸和加权算法后,构建sen-lasso模型所需要的lasso子模型的选择可通过求解如下的优化问题得到:
其中,rmsre(均方根平均相对误差)用于评估集成尺寸为isel的sen模型的泛化性能;iset是用于控制单通道sen-lasso模型结构复杂度的预设定集成尺寸阈值。公式(8)的求解详见文献[1]。
通过上述过程,可得到第jth个通道的sen-lasso模型
通过重复上述过程j次,可获得全部通道机械信号的sen-lasso模型
步骤3、采用sen学习机制对候选单通道sen-lasso模型进行选择和合并,得到多通道机械信号及其频谱特征子集同时优化选择的mlpf模型
由于不同的单通道sen-lasso模型中包含着针对磨机负荷参数的不同的有价值信息,有必要选择性融合这些sen-lasso以构建最终的mlpf模型。与构建基于不同惩罚参数的单通道sen-lasso模型类似,最终mlpf模型的构建是再次采用文献[1]中的sen机制。为保证本发明的完整性,将所最终构建的mlpf模型拟解决的优化问题描述如下:
其中,rmsre用于评估集成尺寸为jsel的最终mlpf模型的泛化性能;jsel是被选择的单通道sen-lasso模型的个数,即最终mlpf模型的集成尺寸;fsen(·)和
最终构建的mlpf模型的预测输出
其中,最终mlpf模型中第jth个单通道sen-lasso模型的加权系数采用下式计算:
其中,σj是第jth个单通道sen-lasso模型的预测输出
实验研究
本实验在直径为602mm和长度为715mm的小型实验磨机上进行,其中磨机筒体的旋转速度为42r/min。实验中,总共有8个同道的机械信号以51200hz的频率进行采集,传感器类型和位置为:(1)固定在磨机筒体表面的2个振动加速度传感器(记为sv1和sv2);(2)与磨机筒体表面相距2mm的2个声传感器(记为sa1和sa2);(3)位于磨机轴承座左侧和右侧的3个振动加速度传感器(记为av1,av2和av3);(4)位于磨机研磨区域下方10mm的1个声传感器(记为milla)。这些机械信号依次被标记为通道1到8,即ch1-ch8。
本次实验是在固定钢球和水负荷、逐渐增加矿石负荷的情况进行的,总共进行了139次实验,即采集了139个样本,其中4/5的样本用做建模的训练和验证数据集,其余的用于模型测试。
单尺度频谱结果
首先,对时域机械振动/振声信号进行滤波处理;然后,采用fft技术将磨机稳定运行过程中每完整周期的数据转换至频域,得到每个通道的多个旋转周期的单尺度频谱;最后,将这些稳定旋转周期的频谱数据进行平均,获得最终的单尺度频谱。总计139次实验的全部8个通道的单尺度频谱如图2(a)至图2(h)所示。
图2(a)至图2(h)表明,不同通道的机械信号频谱具有不同的频率分布特性。对于单一通道的机械信号频谱,不同的频谱特征子集也包含不同的有价值信息。为构建有效的mlpf模型,有必要进行多通道机械信号及其频谱特征子集的选择。
单通道单尺度sen-lasso模型结果
本发明中,预设定单通道sen-lasso模型的集成尺寸为10,并将所采用的候选惩罚参数集合设定为{0.1,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5,10](这些候选惩罚参数的对应编号分别为1-21),即lasso子模型采用不同惩罚参数值时对应不同的回归系数值,进而也对应不同的特征子集。针对磨矿浓度(pd)的众多候选lasso子模型,对不同通道机械信号所选择的频谱特征所对应的回归系数值如图3(a)至图3(h)所示。
图3(a)至图3(h)表明,对不同的单尺度频谱选择不同的特征子集构建单通道sen-lasso模型是必要的。
为了更为简洁的描述,下文中将精度最佳的lasso子模型和加权全部21个lasso子模型的集成模型分别标记为best-sub和enall-sub。
面向pd模型,不同类型的单通道频谱模型的验证数据误差和所包含的子模型惩罚参数编号的统计结果如表1所示。
表1.面向pd模型的不同类型单通道频谱模型的统计结果(rmsre和λji编号)
表1表明:基于通道ch6(av2)的三类mlpf模型均具有佳的预测性能,单通道sen-lasso模型通常具有最小的rmsre;候选集合中前3个较小的惩罚参数(0.1,0.5和1.0,对应的编号为1,2和3),被选择频率很高,表明预设定惩罚参数集合有待于进一步优化确定。
全部单通道sen-lasso模型针对训练、验证和测试数据集的预测误差(rmsres)如图4所示。
基于验证和测试数据的全部单通道sen-lasso模型的预测曲线如图5(a)至图5(d)所示。
上述结果表明,不同通道的机械信号频谱对mlpf模型的贡献不同,有必要对其进行选择性的融合以提高mlpf模型的泛化性能。
选择性融合单通道sen-lasso模型的结果
通过对8个单通道sen-lasso模型进行sen学习机制可获得最终的mlpf模型。基于测试数据,面向pd模型的统计结果(rmsres和通道编号)如表2所示.
表2基于测试数据面向pd模型的统计结果(rmsres和通道编号)
表2表明:采用通道(3,8,1,4,6)所构建的pd模型具有最小的rmsre(0.006431),其集成尺寸为5;同时,当集成尺寸为4时,得到的预测误差为0.006703,其相对于最佳的单通道sen-lasso模型和简单融合全部通道sen-lasso的模型在预测性能提高了1倍多;由表2的备注可知,这4个被选择的通道来源于筒体振动、筒体振声、轴承振动和研磨区域下方振声信号。由此可见,与pd相关的有价值信息蕴含在不同的互补机械信号通道中。
针对cvr,基于本发明所提方法的最小rmsre是0.008145,所选择通道数是ch1和ch6;基于ch6所构建的最佳单通道sen-lasso模型的预测误差是0.008753。针对mbvr,基于本发明所提方法的最小rmsres是0.03322,选择的通道数是ch1和ch6;基于ch6所构建的最佳单通道sen-lasso模型的预测误差是0.0087530.03481。
这些结果表明,在三个磨机负荷参数mbvr、pd和cvr中,pd对机械振动/振声信号具有最佳的灵敏度,该结论与之前的研究相符合。
不同磨机负荷参数的预测曲线如图6(a)至图6(c)所示。
上述结果表明,本发明所提方法对于选择性融合多通道机械振动/振声信号及其单尺度频谱特征子集是有效的。
构建可靠的基于机械振动/振声信号的磨机负荷参数预测(mlpf)软测量模型对实现磨矿过程的运行优化与反馈控制至关重要。本发明提出了基于lasso和sen学习机制的面向多通道机械信号频谱多特征子集的单尺度频谱的建模方法。单通道sen-lasso模型能够选择单尺度频谱特征子集构建lasso子模型并加权融合。最终的mlpf模型通过选择和合并这些单通道的sen-lasso模型得到。基于实验球磨机多通道频谱数据的仿真结果验证了方法的有效性。
参考文献
[1]tang,j.,chai,t.y.,yu,w.,andzhao,l.j.(2013).modelingloadparametersofballmillingrindingprocessbasedonselectiveensemblemultisensorinformation.ieeetransactionsonautomationscience&engineering,10(3),726-740.
[2]tang,j.,tian,f.q.,jia,m.y.,andli,d.l.rotatingmechanicaldeviceloadsoftmeasuringbasedonfrequencyspectraldata-driven.beijing:nationaldefenseindustrypress,2015.(inchinese).
[3]huang,n.e.,shen,z.,andlong,s.r.(1998).theempiricalmodedecompositionandthehilbertspectrumfornon-linearandnonstationarytimeseriesanalysis,proc.royalsoc.londona,454,903-995.
[4]rai,v.k.,andmohanty,a.r.(2007).bearingfaultdiagnosisusingfftofintrinsicmodefunctionsinhilbert-huangtransform.mechanicalsystemsandsignalprocessing,21,2607-2165.
[5]tang,j.,zhao,l.j.,yue,h.,yu,w.,andchai,t.y.(2011).vibrationanalysisbasedonempiricalmodedecompositionandpartialleastsquares,procediaengineering,2011,16,646-652.
[6]wu,z.h.,andhuang,n.e.(2009).ensembleempiricalmodedecompositionforhighfrequencyecgnoisereduction.advancesinadaptivedataanalysis,55(4),193-201.
[7]tang,j.,liu,z.,wu,y.j.,andzhao,l.j.(2014).modelingdifficult-to-measureprocessparametersbasedonintrinsicmodefunctionsfrequencyspectralfeaturesofmechanicalvibrationandacousticalsignals.advancedmaterialsresearch,989-994,3671-3674.
[8]feldman,m.(2006).time-varyingvibrationdecompositionandanalysisbasedonhilberttransform,journalofsoundandvibration,295(3-5),518-530.
[9]jiménez-rodríguez,l.o.,arzuaga-cruz,e.,andvélez-reyes,m.(2007).unsupervisedlinearfeature-extractionmethodsandtheireffectsintheclassificationofhigh-dimensionaldata.ieeetransactionongeoscienceandremotesensing,45(2),469-483.
[10]wang,l.(2008).featureselectionwithkernelclassseparability,ieeetransactionsonpatternanalysisandmachineintelligence,30(9),534-1546.
[11]liu,h.,sun,j.,liu,l.,andzhang,h.(2009).featureselectionwithdynamicmutualinformation.patternrecognition,42(7),1330-1339.
[12]hao,r.h.,yang,z.f.,andwen,x.x.(2009).anefficientgeneselectionalgorithmbasedonmutualinformation.neurocomputing,72(4-6),991-999.
[13]tang,j.,chai,t.y.,cong,q.m.,yuan,b.c.,zhao,l.j.,liu,z.,andyu,w.(2014).softsensorapproachformodelingmillloadparametersbasedonemdandselectiveensemblelearningalgorithm.actaautomaticasinica,40(9),1853-1866.
[14]hawkins,d.m.(2004).theproblemofoverfitting.journalofchemicalinformationandcomputersciences,44(1),1-12.
[15]leardi,r.,seasholtz,m.b.,andpell,r.j.(2002).variableselectionformultivariatecalibrationusingageneticalgorithm:predictionofadditiveconcentrationsinpolymerfilmsfromfouriertransform-infraredspectraldata,analyticachimicaacta,461(2),189-200.
[16]tang,j.,zhao,l.j.,zhou,j.w.,yue,h.,andchai,t.y.(2010).experimentalanalysisofwetmillloadbasedonvibrationsignalsoflaboratory-scaleballmillshell[j].mineralsengineering,23(9),720-730
[17]hastie,t.,tibshirani,r.,andfriedman,j.h.(2001).theelementsofstatisticallearning.springer,august.
[18]tibshirani,r.j.(1996).regressionshrinkageandselectionviathelasso,journaloftheroyalstatisticalsociety,58,267-288.
[19]shevade,s.k.,andkeerthi,s.s.(2003)asimpleandeifficientalgorithmforgeneselectionusingsparselogisticregression.bioinformatics,19(17),2246-2253.
[20]roth,v.(2014).thegeneralizedlasso.ieeetransactiononneuralnetworks,15(1),16-28.
[21]tang,j.,chai,t.y.,yu,w.,andzhao,l.j.(2012).featureextractionandselectionbasedonvibrationspectrumwithapplicationtoestimatetheloadparametersofballmillingrindingprocess,controlengineeringpractice,20(10),991-1004.
[22]tang,j.,chai,t.y.,cong,q.m.,liu,z.,andyu,w.(2015).modelingmillloadparametersbasedonselectivefusionofmulti-scaleshellvibrationfrequencyspectra,controltheoryandapplications,32(12),1582-1591.
[23]tang,j.,qiao,j.f.,liu,z.,chai,t.y.,andyu,w.(2017).modelingmillloadparameterbasedonlassousingmulti-scalehighdimensionalfrequencyspectradata.ieeeinternationalconferenceoninformationandautomation(pp.909-914).ieee,newyork.
[24]tang,j.,qiao,j.f,wu,z.w.,chai,t.y,zhang,j.,andyu,w.(2018).vibrationandacousticfrequencyspectraforindustrialprocessmodelingusingselectivefusionmulti-conditionsamplesandmulti-sourcefeatures.mechanicalsystems&signalprocessing,99,142-168.
[25]lo,y.c.,oblad,a.e.,andherbst,j.a.(1996).costreductioningrindingplantsthroughprocessoptimizationandcontrol.minerals&metallurgicalprocessing,13(1),19-21.
[26]su,z.g.,andwang,p.h.(2009).improvedadaptiveevidentialk-nnruleanditsapplicationformonitoringlevelofcoalpowderfillinginballmill.journalofprocesscontrol,19(10),1751-1762。