李東陽, 常思江, 王中原

(南京理工大學(xué) 能源與動力工程學(xué)院, 江蘇 南京 210094)

0 引言

對于飛行中的彈箭,在攻角δ較小的情況下,可將作用在彈箭上的氣動力和力矩視為攻角的線性函數(shù),并可采用sinδ≈δ、cosδ≈1的線性化假設(shè),據(jù)此形成線性化外彈道理論,該理論曾在很大范圍內(nèi)成功地預(yù)示了彈箭運動,促進了彈箭飛行性能設(shè)計的發(fā)展。然而,彈箭在各種條件下的實際運動,總是具有不同程度的非線性[1]。隨著彈箭外形及其飛行條件的日益復(fù)雜化、多樣化,彈箭運動的非線性特征也愈發(fā)凸顯。如在大攻角情況下,作用于彈箭的氣動力和力矩一般為攻角δ的非線性函數(shù),有的還是彈箭滾轉(zhuǎn)方位角的非線性函數(shù)(如誘導(dǎo)側(cè)向力矩),因而數(shù)學(xué)上無法利用線性理論開展研究;對于一些客觀存在的現(xiàn)象,如“行進軍艦上發(fā)射某尾翼式低速旋轉(zhuǎn)火箭,左舷發(fā)射穩(wěn)定、右舷發(fā)射不穩(wěn)定”、“某些火箭彈飛行中出現(xiàn)錐形運動而致失穩(wěn)”等,線性理論也無法給出合理的解釋。因此,彈箭非線性運動理論應(yīng)運而生并迅速發(fā)展。

20世紀(jì)60年代至本世紀(jì)初,一些學(xué)者就強非線性靜力矩[2]、非對稱再入飛行器的非線性運動[3-4]、轉(zhuǎn)速-攻角閉鎖與災(zāi)變性偏航[5-6]、彈箭非線性動力學(xué)[7-8]、有控彈非線性穩(wěn)定性[9]等問題,利用非線性數(shù)學(xué)工具開展研究,建立了彈箭非線性運動研究的理論方法體系,取得了大量研究成果;Mc Coy等[10]通過開展靶道飛行試驗及氣動系數(shù)辨識,獲得了大量的彈丸非線性氣動力數(shù)據(jù)。近20年來,隨著彈箭技術(shù)的進一步發(fā)展,彈箭非線性運動研究繼續(xù)受到廣泛關(guān)注。Morote 等就無控火箭共振閉鎖[11]、尾翼彈災(zāi)變性偏航控制[12-13]、尾翼彈非線性滾轉(zhuǎn)運動[14]、高階非線性氣動力[15]等開展了一系列深入研究;韓子鵬等[16]對彈箭非線性運動的國內(nèi)外研究成果進行了系統(tǒng)梳理,主要涉及非旋轉(zhuǎn)尾翼彈平面非線性運動、旋轉(zhuǎn)彈箭非線性運動、非線性氣動系數(shù)獲取等方面;李臣明等[17]對遠程彈箭的轉(zhuǎn)速閉鎖現(xiàn)象進行了理論;Li等[18]、Chang等[19]和李東陽等[20]重點研究了彈箭非線性運動方程的解析解,其中文獻[18]在Murphy 等[21]關(guān)于轉(zhuǎn)速-攻角閉鎖解析模型的基礎(chǔ)上,進一步完善了攻角-轉(zhuǎn)速閉鎖的解析解,文獻[19-20]則分別探討了利用同倫分析法、正規(guī)形法求取高精度攻角解析解的可行性,并據(jù)此開展了彈箭穩(wěn)定性分析、參數(shù)設(shè)計等相關(guān)應(yīng)用;此外,Cross等[22]利用數(shù)值仿真對制導(dǎo)炮彈的非線性穩(wěn)定性開展了分析;鐘揚威等[23]、楊杰等[24]利用分岔理論對彈箭非線性角運動等開展了研究。可以預(yù)見,彈箭非線性運動理論將是從根本上提高彈箭性能設(shè)計水平的重要基礎(chǔ)和有力工具。

