張宇靖， 鐘睿

(北京航空航天大學(xué) 宇航學(xué)院, 北京 100083)

繩系衛(wèi)星系統(tǒng)是由系繩連接繩端衛(wèi)星構(gòu)成的空間系統(tǒng)[1]。繩系衛(wèi)星系統(tǒng)表現(xiàn)出了廣闊的應(yīng)用前景，如空間發(fā)電、構(gòu)建空間結(jié)構(gòu)和拖拽離軌等，是近年來(lái)航天研究熱點(diǎn)之一[2]。空間繩系拖拽離軌是借助飛網(wǎng)等機(jī)構(gòu)抓捕，并利用系繩連接主星和目標(biāo)，由主星機(jī)動(dòng)，實(shí)現(xiàn)目標(biāo)拖拽轉(zhuǎn)移的新概念在軌操作技術(shù)[3-5]。在空間碎片主動(dòng)移除技術(shù)中，空間繩系拖拽離軌是一種較為高效和具有應(yīng)用前景的技術(shù)。

在實(shí)際拖拽離軌任務(wù)中，系繩的擺動(dòng)會(huì)引起主星的擾動(dòng)，進(jìn)而影響整個(gè)系統(tǒng)的穩(wěn)定，因此如何抑制離軌過(guò)程中的系繩擺動(dòng)是拖拽離軌任務(wù)的一個(gè)關(guān)鍵問(wèn)題。目前已有一些相關(guān)的理論研究但還不廣泛。文獻(xiàn)[6]進(jìn)行了圓軌道上繩系系統(tǒng)捕捉目標(biāo)的動(dòng)力學(xué)研究。文獻(xiàn)[7]進(jìn)一步進(jìn)行了雙繩系系統(tǒng)進(jìn)行碎片主動(dòng)移除的動(dòng)力學(xué)研究。基于直接算法的最優(yōu)控制理論被廣泛應(yīng)用于系繩回收控制領(lǐng)域，直接算法以其通用性能夠很好地實(shí)現(xiàn)非線性動(dòng)力學(xué)過(guò)程控制而避免對(duì)動(dòng)力學(xué)本身的討論。在這一領(lǐng)域，Williams[8]、Wen等[9]的研究具有代表性。文獻(xiàn)[10]研究得到繩系系統(tǒng)星體姿態(tài)擺動(dòng)將激發(fā)系統(tǒng)高階擺動(dòng)并影響繩系系統(tǒng)的穩(wěn)定，文獻(xiàn)[11]對(duì)此提出保持系數(shù)張力下抑制廢星姿態(tài)擺動(dòng)的控制方法并且文獻(xiàn)[12]對(duì)此方法進(jìn)行了進(jìn)一步研究。

繩系系統(tǒng)的運(yùn)動(dòng)學(xué)具有強(qiáng)非線性，而且拖拽離軌過(guò)程中存在許多約束。故而需要在考慮這些約束因素下設(shè)計(jì)合適的控制律抑制拖拽離軌過(guò)程中系繩擺動(dòng)。從20世紀(jì)70年代開(kāi)始，針對(duì)具有復(fù)雜約束優(yōu)化需求的控制問(wèn)題，多年來(lái)大量工業(yè)實(shí)際中的成功應(yīng)用表現(xiàn)出預(yù)測(cè)控制在解決這類(lèi)問(wèn)題的優(yōu)勢(shì)[13-15]。本文針對(duì)霍曼變軌推力段的繩系系統(tǒng)系繩擺振抑制問(wèn)題，建立了軌道平面內(nèi)繩系系統(tǒng)動(dòng)力學(xué)模型；以最優(yōu)化結(jié)果為參考軌跡、非線性模型預(yù)測(cè)控制方法為基礎(chǔ)設(shè)計(jì)繩系系統(tǒng)拖拽離軌過(guò)程系繩平穩(wěn)控制律，并通過(guò)仿真驗(yàn)證了控制方法的跟蹤能力。

