孫云龍,袁長(zhǎng)清,李政廣

0 引 言

近年來(lái)，利用庫(kù)侖力對(duì)衛(wèi)星進(jìn)行編隊(duì)成為一種新興的編隊(duì)技術(shù).KING等[1]首次提出衛(wèi)星庫(kù)侖力編隊(duì)的概念.通過(guò)衛(wèi)星上的離子噴射器噴射離子，僅消耗數(shù)瓦特的功率使衛(wèi)星表面帶電，利用衛(wèi)星間的靜電引力或斥力形成控制力來(lái)對(duì)編隊(duì)實(shí)現(xiàn)控制，以實(shí)現(xiàn)預(yù)期目標(biāo).與傳統(tǒng)推進(jìn)器的基于動(dòng)量守恒依靠噴射氣體或離子對(duì)衛(wèi)星實(shí)現(xiàn)控制相比，庫(kù)侖力控制不會(huì)產(chǎn)生腐蝕性的羽狀物，不會(huì)影響臨近衛(wèi)星上的敏感元器件，直接延長(zhǎng)了衛(wèi)星的在軌運(yùn)行壽命.另外對(duì)于近距離編隊(duì)，庫(kù)侖控制力的精度可以達(dá)到毫牛的數(shù)量級(jí)，恰滿足編隊(duì)要求.

文獻(xiàn)[2]研究了二星庫(kù)侖編隊(duì)在地月共線平動(dòng)點(diǎn)和三角平動(dòng)點(diǎn)處軌道徑向動(dòng)力學(xué)與軌道徑向線性動(dòng)態(tài)分析.文獻(xiàn)[3]采用Lyapunov反饋控制法探究了考慮太陽(yáng)攝動(dòng)時(shí)二星庫(kù)侖編隊(duì)圓軌道和地月共線平動(dòng)點(diǎn)處軌道徑向穩(wěn)定性問(wèn)題.文獻(xiàn)[4]研究了編隊(duì)反饋控制策略以實(shí)現(xiàn)虛擬結(jié)構(gòu)控制.此控制策略基于推進(jìn)器能力來(lái)實(shí)現(xiàn)對(duì)衛(wèi)星三維運(yùn)動(dòng)的控制.文獻(xiàn)[5]分析了旋轉(zhuǎn)二星編隊(duì)的穩(wěn)定性，結(jié)果表明：徳拜長(zhǎng)度大于分離距離時(shí)，非線性旋轉(zhuǎn)運(yùn)動(dòng)局部穩(wěn)定，否則運(yùn)動(dòng)不穩(wěn)定，在忽略徳拜效應(yīng)時(shí)，受擾動(dòng)的面外運(yùn)動(dòng)總是穩(wěn)定的.文獻(xiàn)[6]研究了針對(duì)平動(dòng)點(diǎn)處二體庫(kù)侖繩系虛擬結(jié)構(gòu)重構(gòu)的非線性動(dòng)力學(xué)和閉環(huán)控制，該控制策略對(duì)針對(duì)執(zhí)行定位任務(wù)和魯棒重構(gòu)任務(wù)的衛(wèi)星編隊(duì)也是有效的.文獻(xiàn)[7]設(shè)計(jì)了針對(duì)一維受限庫(kù)侖編隊(duì)的二階電荷反饋控制策略，并通過(guò)探究系統(tǒng)總能量分析了庫(kù)侖編隊(duì)相對(duì)對(duì)稱運(yùn)動(dòng)條件.文獻(xiàn)[8]通過(guò)研究得出共線三星旋轉(zhuǎn)的解族，設(shè)計(jì)出來(lái)基于線性模型的反饋控制以實(shí)現(xiàn)共線虛擬繩系系統(tǒng)的穩(wěn)定.若系統(tǒng)角動(dòng)量等于預(yù)估的角動(dòng)量，則系統(tǒng)區(qū)域漸進(jìn)穩(wěn)定，據(jù)此可計(jì)算出名義(所需)電荷.文獻(xiàn)[9]探究了共線三星編隊(duì)的不變構(gòu)型解并證明出：對(duì)于任何共線編隊(duì)，存在無(wú)限組開(kāi)環(huán)電荷解使編隊(duì)可保持平衡.文獻(xiàn)[10]還探究了共線三星庫(kù)侖編隊(duì)的解族，分析了不變構(gòu)型解存在的區(qū)域，以確定軌跡可能的范圍.考慮到編隊(duì)構(gòu)型維持和機(jī)動(dòng)時(shí)固有的不穩(wěn)定性的挑戰(zhàn)，文獻(xiàn)[11]證明了對(duì)于徑向相連構(gòu)型，面內(nèi)擾動(dòng)漸進(jìn)穩(wěn)定.文獻(xiàn)[12]等針對(duì)面外運(yùn)動(dòng)進(jìn)行了推廣分析，并設(shè)計(jì)出簡(jiǎn)單控制律，成功消除了面外擾動(dòng)的變化范圍.在文獻(xiàn)[13]中，將傳統(tǒng)電推進(jìn)技術(shù)與庫(kù)侖力結(jié)合設(shè)計(jì)出了混合推進(jìn)策略以適用于近距離編隊(duì)飛行.文獻(xiàn)[14]給出了針對(duì)靜態(tài)四星編隊(duì)電荷解的分析工具.Lee等[15]設(shè)計(jì)了基于Lyapunov法的電荷反饋率來(lái)控制二星編隊(duì)的分離距離和分離率.研究表明徑向三星編隊(duì)也有這些性質(zhì)，電荷反饋控制率可維持編隊(duì)面內(nèi)運(yùn)動(dòng).