彈箭非線性運動的穩(wěn)定性不僅與系統(tǒng)參數(shù)(如氣動參數(shù)、結(jié)構(gòu)參數(shù)等)有關(guān),而且與初始條件密切相關(guān),這是非線性系統(tǒng)與線性系統(tǒng)的不同之處。故要使彈箭達到期望的穩(wěn)定狀態(tài),必須對其初始條件加以限定。定義一個包含初始條件的緊集,若從該集合內(nèi)任一點出發(fā)的軌跡不超出該集合且最終都將收斂到系統(tǒng)的某一個平衡點,則稱該緊集為對應(yīng)穩(wěn)定平衡點的吸引域[25-26]。因此,只有當(dāng)初始運動狀態(tài)在吸引域內(nèi)時,彈箭運動才是穩(wěn)定的,這就使得吸引域的確定變得十分重要。在開展彈箭穩(wěn)定性設(shè)計時,必須通過設(shè)計相關(guān)參數(shù),確保吸引域包含彈箭的可能運動狀態(tài)范圍。然而,一般來說,除一些特殊的系統(tǒng)外,非線性系統(tǒng)準(zhǔn)確的吸引域確定相當(dāng)困難[27],一般通過解析或數(shù)值計算方法進行估計。現(xiàn)有的吸引域估計方法,主要采用李雅普諾夫理論進行問題的描述[28]。對于多項式系統(tǒng),在吸引域問題的求解上普遍采用多項式的平方和方法[25]。該方法的本質(zhì)是將多項式的正定要求松弛為平方和問題,進而利用平方和規(guī)劃方法求解[29]。平方和規(guī)劃包含了多項式半正定問題的求解、平方和優(yōu)化問題與平方和可行性問題的求解,以上問題在具體的吸引域估計問題中都會遇到[26]。盡管該方法具有一定的保守性,且其計算復(fù)雜度隨著系統(tǒng)的維數(shù)和所尋求的李雅普諾夫函數(shù)階次的增大而迅速增大,以及僅適用于多項式系統(tǒng)的局限性,但它使得吸引域的估計變得可行,與基于蒙特卡洛軌跡積分的方法相比,平方和規(guī)劃方法給出的結(jié)果具有可靠性和絕對性[30],并給出了吸引域的解析表達式,這是該方法的優(yōu)勢所在。

因此,本文擬探討平方和方法在彈箭角運動非線性吸引域估計中的應(yīng)用。對此,將針對兩種情形開展吸引域估計:一是無控尾翼彈平面非線性運動,考慮高階非線性靜力矩以及二次方非線性阻尼;二是脈沖末段修正迫彈[31]的空間非線性運動,考慮非線性法向力系數(shù)、馬格努斯力矩系數(shù)、俯仰力矩系數(shù)及赤道阻尼力矩系數(shù)等。由于平方和規(guī)劃僅適用于多項式系統(tǒng),故在角運動方程的建模中需要對幾何非線性等非多項式因素進行多項式近似。本文研究旨在為彈箭非線性運動研究提供新的有效理論工具。

1 基于平方和規(guī)劃的吸引域估計方法

本節(jié)將簡要介紹基于平方和規(guī)劃的非線性吸引域估計方法。對于一個非線性自治多項式:

(1)

式中:x為狀態(tài)變量,x∈n;t為自變量;f(x)是n×1維的多項式向量場;x0為t=0時的狀態(tài)值,即初始條件。不失一般性,設(shè)原點是系統(tǒng)的平衡點,于是f(0)=0,且原點漸近穩(wěn)定,則根據(jù)直接李雅普諾夫定理(定理1),原點的吸引域可定義為

{x0∈n: 若x(0)=x0則

(2)

定理1[25-26]若存在連續(xù)可微標(biāo)量函數(shù)V(x):n→和標(biāo)量γ∈+,使得

V(x)>0, ?x≠0∧V(0)=0

(3)

Ω:={x:V(x)≤γ}

(4)

Ω?{x:(?V(x)/?x)f(x)<0}∪{0}

(5)

則原點漸近穩(wěn)定,Ω是原點的確切吸引域式(2)的子集,可作為吸引域的一個估計,γ為其水平集。

可見,所尋找吸引域的估計,需滿足式(3)～式(5)對V(x)和γ的約束。而式(3)～式(5)中,既含有不等式約束,又含有集合包含約束。利用代數(shù)幾何學(xué),集合包含關(guān)系可以通過著名的推廣S-procedure(定理 2)[32]轉(zhuǎn)化為不等式約束。

定理2(推廣S-procedure)給出多項式g0(x),…,gm(x)∈R[x]和多項式s1(x),…,sm(x)∈Σn,R[x]代表系數(shù)為實數(shù)的多項式集合,Σn代表n個變量組成的平方和多項式集合,若

(6)

則

{x|g1(x),…,gm(x)≥0}?{x|g0(x)≥0}

(7)

因此,問題轉(zhuǎn)化為處理一系列多項式不等式。通常一個多元多項式的非負性是很難確定的,屬于非確定性多項式(NP)困難問題,解決的思路是利用平方和松弛方法[25],將多項式的非負性問題轉(zhuǎn)化為平方和問題,進而利用平方和規(guī)劃方法[25-26]對問題進行求解。例如,驗證一個多項式s(x)∈R(x)是否為正,就等價于驗證是否存在正定矩陣Q,使得

s(x)=ZT(x)QZ(x)∈Σn

(8)

式中:Z(x)為單項式向量。關(guān)于平方和方法的具體內(nèi)容,可參見文獻[26,32]。

綜上所述,吸引域估計問題可以轉(zhuǎn)化為以下平方和問題:

(9)

式中:l1(x)∈R(x)且l1(x)>0,一般取εxTx,ε為大于0的常數(shù)(如取為10-6);s1(x)輔助算子,為一定階數(shù)的平方和多項式。

為了進一步優(yōu)化吸引域的估計,可通過一個稱為形狀函數(shù)的多項式h(x)>0[26, 33],將原問題式(9)改寫為一個雙線性平方和規(guī)劃問題:

(10)

式中:φ為形狀函數(shù)h(x)的水平集,φ∈;l2與l1類似,同為多項式小量;s2與s1類似,亦為平方和輔助算子。

將上述雙線性問題解耦為兩個單線性優(yōu)化問題[26],采用迭代算法(常稱為V-s迭代算法[26-27, 33])進行求解,給定一個已知且可行的初始李雅普諾夫函數(shù)V0和形狀函數(shù)h,并令V(x)=V0,具體步驟簡述如下:

1)固定V(x),利用二分法求得γ的最大值γ*,并記下此時的s2,即

