固體微推力器陣列的推力估計(jì)與分配補(bǔ)償方法

楊 博，趙 旭，苗 峻，劉旭輝，龍 軍

(1.北京航空航天大學(xué)宇航學(xué)院，北京 100191；2.北京控制工程研究所，北京 100190)

0 引 言

微推力器陣列利用微機(jī)電技術(shù)在同一塊芯片上集成大量獨(dú)立可尋址的微推力器，具有高集成度、體積小、質(zhì)量輕、結(jié)構(gòu)簡單、無可移動部件、沒有液體泄漏、制作簡單等優(yōu)點(diǎn)。早在20世紀(jì)90年代，TRW公司就設(shè)計(jì)了這種推力器;此后，法國、韓國和日本的相關(guān)科研機(jī)構(gòu)也設(shè)計(jì)了自己的微推力器陣列;清華大學(xué)、國防科技大學(xué)等少數(shù)高校也制備了原理樣機(jī)進(jìn)行測試[1]。傳統(tǒng)推進(jìn)系統(tǒng)體積質(zhì)量比較大，而微納衛(wèi)星的體積質(zhì)量小、功率電壓低，所以傳統(tǒng)推進(jìn)系統(tǒng)并不適合應(yīng)用在微納衛(wèi)星上，而微推力器陣列克服了傳統(tǒng)推力器體積大質(zhì)量大的劣勢，因此受到了廣泛地關(guān)注和研究[2-3]。

采用微推力器陣列進(jìn)行控制的微納衛(wèi)星，是一個標(biāo)準(zhǔn)的過驅(qū)動系統(tǒng)。過驅(qū)動系統(tǒng)中存在的硬件冗余，使得針對系統(tǒng)的進(jìn)一步優(yōu)化成為可能，為了充分利用這一特點(diǎn)提升航天器性能，引入了推力分配環(huán)節(jié)來對上層控制律給出的需求推力進(jìn)行進(jìn)一步優(yōu)化，即引入推力分配將控制律計(jì)算得到的名義控制力分配給冗余配置的推力器[4]。針對推力分配問題，在許多文獻(xiàn)中已經(jīng)給出一些具體的方法，如查表法[5]、動態(tài)分配法[6]和負(fù)載均衡法[7-8]。然而在實(shí)際應(yīng)用中，僅控制律穩(wěn)定并不能保證引入了推力分配環(huán)節(jié)后系統(tǒng)依然穩(wěn)定，以上推力分配的文獻(xiàn)中均沒有針對這一問題進(jìn)行分析。雖然文獻(xiàn)[9]針對這一問題提出了一種類PD反饋的方法來保障引入推力分配后的閉環(huán)系統(tǒng)仍然穩(wěn)定，但該方法在推力分配環(huán)節(jié)采用了能夠獲得解析解的偽逆分配法，對于不采用偽逆法進(jìn)行推力分配的系統(tǒng)，其穩(wěn)定性仍然少有文獻(xiàn)進(jìn)行分析。

在實(shí)際應(yīng)用中，微推力器的推力測量也是一個難點(diǎn)，其推重比非常小，產(chǎn)生的推力也極小，對測量環(huán)境及測量方法也提出了很高的要求[3]。我國的微小推力器還處于發(fā)展階段，實(shí)現(xiàn)高精度的測量所需技術(shù)要求和成本都很高。因此考慮到受燃料消耗、推力器部件老化以及地面測試不準(zhǔn)確等原因的影響，實(shí)際使用時(shí)推力器產(chǎn)生的推力會偏離標(biāo)稱推力，為了保證控制系統(tǒng)的性能，在推力分配中也應(yīng)該考慮推力偏差的影響。在處理推力偏差問題上，Castaldi等[10]利用微分幾何方法來處理非線性系統(tǒng)中的執(zhí)行器故障估計(jì)問題，在飛行器容錯控制中獲得了較好的效果。Alessro[11]首次在控制分配模塊進(jìn)行推力估計(jì)，但由于配置于傳統(tǒng)的推力分配系統(tǒng)，并沒能證明系統(tǒng)的穩(wěn)定性。Ye[12]利用偏差導(dǎo)數(shù)信息以及自適應(yīng)觀測器的方法重構(gòu)執(zhí)行器乘性效率因子，但該方法需要知道推力偏差的變化率，實(shí)際應(yīng)用價(jià)值受限。

