本發(fā)明涉及航天器攔截交會(huì)控制領(lǐng)域，具體涉及一種快速確定能量最優(yōu)攔截預(yù)測(cè)命中點(diǎn)的數(shù)值優(yōu)化方法。
背景技術(shù)：
：在一般的軌道攔截問(wèn)題中，需要解決的一個(gè)基本問(wèn)題是在預(yù)定時(shí)間內(nèi)，如何從空間中的一點(diǎn)到達(dá)另外一點(diǎn)，這類航天動(dòng)力學(xué)中的兩點(diǎn)邊值問(wèn)題被稱為L(zhǎng)ambert問(wèn)題。數(shù)學(xué)家Gauss給出了Lambert問(wèn)題的經(jīng)典求解，如圖1所示，Gauss問(wèn)題的定義為：給定初始位置矢量r1、終端位置矢量r2，以及從r1到r2的飛行時(shí)間tF和運(yùn)動(dòng)方向，求初始和終端的速度矢量ν1和ν2。在Kepler軌道運(yùn)動(dòng)的前提條件下，位置矢量r1、r2，和速度矢量ν1、ν2共面，所以有r2=fr1+gv1v2=f·r1+g·v1---(1)]]>其中g(shù)=r1r2sinΔνμP=t-a3μ(ΔE-sinΔE)---(3)]]>f·=μPtanΔν2(1-cosΔνP-1r1-1r2)=-μar1r2sinΔE---(4)]]>g·=1-r1P(1-cosΔν)=1-ar2(1-cosΔE)---(5)]]>對(duì)于攔截問(wèn)題，一般知道初始點(diǎn)需要多大的速度，可以使飛行器沿著一條Kepler軌道，滑行至終點(diǎn)。由公式(1)，可以將初始速度表示為v1=r2-fr1g---(6)]]>所以，Gauss問(wèn)題的求解可以簡(jiǎn)化為系數(shù)f、g，以及和的求解。受超越方程的限制，必須用逐次逼近法求解，可將Gauss問(wèn)題的一般解法概括如下：1、給定P、a或ΔE中任意一個(gè)未知數(shù)的初值；2、由公式(2)和(4)計(jì)算另兩個(gè)未知數(shù)的值；3、由公式(3)求解時(shí)間t，并與給定的飛行時(shí)間tF比較，檢驗(yàn)計(jì)算結(jié)果；4、若計(jì)算結(jié)果t與給定時(shí)間tF不一致，修正迭代變量，并重復(fù)步驟2與步驟3。對(duì)于目標(biāo)改變了運(yùn)動(dòng)軌道時(shí)，確定能量最優(yōu)預(yù)測(cè)命中點(diǎn)的數(shù)值優(yōu)化方法問(wèn)題，也就是攔截變軌最小速度修正問(wèn)題，如圖2所示，t1表示動(dòng)能攔截武器E主動(dòng)段關(guān)機(jī)點(diǎn)的時(shí)刻，M表示機(jī)動(dòng)彈頭，T表示打擊目標(biāo)在地面的位置。P1表示如果機(jī)動(dòng)彈頭在t1時(shí)刻不進(jìn)行變軌，則攔截彈會(huì)與彈頭在該點(diǎn)相撞，即P1表示原預(yù)測(cè)命中點(diǎn)。如果彈頭在t1時(shí)刻按一定的突防策略進(jìn)行軌道機(jī)動(dòng)，假設(shè)目標(biāo)M從實(shí)線表示的軌道orb1改變到點(diǎn)劃線表示的軌道orb2，打擊目標(biāo)也隨之修改至T′，則當(dāng)前時(shí)刻t1的攔截器E無(wú)法以當(dāng)前速度v0，在新的目標(biāo)軌道orb2上攔截目標(biāo)M。假設(shè)新的命中點(diǎn)位置是P2，則飛行器的初始速度必須修正為v1，v1可以通過(guò)求解Gauss問(wèn)題獲得。不同于以往的攔截交會(huì)問(wèn)題，當(dāng)目標(biāo)存在軌道機(jī)動(dòng)時(shí)，攔截器選擇不同的預(yù)測(cè)命中點(diǎn)P2將帶來(lái)不同的速度修正Δv＝v1-v0。顯然，對(duì)于燃料有限的攔截器而言，速度修正量Δv＝|Δv|最小的機(jī)動(dòng)策略才是最優(yōu)的，此時(shí)的預(yù)測(cè)命中點(diǎn)P2與飛行時(shí)間tF也分別是最優(yōu)攔截位置與攔截時(shí)間。對(duì)于上述確定能量最優(yōu)預(yù)測(cè)命中點(diǎn)的數(shù)值優(yōu)化方法問(wèn)題，傳統(tǒng)的求解方法是對(duì)時(shí)間t進(jìn)行迭代搜索。P迭代法是假定一個(gè)P的試探值，并由此值計(jì)算出另外兩個(gè)未知數(shù)a和ΔE，然后解出t，并把它與給定的飛行時(shí)間tF比較，以此檢驗(yàn)試探值是否合適。P迭代法可以得出t隨P變化曲線斜率的解析表達(dá)式，因此可用牛頓迭代法修正P的試探值，提高迭代速度。定義三個(gè)常量由Kepler軌道理論，a可以表示為a=mkp(2m-l2)P2+2klP-k2---(8)]]>進(jìn)一步地，由公式(2)、(3)、(4)，可以求出時(shí)間tt=g+a3μ{arccos[1-r1a(1-f)]+r1r2f·μa}---(9)]]>另外，t隨P變化曲線斜率的解析表達(dá)式為dtdP=-g2P-32a(t-g)(k2+(2m-l2)P2mkP2)+a3μ2ksinΔEP(k-lP)---(10)]]>所以，利用P迭代法求解Gauss問(wèn)題的步驟可以表述為1、由公式(7)，計(jì)算常數(shù)k，l和m；2、給定自變量P的試探值；3、由公式(8)，計(jì)算a的值；4、由公式(2)、(3)、(4)，計(jì)算f、g和的值；5、由公式(9)，計(jì)算時(shí)間t，并與給定的飛行時(shí)間tF比較；如果相等，則結(jié)束；不相等，轉(zhuǎn)步驟6；6、結(jié)合公式(10)，并使用牛頓迭代法計(jì)算自變量P的新試探值，并轉(zhuǎn)步驟3。本發(fā)明所要探討的確定能量最優(yōu)預(yù)測(cè)命中點(diǎn)的數(shù)值優(yōu)化方法問(wèn)題，其運(yùn)動(dòng)模型可以簡(jiǎn)化為圖3中O代表地心；r1是攔截彈的初始位置矢量，r2是攔截彈與目標(biāo)的碰撞點(diǎn)，即攔截彈的終端位置矢量；v0是攔截彈機(jī)動(dòng)前的速度矢量，用虛線表示；v1是通過(guò)求解Gauss問(wèn)題所求得的初始速度，用實(shí)線表示；r2與v1的上標(biāo)代表不同的碰撞點(diǎn)以及與之對(duì)應(yīng)的攔截彈初始速度。在已知目標(biāo)軌道信息的條件下，目標(biāo)軌道的升交點(diǎn)赤經(jīng)Ω、軌道傾角i、近地點(diǎn)幅角ω、半長(zhǎng)軸aM和偏心率eM，以及通過(guò)近地點(diǎn)時(shí)刻tpM可以整理為軌道根數(shù)的形式[Ω,i,aM,eM,ω,tpM]終端位置矢量r2受給定飛行時(shí)間tF限制。以下，在不引起誤會(huì)的條件下，將給定飛行時(shí)間簡(jiǎn)稱為攔截時(shí)間t。其中，目標(biāo)運(yùn)動(dòng)的偏近點(diǎn)角EM和真近點(diǎn)角fM隨攔截時(shí)間t的迭代函數(shù)為EM=2·arctan[1-eM1+eMtan(fM2)]EM-eMsin(EM)=2πTM(t-tpM)---(11)]]>其中，目標(biāo)的軌道周期為TM。所以，終端位置矢量r2可以表示為r2=aM(1-eM2)1+eMcos(fM)]]>r2=cos(fM+w)cosΩ-sin(fM+w)cosisinΩcos(fM+w)sinΩ+sin(fM+w)cosicosΩsin(fM+w)sini·r2---(12)]]>至此，已經(jīng)可以通過(guò)傳統(tǒng)方法求解目標(biāo)存在軌道機(jī)動(dòng)的攔截變軌最小速度修正問(wèn)題。傳統(tǒng)方法需要設(shè)計(jì)兩重迭代：內(nèi)層迭代可以通過(guò)P迭代法求解Gauss問(wèn)題，獲得某一預(yù)測(cè)命中點(diǎn)P2所對(duì)應(yīng)的速度修正量Δv；外層迭代則是通過(guò)計(jì)算不同P2對(duì)應(yīng)的Δv，從而確定最小速度修正量Δvmin。計(jì)算流程圖如圖4所示，其中外層迭代使用半分法加速收斂。雖然這類搜索算法思路簡(jiǎn)單，計(jì)算精度也能滿足工程實(shí)際，但是耗時(shí)較長(zhǎng)，即使使用半分法加速搜索，也不利于在線求解。技術(shù)實(shí)現(xiàn)要素：為了解決上述技術(shù)問(wèn)題，本發(fā)明提供一種快速確定能量最優(yōu)攔截預(yù)測(cè)命中點(diǎn)的數(shù)值優(yōu)化方法，其基于最優(yōu)規(guī)劃理論，通過(guò)對(duì)攔截變軌最小速度修正問(wèn)題進(jìn)行數(shù)學(xué)建模，經(jīng)過(guò)大量的公式推導(dǎo)，獲得了Δvmin對(duì)應(yīng)的KKT條件，并通過(guò)設(shè)計(jì)牛頓迭代算法，獲得了解析的梯度信息，實(shí)現(xiàn)了對(duì)Δvmin的快速求解，同時(shí)也求解出最優(yōu)攔截位置與攔截時(shí)間。該方法計(jì)算精度高，速度快，且對(duì)初值不敏感。本發(fā)明采用以下的技術(shù)方案：一種快速確定能量最優(yōu)攔截預(yù)測(cè)命中點(diǎn)的數(shù)值優(yōu)化方法，具體步驟如下：(1)、給定自變量X的試探值，并設(shè)定系數(shù)矩陣u和v；(2)、由公式(22)和公式(24)計(jì)算偏導(dǎo)數(shù)矩陣FF′(X)；(3)、由公式(25)計(jì)算新的自變量取值Xk+1；(4)、判斷是否收斂；如果ΔXX＝XXk+1-XXk小于誤差限，計(jì)算結(jié)束；否則轉(zhuǎn)至步驟(2)，重復(fù)進(jìn)行步驟(2)-(4)；其中，步驟(1)中，所述自變量X＝[tPfM]T，引入兩組等式約束f1=g+a3μ{arccos[1-r1a(1-f)]+r1r2f·μa}-t]]>f2=EM(fM)-eMsin[EM(fM)]-2πTM(t-tpM)]]>使攔截彈的速度修正量最小，描述為數(shù)學(xué)表達(dá)式，則目標(biāo)函數(shù)為J(X)＝(v1x-v0x)2+(v1y-v0y)2+(v1z-v0z)2；至此，把目標(biāo)存在軌道機(jī)動(dòng)的攔截變軌最小速度修正問(wèn)題，轉(zhuǎn)化為一個(gè)含有等式約束的最優(yōu)規(guī)劃問(wèn)題，該最優(yōu)規(guī)劃問(wèn)題的表達(dá)式整理如下：X＝[tPfM]TJ(X)＝(v1x-v0x)2+(v1y-v0y)2+(v1z-v0z)2f1(t,P,fM)=g+a3μ{arccos[1-r1a(1-f)]+r1r2f·μa}-t---(13)]]>f2(t,fM)=EM(fM)-eMsin[EM(fM)]-2πTM(t-tpM)]]>設(shè)帶有Lagrange乘子的增廣目標(biāo)函數(shù)為L(zhǎng)＝J+λ1f1+λ2f2(14)由KKT條件，可得∂L∂t=-λ1-2πTMλ2=0---(15)]]>∂L∂P=∂J∂P+λ1∂f1∂P=0---(16)]]>∂L∂fM=∂J∂fM+λ1∂f1∂fM+λ2∂f2∂fM=0---(17)]]>由公式(15)與(16)，解得λ1=-∂J/∂P∂f1/∂Pλ2=-TM2πλ1=TM2π∂J/∂P∂f1/∂P---(18)]]>將公式(18)代入公式(17)，得到等式f3(P,fM)=∂J∂fM-∂J/∂P∂f1/∂P∂f1∂fM+TM2π∂J/∂P∂f1/∂P∂f2∂fM---(19)]]>所以，KKT條件最終轉(zhuǎn)化為三個(gè)等式約束f1(t,P,fM)=0f2(t,fM)=0f3(P,fM)=0---(20)]]>設(shè)F(X)＝[f1f2f3]T(21)則方程組(21)的解就是攔截變軌最小速度修正問(wèn)題的最優(yōu)解；步驟(2)中，所述公式(22)為方程組F(X)的Jacobi矩陣具體為：F′(X)=∂f1∂t∂f1∂P∂f1∂fM∂f2∂t0∂f2∂fM0∂f3∂P∂f3∂fM---(22)]]>其中，因?yàn)閒2中不顯含P，f3中不顯含t，所以則自變量X的迭代公式可以寫為：XK+1＝XK-[F′(XK)]-1·F(XK)(23)為提高數(shù)值計(jì)算的穩(wěn)定性，公式(22)可以轉(zhuǎn)化為公式(24)，所述公式(24)為：其中，對(duì)角陣u和v為步驟(1)所述的系數(shù)矩陣；所述公式(25)為：如上所述的方法，優(yōu)選地，所述公式(22)中，各項(xiàng)的解析表達(dá)式如下：因?yàn)椴⑶遥@然∂f1∂t=-1∂f2∂t=-2πTM---(27)]]>所以，只需獲得所述J、f1、f2對(duì)P、fM的一階偏導(dǎo)，二階純偏導(dǎo)，以及二階混合偏導(dǎo)，即可獲得公式(22)的解析表達(dá)式；將矢量投影在大地直角坐標(biāo)系，并約定h2xh2yh2z=cos(fM+w)cosΩ-sin(fM+w)cosisinΩcos(fM+w)sinΩ+sin(fM+w)cosicosΩsin(fM+w)sini---(28)]]>則首先，給出J對(duì)P、fM的一階偏導(dǎo)，二階純偏導(dǎo)，以及二階混合偏導(dǎo)，設(shè)JPy1=r2y(fM)g(P,fM)-r1yf(P,fM)g(P,fM)-v0y;JPy2=r2y(fM)-r1y[2-f(P,fM)]2·P·g(P,fM);]]>JPz1=r2z(fM)g(P,fM)-r1zf(P,fM)g(P,fM)-v0z;JPz2=r2z(fM)-r1z[2-f(P,fM)]2·P·g(P,fM);]]>r1h2＝r1xh2x+r1yh2y+r1zh2z；則因?yàn)?part;h2x∂fM∂h2y∂fM∂h2z∂fM=-sin(fM+w)cosΩ-cos(fM+w)cosisinΩ-sin(fM+w)sinΩ+cos(fM+w)cosicosΩcos(fM+w)sini]]>∂g∂fM=1μP·[∂r2∂fM·sr12h2-r2·r1ph2·r1h2sr12h2]]]>∂∂fMf(P,fM)g(P,fM)=∂f∂fM·g-f·∂g∂fMg2]]>所以，設(shè)∂∂fMr2x(fM)g(P,fM)=(∂r2∂fM·h2x+r2·∂h2x∂fM)g-r2·h2x·∂g(P,fM)∂fMg2(P,fM);JfM×2=∂∂fMr2x(fM)g(P,fM)-r1x·∂∂fMf(P,fM)g(P,fM)]]>∂∂fMr2y(fM)g(P,fM)=(∂r2∂fM·h2y+r2·∂h2y∂fM)g-r2·h2y·∂g(P,fM)∂fMg2(P,fM);JfMy2=∂∂fMr2y(fM)g(P,fM)-r1y·f(P,fM)g(P,fM)]]>∂∂fMr2z(fM)g(P,fM)=(∂r2∂fM·h2z+r2·∂h2z∂fM)g-r2·h2z·∂g(P,fM)∂fMg2(P,fM);JfMz2=∂∂fMr2z(fM)g(P,fM)-r1z·∂∂fMf(P,fM)g(P,fM)]]>則∂2J∂P2=-12P∂J∂P+2·(JPx22+JPy22+JPz22)+1Pg∂f∂P·(JPx1·r1x+JPy1·r1y+JPz1·r1z)---(32)]]>因?yàn)樗砸驗(yàn)?part;2r2∂fM2=2e2-e2cos2(fM)+ecos(fM)[1+ecos(fM)]2r2]]>∂2h2x∂fM2∂2h2y∂fM2∂2h2z∂fM2T=-h2xh2yh2zT]]>∂2g∂fM2=1μP[∂2r2∂fM2sr12h2-2·∂r2∂fMr1ph2·r1h2sr12h2-r2·(r1ph22-r1h22)·r12+r1h24sr12h23]]]>∂2f∂fM2=1r1P[2·∂r2∂fM·r1ph2-r2·r1h2-∂2r2∂fM2·(r1-r1h2)]]]>∂2f(P,fM)∂fM2g(P,fM)=∂2f∂fM2g-f∂2g∂fM2-2∂f∂fM∂g∂fM+2fg(∂g∂fM)2g2]]>所以，設(shè)∂2∂fM2r2·h2xg(P,fM)=(∂2r2∂fM2·h2x+2∂r2∂fM∂h2x∂fM-r2h2x)·g-r2·h2x·∂2g∂fM2-∂∂fMr2·h2xg(P,fM)·2g·∂g∂fMg2(P,fM)]]>∂2∂fM2r2·h2yg(P,fM)=(∂2r2∂fM2·h2y+2∂r2∂fM∂h2y∂fM-r2h2y)·g-r2·h2y·∂2g∂fM2-∂∂fMr2·h2yg(P,fM)·2g·∂g∂fMg2(P,fM)]]>∂2∂fM2r2·h2zg(P,fM)=(∂2r2∂fM2·h2z+2∂r2∂fM∂h2z∂fM-r2h2z)·g-r2·h2z·∂2g∂fM2-∂∂fMr2·h2zg(P,fM)·2g·∂g∂fMg2(P,fM)]]>則其次，給出f1對(duì)P、fM的一階偏導(dǎo)，二階純偏導(dǎo)，以及二階混合偏導(dǎo)；因?yàn)棣?＝(2m-l2)P2+2klP-k2；∂a∂P=mk[(l2-2m)P2-k2]χ02;∂2a∂P2=2mk(2m-l2)2P3+3(2m-l2)Pk2+2k3lχ03]]>并且，設(shè)χ1=1-r1a(1-f);χ2=1-χ12;χ3=r1r2tanΔν2·1-cosΔνP;]]>χ4=r1r2tanΔν2(1-cosΔνP-1r1-1r2);∂2∂P21aP=4lm1P3-6km1P4;]]>φ1=arccos[1-r1a(1-f(P,fM))]+χ4aP;]]>φ2=r1r2(1-cosΔν)χ2·∂∂P1aP-χ3P·1ap+χ4·aP2·∂∂P1aP;]]>∂φ2∂P=k·∂2∂P21aP·[2-kaP]-(∂∂P1aP)2·[aP-k][kaP]12·[2-kaP]32+aP2·∂2∂P21aP+χ3·(1aPa·(km1P-lm)P2+12aP·∂∂P1aP)+χ4·141aP(P∂a∂P+a)·∂∂P1aP]]>所以設(shè)χ51＝χ61·2P2+χ62·2·P-2·r2·χ6·χ63；χ6＝r1-r1h2；χ61=(r1h2-r2)·∂r2∂fM+r2·r1ph2;χ63=χ6·∂r2∂fM-r2·r1ph2;]]>χ64=mkχ02;χ65=χ63-r2·χ6a·∂a∂fM;]]>χ7=r1-r1h2r1+r1h2;∂χ7∂fM=-r1·r1ph2(r1+r1h2)·sr12h2;]]>χ8=(r1-r1h2)·r2P-r2-r1;∂χ8∂fM=χ63P-∂r2∂fM;∂χ8∂P=-(r1-r1h2)·r2P2;]]>則∂f1∂fM=∂g∂fM+32aμ∂a∂fMarccos[1-r1a(P,fM)(1-f(P,fM))]+aμχ65Pχ2+∂a∂fM·χ7·χ8+a·∂χ7∂fMχ8+a·χ7·∂χ8∂fMμP---(36)]]>∂2f1∂P2=34gP2+32[121aμ(∂a∂P)2+aμ∂2a∂P2]·φ1+3aμ∂a∂P·φ2+a3μ·∂φ2∂P---(37)]]>因?yàn)樗砸驗(yàn)?chi;9=∂2r2∂fM2χ6-2∂r2∂fMr1ph2+r2r1h2;χ10=1Pχ65χ2;∂χ2∂fM=χ1χ2·χ65aP]]>∂χ5∂fM=[2(∂r2∂fM)2+2r2∂2r2∂fM2]·sr12h22-8r2∂r2∂fMr1h2r1ph2+2r22(r1h22-r1ph22)]]>∂χ51∂fM[2r1ph2∂r2∂fM+(r1h2-r2)·∂2r2∂fM2-r2·r1h2-(∂r2∂fM)2]·2P2-2χ632-2r2χ6χ9+{[2(∂r2∂fM)2+(r1+2r2)·∂2r2∂fM2]·χ6-2(r1+2r2)∂r2∂fMr1ph2+(r1r2+r22)r1h2}·2P]]>∂2a∂fM2=P·(∂χ5∂fM·1χ0-χ64·∂χ51∂fM-2·χ5χ51χ02+2χ64χ512χ0)]]>∂χ10∂fM=1P·{∂2r2∂fM2·χ6-2∂r2∂fMr1ph2+r2r1h2-[r2a∂2a∂fM2+1a∂a∂fM∂r2∂fM-r2a2(∂a∂fM)2]·χ6+r2a∂a∂fMr1ph2}·χ2-χ65·∂χ2∂fMχ22]]>∂2χ7∂fM2=r1r1h2·[r1+(r1xh2x+r1yh2y+r1zh2z)]sr12h2+r1·r1ph2·{r1ph2sr12h2-[r1+(r1xh2x+r1yh2y+r1zh2z)]r1h2r1ph2sr12h2}[r1+(r1xh2x+r1yh2y+r1zh2z)]2sr12h22]]>∂2χ8∂fM2=χ9P-∂2r2∂fM2]]>所以∂2f1∂fM2=∂2g∂fM2+12a∂a∂fM{32aμ∂a∂fMarccos[1-r1a(P,fM)(1-f(P,fM))]+aμχ65Pχ2}+aμ{32∂2a∂fM2arccos[1-r1a(P,fM)(1-f(P,fM))]-32∂a∂fMχ10a+∂χ10∂fM}+∂2a∂fM2χ7χ8+a∂2χ7∂fM2χ8+aχ7∂2χ8∂fM2+2∂a∂fM∂χ7∂fMχ8+2∂a∂fMχ7∂χ8∂fM+2a∂χ7∂fM∂χ8∂fMμP---(39)]]>最后，給出f2對(duì)fM的一階偏導(dǎo)，二階偏導(dǎo)；因?yàn)?part;E∂fM=21+eM1-eM[1+cos(fM)]+1-eM1+eM[1-cos(fM)]]]>∂2E∂fM2=2{1+eM1-eM[1+cos(fM)]+1-eM1+eM[1-cos(fM)]}sin(fM)·(1+eM1-eM-1-eM1+eM)]]>所以∂2f2∂fM2=∂2E∂fM2[1-eMcos(E)]+eMsin(E)·(∂E∂fM)2---(41)]]>至此，由公式(30)至公式(41)，并結(jié)合公式(26)和公式(27)，獲得公式(22)的解析表達(dá)式，算法完整執(zhí)行。本發(fā)明以Kepler軌道理論以及基于P迭代法的Gauss問(wèn)題為基礎(chǔ)，將目標(biāo)存在軌道機(jī)動(dòng)的攔截變軌最小速度修正問(wèn)題，轉(zhuǎn)化為最優(yōu)規(guī)劃問(wèn)題，并通過(guò)最優(yōu)規(guī)劃理論獲得了與Δvmin對(duì)應(yīng)的KKT條件，同時(shí)為快速求解KKT條件，進(jìn)一步地設(shè)計(jì)了牛頓迭代，獲得了解析的梯度信息。相對(duì)于一般的搜索算法，該方法計(jì)算精度高，速度快，且對(duì)初值不敏感。通過(guò)對(duì)比試驗(yàn)，采用一種基于半分法的加速搜索算法，雖然該基于半分法的加速搜索算法相對(duì)于一般的遍歷式搜索算法，其計(jì)算效率得到大幅提高，但依然是本發(fā)明耗時(shí)的30倍。本發(fā)明也和MATLAB自帶的優(yōu)化函數(shù)fmincon進(jìn)行了對(duì)比，結(jié)果也表明無(wú)論是計(jì)算精度還是計(jì)算速度，本文的算法都更具有優(yōu)勢(shì)。附圖說(shuō)明圖1是Gauss問(wèn)題示意圖。圖2是目標(biāo)有軌道機(jī)動(dòng)的，基于最小速度修正量的攔截器交戰(zhàn)示意圖。圖3是簡(jiǎn)化的目標(biāo)存在軌道機(jī)動(dòng)的攔截變軌最小速度修正問(wèn)題示意圖。圖4是采用半分法的傳統(tǒng)求解流程圖。圖5是本發(fā)明中方法的計(jì)算流程圖。具體實(shí)施方式需要說(shuō)明的是，在不沖突的情況下，本申請(qǐng)中的實(shí)施例及實(shí)施例中的特征可以相互組合。下面請(qǐng)參考附圖并結(jié)合實(shí)施例來(lái)詳細(xì)說(shuō)明本發(fā)明。對(duì)于求解目標(biāo)存在軌道機(jī)動(dòng)的快速確定能量最優(yōu)攔截預(yù)測(cè)命中點(diǎn)的數(shù)值優(yōu)化方法，也就是攔截變軌最小速度修正問(wèn)題，一定存在一個(gè)最優(yōu)的攔截時(shí)間t，使得目標(biāo)與攔截彈經(jīng)過(guò)時(shí)間t后，同時(shí)抵達(dá)最優(yōu)攔截位置r2，并且使攔截彈在初始位置r1處的速度改變量最小。所以，本質(zhì)上攔截時(shí)間t是最小速度修正問(wèn)題的唯一自變量，但是公式推導(dǎo)過(guò)于復(fù)雜。在公式(9)與公式(11)，其中的變量P與變量fM，分別是求解Gauss問(wèn)題與求解終端位置矢量r2過(guò)程中的重要迭代變量。將最小速度修正問(wèn)題的自變量擴(kuò)展至X＝[tPfM]T；同時(shí)結(jié)合公式(9)與公式(11)，引入兩組等式約束f1=g+a3μ{arccos[1-r1a(1-f)]+r1r2f·μa}-t]]>f2=EM(fM)-eMsin[EM(fM)]-2πTM(t-tpM)]]>使攔截彈的速度修正量最小，描述為數(shù)學(xué)表達(dá)式，則目標(biāo)函數(shù)為J(X)＝(v1x-v0x)2+(v1y-v0y)2+(v1z-v0z)2至此，把目標(biāo)存在軌道機(jī)動(dòng)的攔截變軌最小速度修正問(wèn)題，轉(zhuǎn)化為一個(gè)含有等式約束的最優(yōu)規(guī)劃問(wèn)題，該最優(yōu)規(guī)劃問(wèn)題的表達(dá)式整理如下X＝[tPfM]TJ(X)＝(v1x-v0x)2+(v1y-v0y)2+(v1z-v0z)2f1(t,P,fM)=g+a3μ{arccos[1-r1a(1-f)]+r1r2f·μa}-t---(13)]]>f2(t,fM)=EM(fM)-eMsin[EM(fM)]-2πTM(t-tpM)]]>設(shè)帶有Lagrange乘子的增廣目標(biāo)函數(shù)為L(zhǎng)＝J+λ1f1+λ2f2(14)由KKT條件，可得∂L∂t=λ1-2πTMλ2=0---(15)]]>∂L∂P=∂J∂P+λ1∂f1∂P=0---(16)]]>∂L∂fM=∂J∂fM+λ1∂f1∂fM+λ2∂f2∂fM=0---(17)]]>由公式(15)與(16)，可以解得λ1=-∂J/∂P∂f1/∂Pλ2=-TM2πλ1=TM2π∂J/∂P∂f1/∂P---(18)]]>將公式(18)代入公式(17)，可以得到等式f3(P,fM)=∂J∂fM-∂J/∂P∂f1/∂P∂f1∂fM+TM2π∂J/∂P∂f1/∂P∂f2∂fM---(19)]]>所以，KKT條件最終轉(zhuǎn)化為三個(gè)等式約束f1(t,P,fM)=0f2(t,fM)=0f3(P,fM)=0---(20)]]>設(shè)F(X)＝[f1f2f3]T(21)則方程組(21)的解就是攔截變軌最小速度修正問(wèn)題的最優(yōu)解。以下介紹如何設(shè)計(jì)牛頓迭代求解方程組(21)。方程組F(X)的Jacobi矩陣是F′(X)=∂f1∂t∂f1∂P∂f1∂fM∂f2∂t0∂f2∂fM0∂f3∂P∂f3∂fM---(22)]]>其中，因?yàn)閒2中不顯含P，f3中不顯含t，所以則自變量X的迭代公式可以寫為XK+1＝XK-[F′(XK)]-1·F(XK)(23)為提高數(shù)值計(jì)算的穩(wěn)定性，可以利用系數(shù)矩陣u和v，將公式(22)轉(zhuǎn)化為公式(24)FF′(X)=abc·F′(X)·def=u·F′(X)·v---(24)]]>其中，對(duì)角陣u和v，分別稱為F(X)量級(jí)統(tǒng)一的系數(shù)矩陣，以及變量X歸一化的系數(shù)矩陣。進(jìn)一步地，迭代公式(23)可以轉(zhuǎn)化為XXK=v-1·XKXXK+1=XXK-[FF′(XK)]-1·u·F(XK)---(25)]]>至此，可以將本發(fā)明的快速確定能量最優(yōu)攔截預(yù)測(cè)命中點(diǎn)的數(shù)值優(yōu)化方法闡述為：1、給定自變量X的試探值，并設(shè)定系數(shù)矩陣u和v；2、由公式(22)和公式(24)計(jì)算偏導(dǎo)數(shù)矩陣FF′(X)；3、由公式(25)計(jì)算新的自變量取值Xk+1。4、判斷是否收斂。如果ΔXX＝XXk+1-XXk小于誤差限，計(jì)算結(jié)束；否則轉(zhuǎn)至步驟2。本發(fā)明的方法實(shí)現(xiàn)了對(duì)Δvmin的快速求解，同時(shí)也求解出最優(yōu)攔截位置與攔截時(shí)間。該算法計(jì)算精度高，速度快，且對(duì)初值不敏感。計(jì)算流程圖如圖5所示。以下整理公式(22)中，各項(xiàng)的解析表達(dá)式。因?yàn)?part;f3∂P=∂2J∂P∂fM+∂2J∂P2·∂f1∂P-∂J∂P·∂2f1∂P2(∂f1∂P)2·(TM2π∂f2∂fM-∂f1∂fM)-∂J/∂P∂f1/∂P·∂2f1∂P∂fM∂f3∂fM=∂2J∂fM2+∂2J∂P∂fM·∂f1∂P-∂J∂P·∂2f1∂P∂fM(∂f1∂P)2·(TM2π∂f2∂fM-∂f1∂fM)∂J/∂P∂f1/∂P·(TM2π∂2f2∂fM2-∂2f1∂fM2)---(26)]]>并且，顯然所以，只需獲得J、f1、f2對(duì)P、fM的一階偏導(dǎo)，二階純偏導(dǎo)，以及二階混合偏導(dǎo)，即可獲得公式(22)的解析表達(dá)式。將矢量投影在大地直角坐標(biāo)系，并約定h2xh2yh2z=cos(fM+w)cosΩ-sin(fM+w)cosisinΩcos(fM+w)sinΩ+sin(fM+w)cosicosΩsin(fM+w)sini---(28)]]>則首先，給出J對(duì)P、fM的一階偏導(dǎo)，二階純偏導(dǎo)，以及二階混合偏導(dǎo)。設(shè)JPy1=r2y(fM)g(P,fM)-r1yf(P,fM)g(P,fM)-v0y;JPy2=r2y(fM)-r1y[2-f(P,fM)]2·P·g(P,fM);]]>JPz1=r2z(fM)g(P,fM)-r1zf(P,fM)g(P,fM)-v0z;JPz2=r2z(fM)-r1z[2-f(P,fM)]2·P·g(P,fM);]]>r1h2＝r1xh2x+r1yh2y+r1zh2z；sr12h2=r12-(r1xh2x+r1yh2y+r1zh2z)2]]>則因?yàn)?part;r2∂fM=r2·eM·sin(fM)1+eM·cos(fM)]]>∂h2x∂fM∂h2y∂fM∂h2z∂fM=-sin(fM+w)cosΩ-cos(fM+w)cosisinΩ-sin(fM+w)sinΩ+cos(fM+w)cosicosΩcos(fM+w)sini]]>∂g∂fM=1μP·[∂r2∂fM·sr12h2-r2·rph2·r1h2sr12h2]]]>∂f∂fM=1P[∂r2∂fM(cosΔν-1)+r2·r1ph2r1]]]>∂∂fMf(P,fM)g(P,fM)=∂f∂fM·g-f·∂g∂fMg2]]>所以，設(shè)∂γMr2x(fM)g(P,fM)=(∂r2∂fM·h2x+r2·∂h2x∂fM)g-r2·h2x·∂g(P,fM)∂fMg2(P,fM);JfMx2=∂∂fMr2x(fM)g(P,fM)-r1x·∂∂fMf(P,fM)g(P,fM)]]>∂∂fMr2y(fM)g(P,fM)=(∂r2∂fM·h2y+r2·∂h2y∂fM)g-r2·h2y·∂g(P,fM)∂fMg2(P,fM);JfMy2=∂∂fMr2y(fM)g(P,fM)-r1y·f(P,fM)g(P,fM)]]>∂∂fMr2z(fM)g(P,fM)=(∂r2∂fM·h2z+r2·∂h2z∂fM)g-r2·h2z·∂g(P,fM)∂fMg2(P,fM);JfMz2=∂∂fMr2z(fM)g(P,fM)-r1z·∂∂fMf(P,fM)g(P,fM)]]>則∂J∂fM=2·(JPx1·JfMx2+JPy1·JfMy2+JPz1·JfMz2)---(31)]]>∂2J∂P2=-12P∂J∂P+2·(JPx22+JPy22+JPz22)+1Pg∂f∂P·(JPx1·r1x+JPy1·r1y+JPz1·r1z)---(32)]]>因?yàn)樗砸驗(yàn)?part;2h2x∂fM2∂2h2y∂fM2∂2h2z∂fM2T=-h2xh2yh2zT]]>∂2g∂fM2=1μP[∂2r2∂fM2sr12h2-2·∂r2∂fMr1ph2·r1h2sr12h2-r2·(r1ph22-r1h22)·r12+r1h24sr12h23]]]>∂2f∂fM2=1r1P[2·∂r2∂fM·r1ph2-r2·r1h2-∂2r2∂fM2·(r1-r1h2)]]]>∂2∂fM2f(P,fM)g(P,fM)=∂2f∂fM2g-f∂2g∂fM2-2∂f∂fM∂g∂fM+2fg(∂g∂fM)2g2]]>所以，設(shè)∂2∂fM2r2·h2xg(P,fM)=(∂2r2∂fM2·h2x+2∂r2∂fM∂h2x∂fM-r2h2x)·g-r2·h2x·∂2g∂fM-∂2∂fMr2·h2xg(P,fM)·2g·∂g∂fMg2(P,fM)]]>∂2∂fM2r2·h2yg(P,fM)=(∂2r2∂fM2·h2y+2∂r2∂fM∂h2y∂fM-r2h2y)·g-r2·h2y·∂2g∂fM2-∂∂fMr2·h2yg(P,fM)·2g·∂g∂fMg2(P,fM)]]>∂2∂fM2r2·h2zg(P,fM)=(∂2r2∂fM2·h2z+2∂r2∂fM∂h2z∂fM-r2h2z)·g-r2·h2z·∂2g∂fM2-∂∂fMr2·h2zg(P,fM)·2g·∂g∂fMg2(P,fM)]]>則∂2J∂fM2=2·JfMx22+JPx1·[∂2∂fM2r2·h2xg(P,fM)-r1x∂2∂fM2f(P,fM)g(P,fM)]+JfMy22+JPy1·[∂2∂fM2r2·h2yg(P,fM)-r1y∂2∂fM2f(P,fM)g(P,fM)]+JfMz22+JPz1·[∂2∂fM2r2·h2zg(P,fM)-r1z∂2∂fM2f(P,fM)g(P,fM)]---(34)]]>其次，給出f1對(duì)P、fM的一階偏導(dǎo)，二階純偏導(dǎo)，以及二階混合偏導(dǎo)。因?yàn)棣?＝(2m-l2)P2+2klP-k2；∂a∂P=mk[(l2-2m)P2-k2]χ02;∂2a∂P2=2mk(2m-l2)2P3+3(2m-l2)Pk2+2k3lχ03]]>并且，設(shè)χ1=1-r1a(1-f);χ2=1-χ12;χ3=r1r2tanΔν2·1-cosΔνP;]]>χ4=r1r2tanΔν2(1-cosΔνP-1r1-1r2);∂2∂P21aP=4lm1P3-6km1P4;]]>φ1=arccos[1-r1a(1-f(P,fM))]+χ4aP;]]>φ2=r1r2(1-cosΔν)χ2·∂∂P1aP-χ3P·1aP+χ4·aP2·∂∂P1aP;]]>∂φ2∂P=k·∂2∂P21aP·[2-kaP]-(∂∂P1aP)2·[aP-k][kaP]12·[2-kaP]32+aP2·∂2∂P21aP+χ3·(1aPa·(km1P-lm)-2P2+12aP·∂∂P1aP)+χ4·141aP(P∂a∂P+a)·∂∂P1aP]]>所以設(shè)χ51＝χ61·2P2+χ62·2·P-2·r2·χ6·χ63；χ6＝r1-r1h2；χ61=(r1h2-r2)·∂r2∂fM+r2·r1ph2;χ63=χ6·∂r2∂fM-r2·r1ph2;]]>χ64=mkχ02;χ65=χ63-r2·χ6a·∂a∂fM;]]>χ7=r1-r1h2r1+r1h2;∂χ7∂fM=-r1·r1ph2(r1+r1h2)·sr12h2;]]>χ8=(r1-r1h2)·r2P-r2-r1;∂χ8∂fM=χ63P-∂r2∂fM;∂χ8∂P=-(r1-r1h2)·r2P2;]]>則，∂f1∂fM=∂g∂fM+32aμ∂a∂fMarccos[1-r1a(P,fM)(1-f(P,fM))]+aμχ65Pχ2+∂a∂fM·χ7·χ8+a·∂χ7∂fMχ8+a·χ7·∂χ8∂fMμP---(36)]]>∂2f1∂P2=34gP2+32[121aμ(∂a∂P)2+aμ∂2a∂P2]·φ1+3aμ∂a∂P·φ2+a3μ·∂φ2∂P---(37)]]>因?yàn)樗?part;2f1∂P∂fM=-12P∂∂fMg(P,fM)-1P1χ2aμ·χ63·(1P+∂χ2∂P1χ2)+12a∂a∂P·{32aμ∂a∂fMarccos[1-r1a(P,fM)(1-f(P,fM))]+aμχ65Pχ2}+32aμ{∂2a∂P∂fMarccos[1-r1a(P,fM)(1-f(P,fM))]+∂a∂fMkχ2·∂∂P1aP}-aμr2χ6[(∂∂P1aP∂a∂fM+1aP∂2a∂P∂fM)-1χ2-1aPχ22∂a∂fM∂χ2∂P]---(38)]]>因?yàn)?chi;9=∂2r2∂fM2χ6-2∂r2∂fMr1ph2+r2r1h2;χ10=1Pχ65χ2;∂χ2∂fM=χ1χ2·χ65aP]]>∂χ5∂fM=[2(∂r2∂fM)2+2r2∂2r2∂fM2]·sr12h22-8r2∂r2∂fMr1h2r1ph2+2r22(r1h22-r1ph22)]]>∂χ51∂fM[2r1ph2∂r2∂fM+(r1h2-r2)·∂2r2∂fM2-r2·r1h2-(∂r2∂fM)2]·2P2-2χ632-2r2χ6χ9+{[2(∂r2∂fM)2+(r1+2r2)·∂2r2∂fM2]·χ6-2(r1+2r2)∂r2∂fMr1ph2+(r1r2+r22)r1h2}·2P]]>∂2a∂fM2=P·(∂χ5∂fM·1χ0-χ64·∂χ51∂fM-2·χ5χ51χ02+2χ64χ512χ0)]]>∂χ10∂fM=1P·{∂2r2∂fM2·χ6-2∂r2∂fMr1ph2+r2r1h2-[r2a∂2a∂fM2+1a∂a∂fM∂r2∂fM-r2a2(∂a∂fM)2]·χ6+r2a∂a∂fMr1ph2}·χ2-χ65·∂χ2∂fMχ22]]>∂2χ7∂fM2=r1r1h2·[r1+(r1xh2x+r1yh2y+r1zh2z)]sr12h2+r1·r1ph2·{r1ph2sr12h2-[r1+(r1xh2x+r1yh2y+r1zh2z)]r1h2r1ph2sr12h2}[r1+(r1xh2x+r1yh2y+r1zh2z)]2sr12h22]]>∂2χ8∂fM2=χ9P-∂2r2∂fM2]]>所以∂2f1∂fM2=∂2g∂fM2+12a∂a∂fM{32aμ∂a∂fMarccos[1-r1a(P,fM)(1-f(P,fM))]+aμχ65Pχ2}+aμ{32∂2a∂fM2arccos[1-r1a(P,fM)(1-f(P,fM))]-32∂a∂fMχ10a+∂χ10∂fM}+∂2a∂fM2χ7χ8+a∂2χ7∂fM2χ8+aχ7∂2χ8∂fM2+2∂a∂fM∂χ7∂fMχ8+2∂a∂fMχ7∂χ8∂fM+2a∂χ7∂fM∂χ8∂fMμP---(39)]]>最后，給出f2對(duì)fM的一階偏導(dǎo)，二階偏導(dǎo)。因?yàn)?part;2E∂fM2=2{1+eM1-eM[1+cos(fM)]+1-eM1+eM[1-cos(fM)]}2sin(fM)·(1+eM1-eM-1-eM1+eM)]]>所以∂2f2∂fM2=∂2E∂fM2[1-eMcos(E)]+eMsin(E)·(∂E∂fM)2---(41)]]>至此，由公式(30)至公式(41)，并結(jié)合公式(26)和公式(27)，可以獲得公式(22)的解析表達(dá)式，方法可以完整執(zhí)行。結(jié)合圖2中的注釋，對(duì)一種可能的交戰(zhàn)情形作如下描述：t1時(shí)刻，攔截器E的位置矢量與速度矢量分別是并且，如果攔截器與目標(biāo)均不進(jìn)行任何機(jī)動(dòng)，二者將在P1點(diǎn)碰撞。此時(shí)，目標(biāo)M進(jìn)行軌道機(jī)動(dòng)，新的軌道根數(shù)為[Ω,i,aM,eM,ω,tpM]＝[1.25442346395840,1.31152621076928,5319700.61120683,0.516254059012830,4.56809019516139,1959.60308392449]；軌道的運(yùn)行周期TM＝3861.3607137952。則攔截器需要將速度v0修正至合適的速度v1，才能在最小燃料消耗，即最小速度修正量的前提下，在新的目標(biāo)軌道上攔截目標(biāo)M。以上是問(wèn)題的初始條件，因?yàn)楸景l(fā)明在推導(dǎo)過(guò)程中，已經(jīng)將該工程問(wèn)題轉(zhuǎn)化為帶有等式約束的最優(yōu)規(guī)劃問(wèn)題，所以在與傳統(tǒng)的搜索方法進(jìn)行對(duì)比的同時(shí)，也將與MATLAB軟件自帶的求解函數(shù)fmincon進(jìn)行對(duì)比。狀態(tài)量X＝[tPfM]T的上下限分別設(shè)為Xmin＝[tminPminfMmin]T＝[1001e6-π]TXmax＝[tmaxPmaxfMmax]T＝[8005e6π]T本發(fā)明和fmincon迭代初值的選取規(guī)則為：t0＝(tmin+tmax)/2；P0沿用攔截器E變軌前的軌道半通徑；fM0通過(guò)將t0代入公式(11)求解獲得。所以，該算例的初值為X0＝[t0P0fM0]＝[4502.4473e6-2.8802]將初值X0代入迭代公式(25)中，完成求解。全部的仿真程序在CPU為I5-2410M，主頻2.3GHz的個(gè)人計(jì)算機(jī)中以及MATLAB2014a仿真環(huán)境下完成，更高效的編譯環(huán)境將提高計(jì)算效率。表1是不同方法的計(jì)算結(jié)果對(duì)比。表1不同方法的計(jì)算精度對(duì)比可以看出，雖然三種方法的計(jì)算結(jié)果以及目標(biāo)函數(shù)的最小值基本都是一致的，但是把計(jì)算結(jié)果代入KKT條件(20)，會(huì)發(fā)現(xiàn)搜索法的結(jié)果不能很好地滿足KKT條件，所以它的計(jì)算精度遠(yuǎn)低于本發(fā)明的方法。其次，我們比較三種算法的計(jì)算效率，基于半分法的加速搜索算法的計(jì)算效率雖然遠(yuǎn)高于fmincon，但是本發(fā)明的計(jì)算方法可以進(jìn)一步將耗時(shí)從13.47毫秒壓縮至0.490毫秒，計(jì)算效率大約提高27.5倍。表2是基于半分法的加速搜索迭代過(guò)程表2搜索法的迭代過(guò)程次數(shù)t(s)P(m)fM(rad)ΔV(m/s)14501.1307e6-2.88021698.384922757.8614e6-2.99123698.53613362.52.9945e6-2.9364780.79294406.251.8565e6-2.9085937.52835384.3752.3607e6-2.9225670.65566373.43752.6593e6-2.9294665.82087378.90632.5057e6-2.9260653.23588376.17192.5814e6-2.9277655.62999377.53912.5433e6-2.9268653.466810378.22272.5244e6-2.9264653.111811378.56452.5150e6-2.9262653.114212378.39362.5197e6-2.9263653.098113378.30812.5221e6-2.9263653.101214378.35082.5209e6-2.9263653.098715378.37222.5203e6-2.9263653.098216378.38292.5200e6-2.9263653.0981表3是本發(fā)明的迭代過(guò)程。表3本發(fā)明的迭代過(guò)程次數(shù)t(s)P(m)fM(rad)ΔV(m/s)14502.4473e6-2.8802333.77812377.39692.5651e6-2.9275655.31643378.38542.5196e6-2.9263653.14844378.38632.5199e6-2.9263653.09815378.38662.5199e6-2.9263653.0981綜上，本發(fā)明基于最優(yōu)規(guī)劃理論，通過(guò)對(duì)攔截變軌最小速度修正問(wèn)題進(jìn)行數(shù)學(xué)建模，經(jīng)過(guò)大量的公式推導(dǎo)，獲得了Δvmin對(duì)應(yīng)的KKT條件，并通過(guò)設(shè)計(jì)牛頓迭代算法，獲得了解析的梯度信息，實(shí)現(xiàn)了對(duì)Δvmin的快速求解，同時(shí)也求解出最優(yōu)攔截位置與攔截時(shí)間。經(jīng)過(guò)比對(duì)驗(yàn)證，該算法計(jì)算精度更高，速度更快。當(dāng)前第1頁(yè)1 2 3