基于可靠性優化的固體火箭姿態控制器設計

嚴 愷 吳燕生 張 兵 鐘 震

1. 北京宇航系統工程研究所，北京 100076；2. 中國航天科技集團公司，北京 100048

0 引言

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

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

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

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

1 固體火箭姿態控制數學描述

1.1 火箭的動力學模型

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

(1)

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

(2)

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

(3)

1.2 性能指標

目前采用可靠性優化方法進行控制器設計的文獻大都采用了LMI來描述系統的穩定性要求，主要是由于LMI具有對設計參數為凸的良好性質，方便優化計算[14]。然而，LMI描述的假設條件過強，且只能描述系統的穩定性條件，因此應用性較差。在實際工程中，控制系統分析與設計主要在頻域進行，除穩定性外，更關心的是系統的穩定裕度。主要包括幅值裕度(Gain Margin, GM)和相位裕度(Phase Margin, PM)，用以表征系統的相對穩定性和動態性能指標。為此，本文將采用GM和PM作為系統的性能指標，計算式為

(4)

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

2 概率分析與設計方法

2.1 可靠性分析

相比確定性設計，概率設計的特點是將系統模型中的不確定性參數描述為滿足一定概率分布的隨機變量，系統的性能指標要求也由確定性要求放寬為概率要求，如式(5)所示

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

(5)

式中x,k分別表示系統的不確定參數和確定性設計變量；g(x,k)表示系統性能的指標函數，稱為性能函數(Performance Function)；R為系統各指標滿足相應要求的概率，即可靠性；R1表示對系統指標滿足的概率要求，當R1=1時，右式與左式等價，概率約束退化為確定性約束。根據概率理論，式(5)中的概率計算表達式為

(6)

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

(7)

此時概率Pi可表示為

(8)

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

(9)

當用此u的線性函數近似指標函數時，由概率論知識有

(10)

由于此時g′(u)服從正態分布，因此有

(11)

由于u*為極限狀態方程上距離原點最近的點，則可得

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

(12)

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

(13)

其中βt表示設計要求的可靠性對應的可靠性指數，當MPP點求得的β*>βt時即認為可靠性滿足要求，如圖1所示。

圖1 RIA法與PMA法示意圖

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

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

(14)

PMA法同樣為一個優化過程，但其約束形式簡單，目標函數復雜，因此收斂性更好，本文將采用PMA法進行系統的可靠性分析，以作為概率優化的基礎。

2.2 可靠性優化設計方法

可靠性優化問題本質上為一個以概率要求為約束的優化問題，可以描述為以下形式

搜索:k
min:J(k)

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

(15)

其中J(k)為優化目標函數，k為設計參數。現有的RBCD文獻大都采用了雙循環優化結構，即將可靠性分析過程直接嵌入到外層參數尋優過程中，如圖2所示。由2.1節可知，基于MPP點的可靠性分析方法本質上也是一個優化過程，因此雙循環法會使兩個優化過程的耦合嵌套，從而導致計算量劇增，降低了該方法的實用性。

圖2 雙循環結構RBCD流程

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

搜索:k

min:J(k)

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

(16)

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

圖3 概率約束的等價描述示意圖

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

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

圖4 SORA法設計流程圖

3 基于可靠性優化的控制器參數設計

根據式(1)所示的火箭動力學模型，得到火箭俯仰通道的開環傳遞函數為

(17)

控制律為

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

(18)

控制系統框圖表示為

圖5 控制系統框圖

系統的開環傳遞函數為

L(s)=(kp+kds)Go

(19)

則設計變量矢量和隨機變量矢量為

(20)

系統指標選取為GM和PM，構造性能函數

(21)

其中PMd,GMd分別為要求的幅值裕度和相位裕度。在相同的姿態角偏差下，控制系統的靜態增益kp與指令發動機擺角的輸出量直接相關，表征了對執行機構的能力需求。kp越大，在同等姿態角偏差下對擺管的擺動能力需求也就越高，且過大kp值對箭體彈性振動的穩定不利。因此，在滿足控制指標的前提下，較小的kp值可以降低控制系統的能力需求。為此優化目標函數可以選取為靜態增益的幅值|kp|。此時控制器參數的概率優化模型為

搜索: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 仿真校驗

本文研究對象為固體火箭的一級飛行段，選取最大動壓點這一典型的特征秒，結合參數的散布情況，所研究模型的剛體系數標稱值和方差如表1所示。

表1 動力學模型參數標稱值

系統的性能指標要求為

(23)

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

表2 極限狀態設計的系統指標

由計算結果可以看出，相位裕度遠大于期望值，因此幅值裕度為主要考慮的指標。按極限拉偏設計的參數保證了系統性能指標的下限滿足給定的指標要求，因此保證了理論上的可靠性為1。同時，可以發現幅值裕度小于9.14dB情況僅占所有情況的5%。此時系統的開環Bode圖及其散布情況如圖6所示。

圖6 極限設計對應Bode圖

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

作為比較，采用概率要求進行設計。設要求的系統可靠性為R1=0.95，即系統在隨機不確定下滿足指標要求概率至少為95%。按此標準進行優化設計得到的控制器參數為k1=[1.58,0.5]，此時系統的各指標如表3所示。

表3 95%可靠性要求設計的系統指標

由計算結果可以發現，雖然此時系統的下極限對應的幅值裕度為6.54dB且小于給定的指標要求，但概率優化設計的結果保證了幅值裕度小于7.83dB的情況占比不超過5%，整體的可靠性為0.959，系統的性能滿足提前給定的概率要求。此時系統Bode圖的散布情況如圖7所示。

圖7 95%可靠性設計對應Bode圖

由圖可知，隨機分布的參數對應的Bode圖的整體散布仍然主要位于5%和95%分位線以內，概率設計能夠滿足預期要求。風攻角為αwp=5°,αwq=0.8°時，計算兩種設計下的擺角，得到

(24)

由上述仿真計算發現，當把系統的可靠性要求降低5%時，兩種方法的控制系統靜態增益kp減小了20%，發動機擺角減小了15%，顯著降低了控制系統的性能需求。

5 結論

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

在未來的工作中，可以引入箭體的彈性振動，并應用概率優化方法調節校正網絡參數，以及基于現代先進控制理論進行控制器設計，進一步提高概率設計方法的實用性。