針對以上問題，本文引入效率因子來描述推力偏差，并利用二次規(guī)劃進(jìn)行效率因子估計(jì)，討論了估計(jì)算法的收斂性。在此基礎(chǔ)上，針對微推力器一次性點(diǎn)火的特性，結(jié)合效率因子估計(jì)值，建立混合整數(shù)規(guī)劃模型進(jìn)行推力分配，使得系統(tǒng)能夠估計(jì)推力偏差并對其進(jìn)行實(shí)時(shí)補(bǔ)償，最后研究了引入推力分配后系統(tǒng)穩(wěn)定性問題。

1 推力分配基本問題描述

對于過驅(qū)動系統(tǒng)，引入推力分配后控制系統(tǒng)的基本結(jié)構(gòu)如圖1所示。

推力分配的目的，即考慮執(zhí)行機(jī)構(gòu)冗余與物理約束，將控制律給出的名義控制力以某種優(yōu)化指標(biāo)分配給具體的執(zhí)行機(jī)構(gòu)，從而實(shí)現(xiàn)系統(tǒng)性能的進(jìn)一步提升。

考慮線性離散系統(tǒng)

x(t+1)=Ax(t)+Bv(t)

(1)

式中：x(t)∈Rn為t時(shí)刻的狀態(tài)向量，A∈Rn×n為系統(tǒng)矩陣，B∈Rn×m為輸入矩陣，v(t)∈Rm為t時(shí)刻的虛擬控制輸入。

對過驅(qū)動系統(tǒng)進(jìn)行推力分配，有

v(t)=Cu(t)

(2)

式中：C∈Rm×k為分配矩陣，u(t)∈Rk為實(shí)際控制力，m

在分配中考慮推力偏差，引入效率因子，則在推力偏差下式(2)可寫為

v(t)=CΨ(t)u(t)

(3)

式中：Ψ(t)=diag(ψ1(t),…,ψk(t))，ψi(t)表示第i個控制力的實(shí)際效率，ψi(t)=1表示該軸實(shí)際控制力與標(biāo)稱控制力相同。

推力分配問題的關(guān)鍵在于如何求解不定方程式(3)以及如何確定效率因子矩陣Ψ(t)。

2 推力分配模型

由于在實(shí)際使用中一般更關(guān)心微推力器陣列上推力器的消耗個數(shù)，將式(3)改寫為

v(t)=C′G(t)z(t)

(4)

式中：G(t)=diag(g1(t),…,gk(t))，gi(t)表示第i個推力器的實(shí)際效率，gi(t)=1表示該推力器實(shí)際推力與標(biāo)稱推力相同。z(t)∈Nk,zi表示第i個推力器陣列上消耗微推力器的個數(shù)，C′∈Rm×k為分配矩陣。

由于推力分配過程總是在一個時(shí)刻進(jìn)行，為了簡便，推力分配建模過程中不再引入t。

推力分配優(yōu)化目標(biāo)選擇燃料最優(yōu)，即選取目標(biāo)函數(shù)為

(5)

將目標(biāo)函數(shù)與約束結(jié)合，則可以寫成如下優(yōu)化問題

s.t.C′Gz=v,z∈Nk

(6)

由于執(zhí)行器不可能完美無誤執(zhí)行上層控制力，引入分配誤差弱化模型約束，取誤差變量為

C′Gz-v=e=e+-e-

(7)

式中：e為分配誤差，e+，e-中元素均非負(fù)。

則可以得到存在分配誤差的優(yōu)化模型