本文以地球同步軌道為衛(wèi)星編隊(duì)的運(yùn)行環(huán)境，建模時(shí)考慮德拜效應(yīng)，首先建立了共線N星庫(kù)侖力編隊(duì)沿徑向方向的相對(duì)運(yùn)動(dòng)動(dòng)力學(xué)模型.然后以四星編隊(duì)為驗(yàn)證算例，針對(duì)變量對(duì)非線性方程進(jìn)行線性化，根據(jù)線性化方程設(shè)計(jì)出LQR控制律并進(jìn)行數(shù)值仿真.考慮到建模誤差，添加積分項(xiàng)設(shè)計(jì)出改進(jìn)型的LQR控制律，并進(jìn)行了數(shù)值仿真驗(yàn)證.

1 動(dòng)力學(xué)建模

1.1 坐標(biāo)系的建立

如圖所示建立Hill(希爾)軌道坐標(biāo)系H:{Or,Oθ,Oh}，Or指向背離地心方向，Oθ指向衛(wèi)星編隊(duì)質(zhì)心的速度方向，Oh指向衛(wèi)星軌道角速度方向，Oh與Or,Oθ構(gòu)成右手系.則第i顆衛(wèi)星在Hill軌道坐標(biāo)系中相對(duì)于編隊(duì)質(zhì)心的位置矢量為ρi，也可表示為ρi=[xiyizi].

編隊(duì)中第i顆衛(wèi)星在地心慣性坐標(biāo)系中的位置矢量表示為ri，質(zhì)心在慣性坐標(biāo)中的位置矢量為rc， 第i顆衛(wèi)星在慣性坐標(biāo)系下的位置矢量可表示為：

(1)

1.2 動(dòng)能與勢(shì)能

根據(jù)速度合成定律，慣性系下第i顆衛(wèi)星的絕對(duì)速度表示為：

vi=vρ+vo

(2)

帶入式(2)得第i顆衛(wèi)星絕對(duì)速度為：

(3)

設(shè)第i顆衛(wèi)星的質(zhì)量為mi，則編隊(duì)系統(tǒng)動(dòng)能為：

(4)

本文中，徑向方向上編隊(duì)構(gòu)型欲保持穩(wěn)定，因此沿徑向方向編隊(duì)可視為剛體結(jié)構(gòu).如圖2、圖3所示建立本體坐標(biāo)系B:{b1,b2,b3}.其中b1方向沿ρn1方向，b1與b2b3構(gòu)成右手系.當(dāng)編隊(duì)與理想構(gòu)型無(wú)速度和位置偏差時(shí)，B系與Hill軌道坐標(biāo)系H完全重合.

