基于可靠性优化的固体火箭姿态控制器设计

严 恺 吴燕生 张 兵 钟 震

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方法，针对幅值裕度和相位裕度指标，对控制器参数进行了概率优化设计。仿真计算表明，传统极限拉偏方法存在极大的保守性，而采用可靠性描述的系统指标要求更符合系统的实际分布情况。通过概率优化得到的设计参数使系统在统计上满足指标要求，降低了系统的保守性，优化了设计方案。

在未来的工作中，可以引入箭体的弹性振动，并应用概率优化方法调节校正网络参数，以及基于现代先进控制理论进行控制器设计，进一步提高概率设计方法的实用性。