基于可靠性優(yōu)化的固體火箭姿態(tài)控制器設(shè)計(jì)

嚴(yán) 愷 吳燕生 張 兵 鐘 震

1. 北京宇航系統(tǒng)工程研究所，北京 100076；2. 中國(guó)航天科技集團(tuán)公司，北京 100048

0 引言

在固體火箭的姿態(tài)控制器設(shè)計(jì)過程中，質(zhì)量偏差、發(fā)動(dòng)機(jī)推力偏差、氣動(dòng)偏差等諸多不確定性因素普遍存在。為了保證火箭在實(shí)際飛行中的穩(wěn)定性和性能，必須在設(shè)計(jì)中對(duì)這些偏差量予以處理。而隨著長(zhǎng)細(xì)比的增加和整流罩尺寸的變大，火箭姿態(tài)控制系統(tǒng)面臨的參數(shù)不確定問題也會(huì)加劇[1]。當(dāng)前工程中最常用的處理方式為根據(jù)經(jīng)驗(yàn)對(duì)這些偏差量進(jìn)行極限拉偏組合，確定系統(tǒng)的極限狀態(tài)，并按照最不利的工作條件進(jìn)行控制器的分析和設(shè)計(jì)[2]。實(shí)際上，這些極端情況出現(xiàn)的概率一般非常低，但為了保證系統(tǒng)在這些極端情況下的性能，需要的代價(jià)卻很大，從而導(dǎo)致設(shè)計(jì)的極度保守，產(chǎn)品效率低下。與之相比，概率設(shè)計(jì)方法充分考慮到不確定參數(shù)的概率分布情況，并以可靠性的方式提出性能要求，得到更加符合工程實(shí)際需求的設(shè)計(jì)，降低了設(shè)計(jì)的保守性[3]。

可靠性(Reliability)為系統(tǒng)指標(biāo)滿足的概率，對(duì)系統(tǒng)可靠性的估計(jì)稱為可靠性分析(Reliability Analyze, RA)，與之對(duì)應(yīng)的控制器設(shè)計(jì)方法稱為基于可靠性的控制器設(shè)計(jì)(Reliability-Based Controller Design, RBCD)[4]。RBCD包含概率估計(jì)和參數(shù)設(shè)計(jì)2個(gè)層次，整體構(gòu)成一個(gè)優(yōu)化問題，求解目標(biāo)為在可衡量的概率水平下獲得一個(gè)最優(yōu)設(shè)計(jì)[5]。文獻(xiàn)[6-8]等將各個(gè)指標(biāo)滿足的概率的加權(quán)和作為優(yōu)化的目標(biāo)函數(shù)，并對(duì)控制器參數(shù)進(jìn)行尋優(yōu)設(shè)計(jì)。這類方法中的可靠性分析通過蒙特卡羅采樣法實(shí)現(xiàn)，整體計(jì)算量極大。此外，指標(biāo)概率加權(quán)和形式的目標(biāo)函數(shù)不能滿足每一單個(gè)指標(biāo)的概率要求。另一類方法將系統(tǒng)參數(shù)的設(shè)計(jì)問題描述為一個(gè)有約束的優(yōu)化問題，而約束為系統(tǒng)的可靠性要求，構(gòu)成可靠性優(yōu)化問題。文獻(xiàn)[9-10]進(jìn)行了魯棒控制器的設(shè)計(jì)。這種方法主要應(yīng)用于結(jié)構(gòu)設(shè)計(jì)領(lǐng)域以及總體參數(shù)設(shè)計(jì)等領(lǐng)域[11],關(guān)于其求解效率和精度的改進(jìn)研究較多，但在控制器設(shè)計(jì)方面的應(yīng)用卻相對(duì)較少。