編隊(duì)質(zhì)心條件為：

m1ρ1+m2ρ2+…+mnρn=0

(5)

由圖2可知：

(6)

聯(lián)立式(5)、(6)可得，在坐標(biāo)系B下，ρ1、ρ2、…、ρn可表示為：

(7)

歐拉角(θ,ψ,φ)分別為俯仰角、偏航角和滾轉(zhuǎn)角，分別表示坐標(biāo)系B與希爾坐標(biāo)系H之間的角度偏差.本文中不考慮衛(wèi)星繞徑向軸的旋轉(zhuǎn)，此時(shí)滾轉(zhuǎn)角φ的值為零.采用3-2-1(ψ-θ-φ)的歐拉角旋轉(zhuǎn)(如圖3所示).此時(shí)，坐標(biāo)系B與Hill坐標(biāo)系H之間的坐標(biāo)轉(zhuǎn)換矩陣為：

(8)

Hill坐標(biāo)系H中，衛(wèi)星的位置向量Hρi可表示成：

(9)

將式(7)、(8)帶入(9)可以得到，在希爾系H中各衛(wèi)星的位置矢量為：

(10)

從上式可以看出，各位置矢量是關(guān)于各衛(wèi)星間間距l(xiāng)i和俯仰角θ、偏航角ψ的函數(shù).對(duì)各個(gè)位置矢量求導(dǎo)可以得到各衛(wèi)星的速度矢量分別為：

(11)

其中：

將式(10)、(11)代入到式(4)中，可以得到徑向共線N星庫(kù)侖編隊(duì)的系統(tǒng)動(dòng)能T為：

(12)

其中：

對(duì)于N星編隊(duì)系統(tǒng)，其重力勢(shì)能可以表示為：

(13)

(14)

其中,u為rc的單位矢量，ti為ρi的單位矢量.

將式(10)代入式(14)中，得到重力勢(shì)能Vg表達(dá)式為：

(15)

衛(wèi)星在太空中運(yùn)行時(shí)，在受到太陽(yáng)光壓、地球非球形攝動(dòng)等因素的影響時(shí)，還受到太空中離子濃度的影響.處于太空的帶電體受其周?chē)鷰喾措姾傻碾x子的屏蔽的影響，其電勢(shì)隨距離長(zhǎng)度的增加呈現(xiàn)出指數(shù)衰減，此長(zhǎng)度用徳拜長(zhǎng)度λd表示，這種現(xiàn)象稱作徳拜效應(yīng).徳拜效應(yīng)對(duì)電勢(shì)的影響可表示為：

式中q為帶電量，L為距離帶電體的距離，λd為徳拜長(zhǎng)度.

太空中等離子密度在不同軌道高度處不同，徳拜長(zhǎng)度也因而不同.近地軌道處，徳拜長(zhǎng)度范圍為2～40 cm，僅適用于衛(wèi)星機(jī)構(gòu)的對(duì)接操作.同步軌道處徳拜長(zhǎng)度范圍為140～1 400 m，恰適用于星間間距在幾十到數(shù)百米的衛(wèi)星編隊(duì).深空中其范圍為8～24 m，也可運(yùn)用于更近距離的編隊(duì)飛行.

對(duì)于N星庫(kù)侖力編隊(duì)系統(tǒng)，每顆衛(wèi)星與其他不相鄰的衛(wèi)星間都具有庫(kù)侖力相互作用.故其庫(kù)侖力勢(shì)能可表示為：

(16)

其中kc為靜電引力常數(shù)，qi和qj分別為第i顆衛(wèi)星與第j顆衛(wèi)星的電荷量，λd為德拜長(zhǎng)度.

1.3 N星編隊(duì)相對(duì)運(yùn)動(dòng)方程

設(shè)系統(tǒng)拉格朗日函數(shù)為[5]：

L=T-(Vg+Ve)

(17)

拉格朗日方程為：

(18)

q=(θ,ψ,l1,l2…ln-1)

其中Qi為第i顆星所受到不包括地球重力和庫(kù)侖力的廣義力.

