曹衛(wèi)東，歐陽(yáng)騁，余 陽(yáng)，李力泓，梁新利，姜博嚴(yán)

(1.河海大學(xué) 物聯(lián)網(wǎng)工程學(xué)院，江蘇 常州 213022;2.中國(guó)船舶工業(yè)綜合技術(shù)經(jīng)濟(jì)研究院，北京 100081)

0 引言

支持向量回歸(Support Vector Regression,SVR)在預(yù)測(cè)領(lǐng)域中取得了可喜成果[1-5]。在車削領(lǐng)域，王興盛等[6]利用最小二乘支持向量機(jī)預(yù)測(cè)鏡片精密車削表面粗糙度，并對(duì)工藝參數(shù)進(jìn)行了分析。結(jié)果表明，刀具圓弧半徑和每圈進(jìn)給量對(duì)表面粗糙度影響較顯著。ALAJMI等[7]使用SVR根據(jù)車削工藝參數(shù)確定刀具磨損，結(jié)果表明使用SVR估計(jì)工藝參數(shù)是可行的。在銑削領(lǐng)域，李堯等[8]使用模糊支持向量回歸進(jìn)行顫振診斷，預(yù)測(cè)準(zhǔn)確率為97.3%。CHARALAMPOUS[9]使用SVR進(jìn)行銑削切削力估算。在鉆削領(lǐng)域，GU等[10]提出了基于SVR的磨損預(yù)測(cè)模型，預(yù)測(cè)結(jié)果與實(shí)驗(yàn)數(shù)據(jù)基本一致。然而，懲罰因子、內(nèi)核類型等超參數(shù)都會(huì)影響到預(yù)測(cè)能力，上述研究大多采用人工經(jīng)驗(yàn)設(shè)置超參數(shù)，需要耗費(fèi)較長(zhǎng)時(shí)間尋找合適的參數(shù)。因此，一些學(xué)者將亨利氣體溶解度算法(Henry Gas Solubility Optimization,HGSO)[11-12]、樽海鞘算法(Salp Swarm Algorithm,SSA)[13-14]、蜻蜓算法[15-16]、蟻獅算法[17-18]、粒子群算法(Particle Swarm Optimization,PSO)[19,20]等智能算法應(yīng)用于SVR的參數(shù)微調(diào)，提高預(yù)測(cè)精度。但是，由于搜索個(gè)體和迭代次數(shù)的增加，運(yùn)行時(shí)間也成倍增加。尤其是在大數(shù)據(jù)預(yù)測(cè)方面，巨大的時(shí)間消耗困擾著許多學(xué)者。雖然使用K-means聚類分析進(jìn)行數(shù)據(jù)預(yù)處理[21-23]可以大大縮短運(yùn)行時(shí)間，但也會(huì)因?yàn)槿狈颖静糠中畔ⅲ档皖A(yù)測(cè)精度。

可以發(fā)現(xiàn)，與智能算法結(jié)合后的SVR具有很好的預(yù)測(cè)能力，預(yù)測(cè)精度較高，但需要更多的運(yùn)行時(shí)間，時(shí)間復(fù)雜度高；K-means聚類分析可以大大降低運(yùn)行時(shí)間，但它需要割舍掉部分樣本數(shù)據(jù)，損害了整體的預(yù)測(cè)性能。受到該問題的啟發(fā)，本文主要對(duì)以下兩個(gè)方面進(jìn)行研究：①K-means聚類對(duì)數(shù)據(jù)集進(jìn)行分析后，尋找一種智能算法優(yōu)化SVR的超參數(shù)，能讓預(yù)測(cè)精度、穩(wěn)定性和時(shí)間復(fù)雜度綜合最優(yōu)，建立一種超參數(shù)自適應(yīng)SVR方法;②使用不同維度的數(shù)據(jù)集以及滾齒工藝參數(shù)預(yù)測(cè)數(shù)據(jù)集對(duì)方法進(jìn)行測(cè)試，顯示其優(yōu)異的綜合性能。

