宋李俊，潘安大，Jing Shi，陳 猛

(1.重慶理工大學(xué) 機(jī)械工程學(xué)院， 重慶 400054；2.辛辛那提大學(xué) 機(jī)械與材料科學(xué)學(xué)院， 辛辛那提 45242)

目前，受到市場需求的波動(dòng)、新技術(shù)的引入以及加工過程中工況復(fù)雜多樣，使得數(shù)控機(jī)床能耗特性變得更加復(fù)雜。針對此現(xiàn)象，許多學(xué)者從不同的角度對機(jī)床能耗進(jìn)行了研究。文獻(xiàn)[1]從普通機(jī)床機(jī)械加工入手對其過程中的功率傳遞和分配、能量來源和消耗、加工路徑等進(jìn)行了分析，提出了普通機(jī)床能耗的物理模型。文獻(xiàn)[2]提出表征機(jī)床自身能量效率的機(jī)床固有能量效率新概念，并基于虛擬工件法，實(shí)現(xiàn)機(jī)床固有能量效率評價(jià)，提升了機(jī)床使用生命周期全過程的綜合能量效率。文獻(xiàn)[3]詳細(xì)研究了機(jī)床加工過程中伺服驅(qū)動(dòng)、主軸傳動(dòng)、液壓以及冷卻潤滑等系統(tǒng)運(yùn)作所產(chǎn)生的電能消耗。文獻(xiàn)[4]基于數(shù)據(jù)驅(qū)動(dòng)，建立了數(shù)控機(jī)床實(shí)時(shí)分布式反饋控制策略，識別和預(yù)測節(jié)能機(jī)會(huì)，提高了當(dāng)前系統(tǒng)的預(yù)測性能、整體利潤和能源效率。文獻(xiàn)[5]采用BP神經(jīng)網(wǎng)絡(luò)建立數(shù)控機(jī)床能耗模型，運(yùn)用遺傳算法對切削工藝參數(shù)進(jìn)行優(yōu)化，達(dá)到簡化傳統(tǒng)經(jīng)驗(yàn)公式計(jì)算過程和降低能耗的目的。上述研究，機(jī)床能耗建模方法大多基于機(jī)床加工物理原理，需要大量的物理參數(shù)，這些參數(shù)通常難以計(jì)算或估計(jì)；同時(shí)，在面對多變的加工環(huán)境下，需要處理大量多樣的數(shù)據(jù)，導(dǎo)致建模困難、能耗預(yù)測不及時(shí)，出現(xiàn)嚴(yán)重的分析和優(yōu)化決策滯后的現(xiàn)象。

高斯過程回歸(gaussian process regression，GPR)是近年發(fā)展起來的基于貝葉斯框架的一種機(jī)器學(xué)習(xí)回歸方法，有著嚴(yán)格的統(tǒng)計(jì)學(xué)習(xí)理論基礎(chǔ)，對處理高維數(shù)、小樣本、非線性等復(fù)雜建模具有很好的適應(yīng)性和泛化能力，同時(shí)能夠提供一套完整的理論、量化預(yù)測值的不確定性、預(yù)測輸出的可靠性判斷，以及能夠以增量的方式更新回歸模型的能力等特點(diǎn)[6]，近幾年逐步得到國內(nèi)學(xué)者的重視，在許多領(lǐng)域得到了成功應(yīng)用。

目前數(shù)控機(jī)床能耗預(yù)測的研究主要采用神經(jīng)網(wǎng)絡(luò)、支持向量機(jī)等方法，相比于GPR具有容易實(shí)現(xiàn)、模型適應(yīng)性、泛化能力好以及超參數(shù)自適應(yīng)獲取等優(yōu)點(diǎn)[6]。因此，采用高斯過程回歸建立數(shù)控機(jī)床能耗預(yù)測模型，以此降低模型的復(fù)雜度，實(shí)現(xiàn)實(shí)時(shí)的能耗預(yù)測，同時(shí)對能耗建模、優(yōu)化研究以及降低工業(yè)能耗、實(shí)現(xiàn)綠色生產(chǎn)等提供一種方法。

1 數(shù)控機(jī)床能耗因素分析

數(shù)控機(jī)床能耗特性、影響因素是數(shù)控機(jī)床能耗建模、預(yù)測及優(yōu)化的基礎(chǔ)。從系統(tǒng)層面來看，數(shù)控機(jī)床各部件的能耗特點(diǎn)可分為兩類：加工動(dòng)力系統(tǒng)能耗源和輔助系統(tǒng)能耗源。加工系統(tǒng)的能量主要用于工件切削加工，可進(jìn)一步分為主傳動(dòng)和進(jìn)給傳動(dòng)能耗系統(tǒng)，在此基礎(chǔ)上，再進(jìn)一步劃分到每個(gè)能耗部件上；液壓系統(tǒng)、潤滑冷卻系統(tǒng)、輔助控制系統(tǒng)、照明指示系統(tǒng)等的能耗與數(shù)控機(jī)床加工動(dòng)力系統(tǒng)能耗相對獨(dú)立，統(tǒng)稱為輔助系統(tǒng)。具體如圖1所示。

