鮑丹， 趙搶搶， 侯保林

(1.南京理工大學(xué) 機械工程學(xué)院，江蘇 南京 210094； 2.江蘇金陵智造研究院有限公司，江蘇 南京 210006)

火炮反后坐裝置安裝在炮身和架體之間，用于提供彈性力和制動力控制后坐部分在火炮射擊時的后坐運動，并使之復(fù)位[1-2]。在此過程中，反后坐裝置減小了火炮架體在射擊時的受力，還把射擊時的全炮后坐運動變?yōu)榭煽氐呐谏砗笞\動，并能自動復(fù)位。反后坐裝置中的參數(shù)包括基本物理參數(shù)、結(jié)構(gòu)參數(shù)和經(jīng)驗系數(shù)。一些參數(shù)基于實驗或經(jīng)驗的系數(shù)不易準(zhǔn)確選定，因而對這些參數(shù)進行辨識具有一定的實用價值。杜中華等[3]采用數(shù)值仿真的方法對某型反后坐裝置的液壓阻力以及密封件工作壓強進行了分析；狄長春等[4]對反后坐裝置后坐復(fù)進過程中復(fù)雜流場性質(zhì)及結(jié)構(gòu)場動態(tài)特性進行了數(shù)值仿真。此外，火炮反后坐隨著機械系統(tǒng)呈現(xiàn)出高速化、大型化和復(fù)雜化的方向發(fā)展，常規(guī)研究手段已經(jīng)無法準(zhǔn)確描述系統(tǒng)中存在的大量不確定性參數(shù)。對于不確定參數(shù)的優(yōu)化方法可分為3類：概率不確定性優(yōu)化方法、非概率不確定性優(yōu)化方法和概率-非概率混合不確定性優(yōu)化方法。工程實際中，往往只能得到有限的樣本數(shù)據(jù)，很難獲得不確定參數(shù)的精確概率分布，所以概率不確定性優(yōu)化方法在工程實際中存在局限性。相對于概率不確定性優(yōu)化方法，非概率不確定性優(yōu)化方法具有對數(shù)據(jù)要求低，計算簡單等優(yōu)點，越來越多的科研人員對此進行了研究。其中采用區(qū)間來描述事物的本質(zhì)和特征的區(qū)間數(shù)理論方法被應(yīng)用于非概率不確定性優(yōu)化。文獻[5-6]分別提出了基于區(qū)間分析的方法進行參數(shù)估計和結(jié)構(gòu)辨識；姜潮[7]提出了基于區(qū)間的不確定性優(yōu)化理論與算法；Jiang等[8]提出一種改進的區(qū)間可能度計算方法提高計算精度；趙搶搶等[9]利用區(qū)間序關(guān)系轉(zhuǎn)換模型將區(qū)間不確定參數(shù)辨識問題轉(zhuǎn)換成確定性優(yōu)化問題，并結(jié)合差分進化算法實現(xiàn)了彈藥協(xié)調(diào)器的參數(shù)辨識；王敏容等[10]以可靠性指標(biāo)為約束條件，提出了一種基于區(qū)間模型的結(jié)構(gòu)非概率穩(wěn)健優(yōu)化設(shè)計方法。當(dāng)用區(qū)間數(shù)描述系統(tǒng)的不確定參數(shù)時，由于目標(biāo)函數(shù)相對于每一個設(shè)計變量的取值構(gòu)成一個區(qū)間，所以辨識問題無法通過傳統(tǒng)的確定性辨識方法求解。針對存在區(qū)間不確定參數(shù)的優(yōu)化問題，可以通過非線性區(qū)間優(yōu)化的數(shù)學(xué)模型轉(zhuǎn)化為確定性問題進行求解。

