本发明涉及一种考虑飞行状态约束的亚轨道飞行器运载性能优化方法，属于亚轨道飞行器设计优化。
背景技术：
1、亚轨道飞行器是一种可以抵达空间边界(约100km)的新型飞行器。亚轨道飞行器由于其可复用性、易维护性及高效费比，已受到旅游业、运输业等行业的广泛关注(chang，2020)。近年来，spaceship one、spaceliner以及lynx等亚轨道飞行器已得到充分发展。然而，亚轨道飞行器在其再入返回地球大气过程中，常常面临强过载、大热流、高动压等多种挑战。为应对上述挑战，国内研究人员对亚轨道飞行器分系统设计优化做了大量研究。wang(2018)使用直接搜索区域法优化飞行器热防护系统，使得飞行过程中的驻点热流最小化。zhang(2014)基于粒子群优化算法，针对亚轨道飞行器的再入返回过程设计了一种最优攻角程序。luo(2023)基于准一维热力学方法建立rbcc发动机性能模型，并基于nsga-ⅱ方法探索比冲增大比和推力-面积比的pareto前沿。然而，亚轨道飞行器总体性能受到几何特征、气动等多个分系统的影响，故亚轨道飞行器设计优化是一类典型的多学科设计优化(multidisciplinary design optimization，mdo)问题。此外，亚轨道飞行器设计优化常常涉及气动力热等高耗时分析模型，若仍采用传统进化类优化算法求解亚轨道飞行器设计优化问题，将导致计算成本剧增。
2、为了提升亚轨道飞行器总体方案阶段的设计效率，十分有必要发展一种考虑飞行状态约束的亚轨道飞行器运载性能优化方法，从而提高设计优化效率、缩短设计周期，从而在总体设计阶段能够快速实现亚轨道飞行器设计方案的快速优选与修改，为亚轨道飞行器系统方案论证与总体设计提供科学的依据与参考。
3、为了更好地说明本发明的技术方案，下面对所涉及的kriging代理模型相关数学基础做一定介绍。
4、kriging代理模型是一种无偏最优估计插值模型，由全局模型和局部偏差项叠加而成，如式(1)所示
5、
6、式中，μ(x)为设计空间内的多项式全局近似模型，当近似对象数值特征未知时，一般取常数μ；局部偏差项z(x)是均值为零、方差为σ2、协方差非零的随机过程，表征在全局近似模型基础上的局部偏差，相关函数与协方差矩阵可分别由式(2)与式(3)表示：
7、
8、cov[z(xi),z(xj)]＝σ2r[r(xi,xj)] (3)
9、式中，r为对称相关矩阵；r(·,·)为相关函数；nv为设计变量维度；xi为训练样本点。μ与σ2的估计值可由最小二乘法根据式(4)得到：
10、
11、式中，nkrg表示构建代理模型的样本点数；i表示nkrg维的单位行向量。通过极大似然法求解式(5)中的子优化问题可确定相关系数θk
12、
13、此外，kriging代理模型可以预估任意点x处预测值的方差s2(x)，从而评估代理模型的近似误差。s2(x)可表示为：
14、
15、式中，:r表示相关函数向量，其具体形式如下
16、
技术实现思路
1、为了解决亚轨道飞行器运载性能提升问题，本发明的目的是提供一种考虑飞行状态约束的亚轨道飞行器运载性能优化方法，该方法在保证亚轨道飞行器方案满足过载、热流、动压等约束条件下，引导亚轨道飞行器设计优化过程快速向最优方案收敛，提升亚轨道飞行器有效载荷。本发明能够适用于不同任务需求的亚轨道飞行器总体方案设计优化问题，提升亚轨道飞行器的性能与设计效率。
2、本发明的目的是通过下述技术方案实现的。
3、本发明公开的考虑飞行状态约束的亚轨道飞行器运载性能优化方法，综合考虑亚轨道飞行器几何、气动等重点学科间耦合关系，建立亚轨道飞行器多学科分析模型与多学科设计优化问题模型。本发明采用基于初始样本扩容机制的近似优化策略，以亚轨道飞行器有效载荷最大为目标函数，并考虑过载、动压等飞行状态约束，对亚轨道飞行器设计变量进行优化。在优化过程中，采用kriging代理模型代替亚轨道飞行器高耗时分析模型，并由约束差分进化算法求解代理模型伪最优解。结合综合改善概率与重点设计空间方法，实现代理模型管理与更新，高效引导亚轨道飞行器运载性能优化过程收敛至全局最优解，即实现考虑飞行状态约束的亚轨道飞行器运载性能优化。本方法克服了传统优化方法求解效率低、全局收敛性差的问题，能够高效挖掘亚轨道飞行器高性能设计方案，并缩短设计周期。
4、本发明公开的考虑飞行状态约束的亚轨道飞行器运载性能优化方法，包括如下步骤：
5、步骤一：确定亚轨道飞行器运载性能优化模型，并确定基于初始样本扩容机制的近似优化策略初始参数。
6、步骤一的具体实现方法如下：
7、以亚轨道飞行器有效载荷质量mpl最大为目标函数，并考虑起飞质量mtakeoff、轴向过载na、法向过载nn、再入段末端速度vend、再入段射程r、再入段末端升阻比cl/cd、飞行动压q、再入段末端速度倾角θtouch、驻点热流密度qstag及机翼前缘热流密度qwing约束，建立亚轨道飞行器亚轨道飞行器运载性能优化模型，如式(8)所示：
8、find x＝[cv,1,cv,2,ch,1,ch,2,lint,lext,broot,btip,χint,χext,lh,wr,fu,fa,kα]
9、min f(x)＝-mpl
10、
11、式中，x表示设计变量，具体包括飞行器头部轮廓垂直方向控制点1坐标cv,1、飞行器头部轮廓垂直方向控制点2坐标cv,2、飞行器头部轮廓水平方向控制点1坐标ch,1、飞行器头部轮廓垂直方向控制点2坐标ch,2、内翼段半展长lint、外翼段半展长lext、翼根弦长broot、翼尖弦长btip、内翼段后掠角χint、外翼段后掠角χext、头部长度lh、后体宽度wr、不可预期燃油质量比fu、附加燃油质量比fa及攻角下降速率kα作为设计变量。xlb与xub分别表示设计变量下界与上界。
12、所述基于初始样本扩容机制的近似优化策略初始参数包括初始采样数ndoe、约束阈值缩放因子η、初始样本扩容数nl、局部搜索阶段新增样本点数na以及最大模型调用次数
13、步骤二：基于亚轨道飞行器设计变量，建立亚轨道飞行器的参数化几何主模型，并获取亚轨道飞行器几何特征。亚轨道飞行器几何特征包括参考面积、头部驻点最大曲率半径、头部驻点最小曲率半径、机翼前缘曲率半径、机翼前缘后掠角、总体积、总表面积、机翼面积及尾翼面积。
14、步骤二的具体实现方法如下：
15、基于如式(9)所示的nurbs方法对亚轨道飞行器头部、机身、机翼、后体及尾翼进行参数化建模，得到参数化几何主模型：
16、
17、式中，di为曲线控制点、ri,p(u)为定义在u∈[0,1]上的有理基函数、ωi为权重系数、ni,p(u)为第i个p阶b样条基函数、n表示节点数量。ni,p(u)可由式(10)、式(11)表示；当节点u与阶数p确定时，b样条基函数可唯一确定
18、
19、式中，ni,p(u)为第i个p阶b样条基函数、ui为第i个节点。
20、由式(9)-(11)生成头部、机身、机翼、后体及尾翼的参数化外形；并测量几何特征(参考面积、头部驻点最大曲率半径、头部驻点最小曲率半径、机翼前缘曲率半径、机翼前缘后掠角、总体积、总表面积、机翼面积及尾翼面积)，作为获取升力系数、阻力系数、头部驻点热流密度、机翼内翼段前缘热流密度、机翼外翼段前缘热流密度、再入返回质量及弹道模型的输入。
21、步骤三：依据步骤二中的亚轨道飞行器几何主模型及参考面积，获取亚轨道飞行器升力系数及阻力系数。
22、步骤三的具体实现方法如下：
23、①.基于步骤二中的参数化几何主模型及计算域，划分非结构化网格；并于头部驻点及机翼处添加局部细化处理。
24、于头部及机翼处添加边界层网格，其首层网格高度如式(12)所示
25、
26、式中，无量纲壁面距离y+≈1、lr为飞行器参考长度、re为雷诺数，由式(13)表示
27、
28、式中，ρ∞为大气密度、v∞为飞行工况中的来流速度、μ为大气粘度；
29、经过局部细化处理与添加边界层网格后，得到网格a。
30、②.基于步骤二测量得到的亚轨道飞行器的参考面积及①中的网格a，由cfd方法求解亚轨道飞行器升力系数及阻力系数。在cfd求解中采用密度基求解器及能量方程，并将来流设置为压力远场边界条件；使用对流迎风分裂格式进行空间离散，并采用green-gauss节点法计算梯度；此外，在开始求解前使用3级fmg初始化以提升求解收敛速度。
31、步骤四：根据步骤二中的亚轨道飞行器头部驻点最大曲率半径、头部驻点最小曲率半径、机翼前缘曲率半径、机翼前缘后掠角，以及飞行工况，获取亚轨道飞行器头部驻点热流密度、机翼内翼段前缘热流密度及机翼外翼段前缘热流密度。
32、步骤四的具体实现方法如下：
33、①.基于式(14)、式(15)计算驻点气体显焓及气体壁面焓值
34、
35、式中，cp为空气比热比、t∞为来流温度、r为理想气体常数、w为大气常数、tw为壁面温度；
36、②.基于式(16)，计算驻点热流密度值
37、
38、式中，rmin为最小头部驻点曲率半径、rmax为最大头部驻点曲率半径、ρsea为海平面大气密度、ρ∞为来流密度；
39、③.基于式(17)、式(18)，计算内、外翼段前缘热流密度值
40、
41、式中，qsp为等半径球热流、χint为机翼内翼段前缘后掠角、χext为机翼外翼段前缘后掠角、na为修正指数、α为飞行工况中的攻角；qsp与na的具体值由式(19)、式(20)给出
42、
43、式中，rl为机翼前缘曲率半径、h300k为温度为300k的空气焓值、χ为机翼后掠角。
44、步骤五：依据步骤二中的亚轨道飞行器总体积、总表面积、机翼面积、尾翼面积，获取亚轨道飞行器再入返回质量。
45、步骤五的具体实现方法如下：
46、①.基于式(21)，计算亚轨道飞行器干重
47、mdry＝mext+mp+mtps+mpl+muc+mes+mof (21)
48、式中，mext为外部结构质量、mp为推进系统质量、mtps为热防护系统质量、muc为起落架质量、mes为机上电子系统质量、mof为其他液体质量。其中，mext由式(22)表示
49、
50、式中，mfu为机身质量、mwing为机翼质量、mvt为尾翼质量、wfu为机身单位质量因子、wwing为机翼单位质量因子、wvt为尾翼单位质量因子、sfu为机身表面积、swing为机翼表面积、svt为尾翼表面积、g0为海平面重力加速度；
51、推进系统质量mp、热防护系统质量mtps、燃料质量mfuel及其他液体质量mof由式(23)-式(26)表示
52、mp＝np×wp×sin (23)
53、mtps＝wtps×swing (24)
54、mfuel＝ρfuel×vtank×kfuel (25)
55、mof＝mfuel×kof (26)
56、式中，np为发动机数量、wp为发动机单位质量因子、sin为发动机进气道面积、wtps为热防护系统单位质量因子、ρfuel为燃油密度、vtank为燃油储箱体积、kfuel为燃油装载系数、kof为燃油-其他液体比例系数；
57、②.求解式(27)中的线性方程组，获取亚轨道飞行器起飞质量mtakeoff及干重mdry
58、
59、③.基于式(28)计算亚轨道飞行器再入返回质量
60、mentry＝mdry+(fu+fa)·mfuel (28)
61、步骤六：依据步骤二中的亚轨道飞行器参考面积、步骤三中的升力系数及阻力系数、亚轨道飞行器再入返回质量以及飞行工况，建立亚轨道飞行器的弹道模型。
62、步骤六的具体实现方法如下：
63、①.基于式(29)，确定亚轨道飞行器再入返回攻角：
64、
65、式中，v(t)为飞行速度；αmax为初始再入攻角、αmin为最终再入攻角、kα为攻角下降速率、vmax为气动减速开始前的最大速度、t1为达到最大速度所需时间、t2为速度变化量为δv时经历的时间；δv由式(30)表示
66、
67、②.求解式(31)中的偏微分方程组，获取亚轨道飞行器再入返回过程中的速度、速度倾角、航向角、纬度、经度及地心距
68、
69、式中，v为飞行速度、θ为速度倾角、σ为偏航角、r为地心距、φ为纬度、λ为经度、m为为飞行器再入返回质量、gr为重力加速度地心方向分量、ωe为地球自转角速度、gω为重力加速度地球自转方向分量、ν为亚轨道飞行器倾侧角、l(v,α)为升力、d(v,α)为阻力。
70、步骤七：基于步骤一中的亚轨道飞行器运载性能优化模型，以亚轨道飞行器有效载荷质量最大为优化目标，采用初始样本扩容机制对亚轨道飞行器样本进行扩容，利用样本数据库中的所有样本信息构建目标函数与约束函数的kriging代理模型，采用kriging代理模型代替亚轨道飞行器高耗时分析模型，并由约束差分进化算法求解代理模型伪最优解；结合综合改善概率与重点设计空间方法，在满足亚轨道飞行器各学科设计要求下实现代理模型管理与更新，高效引导亚轨道飞行器运载性能优化过程收敛至全局最优解。
71、步骤七的具体实现方法如下：
72、①.确定基于初始样本扩容机制的近似优化策略初始参数，包括初始采样数ndoe、约束阈值缩放因子η、初始样本扩容数nl、局部搜索阶段新增样本点数na以及最大模型调用次数
73、②.采用拉丁超方试验设计方法在设计空间内构造ndoe个初始样本点，并调用真实分析模型计算初始样本点处的模型响应值；将所有样本点及响应值加入样本数据库中；设置优化迭代次数iter为1；
74、③.基于初始样本扩容机制，获取nl个样本点；
75、(1).基于步骤②的样本数据库，建立如式(32)所示的随机森林分类器
76、
77、式中，nrf为构建随机森林分类器的样本点数、为第i个样本点的分类值、ci为第i个样本点的约束违背度、为第j个约束函数的归一化取值、为第j个约束函数的约束阈值、pthresh为约束违背阈值；pthresh表示为
78、pthresh＝ci,min+η·(ci,max-ci,min) (33)
79、式中，ci,min为最小约束违背度、ci,max为最大约束违背度、η为约束阈值缩放因子；
80、在设计空间中生成大量简单样本点，并由构建的随机森林分类器对简单样本点进行分类；若分类值为1，则记录该简单样本点为优质样本点；
81、(2).统计优质样本点数；若无优质样本点，则将约束阈值上调10％；否则由k-means方法获取优质样本点聚类中心；
82、(3).在优质样本点聚类中心与步骤②的样本数据库中具有最小目标函数响应值的样本点之间，随机生成一个样本点；调用真实分析模型计算该样本点真实模型响应值，并将该样本信息加入步骤②的样本数据库中；检查新增样本点数是否达到nl，若是，扩容机制终止；否则返回步骤(1)，扩容过程继续；
83、④.利用样本数据库中的所有样本信息构建目标函数与约束函数的kriging代理模型，并由约束差分进化方法对构造的kriging代理模型进行优化，得到伪最优解调用真实分析模型计算伪最优解真实模型响应值，并将伪最优解样本信息加入样本数据库中；
84、⑤.求解式(34)所示的子优化问题，在设计空间内获取综合改善概率最大的样本点xpi
85、
86、式中，pi(·)为目标函数改善概率、p(gi(·)≤0)为第i个约束满足概率、dm(x,xp-min)为样本点与当前样本数据库中最优综合改善概率样本点的曼哈顿距离；pi(·)、p(gi(·)≤0)及dm(x,xp-min)分别由式(35)、式(36)、式(37)表示
87、
88、式中，φ(·)为标准正态分布的分布函数、fmin为样本数据库中目标函数最小值、为目标函数的kriging代理模型预测值、为第i个约束函数的kriging代理模型预测值、s(·)为目标函数的kriging代理模型预测方差、sg,i(·)为第i个约束函数的kriging代理模型预测方差；调用真实分析模型计算该综合改善概率最大的样本点xpi真实模型响应值，并将xpi样本信息加入样本数据库中；
89、⑥.基于当前样本数据库，构建重点采样空间；重点采样空间半径如式(38)所示
90、
91、式中，xopt,k为第k次迭代最优可行解、xopt,k-1为第k-1次迭代最优可行解、为迭代次数为1时由逐一校验法所得误差最大的样本点；若xopt,k的最优性相较于xopt,k-1有所改善，则将xopt,k作为重点设计空间中心，否则将xopt,k-1作为重点设计空间中心；
92、在重点采样空间内随机生成一个样本点xsss，调用真实分析模型计算xsss真实模型响应值，并将xsss样本信息加入样本数据库中；
93、⑦.判断当前样本数据库中样本点数是否达到若达到则停止优化；否则令优化迭代次数iter为iter+1，转到步骤七④，更新目标函数与约束函数的kriging代理模型，优化过程继续。
94、步骤八：输出当前样本数据库中的可行最优解，作为亚轨道飞行器优化结果，即实现考虑飞行状态约束的亚轨道飞行器运载性能优化。
95、还包括步骤九：将步骤一至步骤八所述的考虑飞行状态约束的亚轨道飞行器运载性能优化方法，应用于亚轨道飞行器运载性能优化，提升亚轨道飞行器的性能与设计效率；根据步骤八得到的亚轨道飞行器优化结果，使亚轨道飞行器在起飞质量、最大轴向过载、最大法向过载、再入段末端速度、再入段射程、再入段末端升阻比、最大飞行动压、再入段末端速度倾角、最大驻点热流密度、最大机翼前缘热流密度约束下实现运载性能最大化。所述的亚轨道飞行器运载性能优化包括单级水平起降亚轨道飞行器运载性能优化、两级水平起降亚轨道飞行器运载性能优化、垂直起降亚轨道飞行器运载性能优化等领域。
96、有益效果：
97、1、针对多学科分析模型计算耗时，且有限计算成本下数据挖掘不充分问题，本发明公开的考虑飞行状态约束的亚轨道飞行器运载性能优化方法，综合考虑亚轨道飞行器几何主模型、气动力热模型、质量估算模型及弹道模型之间的参数传递关系，建立了考虑过载、动压、热流等飞行状态约束的亚轨道飞行器多学科分析模型。采用基于初始样本扩容机制的近似优化策略，实现亚轨道飞行器优质设计方案的高效挖掘。本发明公开的考虑飞行状态约束的亚轨道飞行器运载性能优化方法结合一种基于初始样本扩容机制的近似优化策略，对亚轨道飞行器设计变量进行扩容，利用样本数据库中的所有样本信息构建目标函数与约束函数的kriging代理模型，采用kriging代理模型代替亚轨道飞行器高耗时分析模型，并由约束差分进化算法求解代理模型伪最优解；结合综合改善概率与重点设计空间方法，在满足亚轨道飞行器各学科设计要求下实现代理模型管理与更新，高效引导亚轨道飞行器运载性能优化过程收敛至全局最优解。由该方法所得优化方案的有效载荷相较于初始方案可提升171kg以上。
98、2、本发明公开的考虑飞行状态约束的亚轨道飞行器运载性能优化方法，采用基于初始样本扩容机制的近似优化策略对亚轨道飞行器设计变量进行扩容，利用样本数据库中的所有样本信息构建目标函数与约束函数的kriging代理模型，采用kriging代理模型代替亚轨道飞行器高耗时分析模型，并由约束差分进化算法求解代理模型伪最优解；结合综合改善概率与重点设计空间方法，在满足亚轨道飞行器各学科设计要求下实现代理模型管理与更新，高效引导亚轨道飞行器运载性能优化过程收敛至全局最优解，进一步提高亚轨道飞行器运载性能优化的效率。所述全局最优解对应亚轨道飞行器优化结果，根据亚轨道飞行器优化结果使亚轨道飞行器在飞行状态约束下实现运载性能最大化。
1.考虑飞行状态约束的亚轨道飞行器运载性能优化方法，其特征在于：包括如下步骤，
2.如权利要求1所述的方法，其特征在于：步骤一的具体实现方法如下：
3.如权利要求1所述的方法，其特征在于：步骤二的具体实现方法如下：
4.如权利要求1所述的方法，其特征在于：步骤三的具体实现方法如下：
5.如权利要求1所述的方法，其特征在于：步骤四的具体实现方法如下：
6.如权利要求1所述的方法，其特征在于：步骤五的具体实现方法如下：
7.如权利要求1所述的方法，其特征在于：步骤六的具体实现方法如下：
8.如权利要求1所述的方法，其特征在于：步骤七的具体实现方法如下：
9.如权利要求8所述的方法，其特征在于：步骤七所述目标函数kriging代理模型为：有效载荷质量mpl的kriging代理模型；步骤七所述约束函数kriging代理模型为：起飞质量mtakeoff的kriging代理模型、轴向过载na的kriging代理模型、法向过载nn的kriging代理模型、再入段末端速度vend的kriging代理模型、再入段射程r的kriging代理模型、再入段末端升阻比cl/cd的kriging代理模型、飞行动压q的kriging代理模型、再入段末端速度倾角θtouch的kriging代理模型、头部驻点热流密度qstag的kriging代理模型以及机翼前缘热流密度qwing的kriging代理模型。
10.如权利要求1、2、3、4、5、6、7、8或9所述的方法，其特征在于：还包括步骤九：将步骤一至步骤八所述的考虑飞行状态约束的亚轨道飞行器运载性能优化方法，应用于亚轨道飞行器运载性能优化，提升亚轨道飞行器的性能与设计效率；根据步骤八得到的亚轨道飞行器优化结果，使亚轨道飞行器在起飞质量、最大轴向过载、最大法向过载、再入段末端速度、再入段射程、再入段末端升阻比、最大飞行动压、再入段末端速度倾角、最大驻点热流密度、最大机翼前缘热流密度约束下实现运载性能最大化。所述的亚轨道飞行器运载性能优化包括单级水平起降亚轨道飞行器运载性能优化、两级水平起降亚轨道飞行器运载性能优化、垂直起降亚轨道飞行器运载性能优化等领域。
