撓性航天器多目標(biāo)魯棒姿態(tài)控制的DPSO算法實(shí)現(xiàn)

王夢菲，張軍,*

1. 北京控制工程研究所, 北京 100190 2. 空間智能控制技術(shù)國家級重點(diǎn)實(shí)驗(yàn)室，北京 100190

由于未建模動態(tài)和建模的不確定性，標(biāo)稱模型下設(shè)計(jì)的控制器在實(shí)際中性能會出現(xiàn)偏差。在軌時航天器常受到內(nèi)外干擾作用，特別是帶有大型撓性子結(jié)構(gòu)的現(xiàn)代航天器，其撓性結(jié)構(gòu)易受干擾激發(fā)產(chǎn)生振動，嚴(yán)重影響姿態(tài)指向精度和穩(wěn)定度，甚至導(dǎo)致航天器失穩(wěn)，如美國的Explore I[1]。因此，大型撓性航天器高精度、高穩(wěn)定度姿態(tài)控制是一個熱點(diǎn)問題也是控制界的難題。

文獻(xiàn)[2-4]研究了經(jīng)典PID控制方法在撓性航天器姿態(tài)控制中的應(yīng)用；當(dāng)模型存在不確定性、未建模模態(tài)發(fā)生振動以及耦合振動發(fā)生變化時，其控制性能下降，且PID控制對外部干擾的魯棒性也有待加強(qiáng)。文獻(xiàn)[5]針對結(jié)構(gòu)振動，基于狀態(tài)反饋設(shè)計(jì)了時間最優(yōu)的LQR控制器，但不適合高頻范圍內(nèi)未建模模態(tài)的振動抑制，且不能滿足模型不確定性時的魯棒控制。文獻(xiàn)[6]改進(jìn)了自適應(yīng)滑?？刂品椒?，最大限度地減少了不確定性、擾動和柔性動態(tài)坐標(biāo)測量等帶來的困難，并在一定程度上改善了滑?？刂乒逃械亩墩瘳F(xiàn)象。文獻(xiàn)[7]證明了基于姿態(tài)角和角速度的反饋對于大范圍攝動的穩(wěn)定性，并針對含參數(shù)不確定性和未建模動態(tài)的姿態(tài)控制系統(tǒng)，設(shè)計(jì)了基于線性矩陣不等式(linear matrix inequality, LMI)的魯棒H∞控制器。文獻(xiàn)[8]針對撓性航天器姿態(tài)機(jī)動問題，以姿態(tài)達(dá)到目標(biāo)角度且保持穩(wěn)定的時間為適應(yīng)度值，將粒子群優(yōu)化(particle swarm optimization, PSO)算法與黃金分割控制等結(jié)合，設(shè)計(jì)了一種離線控制器；相對于常規(guī)的控制方法，該方法的優(yōu)點(diǎn)在于：可根據(jù)系統(tǒng)實(shí)際動力學(xué)特性和控制要求選擇最優(yōu)機(jī)動路徑，實(shí)現(xiàn)姿態(tài)快速機(jī)動和穩(wěn)定。文獻(xiàn)[9-11]將神經(jīng)網(wǎng)絡(luò)、模糊控制理論應(yīng)用于姿態(tài)控制系統(tǒng)設(shè)計(jì)；但模糊控制效果往往依賴于操作者經(jīng)驗(yàn)；為提高可靠性，神經(jīng)網(wǎng)絡(luò)一般需在線學(xué)習(xí)，受星載計(jì)算機(jī)限制。

上述撓性航天器姿態(tài)控制方法大多是針對單一控制目標(biāo)而設(shè)計(jì)的，但在工程實(shí)際中，控制系統(tǒng)的設(shè)計(jì)需考慮魯棒性、快速性、精度和執(zhí)行機(jī)構(gòu)輸出飽和等多目標(biāo)要求，這就屬于典型的多目標(biāo)設(shè)計(jì)問題，關(guān)于此類問題的研究較少。