圖1 數(shù)控機(jī)床能耗系統(tǒng)構(gòu)成框架

上述分析中數(shù)控機(jī)床能量源眾多，能耗形式復(fù)雜多樣，最終作用于切削加工；在切削區(qū)域，工件切削層發(fā)生塑性變形主要受到切削力、切削熱、切削溫度、刀具角度、加工材料等因素以及它們之間相互作用的影響，其中切削力決定著切削熱的產(chǎn)生，并影響刀具的磨損和已加工表面的質(zhì)量[7]。此外，根據(jù)切削力Fc經(jīng)驗(yàn)公式、切削功率Pc公式，本文選取切削速度、進(jìn)給量、切削深度作為數(shù)控機(jī)床切削能耗預(yù)測模型的輸入變量。

Fc=CFcapxFcfyFcKFc

(1)

Pc=Fcvc10-3

(2)

式中：CFc為工件材料、切削條件所決定的系數(shù)；xFc、yFc為指數(shù)常數(shù)；KFc為實(shí)際與經(jīng)驗(yàn)切削力不符時(shí)，各種因素對切削力的修正系數(shù)之積。

2 基于高斯過程回歸的數(shù)控機(jī)床切削能耗預(yù)測模型

2.1 高斯過程預(yù)測

高斯過程具有嚴(yán)格的統(tǒng)計(jì)學(xué)習(xí)理論基礎(chǔ)，性質(zhì)完全由它的均值函數(shù)m(x)和協(xié)方差函數(shù)k(x，x′)確定，公式如下：

f(x)～GP(m(x))，k(x，x′))

(3)

其中x，x′∈Rl為任意隨機(jī)變量。

(4)

由式(4)推導(dǎo)出觀測y和預(yù)測值y*的聯(lián)合先驗(yàn)分布：

(5)

式中：K(X，X)為n×n階對稱正定協(xié)方差矩陣；K(X，x*)為測試點(diǎn)x*與訓(xùn)練集X中各輸入點(diǎn)的n×1階協(xié)方差矩陣；k(x*，x*)為測試點(diǎn)x*自身的協(xié)方差矩陣；In為n階單位矩陣。

根據(jù)貝葉斯后驗(yàn)概率公式，推導(dǎo)出預(yù)測值y*的后驗(yàn)分布

(6)

式中：

(7)

cov(y*)=k(x*，x*)-K(x*，X)K(X，x*)×

(8)

2.2 高斯過程訓(xùn)練

考慮目標(biāo)函數(shù)的特征以及能量消耗的變化呈平穩(wěn)性[8]，本文選擇自動(dòng)相關(guān)性測性型核函數(shù)，即平方指數(shù)核函數(shù)，公式如下：

(9)

對于超參數(shù)最優(yōu)解的求取，一般通過對訓(xùn)練樣本條件概率負(fù)對數(shù)似然函數(shù)L(θ)=-logp(y|X，θ)的超參數(shù)θ求偏導(dǎo)，采用共軛梯度法對偏導(dǎo)數(shù)進(jìn)行最小化求解獲得。求得到超參數(shù)最優(yōu)解，再利用式(7)和(8)，則可得到測試點(diǎn)x*對應(yīng)y*及其方差σy*。

3 差分進(jìn)化-高斯回歸過程回歸模型

高斯過程回歸超參數(shù)集θ對預(yù)測結(jié)果影響很大，當(dāng)前最優(yōu)參數(shù)常采用共軛梯度算法來確定，但共軛梯度法存在優(yōu)化效果過于依賴初始值、迭代次數(shù)難定以及難以獲得全局最優(yōu)解等缺陷，不能建立整體性能最優(yōu)的預(yù)測模型。針對上述問題，在最優(yōu)超參數(shù)求解過程中，引入DE算法。

3.1 差分進(jìn)化算法DE

差分進(jìn)化算法(differential evolution，簡稱DE)，具有收斂快、控制參數(shù)少、優(yōu)化結(jié)果穩(wěn)健等優(yōu)點(diǎn)。其流程如下：

1) 初始化種群

隨機(jī)產(chǎn)生初始化種群：

X={xi(0)=[xi1，xi2，…，xiD]}

xij=aj+rand×(bj-aj)