2)固定V(x)、變換γ*,利用二分法求得φ的最大值φ*,并記下此時的s1,即

3)固定s1、s2、φ*、γ*,求解新的滿足以下約束的李雅普諾夫函數(shù)V(x),這里可根據(jù)需要通過待定系數(shù)設(shè)定V(x)的階次(如2次、4次或6次李雅普諾夫函數(shù)):

-[(?V(x)/?x)f+l2]-(γ*-V(x))s2∈Σn
(γ*-V(x))-(φ*-h)s1∈Σn
V(x)-l1∈Σn
V(0)=0

4)V(x)=V(x)/γ*;

5)重復(fù)以上步驟,得到最終的李雅普諾夫函數(shù)V(x)=V(x)/γ*和相應(yīng)的吸引域估計Ω:={x∈n:V(x)<1}。停止條件可為以下任意一種:

①遇到數(shù)值計算不可行問題;

②相鄰兩次φ的最優(yōu)值小于既定容許度,如0.001;

③預(yù)定的迭代次數(shù)i。

值得注意的是,在尋找新的李雅普諾夫函數(shù)這一步(步驟3)中,s1、s2為已知且固定,與李雅普諾夫函數(shù)的階次無直接關(guān)系。

2 考慮高階氣動力的無控彈箭角運動吸引域

引起彈箭非線性運動的因素包括氣動力非線性、結(jié)構(gòu)非線性及幾何非線性等,其中又以氣動力非線性占主導(dǎo)地位[1]。Morote等[12]的研究結(jié)果表明,對于某彈箭,當(dāng)靜力矩系數(shù)關(guān)于攻角的階數(shù)達到7次時才能準(zhǔn)確地擬合試驗數(shù)據(jù);Liano等[15]也發(fā)現(xiàn),為準(zhǔn)確預(yù)測轉(zhuǎn)速閉鎖,當(dāng)攻角為12°和20°時,與滾轉(zhuǎn)角相關(guān)的非線性力矩關(guān)于攻角的階數(shù)分別不能低于 5次和7次。因此,對于彈箭非線性運動研究,有必要引入高階氣動力系數(shù)。將利用第1節(jié)中的方法,研究無控尾翼彈平面非線性角運動的吸引域估計。

2.1 無控尾翼彈非線性平面角運動方程

引入高階氣動力,考慮5次方和7次方非線性靜力矩以及2次非線性阻尼對角運動的影響。無控尾翼彈的平面非線性角運動方程可描述為

δ″+(H0+H2δ2)δ′-(M0+M2δ2+M4δ4+M6δ6)δ=0

(11)

2.2 平衡點的穩(wěn)定性

若該角運動方程的平衡點記為δ*,其滿足:

(12)

因此,在高階非線性靜力矩系數(shù)的影響下,角運動可能存在除原點以外的其他平衡點。當(dāng)不考慮阻尼項時,平衡點對應(yīng)的特征根λ滿足

λ2=M0+3M2δ*2+5M4δ*4+7M6δ*6

(13)

故原點對應(yīng)的特征根滿足λ2=M0<0,則平衡點為中心,穩(wěn)定性無法由線性理論確定。以3次非線性靜力矩為例,當(dāng)M0、M2同號時,不存在其他平衡點,理論上不論初始條件如何,攻角都將做幅值為δ0的擺動;當(dāng)M0、M2異號時,存在一對實平衡點δ*2=-M0/M2,且當(dāng)M0<0時,δ*為不穩(wěn)定鞍點。

當(dāng)考慮阻尼項時,阻尼項的存在雖然不改變角運動方程式(11)的平衡點位置,但影響其穩(wěn)定性,此時特征根滿足

λ2+(H0+H2δ*2)λ-(M0+3M2δ*2+5M4δ*4+7M6δ*6)=0

(14)