1 繩系系統(tǒng)動(dòng)力學(xué)建模

1.1 建模假設(shè)和坐標(biāo)系定義

為突出關(guān)鍵問(wèn)題，做如下假設(shè)，主要考慮系繩和主星的面內(nèi)姿態(tài)運(yùn)動(dòng)，即假設(shè)主星的滾轉(zhuǎn)角、偏航角，以及系繩的面外擺角為小角度，目標(biāo)視為質(zhì)點(diǎn)；不考慮系繩的彈性和質(zhì)量；忽略除重力外的其余空間干擾力，主動(dòng)力僅考慮推力，作用于主星質(zhì)心；忽略除重力梯度力矩外干擾力矩，控制力矩作用于主星。

首先建立軌道坐標(biāo)系So，zo軸沿地心連線方向由地球指向繩系衛(wèi)星系統(tǒng)(TSS)質(zhì)心，xo軸位于TSS質(zhì)心軌道平面內(nèi)與zo軸垂直并指向軌道運(yùn)動(dòng)方向，yo軸遵循右手準(zhǔn)則，三軸單位矢量依次為eox、eoy和eoz。另外定義TSS本體坐標(biāo)系Sb，zb軸沿著系繩指向主星，xb軸垂直系繩且在軌道平面內(nèi)，三軸單位矢量依次為ebx、eby和ebz。為描述主星姿態(tài)運(yùn)動(dòng)，定義主星本體坐標(biāo)系Sb1(三軸單位矢量依次為eb1x、eb1y和eb1z)，記Sb1相對(duì)軌道坐標(biāo)系So的滾動(dòng)、俯仰和偏航角分別為φ1、θ1、ψ1(3—1—2旋轉(zhuǎn)順序)。繩系系統(tǒng)坐標(biāo)系如圖1所示。

圖1 繩系系統(tǒng)坐標(biāo)系示意圖Fig.1 Schematic of coordinate system of tethered system

1.2 簡(jiǎn)化動(dòng)力學(xué)方程

利用上述假設(shè)，整個(gè)主星拖拽系統(tǒng)可以由包含系繩長(zhǎng)度、系繩面內(nèi)外擺角和主星姿態(tài)角，將動(dòng)力學(xué)方程線性化，消去高階項(xiàng)可以得到系統(tǒng)的動(dòng)力學(xué)方程[16]。為簡(jiǎn)化研究，忽略系統(tǒng)面外運(yùn)動(dòng)，本文研究霍曼變軌的第一次推力過(guò)程，主星初始位于同步軌道，變軌目標(biāo)為墳?zāi)管壍溃虼讼到y(tǒng)運(yùn)行的橢圓軌道偏心率很小，可以認(rèn)為是圓軌道。在偏置耦合力當(dāng)中引入上述假設(shè)，并忽略包括二階導(dǎo)數(shù)在內(nèi)的高階小量。

簡(jiǎn)化繩系系統(tǒng)的動(dòng)力學(xué)方程式為[17]

(1)

(2)

(d1zsin(θ1-φ)+d1xcos(θ1-φ))Fte+T1y

(3)

式中：φ為系繩面內(nèi)擺角(軌道坐標(biāo)系So繞其y軸旋轉(zhuǎn)φ得到TSS本體坐標(biāo)系Sb)；I1x、I1y和I1z為主星三軸慣量；T1y為主星y軸控制力矩大?。籉te為系繩張力；m1和m2分別為主星、子星的質(zhì)量；m為系統(tǒng)總質(zhì)量；l為系繩原長(zhǎng)；ωo為軌道角速度；d1x、d1z分別為主星上x(chóng)、z方向上的繩系點(diǎn)偏置大??；Fp為主星發(fā)動(dòng)機(jī)推力大小；αp為主星發(fā)動(dòng)機(jī)推力相對(duì)于軌道運(yùn)動(dòng)方向的夾角。動(dòng)力學(xué)方程式(1)～式(3)形式較為簡(jiǎn)單，利用了小角度假設(shè)，因此僅供穩(wěn)定控制器設(shè)計(jì)所用。