將式(12)、(15)、(16)代入式(18)得共線N星庫(kù)侖編隊(duì)相對(duì)運(yùn)動(dòng)方程為：

Ω2(3cos2θcos2ψ-1)]-

Ω2(3cos2θcos2ψ-1)]-

……

(19)

此處Qij為衛(wèi)星間電荷乘積 ，即Qij=qiqj.上述的N+1式具有很強(qiáng)的非線性，且衛(wèi)星之間的距離li(i=1,2…n-1)與俯仰角θ和偏航角ψ高度耦合，這對(duì)于控制器的設(shè)計(jì)及仿真校驗(yàn)都是不小的考驗(yàn).

本文選取地球同步軌道為庫(kù)倫編隊(duì)運(yùn)行環(huán)境.在這一軌道高度處，徳拜效應(yīng)對(duì)于編隊(duì)中各衛(wèi)星間庫(kù)侖力的作用比較有限.故下文在進(jìn)行建模時(shí)將略去徳拜效應(yīng)的影響.對(duì)于多星編隊(duì)的建模及仿真，仍不失一般性.

1.4 四星編隊(duì)相對(duì)運(yùn)動(dòng)方程

上文求得了共線N星庫(kù)侖編隊(duì)的相對(duì)運(yùn)動(dòng)方程組.由于其中各方程具有很強(qiáng)的非線性，各方程間相應(yīng)變量又高度耦合，因此完全對(duì)N星庫(kù)侖編隊(duì)進(jìn)行仿真校驗(yàn)是不實(shí)際的.不失一般性地，為了驗(yàn)證動(dòng)力學(xué)方程的正確，進(jìn)行有效的數(shù)值模擬，下文將進(jìn)行針對(duì)四星編隊(duì)的動(dòng)力學(xué)模型建立及數(shù)值仿真.

令a1=m1(m2+m3+m4)，a2=(m1+m2)(m3+m4)，a3=(m1+m2+m3)m4，a4=m1(m3+m4)，a5=m1m4，a6=(m1+m2)m4.

(20)

共線四星庫(kù)侖編隊(duì)在徑向方向相連，當(dāng)該編隊(duì)不受到其他攝動(dòng)因素的影響而達(dá)到理想構(gòu)型時(shí)(編隊(duì)相對(duì)于Hill軌道坐標(biāo)系保持靜態(tài)穩(wěn)定)，此時(shí)相對(duì)運(yùn)動(dòng)方程滿足以下平衡條件：

l1=l1f,l2=l2f,l3=l3f

(21)

式中l(wèi)1f為星1與星2之間的期望距離，l2f為星2與衛(wèi)星3之間的期望距離，l3f為星3與星4之間的期望距離.將上式分別帶入到運(yùn)動(dòng)方程，可得到保持共線四星庫(kù)侖靜態(tài)編隊(duì)平衡所需的參考電荷積為：

(22)

Q12f、Q13f、Q14f、Q23f、Q24f、Q34f為期望距離所對(duì)應(yīng)的參考電荷積.對(duì)于所期望的共線四星庫(kù)侖編隊(duì)來(lái)說(shuō)，為使庫(kù)侖力作為控制力，以保持衛(wèi)星編隊(duì)的靜態(tài)構(gòu)型達(dá)到期望，各衛(wèi)星的帶電量需要滿足以上各式.

對(duì)于俯仰和偏航運(yùn)動(dòng)方程分別以小俯仰角和小偏航角進(jìn)行線性化處理；式中的衛(wèi)星間距離l1、l2、l3以間距誤差δl1、δl2、δl3進(jìn)行線性化處理，電荷積Q12、Q13、Q14、Q23、Q24、Q34以電荷積誤差δQ12、δQ13、δQ14、δQ23、δQ24、δQ34進(jìn)行線性化處理，如下所示：

l1=l1f+δl1,l2=l2f+δl2,l3=l3f+δl3,

Q12=Q12f+δQ12,Q13=Q13f+δQ13,

Q14=Q14f+δQ14,Q23=Q23f+δQ23,

Q24=Q24f+δQ24,Q34=Q34f+δQ34

(23)

