航天器姿態(tài)控制的區(qū)間參數(shù)化方法*

汪 洋，龔 薇*，張江偉，余 琪，朱勃帆

0 引 言

航天器姿態(tài)控制是航天器執(zhí)行空間任務(wù)中重要組成成分，由于航天任務(wù)的多樣性，在空間任務(wù)執(zhí)行過(guò)程中，實(shí)現(xiàn)快速有效的控制顯得格外重要.通過(guò)給定驅(qū)動(dòng)力矩，使得航天器姿態(tài)達(dá)到期望姿態(tài)是航天器姿態(tài)控制的主要任務(wù).

在處理航天器最優(yōu)控制問(wèn)題上，一般被分為間接法和直接法，間接法是基于一階最優(yōu)性條件來(lái)解決問(wèn)題，直接法是解決最優(yōu)問(wèn)題的近似問(wèn)題[1].在解最優(yōu)控制問(wèn)題的近似問(wèn)題的方法中，有兩種常見(jiàn)的方法，一種是將外部控制表示為時(shí)間和一系列參數(shù)的顯示解析函數(shù)，通過(guò)逆動(dòng)力學(xué)粒子群優(yōu)化技術(shù)(IPSO)來(lái)解決最小時(shí)間機(jī)動(dòng)問(wèn)題[2-4]，另一種是將控制按照一定規(guī)律的近似值處理,如基于粒子群優(yōu)化[5,6]，偽譜法優(yōu)化[7]，混合二階段法[8].

解決最優(yōu)控制問(wèn)題的近似問(wèn)題時(shí)，最優(yōu)控制問(wèn)題的數(shù)值解都要從無(wú)限維的優(yōu)化問(wèn)題轉(zhuǎn)換到有限的參數(shù)優(yōu)化問(wèn)題[9].雖然文獻(xiàn)[5-8]都給出了滿足約束的航天器的最優(yōu)控制的數(shù)值解，但是最優(yōu)控制問(wèn)題的數(shù)值解切換時(shí)間點(diǎn)的問(wèn)題并沒(méi)有考慮到.時(shí)間變換的方法[10,11]可以解決這樣的問(wèn)題，但是過(guò)多的時(shí)間切換節(jié)點(diǎn)將導(dǎo)致計(jì)算量的迅速攀升.并且，在解決最優(yōu)控制問(wèn)題中，無(wú)效切換時(shí)間點(diǎn)是需要避免的.

本文考慮了航天器姿態(tài)的最優(yōu)控制問(wèn)題，提出一種新的區(qū)間參數(shù)化的計(jì)算方法，此方法如下所述：首先將航天器姿態(tài)控制的能耗，以及系統(tǒng)的收斂速度等性能指標(biāo)建立為目標(biāo)函數(shù)，考慮實(shí)際過(guò)程中的物理限制，例如陀螺儀測(cè)速限制的角速度約束，以及反作用輪的實(shí)際輸出限制等作為有界輸入和有界狀態(tài)約束.在解決這樣的最優(yōu)控制問(wèn)題的時(shí)候，一種控制參數(shù)化的方法被使用[12-16]，控制值將由分段常值函數(shù)表示，解決該問(wèn)題的近似問(wèn)題.在解決連續(xù)狀態(tài)不等式約束的時(shí)候，采用約束轉(zhuǎn)錄法[17]與局部平滑法相結(jié)合來(lái)處理連續(xù)狀態(tài)不等式約束． 經(jīng)過(guò)這樣的處理，原約束變?yōu)橐粋€(gè)標(biāo)準(zhǔn)約束的形式.同時(shí)問(wèn)題的可行域也會(huì)隨之變化，通過(guò)[9]提出的漸近調(diào)整相關(guān)參數(shù)的方法來(lái)解決.

在解決切換時(shí)間的問(wèn)題上，一種新的區(qū)間時(shí)間變換的方法被提出，通過(guò)改變各個(gè)子區(qū)間長(zhǎng)度，間接的使得切換時(shí)間點(diǎn)線性變換，從而解決最優(yōu)問(wèn)題切換時(shí)間點(diǎn)確定的難題.在經(jīng)過(guò)上述處理以后，此問(wèn)題可以作為數(shù)學(xué)規(guī)劃問(wèn)題，由任意一種優(yōu)化方法解決.一種基于梯度的miser3.2[18,21]軟件被用來(lái)解決這樣的問(wèn)題，并給出了相關(guān)的代價(jià)函數(shù)以及約束的梯度公式.數(shù)值仿真表明，提出的方法是有效的.