為處理方便，采用的無(wú)量綱時(shí)間和無(wú)量綱長(zhǎng)度分別為

(4)

式中：lt為單位繩長(zhǎng)，可取為變軌前的TSS繩長(zhǎng)。

做如下定義:

(5)

則可將式(5)轉(zhuǎn)化為狀態(tài)空間形式的系統(tǒng)動(dòng)力學(xué)方程，x′=f(x,u)展開(kāi)如下:

(6)

(7)

2 基于最優(yōu)化求參考軌跡

參考軌跡采用最優(yōu)化的方法求得。最優(yōu)化方法一般描述為：尋找狀態(tài)和控制，使得性能指標(biāo)達(dá)到最小，同時(shí)滿足包括系統(tǒng)狀態(tài)方程，一系列始末約束和過(guò)程約束在內(nèi)的所有約束條件。下面針對(duì)由主星俯仰運(yùn)動(dòng)和系繩面內(nèi)擺振組成的主星系統(tǒng)面內(nèi)運(yùn)動(dòng)，具體分析應(yīng)采用的性能指標(biāo)和應(yīng)滿足的約束條件，以建立最優(yōu)問(wèn)題。

2.1 性能函數(shù)選取

希望在最短的時(shí)間內(nèi)完成所需的脈沖，因此將問(wèn)題視為時(shí)間最短的最優(yōu)問(wèn)題，其性能函數(shù)為

(8)

2.2 控制約束

1) 始末約束

在開(kāi)始變軌推力之前，繩系系統(tǒng)初值為

(9)

式中：φ0為系繩面內(nèi)擺角初值，π/2-φmax≤φ0≤π/2+φmax，φmax為系繩面內(nèi)擺角所允許的最大幅值。

變軌推力結(jié)束之后，繩系系統(tǒng)回到水平狀態(tài)，且繩長(zhǎng)變化在可接受的范圍內(nèi)，終值狀態(tài)如下：

(10)

式中：Λmin為出于防碰撞安全考慮所定義的最短無(wú)量綱繩長(zhǎng)；Λmax為硬件所限最大的無(wú)量綱繩長(zhǎng)。

2) 過(guò)程約束

系統(tǒng)的狀態(tài)變量應(yīng)滿足系統(tǒng)動(dòng)力學(xué)方程。對(duì)于控制量實(shí)際上也必然是有一個(gè)限度的。對(duì)控制量的約束這里主要考慮主星發(fā)動(dòng)機(jī)推力方向φp的約束。

為了保證安全性，約束如下：

系繩面內(nèi)擺角π/2-φmax≤φ≤π/2+φmax；主星俯仰角π/2-θ1max≤θ1≤π/2+θ1max；主星發(fā)動(dòng)機(jī)推力方向αp,min≤αp≤αp,max；優(yōu)化過(guò)程中對(duì)其他工程中重要性更低的量不做約束。其中φmax為系繩面內(nèi)擺角所允許的最大幅值；?1max為主星俯仰運(yùn)動(dòng)所允許的最大幅值；αp,min為主星發(fā)動(dòng)機(jī)推力方向所能達(dá)到的最小值，αp,max為主星發(fā)動(dòng)機(jī)推力方向所能達(dá)到的最大值。

最后，整個(gè)推力段為避免系繩松弛，系繩的張力為非負(fù)。

(11)

3) 時(shí)間約束

顯然推力時(shí)間不可能是不受限的，最短時(shí)間應(yīng)為推力方向一直維持軌道速度方向時(shí)，達(dá)到所需速度增量的時(shí)間，其表達(dá)式為

(12)

式中：ΔVp1為主星霍曼變軌第一段所需速度增量，其表達(dá)式為

(13)

其中：r1和r2分別為變軌前后的軌道半徑。

最大時(shí)間可人為設(shè)定，無(wú)量綱化即有

(14)

4) 推力約束

