沈海东，佘智勇，曹 瑞，刘燕斌，陆宇平

(1. 南京航空航天大学航天学院，南京 211106；2. 北京空天技术研究所，北京 100074；3. 南京航空航天大学自动化学院，南京 211106)

0 引 言

凭借着发射成本低、任务可靠性高、部署方便灵活的优势，吸气式水平起降可重复使用空天飞行器成为未来空天运载系统发展的主要方向，在军事和民用领域都具有广阔的应用前景[1]。因此，其关键技术的研究吸引了国内外学者的广泛关注[2-3]。

有关水平起降空天飞行器的研究始于20世纪80年代，目前已经公开的方案均处于任务可行性论证阶段，典型代表包括：美国的国家空天飞行器计划[4](National aerospace plane, NASP)，英国的“霍托尔”计划[5](Horizontal takeoff and landing, HOTOL)、“云霄塔”计划[6](Skylon)，德国的“桑格尔”计划[7](Sänger II)等。受当时技术水平、研制难度、经费及政策的影响，这些计划先后被搁置，但世界各主要大国在关键技术攻关上的努力并未终止[8]。

进入21世纪后，世界各主要大国陆续在高超声速技术方面取得突破。总体设计方面，由于空天飞行器各系统间存在严重的耦合效应，传统“递进式”的飞行器设计方法已无法满足要求，因此多学科优化的设计方法引起了广泛重视[9]。推进系统是决定空天飞行器能否实现大空域、宽速域飞行的关键因素。单一模态的发动机显然无法满足该任务需求，因此国内外对包括涡轮基组合循环发动机[10]、火箭基组合循环发动机[11]、协同吸气式火箭发动机[12]等在内的多种组合循环推进系统开展了研究。

近年来，随着美国洛克希德·马丁公司SR-72[13]和波音公司Valkyrie II[14]概念方案的提出，有关水平起降空天飞行器的研究再次升温。为满足空天飞行器水平起降、高超声速巡航的任务要求，这两种方案高度类似，均采用了大后掠三角翼无尾式气动布局和并联式涡轮基组合循环发动机。

无尾布局飞行器一般采用升降副翼同时控制俯仰、滚转通道，这种特殊的配置方式增加了舵面设计的复杂度。首先，由于缺乏水平尾翼，升降副翼既需要完成飞行器的配平，又要实现纵向、横侧向的动态操控，即升降副翼的设计需要兼顾飞行器横侧向静态、动态性能。其次，相较于传统升降舵，机翼后缘升降副翼的纵向力臂更短，导致其低速阶段纵向控制能力较弱。

为增强飞行器纵向控制效率，最直接的方法就是增加控制舵面的面积。但如果舵面面积过大，一方面舵面偏转时引起的气动性能参数变化量过大，另一方面维持舵面偏转角所需的铰链力矩也会大幅增加[15]。另一种解决方法是通过向后移动质心位置，放宽飞行器纵向静稳定度。但这种方式会降低飞行器的静稳定性，实际飞行中需要引入增稳系统以改善飞行品质。随着飞行器静不稳定度的增强，增稳系统抵抗外界扰动时所需的舵面偏转角速率增大，这对作动器造成了严峻的挑战。此外，为了确保放宽静稳定性后的飞行器仍然可控，需要设计增稳控制器。但控制器的设计需要依赖具体的对象，这就造成了概念设计阶段控制器设计与舵面设计间的相互耦合。

这种问题在控制领域被称为本体—控制一体化优化，可以通过顺序、迭代、嵌套、同步等不同策略求解，其中，顺序和迭代策略无法保证所得解的全局最优性[16]。文献[17]通过将控制模块融入到了多学科优化(Multidisciplinary optimization, MDO)过程中，运用协同优化算法实现了民航客机的多学科并行分析和设计。仿真结果表明，放宽静稳定性后飞行器升降舵面积显著减小。文献[18]提出了一种将控制参数和模型参数同步优化的算法，基于H∞控制器，将飞行品质、机动性等性能指标表示成线性矩阵不等式(Linear matrix inequality, LMI)的形式进行数值求解。该方法在文献[19]中得到进一步发展，相关成果被应用到了高超声速飞行器控制一体化设计中[20]。