minJ=∑ki=1zis．t．e+-e-=C′Gz-vz∈Nk，e+≤emax，e-≤emax{

(8)

式中：emax為容許最大推力分配誤差。

針對該模型的求解問題是一個混合整數(shù)規(guī)劃問題，為便于求解將其轉(zhuǎn)化為標(biāo)準(zhǔn)優(yōu)化模型的形式

minlTys．t．Dy=-vy≥0，z∈Nk，e+≤emax，e-≤emax{

(9)

這樣就得到了一個標(biāo)準(zhǔn)的混合整數(shù)規(guī)劃模型，利用現(xiàn)已成熟的數(shù)學(xué)優(yōu)化算法(如分支定界法，內(nèi)點(diǎn)法等)就可以進(jìn)行優(yōu)化求解，進(jìn)行在容許分配誤差下的最省燃料分配。

注2. 無論取何種式(9)的變形體，約束中需有e+≤emax,e-≤emax，即分配誤差有界，該條件是穩(wěn)定性分析中的重要條件。為了保證式(9)有解，同時(shí)燃料盡可能少消耗，在已知G的情況下一般選擇emax為微推力陣列所能產(chǎn)生的單位推力。

注3. 假設(shè)一個推力器陣列上的微推力器推力均相等，陣列平貼于立方星表面，由于每個微推力器產(chǎn)生推力效應(yīng)相同，在不考慮姿態(tài)的情況下，可建模為式(9)的形式，只關(guān)心每個陣列上點(diǎn)火幾個微推力器。

注4. 在考慮小衛(wèi)星姿態(tài)的情況下，有時(shí)候需要針對推力器陣列上的每一個推力器進(jìn)行分配建模，即將式(9)的混合整數(shù)規(guī)劃模型變形為0-1規(guī)劃模型(混合整數(shù)規(guī)劃的特例)，1代表該微推力器點(diǎn)火，即有

minlTy′s．t．D′y′=-v，y′≥0z′∈Nk，0≤zi≤1，i=1，2，…，ke+≤emax，e-≤emaxì?í????

(10)

由于微推力器陣列上的微推力器一次性使用的特性，每一次分配結(jié)束后需要將C″矩陣中對應(yīng)該微推力器的列置0或剔除，防止下一次分配中再使用該微推力器。

3 效率因子估計(jì)算法

3.1 基于二次規(guī)劃的估計(jì)算法

從推力分配模型中可以看出，在已知效率因子的情況下，通過推力分配能夠在不修改上層控制律的情況下很好地補(bǔ)償實(shí)際推力器偏差，然而實(shí)際中，效率因子真值很難獲得，只能通過對系統(tǒng)狀態(tài)的測量來對效率因子進(jìn)行一個估計(jì)，而估計(jì)的優(yōu)劣程度則代表了最后推力分配中補(bǔ)償?shù)膬?yōu)劣。

考慮實(shí)際工藝，在制作一片微推力器陣列時(shí)，其每一個微推力器具有相同的推力大小，即每一片微推力器陣列上的推力器具有相同的效率因子，那么在估計(jì)效率因子時(shí)，可用該片已經(jīng)使用的推力器來估計(jì)其他推力器的推力。

考慮到二次規(guī)劃問題在優(yōu)化領(lǐng)域中研究得比較深入，只需要建立其數(shù)學(xué)估計(jì)模型，有非常成熟的優(yōu)化算法可以對其進(jìn)行求解。這里采用一種二次規(guī)劃進(jìn)行效率因子估計(jì)的方法。該估計(jì)算法通過對過去p個采樣點(diǎn)系統(tǒng)狀態(tài)信息的處理，得到在這p個點(diǎn)中效率因子增量的一個優(yōu)化估計(jì)值，根據(jù)效率因子增量估計(jì)值可以計(jì)算得到效率因子的估計(jì)值。

將推力分配式(4)代入模型(1),可得

x(t+1)=Ax(t)+BuG(t)z(t)

(11)

式中：Bu=BC′。

為了能夠進(jìn)行效率矩陣估計(jì)，首先定義效率矩陣增量如下

Φ=G(t)-G(t-1)

(12)

式中：Φ=diag(φ1，…,φk)，φi=gi(t)-gi(t-1)。

第一步：求解以下二次規(guī)劃問題

min ∑pi=1＾sTi＾sis．t． x(t)-Ax(t-1)-Bu(＾G(t-1)+＾Φ)· z(t-1)=＾s1x(t-1)-Ax(t-2)-Bu(＾G(t-1)+＾Φ)· z(t-2)=＾s2?x(t-p+1)-Ax(t-p)-Bu(＾G(t-1)+＾Φ)· z(t-p)=＾spì?í????????????

(13)

(14)

(15)

(16)

(17)

3.2 算法性能分析

為了討論二次規(guī)劃算法對于慢時(shí)變效率因子的估計(jì)能力，設(shè)t′為推力發(fā)生偏差的時(shí)刻，同時(shí)假設(shè)在發(fā)生偏差后一段時(shí)間內(nèi)，效率因子保持不變。

