本發(fā)明涉及制導(dǎo)技術(shù)領(lǐng)域，具體涉及一種基于模型預(yù)測(cè)控制技術(shù)的多導(dǎo)彈協(xié)同作戰(zhàn)制導(dǎo)方法。

背景技術(shù)：
隨著科技的進(jìn)步，導(dǎo)彈防御系統(tǒng)正逐步完善，導(dǎo)彈的突防能力受到嚴(yán)重威脅。多枚導(dǎo)彈對(duì)目標(biāo)進(jìn)行飽和攻擊是提高導(dǎo)彈的戰(zhàn)場(chǎng)生存能力及打擊能力的重要手段，正受到越來(lái)越多的關(guān)注。因此，在考慮戰(zhàn)場(chǎng)實(shí)際情況的前提下，研究能使多枚導(dǎo)彈同時(shí)到達(dá)目標(biāo)、具有攻擊時(shí)間約束的協(xié)同制導(dǎo)方法具有重要意義。目前現(xiàn)有技術(shù)在解決在具有攻擊時(shí)間約束的條件下實(shí)現(xiàn)多枚導(dǎo)彈同時(shí)到達(dá)目標(biāo)的協(xié)同制導(dǎo)技術(shù)問(wèn)題中，會(huì)出現(xiàn)以下問(wèn)題：1、一般都會(huì)假定導(dǎo)彈的速度為常值，然而實(shí)際過(guò)程中，很難保證導(dǎo)彈速度為常值，導(dǎo)彈速度會(huì)因推力、飛行高度的變化或外界擾動(dòng)而發(fā)生變化。2、需對(duì)導(dǎo)彈的剩余飛行時(shí)間進(jìn)行估計(jì)的問(wèn)題。而在實(shí)際飛行過(guò)程中，導(dǎo)彈的剩余飛行時(shí)間也很難進(jìn)行準(zhǔn)確的估計(jì)，尤其當(dāng)導(dǎo)彈的速度時(shí)變時(shí)，剩余飛行時(shí)間的精確估計(jì)變得更加困難。3、導(dǎo)彈飛行過(guò)程中，常會(huì)受到如隨機(jī)風(fēng)等各種干擾，而現(xiàn)有技術(shù)卻多考慮的是理想情況下的協(xié)同制導(dǎo)。4、導(dǎo)彈的舵面偏轉(zhuǎn)有限，即導(dǎo)彈控制系統(tǒng)的輸入具有非線性飽和性，進(jìn)而導(dǎo)致攻角變化量的限制。而現(xiàn)有技術(shù)并未考慮。綜上所述，現(xiàn)有技術(shù)多考慮在理想(無(wú)干擾)、較簡(jiǎn)單(速度為常值)的情況下的多導(dǎo)彈協(xié)同作戰(zhàn)制導(dǎo)。