假定推力時(shí)間相對(duì)軌道周期而言很小(這個(gè)假設(shè)在同步軌道是成立的)，有限推力段在整個(gè)變軌過(guò)程中可以等效為脈沖推力，則推力提供的沖量在軌道方向應(yīng)提供霍曼變軌所計(jì)算的速度增量，而在垂直軌道方向?yàn)?，即

(15)

(16)

至此，系統(tǒng)面內(nèi)運(yùn)動(dòng)的最優(yōu)控制問(wèn)題已經(jīng)建立：尋找滿足控制約束的狀態(tài)和控制集合，使得如式(8)所示的性能指標(biāo)最小。

2.3 最優(yōu)問(wèn)題求解

本研究采用Legendre-Gauss-Lobatto(LGL)偽譜法，直接LGL偽譜法是目前應(yīng)用最廣的偽譜方法之一[18]，具有較高的收斂速度并且能方便地處理高階導(dǎo)數(shù)，本文僅將此方法作為求解參考軌跡的方法，優(yōu)化方法本身不是本文的重點(diǎn)。通過(guò)偽譜法將問(wèn)題化為一個(gè)帶約束的非線性規(guī)劃問(wèn)題，有成熟的尋優(yōu)方法可以利用，本文采用MATLAB軟件自嵌的函數(shù)fmincon進(jìn)行求解。通過(guò)求解這個(gè)最優(yōu)控制問(wèn)題得到參考軌跡。

3 基于非線性模型預(yù)測(cè)控制律設(shè)計(jì)

3.1 模型離散化

為了設(shè)計(jì)控制器，首先要將連續(xù)時(shí)間系統(tǒng)模型轉(zhuǎn)化為離散時(shí)間系統(tǒng)。將第1節(jié)得到的非線性連續(xù)時(shí)間系統(tǒng)x′=f(x,u)在(x(t),u(t-1))處線性化得到如下線性時(shí)變系統(tǒng)[19]：

Bt(u(t)-ur(t))

(17)

式中：xr(t)為參考狀態(tài)；ur(t)為參考控制量，由第1節(jié)給出;At和Bt的計(jì)算公式如下：

(18)

設(shè)定系統(tǒng)的采樣周期為T(mén)，得到離散化系統(tǒng)如下[20]：

x(k+1)=fd(x(k),u(k))=

Ad,tx(k)+Bd,tu(k)+dd,t

(19)

式中：

(20)

dd,t=xr(k+1)-Ad,txr(k)-Bd,tur(k)

(21)

其中：I為單位向量。

3.2 性能函數(shù)選取

模型預(yù)測(cè)控制的目標(biāo)就是跟蹤期望的主星的三軸姿態(tài)和系繩擺角，即跟蹤參考軌跡。在離軌拖拽過(guò)程中，對(duì)于給定的主星參考姿態(tài)及系繩參考擺角，使得在空間拖拽過(guò)程中，主星姿態(tài)和系繩擺角的模型預(yù)測(cè)控制下的軌跡與最優(yōu)化得到的參考軌跡相同。因此以主星姿態(tài)和系繩擺角速度跟蹤誤差、控制力矩和終端誤差作為優(yōu)化目標(biāo)，建立性能函數(shù)如下[15]：

e(Np)TPe(Np)

(22)

式中：e(k)為跟蹤誤差；Q、R和P為相應(yīng)的跟蹤誤差、控制量和終端加權(quán)矩陣；Np和Nu分別為預(yù)測(cè)和控制時(shí)域。

3.3 主星姿態(tài)和系繩擺角的預(yù)測(cè)

利用繩系系統(tǒng)當(dāng)前時(shí)刻k的信息x(k)，包括主星姿態(tài)、系繩繩長(zhǎng)和繩長(zhǎng)變化，得到k時(shí)刻非線性離散狀態(tài)空間模型如下[20]：

(23)