可靠性分析的方法和優(yōu)化與可靠性分析關(guān)系的處理是可靠性優(yōu)化問題的2個(gè)主要問題。在可靠性分析方法中，性能測(cè)度法[12](Performance Measure Approach, PMA)由于計(jì)算效率高、收斂性好而在結(jié)構(gòu)設(shè)計(jì)中得到了廣泛的應(yīng)用。對(duì)于可靠性分析問題與優(yōu)化問題關(guān)系的處理，序列解耦與可靠性評(píng)估(Sequential Optimization and Reliability Assessment, SORA)算法[13]通過將可靠性分析過程與參數(shù)尋優(yōu)過程解耦，將嵌套結(jié)構(gòu)轉(zhuǎn)化為順序結(jié)構(gòu)，進(jìn)而提高問題的求解效率。當(dāng)前，RBCD領(lǐng)域鮮有文獻(xiàn)采用這種高效求解算法進(jìn)行固體火箭的控制器參數(shù)設(shè)計(jì)，因此有必要對(duì)此進(jìn)行研究。

本文針對(duì)固體火箭的姿態(tài)控制問題，采用概率優(yōu)化方法設(shè)計(jì)火箭的控制器參數(shù)。首先建立火箭的縱向通道線性化模型，并使用概率描述控制系統(tǒng)的性能指標(biāo)；隨后，采用高效的PMA方法，近似估計(jì)指標(biāo)滿足的概率，采用SORA方法，解耦參數(shù)尋優(yōu)過程，極大地提高優(yōu)化問題的求解效率；最后，設(shè)計(jì)火箭縱向通道模型的PD控制器參數(shù)，并與傳統(tǒng)極限拉偏方法的設(shè)計(jì)結(jié)果進(jìn)行對(duì)比分析。結(jié)果表明，采用概率優(yōu)化方法設(shè)計(jì)的控制器能夠在保證給定的系統(tǒng)可靠性的基礎(chǔ)上，降低系統(tǒng)設(shè)計(jì)的保守性。

1 固體火箭姿態(tài)控制數(shù)學(xué)描述

1.1 火箭的動(dòng)力學(xué)模型

由于固體火箭的動(dòng)力學(xué)模型為一高度非線性時(shí)變微分方程組，因此在進(jìn)行控制器設(shè)計(jì)時(shí)一般選取典型特征秒進(jìn)行系數(shù)凍結(jié)，并采用小擾動(dòng)假設(shè)得到線性化的微分方程組。不失一般性，本文以火箭俯仰通道為研究對(duì)象，忽略發(fā)動(dòng)機(jī)擺管慣性和箭體彈性，得到線性化的剛體小偏差運(yùn)動(dòng)方程：

(1)

式中：θ,ω,φ和δφ分別為速度傾角、俯仰角速度、俯仰角和噴管等效擺角，c1～c3,b1～b3為剛體系數(shù)，表達(dá)式為

(2)

ci～N(μci,σci),
bi～N(μbi,σbi)
i=1,2,3

(3)

1.2 性能指標(biāo)

目前采用可靠性優(yōu)化方法進(jìn)行控制器設(shè)計(jì)的文獻(xiàn)大都采用了LMI來描述系統(tǒng)的穩(wěn)定性要求，主要是由于LMI具有對(duì)設(shè)計(jì)參數(shù)為凸的良好性質(zhì)，方便優(yōu)化計(jì)算[14]。然而，LMI描述的假設(shè)條件過強(qiáng)，且只能描述系統(tǒng)的穩(wěn)定性條件，因此應(yīng)用性較差。在實(shí)際工程中，控制系統(tǒng)分析與設(shè)計(jì)主要在頻域進(jìn)行，除穩(wěn)定性外，更關(guān)心的是系統(tǒng)的穩(wěn)定裕度。主要包括幅值裕度(Gain Margin, GM)和相位裕度(Phase Margin, PM)，用以表征系統(tǒng)的相對(duì)穩(wěn)定性和動(dòng)態(tài)性能指標(biāo)。為此，本文將采用GM和PM作為系統(tǒng)的性能指標(biāo)，計(jì)算式為