i=1，2，…，NP

j=1，2，…，D

(10)

2) 變異操作

從第g代種群中隨機(jī)抽取3個(gè)不同的個(gè)體xr(g)、xs(g)、xt(g)，對種群中的每個(gè)個(gè)體xi(g)進(jìn)行如下變異操作，產(chǎn)生新的個(gè)體：

ui(g+1)=xr(g)+F·(xs(g)+xt(g))

(11)

式中，F(xiàn)∈[0，1]為縮放因子。

3) 交叉操作

對第g代種群X(g)中的每個(gè)向量xi(g)與其進(jìn)行變異操作后得到的新個(gè)體ui(g+1)進(jìn)行交叉操作，產(chǎn)生第g+1代的個(gè)體vi(g+1)。為了能夠充分利用各代種群間的信息，本文采用二項(xiàng)交叉方式。

(12)

式中：CR為交叉概率，取值范圍(0,1);jrand為隨機(jī)整數(shù)，取值范圍區(qū)間[1，D]。

最后，對vij(g+1)的有效性進(jìn)行判斷，若滿足則保留。公式如下：

vij(g+1)=

(13)

4) 選擇操作

DE算法選擇進(jìn)入g+1代種群中的個(gè)體，采用貪婪策略，即根據(jù)目標(biāo)函數(shù)值的大小來選擇：

(14)

5) 終止條件

如果迭代次數(shù)g超過最大迭代次數(shù)或者求解精度達(dá)到要求時(shí)，則停止搜索；否則，對種群再次進(jìn)行變異、交叉和選擇操作，直到滿足條件為止。

3.2 基于差分進(jìn)化算法高斯過程回歸預(yù)測框架設(shè)計(jì)

基于Matlab平臺，采用DE算法對樣本訓(xùn)練過程的GPR進(jìn)行最優(yōu)超參數(shù)自動(dòng)搜索，形成差分進(jìn)化—高斯過程回歸模型(DE-GPR)，流程如圖2所示。該預(yù)測框架步驟如下：

步驟1 將樣本隨機(jī)分為學(xué)習(xí)樣本和測試樣本。

步驟2 確定DE進(jìn)化模式和初始化控制參數(shù)。本文采取DE/rand/1/bin進(jìn)化模式、交叉概率CR=0.2、種群數(shù)量Np=30、變異算子中縮放因子F=0.8、進(jìn)化次數(shù)200。

步驟3 DE-GPR模型讀入學(xué)習(xí)樣本、測試樣本以及GPR模型超參數(shù)種群；對模型進(jìn)行訓(xùn)練，獲得預(yù)測能力之后，再對預(yù)測樣本進(jìn)行預(yù)測。

步驟4 根據(jù)DE算法適應(yīng)度函數(shù)計(jì)算每個(gè)個(gè)體的適應(yīng)度。這里適應(yīng)度函數(shù)為：

(15)

式中：f(xi)為第i個(gè)測試樣本的預(yù)測值；yi為第i個(gè)測試樣本的樣本值;M為測試樣本個(gè)數(shù)。

步驟5 判斷迭代次數(shù)g是否達(dá)到最大值。如達(dá)到，則算法結(jié)束，返回當(dāng)前適應(yīng)度最高的個(gè)體；反之，進(jìn)入下一循環(huán)。

步驟6 進(jìn)行變異、交叉和選擇操作，得到適應(yīng)度較高的個(gè)體，形成與初代種群相同規(guī)模的新一代種群，計(jì)數(shù)器g=g+1；進(jìn)入步驟3。

步驟7 重復(fù)步驟3～6，當(dāng)?shù)螖?shù)達(dá)到最大值時(shí)，算法結(jié)束，返回GPR模型最優(yōu)超參數(shù)，并確定GPR模型。

圖2 DE-GPR預(yù)測流程

4 實(shí)驗(yàn)仿真

4.1 實(shí)驗(yàn)方案

本文機(jī)床選用CKA6136P數(shù)控機(jī)床，主軸轉(zhuǎn)速32～2 000 r/min、功率4.5 kW；加工工件材料40Cr、長50 mm、直徑40 mm；切削速度vc、進(jìn)給量f、背吃刀量ap作為模型的輸入變量；為了更加準(zhǔn)確地反映數(shù)控機(jī)床切削加工過程中的能耗，采用比能(specific energy consumption，SEC)作為能耗數(shù)控機(jī)床能耗模型的輸出變量，公式如下

(16)

式中：E為加工能耗；V為去除材料體積。