本文為解決火炮后坐復(fù)進過程中反后坐裝置的參數(shù)辨識問題，建立了火炮后坐復(fù)進過程的解析模型。為了提高辨識效率，同時建立了神經(jīng)網(wǎng)絡(luò)模型，并根據(jù)不確定參數(shù)的變化特性對相關(guān)性能參數(shù)進行了分類以及用區(qū)間數(shù)進行描述。針對在一次動作過程中可通過參數(shù)辨識的方法得到其確定值的第1類參數(shù)，利用區(qū)間可能度轉(zhuǎn)換模型，將區(qū)間不確定性問題轉(zhuǎn)換為確定性問題進行求解。以反映辨識結(jié)果與測試結(jié)果的相似程度的時間序列數(shù)據(jù)相似度作為優(yōu)化目標(biāo)，利用粒子群算法實現(xiàn)參數(shù)辨識。針對只能確定其變化區(qū)間的參數(shù)，利用粒子群算法結(jié)合區(qū)間優(yōu)化方法優(yōu)化其區(qū)間。

1 火炮反后坐裝置運動過程數(shù)學(xué)建模

1.1 火炮復(fù)進機和制退機數(shù)學(xué)方程

火炮在后坐運動結(jié)束后，后坐部分在復(fù)進機力作用下回到待發(fā)射位置。取后坐部分為受力體,后坐過程中的運動微分方程為[1]：

mdu/dt=m(d2x)/(dt2)=Fpt-FR

(1)

Fpt=Fφh+Ff+F+FT-mgsinφ

(2)

式中：Fpt為炮膛合力；FR為后坐阻力；Fφh為制退機力；Ff為復(fù)進機力；F為反后坐裝置密封裝置的摩擦力；FT為導(dǎo)軌的摩擦力；u為后坐部分后坐速度；x為后坐部分后坐位移；m為后坐部分質(zhì)量；φ為火炮射角。其中，F(xiàn)f、Fφh為：

Ff=AfPf 0(V0/(V0-Afx))n

(3)

(K2ρAfj3)/(2A12)u2

(4)

火炮復(fù)進過程中的運動微分方程為：

mdu/dt=m(d2x)/(dt2)=

Ff-Fφh-(F+FT+mgsinφ)

(5)

Ff=AfPf 0(V0/(V0-Afx))n

(6)

式中:Af為復(fù)進機活塞工作面積；Pf 0為復(fù)進機氣體的初壓力；V0為復(fù)進機氣體的初體積；n為復(fù)進機氣體非線性指數(shù)。后坐時，由于制退桿從制退機內(nèi)抽出而加大了制退機內(nèi)的空間體積，使之形成制退機非工作腔內(nèi)的真空。復(fù)進時，在制退機非工作腔真空消失以前，總的復(fù)進液壓阻力中不含有制退機復(fù)進液壓阻力；當(dāng)制退機非工作腔的真空消失之后，總的液壓阻力包含制退機復(fù)進液壓阻力。

(7)

(8)

(9)

式中：P1為制退機P1腔室的壓力；Ao為制退桿活塞的工作面積；Ap為節(jié)制環(huán)孔面積。

圖1 制退機結(jié)構(gòu)及壓力測試腔示意Fig.1 Recoil brake structure and pressure test chamber figure

1.2 火炮反后坐裝置待辨識參數(shù)分類

由1.1節(jié)可知，火炮反后坐裝置系統(tǒng)包含參數(shù)眾多。待辨識的參數(shù)有3個，分別是復(fù)進機氣體初壓Pf 0、復(fù)進機氣體非線性指數(shù)n和制退機主流液壓阻力系數(shù)K1。根據(jù)3個參數(shù)在發(fā)射過程中的變化特性分為2類。第1類參數(shù)是難以通過測量的方法獲得，在動作過程中可近似看作確定的，可通過辨識的方法獲得，即復(fù)進機氣體初壓Pf 0。Pf 0是復(fù)進機的基本參數(shù)，由于密封裝置的泄露，結(jié)構(gòu)尺寸誤差等各方面原因，使得在發(fā)射過程中很難測得其準(zhǔn)確值。本文根據(jù)在經(jīng)驗值的基礎(chǔ)上適當(dāng)增大范圍的原則確定Pf 0的初始區(qū)間。

