楊 博，魏 翔，于 賀，劉超凡

(北京航空航天大學(xué)宇航學(xué)院，北京 100191)

0 引 言

火星探測器在火星表面采集土壤和巖石樣品返回地球，可以用于探究火星的地質(zhì)特性[1]，為人類實現(xiàn)火星“移民”的目標(biāo)提供理論依據(jù)。由于火星與地球通信延遲大、火星大氣環(huán)境復(fù)雜等因素[2]，火星返回上升器的姿態(tài)控制系統(tǒng)需要考慮氣動力、液體晃動力等諸多干擾力，同時應(yīng)該具備自主容錯能力[3]。我國天問一號火星探測器的軟著陸GNC系統(tǒng)設(shè)計有干擾力和干擾力矩的快速識別算法，以及推力器故障快速診斷和重構(gòu)算法，使探測器具備了在著陸過程中的高容錯控制能力，提高了此過程姿態(tài)控制的魯棒性[4]。近年來，隨著火星返回任務(wù)研究大力開展，使用冗余執(zhí)行器的過驅(qū)動容錯控制設(shè)計思想，成為當(dāng)前容錯控制技術(shù)應(yīng)用的主要研究方向之一。

對于過驅(qū)動系統(tǒng)，控制分配是容錯技術(shù)應(yīng)用的典型方法[5-6]。在保持頂層控制律不變的條件下，通過對推力器進(jìn)行控制指令再分配，實現(xiàn)故障狀態(tài)下的容錯控制[7-8]。常見的控制分配方法有直接分配法、daisy chain法[9]，在此基礎(chǔ)上Zhou等[10]提出了兩種控制分配方法用于分析在推力器卡死故障下的容錯性能；Li等[11]利用基于優(yōu)化求解的控制分配方法提出了一種推力器卡死故障容錯控制方法；Wang等[12]提出基于自適應(yīng)滑模技術(shù)的故障條件下控制指令再分配方法；Galvan-Guerra等[13]提出利用積分滑?？刂萍夹g(shù)補(bǔ)償故障產(chǎn)生的干擾的容錯控制方法。上述學(xué)者所提的方法對于火星上升器推力器的故障均能產(chǎn)生一定的容錯控制效果，但都忽略了執(zhí)行器自身的動態(tài)特性對姿態(tài)閉環(huán)系統(tǒng)造成的影響。Yan等[14]研究表明執(zhí)行器動態(tài)特性能夠影響閉環(huán)系統(tǒng)的控制性能。董哲等[15]考慮執(zhí)行器工作特性時僅將響應(yīng)速度作為優(yōu)化指標(biāo)，忽略了執(zhí)行器動態(tài)特性引起的分配誤差。Bolling[16]通過對執(zhí)行器施加過驅(qū)動指令力圖抵消動態(tài)特性的影響，但是這種方法并沒有充分利用動態(tài)模型的信息，故而造成補(bǔ)償效果不理想。Hanger等[17]提出的基于模型預(yù)測的控制分配方法，在簡化意義上解決了考慮執(zhí)行器動態(tài)特性情況下的控制分配問題。

本文基于模型預(yù)測的基本控制分配思想，對火星上升器的推力器故障提出一種基于模型預(yù)測的推力再分配容錯控制新方法，其突出特點是依據(jù)推力器各自動態(tài)特性，綜合考慮分配誤差和燃料消耗設(shè)計優(yōu)化指標(biāo)，建立推力器故障模型作為優(yōu)化求解器的輸出約束，利用模型預(yù)測控制分配方法，實現(xiàn)故障狀態(tài)下的姿態(tài)控制系統(tǒng)快速鎮(zhèn)定。

1 火星上升器的控制系統(tǒng)推力分配模型

火星返回上升器姿態(tài)控制系統(tǒng)模型為仿射非線性模型：

(1)

式中：θ為上升器姿態(tài)控制模型中的狀態(tài)變量；虛擬控制量v為控制力及控制力矩沿上升器本體坐標(biāo)系的三軸分量。假設(shè)上升器帶有N個推力器，li=[li,x,li,y,li,z]T,i=1,2,…,N，代表推力器到上升器本體坐標(biāo)系的位置量，推力器噴管擺動使推力方向變化。設(shè)ob,xb,yb,zb為上升器本體坐標(biāo)系，oi,xi,yi,zi為推力器坐標(biāo)系，其關(guān)系如圖1所示。

圖1 推力器幾何位置示意圖Fig.1 Geometric position of the thruster