針對(duì)上述問題，本文旨在尋找一種智能算法,以平衡預(yù)測(cè)能力與運(yùn)行時(shí)間。哈里斯鷹算法(Harris Hawks Optimization, HHO)[24]是一種新穎的元啟發(fā)式智能算法，具有非常優(yōu)異的性能，但有落入局部最優(yōu)解的缺陷。本文引入邏輯混沌映射，降低落入局部最優(yōu)解的可能，以形成混沌哈里斯鷹算法(Chaos HHO, CHHO)，將該方法用于K-means聚類與SVR混合方法中，形成超參數(shù)自適應(yīng)SVR方法，即Kmeans-CHHO-SVR。最后,在測(cè)試集上進(jìn)行仿真實(shí)驗(yàn)，驗(yàn)證了算法的優(yōu)越性。隨后，將Kmeans-CHHO-SVR應(yīng)用到滾齒工藝參數(shù)預(yù)測(cè)領(lǐng)域，驗(yàn)證其工程應(yīng)用效果。

1 超參數(shù)自適應(yīng)支持向量回歸方法

Kmeans-CHHO-SVR的主要框架如下：①使用訓(xùn)練集和K-means聚類分析獲取聚類中心CCs;②根據(jù)SVR超參數(shù)集群和CCs獲得驗(yàn)證集的預(yù)測(cè)結(jié)果;③使用CHHO超參數(shù)集群，并記錄最佳SVR超參數(shù)。重復(fù)步驟②和步驟③，直至達(dá)到終止條件?；谧罴裇VR超參數(shù)，得到測(cè)試集的預(yù)測(cè)結(jié)果，預(yù)測(cè)過程如圖1所示。SVR超參數(shù)集群X由SVR參數(shù)個(gè)體Xi組成。X={X1、X2、…、Xn}，其中n是個(gè)體數(shù)。

Kmeans-CHHO-SVR具體步驟如下：

步驟1設(shè)置最大迭代次數(shù)MaxIter，迭代次數(shù)t，t=1。

步驟2將數(shù)據(jù)歸一化到-1～1之間，將其分為訓(xùn)練集集、驗(yàn)證集和測(cè)試集。

步驟3基于訓(xùn)練集，使用K-means聚類分析來獲取聚類中心CCs。

步驟4輸入SVR參數(shù)的下界LB和上界UB，根據(jù)LB和UB隨機(jī)初始化X。

步驟5基于X和CCs，使用SVR獲取訓(xùn)練模型。

步驟6根據(jù)訓(xùn)練模型和驗(yàn)證集獲取驗(yàn)證集的均方誤差MSEV。

步驟7找出第t代MSEV的最小值和最佳SVR參數(shù)Xopt。

步驟8若t>MaxIter，轉(zhuǎn)步驟11；否則，轉(zhuǎn)步驟9。

步驟9使用Xopt和CHHO的核心步驟來更新X。

步驟10t=t+1，轉(zhuǎn)步驟5。

步驟11根據(jù)Xopt使用SVR獲取測(cè)試集的均方誤差MSET。

1.1 混沌哈里斯鷹優(yōu)化器(CHHO)

哈里斯鷹算法(HHO)是HEIDARI等提出的新型啟發(fā)式算法[24]，靈感來自哈里斯鷹的合作行為和追逐風(fēng)格。哈里斯鷹的捕獵過程描述如下：

(1)探索階段 使用式(1)模擬哈里斯鷹探索獵物的行為。

(1)

(2)開發(fā)階段 該階段哈里斯鷹根據(jù)兔子逃跑方式和環(huán)境影響在軟圍攻、硬圍攻、漸進(jìn)式快速俯沖的軟圍攻和漸進(jìn)式快速俯沖的硬圍攻4種追逐方式中選擇一種去突襲兔子。

1)軟圍攻：

X(t+1)=△X(t)-E|J·Xopt(t)-X(t)|，

(2)