然而，基于传统H∞理论设计得到的是没有明确结构的全阶控制器，很难应用到工程中。考虑控制器结构约束时，原来的线性矩阵不等式约束将变为双线性矩阵不等式(Bilinear matrix inequality, BMI)，对于该问题的求解属于非确定多项式难题(Non-deterministic polynomial hard, NP-hard)，传统的基于梯度的优化算法不能直接应用于该问题。文献[21]首次提出了基于非光滑优化的结构化H∞设计方法，通过Clarke次微分理论计算性能指标的局部最优值。为保证所得到的局部最优值更接近于全局最优，需要选择不同的初始值进行多次优化。该方法具有较强的通用性，在控制器参数调节、本体—控制一体化设计等领域得到了广泛应用[22-23]。

目前，对空天飞行器的研究主要集中在高速段，对低速阶段的研究较少。文献[24]提出了一种面向控制的空天飞行器低速动力学建模方法。通过参数化建模方法实现了飞行器几何构型的快速调整，基于势流理论和0维混合排气涡扇发动机原理，分析了低速阶段的飞行器气动/推进特性。在此基础上，本文构建了包含总体参数的气动特性代理模型，提出了一种基于非光滑优化算法的空天飞行器控制舵面—控制参数一体化设计方法。

1 空天飞行器面向控制舵面优化的参数化建模

1.1 空天飞行器动力学模型

如图1所示，为实现大空域、宽速域飞行，目标飞行器采用大后掠无尾三角翼布局，依靠涡轮基组合循环发动机提供动力，其具体特征参数可参考文献[24]。

图1 目标飞行器几何构型图[24]Fig.1 Object vehicle configuration[24]

空天飞行器机体/发动机一体化的设计使得发动机对飞行器的姿态十分敏感，为保证发动机的稳定工作，实际飞行中应当避免出现较大的横侧向机动。因此，本文主要研究空天飞行器的纵向运动。

忽略地球自转的影响，根据拉格朗日方程[25]可推导得到飞行器纵向刚体动力学方程：

(1)

式中:m为飞行器的质量，Iy为绕机体轴y轴的转动惯量，状态量x=[V,α,h,q,θ]T分别为飞行速度、迎角、高度、俯仰角速度和俯仰角；L,D,T,M分别为飞行器所受升力、阻力、推力和俯仰力矩，它们的表达式如式(12)所示。

(2)

式中:c为参考弦长，ρ为大气密度，s为飞行器参考面积，CL,CD,CT,CmC分别表示升力系数、阻力系数、推力系数及俯仰力矩系数。

飞行器气动系数由两部分组成：洁净体气动系数Ci(clean)及舵面偏转引起的气动系数增量ΔCi，即：

Ci=Ci(clean)+ΔCi

(3)

1.2 升降副翼定义及气动系数增量求解

定义升降副翼的面积缩放系数η为：

(4)

式中:Sini为升降副翼初始面积。

对于全展长升降副翼，η可表示为：

(5)

式中:cini为初始舵面弦长(25%机翼弦长)，cele为调整后的弦长，不同η值对应的升降副翼如图2所示。

图2 不同缩放系数η下的升降副翼(阴影部分)Fig.2 Elevon configuration for different η

如图3所示，对于不同η值，可以根据几何外形参数化方法获得对应的升降副翼构型，结合面元法即可获得不同飞行条件下升降副翼偏转引起的气动系数增量。为增强所得数据的可靠性，可进一步结合少量CFD数据对获得的数据进行校正。

图3 不同尺寸升降副翼气动参数计算Fig.3 Aerodynamic modeling for different elevon sizes