為了確定攻角運動的穩(wěn)定性,需綜合考慮系統(tǒng)平衡點式(12)和阻尼項H0、H2的影響。原點的穩(wěn)定性將由線性靜力矩系數(shù)M0和線性阻尼系數(shù)H0共同確定,即當(dāng)M0<0時,若H0<0,則原點為不穩(wěn)定平衡點,若H0>0,則原點為穩(wěn)定平衡點。此外,文獻[20]表明,考慮阻尼時,符號相反的H0、H2使得非線性的角運動系統(tǒng)出現(xiàn)了極限環(huán)。極限環(huán)的穩(wěn)定性與原點穩(wěn)定性的共同作用決定了攻角幅值的變化規(guī)律,如表1所示。詳細討論參見文獻[20]。

表1 系統(tǒng)穩(wěn)定性分析

由表1可知,對于不同符號H0、H2的組合,可能存在穩(wěn)定或不穩(wěn)定的極限環(huán)。對于原點為穩(wěn)定平衡點的系統(tǒng),設(shè)計者感興趣的是彈箭在什么初始條件下角運動可以收斂到原點。

2.3 吸引域估計結(jié)果

采用文獻[19]中的平面角運動參數(shù)開展具體的吸引域估計,其中馬赫數(shù)Ma=0.6,靜力矩系數(shù)M0=-5.0×10-5,M2=-4.5×10-4,M4=-8.0×10-3,M6=0(由于缺乏數(shù)據(jù),這里取為0值),阻尼力矩系數(shù)H0=2.0×10-3,H2=-0.4。如表1所示,對于該H0、H2的具體值,原點為穩(wěn)定平衡點,且存在一個不穩(wěn)定極限環(huán),相軌線圖如圖1所示。

圖1 不同初始條件下攻角運動相平面圖

當(dāng)角運動初始條件位于此極限環(huán)內(nèi)時,攻角運動收斂到終點;角運動初始條件在此極限環(huán)外時,攻角運動發(fā)散。因此,穩(wěn)定原點的確切吸引域由該極限環(huán)表示(見圖1中的陰影區(qū)域),該極限環(huán)可通過數(shù)值計算并逆向積分得到。需要說明的是,并非任何系統(tǒng)都可以采用逆向數(shù)值積分方法獲得其確切吸引域。如Van de Pol系統(tǒng),由于其形成了不穩(wěn)定極限環(huán),從環(huán)內(nèi)所有點出發(fā)的軌跡均收斂于環(huán)內(nèi)平衡點,故該平衡點的吸引域就是該極限環(huán)包圍的區(qū)域,進而可以通過逆向積分,使軌跡最終收斂于極限環(huán),得到確切吸引域,而角運動方程式(11)所示系統(tǒng),本質(zhì)上恰為一Van de Pol系統(tǒng)。

取x1=δ、x2=δ′,則角運動方程式(11)可表示為

(15)

利用V-s迭代算法對原點的吸引域進行估計,選取初始李雅普諾夫函數(shù)

V0=xTPx

(16)

式(16)中的P由李雅普諾夫方程ATP+PA=-I計算得到,A=(?f/?x)|x=0,I為單位矩陣;形狀函數(shù)取為h(x)=xTx;利用V-s迭代算法,分別尋找 2次、4次和6次李雅普諾夫函數(shù)對吸引域進行估計,最終得到的李雅普諾夫函數(shù)分別記為V2、V4、V6:

V2=0.031 708δ2+134.146 533δδ′+585 313 523.361 843δ′2

(17)

V4=1.341 471e-5δ4-0.090 763δ3δ′+453 519.717 597δ2δ′2+228 473 313 060.431 2δδ′3+2.277 762e16δ′4-3.917 141e-7δ3+0.085 829δ2δ′+776.309 701δδ′2+2 696 018 997.682 03δ′3+0.027 334δ2+481.738 113δδ′+481 989 806.810 080δ′2

(18)

V6=6.535 639e9δ6-0.002 437δ5δ′-20.835 500δ4δ′2+5.912 217e8δ3δ′3+4.210 891e14δ2δ′4+0.000 19δ4δ′+1.829 569e19δδ′5+2.396 202e24δ′6-3.296 192e-10δ5-32.114 797δ3δ′2+9.010 24e6δ2δ′3-4.976 216e11δδ′4-1.172 415 2e16δ′5+1.818 161e-5δ4-6.580 880δ3δ′-7.076 783e4δ2δ′2+1.654 866e11δδ′3-8.684 899e14δ′4+2.375 432e-8δ3-0.010 868δ2δ′+1 825.386 293δδ′2+37 720 324.123 641 36δ′3+0.014 637 006 762 577 65δ2+429.817 180 706 390 3δδ′+285 374 579.309 986 1δ′2

(19)

吸引域估計結(jié)果如圖2所示。

圖2 不同階次李雅普諾夫函數(shù)估計下的吸引域

