柴 勁,張鵬飛,張 意,齊竹昌,王天明,韓旭東
(西安現(xiàn)代控制技術(shù)研究所, 西安 710000)
升阻特性作為飛行器的基本氣動特性,對總體射程指標(biāo)以及滑翔方案彈道設(shè)計(jì)等影響較大,因此在項(xiàng)目前期論證以及飛行試驗(yàn)中被重點(diǎn)關(guān)注[1-4]。在得到飛行試驗(yàn)數(shù)據(jù)之后,工程實(shí)踐中常用2種方法進(jìn)行參數(shù)辨識:一種方法是彈道復(fù)合法,即在彈道程序中對氣動插值結(jié)果乘以修正系數(shù),以使速度、彈道傾角、彈道等外彈道參數(shù)接近實(shí)測值,這種方法工作量較大且過程繁瑣,辨識過程依賴工程師的經(jīng)驗(yàn)判斷,辨識結(jié)果精度有限;另一種辨識方法是通過外彈道數(shù)據(jù)求得彈體的受力信息,根據(jù)質(zhì)心動力學(xué)方程直接求解升力阻力,從而得到阻力系數(shù)與升力系數(shù)的辨識值,該方法簡單方便,但需要對觀測數(shù)據(jù)直接進(jìn)行微分,辨識噪聲及辨識誤差較大。
目前常見的氣動辨識算法有最小二乘法[5-7]、極大似然估計(jì)[8-9]、卡爾曼濾波[10-11]以及神經(jīng)網(wǎng)絡(luò)[12-13]、蜂群算法[14]、基于學(xué)生心理的啟發(fā)式優(yōu)化算法[15]等智能算法。這些智能算法模型復(fù)雜且對觀測量及其誤差分布等信息需求較多,不利于工程應(yīng)用。
本文中的辨識算法受模型預(yù)測靜態(tài)規(guī)劃(MPSP)[16-17]啟發(fā),該算法多用于帶有終端約束的制導(dǎo)律設(shè)計(jì),通過狀態(tài)靈敏度矩陣建立了控制變量與終端誤差的聯(lián)系。本文中將待辨識氣動參數(shù)類比制導(dǎo)律中的控制量,將辨識誤差類比制導(dǎo)律中的終端約束誤差,建立起辨識氣動參數(shù)與外彈道數(shù)據(jù)誤差之間的聯(lián)系,通過算法自動調(diào)整氣動參數(shù),從而實(shí)現(xiàn)仿真彈道與外彈道數(shù)據(jù)誤差小于指定精度。本文中提出的氣動辨識方法簡單新穎,可通過編制程序自動執(zhí)行,無需通過人工手動復(fù)合彈道,減少辨識工作量,提高辨識精度,為飛行試驗(yàn)的參數(shù)辨識提供了一種新思路。
暫時不考慮有助推發(fā)動機(jī)的情況,使用簡化的縱向彈體質(zhì)心動力學(xué)方程作為辨識模型,即:
(1)
式中:V、θ、x、y、q、S、cx、cy、m、g分別為速度、彈道傾角、發(fā)射系射距、發(fā)射系高度、動壓、參考面積、阻力系數(shù)、升力系數(shù)、質(zhì)量、重力加速度。
由于狀態(tài)變量存在量級差異,先進(jìn)行歸一化處理。使用最大、最小標(biāo)準(zhǔn)化處理,歸一化處理算子為:
(2)
將歸一化后的變量代入式(1)中,即得到歸一化后的縱向質(zhì)心運(yùn)動模型。
首先確定辨識步長,在每個辨識步長內(nèi)應(yīng)用MPSP算法進(jìn)行迭代計(jì)算,使?fàn)顟B(tài)誤差收斂至誤差閾值以下,即得到待辨識參數(shù)的估計(jì)值。算法流程圖如圖1所示。
圖1 算法流程Fig.1 Algorithm flow chart
提供外彈道數(shù)據(jù)、理論氣動參數(shù)以及彈體結(jié)構(gòu)參數(shù)等辨識輸入。
辨識步長設(shè)置要根據(jù)飛行狀態(tài)確定。若彈體全程飛行狀態(tài)變化范圍不大,辨識步長可以適當(dāng)取長;若飛行狀態(tài)變化劇烈,辨識步長需減小。
(3)
采用歐拉法對式(3)進(jìn)行離散處理(為簡化表示,忽略步長上標(biāo)p,不失一般性,下同),得:
Xk=Xk-1+hf(Xk-1,λk-1)=Fk-1(Xk-1,λk-1)
(4)
式中:h為數(shù)值積分步長;k=1,2,…,N為離散點(diǎn)。
考慮式(4)全微分形式為:
(5)
式中,
(6)
(7)
(8)
(9)
令k=N,N-1,…,2,將式(5)循環(huán)代入式(9)中可得終端誤差近似為:
(10)
式(10)中,
(11)
由于初始狀態(tài)量為常量,故dX1=0。式(10)可簡化為:
(12)
式(12)通過靈敏度矩陣Bk建立了終端誤差ΔYN與待辨識參數(shù)變化量dλk之間的聯(lián)系。其物理含義是,當(dāng)實(shí)際外彈道數(shù)據(jù)與辨識模型的狀態(tài)量間存在誤差時,可通過調(diào)整辨識模型中的待辨識參數(shù)去減小該誤差,直到其滿足精度要求。理論上滿足誤差要求的辨識參數(shù)組合并不唯一。但存在一組既滿足誤差要求,同時與辨識初值相比變化量又最小的一組解作為“較優(yōu)”的辨識結(jié)果。因此可以通過構(gòu)造帶有誤差要求約束、最小化辨識參數(shù)變化量的優(yōu)化問題。其目標(biāo)函數(shù)為:
(13)
式中,Rk為權(quán)重矩陣,可對辨識參數(shù)的調(diào)整量影響進(jìn)行加權(quán)。對式(12)和式(13)組成的優(yōu)化問題,可通過拉格朗日乘子法轉(zhuǎn)化為無約束優(yōu)化問題求解。省去求解過程,直接給出最優(yōu)解的解析形式為:
(14)
待辨識參數(shù)更新為:
(15)
至此,基于MPSP的氣動辨識算法推導(dǎo)完畢。將其整理為:
步驟1提供辨識算法需要的外彈道數(shù)據(jù)、初始?xì)鈩訁?shù)及彈體相關(guān)參數(shù)作為辨識輸入。
步驟2設(shè)置辨識步長。
步驟4計(jì)算矩陣Bk、Aλ、Bk。
以某型155 mm制導(dǎo)炮彈為例,驗(yàn)證本文中辨識算法的有效性。主要彈體參數(shù)及仿真參數(shù)設(shè)置如表1所示。其他仿真條件設(shè)置為:大氣模型使用理想大氣。使用真實(shí)阻力系數(shù)和升力系數(shù)進(jìn)行彈道仿真,生成的速度、彈道傾角、射程、彈道高作為辨識輸入。阻力系數(shù)初值使用零攻角和馬赫數(shù)進(jìn)行插值計(jì)算,并將插值結(jié)果拉偏+30%作為辨識初值。升力系數(shù)的初始迭代值設(shè)為零。權(quán)重矩陣Rk設(shè)為:
表1 參數(shù)設(shè)置
3.2.1辨識精度分析
根據(jù)以上設(shè)置的仿真條件進(jìn)行參數(shù)辨識。辨識步長與每個辨識步長對應(yīng)的迭代次數(shù)如圖2所示??偙孀R步長為47步。每個步長內(nèi)迭代次數(shù)在4~6次,且迭代5次的情況占多數(shù)。可見算法在彈道全段均具有較好的收斂性與一致性。
圖2 辨識步長及每個步長內(nèi)的迭代次數(shù)
辨識后的彈道狀態(tài)量與實(shí)際彈道狀態(tài)量對比如圖3—圖6所示??梢姳孀R后的彈道速度、彈道傾角、射距、彈道高與真實(shí)實(shí)際值基本重合,收斂誤差小于指定精度要求。
圖3 真實(shí)速度與參數(shù)辨識后的速度對比
圖4 真實(shí)彈道傾角與參數(shù)辨識后的彈道傾角對比
圖5 真實(shí)射距與參數(shù)辨識后的射距對比
圖6 真實(shí)射高與參數(shù)辨識后的射高對比
阻力系數(shù)與升力系數(shù)辨識結(jié)果如圖7、圖8所示。
圖7 阻力系數(shù)(真實(shí)、初始值以及辨識值)對比
圖8 升力系數(shù)(真實(shí)、初始值以及辨識值)對比
由圖7、圖8可見,在初始阻力系數(shù)、升力系數(shù)存在誤差的情況下,經(jīng)過迭代,辨識結(jié)果可以較好地收斂至真值。辨識誤差如圖9所示,無控段阻力系數(shù)誤差除,在開舵等個別點(diǎn)達(dá)到2%,其余絕大部分飛行段誤差均在1%以內(nèi)。有控段阻力系數(shù)誤差在0.5%以內(nèi)甚至更小。升力系數(shù)誤差基本在0.5%~1%以內(nèi)。
圖9 阻力系數(shù)、升力系數(shù)辨識誤差
3.2.2算法在線計(jì)算能力分析
從算法復(fù)雜度的角度來看,假設(shè)每個辨識步長內(nèi)算法迭代次數(shù)為n,每個辨識步長內(nèi)的節(jié)點(diǎn)數(shù)記為N,則本文中算法的復(fù)雜度為O(n·N)。由于MPSP算法收斂性較快,通過幾次迭代一般即可滿足要求,因此迭代次數(shù)n通常較小。本文中的仿真案例對阻力系數(shù)初始值進(jìn)行了30%的拉偏,升力系數(shù)初始值直接賦值為零,此時迭代次數(shù)為4~6次??紤]實(shí)際情況氣動參數(shù)初值誤差一般小于文中的仿真設(shè)置,因此迭代次數(shù)會更少,算法復(fù)雜度可降至O(N),具備在線計(jì)算的潛力。
從算法耗時的角度來看,在CPU為AMD Ryzen 7 4800H@2.9Hz、內(nèi)存為16G的計(jì)算機(jī)上,采用Matlab 2019a編制程序?qū)Ρ疚闹兴惴ǖ暮臅r進(jìn)行了統(tǒng)計(jì)。結(jié)果如圖10所示。橫軸表示全彈道的辨識步長序號,總共47步。其中前10步為無控段辨識,每個步長為4 s,平均耗時約0.34 s;后37步為有控段辨識,每個步長為2 s,平均耗時約0.17 s。可見算法耗時較低。
圖10 每個步長的耗時
本文中提出的基于MPSP算法的氣動辨識方法,是一種類比帶有終端約束制導(dǎo)律設(shè)計(jì)思想的一種辨識新思路。通過在某型155mm制導(dǎo)炮彈的升阻力系數(shù)的辨識應(yīng)用中證實(shí)了其有效性,并有如下結(jié)論:
1) 在僅依賴輸入外彈道數(shù)據(jù)及彈體質(zhì)量等必要信息的前提下,本文中算法可對升、阻力系數(shù)進(jìn)行辨識,且辨識精度較高,阻力系數(shù)辨識誤差小于2%,升力系數(shù)辨識誤差小于1%。
2) 本文算法對待辨識氣動參數(shù)初值精度要求不高。在阻力系數(shù)拉偏30%、升力系數(shù)未知的情況下,通過數(shù)次迭代即可收斂至真值附近。
3) 本文辨識算法設(shè)計(jì)簡單、運(yùn)算速度快,有利于工程實(shí)現(xiàn)。