利用式(23)可得到線性化動(dòng)力學(xué)方程：

3Ω2(a1δl1+a4δl2+a5δl3)-

3Ω2(a2δl2+a4δl1+a6δl3)-

3Ω2(a3δl3+a5δl1+a6δl2)-

(24)

(25)

式(25)中前兩式為線性姿態(tài)動(dòng)力學(xué)方程，δL與ψ相互耦合.后3個(gè)式子為編隊(duì)關(guān)于徑向理想構(gòu)型間距誤差的線性微分方程.可以通過(guò)控制電荷使得δL、ψ方程漸近穩(wěn)定.

2 控制器設(shè)計(jì)

本節(jié)將基于線性化相對(duì)運(yùn)動(dòng)動(dòng)力學(xué)模型，僅利用庫(kù)侖力作為控制力，設(shè)計(jì)LQR控制律使徑向共線四星庫(kù)侖編隊(duì)保持穩(wěn)定構(gòu)型.考慮到模型誤差，設(shè)計(jì)了改進(jìn)型的LQR控制律，以提高傳統(tǒng)LQR控制器在控制編隊(duì)衛(wèi)星構(gòu)型時(shí)的魯棒性.

2.1 傳統(tǒng)LQR控制器

衛(wèi)星編隊(duì)的俯仰角θ不與偏航角ψ和間距誤差δl1、δl2、δl3耦合，不受庫(kù)侖力控制，所以不考慮俯仰運(yùn)動(dòng).下面僅利用衛(wèi)星間相互作用的庫(kù)侖力對(duì)編隊(duì)的面內(nèi)運(yùn)動(dòng)(即間距誤差δl1、δl2、δl3和偏航角ψ)進(jìn)行控制[17].首先，定義狀態(tài)向量x為：

控制器期望實(shí)現(xiàn)的目標(biāo)為x→0.

將線性化的共線四星庫(kù)侖編隊(duì)構(gòu)型相對(duì)運(yùn)動(dòng)動(dòng)力學(xué)模型表示為狀態(tài)空間形式為：

(26)

其中矩陣A1為8×8的狀態(tài)矩陣，矩陣B1為8×6控制矩陣，矩陣C1為4×8矩陣.

定義輸入量：u=[δQ12δQ23δQ34δQ13δQ24δQ14]T,通過(guò)設(shè)計(jì)輸入量u,使系統(tǒng)滿足性能指標(biāo)：

(27)

其中Q1和R1分別為狀態(tài)變量和輸入變量的加權(quán)矩陣.

當(dāng)控制輸入為：

(28)

性能指標(biāo)泛函數(shù)最小.式中矩陣P滿足Riccati微分方程，即：PA+ATP-PBR-1BTP+Q=0.

2.2 改進(jìn)LQR控制器

針對(duì)地球同步軌道處四星庫(kù)侖力衛(wèi)星編隊(duì)運(yùn)行時(shí)力學(xué)環(huán)境復(fù)雜、非線性強(qiáng)、控制變量耦合性 強(qiáng)、傳統(tǒng)LQR控制器魯棒性不足等現(xiàn)象，在傳統(tǒng)LQR控制器的基礎(chǔ)上，考慮到未建模攝動(dòng)力的影響及線性化產(chǎn)生的誤差，設(shè)計(jì)出一種改進(jìn)型的LQR控制器，以提高控制的魯棒性.

考慮如下式所示的帶有不確定項(xiàng)的動(dòng)力學(xué)系統(tǒng)的狀態(tài)空間表達(dá)式：

(29)

式中ω為未建模攝動(dòng)力.令ω=Dv，v為解決不確定項(xiàng)的輔助擴(kuò)展控制向量，取D=I-BB+.

為使系統(tǒng)在不確定擾動(dòng)ω(x)的作用下漸進(jìn)穩(wěn)定，可通過(guò)尋找反饋控制律[uv]T,使系統(tǒng)滿足下列目標(biāo)泛函數(shù).

uTHu+ρ2‖v‖2)dt

(30)

式中F和H均為對(duì)角矩陣；α≥0，β≥0，ρ≥0;且2ρ2‖v‖2≤β2‖x‖.