由圖2可見,隨著李雅普諾夫函數(shù)階次的增加,吸引域的近似效果逐漸提升。6階李雅普諾夫函數(shù)的結(jié)果已無限接近確切吸引域。這一結(jié)果正如文獻[34]所得結(jié)論,李雅普諾夫函數(shù)的階次越高,可以越準(zhǔn)確地接近確切吸引域。主要原因在于,當(dāng)李雅普諾夫函數(shù)的階次越高,可表達成的形狀也越為多樣,故高階李雅普諾夫函數(shù)有能力表示形狀復(fù)雜的吸引域。因此,基于平方和規(guī)劃的非線性吸引域估計方法可用于分析彈箭角運動的穩(wěn)定性,能夠定量地確定角運動的初始穩(wěn)定范圍。

3 脈沖末修迫彈的非線性角運動吸引域

隨著彈箭技術(shù)的發(fā)展,彈箭非線性運動研究已不局限于無控彈箭,對于有控彈箭的非線性分析也已開展過相關(guān)研究[9, 22],但限于研究工具的缺乏,這些工作尚不深入。將以脈沖控制的末修迫彈[31]為對象,開展非線性吸引域估計,根據(jù)所得吸引域反求脈沖參數(shù)設(shè)計應(yīng)滿足的約束,據(jù)此提出脈沖控制參數(shù)的設(shè)計方法。

3.1 非線性角運動方程及其多項式近似

為便于描述,定義彈軸坐標(biāo)系(簡稱彈軸系),其原點O位于彈質(zhì)心,x軸沿彈軸指向彈體頭部為正,y軸沿水平方向,z軸在鉛錘面內(nèi)垂直于y軸指向下為正。在彈軸系內(nèi),彈體的速度和角速度分別表示為U=[uvw]T和ω=[pqr]T,u、v、w分別為彈體速度在彈軸坐標(biāo)系三軸上的分量,p、q、r分別為彈體角速度在彈軸坐標(biāo)系三軸上的分量;彈軸系相對于慣性坐標(biāo)系的角速度在彈軸系內(nèi)的投影可表示為Ωx=[pxqr]T,px表示該角速度的軸向分量。作用在彈體的橫向(即沿y軸和z軸)氣動力和氣動力矩可以表示為如下復(fù)數(shù)形式:

(20)

彈體的橫向運動在彈軸系下可表示為

(21)

(22)

不考慮馬格努斯力,氣動力和氣動力矩系數(shù)可表示為

Cy+iCz=-CNαξ

(23)

(24)

其中:CNα為法向力系數(shù);CMpα為馬格努斯力矩系數(shù);CMα為俯仰力矩系數(shù)。

考慮氣動力非線性,則有

(25)

式中:CNα0、CNα2分別為線性和三次方法向力系數(shù);CMpα0、CMpα2分別為線性和三次方馬格努斯力矩系數(shù)。

將式(23)和式(24)代入橫向運動方程組式(21),且彈軸系下角速度px≈0,并將自變量從時間t換為無量綱弧長s。為方便起見,選擇變量

(26)

則角運動方程可寫為

(27)

(28)

式中:η為幾何非線性項,η=u/U。

由幾何關(guān)系可知

(29)

則有

(30)

由氣動力和重力沿速度方向的分量可得速度方程

(31)

其中:CD為阻力系數(shù)。

由于平方和計算工具(如SOSTOOLs、SOSOPTs、SeDuMi等[26-27, 33])只能應(yīng)用于多項式,故式(27)需表示為多項式的形式,即對幾何非線性η進行多項式近似。

為此,將式(29)和式(30)泰勒展開并取1次近似,可得

(32)

式中:O(·)為高階小量表示符號。

為驗證上述多項式近似系統(tǒng)代替原系統(tǒng)進行角運動穩(wěn)定性分析的準(zhǔn)確性,首先要分析二者的平衡點和平衡點穩(wěn)定性是否一致,然后通過數(shù)值計算對比二者所得軌跡的一致性。

3.2 脈沖未作用時的吸引域估計結(jié)果

以文獻[10]中提供的某120 mm迫彈作為無控彈平臺,結(jié)構(gòu)參數(shù)如表2所示。

表2 某120 mm迫彈結(jié)構(gòu)參數(shù)[10]

由于末修迫彈是在彈道末段(如距離目標(biāo)800～1 000 m)進行修正控制,其大部分彈道為無控彈道,末段啟控點參數(shù)可作為末段運動的初始條件。假設(shè)彈體不滾轉(zhuǎn)(P=0),以初速318 m/s、射角45°對上述120 mm迫彈進行無控彈道計算,其至目標(biāo)斜距800 m處的彈道參數(shù)如表3所示,對應(yīng)的角運動方程參數(shù)如表4所示。

表3 彈道末段起始位置的彈道參數(shù)

表4 彈道末段起始位置的角運動方程參數(shù)