推力器i的推力Ti(t)映射到上升器本體坐標(biāo)系中的力Tb,i(t)和力矩Mb,i(t)表示為：

(2)

(3)

(4)

記N個推力器推力T為

(5)

根據(jù)式(3)，N個推力器產(chǎn)生的總的力和力矩為

(6)

(7)

映射矩陣G為

(8)

映射矩陣H為

(9)

根據(jù)效率矩陣K的定義,有

(10)

1.1 推力器線性二階環(huán)節(jié)

將火星上升器的推力器動態(tài)模型建模為線性二階環(huán)節(jié)，如下：

(11)

式中：xi,a表示第i個執(zhí)行器的動態(tài)模型的狀態(tài)變量；ui為控制輸入；yi表示第i個推力器的推力沿推力器坐標(biāo)系的三軸分量。假設(shè)各個執(zhí)行器的動態(tài)特性相互獨立，對N個推力器則有

(12)

1.2 推力器模型修正

若火星上升器推力器的動態(tài)模型存在建模參數(shù)誤差，則對式(12)進(jìn)行如下修正

(13)

假設(shè)誤差項ΔAa和ΔBa可以由一組特定的未知參數(shù)表示，如下所示

ΔAaxa+ΔBau=d(xa,u)?

(14)

假設(shè)：推力器負(fù)載動態(tài)模型的狀態(tài)量xa以及映射矩陣d(xa,u)具體表達(dá)式已知。

將式(14)代入式(13)中可以得到

(15)

用自適應(yīng)估計算法對式(15)模型參數(shù)進(jìn)行估計。定義如下的狀態(tài)估計方程：

(16)

(17)

(18)

(19)

根據(jù)如下的動態(tài)方程定義矩陣Qe∈p×p和Ce∈p:

(20)

(21)

假設(shè)存在時間點tc使得Qe(tc)是可逆的，那么對于?t≥tc有

(22)

式中：?c為不確定參數(shù)?的估計輸出。根據(jù)式(20)有

(23)

(24)

(25)

(26)

假設(shè)(Ak,Bk)是可控的，(Ak,C)是可觀的?？刂品峙淙蝿?wù)即求解輸入uk-1使得Kyk與虛擬控制量vk相等。

2 火星上升器推力控制分配

控制分配任務(wù)的總體框圖如圖2所示：

圖2 動態(tài)控制分配總體方案圖Fig.2 Scheme diagram of the dynamic control allocation

2.1 基于優(yōu)化求解的參考控制輸入規(guī)劃

由式(26)得到

(27)

(28)

式中：

(29)

寫為優(yōu)化問題，定義

(30)

(31)

式中：ω1,ω2,ω3為調(diào)節(jié)參數(shù)；ω1為控制量的模，表示燃料消耗最??；ω2表示分配誤差最小；ω3表示松弛因子最小。則約束條件為:

(32)

式中：Ωy表示輸出約束。

2.2 基于模型預(yù)測的控制分配

令性能指標(biāo)函數(shù)為：

(33)

式中：Q>0,R>0,Ψ>0,xk+i|k是在時刻k時的預(yù)測狀態(tài)量,xk|k=xk。

定義跟蹤誤差為：

(34)

將式(33)寫為:

(35)

由式(26)和式(34)，i=0,1,…,P-1，得到如下的約束方程:

(36)

系統(tǒng)參數(shù)Ak和Bk為自適應(yīng)估計方法得到的當(dāng)前時刻的估計值，則最優(yōu)控制為:

(37)

3 推力器容錯控制再分配設(shè)計

以推力器卡死故障為例，故障約束可以表示為推力沿三軸分量的線性約束，即存在矩陣M使得

Myk=0

(38)

說明：假設(shè)共有四臺推力器，推力器1卡死時，輸出推力T1存在如下所示的約束:

(39)

式中：T1x,T1y,T1z分別為T1沿推力器坐標(biāo)系的三軸分力;α*,β*,γ*分別為T1與推力器坐標(biāo)系三軸的夾角。記約束矩陣為

(40)

在預(yù)測時刻i, 1≤i≤P，推力器預(yù)測輸出yk+i|k的表達(dá)式為：

(41)

(42)

Skz+Tk=yk

(43)

式中：

(44)

(45)

若卡死故障發(fā)生，則將式(38)代入式(43)作為新增約束：

M(Skz+Tk)=0

(46)

模型預(yù)測跟蹤控制求解過程中，對于輸出yk+i|k的速率約束可以通過相鄰控制周期的時間間隔Ts轉(zhuǎn)換為位置約束，速率約束如下，0≤i≤P:

(47)

由式(47)和式(43)整理得到：

(48)

式中：

(49)

根據(jù)式(48)可以得到關(guān)于優(yōu)化變量z的線性約束:

(50)

4 仿真校驗

4.1 上升器姿態(tài)動力學(xué)模型

假設(shè)上升器在射面內(nèi)運動，帶液體晃動[19-20]，則火星上升器姿態(tài)動力學(xué)方程如下所示，沿xb軸方向:

T-FX

(51)

沿zb軸方向:

F-FZ

(52)

繞yb軸轉(zhuǎn)動:

mfgasin(θ+φ)=M+(F-FZ)b-MY

(53)

液體晃動等效單擺動力學(xué)方程:

(54)

火星大氣密度隨海拔高度變化可以建模為一維指數(shù)模型

ρh=ρ0e-(h-h0)/hs

(55)

假定姿態(tài)動力學(xué)方程(51)中推力T較大且為恒定常值，并且不考慮液體燃料的消耗，在上升器機(jī)動時，姿態(tài)角變化較小，晃動為微幅晃動，則液體晃動和氣動力對xb軸方向的加速度影響較小，因此可以用下面的近似方程代替

(56)

定義控制指令(u1,u2)如下所示

(57)

式中：

(58)

mgbsinθ-MY

(59)

(60)

則可以得到簡化的姿態(tài)動力學(xué)方程

(61)

式中：

(62)

(63)

式中：

(64)

利用反饋線性化設(shè)計狀態(tài)反饋控制器[21]：令S=[05×1,I5×5]，定義映射如下所示，

H0(x)=Sx

(65)

使用李導(dǎo)數(shù)記號，有下列表達(dá)式：

(66)

反饋輸入項為

(67)

v=u*+A(x)FS0x

(68)

其中反饋矩陣F為常值矩陣。

4.2 仿真條件

(1)仿真實驗中，推力器數(shù)目N=4，推力器動態(tài)模型為:

(69)

式中：自然頻率ω=40 rad/s；阻尼η=0.05。模型參數(shù)誤差不確定項映射矩陣d(xa,u)記為d:

(70)

不確定參數(shù)?維度p為1,真實值為2。

(2)上升器系統(tǒng)參數(shù)：m=600 kg,I=720 kg·m-2,mf=100 kg,a=0.32 m,b=0.25 m,T=3500 N,ε=0.19 kg·m2·s-1,If=90 kg·m-2；火星大氣參數(shù)：ρ0=2×10-4kg·m-3,h0=40000 m,hs=7500 m,h=5000 m；氣動參數(shù)：Cz=0.04,xd=0.75 m,xz=0.62 m,Sm=0.23 m2；控制器反饋矩陣取值為：

(71)

l=[l1,l2,l3,l4]=

(72)

四個推力器模型輸出約束為：

(73)

式中：推力最大值Tmax=2000 N，擺角最大值αmax=10°；推力變化速率最大為150 N/s。

4.3 仿真驗證及結(jié)果分析

4.3.1推力器為二階模型的控制分配

假定推力器為線性二階環(huán)節(jié)，見1.1節(jié)，得到如下仿真結(jié)果。

圖3所示為姿態(tài)角θ和液體晃動等效擺角φ的變化規(guī)律，相應(yīng)的角速度如圖4所示。

圖3 上升器姿態(tài)角和等效單擺擺角Fig.3 Attitude angle and equivalent single pendulum deflection angle of the ascent vehicle

圖4 上升器姿態(tài)角速度和等效單擺擺角速度Fig.4 Attitude angular velocity and equivalent single pendulum deflection angle velocity of the ascent vehicle

從圖3和圖4可以看出，上升器的姿態(tài)角和液體晃動的等效擺角在5 s后均趨于零，說明若推力器為二階系統(tǒng)模型時，實施推力控制分配，可以快速抑制由液體燃料晃動造成的姿態(tài)角誤差，提高了系統(tǒng)可控性。

四臺推力器的推力如圖5所示，擺角如圖6所示，描述了各推力器在參與推力動態(tài)控制分配過程中的變化規(guī)律。

圖5 四臺推力器推力Fig.5 Thrust of the four engines

圖6 四臺推力器擺角Fig.6 Swing angle of the four engines

各推力器推力最大值限制在2000 N以下，輸出推力趨于875 N，圖6中推力器擺角最大值為10°，滿足預(yù)先假定的推力限制條件，結(jié)果說明，四臺推力器不但可以將上升器控制在彈面內(nèi)，同時還可以沿xb軸方向使上升器上升獲得所需求的3500 N額定推力。