1 航天器姿態(tài)控制的模型以及問(wèn)題描述

在相對(duì)于慣性坐標(biāo)系的本體坐標(biāo)系下，航天器的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)模型可以由以下的公式表述[19]：

設(shè)x1(t)=q1(t)，x2(t)=q2(t)，x3(t)=q3(t)，x4(t)=q0(t)，x5(t)=ω1(t)，x6(t)=ω2(t)，x7(t)=ω3(t).

則衛(wèi)星姿態(tài)控制的動(dòng)力學(xué)方程為：

考慮如下的二次型性能指標(biāo)的代價(jià)函數(shù)：

(3)

并且在實(shí)際的空間任務(wù)中，一般需要限制控制轉(zhuǎn)矩和角速度完成衛(wèi)星姿態(tài)機(jī)動(dòng)過(guò)程[20].

考慮航天器反作用輪的輸出物理限制，限制反作用輪輸出力矩在0.4 N·m，引入如下控制約束：

另外考慮到陀螺儀的測(cè)速飽和，限制航天器的轉(zhuǎn)動(dòng)角速0.2 rad/s,引入如下的狀態(tài)約束：

因此航天器姿態(tài)控制的最優(yōu)化問(wèn)題可以表述成問(wèn)題P的形式：

(4)

2 控制參數(shù)化

問(wèn)題P是一個(gè)帶有連續(xù)狀態(tài)不等式約束以及終端約束非線性最優(yōu)控制問(wèn)題，傳統(tǒng)的最優(yōu)控制理論解決這樣的問(wèn)題是非常困難的.一種控制參數(shù)化的方法被提出用于解決這樣的帶有連續(xù)狀態(tài)不等式約束的最優(yōu)控制問(wèn)題.控制參數(shù)化將一個(gè)無(wú)限維的最優(yōu)控制問(wèn)題轉(zhuǎn)化為一個(gè)有限維最優(yōu)參數(shù)選擇問(wèn)題，后者更容易解決[12-16].

對(duì)于k=1,2,…,p一個(gè)分段常值函數(shù)將作為控制的p段組成成分，控制值切換時(shí)間序列為：

τ=[τ0τ1…τp]T

其中，τ0=0，τp=tf

控制將由如下函數(shù)的表達(dá)式近似：

其中，σk表示分段控制值，χI(t)是指示函數(shù)，其定義為：

其中，I表示σk對(duì)應(yīng)的時(shí)間區(qū)間.

經(jīng)過(guò)如上的變換原問(wèn)題P轉(zhuǎn)換成一個(gè)最優(yōu)參數(shù)選擇問(wèn)題Pc.

盡管控制參數(shù)化的方法是可行的，但是尋求的最優(yōu)控制是一個(gè)分段常值函數(shù).顯然，對(duì)這個(gè)分段常值函數(shù)來(lái)說(shuō)，切換時(shí)間點(diǎn)是未知的，如何精確的確定切換時(shí)間點(diǎn)是一件十分麻煩的事情.切換時(shí)間點(diǎn)變換的方法有效的解決了這個(gè)問(wèn)題，將切換時(shí)間點(diǎn)作為決策變量進(jìn)行求解，能夠有效的求出實(shí)際控制的最佳切換時(shí)間點(diǎn).

3 切換時(shí)間點(diǎn)變換

切換時(shí)間點(diǎn)變換法[10,11]是將切換時(shí)間點(diǎn)τi,1≤i≤p-1作為待優(yōu)化的控制參數(shù)進(jìn)行優(yōu)化.

為此，一個(gè)新的時(shí)間映射s∈[0,1],構(gòu)造一個(gè)從時(shí)域t∈[0,tf]的固定切換時(shí)間點(diǎn)k/p,k=1，2…,p到的s∈[0,1]可變時(shí)間切換點(diǎn)的變換.

定義如下的微分等式進(jìn)行變換：

初始條件：

t(0)=0

其中，θi≥0,i=1,2,…，p.

通過(guò)對(duì)式子以及初始條件積分，對(duì)于s∈[(k-1)/p,k/p),k=1,2,…,p有：

4 區(qū)間時(shí)間變換

值得注意的是，隨著控制段數(shù)的增加，切換時(shí)間點(diǎn)數(shù)量隨之增加，這將導(dǎo)致這個(gè)問(wèn)題的計(jì)算量迅速增加.因此，區(qū)間時(shí)間變換的方法被提出用于有效減少計(jì)算切換時(shí)間點(diǎn).