此時,角運動方程式(27)具有唯一平衡點x*=[0.000 2 rad 0 0 0]T。在該平衡點進行線性化,線性化系統(tǒng)特征方程的兩對特征根均具有負實部,則根據(jù)線性化理論,角運動在平衡點x*的鄰域內(nèi)可穩(wěn)定收斂到x*。狀態(tài)變量的數(shù)值積分結(jié)果如圖3所示。

圖3 原模型和多項式近似模型的相軌跡比較

如圖3所示,多項式近似系統(tǒng)和原系統(tǒng)軌跡是一致的,驗證了上述多項式近似系統(tǒng)代替原系統(tǒng)進行角運動穩(wěn)定性分析的準(zhǔn)確性。

對于彈箭的實際角運動而言,所允許的攻角和角速度是有范圍的,故需要找到具有實際意義的攻角和角速度穩(wěn)定范圍,即開展吸引域估計。

假設(shè)在選定點附近一段時間內(nèi),彈箭速度、飛行高度、彈道傾角的變化忽略不計,則氣動系數(shù)保持不變。將非線性氣動力式(25)和幾何非線性近似式(32)代入角運動方程式(27),并通過坐標(biāo)轉(zhuǎn)換將平衡點移至原點,之后可利用平方和規(guī)劃估計其吸引域。

利用V-s迭代算法,分別尋找二次和四次李雅普諾夫函數(shù),最終所得函數(shù)記為V2和V4:

V2=24.063 268α2+344.850 517αα′+24.063 267β2+344.850 514ββ′+279 665.08α′2+279 665.07β′2-0.009 189α-0.065 849α′

(33)

V4=0.291 549α4+8.517 460α3α′+0.583 178α2β2+8.467 055α2ββ′+1 204.838 384α2α′2+2 960.584 841α2β′2+8.466 879αβ2α′-3 496.704 7αβα′β′+4 404.953 287αα′3-0.000 430αα′2β′+4 305.759 8αα′β′2-0.000 508αβ′3+0.291 549β4+8.517 493β3β′+2 960.576 6β2α′2-1.018 918e-6β2α′β′+1 204.859 8β2β′2+0.000 402βα′3+4 309.150 272βα′2β′+0.000 544βα′β′2+4 404.883 9ββ′3+2 692 175.590 3α′4-8.036 145e-5α′3β′+5 401 623.866 9α′2β′2+0.000 749α′β′3+2 692 218.734 4β′4-0.000 242α3+0.018 099α2α′-0.000 243αβ2+0.028 630 1αββ′-0.074 852αα′2-6.852 724e-6αα′β′-2.955 652αβ′2-0.010 157β2α′

(34)

上述李雅普諾夫函數(shù)V2、V4等,必然滿足V(x)>0。這是由于V(x)>0為多項式正定的條件,在平方和規(guī)劃中已將其轉(zhuǎn)化成平方和約束,即V(x)-l1(x)∈Σn(體現(xiàn)在V-s算法的第3步),故V(x)>0本身即為問題的約束條件,求得的V2、V4等自然滿足該約束。計算所得吸引域分別為

(35)

其空間表達如圖4所示,其在各相平面的軌跡(剖面圖)如圖5所示。

圖4 吸引域估計結(jié)果的空間表示

圖5 吸引域估計結(jié)果的剖面圖

如圖4和圖5所示,基于平方和規(guī)劃的吸引域估計給出了很好的結(jié)果。需要說明的是,上述計算是在假設(shè)彈箭飛行環(huán)境不變的條件下得到,因而只在相對較短的一段彈道上成立。上述計算結(jié)果表明,平方和規(guī)劃方法用于確定角運動穩(wěn)定范圍是有效的。值得注意的是,由于式(26)表示的氣動力僅考慮到立方次,如果能夠建立更為準(zhǔn)確的高階氣動力模型,將得到更加準(zhǔn)確的結(jié)果。

3.3 脈沖作用下彈體穩(wěn)定的初始條件范圍

對于脈沖末修迫彈,發(fā)射后先做無控飛行,當(dāng)飛入預(yù)定區(qū)域后(彈道末端,如距離目標(biāo)800 m斜距處)啟控,進行彈道修正以提高對目標(biāo)的射擊精度。由于采用的控制執(zhí)行機構(gòu)為脈沖發(fā)動機,在設(shè)計脈沖控制參數(shù)時,需要確定脈沖大小J、脈沖作用的軸向位置Lx(距質(zhì)心,作用在質(zhì)心前為正)以及脈沖作用的周向位置φJ(rèn)(與彈軸系y軸的夾角,即脈沖方位角)。由于脈沖作用的時間極短(通常約幾毫秒至幾十毫秒),可認(rèn)為其作用前后僅引起攻角和攻角速度的變化。