針對上述多目標(biāo)魯棒控制問題，本文在自適應(yīng)PSO算法中加入變異操作，提出一種差分粒子群優(yōu)化(differential particle swarm optimization, DPSO)算法，以提高粒子的多樣性和搜索性能；同時將該算法與魯棒D-穩(wěn)定極點(diǎn)約束以及Pareto最優(yōu)結(jié)合，涉及的數(shù)學(xué)公式較少，設(shè)計(jì)靈活。相比于傳統(tǒng)帶極點(diǎn)配置的LMI方法，本文方法不受特定指標(biāo)和約束模式的限制，可避免凸優(yōu)化過程，并減小求解的保守性[12]。相對于文獻(xiàn)[13-14]提到的借助特征結(jié)構(gòu)配置的參數(shù)化方法，本文優(yōu)化過程無需將多目標(biāo)轉(zhuǎn)化為一個加權(quán)目標(biāo)函數(shù)，可避免加權(quán)系數(shù)的選擇困難，特別是多目標(biāo)存在矛盾關(guān)系的情況；而且文獻(xiàn)中方法需將指標(biāo)轉(zhuǎn)化為特定的參數(shù)化形式，對參數(shù)攝動敏感度的優(yōu)化無法保證系統(tǒng)魯棒穩(wěn)定。本文方法適用的指標(biāo)形式更廣，能實(shí)現(xiàn)所允許的參數(shù)不確定范圍內(nèi)的魯棒D-穩(wěn)定；可為大型撓性航天器魯棒控制器設(shè)計(jì)提供一種思路。

1 撓性航天器動力學(xué)模型

帶有一對太陽翼的航天器姿態(tài)動力學(xué)方程可表示為[15]：

(1)

本文僅考慮航天器姿態(tài)定向保持和姿態(tài)鎮(zhèn)定情況，假設(shè)姿態(tài)為小角度，此時姿態(tài)運(yùn)動學(xué)為

(2)

Tg=Agθ

(3)

(4)

(5)

考慮模型不確定性，矩陣A、B可表示為[16]：

式中：下標(biāo)0表示為標(biāo)稱系統(tǒng)下的矩陣；D、F1和F2為適當(dāng)維數(shù)的常矩陣；Δ為未知矩陣且滿足ΔTΔ≤I。忽略二階不確定項(xiàng)并只考慮轉(zhuǎn)動慣量的不確定性可得：

(6)

式中：d∈[0,1)，當(dāng)d=0時系統(tǒng)不存在不確定性；

diag[4(J22+J33),3(J11+J33),J11+J22]

當(dāng)不考慮撓性模態(tài)時，e12、e13、e21、e31均為零矩陣。

證明：首先推導(dǎo)矩陣相加求逆結(jié)果，如對于矩陣M有(M+ΔM)-1=Μ-1+N，其中ΔM為不確定項(xiàng)，N為待求解的矩陣。將等式兩邊分別左乘M+ΔM得：

于是

(M+ΔM)-1=Μ-1-(MΔM-1Μ+Μ)-1

同理可推導(dǎo)：

保留ΔM的一次項(xiàng)，近似可得

(M+ΔM)-1=Μ-1-M-1ΔMM-1

那么忽略二階不確定項(xiàng)，考慮轉(zhuǎn)動慣量的不確定性有：

記

(7)

同理可得：

(8)

(9)

由式(7)～(9)可得式(6)。得證。

2 魯棒控制器智能優(yōu)化設(shè)計(jì)

2.1 智能優(yōu)化算法介紹

(1)PSO算法

PSO算法是一種受鳥群覓食啟發(fā)的隨機(jī)搜索算法[17]，具有結(jié)構(gòu)簡單、參數(shù)少和搜索速度快等優(yōu)點(diǎn)，已廣泛應(yīng)用到自然科學(xué)和工程科學(xué)多個領(lǐng)域。該算法中速度和位置更新如下。