設(shè)定預(yù)測(cè)時(shí)域?yàn)镹p，則可以根據(jù)繩系系統(tǒng)的當(dāng)前狀態(tài)信息，利用離散化的繩系系統(tǒng)動(dòng)力學(xué)及運(yùn)動(dòng)學(xué)方程，經(jīng)過(guò)迭代計(jì)算，得到預(yù)測(cè)時(shí)域Np內(nèi)的繩系系統(tǒng)狀態(tài)信息：

(24)

式中：x(k|k)=x(k)；y(k+i|k)為k時(shí)刻對(duì)未來(lái)k+i時(shí)刻系統(tǒng)輸出y(k+i)的預(yù)測(cè)值；u(k+i|k)為k時(shí)刻計(jì)算k+i時(shí)刻作用于系統(tǒng)的控制量。

因此得到預(yù)測(cè)時(shí)域Np內(nèi)的控制量{u(k),u(k+1),…,u(k+Np-1)}。同時(shí)假設(shè)在區(qū)間[Np，Nu]內(nèi)控制量保持不變，則可以得到

u(k+Nu-1)=u(k+Nu)=…=

u(k+Np-1)

(25)

因此利用繩系系統(tǒng)狀態(tài)模型和系統(tǒng)過(guò)去的控制量或輸出量，通過(guò)以上迭代計(jì)算出系統(tǒng)Np步預(yù)測(cè)輸出為

y(k)=(y(k+1|k),y(k+2|k),…,

y(k+Np|k))

(26)

3.4 控制約束

1) 始末約束

在開(kāi)始變軌推力之前，繩系系統(tǒng)處于水平狀態(tài)，變軌推力結(jié)束之后，繩系系統(tǒng)回到水平狀態(tài)，這里約束和求解最優(yōu)問(wèn)題中的相同。

2) 過(guò)程約束

系統(tǒng)的狀態(tài)變量應(yīng)滿足系統(tǒng)動(dòng)力學(xué)方程。

模型預(yù)測(cè)控制過(guò)程約束包含跟蹤誤差的約束，對(duì)于狀態(tài)過(guò)程跟蹤誤差約束如下：

系繩面內(nèi)擺角跟蹤誤差-Δφmax≤Δφ≤Δφmax；主星俯仰運(yùn)動(dòng)跟蹤誤差-Δθ1max≤Δθ1≤Δθ1max。

優(yōu)化過(guò)程中對(duì)其他工程中重要性更低的量不做約束。其中Δφmax為系繩面內(nèi)擺角跟蹤誤差所允許的最大幅值；Δθ1max為主星俯仰運(yùn)動(dòng)跟蹤誤差所允許的最大幅值。

系繩張力約束與最優(yōu)問(wèn)題中約束相同。

3) 推力約束

這里約束和求解最優(yōu)問(wèn)題中的相同。

3.5 非線性模型預(yù)測(cè)控制律實(shí)現(xiàn)

為了方便運(yùn)算，將跟蹤誤差加權(quán)矩陣、控制力矩加權(quán)矩陣和終端加權(quán)矩陣設(shè)計(jì)為單位矩陣。并由式(19)和以上步驟可以形成如下滿足目標(biāo)所要求的約束條件的繩系系繩擺振抑制控制的優(yōu)化問(wèn)題[15]。

(27)

同樣通過(guò)偽譜法將問(wèn)題化為一個(gè)帶約束的非線性規(guī)劃問(wèn)題，并采用MATLAB軟件自嵌的函數(shù)fmincon進(jìn)行求解。

4 基于非線性模型預(yù)測(cè)控制方法仿真驗(yàn)證

取仿真參數(shù)如下：主星初始俯仰角為90°(90°為水平)，俯仰角速度為0；系繩初始擺角85°(90°為水平)，初始擺角速度為0。且繩系系統(tǒng)參數(shù)如表1所示。

設(shè)定拖拽約束要求如下：最短繩長(zhǎng)10 m，最大繩長(zhǎng)1 000 m；主星的姿態(tài)角擺動(dòng)幅度不大于2°，姿態(tài)角速度不大于0.2 (°)/s；系繩擺動(dòng)幅度不大于5 (°)/s，擺動(dòng)角速度不大于0.2 (°)/s。主星推力方向最小值-π/4，主星推力方向最大值π/4。