設(shè)脈沖作用前(記為t=0-)攻角狀態(tài)為x0-=[α0,α′0,β0,β′0]T,脈沖作用后(記為t=0+)x0+=[α,α′,β,β′]T,則有

x0+=x0-+xJ

(36)

式中:xJ為脈沖作用引起的角運動狀態(tài)增量,

(37)

若要保證脈沖作用前后的角運動穩(wěn)定,需同時滿足:

(38)

則穩(wěn)定的角運動范圍為上述兩個集合的交集Ωx0-∩Ωx0+,設(shè)其仍具有吸引域所對應(yīng)的李雅普諾夫函數(shù)之形式,則該交集可近似為

Ω∩={x0-|V((x0-+x0+)/2)<γ∩}

(39)

式中:γ∩為函數(shù)V((x0-+x0+)/2)的水平集。因此,可以方便地利用定理2處理集合包含關(guān)系,并利用平方和規(guī)劃得到盡可能大的水平集γ∩。

考慮算例:脈沖沖量大小為J=60 N·s,作用軸向位置為Lx=0.083d(即質(zhì)心前0.083d)。將上述脈沖參數(shù)值代入式(38),并利用平方和規(guī)劃通過式(39)計算不同脈沖方位角φJ(rèn)所對應(yīng)的角運動穩(wěn)定范圍,即交集Ω∩的水平集γ∩。計算結(jié)果表明,對于不同的φJ(rèn),水平集γ∩的大小相同,但穩(wěn)定區(qū)域Ω∩的位置隨φJ(rèn)的變化而變化,如圖6所示。

圖6 交集估計Ω∩在α′-β′平面上的位置變化

在圖6所示α′-β′平面上,當(dāng)φJ(rèn)以45°為間隔從0～360°取值時,紅色曲線為不同φJ(rèn)所對應(yīng)的Ωx0+,綠色曲線為穩(wěn)定角運動初始條件的估計Ω∩,顯然,Ωx0+的位置隨著φJ(rèn)的變化而變化,進而Ω∩也隨著交集的變化而變化。

從圖6還可看出,存在一定的區(qū)域,在任意脈沖方位角下,彈箭角運動均穩(wěn)定。不妨將這樣的區(qū)域記為Ω?,且設(shè)其仍具有所對應(yīng)的李雅普諾夫函數(shù)之形式,即

Ω?={x0-|V(x0-)<γ?}

(40)

則仍然可以通過集合包含問題優(yōu)化計算得到水平集γ?的盡可能最大值。

由于上述問題沒有精確解,為驗證平方和規(guī)劃所得結(jié)果的正確性,下面引入另一種計算方法。由于二次李雅普諾夫函數(shù)給出的吸引域恰為橢球體,故γ∩和γ?也可通過一定的幾何關(guān)系得到。設(shè)二次李雅普諾夫函數(shù)給出的吸引域Ω2D的半徑(即半長軸長)為rΩ,并定義

(41)

為脈沖增量xJ對應(yīng)的半徑,則易知脈沖作用下穩(wěn)定角運動范圍Ω∩的存在條件為

max (rJk/rΩk)<1,k=1,2,3,4

(42)

方便起見,下面將max (rJk/rΩk)簡記為max (rJ/rΩ)。

對于二次李雅普諾夫函數(shù),設(shè)橢球體Ω∩和Ω?對應(yīng)的半徑分別為r∩和r?,則根據(jù)幾何關(guān)系,可得到半徑與水平集之間的關(guān)系為

γ∩=max2(r∩/rΩ),γ?=max2(r?/rΩ)

(43)

各半徑之間的關(guān)系為

(44)

r?=2r∩-rΩ

(45)

同時可得到

(46)

且

(47)

上述關(guān)系意味著,可以通過優(yōu)化計算得到γ∩后,直接由式(46)計算得到水平集γ?;另一種方法是通過已知的rJ和rΩ直接計算γ∩和γ?。

根據(jù)式(44)第1式可得,若r∩存在,也即脈沖在某個方位角φJ(rèn)下作用時,存在穩(wěn)定角運動范圍的條件為

max (rJ/rΩ)<2

(48)

根據(jù)式(44)第2式可得,若r?存在,即存在任意方位角下均能確保彈箭角運動穩(wěn)定的初始條件,那么

max (rJ/rΩ)<1

(49)

對于本節(jié)算例中的二次李雅普諾夫函數(shù)V2,其半徑為

(50)

當(dāng)脈沖沖量為J=60 N·s、作用位置Lx=0.083d時,脈沖增量對應(yīng)的半徑為

(51)

利用平方和規(guī)劃方法處理集合包含關(guān)系得到的γ∩,進而根據(jù)式(46)計算得到的γ?如表5所示,同時給出了根據(jù)式(47)直接計算的γ∩、γ?值。

表5 吸引域半徑計算結(jié)果