式中：k(k=1,2,…,G)為當(dāng)前迭代次數(shù)，G為最大迭代次數(shù)；vi和xi分別為第i個粒子的速度和位置向量，取值區(qū)間分別為[vmin,vmax]和[xmin,xmax]，i∈{1,…,N}，N為群體個體總數(shù)；pi和g分別為個體及全局極值，pbest(i)、gbest為對應(yīng)的適應(yīng)度值；c1、c2為加速常數(shù)即學(xué)習(xí)率；r1、r2為[0，1]范圍隨機(jī)數(shù)；w(k)=wmax-k(wmax-wmin)/G為慣性權(quán)重，取值區(qū)間為[wmin,wmax]。

(2)差分進(jìn)化(differential evolution, DE)算法

DE算法是一種魯棒性較強(qiáng)的隨機(jī)搜索算法，其基本思想來源于遺傳算法，這里選取差分進(jìn)化算法DE/rand/1/bin策略，種群參數(shù)定義同上述PSO算法，關(guān)鍵步驟如下[18-19]。

1)變異,計(jì)算變異向量:

υi(k+1)=xl1(k)+F(k)·[xl2(k)-xl3(k)]

2)交叉，計(jì)算試驗(yàn)變量：

式中：交叉算子CR=0.3×[1+rand(0,1)]；randb(j)表示產(chǎn)生[0,1]之間隨機(jī)數(shù)發(fā)生器的第j個估計(jì)值；rnbr(i)∈(1,2,…,D)表示一個隨機(jī)選擇的序列，以確保ui至少從υi中獲得一個參數(shù)。

3)選擇，DE算法按照貪婪準(zhǔn)則將試驗(yàn)向量ui(k+1)與當(dāng)前種群中的目標(biāo)向量xi(k)進(jìn)行比較，較優(yōu)的向量將在下一代種群中出現(xiàn)。

(3)DPSO算法

DPSO算法在PSO算法的基礎(chǔ)上加入DE算法的變異操作，在保證粒子群朝著好的方向發(fā)展的同時增加了粒子的多樣性，提高粒子的全局搜索能力，有效避免了早熟至局部最優(yōu)解的情況[20]。算法流程和部分參數(shù)定義同PSO算法，其速度和位置更新如下。

(10)

F(k)定義同DE算法，r3為[0，1]范圍內(nèi)的隨機(jī)數(shù)，xU(k)和xL(k)為當(dāng)代粒子中的較優(yōu)和較劣解，以保證粒子向有益的方向發(fā)展。

2.2 多目標(biāo)約束描述

由于撓性模態(tài)難以觀測或直接測量，本節(jié)參考文獻(xiàn)[7]，采用姿態(tài)角和角速度反饋設(shè)計(jì)控制器，下面關(guān)于多目標(biāo)約束的描述建立在此基礎(chǔ)上。

考慮轉(zhuǎn)動慣量的不確定性和未建模動態(tài)、各種干擾對系統(tǒng)的影響，以及避免控制飽和的現(xiàn)象，本節(jié)在控制器設(shè)計(jì)時加入魯棒D-穩(wěn)定約束，并對控制能量和內(nèi)外干擾抑制相關(guān)指標(biāo)進(jìn)行優(yōu)化。極點(diǎn)約束和兩個典型優(yōu)化指標(biāo)的顯示表達(dá)以及Pareto最優(yōu)理論描述如下。

2.2.1 魯棒D-穩(wěn)定

(1)含撓性模態(tài)

設(shè)穩(wěn)定裕度參數(shù)為h(h<0)，那么式(5)系統(tǒng)對應(yīng)如圖1所示，魯棒D-穩(wěn)定實(shí)現(xiàn)為：

(11)

式中：K∈Rm×n為相應(yīng)維數(shù)的反饋增益矩陣，在優(yōu)化的過程中視為已知，此時控制輸入u=BKy；ε>0為標(biāo)量。優(yōu)化時可將標(biāo)稱系統(tǒng)下的撓性模態(tài)極點(diǎn)限制在一定變化范圍內(nèi)，如實(shí)部區(qū)間為[a1,a2](a1,a2<0)，虛部絕對值區(qū)間為[b1,b2](b1,b2>0)。

圖1 極點(diǎn)配置區(qū)域(1)Fig.1 Pole assignment region(1)