區(qū)間時(shí)間變換方法，將控制切換的時(shí)間序列分為不同的子區(qū)間，通過(guò)改變子區(qū)間的長(zhǎng)度，使得子區(qū)間內(nèi)的控制切換時(shí)間隨區(qū)間長(zhǎng)度線性變化，從而達(dá)到改變各個(gè)子區(qū)間控制值的切換時(shí)間的目的.在減小計(jì)算量的同時(shí)，優(yōu)化控制值切換時(shí)間.如圖1所示，傳統(tǒng)的時(shí)間切換點(diǎn)變化技術(shù)，是將每一個(gè)時(shí)間切換點(diǎn)進(jìn)行計(jì)算,如圖1紅線所示，而區(qū)間時(shí)間變換方法是將時(shí)間分為若干個(gè)子集(子集中的切換時(shí)間點(diǎn)等分)，優(yōu)化每個(gè)子集的所占的時(shí)間長(zhǎng)度，然后對(duì)單個(gè)集合等分延伸控制時(shí)間節(jié)點(diǎn)，如圖二黑線所示.

圖1 時(shí)間變換方法Fig.1 Time transformation method

因此，以控制切換時(shí)間序列分為三個(gè)子區(qū)間為例，將控制值進(jìn)行分區(qū)間處理，

定義如下的區(qū)間：

U1=[0,p1/p),U2=[p1/p,(p2-p1)/p),

U3=[(p2-p1)/p,p3/p],

V=[(k-1)/p,k/p),k=1,2,…,p

定義如下的微分等式進(jìn)行變換，

(5)

顯然可以得到，當(dāng)s∈V時(shí)，有：

利用上述的微分等式(3)，可以對(duì)上述控制系統(tǒng)f(x,u)進(jìn)行一個(gè)新的時(shí)間變換：

=vp(s)f(x(s),u(s))

可以證明，經(jīng)過(guò)如上的變換，Pc變換成問(wèn)題Pt.

可以證明如上的變換的滿足最優(yōu)解的收斂性定理.其證明過(guò)程類似文獻(xiàn)[11]中的定理一和定理二.

5 連續(xù)狀態(tài)不等式約束的計(jì)算方法

在本節(jié)中，應(yīng)用約束轉(zhuǎn)錄結(jié)合一種局部光滑的方法來(lái)處理不等式約束.

首先，將連續(xù)狀態(tài)不等式約束等效為如下的等式約束：

Gi(σp,θ1,θ2,θ3)=

i=1,2,…,6

由于等效的約束函數(shù)不可微，因此對(duì)于每個(gè)i=1,2,…,6，都用光滑函數(shù)(4)來(lái)近似.

式中給出的ε>0為可調(diào)參數(shù)，所以式子可微，但是可行域不同，因此需要應(yīng)用參數(shù)調(diào)整漸進(jìn)的方法來(lái)解決.具體而言，對(duì)于每一個(gè)i=1,2,…,6,設(shè)

(6)

Gi,ε(σp,θ1,θ2,θ3)=

于是將其替換后得到問(wèn)題Pε

利用數(shù)值法求解問(wèn)題Pε，得到的最優(yōu)解不滿足(2)的約束規(guī)范.于是引入如下等式：

Gi,ε,γ(σp,θ1,θ2,θ3)=-γ+Gi,ε(σp,θ1,θ2,θ3)≤0

替換得到問(wèn)題Pε,γ

可以證明上述的變換滿足最優(yōu)解收斂定理，其證明過(guò)程類似文獻(xiàn)[9]中的定理8.3.1和定理8.3.3，給出如下算法近似地解問(wèn)題P

對(duì)于任意的γ>0，ε>0，在此問(wèn)題中取ε=10-1，γ=5ε/16，εmin=10-3；

第三步：若第二步滿足轉(zhuǎn)至第五步，否則轉(zhuǎn)至第四步；

第四步：令γ=γ/2，轉(zhuǎn)至第一步；

第五步：令ε=ε/10，γ=γ/10，轉(zhuǎn)至第一步，當(dāng)ε≤εmin時(shí)，算法停止.

6 代價(jià)函數(shù)和約束的梯度公式

提出一種通過(guò)代價(jià)函數(shù)和連續(xù)狀態(tài)不等式約束對(duì)u的梯度求解問(wèn)題Pε,γ的計(jì)算方法，為了求解每一個(gè)u對(duì)應(yīng)的代價(jià)函數(shù)和連續(xù)狀態(tài)不等式約束的梯度，給出以下兩個(gè)定理.證明過(guò)程均類似文獻(xiàn)[9]中的定理5.2.1