△X(t)=Xopt(t)-X(t)。

(3)

2)硬圍攻：

X(t+1)=Xopt(t)-E|△X(t)|。

(4)

3)漸進(jìn)式快速俯沖的軟圍攻：

Y=Xopt(t)-E|J·Xopt(t)-X(t)|。

(5)

哈里斯鷹根據(jù)萊維飛行進(jìn)行俯沖：

Z=Y+S×LF(D)。

(6)

式中：D表示問題的維數(shù)，S表示1×D大小的隨機(jī)向量，LF是萊維飛行函數(shù)，可以使用式(7)得到

(7)

其中μ和θ表示0～1之間的隨機(jī)值。

鷹在軟圍攻階段的位置更新可以總結(jié)為：

(8)

式中Y和Z可以使用式(5)和式(6)求得。

4)漸進(jìn)式快速俯沖的硬圍攻：

(9)

Y=Xopt(t)-E|J·Xopt(t)-Xm(t)|，

(10)

Z=Y+S×LF(D)。

(11)

CHHO核心操作的偽代碼如圖2所示，混沌值rC決定了算法進(jìn)入了哪個(gè)階段。rC是混沌隨機(jī)值[25]，使用式(12)初始化rC：

科文學(xué)院新辦公樓主要是行政人員,網(wǎng)絡(luò)用戶主要使用DB數(shù)據(jù)庫(kù)訪問、E-mail郵件和WWW訪問,上網(wǎng)人數(shù)在60人左右;辦公樓(老校區(qū))主要是教學(xué)人員,網(wǎng)絡(luò)用戶主要使用FTP文件傳輸、E-mail郵件和WWW訪問,上網(wǎng)人數(shù)在30人左右;新教學(xué)樓北區(qū)2～4層、教學(xué)樓部分樓層是學(xué)生機(jī)房,網(wǎng)絡(luò)用戶主要使用FTP文件傳輸、DB查詢和WWW訪問,上網(wǎng)學(xué)生數(shù)峰值在640人左右．網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)模型如圖2所示,其中,NIC為科文學(xué)院網(wǎng)絡(luò)中心,Student為學(xué)生子網(wǎng),Teach為教工子網(wǎng),Office為辦公子網(wǎng).Student子網(wǎng)拓?fù)浣Y(jié)構(gòu)如圖3所示.

rC(t+1)=4rC(t)(1-rC(t))。

(12)

1.2 K-means聚類分析

K-means聚類應(yīng)用十分廣泛，且十分有效[26]。K-means聚類算法的步驟如下。

步驟1設(shè)置初始聚類中心個(gè)數(shù)cNum并初始化聚類中心。

步驟2計(jì)算樣本到每個(gè)聚類中心的距離,并將其分配給與其距離最小的聚類中心所對(duì)應(yīng)的類中。

步驟3根據(jù)類中成員重新計(jì)算每類的聚類中心。

步驟4重復(fù)步驟2和步驟3，直到滿足沒有對(duì)象被重新分配給不同的聚類。

本文選擇歐式距離作為距離度量，通過隨機(jī)操作對(duì)聚類中心進(jìn)行初始化。參數(shù)cNum的設(shè)置較為關(guān)鍵，實(shí)驗(yàn)1具體描述了它的選擇過程和結(jié)果。

(1)RBF核函數(shù)：

KT(·)=exp(-γ|u-v|2)。

(2)Sigmod核函數(shù)：

KT(·)=tanh(γ(uTv)+coef)。

2 實(shí)例驗(yàn)證

驗(yàn)證基礎(chǔ)條件包括:①計(jì)算機(jī)，CPU為Intel Xeon W-2123，內(nèi)存為16 GB，WIN 10系統(tǒng);②MATLAB R2018a;③LIBSVM;④大數(shù)據(jù)集，從UCI機(jī)器學(xué)習(xí)共享庫(kù)[30]中獲取。5個(gè)大數(shù)據(jù)集如表1所示。

表1 數(shù)據(jù)集的基本信息