證明：參考文獻(xiàn)[16]，考慮不確定性和保證系統(tǒng)一定的魯棒穩(wěn)定裕度有：

P[(A0+ΔA)+(B0+ΔB)KC]+

[(A0+ΔA)+(B0+ΔB)KC]TP-hI<0

式中：P為對稱正定矩陣。記Y=P(A0+B0KC)+(A0+B0KC)TP-hI，則有

Y+PDΔ(F1+F2K)+

(F1+F2K)TΔT(PD)T<0

上述矩陣不等式對所有滿足ΔTΔ≤I的不確定矩陣Δ成立，當(dāng)且僅當(dāng)存在ε，使得

Y+εPDDTP+ε-1(F1+F2K)T(F1+F2K)<0

應(yīng)用矩陣的Shur補(bǔ)性質(zhì)，可以等價為：

上式左右同乘矩陣diag{P-1,I}，并記X1=P-1，便可得到式(11)，得證。

(2)不含撓性模態(tài)

針對不確定系統(tǒng)，對給定區(qū)域D，若存在對稱矩陣X1、X2和標(biāo)量ε1、ε2>0同時滿足：

圖2 極點(diǎn)配置區(qū)域(2)Fig.2 Pole assignment region(2)

(12)

(13)

式中：A0、B0、C不含撓性模態(tài)信息且此時C為是單位陣；此處X1可不等于X2，ε1可不等于ε2，以減小保守性。則該不確定系統(tǒng)為魯棒D-穩(wěn)定。證明過程同式(11)。

2.2.2 (內(nèi)外)干擾抑制能力

式中：Ac=A+BKC。建立從w到y(tǒng)的傳遞函數(shù)矩陣Gyw(s)，那么根據(jù)H2控制理論思想，用‖Gyw(s)‖2定義干擾對輸出影響的大小[13]：

2.2.3 控制能量

現(xiàn)代航天器多以動量輪或控制力矩陀螺為姿態(tài)控制執(zhí)行機(jī)構(gòu)，為防止因執(zhí)行機(jī)構(gòu)飽和而導(dǎo)致的系統(tǒng)性能變差或失穩(wěn)現(xiàn)象，需加入控制能量約束。直接將飽和非線性因素加入控制器設(shè)計(jì)會帶來分析困難，所以參考文獻(xiàn)[13]，基于航天器姿態(tài)定向保持時姿態(tài)角偏差和姿態(tài)角速度均為小量的特點(diǎn)，通過優(yōu)化‖K‖2的值來盡可能減少控制能量‖u‖。記

f2(K)=‖K‖2

2.2.4 Pareto最優(yōu)原則

一般優(yōu)化指標(biāo)之間的關(guān)系難以確定，且不乏相互沖突的情況，因此無法找到多目標(biāo)均為最優(yōu)的唯一解，只能找到在多目標(biāo)間平衡的次優(yōu)解，此時需應(yīng)用Pareto最優(yōu)解的概念。假設(shè)存在性能指標(biāo)函數(shù)f1(x),…，fs(x)(s≥2)，其中x為搜索空間φ中的向量，與K對應(yīng)，有如下定義：

x=[K(1,1),…,K(1,n),K(2,1),…,

K(2,n),…,K(m,n)]

Pareto最優(yōu)解：如果不存在同時滿足下式的x∈φ，那么x*∈φ則為Pareto最優(yōu)解或非劣解(此解不唯一)[21]：

(14)

(15)

如果存在滿足上述條件的解x，那么則稱該解Pareto占優(yōu)于解x*。

創(chuàng)建由粒子每代更新產(chǎn)生的非劣解和占優(yōu)解組合而成的外部解集；因?yàn)橄氲玫矫總€指標(biāo)都相對較優(yōu)而不是某單個指標(biāo)最優(yōu)的解，所以為了提高算法的快速性，不按擁擠度法[12]而是按距離評價法更新該解集，保留距離較小的解。做歸一化處理，對應(yīng)x的距離d(x)定義如下：

式中：maxfs和minfs分別表示第s個目標(biāo)函數(shù)在外部解集中取得的最大和最小值。以s=2為例，如圖3所示。