1.3 可变构型空天飞行器气动特性代理建模

为提高后续优化效率，需要对空天飞行器气动特性进行代理建模，即构建一个高效的气动特性解析模型。

传统的多项式代理模型，大多基于预先已知的模型结构，将模型辨识问题转换为最小二乘求解。对于构型可变的空天飞行器，传统的气动系数模型结构无法反映出结构参数的影响。因此本文提出了一种基于鸽群算法[26]的可变结构飞行器气动特性代理模型自动获取方法。

该算法的流程如图4所示，具体步骤如下：

图4 基于鸽群算法的气动特性代理建模Fig.4 Flowchart of aerodynamic surrogate modeling based on PIO algorithm

1) 读取自变量x1∈Ro，因变量y1；

3) 构建当前阶次对应的完整代理模型形式：

4) 根据最小二乘法，求取代理模型多项式各项对应的系数a1,a2, …,ar0；

7) 基于高斯变异鸽群算法在所构建代理模型优化初始集合中进行搜索；

9) 判断搜索次数k

10) 输出最高次数为n，最大项数为r的最简代理模型：

12) 输出不同尺寸升降副翼对应的多项式系数，并进一步将其写作升降副翼缩放系数η的表达式，即可获得包含总体参数的气动特性代理模型：

(6)

(7)

2 控制一体化设计

本节介绍了一种综合考虑飞行品质、鲁棒性、机动性等多种约束下的空天飞行器升降副翼缩放系数η与控制器参数K一体化优化设计方法。其主要思想为在给定的控制结构下，将飞行器控制参数调整、升降副翼尺寸设计转化为多约束下的多目标优化问题，然后运用非光滑优化算法进行求解，以同时获得满足约束条件的控制器参数和模型对象参数。

2.1 模型参考跟踪控制

本文采用的纵向模型参考跟踪控制结构如图5所示，主要包括C*控制器、飞行器短周期模型、参考模型和执行器模型等四个部分。图中Nzc表示过载指令，nzref为期望指令跟踪信号，nz为实际过载跟踪信号，z1表示跟踪误差，Δδe为升降副翼偏转角，uact为执行器输出信号。

图5 基于C*结构的控制舵面—控制参数一体化设计Fig.5 Integrated design and optimization based on C* structure

1)C*控制器

为实现法向过载的跟踪，使用了工程中经典的C*控制器[24]，即通过俯仰角速率和过载反馈信号来控制俯仰通道稳定。

基于该结构进行飞行器纵向控制时，共有四个控制参数需要调节，分别为过载反馈增益Knz，俯仰角速度反馈增益Kq，指令信号前馈增益Knzc，积分器增益Ki，即K=[Knz,Kq,Knzc,Ki]。则C*控制器总输出为升降副翼偏转角的增量：

(8)

2)短周期模型

将式(1)在平衡点x0处进行泰勒展开，仅保留线性项，可以得到飞行器短周期小扰动线性方程：

(9)

(10)

式(9)、(10)中气动力/力矩量纲导数的具体定义见文献[25]。

3)参考模型

期望的闭环系统响应通过二阶系统Gref表示：

(11)

式中:ξref和ωref分别为期望闭环响应模型的阻尼比和自然频率，本文取ξref=0.7，ωref=1 rad·s-1。

4)执行器模型

控制信号Δδe通过执行器后输入到飞行器动力学模型，本文将执行器模型等价为一个二阶系统，其传递函数为：

(12)

式中:自然频率ωact=8.8 rad·s-1(带宽1.4 Hz)，阻尼比ξact=0.8。

2.2 性能指标

总的来说，可将系统性能指标分为两大类，即时域性能指标和频域性能指标。

1)时域性能指标

在大多数情况下，性能约束根据时域指标给出，如最大超调量、上升时间、幅值限制等。这些约束可以通过闭环系统对固定测试输入信号的响应来表示，包括衰减系数、阻尼比、自然频率等。