(4)

其中GM,PM分別表示GM和PM。

2 概率分析與設(shè)計(jì)方法

2.1 可靠性分析

相比確定性設(shè)計(jì)，概率設(shè)計(jì)的特點(diǎn)是將系統(tǒng)模型中的不確定性參數(shù)描述為滿足一定概率分布的隨機(jī)變量，系統(tǒng)的性能指標(biāo)要求也由確定性要求放寬為概率要求，如式(5)所示

g(x,k)>0?P(g(x,k)>0)=R≥R1

(5)

式中x,k分別表示系統(tǒng)的不確定參數(shù)和確定性設(shè)計(jì)變量；g(x,k)表示系統(tǒng)性能的指標(biāo)函數(shù)，稱為性能函數(shù)(Performance Function)；R為系統(tǒng)各指標(biāo)滿足相應(yīng)要求的概率，即可靠性；R1表示對(duì)系統(tǒng)指標(biāo)滿足的概率要求，當(dāng)R1=1時(shí)，右式與左式等價(jià)，概率約束退化為確定性約束。根據(jù)概率理論，式(5)中的概率計(jì)算表達(dá)式為

(6)

P的值需要通過對(duì)概率密度函數(shù)fg的多維積分才能進(jìn)行準(zhǔn)確的求取。而實(shí)際上，由于函數(shù)g(·,·)為各隨機(jī)變量的非線性組合，f的具體形式一般很難得到，而且多維積分本身也需要較高的計(jì)算成本，因此一般采用近似方法進(jìn)行估計(jì)，目前應(yīng)用最廣泛的方法為以FORM為代表的可靠性分析方法。

(7)

此時(shí)概率Pi可表示為

(8)

上式中導(dǎo)致積分難以求解的部分主要為當(dāng)指標(biāo)函數(shù)非線性較強(qiáng)時(shí)，滿足g>0的積分域難以確定。為此FORM法將g′(u)在u*處進(jìn)行一階Taylor展開，得到

(9)

當(dāng)用此u的線性函數(shù)近似指標(biāo)函數(shù)時(shí)，由概率論知識(shí)有

(10)

由于此時(shí)g′(u)服從正態(tài)分布，因此有

(11)

由于u*為極限狀態(tài)方程上距離原點(diǎn)最近的點(diǎn)，則可得

Pi=Φ(β*)=Φ(‖u*‖)

(12)

進(jìn)而得到概率Pi的近似解。FORM法中的核心步驟為MPP點(diǎn)的搜索。上述流程即為可靠性指數(shù)法(RIA)，其中由MPP點(diǎn)求得的β*稱為可靠性指數(shù)，表征可靠性的大小。RIA法進(jìn)行可靠性分析可以描述為如下形式

(13)

其中βt表示設(shè)計(jì)要求的可靠性對(duì)應(yīng)的可靠性指數(shù)，當(dāng)MPP點(diǎn)求得的β*>βt時(shí)即認(rèn)為可靠性滿足要求，如圖1所示。

圖1 RIA法與PMA法示意圖

式(13)實(shí)際上為一個(gè)優(yōu)化問題，其目標(biāo)函數(shù)形式簡(jiǎn)單，約束卻很復(fù)雜，因此導(dǎo)致優(yōu)化過程可能出現(xiàn)不穩(wěn)定問題，求解效率較低。PMA法是與RIA法互為逆問題的一種方法，其基本思想為在βt圓上進(jìn)行MPP點(diǎn)搜索，搜索得到的逆MPP點(diǎn)(IMPP)如果滿足g′>0，則說明指標(biāo)滿足區(qū)域存在滿足指定可靠性要求的點(diǎn)，反之則不存在，如圖1所示。具體求解過程可以總結(jié)如下

min:g
s.t.:‖u‖=βt,g′*(uMPP)>0

(14)