第2類參數(shù)是區(qū)間不確定參數(shù)，無法確定其準(zhǔn)確的分布規(guī)律，只能確定其變化的區(qū)間。對于火炮反后坐裝置而言，2個區(qū)間不確定的參數(shù)分別為：復(fù)進機氣體非線性指數(shù)n和制退機主流液壓阻力系數(shù)K1。復(fù)進機氣體非線性指數(shù)n取決于復(fù)進機的散熱條件和活塞運動速度。制退機液壓阻力系數(shù)是反映制退機內(nèi)液體真實流動因素的符合系數(shù)，這些因素有液體流動的沿程損失、流動的局部損失、流液孔處液流截面的收縮、液體可壓縮性和流動的非定常流動等。本文在經(jīng)驗區(qū)間的基礎(chǔ)上適當(dāng)增大區(qū)間范圍的原則確定K1和n的初始區(qū)間，如表1所示。

表1 待辨識參數(shù)初始區(qū)間Table 1 Initial intervals of parameters to be identified

1.3 神經(jīng)網(wǎng)絡(luò)模型的建立

火炮后坐復(fù)進過程是復(fù)雜的運動過程，而辨識過程需要大量的循環(huán)計算，為了提高辨識效率，利用前述數(shù)學(xué)模型，建立神經(jīng)網(wǎng)絡(luò)模型。反向傳播神經(jīng)網(wǎng)絡(luò)(BP神經(jīng)網(wǎng)絡(luò))是一種多層前向反饋神經(jīng)神經(jīng)網(wǎng)絡(luò)，具有強大的自適應(yīng)、泛化和非線性逼近能力。BP神經(jīng)網(wǎng)絡(luò)通常采用梯度下降搜索算法，不斷調(diào)整網(wǎng)絡(luò)權(quán)值和閾值，使得網(wǎng)絡(luò)的輸出值與實際輸出值的均方誤差最小。文章在3個待辨識參數(shù)的初始區(qū)間進行拉丁超立法抽樣1 000組，代入數(shù)學(xué)模型，得到1 000組P1腔壓力曲線，并計算1 000組仿真曲線與測試壓力曲線的曲線相似度。以3個變量為輸入，曲線相似度為輸出，建立神經(jīng)網(wǎng)絡(luò)的1 000組樣本，神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)如圖2所示。

圖2 神經(jīng)網(wǎng)絡(luò)模型的結(jié)構(gòu)Fig.2 The structure of neural network model

2 參數(shù)辨識問題的描述及求解方法

2.1 參數(shù)辨識問題描述

整個參數(shù)辨識過程分為2個步驟。1)對第1類參數(shù)的辨識可利用尋優(yōu)的方法求解，即在主流液壓阻力系數(shù)K1和復(fù)進機氣體非線性指數(shù)n在初始區(qū)間范圍內(nèi)變化的情況下，尋找出最優(yōu)的復(fù)進機氣體初壓Pf 0，使得辨識結(jié)果與測試結(jié)果的時間序列數(shù)據(jù)相似度區(qū)間的可能度最大。其數(shù)學(xué)模型為：

(10)

式中：X=Pf 0為設(shè)計變量，D=(K1,n) 為不確定變量，由于D的取值為一區(qū)間數(shù)，因此對于任何一個確定的設(shè)計變量X，相似度fI(Xi,D) 的取值構(gòu)成一個區(qū)間數(shù)，且均在區(qū)間[0,1]。將相似度區(qū)間轉(zhuǎn)換成區(qū)間可能度P(fI(Xi,D)≥VI)，VI為性能區(qū)間，求解式后，可以得到一個最優(yōu)設(shè)計向量X*，使得目標(biāo)函數(shù)的可能度取最大值Pmax=P(fI(X*,D)≥VI)。傳統(tǒng)的確定性優(yōu)化方法中，決策的判定都是基于目標(biāo)函數(shù)在各個設(shè)計向量處的具體數(shù)值。針對此類目標(biāo)函數(shù)為區(qū)間的參數(shù)辨識問題，本文提出了一種基于區(qū)間可能度數(shù)學(xué)轉(zhuǎn)換模型將區(qū)間不確定優(yōu)化問題轉(zhuǎn)換為確定性問題的方法，再結(jié)合粒子群算法進行辨識。