时域性能指标与参考模型约束效果上类似，实际使用中可根据情况自由选择。

2)频域性能指标

另一方面，系统对外界干扰的抑制能力或对结构不确定性的鲁棒性更容易通过相应闭环系统传递函数的H∞范数表示。本文考虑了如下两种频域性能指标：

(1)闭环系统鲁棒性能指标

闭环系统的性能通过传递函数的H∞范数来衡量，则以闭环系统鲁棒性能为指标的参考模型跟踪问题可描述为:

(13)

式中：TNzc→z1表示跟踪误差信号z1到过载指令信号Nzc的传递函数，P为开环系统，K为待优化控制器参数。式(13)表示通过调整控制参数K使得系统P内部稳定，且传递函数TNzc→z1的H∞范数最小。

该方案从H∞范数的角度出发，力求将参考动力学模型与闭环飞行器之间的差异最小化。如果闭环系统在整个频域内与参考模型完全匹配，则最优H∞范数J的值为零。然而，由于物理限制，完全匹配在实际中是不可行的，因此J的值总大于零。

(2)拉起机动下控制舵面最大偏转角及角速率约束

(14)

2.3 优化问题构建及求解

控制舵面一体化设计旨在寻找同时满足：1)闭环系统飞行品质约束；2)闭环系统鲁棒性能约束；3) 最大过载机动时舵面偏转角度及偏转角速率约束的最小舵面尺寸ηopt及其对应的控制参数K。

(15)

H∞范数随系统参数的变换呈现出非凸、非光滑的特性，即优化问题是一个非凸、非光滑问题，需要采用非光滑优化算法进行求解。本节将对非光滑优化算法进行简单介绍，所涉及的具体理论可参考文献[21]。

为方便描述，可将优化问题一般化表述为：

minf(K)
s.t.g(K)≤0

式中:目标函数f、约束条件g可以是时域、频域约束的综合。

为求解该约束问题，引入如下进度函数：

F(K+,K)=

max{f(K+)-f(K)-μg(K)+;

g(K+)-g(K)+}

式中:μ为大于0的固定常数，g( )+=max{g,0}。K表示当前迭代，K+为下一步迭代。

根据文献[22]可知，进度函数F的临界值K*也是原始优化问题的临界值，其满足条件0∈∂1F(K*,K*)，其中∂1F(K*,K*)表示函数F在点K*处的Clarke次微分。该临界点可由迭代下降方法求解。首先，如果当前迭代不是进度函数的临界值，即0∉∂1F(K,K)，则在K的邻域内存在一点K+，满足F(K+,K)

迭代求解过程中，点K处的下降步长ΔK通过求解如下最小值问题得到：

(16)

根据求解得到的下降步长，进行下一轮迭代。令K+=K+ΔK，并检查K+是否在可行范围内。若此时k+已超出了可行解范围，则需进行回归搜索直至K+=K+l·ΔK满足如下Armijo条件[21]：

F(K+lΔK,K)-F(K,K)<γ0lF′(·,K)(K;ΔK)

(17)

其中，0<γ0<1，0

综上所述，可将非光滑优化算法流程归纳如下：

1) 初始化，选择一个闭环稳定控制器K1；

3) 求解式(6)获得下一步迭代方向ΔK；

4) 线性搜索，寻找满足约束条件的参数l；

5) 令Kj+1=Kj+l·ΔK，j的数值加1，并跳转至步骤2)。

需要注意的是，非光滑优化算法求解的是局部最小值，在仿真过程中可通过随机选取几组初值的方法，减弱算法对初值的敏感性。

3 仿真分析

3.1 气动系数代理模型

分别令η=0.2,0.4,0.6,0.8,1，构建不同升降副翼尺寸下的空天飞行器低速气动数据库。并基于1.3节提出的代理模型自动获取算法，得到飞行器气动特性拟合表达式(拟合优度阈值为0.99)。其中，洁净体飞行器气动系数Ci(clean)的拟合表达式为：