定義估計(jì)誤差如下

(18)

(19)

證.實(shí)際上可以假設(shè)t′時(shí)刻后真實(shí)效率因子矩陣為G′，那么在t>t′時(shí)

x(t+1)=Ax(t)+BuG′z(t)

(20)

在最小擬合誤差均為0的情況下，可以認(rèn)為存在多個增量矩陣Φ滿足式(14)中的約束條件。即

(21)

式中：Q代表優(yōu)化問題式(14)可行解的集合，

式(13)有解，則可行解約束方程

(22)

從定理1可以看出，要求對所有推力器的效率因子都能準(zhǔn)確估計(jì)，需要在過去p個采樣點(diǎn)中，所有推力器都有獨(dú)立的推力輸出，即對于p個推力器陣列，至少需要p個采樣點(diǎn)信息才能保證算法收斂。

根據(jù)推力分配模型以及推力估計(jì)算法，可以得到推力估計(jì)與分配補(bǔ)償方法的流程圖(見圖2)。

4 系統(tǒng)穩(wěn)定性分析

一般研究推力分配的文獻(xiàn)中很少討論引入分配之后系統(tǒng)的穩(wěn)定性問題，但實(shí)際上，由于引入了推力分配，同時(shí)采用了二次規(guī)劃的方法對效率矩陣進(jìn)行估計(jì)，那么控制系統(tǒng)的穩(wěn)定性必然會受分配誤差的影響，而分配誤差又主要受到推力估計(jì)誤差的影響。本文將對該部分影響進(jìn)行一個理論分析。

假設(shè)v′為推力分配給出的實(shí)際控制力，代入式(1)，可得

x(t+1)=Ax(t)+Bv′(t)

(23)

考慮狀態(tài)反饋控制

v(t)=-Kx(t)

(24)

則有

x(t+1)=Ax(t)+Bv′(t)+B(-Kx(t))-

B(-Kx(t))=(A-BK)x(t)+

B(v′(t)-(-Kx(t)))=(A-

BK)x(t)+B(v′(t)-v(t))

(25)

實(shí)際上

v′(t)-v(t)=C′Gz-v(t)=e

(26)

所以

x(t+1)=(A-BK)x(t)+Be

(27)

從式(9)的約束條件中可以得到

|e|≤emax

(28)

式中：|*|代表對每一個元素取絕對值。

式(27)所描述的運(yùn)動穩(wěn)定，首先需要(A-BK)的特征值均在單位圓內(nèi)，同時(shí)系統(tǒng)存在一個有界外擾動，針對這樣的系統(tǒng)，文獻(xiàn)[14]證明是一個有界穩(wěn)定系統(tǒng)，其最終收斂于集合C。

|V-1B|emax+|V|ζ}

(29)

以上推導(dǎo)基于式(9)存在解，為了保證式(9)總是有解，取emax中的最小分量為μ，需保證

(30)

注6.由于該收斂界的計(jì)算基于最大推力分配誤差，然而實(shí)際中并不會每次分配都以分配誤差最大的情況出現(xiàn)，所以該收斂界的計(jì)算具有一定的保守性。

5 仿真校驗(yàn)

考慮編隊(duì)近距離和主星圓軌道的情況，從星在Hill坐標(biāo)系下的相對運(yùn)動方程可以簡化成Clohessy-Wiltshire方程模型[15]：

(31)

式中：ωt為主星軌道平均角速度，f為控制加速度。

則有

(32)

式中：

取采樣周期T=0.1 s,對模型離散化有

x(t+1)=Ax(t)+Bv(t)

(33)

與文獻(xiàn)[1]類似，取編隊(duì)構(gòu)形為等邊三角形，三顆從星分別位于三角形的三個頂點(diǎn)，一顆主星位于三角形的中心，構(gòu)形尺度為1.5 m，四顆衛(wèi)星均為立方星，尺寸為10 cm×10 cm×10 cm，質(zhì)量為1 kg。主星軌道選擇低軌太陽同步軌道，各星初始時(shí)標(biāo)準(zhǔn)軌道參數(shù)如表1所示。