根據(jù)表5所示結(jié)果可見,優(yōu)化計算所得估計結(jié)果與幾何關(guān)系式(44)給出的結(jié)果非常接近,由此認(rèn)為二者互相驗證了各自的準(zhǔn)確性。

在此基礎(chǔ)上,對于已經(jīng)過驗證的平衡點吸引域Ω及其半徑rΩ,若給定脈沖大小J和作用位置Lx,根據(jù)式(48)和式(49),即可判斷此時是否存在穩(wěn)定的角運動范圍以及對任意脈沖方位角均穩(wěn)定的角運動范圍,并可根據(jù)式(47)計算出穩(wěn)定角運動范圍的水平集,進而由式(39)和式(40)確定出角運動穩(wěn)定范圍的具體位置。

3.4 應(yīng)用:脈沖參數(shù)設(shè)計

上述結(jié)論可用于脈沖參數(shù)的設(shè)計,即根據(jù)吸引域估計結(jié)果反向確定出脈沖參數(shù)的設(shè)計范圍,該范圍內(nèi)的任意脈沖作用均可使彈箭保持穩(wěn)定的角運動。

將脈沖增量對應(yīng)的吸引域半徑表達式(41)代入式(47),可得

(52)

即

(53)

根據(jù)表5,有γ?=0.068 6,代入式(53)即可得出滿足水平集γ?=0.068 6的脈沖參數(shù)J,Lx的所有可能組合(稱為有效脈沖參數(shù)組合),如圖7曲線與縱軸包圍的區(qū)域所示,縱軸為脈沖作用在彈軸上的位置。

圖7 滿足水平集γ?=0.068 6的脈沖參數(shù)組合

圖7中曲線與縱軸包圍的區(qū)域即為脈沖參數(shù)設(shè)計的有效范圍。由此可見,對于一定的脈沖軸向作用位置,均對應(yīng)一定的脈沖大小。這些脈沖參數(shù)的組合可保證,在任意脈沖方位角φJ(rèn)作用下,只要初始角運動狀態(tài)在Ω?范圍內(nèi),受控后的彈箭角運動仍可保持穩(wěn)定。顯然,上述有效脈沖參數(shù)組合,可直接用于非線性條件下的脈沖參數(shù)設(shè)計。

此外,當(dāng)γ?∈(0,1]時,有效脈沖參數(shù)組合如圖8所示。豎線左邊給出了滿足式(53)第1式的脈沖范圍,對應(yīng)顏色的兩條曲線所包圍的區(qū)域給出了滿足式(53)第2式的脈沖參數(shù)組合范圍,二者的重合區(qū)域即為γ?所對應(yīng)的有效脈沖參數(shù)組合范圍。由圖8可見,當(dāng)γ?越小,有效脈沖參數(shù)組合所對應(yīng)的范圍越大。

圖8 不同γ?下的有效脈沖參數(shù)組合

將已知的rΩ(即式(50))代入式(43),計算可得不同γ?對應(yīng)的Ω?的半徑r?,進而確定了α、β、α′、β′各自可取到的最大值,如圖9所示。由于角運動方程式(27)的對稱性,α和β的最大取值相等,α′和β′的最大取值相等。對于本節(jié)算例,當(dāng)γ?=0.068 6,α、β的最大值為3.06°,α′、β′的最大值為0.028°/cal。

圖9 γ?與Ω?中攻角和攻角速度最大值的對應(yīng)關(guān)系

4 結(jié)論

本文探討了平方和規(guī)劃方法在彈箭角運動非線性吸引域估計中的應(yīng)用。得到以下主要結(jié)論:

1)對于考慮高階氣動力系數(shù)的無控尾翼彈非線性平面角運動吸引域估計,通過構(gòu)造李雅普諾夫函數(shù),得到了高精度的吸引域估計結(jié)果,且李雅普諾夫函數(shù)的階數(shù)越高,吸引域估計值與確切值越接近,驗證了方法的可行性和有效性。

2)針對考慮了諸非線性氣動力系數(shù)的脈沖末修迫彈空間角運動吸引域估計,通過構(gòu)建高精度多項式近似模型開展平方和規(guī)劃,對于采用二次李雅普諾夫函數(shù)的情形,吸引域估計結(jié)果與直接利用橢球幾何關(guān)系所得結(jié)果一致,驗證了方法的準(zhǔn)確性。

3)由于脈沖作用主要影響彈箭的初始運動狀態(tài),故根據(jù)脈沖末修迫彈角運動的吸引域估計,可反算并確定出脈沖沖量及其軸向作用位置的有效組合范圍(滿足非線性穩(wěn)定要求),由此得到非線性條件下脈沖參數(shù)設(shè)計的有效方法。

本文研究結(jié)果表明,平方和規(guī)劃方法為彈箭角運動吸引域估計、非線性條件下的彈箭參數(shù)設(shè)計等提供了一個有力的工具,為后續(xù)更為深入的應(yīng)用研究奠定了基礎(chǔ)。