高超声速飞行器鲁棒控制系统的设计
Christopher I. Marrison and Robert F. Stengel
Princeton University, Princeton, New Jersey 08544
本文设计了高超声速飞行器纵向平面鲁棒控制系统。飞行器纵向平面的非线性数学模型包含了28个不确定参数。利用遗传算法搜索每个控制器的系数设计空间；利用蒙特卡洛算法检验每个搜索点处的稳定性和鲁棒性。补偿器的鲁棒性用概率函数来表示，该函数表示在参数可能变动范围内，闭环系统的稳定性等性能指标落入允许范围的概率。设计了一性能指标函数，使其最小，从而产生可能控制器系数空间。这种设计方法综合考虑了不同的设计目标，辨识了鲁棒性指标下的系数的不确定性。这种方法有效利用了计算工具，广泛考虑了工程知识，设计出了能够应用于实际的控制系统。
本文中用到的符号：
a ——声速，ft/s
D C ——阻力系数
L C ——升力系数
()M C q ——俯仰角速率引起的俯仰力矩系数
()M C α——攻角引起的俯仰力矩系数
()M C E δ——舵偏引起的俯仰力矩系数
T C ——发动机推力系数
c ——参考长度，80ft
D ——阻力，lbf
h ——高度，ft
yy I ——俯仰转动惯量，6710?slug-ft 2
L ——升力，lbf
M ——马赫数
yy M ——绕俯仰轴的转动力矩，lbf-ft
m ——质量，9375slugs
q ——俯仰速率，rad/s
E R ——地球半径，20 903 500 ft
r ——距地心距离，ft
S ——参考面积，3603ft 2
T ——推力，lbf
V ——速度，ft/s
α——攻角，rad
β——喷管开度
γ——弹道倾角，rad
E δ——舵偏角，rad
μ——地球引力常数，161.3910?ft 3/s 2
ρ——空气密度，slugs/ft 3
引言
高超音速飞行给控制系统的设计提出了新的挑战。极高速度使得飞行器对飞行条件的变化非常敏感，例如，高度110 000ft ，速度15Mach （15 060 ft/s ），攻角增加1°就会产生11.5ft/s 2的法向加速度，即大概3g 的过载。对于传统的飞行器配置方案，长周期振动有可能轻微发散或者不稳定，短周期振动可能不稳定，推力幅值也易受到攻角变化的影响。大气特性以及气动参数测量上的困难加剧了这个问题。但是，高超声速飞行器的飞行控制系统必须保证飞行器总是处于稳定状态，对于飞行员及自动驾驶仪的指令有满意的响应特性，外界的干扰对飞行器不会产生难以接受的附加运动。本论文设计了满足上述要求的控制系统，重点考虑了参数不确定性条件下的鲁棒性。
随机鲁棒性分析用概率P 量化了补偿器?欠鲁棒性的程度，补偿器?模型参数的变化将会引起闭环系统不可接受的响应特性。概率函数P 的一种简单的设计方案是将某个指标函数在不确定参数变动空间上积分，即
[(),]pr()P I H d νν?νν=? (1) H 表示飞行器的数学模型，ν是不确定参数的可能变动空间；pr()ν是ν的概率密度函数。[]I 是双指标函数，如果()H ν和?构成一个可以接受的系统其值为零，反之为1。Fig.1描述了空间划分为可接受和不可接受区域。
性能不可接受的标准由设计者制定，例如，可以指定为不确定性，不满足响应包线，超过了控制器的阈值限制。这种概率设计原则可以应用于线性或非线性、时变或非时变、连续或者离散系统的设计。
在多数情况下，方程(1)并不能直接进行数值积分。实际中常用蒙特卡洛方法1替代pr()ν，它利用了不确定参数ν的随机采样值。每次试验的个体选择记为k ν。对于每个k ν，
检验闭环系统的响应特性，以检验其是否满足可接受条件。例如，线性系统的稳定性指标函数其值为1，如果其闭环系统的特征值存在实部为正的情况；反之其值为0。性能指标函数，例如阶跃响应的超调量或者稳定时间超过设计允许范围，其值为1；反之为0。从空间ν中N 次采样，重复进行蒙特卡洛试验。N 值取决于所需要的精度。方程(1)可以用下式进行估计： 11[(),]N k P I H N ν?==∑
(2)
当N →∞时上式逼近精确值。性能特性和稳定性的概率估计综合在一起考虑不仅可以用来描述控制系统的鲁棒性，而且可以用来设计鲁棒控制器2-10。
对于特定应用的补偿器的设计是一个主观过程，需要对很多指标进行权衡考虑，有时这些性能指标根本不能比较。这种权衡可以考虑将各个概率赋以主观权重，用标量性能函数J 来表示。如果补偿器?中的参数能够用常值向量d 来表示，那么P 和J 是d 的函数，鲁棒控制系统的设计是要使下面的性能指标函数最小：
21??()[()]M m m m J w p ==∑d d
(3)
其中，M 是稳定性和性能指标的数量，?()m p
d 是第m 个准则函数，m w 是第m 个准则函数的权重。参考文献[10]针对线性控制问题利用线搜索（lin
e search ）设计了鲁棒补偿器，文献
[11]利用遗传算法来估计?()J
d 的最小值。上述方法得到的补偿器鲁棒性非常强，控制失败的概率远低于按其它方法设计的控制器8。文献[10]中的线搜索是一种有效但费时的方法；遗传算法相比下减小了一个数量级的计算量11。
随机鲁棒性分析和设计方法对传统设计方法进行了补充，它的稳定性和性能指标函数都是基于传统设计准则。如果关注的是频率特性，很容易将频域性能指标加入到指标函数中去。但是需要注意的是，幅值增益和相位裕度并不是鲁棒性的可靠指标9。理由很简单：实际参数变动时，并不只是引起幅值增益和相位裕度的变化，也会影响那奎斯特图的形状。最大奇异值会引起增益裕度的多元化，并影响鲁棒性，原因也是如此7。
本文利用蒙特卡洛方法和遗传算法12-15解决复杂问题：在高度110 000ft ，马赫数15的飞行条件下进行高超音速巡航飞行器的飞行控制系统的设计。该飞行器的数学模型为纵向平面内的模型，包含了重要的非线性影响因素。我们考虑了39个稳定性和性能鲁棒指标准则，因此较其它鲁棒性设计方法的适应性更强。该设计方法利用一种巡航飞行条件进行了检验，并容易扩展到其它飞行条件，利用增益调度方法将结果扩展到全部飞行包线。
高超音速飞行器的数学模型
纵向平面的数学模型包含了重力模型的平方反比定律，以及由弯曲的飞行轨道引起的向心加速度。在高度110 000ft ，马赫数15飞行条件下，非旋转地球的向心加速度大小为10.8ft/s 2。速度，飞行倾角，高度，攻角，和俯仰角速度的数学方程分别如下： 2cos sin T D V m r αμγ-=- (4) 22sin ()cos L T V r mV Vr αμγγ+-=- (5)
sin h V γ= (6)
q αγ=- (7) /yy yy q M I =
(8) 升力、阻力、推力、俯仰力矩、距地心距离分别为：
212L L V SC ρ= (9)
212D D V SC ρ= (10)
212T T V SC ρ= (11)
21[()()()]2yy M M M M V Sc C C E C q ραδ=++ (12) E r h R =+
(13) 空气动力系数和大气参数是飞行器状态和控制量的函数，在仿真中通过查表插值或者曲线拟合得到。此处利用相对简单的函数来表示参考巡航点处的气动参数。大气密度和声速模型来源于参考文献[16]，气动参数借鉴NASA Langley 高超声速飞行器仿真模型17。
此处设计中，28个与惯量、推力、气动等有关的参数假设为不确定的，v 表示不确定向量ν中的一个元素。这些参数为：
10m v m = (14)
20yy I v I = (15)
30S v S = (16) 40c v c =
(17) 80.00238exp()24,000h v ρ-= (18)
924567(8.99109.1610996)a v v h v h --=?-?+ (19)
/M V a = (20)
109 1.91(0.493)L v C v M α=+ (21)
2211121314150.0082(171 1.152)(0.00120.0541)D C v v v v M v M αα=++?-+ (22) 218161701738[1164()](1)v k v v M
αα=--+ 1919(10.15),1(10.15),1
T k v if C k v if ββββ+? (23)
21342202223()10[0.06][21201]v M M C v e v v ααα--=--+- (24)
224252627()(/2)(0.025 1.37)(0.00210.00530.23)M C q c V qv v M v v αα=-+?-+- (25) 28()0.00051()M C E v E δδα=- (26)
方程(4)~(26)定义了飞行器的数学模型()H ν。此处每个v 假定服从正态分布，均值为1，标准偏差为0.1，这样就定义了pr()ν。常用的变量如机翼面积和平均气动弦长此处假定未知。在实际中，应考虑变量的不确定性，但不确定性很小时，可作为确定值处理。将变量的不确定性假设为正态分布，这样能够保证不稳定性的概率永远不为零，这是由于变量两侧的分布趋向于正或负无穷。将每个v 取值为1，得到标称系统。在蒙特卡洛试验中，每个v 值由随机数发生器产生，这种方法并不仅限用于正态分布情况。不确定性可以用恰当的概率分布来表示，例如均匀分布，多峰分布，离散分布。
将飞行器纵向平面数学模型线性化，得到巡航飞行条件下特征点处（M =15，h =110 000ft ，q =0deg/s ，0γ=deg ）的开环特征值-0.895，0.784，-0.00021±0.0362j ，0.00011。特征点（平衡点）利用非线性模型计算得到。前两个特征值表示不稳定的短周期运动模态，复数特
征值表示轻微震荡的长周期运动模态，最后一个特征值表示高度模态存在轻微不稳定18。这
样，必须设计稳定的反馈控制系统限制巡航条件下姿态和高度的发散。
巡航特征点处的非线性模型的响应特性由攻角、推力、舵偏等确定，这些变量提供力和力矩，保证飞行器的稳定、水平飞行。确定特征点后，能够得到非线性模型的阶跃响应。利用鲁棒性能准则评估变量的响应特性。
稳定性和性能准则
此处控制系统的设计中主要关注三个方面的鲁棒性：稳定性，速度响应特性和高度响应特性。它们用39个指标函数来表示。特征点及其邻域范围内的稳定性利用相应的线性化模型来判断（准则1），稳定性由闭环系统的特征值决定19。有19条准则（准则2~20）根据施加在非线性模型上的速度阶跃响应（速度变化100ft/s ）决定，见表1。有两条并行准则的（例如,25V Ts I 和,50V Ts I ）则是使遗传算法寻优时区分补偿器，某些情况下补偿器会严重违背高要求准则，而较少违背低要求准则。
剩下的19条准则（21~39）建立在高度的阶跃响应基础上（高度变化2000ft ）。高度准则类似于速度准则（将表1中的下标V 用h 来代替），不同之处在于时间延长一倍。
设计中采用如下的性能指标函数：
3921??()[()]m m m J d w p d ==∑
(27)
其中，稳定性准则的权重取为10，基本的性能准则权重为1高要求准则权重为0.1。
对于一般函数的全局搜索方法，在有限的试验中不能保证得到全局最小值。但有一些能够估计全局最小值的方法，方便设计者进行判断。很早人们就知道，随机搜索在利用很多的采样数据情况下，搜索的最终结果（final J ） 服从Weibull 分布20。Weibull 分布函数为 {}Pr()1exp [()/]c final J x x a b ≤=--- (28) 其中0x a ≥≥，0b >，0c >。a 是位置变量，b 是规模变量，c 是形状变量。a 是final J 的可能最小值。在寻找最小值的过程中，*()a J d =。如果搜索算法运行n 次，得到n 个final J 值，从而利用这些值得到final J 的经验分布。a ，b ，c 的取值根据Weibull 分布和经验分布的偏差利用最小二乘算法得到。
文献[21]中，随机启发式搜索算法的结果同单纯随机搜索算法一样，也服从Weibull 分布。Golden 21利用Weibull 分布来估计*()J d ，解决了旅行商问题。文献[21]得到了置信区间 *,,Pr[()]1n final n final n J b J d J e --≤≤=-
(29) 其中,final n J 为n 次搜索中的最优值。Cooke 22发展了另一种估计*()J d 的方法，并不需要事
先假设特定的概率分布。Cooke 考虑了i Y 的n 次随机采样值，以及概率分布函数Pr()Y x <。按照由小到大排序（1i i Y Y +<），最大值为n Y ，它的概率分布函数定义为Pr()n Y x <，其值
等于[Pr()]n Y x <。n 次采样最大值的期望值为
max
()pr ()n min Y n Y x Y E Y x x dx ==? (30)
其中，pr ()n Y x 是概率密度函数，即Pr()n Y x <。通过部分积分，Cooke 得到了估计器 1max 0(1)?2(1)[1]n n n n n i i i i Y Y Y n n --=+??=----????∑ (31)
利用关系式Y J =-，得到*
()J d 的估计器 1
*0(1)()2(1)[1]n n n n n i i i i J d J J n n --=+??=----????∑ (32)
其中1i i J J +<。这种方法计算简单，同时能够得到a 的估计值。
补偿器结构
控制量由三部分信号构成：根据特征点处得到的标称控制量0u ；从特征点改变到另一
种飞行状态的控制增量*u ?；保证稳定性和响应特性的动态控制量d u ?。
将(4)~(8)式在特征点0x 处进行一阶泰勒展开，得到*u ?、d u ?的表达式。
(,)x f x u = (33)
[,,,,]T x V h q γα= (34)
[,]T u E βδ= (35)
0x x x =+ (36)
*00d u u u u u u =+=++ (37) 00(,)x f x u F x G u F x G u ≈++=+
(38) 雅可比矩阵F f x =??，G f u =??的值在特征点处取得。
如果*y 是指令向量，忽略干扰，得到平衡状态*x 和控制量*u 1***0x u F G x H H y u -??????=???????????? (39)
其中***x u y H x H u =+。此处得到的控制律为连续时间形式，线性二次型最优调节器
（LQR ）能够起到比例-积分（PI ）补偿的作用2。积分补偿能够消除稳态误差，这种误差由实际系统与此处用到的系统(39)之间的偏差产生。动态响应特性通过前馈补偿得到改善，需要用到*x 和*u 。对于非线性系统，将控制误差限制得很小，得到状态和控制的参考值非常重要。这点在分析、设计及试验中得到了广泛地验证（见文献[2]中的6.2）。
对于高超声速飞行器，需要控制的量为速度和高度，因此将状态变量扩展为
[,,,,,,]T PI v h x V h q γαζζ= (40)
*()v V V dt ζ=-? (41)
*()h h h dt ζ=-? (42) 最优控制律为
PI PI d u C x =- (43) 其中，PI C 为增益矩阵。PI C 由最优算法中性能指标函数中的权重决定，
2PI PI 0PI PI T T T T g d d L x Q x k V x F F x u R u γγ=++
(44) Q 、R 分别为状态和控制的加权矩阵，为简化起见，均定义为对角矩阵。20PI PI
T T g k V x F F x γγ项决定法向加速度的比重，近似等于0V γ。定义F γ为矩阵F 中与γ有关的行向量，法向加
速度等于0PI V F x γ。g k 为标量。
将控制增益乘以标量常值C k 有可能会改进鲁棒性，尽管一定程度上降低了最优性23。最终的反馈控制律为
PI PI d C u k C x =- (45)
有10个设计参数需要用遗传算法进行寻优。g k 、C k 、R 中的对角元素、Q 中的六个对角元素（速度权重Q V 已确定）
[,,,,,,,,,]h q v h g E C d Q Q Q Q Q Q k R R k γαζζβδ= (46)
利用遗传算法对Q 、R 进行寻优，从而将设计者从繁重的计算工作中解放出来。在[稳定性和性能准则]一节中看到，很方便寻找概率性能指标函数中的39个权重（方程(3)）。极端重要的准则赋以高权重，例如给稳定性准则赋以权重10；重要的准则赋以中间权重，例如1；不重要的赋以低权重，例如0.1。在此处，Q 、R 为对角矩阵，只寻优对角元素，非对角元素和状态-控制耦合元素可以通过扩展搜索算法或者利用本节末处提到的隐式模型跟随结构（implicit-model-following ）。矩阵的正则性则由搜索范围d 决定，不允许搜索非正则数值。如果Q 、R 取值不理想，则将其丢掉。由遗传算法确定的最优补偿器具有满意的Q 、R 。
d 的初始值由x 和u 的期望偏差决定，将Q 、R 中的对应元素取为期望值倒数的平方24。例如，高度h 的期望偏差为2000ft ，则将h Q 设为2
(12000)。系统的响应由这些权重决
定。不断调整权重，得到仅违背少数准则的响应特性。这些值为搜索提供了初始值。
此处采用LQ-PI 控制器结构来验证概率鲁棒控制方法，也可以利用其它控制结构。基于比例-积分-滤波方法的随机设计方法已经用于高超声速巡航问题25，隐式模型跟随算法已经采用概率方法研究3,4，文献[26]发展了一种确定单输入单输出传递函数系数的概率搜索算法。
仿真结果
在特征点处（M =15，h =110 000ft ，q =0deg/s ，0γ=deg ）对开环和闭环系统进行了随机鲁棒性分析。利用(14)~(25)式给出的参数不确定性，开环系统不稳定的概率为0.816，违背其它准则的概率在0.816~1之间。
利用基本的补偿器，开环系统的不稳定概率得到降低（0.014），违背其它准则的概率也
得到降低。违背了稳定时间准则，对于速度响应，,25P 0.966V Ts =；对于高度响应，,50P 1h Ts =。
利用基本的补偿器，遗传算法基于(27)式对基本补偿器的参数进行寻优。由于稳定性的鲁棒性非常重要，在P i 中的权重是最大的（权重为10）。给低要求的性能准则赋以权重1，例如要求稳定时间大于50s 。给中等要求的性能准则赋以权重0.1，例如妖气稳定时间大于25s 。
优化开始时，先进行随机搜索。然后优化群体，直至满足结束条件。例如，结束条件为
达到20000次蒙特卡洛试验。利用(32)式，运行8次来估计全局最小值，?J
的8次结果分别为2.09，1.86，2.16，2.71，1.99，1.83，1.72，1.74。这样，第七次给出了最优结果。
全局最小值的估计方法可以由下面两种方法获得：利用得到的结果确定Weibull 分布（式
(28)），或者利用Cooke 估计器（式(32)）。此处两种方法全部采用。利用前面得到?J
的8次结果，(32)式得到全局最小值为1.703。为了证实这一结果，Weibull 分布中的a 取为此值。这8个观测数值用来估计遗传算法搜索结果的概率分布，b 、c 的取值应使(27)式符合观测结果，这样得到a =1.703，b =0.28，c =0.8。8次结果及相应的Weibull 分布见图2所示。可见，全局最小值的估计值取为1.703是可行的。
如果假设b 的取值也是准确的，那么(29)式可以写为
*8Pr[1.720.28() 1.72]10.9997J d e --≤≤=-= (47) 可见全局最小值有99.97%的概率落在1.44~1.72之间。进一步，接受该全局最小值，采用给定的补偿器结构，J 的最优值为1.703，第七次运行时补偿器得到的值为1.72，误差仅为1%。 这种分析和计算方法有很大的计算量，尽管已经充分利用了计算机。在Silicon Graphics Indy 工作站上运行matlab ，采用本文的方法设计补偿器所需要的时间为24h 。转换成C 或者Fortran 语言，并利用最新的并行计算机，所需时间以分或者秒计。
一般情况下鲁棒控制问题所需的计算量很大。精确求解所需要的运行次数与变量数量是指数级的关系，而且对于一个简单的、连续时间的不确定参数，变量的数量也是无穷的。所有的鲁棒设计问题需要在一定程度上逼近所求问题。此处利用蒙特卡洛方法评估和数值搜索。利用蒙特卡洛方法求解特征值和响应特性所需的运算次数与状态变量、控制量是多项式的关系，与响应时间、不确定变量的数量是线性关系。遗传算法所要的运行时间与设计参数是线性关系。
有很多方法将所得到的控制系统扩展到整个飞行空域。此处利用LQ 设计方法，最合适的方法是增益调度。例如，本文的第二作者和他的同事在上世纪70年代设计了集中了增益调度、LQG 、数字方法的直升机控制器，飞行试验证实了良好的实时操作性能（盘旋，后飞，爬高，转弯，水平巡航）27,28。
对最优补偿器进行分析
图3比较了基本补偿器与优化后的补偿器的鲁棒性。不稳定的概率由0.014减小至0.001，违背稳定时间准则的概率改善了4.5倍（,50P V Ts ）和1.02倍（,100P h Ts ），,50P h Ts 仍为
1。除了三条，其它准则得到改善：响应速度指令耗费超过50%燃油的概率（,50P V T δ），响应高度指令过载大于2或4的概率（,2P h g ，,4P h g ）。
存在0.1%的不稳定性也是飞行器的隐患。降低这个概率可以通过增大(27)式中的P i 的权重，允许LQ 性能指标函数中存在权重的交叉乘积项，可以改变控制器的结构。如果变量服从正态分布，P i 不能取为0。这不是补偿器的特点，而是pr()ν的内在属性。为了描述这一点，利用有界分布重新对补偿器进行分析：正态分布限制在±0.2（之前的分析中标准差为0.1）。这种情况下，不稳定的概率变为零。对于具体的应用了解变量的变动范围非常重要。
利用标称参数，带有基本补偿器的闭环系统的特征值为-1.1956±0.9184i ，-0.0210±