由于三顆從星控制效應(yīng)類似，為了節(jié)省篇幅，本文主要以從星1的結(jié)果來說明。取從星1的初始偏差e1=[2, -2.5, 2, 0.5, -1.2, 1]T，即從星1與其標(biāo)準(zhǔn)位置/速度的偏差?？紤]推力器陣列單位推力為0.01 N，產(chǎn)生的推力為0.01 N的整數(shù)倍，則分配矩陣有

(34)

取各推力器效率因子初值為1，推力器的真實(shí)效率因子為g=[1.2 1.3 0.9 0.7 0.8 1.1]T，p=6 ，假設(shè)同一個陣列上的微推力器具有相同的效率，且推力偏差初始時(shí)刻即存在，采用本文所提推力估計(jì)方法進(jìn)行在線推力估計(jì)，6個推力器的效率估計(jì)如圖3～4所示。

表1 編隊(duì)衛(wèi)星初始標(biāo)準(zhǔn)軌道要素Table 1 Initial orbital elements of formation satellites

從圖3～4可以看出，在初始時(shí)便出現(xiàn)效率因子偏差的情況下，利用二次規(guī)劃算法可以很好地估計(jì)常值效率因子，6個推力器的效率因子估計(jì)值基本都收斂到真值，即獲得了推力器的真實(shí)推力。

同時(shí)考察從星1的位置誤差變化過程，得到圖5～7所示的位置偏差曲線。其中，理想控制指未出現(xiàn)推力偏差的情況；實(shí)際控制指出現(xiàn)推力偏差但未采用本文方法進(jìn)行估計(jì)補(bǔ)償；本文方法指采用推力估計(jì)與分配方法進(jìn)行推力偏差實(shí)時(shí)估計(jì)與補(bǔ)償。

從圖5～7可以得出以下結(jié)論：

1)利用推力估計(jì)與分配算法，能夠在推力器出現(xiàn)推力偏差的情況下，讓系統(tǒng)的狀態(tài)響應(yīng)接近無推力偏差的情況。

2)不采用該推力估計(jì)與分配算法，由于推力偏差的存在，使得系統(tǒng)動態(tài)響應(yīng)過程與標(biāo)準(zhǔn)響應(yīng)有所偏離，如圖5中穩(wěn)定后曲線完全偏離標(biāo)準(zhǔn)響應(yīng)，而圖6、圖7則出現(xiàn)了高于標(biāo)準(zhǔn)控制曲線的超調(diào)。

3)若K矩陣根據(jù)實(shí)際要求進(jìn)行極點(diǎn)配置或者采用最優(yōu)理論設(shè)計(jì)，那么不采用該分配算法，則可能在實(shí)際使用中出現(xiàn)推力偏差時(shí)喪失最優(yōu)性或者動態(tài)性能，從而出現(xiàn)非預(yù)期的效果。

4)該方法使得推力器在出現(xiàn)推力偏差后系統(tǒng)的動態(tài)、穩(wěn)態(tài)特性仍能保持近似不變，提升微推力器的工作能力，保障實(shí)際控制性能，減輕地面測試中對推力測量的要求，降低測試成本，為微推力器的在軌工程應(yīng)用打下基礎(chǔ)。

6 結(jié) 論

針對實(shí)際使用中微推力器推力輸出無法準(zhǔn)確測量的問題，本文采用一種二次規(guī)劃方法對實(shí)際推力進(jìn)行估計(jì)，并利用混合整數(shù)規(guī)劃算法進(jìn)行微推力器的推力分配，實(shí)現(xiàn)推力偏差的實(shí)時(shí)估計(jì)與補(bǔ)償。仿真結(jié)果表明,該方法能夠很好地估計(jì)常值效率因子，同時(shí)配合推力分配策略，可以保證推力器在出現(xiàn)推力偏差時(shí)系統(tǒng)的狀態(tài)響應(yīng)仍然與未出現(xiàn)推力偏差時(shí)的標(biāo)準(zhǔn)狀態(tài)響應(yīng)近似相同，保障了系統(tǒng)在出現(xiàn)推力偏差時(shí)的控制性能，并證明了采用該分配算法的系統(tǒng)有界穩(wěn)定。