(18)

式中：多项式系数a*，b*，c*如表1所示。

表1 洁净体气动特性代理模型系数Table 1 Coefficients of aerodynamic surrogate model (clean configuration)

同理，可获得不同尺寸舵面偏转引起的气动特性增量拟合表达式为：

(19)

式中:不同尺寸对应的拟合表达式系数如表2所示。

表2 不同控制面气动系数Table 2 Coefficients of aerodynamic surrogate model (control surface incremental)

3.2 静稳定性分析

不同质心位置下的飞行器俯仰力矩CmC可根据如下公式计算得到：

CmC=CmR+(xR-xC)(CDsinα+CLcosα)/c+(zR-zC)(-CDcosα+CLsinα)/c

(20)

式中:(xR,0,zR)为力矩参考点坐标，(xC,0,zC)为质心位置坐标。

定义质心位置沿机体轴方向的坐标xC与机身长度Lf的比值为归一化质心位置参数：

(21)

如图6所示，随着质心位置的后移，配平舵面和配平攻角均随之减小。进一步后移质心，飞行器变为静不稳定，配平舵偏为正且随着质心后移而增大，配平攻角减小。

图6 起飞状态随质心位置变化关系Fig.6 Take-off states for varying center of gravity

3.3 控制舵面一体化优化

选取(Ma0.4,h=9 km)作为低动压典型工作点，采用如图5所示的参考模型跟踪控制方案进行控制舵面一体化设计。

1)确定质心位置

(22)

此时的优化问题转换为在不同质心位置下，求解满足约束条件(13)和(14)的最小控制舵面尺寸ηopt及对应的控制参数。

仿真过程中所涉及的约束值如下：

(1)最大舵面偏转角：δemax=30°，需要注意的是该最大舵面偏转角包含配平舵偏，即如果配平舵面角较大，则可用于机动的舵面偏转角余量就偏小;

运用2.3节所介绍的非光滑优化算法分别对三个质心位置下的一体化优化问题进行求解，结果如表3所示。

表3 不同质心位置下的舵面优化结果Table 3 Optimization results for different

图7 不同质心位置下的跟踪结果Fig.7 Tracking results for varying center of gravity

2)包线内的控制舵面一体化设计

确定质心后，式(9)中的系数可写作飞行状态、控制舵面尺寸的函数，即：

(23)

此时，优化问题转换为：寻找最小控制舵面尺寸，使得对于包线范围内任意点，均存在符合约束条件的稳定控制器。

包线内的优化结果如图8所示，飞行器在各工作点满足稳定性、鲁棒性、偏转角度及偏转角速率的最小舵面尺寸的最大值在(Ma0.4,h=9 km)处取到，所以空天飞行器舵面优化的结果为ηopt=0.79。最后，对优化后的飞行器构型进行起飞性能验证可知，飞行器起飞配平攻角为7°，配平舵面偏转角为2.4°，即优化后的空天飞行器能够满足水平起飞的任务需求。

图8 不同状态点的舵面优化结果Fig.8 Optimization results for varying operation points

4 结 论

本文提出了一种基于代理模型的空天飞行器控制舵面设计及优化算法。首先，基于鸽群算法实现了可变构型空天飞行器气动特性参数代理模型的自动获取。然后，构建了飞行器纵向参考模型跟踪控制器。将飞行品质约束下的飞行器控制舵面设计问题转换成多约束多参数优化问题，并运用非光滑优化算法对该优化问题进行求解。最后，将本文所提的方法用于目标空天飞行器升降舵面设计。仿真结果表明，方算法能够实现给定控制结构下的空天飞行器控制舵面—控制参数一体化优化。相较于按经验设计的初始升降副翼构型，优化后的升降副翼面积减小了约21%。综上所述，本文提出的基于代理模型的飞行器控制舵面一体化优化方法能够与工程中常用的控制结构有效结合，具有较强的应用价值。