圖3 距離評價法Fig.3 Distance evaluation

2.3 多目標(biāo)魯棒控制器設(shè)計(jì)流程

基于DPSO算法的控制器設(shè)計(jì)步驟如下。

步驟1：初始化粒子群。包括算法參數(shù)和位置、速度向量，其中初始位置向量由PID控制下的增益矩陣和LMI方法下求得的幾組反饋增益矩陣定義，初始速度向量為零向量；初始化粒子i(i=1,2,…，N)的歷史最優(yōu)位置pi和最優(yōu)值p1best(i)和p2best(i)并選取第一個粒子為全局最優(yōu)粒子，對應(yīng)位置g和最優(yōu)值g1best和g2best。迭代次數(shù)k=1，且xU(k)=xL(k)=x1(k)。

步驟2：更新粒子速度和位置。k=k+1；按式(10)更新速度和位置，并做邊界條件處理；判斷是否滿足魯棒-D穩(wěn)定極點(diǎn)約束，若滿足進(jìn)入步驟3，否則進(jìn)入步驟6。

步驟3：更新個體粒子最優(yōu)解和全局最優(yōu)解。計(jì)算個體粒子適應(yīng)度值(性能指標(biāo))f1(xi),f2(xi)，按Pareto占優(yōu)更新最優(yōu)解；將每個粒子當(dāng)代更新的非劣解和占優(yōu)解存入外部解集。

步驟4：更新外部解集。計(jì)算解集中每個解對應(yīng)的距離，找到距離最短的解，然后將其他解對應(yīng)的適應(yīng)度值與該解比較，保留下非劣解和占優(yōu)解；判斷保留下的解的個數(shù)是否大于集合上限N，若是，進(jìn)入步驟5，否則剔除距離較大的部分解。

步驟5：計(jì)算xU(k)和xL(k)。計(jì)算各粒子的適應(yīng)度值f1(xi),f2(xi)，并求取所有粒子適應(yīng)度值的平均值；將每個粒子對應(yīng)的適應(yīng)度值與平均值比較，在Pareto占優(yōu)的粒子中任選一個將其解作為xU(k)，在剩余粒子中任選一個將其解作為xL(k)；當(dāng)不存在Pareto占優(yōu)的粒子時，xU(k)=xL(k)，取任意粒子位置向量。

步驟6：判斷是否滿足k≤G，如滿足返回步驟2，否則進(jìn)入步驟7。

步驟7：畫圖，并定義仿真參數(shù)。

3 仿真校驗(yàn)

以某大型航天器為例，取其一對翼板的前兩階撓性模態(tài)進(jìn)行控制器設(shè)計(jì)和仿真。標(biāo)稱系統(tǒng)轉(zhuǎn)動慣量、耦合系數(shù)、帆板撓性約束模態(tài)角頻率對角陣和阻尼系數(shù)分別為：

圖2中參數(shù)α=50°，h1=-0.2，h2=-0.02；撓性模態(tài)極點(diǎn)相關(guān)參數(shù)h=-0.001，a1=-0.03，a2=-0.001 5，b1=0.25，b2=0.35。優(yōu)化算法參數(shù)取值如表1所示；以控制能量為適應(yīng)度值為例，PSO、DE和DPSO三種算法的對比優(yōu)化結(jié)果見圖4。

從圖4可看出DE算法進(jìn)化較慢，PSO算法過早收斂，容易陷入局部極值，本文提出的DPSO算法尋優(yōu)能力最好，因此采用DPSO算法進(jìn)行多目標(biāo)優(yōu)化。圖5為基于Pareto最優(yōu)原則的多目標(biāo)優(yōu)化結(jié)果，可以看出在距離評價法下，外部解集中的非劣解集合接近Pareto前沿[12]。

表1 三種算法下參數(shù)設(shè)置

圖4 三種算法下進(jìn)化曲線Fig.4 Evolution curve under three algorithms

圖5 外部解集Fig.5 Set of external solutions

為驗(yàn)證本文控制器的魯棒性能、干擾抑制能力和控制能量優(yōu)化效果，仿真時不妨取轉(zhuǎn)動慣量(不同于J)和外部干擾為：