設(shè)定模型預(yù)測(cè)控制誤差約束：系繩面內(nèi)擺角跟蹤誤差所允許的最大幅值0.1°；主星俯仰運(yùn)動(dòng)跟蹤誤差所允許的最大幅值0.1°。

通過(guò)最優(yōu)化方法求的最優(yōu)軌跡即參考軌跡如圖2所示。使用模型預(yù)測(cè)控制得到的跟蹤路徑和最優(yōu)化得到的參考軌跡的對(duì)比，如圖3所示。

從圖2中不難發(fā)現(xiàn)：采用最優(yōu)化方法得到的參考軌跡能使實(shí)際離軌拖拽繩系系統(tǒng)回到水平位置，并使系繩面內(nèi)擺角和主星俯仰角滿足約束要求，且系繩出于緊繃狀態(tài)。而且過(guò)程比較平滑。說(shuō)明最優(yōu)化得到的參考軌跡是可用的。

從圖3中可以看到：模型預(yù)測(cè)控制下的跟蹤路徑的各狀態(tài)很好地跟蹤了最優(yōu)路徑即最優(yōu)化得到的參考路徑。并且繩系系統(tǒng)的狀態(tài)控制在了約束范圍內(nèi)，即系繩面內(nèi)擺角和主星俯仰角滿足約束要求且系繩出于緊繃狀態(tài)。而且明顯可以通過(guò)縮小誤差約束范圍來(lái)提高跟蹤效果，但這樣又勢(shì)必提高計(jì)算量。又能從圖3看到，模型預(yù)測(cè)控制下的除了系繩面內(nèi)擺角速度外，其他跟蹤路徑比較平滑，但系繩面內(nèi)擺角速度抖動(dòng)比較厲害，不夠平滑，如工程需要，須在跟蹤控制中對(duì)擺角速度誤差適當(dāng)約束。模型預(yù)測(cè)算法具有在線優(yōu)化的特點(diǎn)，魯棒性較好，但是其計(jì)算量較大的問(wèn)題需要在今后進(jìn)一步研究。

表1 繩系系統(tǒng)參數(shù)

圖2 系繩面內(nèi)擺角、系繩面內(nèi)擺角速度、主星俯仰角、主星俯仰角速度和繩速的實(shí)際路徑和最優(yōu)路徑Fig.2 Actual path and optimal path of tether’s in-plane swing angle, tether’s in-plane swing angular velocity, main satellite pitching angle, main satellite pitching angular velocity and tether speed

圖3 系繩面內(nèi)擺角、系繩面內(nèi)擺角速度、主星俯仰角、主星俯仰角速度和繩速的跟蹤路徑和最優(yōu)路徑Fig.3 Tracking path and optimal path of tether’s in-plane swing angle, tether’s in-plane swing angular velocity, main satellite pitching angle, main satellite pitching angular velocity and tether speed

5 結(jié) 論

1) 考慮主星姿態(tài)、系繩擺角穩(wěn)定、系繩張力為正、系繩長(zhǎng)度變化有限等約束，針對(duì)繩系系統(tǒng)拖拽離軌推力作用過(guò)程中系繩擺振的抑制問(wèn)題，進(jìn)行系統(tǒng)狀態(tài)的最優(yōu)軌跡求解，最優(yōu)軌跡在滿足變軌速度增量要求的同時(shí)實(shí)現(xiàn)了推力時(shí)間最短。

2) 仿真結(jié)果表明，所設(shè)計(jì)的非線性模型預(yù)測(cè)控制方法很好地跟蹤了最優(yōu)軌跡，并且滿足控制約束條件；在變軌推力作用下主星姿態(tài)和系繩擺動(dòng)穩(wěn)定，且動(dòng)態(tài)過(guò)程平滑。