對比未使用1.2節(jié)模型預(yù)測跟蹤控制方法的推力變化曲線如圖7所示，明顯看出初始階段推力變化幅度大，而且實施時間長，倘若出現(xiàn)更大的隨機(jī)擾動，將很難保證推力分配使系統(tǒng)快速穩(wěn)定，甚至可能導(dǎo)致上升失敗。

圖7 未使用模型預(yù)測跟蹤控制時的推力Fig.7 Thrust without model predictive tracking control

使用偽逆法和加權(quán)偽逆法[22]結(jié)合線性狀態(tài)反饋跟蹤控制作為仿真對比，得到如下所示的實際控制力和控制力矩的結(jié)果:

根據(jù)圖8和圖9可以得出，本文提出的基于模型預(yù)測的控制分配方法能夠有效避免由于推力約束引起的控制量抖動，相對于偽逆法、加權(quán)偽逆法輸出的實際控制力和力矩更加平穩(wěn)，這是因為優(yōu)化求解的過程中不僅滿足推力大小和方向約束，同時考慮了預(yù)測時間域內(nèi)的分配誤差和燃料消耗指標(biāo)。

圖8 三種控制分配方法下的實際控制力Fig.8 Actual control force under the three control methods

圖9 三種控制分配方法下的實際控制力矩Fig.9 Actual control torque under the three control methods

4.3.2推力自適應(yīng)動態(tài)控制分配

根據(jù)1.2節(jié)做推力器模型修正，即對二階模型作參數(shù)估計，用2.1節(jié)的自適應(yīng)動態(tài)分配方法得到如下結(jié)果:

圖10 推力器動態(tài)模型參數(shù)估計Fig.10 Parameter estimation of the thrusters’ dynamic model

圖11 使用和未使用自適應(yīng)分配算法的推力偏差對比Fig.11 Comparison of thrust deviations with and without the adaptive allocation algorithm

由圖11可以得出，在1～2 s期間，未使用動態(tài)自適應(yīng)分配算法時推力發(fā)生明顯偏差，約130 N，這是由于此時推力指令變化幅度較大，如1.8 s時4號推力器的推力從1300 N下降到800 N以下。

而在使用動態(tài)自適應(yīng)分配方法后，推力與基準(zhǔn)值的偏差最大不到50 N，可以認(rèn)為，使用模型參數(shù)估計可以對推力指令實施有效跟蹤控制，可以有平滑輸出減少超調(diào)，在一定范圍使推力偏差減小60%以上。

4.3.3推力器容錯控制再分配

假定推力器發(fā)生卡死故障，對推力器卡死故障實施容錯控制再分配。

在1.95 s時，令1號推力器擺角卡死在-6°，此時式(38)中的約束矩陣為：

(74)

在2.45 s時，使用容錯控制再分配方法對推力進(jìn)行再分配，仿真結(jié)果如圖12所示。

圖12 上升器姿態(tài)角和等效單擺偏轉(zhuǎn)角Fig.12 Attitude angle and equivalent single pendulum deflection angle of the ascent vehicle

可以看出，收斂時間、穩(wěn)定誤差、超調(diào)等各項指標(biāo)均較優(yōu)越。

圖13和圖14所示為四臺推力器在容錯控制過程中各自的推力和偏角的變化過程，可以看出，本文提出的推力器容錯控制再分配方法對解決火星上升器推力器卡死故障是非常有效的。

圖13 四臺推力器推力Fig.13 Thrust of the four engines

圖14 四臺推力器擺角Fig.14 Swing angle of the four engines

5 結(jié) 論

本文針對火星返回上升器的推力器故障容錯技術(shù)研究，提出基于模型預(yù)測的容錯控制再分配方法，結(jié)論如下：根據(jù)火星返回器的推力器動態(tài)特性獲得基礎(chǔ)模型，并采用參數(shù)估計方法實時修正逼近真實模型，使推力器輸出誤差能夠減小60%以上；使用推力再分配技術(shù)，實現(xiàn)推力器容錯控制再分配，3～5 s即可實現(xiàn)對控制指令的穩(wěn)定跟蹤，姿態(tài)控制系統(tǒng)快速鎮(zhèn)定；用推力器故障狀態(tài)下的輸出約束調(diào)整優(yōu)化求解器的約束域，實現(xiàn)最小化分配誤差和最小燃料消耗意義下的最優(yōu)控制再分配。