圖8所示為三軸姿態(tài)受到的干擾作用，圖9～圖11為閉環(huán)系統(tǒng)仿真結(jié)果。由圖9知航天器姿態(tài)穩(wěn)定時間相對于PID控制減小，基本不超過100 s；在300 s后，PID控制姿態(tài)誤差在1.35×10-4(°)以內(nèi)，而本文方法下姿態(tài)誤差在6.2×10-5(°)以內(nèi)，減小了約54%，說明對干擾進(jìn)行了有效抑制；圖10顯示在本文仿真算例下，兩種方法控制力矩相差不大，本文方法控制力矩峰值約為6.99N·m，PID控制峰值約為8.42N·m；由圖11知，本文控制方法下的航天器前兩階撓性模態(tài)振動幅度與PID控制相當(dāng)，且處于較小的量級，其中與姿態(tài)耦合作用較強(qiáng)的第一階撓性模態(tài)振動衰減更快，在100 s左右接近穩(wěn)定(比PID快約250 s)。綜上，本文方法對航天器高精度高穩(wěn)定度魯棒姿態(tài)控制具有良好的效果。

圖 6 極點(diǎn)配置結(jié)果(圖中線條與圖2對應(yīng))Fig.6 Result of pole assignment

圖7 開閉環(huán)系統(tǒng)奇異值曲線Fig.7 Singular value of open-loop and closed-loop system

圖8 姿態(tài)干擾Fig.8 Attitude interference

圖9 姿態(tài)角Fig.9 Attitude angle

圖10 控制力矩Fig10 Control torque

圖11 撓性模態(tài)位移Fig.11 Flexible modal displacement

4 結(jié)束語

本文針對大型撓性航天器的高精度高穩(wěn)定度姿態(tài)控制問題，將改進(jìn)智能算法、現(xiàn)代控制理論和多目標(biāo)Pareto最優(yōu)思想結(jié)合，提出一種基于輸出反饋和DPSO算法的多目標(biāo)魯棒控制器設(shè)計(jì)方法，并推導(dǎo)了含轉(zhuǎn)動慣量不確定性的系統(tǒng)模型和魯棒D-穩(wěn)定的LMI表達(dá)。主要結(jié)論如下：

1)相較于DE算法和PSO算法，本文提出的DPSO算法尋優(yōu)效果更好。在魯棒區(qū)域極點(diǎn)約束下，利用DPSO算法和Pareto最優(yōu)原則，對控制能量和魯棒性能進(jìn)行優(yōu)化得到的控制器實(shí)現(xiàn)了閉環(huán)系統(tǒng)的魯棒D-穩(wěn)定，且起到了一定的振動抑制作用。

2)仿真驗(yàn)證了本文方法的有效性。在本算例中，相對于PID控制，本文方法下系統(tǒng)響應(yīng)速度較快，其中第1階撓性模態(tài)振動衰減時間減少約250 s，控制力矩峰值減小約1.43N·m，姿態(tài)穩(wěn)態(tài)誤差減小約54%。

3)將智能優(yōu)化與現(xiàn)代控制理論結(jié)合，涉及的數(shù)學(xué)公式和參數(shù)較少，靈活性較強(qiáng)。本文方法彌補(bǔ)了經(jīng)典PID控制在不確定系統(tǒng)模型下魯棒性能以及抗干擾能力的不足；減小了帶極點(diǎn)配置的LMI方法在多指標(biāo)約束凸優(yōu)化問題中的求解保守性；也能避免傳統(tǒng)方法將多目標(biāo)轉(zhuǎn)化為加權(quán)指標(biāo)函數(shù)時，由于多目標(biāo)關(guān)系難以確定而導(dǎo)致的加權(quán)系數(shù)選擇困難。

4)本文智能優(yōu)化方法還可應(yīng)用到更多復(fù)雜以及非線性系統(tǒng)控制參數(shù)優(yōu)化問題上，為大型撓性航天器魯棒控制器設(shè)計(jì)提供一種思路。