根據(jù)CKA6136P數(shù)控機(jī)床主要參數(shù)，設(shè)定主軸轉(zhuǎn)速、進(jìn)給速度，通過機(jī)床倍率參數(shù)來改變數(shù)值；切削深度在0.5～2 mm進(jìn)行逐級變化，編寫數(shù)控程序進(jìn)行循環(huán)執(zhí)行；切削速度根據(jù)實(shí)際的工件直徑與主軸轉(zhuǎn)速計(jì)算得出；在機(jī)床總電源處獲取總電壓和總電流，并記錄完成每道工序的時(shí)間，以此計(jì)算出每道工序的加工能耗；記錄每道工序加工前后的工件半徑，結(jié)合該道工序的進(jìn)給量得到去除材料的體積。經(jīng)整理，部分?jǐn)?shù)據(jù)見表1。

為了驗(yàn)證DE-GPR模型的預(yù)測結(jié)果的有效性與可行性，使用傳統(tǒng)共軛梯度方法計(jì)算超參數(shù)的GPR模型，BP神經(jīng)網(wǎng)絡(luò)(back propagation neural network，BPNN)模型作為對比；采用均方根誤差(root mean square error，RMSE)和平均相對百分比誤差( Average Relative Percentage Error，MPAE)作為模型性能的評價(jià)指標(biāo)，公式如下:

(17)

(18)

表1 部分實(shí)驗(yàn)數(shù)據(jù)

4.2 實(shí)驗(yàn)結(jié)果分析

采用DE-GPR方法建立數(shù)控機(jī)床能耗預(yù)測模型，基于Matlab2016b編程實(shí)現(xiàn)模型仿真驗(yàn)證和對比，在運(yùn)行平臺Windows10、Interi5 CPU2.6 GHz、4GB RAM上對模型進(jìn)行實(shí)驗(yàn)。BP、GPR、DE-GPR模型的預(yù)測結(jié)果和預(yù)測誤差變化曲線分別如圖3、4所示。

圖3 BPNN、GPR、DE-GPR模型預(yù)測結(jié)果

圖4 BPNN、GPR、DE-GPR模型預(yù)測誤差變化曲線

從圖3可以看出：在對小樣本能耗數(shù)據(jù)進(jìn)行處理時(shí)，相比BPNN，GPR、DE-GPR對能耗的變化能夠很好地跟蹤預(yù)測值，與實(shí)際值比較接近。從圖4中可以看出：GPR 、DE-GPR預(yù)測誤差波動(dòng)范圍相比于BP小，預(yù)測更加精確、穩(wěn)定。

BPNN、GPR、DE-GPR的預(yù)測結(jié)果和誤差見表2、3。對預(yù)測誤差表3進(jìn)行分析，可以得到以下結(jié)論：

1) 與BPNN模型相比，在小樣本數(shù)據(jù)條件下，采用GPR、DE-GPR模型的預(yù)測精度有所提高、預(yù)測誤差有所下降。這表明在數(shù)控機(jī)床加工過程中，采用GPR建模的能耗預(yù)測效果比采用BP神經(jīng)網(wǎng)絡(luò)更好。

2) 相比GPR模型，DE-GPR模型采用了DE算法對GRP超參數(shù)進(jìn)行了優(yōu)化，獲得了全局最優(yōu)解，克服了采用共軛梯度法對初始值敏感、易限入局部最優(yōu)解的缺陷，使其預(yù)測誤差有所下降。

表2 GPR、BPNN、DE-GPR模型預(yù)測結(jié)果對比

表3 GPR、BPNN、DE-GPR模型的誤差對比

5 結(jié)論

采用高斯過程回歸方法建立數(shù)控機(jī)床預(yù)測模型取得了以下成果：

1) 降低建立數(shù)控機(jī)床預(yù)測模型復(fù)雜性，實(shí)現(xiàn)及時(shí)的能耗預(yù)測；

2) 引入差分進(jìn)化算法對所建的數(shù)控機(jī)床預(yù)測模型進(jìn)行優(yōu)化，提高了預(yù)測的精度，并通過實(shí)驗(yàn)仿真，驗(yàn)證了所建模型的有效性和可行性；

3) 為數(shù)控機(jī)床的能耗建模、預(yù)測、節(jié)能優(yōu)化研究提供一種新的方法，同時(shí)也為解決工業(yè)能耗、實(shí)現(xiàn)綠色生產(chǎn)等實(shí)際問題建立基礎(chǔ)。

目前，本文所采用的數(shù)控機(jī)床能耗預(yù)測模型，只從工藝參數(shù)與能耗關(guān)系的角度出發(fā)，未能更加全面精確地預(yù)測機(jī)床能耗。在未來，計(jì)劃進(jìn)行進(jìn)一步的實(shí)驗(yàn)研究，以提高能耗預(yù)測模型的穩(wěn)健性和普遍性。