[1] 范林東，楊博，苗峻，等. 基于SiC MEMS陣列的高精度微納衛(wèi)星編隊(duì)保持[J]. 中國空間科學(xué)技術(shù), 2016, 36(2):37-45. [Fan Lin-dong，Yang Bo，Miao Jun，et al. High precision micro-nano satellite formation keeping based on SiC MEMS micro thrust array[J]. Chinese Space Science and Technology, 2016, 36(2):37-45.]

[2] 楊靈芝, 魏延明, 劉旭輝. MEMS固體微推力器陣列發(fā)展研究[J]. 空間控制技術(shù)與應(yīng)用, 2016, 42(1):13-19. [Yang Ling-zhi, Wei Yan-ming, Liu Xu-hui. Development of MEMS solid propellant micro-thruster array[J]. Chinese Space Science and Technology, 2016, 42(1):13-19.]

[3] 劉旭輝，方蜀州，馬紅鵬，等. 基于固體微推力器陣列的衛(wèi)星控制一體化算法[J]. 固體火箭技術(shù), 2012, 35(1):17-23. [Liu Xu-hui，F(xiàn)ang Shu-zhou，Ma Hong-peng，et al. Integrated control algorithm based on solid propellant micro-thruster array[J]. Journal of Solid Rocket Technology, 2012, 35(1):17-23.]

[4] Johansen T A, Fossen T I. Control allocation-a survey[J]. Automatica, 2013, 49(5):1087-1103.

[5] 王敏, 解永春. 考慮推力器推力上界及故障情況的航天器實(shí)時(shí)指令分配最優(yōu)查表法[J]. 宇航學(xué)報(bào), 2010, 31(6):1540-1546. [Wang Min，Xie Yong-chun. Spacecraft thrusters real-time command allocation algorithm in consideration of thrust upper bounds and thruster failure[J]. Journal of Astronautics, 2010, 31(6):1540-1546.]

[6] 唐生勇，張世杰，陳閩，等. 交會對接航天器推力器分配算法研究[J].宇航學(xué)報(bào), 2008, 29(4):1120-1125. [Tang Sheng-yong，Zhang Shi-jie，Chen Min，et al. Research on a thrust allocation algorithm of spacecraft in RVD[J]. Journal of Astronautics， 2008， 29(4):1120-1125.]

[7] 張世杰，段晨陽，趙亞飛. 考慮負(fù)載均衡的過驅(qū)動航天器推力器分配方法[J]. 宇航學(xué)報(bào), 2015, 36(7):826-832. [Zhang Shi-jie, Duan Chen-yang, Zhao Ya-fei. Thruster allocation for over-actuated spacecraft with load balancing[J]. Journal of Astronautics, 2015, 36(7):826-832.]

[8] Bodson M, Frost S A. Load balancing in control allocation[J]. Journal of Guidance Control & Dynamics, 2012, 34(2):380-387.

[9] Hu Q, Li B, Zhang Y. Nonlinear proportional- derivative control incorporating closed-loop control allocation for spacecraft[J]. Journal of Guidance Control & Dynamics, 2014, 37(3):799-812.

[10] Castaldi P, Mimmo N, Simani S. Differential geometry based active fault tolerant control for aircraft [J]. Control Engineering Practice, 2014, 32:227-235.

[11] Alessro C, Emanuele G. Fault‐tolerant adaptive control allocation schemes for over-actuated systems[J]. International Journal of Robust & Nonlinear Control, 2010, 20(17):1958-1980.

[12] Ye D, Park J H, Fan Q. Adaptive robust actuator fault compensation for linear systems using a novel fault estimation mechanism[J]. International Journal of Robust & Nonlinear Control, 2016, 26(8):1597-1614.

[13] Seron M M, De D, Jose A. Robust actuator fault compensation accounting for uncertainty in the fault estimation[J]. International Journal of Adaptive Control & Signal Processing, 2014, 28(12):1440-1453.

[14] Kofman E, Haimovich H, Seron M M. A systematic method to obtain ultimate bounds for perturbed systems[J]. International Journal of Control, 2007, 80(2):167-178.

[15] 張世杰, 段廣仁. 分布式衛(wèi)星編隊(duì)飛行隊(duì)形保持協(xié)同控制[J]. 宇航學(xué)報(bào), 2011, 32(10):2140-2145. [Zhang Shi-jie, Duan Guang-ren. Cooperative control for distributed satellites formation keeping[J]. Journal of Astronautics,2011, 32(10):2140-2145.]