PMA法同樣為一個(gè)優(yōu)化過程，但其約束形式簡(jiǎn)單，目標(biāo)函數(shù)復(fù)雜，因此收斂性更好，本文將采用PMA法進(jìn)行系統(tǒng)的可靠性分析，以作為概率優(yōu)化的基礎(chǔ)。

2.2 可靠性優(yōu)化設(shè)計(jì)方法

可靠性優(yōu)化問題本質(zhì)上為一個(gè)以概率要求為約束的優(yōu)化問題，可以描述為以下形式

搜索:k
min:J(k)

s.t.:Pi(gi(x,k)>0)≥Ri,i=1,2,…,N

(15)

其中J(k)為優(yōu)化目標(biāo)函數(shù)，k為設(shè)計(jì)參數(shù)?，F(xiàn)有的RBCD文獻(xiàn)大都采用了雙循環(huán)優(yōu)化結(jié)構(gòu)，即將可靠性分析過程直接嵌入到外層參數(shù)尋優(yōu)過程中，如圖2所示。由2.1節(jié)可知，基于MPP點(diǎn)的可靠性分析方法本質(zhì)上也是一個(gè)優(yōu)化過程，因此雙循環(huán)法會(huì)使兩個(gè)優(yōu)化過程的耦合嵌套，從而導(dǎo)致計(jì)算量劇增，降低了該方法的實(shí)用性。

圖2 雙循環(huán)結(jié)構(gòu)RBCD流程

為了解決雙循環(huán)法的缺陷，本文采用SORA法提高求解效率，并針對(duì)控制器參數(shù)設(shè)計(jì)問題的特殊性對(duì)其進(jìn)行調(diào)整。SORA方法的基本思想為將優(yōu)化結(jié)構(gòu)解耦，將概率約束轉(zhuǎn)化為確定性約束，采用PMA算法在指定的可靠度水平進(jìn)行可靠性分析，通過移動(dòng)確定性約束的邊界，使之與等效的概率約束邊界重合，達(dá)到優(yōu)化問題和概率約束的同步。根據(jù)以上的描述，式(15)可以改寫為以下形式

搜索:k

min:J(k)

s.t.:gi(xMPPi,k)>0,i=1,2,…,N

(16)

其中xMPPi為上一步可靠性分析得到的MPP點(diǎn)。在PMA方法中，可以由MPP點(diǎn)處指標(biāo)函數(shù)值的符號(hào)等價(jià)地描述指標(biāo)的可靠性要求，通過式(14)求取MPP點(diǎn)。以二維隨機(jī)變量為例，原空間的概率約束與確定性約束的等價(jià)表示如圖3所示。

圖3 概率約束的等價(jià)描述示意圖

圖中，曲線A表示確定性約束邊界g(x,k)=0，由于x是隨機(jī)變量，因此均值位于此邊界線上的設(shè)計(jì)可靠性一般較低。曲線B表示概率約束邊界P[g(x,k)]=R，即均值位于曲線B上的設(shè)計(jì)點(diǎn)，其MPP點(diǎn)恰好位于曲線A上。由此可見，概率約束的確定性等價(jià)描述為gi(xMPPi,k)>0，即設(shè)計(jì)點(diǎn)的MPP點(diǎn)須位于可行區(qū)域g(x,k)>0中。

在經(jīng)典的SORA法中，隨機(jī)變量x的均值同樣為設(shè)計(jì)變量，對(duì)應(yīng)于結(jié)構(gòu)設(shè)計(jì)中設(shè)計(jì)參數(shù)與制造得到的實(shí)際參數(shù)亦會(huì)存在偏差的事實(shí)。在控制器參數(shù)設(shè)計(jì)中，隨機(jī)變量的均值是提前給定的標(biāo)稱參數(shù)，設(shè)計(jì)變量只有確定性參數(shù)即控制器的增益k。然而這種情況下SORA法同樣適用，只需取消平移隨機(jī)變量的步驟，而通過調(diào)整確定性設(shè)計(jì)變量k的方式對(duì)約束邊界予以調(diào)整，對(duì)應(yīng)的設(shè)計(jì)流程如圖4所示。