2)在辨識出第1類參數(shù)的基礎(chǔ)上，對第2類參數(shù)的區(qū)間進行優(yōu)化。每一組(K1,n)對應(yīng)一個相似度值f(K1,n)，如圖3所示。

圖3 參數(shù)值與相似度值對應(yīng)關(guān)系Fig.3 Correspondence between parameters and similarity values

需要找到待辨識參數(shù)制退機的主流液壓阻力系數(shù)K1和復(fù)進機氣體非線性指數(shù)n的最大區(qū)間，使得該區(qū)間內(nèi)所有的相似度值均大于μ。其數(shù)學(xué)模型為：

(11)

(12)

(13)

針對此類區(qū)間優(yōu)化問題，本文采用式(11)的數(shù)學(xué)模型結(jié)合粒子群算法進行優(yōu)化。

2.2 時間序列相似度計算

本文采用離散序列的一致性度量方法，即動態(tài)時間規(guī)整(DTW)計算仿真曲線與測試曲線的相似度來評價辨識結(jié)果。試驗所得的壓力數(shù)據(jù)序列長度與仿真所得的壓力數(shù)據(jù)序列長度不一致，并且2組數(shù)據(jù)序列的時間軸無法完全對齊，所以一般的歐氏距離法存在不足。DTW方法用滿足一定條件的時間規(guī)整函數(shù)描述兩者之間的時間對應(yīng)關(guān)系，它允許時間軸上的漂移。將DTW距離進行歸一化處理，所得的值即為2組時間序列數(shù)據(jù)之間的相似度[11-13]。動態(tài)時間規(guī)劃的具體步驟為：

1)假設(shè)2組時間序列分別為n維和m維，構(gòu)造一個n×m的矩陣，用于存放兩序列點對點之間的距離(一般可采用歐氏距離)，距離越小表示兩點之間的相似度越高。

2)把矩陣看成一個網(wǎng)格，算法的目的為找到一條通過此矩陣網(wǎng)格的最有路徑，改路徑通過的網(wǎng)格點即為兩組時間序列經(jīng)過對齊后的點對；

3)DTW算法定義一個歸整路徑距離，找到最優(yōu)路徑后，將所有相似點之間距離和進行歸一化處理，來衡量2組時間序列之間的相似性。

2.3 區(qū)間可能度轉(zhuǎn)換模型

根據(jù)區(qū)間數(shù)學(xué)，區(qū)間數(shù)被定義為一對有序的實數(shù)[14]：

AI=[AL,AR]={xAL≤x≤AR,x∈R}

(14)

式中上標(biāo)I、L、R分別表示區(qū)間、區(qū)間下界和區(qū)間上界。區(qū)間的寬度定義為：

L(AI)=AR-AL

(15)

區(qū)間可能度方法是基于模糊集來構(gòu)造區(qū)間可能度，是一種定量比較區(qū)間之間的優(yōu)劣程度[15]。將區(qū)間B視為文章辨識的性能區(qū)間，則區(qū)間A與區(qū)間B的位置關(guān)系如圖4所示。

圖4 區(qū)間AI和BI的3種位置關(guān)系Fig.4 Three kinds of positional relationships between interval AI and BI

記P(AI≥BI)為區(qū)間AI≥BI的可能度。本文采用了2種不同區(qū)間可能度的構(gòu)造方法:

構(gòu)造方法(Ⅰ):

P(AI≥BI)=max[0,L(AI)+L(BI)-max(0,BR-AL)](L(AI)+L(BI))-1

(16)

構(gòu)造方法(Ⅱ)：

(17)

2.4 粒子群算法

本文對2類參數(shù)的辨識均采用了粒子群算法[16-18]，該算法是智能優(yōu)化算法的一種，有著精度高、收斂快的特點。粒子群算法優(yōu)化的原理是通過不斷更新粒子速度和位置，并計算其適應(yīng)度值，直到適應(yīng)度值符合要求為止。粒子群算法的流程如圖5所示。

圖5 粒子群優(yōu)化算法的流程Fig.5 The process of particle swarm optimization algorithm

根據(jù)2類參數(shù)的不同辨識要求，設(shè)置粒子群算法的參數(shù)，包括種群規(guī)模，最大迭代次數(shù)、學(xué)習(xí)因子、適應(yīng)度函數(shù)以及粒子位置更新準(zhǔn)則。

3 火炮發(fā)射試驗與辨識結(jié)果分析

3.1 火炮發(fā)射試驗

文中目標(biāo)函數(shù)是通過數(shù)學(xué)解析模型所得的仿真曲線與試驗曲線之間的相似度。測試曲線是在火炮發(fā)射試驗過程中，測得的P1腔室的壓力曲線。圖6為試驗過程示意圖，圖7為測得的壓力曲線。

圖6 發(fā)射試驗Fig.6 The launch test

圖7 P1腔室的壓力測試曲線Fig.7 The pressure test curves of P1 cavity

3.2 辨識結(jié)果分析

3.2.1 第1類參數(shù)辨識過程及結(jié)果

設(shè)置粒子群智能算法的相關(guān)參數(shù)：種群規(guī)模為100;最大迭代次數(shù)為50；適應(yīng)度函數(shù)即為區(qū)間可能度；學(xué)習(xí)因子c1=0.5,c2=0.5；更新權(quán)重w=0.8。粒子速度及位置更新準(zhǔn)則為：

vid(k+1)=wvid(k)+c1r1(pid-xid(k))+
c2r2(pgd-xid(k))

(18)

xid(k+1)=xid(k+1)+vid(k+1)

(19)

式中:k是迭代次數(shù)；r1、r2均為[-1,1]的均勻隨機數(shù)；pid、pgd分別為個體最優(yōu)值和全局最優(yōu)值。分別采用2種區(qū)間可能度的構(gòu)造方法對第1類參數(shù)進行辨識，性能區(qū)間VI為[0.8,1.0]。當(dāng)采用第1種可能度構(gòu)造方法時，可得復(fù)進機氣體初壓Pf 0的最優(yōu)值為7.588 2 MPa,所用時間為6.708 s；當(dāng)采用第2種區(qū)間可能度構(gòu)造方法時，可得復(fù)進機氣體初壓Pf 0的最優(yōu)值為7.342 0 MPa，所用時間為6.563 s。然后，K1和n在初始區(qū)間采樣，分別同2種可能度構(gòu)造法下得到的復(fù)進機氣體初壓Pf 0最優(yōu)值代入模型，可得P1腔室壓力曲線簇與測試曲線的對比，即如圖8和圖9所示。分別采用2種區(qū)間可能度的構(gòu)造方法辨識出復(fù)進機氣體初壓值。由圖8和圖9可得，采用第2種可能度構(gòu)造方法得到的復(fù)進機氣體初壓值代入神經(jīng)網(wǎng)絡(luò)模型，得到的P1腔室壓力曲線簇更加穩(wěn)定，且與測試曲線的整體相似度更高。因此，復(fù)進機氣體初壓Pf 0的最優(yōu)值取7.342 0 MPa。以此為基礎(chǔ)，進行第2類參數(shù)的區(qū)間優(yōu)化。圖8中，該結(jié)果是在采用第1種可能度構(gòu)造法的情況下獲得。黑色曲線為測試曲線，其他不同顏色曲線為K1和n在不同采樣樣本下得到的辨識曲線。圖9中，該結(jié)果是在采用第2種可能度構(gòu)造法的情況下獲得。黑色曲線為測試曲線，其他不同顏色曲線為K1和n在不同采樣樣本下得到的辨識曲線。