技術(shù)實(shí)現(xiàn)要素：
有鑒于此，本發(fā)明提供了一種基于模型預(yù)測(cè)控制技術(shù)的多導(dǎo)彈協(xié)同制導(dǎo)方法，能夠?qū)崿F(xiàn)對(duì)精準(zhǔn)的協(xié)同制導(dǎo)。一種基于模型預(yù)測(cè)控制技術(shù)的多導(dǎo)彈協(xié)同作戰(zhàn)制導(dǎo)方法，針對(duì)每枚導(dǎo)彈，采用如下具體步驟：步驟一、設(shè)定擊中目標(biāo)的指定攻擊時(shí)間Td，并根據(jù)指定攻擊時(shí)間獲得指定彈目距離變化率其中，r0為初始狀態(tài)下導(dǎo)彈與目標(biāo)之間的距離；步驟二、令理論彈目距離變化率與指定彈目距離變化率相等，根據(jù)彈目相對(duì)運(yùn)動(dòng)方程、角度間的幾何關(guān)系和導(dǎo)彈的質(zhì)心動(dòng)力學(xué)模型，推導(dǎo)導(dǎo)彈在指定攻擊時(shí)間命中目標(biāo)的理論飛行軌跡；步驟三、基于彈目相對(duì)運(yùn)動(dòng)方程與導(dǎo)彈的質(zhì)心動(dòng)力學(xué)方程，建立用于描述實(shí)際彈目距離變化的預(yù)測(cè)模型；步驟四、對(duì)預(yù)測(cè)模型進(jìn)行線性化、離散化和標(biāo)準(zhǔn)化處理，獲得預(yù)測(cè)模型的標(biāo)準(zhǔn)形式；步驟五、根據(jù)性能指標(biāo)標(biāo)準(zhǔn)表達(dá)式，建立反應(yīng)實(shí)際飛行軌跡對(duì)理論飛行軌跡跟蹤效果，同時(shí)又反應(yīng)控制量變化大小的性能指標(biāo)函數(shù)；步驟六、以攻角的最大幅值umax作為控制量約束，以步驟四中獲得的預(yù)測(cè)模型的標(biāo)準(zhǔn)形式作為軌跡約束，以步驟五中獲得的函數(shù)作為性能指標(biāo)函數(shù)，利用凸優(yōu)化工具，獲得一系列在滿足最小性能指標(biāo)函數(shù)值的條件下的最優(yōu)控制序列；并取最優(yōu)序列中的第一項(xiàng)代入彈目相對(duì)運(yùn)動(dòng)方程組、彈道傾角變化率表達(dá)式、導(dǎo)彈的質(zhì)心動(dòng)力學(xué)模型、導(dǎo)彈阻力表達(dá)式和導(dǎo)彈升力表達(dá)式，獲得實(shí)際彈目距離、彈道傾角、導(dǎo)彈的速度、導(dǎo)彈的速度前置角和視線角，并將獲得的實(shí)際彈目距離、彈道傾角、導(dǎo)彈的速度、導(dǎo)彈的速度前置角和視線角作為下一時(shí)刻的初值，再進(jìn)行下一時(shí)刻優(yōu)化問(wèn)題的求解；直至實(shí)際彈目距離r＜1000m時(shí)，轉(zhuǎn)入比例導(dǎo)引進(jìn)行控制。特別地，獲得理論飛行軌跡的具體方法為：第21步：建立彈目相對(duì)運(yùn)動(dòng)方程組，并根據(jù)指定彈目距離變化率表達(dá)式，獲得導(dǎo)彈速度前置角表達(dá)式；第22步：根據(jù)設(shè)定的反余弦三角函數(shù)的定義域及最大框架角的限制，對(duì)導(dǎo)彈速度前置角表達(dá)式進(jìn)行分段描述，獲得導(dǎo)彈速度前置角的分段函數(shù)表達(dá)式；第23步：采用低通濾波器對(duì)導(dǎo)彈速度前置角進(jìn)行低通濾波，獲得低通濾波器輸出的速度前置角ηd；第24步：根據(jù)彈目相對(duì)運(yùn)動(dòng)關(guān)系，獲得彈道傾角變化率表達(dá)式；第25步：建立導(dǎo)彈的質(zhì)心動(dòng)力學(xué)模型和升力表達(dá)式，獲得攻角表達(dá)式；第26步：將獲得的彈目相對(duì)運(yùn)動(dòng)方程、導(dǎo)彈速度前置角表達(dá)式、低通濾波器輸出的速度前置角表達(dá)式、彈道傾角變化率表達(dá)式、質(zhì)心動(dòng)力學(xué)模型、導(dǎo)彈的升力表達(dá)式和攻角表達(dá)式聯(lián)立，獲得理論飛行軌跡。有益效果：本發(fā)明在考慮導(dǎo)彈速度可變的前提下，通過(guò)將導(dǎo)彈的協(xié)同作戰(zhàn)問(wèn)題轉(zhuǎn)化為理論軌跡跟蹤問(wèn)題，避免了估算剩余飛行時(shí)間，提高了協(xié)同作戰(zhàn)的準(zhǔn)確性和可實(shí)現(xiàn)性。而且，采用模型預(yù)測(cè)控制技術(shù)，能夠在有外界干擾且攻角存在有限值的情況下，實(shí)現(xiàn)對(duì)理論軌跡進(jìn)行快速精確的跟蹤。附圖說(shuō)明圖1為本發(fā)明系統(tǒng)流程圖。圖2為三枚導(dǎo)彈協(xié)同攻擊目標(biāo)彈道圖。圖3為彈目相對(duì)運(yùn)動(dòng)圖。圖4導(dǎo)彈1彈目距離變化圖；圖5為導(dǎo)彈1協(xié)同段實(shí)際彈目距離與期望彈目距離差值變化圖；圖6為導(dǎo)彈1速度前置角變化圖；圖7為導(dǎo)彈1速度變化圖；圖8為導(dǎo)彈1攻角變化圖。具體實(shí)施方式下面結(jié)合附圖并舉實(shí)施例，對(duì)本發(fā)明進(jìn)行詳細(xì)描述。本發(fā)明提供了一種基于模型預(yù)測(cè)控制技術(shù)的多導(dǎo)彈協(xié)同制導(dǎo)方法；如圖2所示，其基本思想在于：設(shè)定一指定攻擊時(shí)間，通過(guò)讓理論情況下的彈目距離變化率等于指定的彈目距離變化率，使得理論情況下的導(dǎo)彈能夠在該指定攻擊時(shí)間，多枚導(dǎo)彈同時(shí)命中目標(biāo)，實(shí)現(xiàn)協(xié)同作戰(zhàn)。其中，本文中所指的彈目距離即為導(dǎo)彈到目標(biāo)之間的距離。這樣一來(lái)，就避免了現(xiàn)有技術(shù)中需要預(yù)估剩余時(shí)間的問(wèn)題。此外，實(shí)際導(dǎo)彈的飛行會(huì)受外界干擾，本發(fā)明采用模型預(yù)測(cè)控制技術(shù)，實(shí)現(xiàn)實(shí)際飛行軌跡追蹤理論飛行軌跡，從而使得導(dǎo)彈在無(wú)論是否受到外界干擾的情況下，都能在指定攻擊時(shí)間到達(dá)目標(biāo)，實(shí)現(xiàn)協(xié)同作戰(zhàn)。針對(duì)每枚導(dǎo)彈，其具體方法為：如圖1所示：步驟一、設(shè)定擊中目標(biāo)的指定攻擊時(shí)間Td，并根據(jù)指定攻擊時(shí)間獲得指定彈目距離變化率其中，r0為初始狀態(tài)下的彈目距離；其中，為實(shí)現(xiàn)各導(dǎo)彈協(xié)同，每枚導(dǎo)彈的指定攻擊時(shí)間相同。步驟二、獲得各導(dǎo)彈在指定攻擊時(shí)間命中目標(biāo)的理論飛行軌跡；如圖3所示：在目標(biāo)的速度和速度前置角已知的情況下，若想獲得導(dǎo)彈的理論飛行軌跡，首先需要獲得有關(guān)導(dǎo)彈自身以及導(dǎo)彈與目標(biāo)之間關(guān)系的相關(guān)參數(shù)表達(dá)式；再根據(jù)各相關(guān)參數(shù)表達(dá)式，聯(lián)立方程求解獲得理論彈目距離rep，即理論飛行軌跡。為此，本發(fā)明根據(jù)所需要獲得的相關(guān)參數(shù)表達(dá)式，利用現(xiàn)有理論知識(shí)，通過(guò)聯(lián)立方程組的形式，獲得最終的理論飛行軌跡表達(dá)式；具體為：第21步：建立彈目相對(duì)運(yùn)動(dòng)方程組：其中，為理論彈目距離變化率；vT為目標(biāo)T的速度；ηT為目標(biāo)T的速度前置角；vM為導(dǎo)彈M的速度；ηM為導(dǎo)彈M的速度前置角；q為視線角；第22步：為了能夠保證沿理論飛行軌跡飛行的導(dǎo)彈在指定攻擊時(shí)間到達(dá)目標(biāo)，故令理論彈目距離變化率與指定彈目距離變化率相等；聯(lián)立方程(1)和方程組(2)，獲得導(dǎo)彈M的速度前置角的表達(dá)式：令根據(jù)反余弦三角函數(shù)的定義域及最大框架角的限制，將式(3)改為分段函數(shù)：其中，η0為導(dǎo)彈初始速度的前置角，ηmax為考慮到最大框架角限制的速度前置角。第23步：由于考慮到獲得的導(dǎo)彈M的速度前置角及其導(dǎo)數(shù)不滿足連續(xù)性要求，本發(fā)明采用低通濾波器對(duì)其進(jìn)行濾波，表示為：其中，ηd為低通濾波器輸出的速度前置角，τ為濾波器時(shí)間常數(shù)。第24步：根據(jù)彈目相對(duì)運(yùn)動(dòng)中角度間的幾何關(guān)系：θ＝q-ηM(6)其中，θ為導(dǎo)彈的彈道傾角；故結(jié)合公式(5)，獲得導(dǎo)彈的彈道傾角變化率為：第25步：建立導(dǎo)彈的質(zhì)心動(dòng)力學(xué)模型：其中，m為導(dǎo)彈的質(zhì)量；α為導(dǎo)彈的攻角，P為導(dǎo)彈受到的推力；X為導(dǎo)彈受到阻力；其中，X0為零升阻力，為阻力對(duì)攻角平方的導(dǎo)數(shù)。Y為導(dǎo)彈受到的升力，g為重力加速度。由于舵偏角引起的升力較小，所以導(dǎo)彈的升力Y≈Yα·α(10)當(dāng)攻角較小時(shí)，令sinα＝α，此時(shí)結(jié)合公式(10)，公式(8)則可變?yōu)椋?![CDATA[α=mvMθ·+mgcosθP+Yα---(11)]]>根據(jù)公式(2)至公式(11)，此時(shí)，便得到關(guān)于有關(guān)導(dǎo)彈的速度前置角ηM、彈道傾角的變化率導(dǎo)彈的速度變化率和導(dǎo)彈的攻角α，還得到的關(guān)于導(dǎo)彈與目標(biāo)之間的理論彈目距離變化率以及視線角的變化率第26步：將公式(2)、公式(3)、公式(5)、公式(8)、公式(10)和公式(11)聯(lián)立，進(jìn)而獲得理論的彈目彈目距離rep，即理論飛行軌跡。步驟三：由于理論飛行軌跡能夠?qū)崿F(xiàn)在指定攻擊時(shí)間到達(dá)目標(biāo)位置，為此，只要使實(shí)際飛行軌跡跟蹤理論飛行軌跡即可實(shí)現(xiàn)協(xié)同作戰(zhàn)。故本發(fā)明采用預(yù)測(cè)控制方法，建立用于描述導(dǎo)彈實(shí)際彈目距離變化率的預(yù)測(cè)模型；實(shí)現(xiàn)實(shí)際導(dǎo)彈跟蹤理論飛行軌跡。根據(jù)公式(2)中所建立的彈目相對(duì)運(yùn)動(dòng)方程組，對(duì)公式(2)中的第一個(gè)表達(dá)式兩端求導(dǎo)；假設(shè)不考慮目標(biāo)的運(yùn)動(dòng)，即vT為零。同時(shí)，導(dǎo)彈的實(shí)際飛行與理論飛行過(guò)程中彈目相對(duì)運(yùn)動(dòng)方程形式相同，即僅需將rep替換為r即可。將公式(8)代入求導(dǎo)后的結(jié)果，得到：為此，建立實(shí)際彈目距離的預(yù)測(cè)模型：只要導(dǎo)彈實(shí)際飛行軌跡能夠跟蹤理論飛行軌跡，導(dǎo)彈即可實(shí)現(xiàn)協(xié)同作戰(zhàn)，因此，此處建立公式(13)表達(dá)的關(guān)于彈目距離的一階和二階方程。步驟四、由于模型預(yù)測(cè)一般應(yīng)用于慢時(shí)變系統(tǒng)，其每步求解優(yōu)化問(wèn)題造成的大計(jì)算量，為了能夠解決模型預(yù)測(cè)控制求解時(shí)間長(zhǎng)、難以應(yīng)用于快時(shí)變系統(tǒng)的缺點(diǎn)，本發(fā)明采用凸優(yōu)化工具CVXGEN，以解決上述問(wèn)題。為了能夠采用凸優(yōu)化工具對(duì)實(shí)際彈目距離的預(yù)測(cè)模型求解，需預(yù)先對(duì)預(yù)測(cè)模型進(jìn)行線性化、離散化和標(biāo)準(zhǔn)化處理，即進(jìn)行以下常規(guī)操作：第41步：定義狀態(tài)變量輸入u＝α，則此時(shí)，系統(tǒng)的狀態(tài)空間表達(dá)式可以寫(xiě)為：第42步：對(duì)公式(14)在理論飛行軌跡(xep,uep)處線性化處理，獲得其中，矩陣中各元素表達(dá)式為：其中，下標(biāo)“ep”表示理論參數(shù)即矩陣A、B、C的元素中所涉及到的參數(shù)vm、η、r和α為步驟二求得的理論軌跡參數(shù)。第43步：對(duì)公式(15)利用歐拉法進(jìn)行離散化處理，得到：x(k+1)＝A1x(k)+B1u(k)+C1(16)其中，k代表當(dāng)前采樣時(shí)刻，矩陣A1＝I+AΔt，B1＝BΔt，C1＝CΔt，I為2×2單位矩陣，Δt為采樣時(shí)間。第44步：將公式(16)寫(xiě)為預(yù)測(cè)控制的標(biāo)準(zhǔn)形式：x(k+i+1|k)＝A1(k+i|k)x(k+i|k)+B1(k+i|k)u(k+i|k)+C1(k+i|k)(17)其中，i＝0…P，P為模型預(yù)測(cè)控制的控制時(shí)域，x(k+i|k)代表k時(shí)刻對(duì)k+i時(shí)刻變量的預(yù)測(cè)。步驟五、在考慮到導(dǎo)彈在飛行過(guò)程中，常會(huì)受到如隨機(jī)風(fēng)等各種干擾，且導(dǎo)彈的舵面偏轉(zhuǎn)有限的情況，采用預(yù)測(cè)控制與凸優(yōu)化相結(jié)合的技術(shù)實(shí)現(xiàn)導(dǎo)彈實(shí)際彈目軌跡對(duì)理論飛行軌跡的跟蹤，從而實(shí)現(xiàn)協(xié)同作戰(zhàn)。為實(shí)現(xiàn)對(duì)理論飛行軌跡的跟蹤，故性能指標(biāo)中需要有跟蹤誤差項(xiàng)；為了避免控制量過(guò)快變化，性能指標(biāo)中還需加入控制量變化項(xiàng)。同時(shí)，根據(jù)性能指標(biāo)函數(shù)標(biāo)準(zhǔn)表達(dá)式，將性能指標(biāo)選取如下：其中，x1為導(dǎo)彈的實(shí)際彈目距離，x1ep為理論彈目距離，Q與Rd分別為相應(yīng)的權(quán)重矩陣，均為正定矩陣。性能指標(biāo)函數(shù)中的(x1(k+i|k)-xep(k+i|k))TQ(x1(k+i|k)-xep(k+i|k))代表實(shí)際飛行軌跡對(duì)理論飛行軌跡的跟蹤效果，(u(k+i|k)-u(k+i-1|k))TRd(u(k+i|k)-u(k+i-1|k))代表控制量變化的大小。步驟六、考慮到導(dǎo)彈的舵面偏轉(zhuǎn)有限導(dǎo)致攻角的限制的情況，對(duì)控制量u進(jìn)行限幅，即令|u|≤umax(19)其中，umax大于零，為控制量的幅值。根據(jù)公式(17)、(18)和(19)，將模型預(yù)測(cè)控制的標(biāo)準(zhǔn)形式寫(xiě)為：其中，s.t.代表約束條件，公式(20)的物理含義為：在給定Q與Rd之后，以x(k+i+1|k)＝A1(k+i|k)x(k+i|k)+B1(k+i|k)u(k+i|k)+C1(k+i|k)和|u|≤umax作為約束條件，采用CVXGEN工具求解式(20)所示的優(yōu)化問(wèn)題，獲得能夠滿足J為最小值的最優(yōu)控制序列[u(k),u(k+1),…,u(k+P)]T。其中，x(k+i+1|k)＝A1(k+i|k)x(k+i|k)+B1(k+i|k)u(k+i|k)+C1(k+i|k)中的A1、B1、C1矩陣根據(jù)理論飛行軌跡參數(shù)得到，表征動(dòng)態(tài)約束。將公式(17)中包含的關(guān)于x1的關(guān)系表達(dá)式代入公式(20)，起到的是理論飛行軌跡和指定攻擊時(shí)間的限制。而|u|≤umax則是對(duì)攻角的限制。獲得最優(yōu)控制序列[u(k),u(k+1),…,u(k+P)]T；為了實(shí)現(xiàn)對(duì)導(dǎo)彈實(shí)際彈目距離的控制，需得到第k+1時(shí)刻的彈目距離、彈道傾角、導(dǎo)彈的速度、導(dǎo)彈的速度前置角和視線角；為此，選取最優(yōu)控制序列[u(k),u(k+1),…,u(k+P)]T的首項(xiàng)作為輸入量，在給定k時(shí)刻的導(dǎo)彈運(yùn)動(dòng)參數(shù)即彈目距離r(k)、彈道傾角θ(k)、導(dǎo)彈的速度vM(k)、導(dǎo)彈的速度前置角ηM(k)和視線角q(k)情況下，將其代入由公式(2)、公式(8)、公式(9)和公式(10)組成的方程組中，通過(guò)利用數(shù)值積分方法，獲得第k+1時(shí)刻的彈目距離r(k+1)、彈道傾角θ(k+1)、導(dǎo)彈的速度vM(k+1)、導(dǎo)彈的速度前置角ηM(k+1)和視線角q(k+1)，并將獲得的實(shí)際彈目距離、彈道傾角、導(dǎo)彈的速度、導(dǎo)彈的速度前置角和視線角作為下一時(shí)刻的初值，再進(jìn)行下一時(shí)刻優(yōu)化問(wèn)題的求解。直到實(shí)際彈目距離r＜1000m時(shí)，轉(zhuǎn)入比例導(dǎo)引進(jìn)行控制。實(shí)施例：如圖4至8所示，假設(shè)三枚同種類型導(dǎo)彈協(xié)同攻擊一艘靜止的艦艇，艦艇的位置為(xt,yt)＝(10000m,0)。導(dǎo)彈Mi(i＝1,2,3)的初始位置(xm0,ym0)、初始速度Vm0、初始彈道傾角θm0如表1所示。導(dǎo)彈的結(jié)構(gòu)與氣動(dòng)數(shù)據(jù)如表2所示。表1導(dǎo)彈的初始運(yùn)動(dòng)參數(shù)導(dǎo)彈(xm0,ym0)/mVm0/ms-1θm0/°導(dǎo)彈1(0,0)23030導(dǎo)彈2(50,0)22025導(dǎo)彈3(100,0)21020表2導(dǎo)彈的結(jié)構(gòu)與氣動(dòng)數(shù)據(jù)表2中，S、L分別為導(dǎo)彈的特征面積與特征長(zhǎng)度，cx0、分別為零升阻力系數(shù)、阻力系數(shù)對(duì)攻角平方的導(dǎo)數(shù)與升力系數(shù)對(duì)攻角的導(dǎo)數(shù)。導(dǎo)彈的指定攻擊時(shí)間Td＝60s，低通濾波器時(shí)間常數(shù)τ＝0.1s，模型預(yù)測(cè)控制的控制時(shí)域P＝30，性能指標(biāo)中的權(quán)重Q＝10，Rd＝0.1，控制量幅值umax＝10°。為驗(yàn)證控制系統(tǒng)的抗干擾能力，假設(shè)導(dǎo)彈飛行過(guò)程中氣動(dòng)參數(shù)有20％的攝動(dòng)；為驗(yàn)證控制系統(tǒng)對(duì)輸入的限幅作用，假設(shè)導(dǎo)彈1由于發(fā)射擾動(dòng)，初速變?yōu)?35m/s。三枚導(dǎo)彈協(xié)同作戰(zhàn)時(shí)的彈道及其它運(yùn)動(dòng)參數(shù)變化如圖4至8所示。綜上所述，以上僅為本發(fā)明的較佳實(shí)施例而已，并非用于限定本發(fā)明的保護(hù)范圍。凡在本發(fā)明的精神和原則之內(nèi)，所作的任何修改、等同替換、改進(jìn)等，均應(yīng)包含在本發(fā)明的保護(hù)范圍之內(nèi)。