定理1.代價(jià)函數(shù)的梯度公式為：

k=1,2,3

其中，H0(x,σp,θ1,θ2,θ3,λ0(s))為哈密爾頓函數(shù)

H0(x,σp,θ1,θ2,θ3,λ0(s))=

協(xié)態(tài)方程的解：

邊界條件為：

(λ0(1))T=0

定理2.連續(xù)狀態(tài)不等式約束函數(shù)的梯度公式如(5)，(6)所示.

(7)

k=1,2,3

(8)

其中，Hi(x,σp,θ1,θ2,θ3,λi(s))為哈密爾頓函數(shù).

Hi(x,σp,θ1,θ2,θ3,λi(s))=Gi(x,σp,θ1,θ2,θ3)+

協(xié)態(tài)方程的解：

邊界條件為：

(λi(1))T=0

7 計(jì)算結(jié)果與仿真

本文選用反作用輪為動(dòng)力機(jī)構(gòu)的航天器姿態(tài)控制問(wèn)題，考慮實(shí)際過(guò)程中的輸入有界約束，以及角速度有界約束，驗(yàn)證算法的有效性，?。?/p>

選取初始點(diǎn)θ=[100 100 100]，u=[0 0 0]，初始姿態(tài)Q0(t)=[0.512 0.765 0.14 0.3647]T，期望姿態(tài)，Qd(t)=[0 0 0 1]T.在(1)中的權(quán)值矩陣為Q=diag{100,100,100,100,10,10,10}，R=E進(jìn)行仿真測(cè)試.

最終得到最優(yōu)區(qū)間映射值，θ=[75.7786 151.905 83.093]，最優(yōu)值2.972.最優(yōu)控制輸入值如圖4所示.其狀態(tài)響應(yīng)曲線如圖2，圖3所示，角速度值和控制值均滿足實(shí)際的物理約束.仿真采取30段控制值近似，若通過(guò)全時(shí)間切換節(jié)點(diǎn)計(jì)算，將要計(jì)算30個(gè)切換時(shí)間節(jié)點(diǎn)變量，通過(guò)區(qū)間變換的方法，只需要計(jì)算3個(gè)區(qū)間長(zhǎng)度變量，計(jì)算的決策變量數(shù)，只需要占到10%.其控制響應(yīng)曲3如圖所示，控制值變化區(qū)間，第一個(gè)區(qū)間內(nèi)，變換時(shí)間間隔小，此時(shí)變化幅度較大，在第二個(gè)區(qū)間，控制轉(zhuǎn)矩的變換幅度較小，或基本不變，此區(qū)間較長(zhǎng)，在第三個(gè)區(qū)間，控制轉(zhuǎn)矩小幅頻繁變化，變換時(shí)間間隔較小，實(shí)現(xiàn)了區(qū)間時(shí)間變化，所得到的子區(qū)間長(zhǎng)度是最優(yōu)的，驗(yàn)證了方法的有效性.

圖2 最優(yōu)控制下的四元數(shù)響應(yīng)曲線圖Fig.2 Quaternion response graphs under optimal control

圖3 最優(yōu)控制下的角速度響應(yīng)曲線Fig.3 The angular velocity response curves under optimal control

圖4 最優(yōu)控制轉(zhuǎn)矩Fig.4 Optimal control torques

8 結(jié) 論

本文研究了一類航天器姿態(tài)控制問(wèn)題，通過(guò)區(qū)間化的參數(shù)化方法，最優(yōu)控制解的數(shù)值解由一系列不固定的切換時(shí)間點(diǎn)的分段常值函數(shù)近似，提出的區(qū)間時(shí)間變換方法，相比傳統(tǒng)的定區(qū)間長(zhǎng)度，以及時(shí)間變換方法，本方法在確定最優(yōu)切換時(shí)間基礎(chǔ)上，減小了切換時(shí)間點(diǎn)的計(jì)算量.并且考慮了實(shí)際執(zhí)行過(guò)程中的物理約束，運(yùn)用了約束轉(zhuǎn)錄和局部平滑的方法處理了連續(xù)狀態(tài)不等式約束.得到的最優(yōu)控制問(wèn)題的近似解接近最優(yōu)解并且最優(yōu)解滿足實(shí)際工作情況，因此在實(shí)際執(zhí)行過(guò)程中，可以當(dāng)作控制的初始估計(jì)，或者作為參考值，具有較大的使用價(jià)值.