圖4 SORA法設(shè)計(jì)流程圖

3 基于可靠性優(yōu)化的控制器參數(shù)設(shè)計(jì)

根據(jù)式(1)所示的火箭動(dòng)力學(xué)模型，得到火箭俯仰通道的開環(huán)傳遞函數(shù)為

(17)

控制律為

δ(s)=(kp+kds)Δφ

(18)

控制系統(tǒng)框圖表示為

圖5 控制系統(tǒng)框圖

系統(tǒng)的開環(huán)傳遞函數(shù)為

L(s)=(kp+kds)Go

(19)

則設(shè)計(jì)變量矢量和隨機(jī)變量矢量為

(20)

系統(tǒng)指標(biāo)選取為GM和PM，構(gòu)造性能函數(shù)

(21)

其中PMd,GMd分別為要求的幅值裕度和相位裕度。在相同的姿態(tài)角偏差下，控制系統(tǒng)的靜態(tài)增益kp與指令發(fā)動(dòng)機(jī)擺角的輸出量直接相關(guān)，表征了對(duì)執(zhí)行機(jī)構(gòu)的能力需求。kp越大，在同等姿態(tài)角偏差下對(duì)擺管的擺動(dòng)能力需求也就越高，且過大kp值對(duì)箭體彈性振動(dòng)的穩(wěn)定不利。因此，在滿足控制指標(biāo)的前提下，較小的kp值可以降低控制系統(tǒng)的能力需求。為此優(yōu)化目標(biāo)函數(shù)可以選取為靜態(tài)增益的幅值|kp|。此時(shí)控制器參數(shù)的概率優(yōu)化模型為

搜索:kp,kd
min:|kp|

s.t.:

P(g1(x,k)>0)≥R1,
P(g2(x,k)>0)≥R1，
kpl≤kp≤kpu,kdl≤kd≤kdu

(22)

4 仿真校驗(yàn)

本文研究對(duì)象為固體火箭的一級(jí)飛行段，選取最大動(dòng)壓點(diǎn)這一典型的特征秒，結(jié)合參數(shù)的散布情況，所研究模型的剛體系數(shù)標(biāo)稱值和方差如表1所示。

表1 動(dòng)力學(xué)模型參數(shù)標(biāo)稱值

系統(tǒng)的性能指標(biāo)要求為

(23)

由2.1節(jié)可知，極限拉偏包絡(luò)設(shè)計(jì)問題可以視為式(22)中R1=1的概率設(shè)計(jì)，此時(shí)要求系統(tǒng)的最壞情況指標(biāo)滿足性能要求。進(jìn)行SORA優(yōu)化，得到此時(shí)控制器參數(shù)k1=[1.89,0.48]，對(duì)設(shè)計(jì)結(jié)果進(jìn)行MCS仿真驗(yàn)證，仿真500次，得到系統(tǒng)上下極限拉偏對(duì)應(yīng)的性能指標(biāo)，5%、95%分位線對(duì)應(yīng)的性能指標(biāo)，和此時(shí)系統(tǒng)的可靠性R如表2所示。

表2 極限狀態(tài)設(shè)計(jì)的系統(tǒng)指標(biāo)

由計(jì)算結(jié)果可以看出，相位裕度遠(yuǎn)大于期望值，因此幅值裕度為主要考慮的指標(biāo)。按極限拉偏設(shè)計(jì)的參數(shù)保證了系統(tǒng)性能指標(biāo)的下限滿足給定的指標(biāo)要求，因此保證了理論上的可靠性為1。同時(shí)，可以發(fā)現(xiàn)幅值裕度小于9.14dB情況僅占所有情況的5%。此時(shí)系統(tǒng)的開環(huán)Bode圖及其散布情況如圖6所示。