圖8 辨識曲線簇與測試曲線的對比ⅠFig.8 Comparison between identification curve cluster and test curve(Ⅰ)

圖9 辨識曲線簇與測試曲線的對比ⅡFig.9 Comparison between identification curve cluster and test curve(Ⅱ)

3.2.2 第2類參數(shù)辨識過程及結(jié)果

設(shè)置粒子群智能算法的相關(guān)參數(shù)：種群規(guī)模為100；最大迭代次數(shù)為200；c1=1.2,c2=1.3；適應(yīng)度函數(shù)為Lij→rk(x,y)=(xr-xi)2+(yk-yj)2；速度位置更新準(zhǔn)則如式(18)和式(19)。μ分別取0.50、0.80、0.95進行優(yōu)化，得到如表2結(jié)果。

為驗證該算法的優(yōu)化精度，當(dāng)μ取0.50時，參數(shù)K1和n在優(yōu)化后所得區(qū)間內(nèi)隨機采樣1 000組，代入數(shù)學(xué)模型中，得到P1腔室壓力曲線簇，并計算相似度，得到相似度在優(yōu)化后相似度區(qū)間范圍內(nèi)的概率為97.4%。同理，當(dāng)μ分別取0.80、0.95時，參數(shù)K1和n在優(yōu)化后所得區(qū)間內(nèi)隨機采樣1 000組，帶入數(shù)學(xué)模型中，并計算P1腔室壓力曲線簇的相似度，可得相似度值在優(yōu)化后相似度區(qū)間內(nèi)的概率分別為97.2%和98.0%。3種情況相似度均較高，可得本文的算法具有較高的優(yōu)化精度。

表2第2類參數(shù)辨識結(jié)果

Table2Thesecondkindofparameteridentificationresults

μ參數(shù)K區(qū)間參數(shù)n區(qū)間相似度區(qū)間0.50[1.489,2.000][1.045,1.826][0.500,0.922]0.80[1.783,1.788][1.083,1.801][0.800,0.883]0.951.621[1.530,1.767][0.950,0.985]

4 結(jié)論

1)為了辨識火炮反后坐裝置不確定性參數(shù)，文章建立了火炮發(fā)射過程中，反后坐裝置的數(shù)學(xué)解析模型，并為提高辨識效率，建立神經(jīng)網(wǎng)絡(luò)模型。為了得到P1腔室的壓力曲線作為辨識過程中的目標(biāo)曲線，進行了火炮發(fā)射試驗。

2)火炮復(fù)進機的氣體初壓屬于第1類參數(shù)，針對此類不確定參數(shù)，利用區(qū)間可能度數(shù)學(xué)轉(zhuǎn)換模型，將區(qū)間不確定優(yōu)化問題轉(zhuǎn)化為確定性優(yōu)化問題進行求解。以測試曲線和數(shù)值仿真曲線的時間序列數(shù)據(jù)相似度作為優(yōu)化目標(biāo)，利用粒子群算法實現(xiàn)參數(shù)辨識。

3)火炮制退機的主流液壓阻力系數(shù)和復(fù)進機的非線性指數(shù)屬于第2類參數(shù)。針對此類無法確定其準(zhǔn)確的分布規(guī)律，只能確定其變化的區(qū)間，文章提出了一種基于粒子群算法結(jié)合區(qū)間優(yōu)化算法優(yōu)化其區(qū)間范圍。辨識結(jié)果表明文章提出的同時針對2種不同類型參數(shù)的辨識方法具有可行性和有效性。