為了探究Kmeans-CHHO-SVR的最佳設(shè)置和算法性能，實(shí)驗(yàn)由以下3個(gè)部分組成：首先進(jìn)行了cNum設(shè)置實(shí)驗(yàn)，探究cNum的最佳取值;實(shí)驗(yàn)2比較了Kmeans-CHHO-SVR(CA0), Kmeans-SVR(CA1), Kmeans-HGSO-SVR(CA2), Kmeans-SSA-SVR(CA3), Kmeans-DA-SVR(CA4), Kmeans-ALO-SVR(CA5), Kmeans-PSO-SVR(CA6)和 Kmeans-HHO-SVR(CA7)在各個(gè)數(shù)據(jù)集下的預(yù)測(cè)精度和穩(wěn)定性，分析了Kmeans-CHHO-SVR的性能；第3個(gè)實(shí)驗(yàn)研究了每個(gè)方法的時(shí)間復(fù)雜度。表2給出了各方法的初始參數(shù)設(shè)置。

表2 參數(shù)設(shè)置

2.1 實(shí)驗(yàn)1：cNum設(shè)置

本實(shí)驗(yàn)研究了不同數(shù)據(jù)集下，Kmeans-CHHO-SVR的cNum設(shè)置對(duì)預(yù)測(cè)精度和運(yùn)行時(shí)間的影響。圖3給出了cNum取[20,200]下的MSET和運(yùn)行時(shí)間結(jié)果，共計(jì)20次實(shí)驗(yàn)。可以看出，隨著cNum的增加，MSET值在數(shù)據(jù)集D1～D4上相對(duì)穩(wěn)定，在數(shù)據(jù)集D5上劇烈波動(dòng)；運(yùn)行時(shí)間在D1和D3上呈總體上升趨勢(shì)，在另外3個(gè)數(shù)據(jù)集上出現(xiàn)了波動(dòng)和上升，并最終有所下降。

為了選擇合適的cNum，定義CMT=β·NMSET+δ·NTC，用來衡量cNum選取的優(yōu)劣，CMT越小越好，其中β=0.8，δ=0.2，NMSET和NTC是歸一化后的MSET值和運(yùn)行時(shí)間。

表3給出了Kmeans-CHHO-SVR在5個(gè)數(shù)據(jù)集下的CMT值，OptC表示cNum的合適值。舉例說明：在數(shù)據(jù)集D1下，cNum選擇60最合適，對(duì)應(yīng)的CMT為0.023。

表3 CMT和OptC結(jié)果(CMT最小值用下劃線標(biāo)示)

2.2 實(shí)驗(yàn)2：預(yù)測(cè)精度和穩(wěn)定性

為了檢驗(yàn)Kmeans-CHHO-SVR的預(yù)測(cè)精度和穩(wěn)定性，將其分別與CA1, CA2, CA3, CA4, CA5, CA6相比，驗(yàn)證其在數(shù)據(jù)集D1-D5上的性能。初始參數(shù)如表2和表3所示。實(shí)驗(yàn)選用十倍交叉驗(yàn)證，訓(xùn)練集、驗(yàn)證集、測(cè)試集的比例為98∶1∶1。根據(jù)第1章的方法步驟和數(shù)據(jù)集獲取MSEV和MSET，重復(fù)該過程10次。

表4給出了預(yù)測(cè)精度(MSEV和MSET)的平均值(Avg)和標(biāo)準(zhǔn)差(Std)，圖4和圖5給出了MSEV和MSET的箱型圖。結(jié)果顯示，Kmeans-CHHO-SVR在數(shù)據(jù)集D1-D5上的MSEV和D2-D5上的MSET上表現(xiàn)最佳(聯(lián)合最佳)。如圖6所示，考慮驗(yàn)證集的標(biāo)準(zhǔn)差StdV，Kmeans-CHHO-SVR在D1、D3和D4上的輸出最穩(wěn)定；考慮測(cè)試集的標(biāo)準(zhǔn)差StdT，Kmeans-CHHO-SVR在D3和D4上的輸出最穩(wěn)定。