圖6 極限設(shè)計(jì)對(duì)應(yīng)Bode圖

由圖6可以看出，絕大部分情況對(duì)應(yīng)的曲線位于5%和95%分位線定義的包絡(luò)以內(nèi)，而按極限拉偏處理產(chǎn)生的包絡(luò)遠(yuǎn)大于實(shí)際可能出現(xiàn)的范圍。從物理意義上講，這是由于極限包絡(luò)是6個(gè)隨機(jī)參數(shù)同時(shí)進(jìn)行3σ拉偏，且需要滿足特定的拉偏方向的情況下得到的，因此同時(shí)出現(xiàn)的概率極低。這充分說明了極限設(shè)計(jì)帶來的保守性。

作為比較，采用概率要求進(jìn)行設(shè)計(jì)。設(shè)要求的系統(tǒng)可靠性為R1=0.95，即系統(tǒng)在隨機(jī)不確定下滿足指標(biāo)要求概率至少為95%。按此標(biāo)準(zhǔn)進(jìn)行優(yōu)化設(shè)計(jì)得到的控制器參數(shù)為k1=[1.58,0.5]，此時(shí)系統(tǒng)的各指標(biāo)如表3所示。

表3 95%可靠性要求設(shè)計(jì)的系統(tǒng)指標(biāo)

由計(jì)算結(jié)果可以發(fā)現(xiàn)，雖然此時(shí)系統(tǒng)的下極限對(duì)應(yīng)的幅值裕度為6.54dB且小于給定的指標(biāo)要求，但概率優(yōu)化設(shè)計(jì)的結(jié)果保證了幅值裕度小于7.83dB的情況占比不超過5%，整體的可靠性為0.959，系統(tǒng)的性能滿足提前給定的概率要求。此時(shí)系統(tǒng)Bode圖的散布情況如圖7所示。

圖7 95%可靠性設(shè)計(jì)對(duì)應(yīng)Bode圖

由圖可知，隨機(jī)分布的參數(shù)對(duì)應(yīng)的Bode圖的整體散布仍然主要位于5%和95%分位線以內(nèi)，概率設(shè)計(jì)能夠滿足預(yù)期要求。風(fēng)攻角為αwp=5°,αwq=0.8°時(shí)，計(jì)算兩種設(shè)計(jì)下的擺角，得到

(24)

由上述仿真計(jì)算發(fā)現(xiàn)，當(dāng)把系統(tǒng)的可靠性要求降低5%時(shí)，兩種方法的控制系統(tǒng)靜態(tài)增益kp減小了20%，發(fā)動(dòng)機(jī)擺角減小了15%，顯著降低了控制系統(tǒng)的性能需求。

5 結(jié)論

本文對(duì)火箭俯仰通道的不確定參數(shù)進(jìn)行了概率描述，基于可靠性分析理論和高效的SORA方法，針對(duì)幅值裕度和相位裕度指標(biāo)，對(duì)控制器參數(shù)進(jìn)行了概率優(yōu)化設(shè)計(jì)。仿真計(jì)算表明，傳統(tǒng)極限拉偏方法存在極大的保守性，而采用可靠性描述的系統(tǒng)指標(biāo)要求更符合系統(tǒng)的實(shí)際分布情況。通過概率優(yōu)化得到的設(shè)計(jì)參數(shù)使系統(tǒng)在統(tǒng)計(jì)上滿足指標(biāo)要求，降低了系統(tǒng)的保守性，優(yōu)化了設(shè)計(jì)方案。

在未來的工作中，可以引入箭體的彈性振動(dòng)，并應(yīng)用概率優(yōu)化方法調(diào)節(jié)校正網(wǎng)絡(luò)參數(shù)，以及基于現(xiàn)代先進(jìn)控制理論進(jìn)行控制器設(shè)計(jì)，進(jìn)一步提高概率設(shè)計(jì)方法的實(shí)用性。