當(dāng)α=ρ=0，β=1，且F=Q1，H=R1時(shí)，即轉(zhuǎn)換為式(27)所示的傳統(tǒng)LQR控制.

為了使性能指標(biāo)泛函數(shù)J1最小，取控制輸入為：

(31)

式中矩陣P滿足Riccati微分方程，且

3 數(shù)值仿真

3.1 初始條件

本節(jié)針對(duì)徑向共線四星庫(kù)侖衛(wèi)星編隊(duì)飛行隊(duì)形保持問(wèn)題作為仿真算例.在Matlab/Simulink環(huán)境下進(jìn)行數(shù)值仿真，對(duì)共線四星庫(kù)侖力相對(duì)運(yùn)動(dòng)方程和控制算法的有效性進(jìn)行檢驗(yàn).航天器參數(shù)取值[16]為：m1=m2=m3=m4=50 kg,l1f=l2f=l3f=20 m,仿真參數(shù)的取值為：kc=8.99×109N·m2/C2,Ω=2.66×10-6rad/s,δψ(0)=0.1 rad，δl1(0)=δl2(0)=δl3(0)=0.6 m，q1=1.869 7 mC,q2=42.445 mC，q3=-46.569 mC，q4=-195.103 mC.

傳統(tǒng)LQR控制中，狀態(tài)變量的加權(quán)矩陣Q1為：

Q1=diag{0.001,0.1,0.001,0.01,0.01,0.001,0.1,0.001}

輸入變量的加權(quán)矩陣R1為：

R1=diag{108,108,108,107,107,106}

在改進(jìn)LQR的魯棒控制中，取α=ρ=β=1.狀態(tài)變量的加權(quán)矩陣Q2為：

Q2=diag{1,10,30,0.5,0.01,0.01,0.1,0.01,0.01,0.1,0.001,0.1}

輸入變量的加權(quán)矩陣R2為：

R2=diag{107,108,107,107,107,107}

3.2 仿真結(jié)果及分析

在MATLAB中進(jìn)行數(shù)值模擬.仿真結(jié)果如圖4所示.

圖4為L(zhǎng)QR控制策略下有無(wú)積分項(xiàng)姿態(tài)角ψ隨軌道時(shí)間的變化情況.可以看出，與無(wú)積分項(xiàng)相比，偏航角誤差的仿真趨勢(shì)與其基本相同，收斂略早.圖5、圖6及圖7分別表示LQR控制策略下有無(wú)積分項(xiàng)衛(wèi)星間距l(xiāng)1、l2和l3的變化情況.可以看出，加入積分項(xiàng)的控制有效的減少了間距收斂至穩(wěn)定的時(shí)間，減小了間距變化的振蕩幅度.在相同軌道時(shí)間時(shí)，帶積分項(xiàng)控制策略下間隔距離誤差優(yōu)于無(wú)積分項(xiàng).由于積分項(xiàng)的補(bǔ)償作用，在控制穩(wěn)定階段的穩(wěn)定誤差大大減少.

4 結(jié) 論

本文研究了徑向共線多體庫(kù)侖衛(wèi)星編隊(duì)在地球同步軌道處軌道動(dòng)力學(xué)與控制問(wèn)題.通過(guò)拉格朗日函數(shù)推導(dǎo)出了共線N星庫(kù)侖編隊(duì)在同步軌道處的一般動(dòng)力學(xué)方程.以四星庫(kù)侖編隊(duì)作為驗(yàn)證算例，針對(duì)共線四星衛(wèi)星編隊(duì)設(shè)計(jì)了LQR控制律和加入積分項(xiàng)的LQR控制律.通過(guò)Matlab/Simulink數(shù)值仿真分析(仿真過(guò)程中忽略德拜效應(yīng))，驗(yàn)證了改進(jìn)LQR控制律的有效性.仿真結(jié)果表明，改進(jìn)型LQR控制律能夠更快達(dá)到期望構(gòu)型，該控制可行、有效.基于本文研究，可進(jìn)一步研究多星庫(kù)侖編隊(duì)在同步軌道處的切向或法向軌道動(dòng)力學(xué).