表4 MSEV和MSET結(jié)果

2.3 實(shí)驗(yàn)3：時(shí)間復(fù)雜度

時(shí)間復(fù)雜度描述了某個(gè)方法的運(yùn)行時(shí)間。通常使用O進(jìn)行表述。測(cè)算依據(jù)：本文方法屬于組合方法，它的時(shí)間復(fù)雜度隨著樣本數(shù)目、樣本屬性數(shù)目、采用哪種核函數(shù)等參數(shù)而變化，采用O表述不能直觀地與其他方法進(jìn)行比較。故本文方法運(yùn)行10次耗費(fèi)的平均時(shí)間作為時(shí)間復(fù)雜度的測(cè)算依據(jù)。本實(shí)驗(yàn)通過比較相同的條件下，在數(shù)據(jù)集D1～D5上預(yù)測(cè)消耗的平均時(shí)間衡量復(fù)雜度，消耗時(shí)間越少，時(shí)間復(fù)雜度越低，測(cè)試過程與實(shí)驗(yàn)1幾乎相同。定義CCRk評(píng)價(jià)時(shí)間復(fù)雜度如下：

(13)

(14)

式中:Crate表示歸一化后的時(shí)間消耗，Cold表示時(shí)間消耗，CMin(CMax)表示最小(最大)時(shí)間消耗。nd=10。

表5給出了實(shí)驗(yàn)結(jié)果，可以看出，Kmeans-CHHO-SVR的表現(xiàn)不理想，但與其他算法的差距不大。

表5 各個(gè)算法的平均消耗時(shí)間

3 實(shí)驗(yàn)討論

表6匯總了所有結(jié)果，其中“W”表示獲勝，“T”并列，“L”表示失敗。定義Rank評(píng)估算法綜合表現(xiàn)，Rank越低越好。

(15)

表6顯示，Kmeans-CHHO-SVR模型的Rank值為0.226，綜合表現(xiàn)最佳。可以看到，由于CHHO的存在，Kmeans-CHHO-SVR模型在許多數(shù)據(jù)集上預(yù)測(cè)精度最佳，且在D3和D4上擁有最好的穩(wěn)定性，雖然混沌搜索增加了時(shí)間消耗，但模型取得了更佳的性能。因此，相比其他方法，Kmeans-CHHO-SVR效率更高，在處理大數(shù)據(jù)方面具有較大優(yōu)勢(shì)。

表6 各個(gè)方法的綜合結(jié)果

4 工程驗(yàn)證

使用本文方法對(duì)滾齒工藝參數(shù)進(jìn)行預(yù)測(cè)[31]，并在實(shí)際機(jī)床上進(jìn)行加工。工程條件：①數(shù)控高速滾齒機(jī)YS3120CNC6-S，該機(jī)床為重慶某機(jī)床制造公司研發(fā)的六軸四聯(lián)動(dòng)智能高速滾齒機(jī)。配有西門子840Dsl數(shù)控系統(tǒng)，加工模數(shù)最大為6 mm; ②齒坯; ③計(jì)算機(jī)及MATLAB軟件；④三豐高精度數(shù)顯千分尺。實(shí)驗(yàn)基本步驟：將本文方法在MATLAB軟件上實(shí)現(xiàn)，基于歷史工藝參數(shù)決策數(shù)據(jù)生成滾齒參數(shù)，經(jīng)工藝人員確認(rèn)后進(jìn)行編程并加工，使用三豐高精度數(shù)顯千分尺測(cè)量跨棒距，以驗(yàn)證工件是否合格,加工設(shè)備如圖7所示。

歷史工藝參數(shù)決策數(shù)據(jù)如表7所示，種類為漸開線圓柱齒輪(包含非標(biāo)齒輪)，材料為20CrMoTi，刀具為TiAlN涂層滾刀，加工精度為7級(jí)，采用逆滾軸向進(jìn)給方式加工，一次走刀，走刀路徑如圖8所示。待決策問題Q={f1,f2,f3,f4,f5,f6,f7,f8,f9}={2, 0.349, 41, 0.456, 94.011, 13,71, 3, 4.005}，跨棒距要求為95.09±0.025 mm。使用本文方法決策出Q的工藝參數(shù)(382.6 r/min,1.84 mm/r)，決策耗時(shí)4.2 s。工藝人員依據(jù)該工藝參數(shù)編制數(shù)控程序，操作人員依據(jù)數(shù)控程序操作數(shù)控系統(tǒng)對(duì)工件進(jìn)行加工，用三豐千分尺測(cè)量跨棒距為95.106 mm，工件滾齒合格，可以進(jìn)入熱處理及磨齒等后續(xù)工序。

表7 歷史工藝參數(shù)決策數(shù)據(jù)集

將本文方法(CA0)與CA1、CA2、CA3、CA4、CA5、CA6相比，由于沒有Q的歷史決策結(jié)果，無法比較預(yù)測(cè)精度，故本文在工件是否合格、決策耗時(shí)、切削耗時(shí)上進(jìn)行比較。結(jié)果如表8所示，可以發(fā)現(xiàn)，本文方法在決策耗時(shí)上排名第3。在保證質(zhì)量的情況下，切削耗時(shí)是最少的，排名第1。本文方法可以很好且快速地預(yù)測(cè)滾齒工藝參數(shù)，保證加工質(zhì)量，且切削耗時(shí)少，與其他方法相比，具備一定優(yōu)勢(shì),驗(yàn)證了本文方法的工程效果。

表8 各個(gè)算法應(yīng)用效果比較

5 結(jié)束語(yǔ)

本文提出一種利用K-means聚類分析和CHHO優(yōu)化SVR超參數(shù)的新方法，并進(jìn)行了實(shí)驗(yàn)驗(yàn)證與工程驗(yàn)證。K-means聚類用于提高該算法的速度，CHHO用于提高預(yù)測(cè)精度。該方法能夠?qū)崿F(xiàn)預(yù)測(cè)問題中，特別是在大數(shù)據(jù)預(yù)測(cè)問題中的SVR超參數(shù)自適應(yīng)調(diào)整，且能實(shí)現(xiàn)預(yù)測(cè)精度的提高，穩(wěn)定性維持在較高水平。實(shí)驗(yàn)將Kmeans-CHHO-SVR與Kmeans-SVR、Kmeans-HGSO-SVR、Kmeans-SSA-SVR、Kmeans-DA-SVR、Kmeans-ALO-SVR、Kmeans-PSO-SVR、Kmeans-HHO-SVR比較，結(jié)果顯示與其他成熟方法相比，Kmeans-CHHO-SVR的預(yù)測(cè)性能最佳，獲得了最佳綜合得分0.226。針對(duì)MSEV指標(biāo)，本文方法在5組數(shù)據(jù)集上全面戰(zhàn)勝其他方法，為最優(yōu)；針對(duì)MSET指標(biāo)，本文方法在4組數(shù)據(jù)集上戰(zhàn)勝其他方法，為最優(yōu)。針對(duì)StdVRK指標(biāo)，本文方法排名第2；針對(duì)StdTRK指標(biāo)，本文方法排名第3。綜合評(píng)價(jià)顯示本文方法綜合性能最佳。此外，在滾齒工藝參數(shù)預(yù)測(cè)領(lǐng)域進(jìn)行了工程驗(yàn)證，證實(shí)了本文方法的實(shí)際應(yīng)用效果。綜上，可以認(rèn)為Kmeans-CHHO-SVR在預(yù)測(cè)領(lǐng)域擁有巨大潛力。未來將進(jìn)一步提升本文方法的穩(wěn)定性和降低方法運(yùn)行時(shí)間，爭(zhēng)取應(yīng)用于更多更廣的工業(yè)領(lǐng)域。