一种针对高耗时约束的飞行器快速近似优化方法与流程

本发明涉及一种针对高耗时约束的飞行器快速近似优化方法，属于飞行器设计中的多学科优化技术领域。

背景技术：

随着科技的进步，工程优化问题越来越复杂，许多分析和仿真软件都被应用于设计和研究中，但是这些分析和仿真问题大多都是高精度分析模型，例如结构分析中采用的有限元分析(fea)模型、气动分析中使用的计算流体力学(cfd)分析模型等。高精度分析模型在提高分析精度和可信度的同时也带来了计算耗时的问题，虽然当今计算机软硬件技术已经取得了很大的进步，但是调用高精度分析模型完成一次分析仍然极其耗时，例如使用cfd模型完成一次气动仿真分析仍需要数小时甚至数十小时；其次，现代工程设计问题往往涉及多个相互耦合的学科。譬如，航天器设计涉及结构、姿态控制、轨道转移、位置保持、电源等学科，各学科相互影响，相互制约，飞行器的性能是各学科耦合的综合体现。由于学科之间的耦合关系，工程设计问题的系统分析表现为多学科分析。本质上，多学科分析过程是一个典型的非线性求解过程，每次多学科分析都需要进行多次迭代，计算耗时，如果各学科都采用高精度分析模型，计算量将非常庞大。再次，工程优化设计过程中需要通过反复迭代才能收敛到局部或全局最优解，而每次迭代都需要进行多次工程设计问题的多学科分析，计算成本将进一步增加。传统的梯度算法往往只能找到分析问题的局部最优解，不具备全局搜索能力。

为了获得工程设计问题的全局最优解，常常直接使用具有全局搜索能力的优化算法，例如遗传算法(ga)、模拟退火算法(sa)等。但是，与传统的梯度算法相比，全局优化算法所需计算量更大。例如，使用遗传算法对一个分析模型进行优化通常需要调用成千上万次分析系统。对于大量采用高精度学科分析模型的复杂现代工程设计问题，传统的全局优化策略的计算成本过大甚至难以接受。此外，目前多数高精度分析模型都采用商业cae软件建立的黑盒模型(black-boxmodel)，与优化算法(优化器)的接口相当困难。

为了降低传统的全局优化算法对现在工程设计问题优化过程中所涉及的高计算量，基于代理模型的优化方法被逐渐的为人们所研究。其本质上就是构造与高精度分析模型近似，但计算成本更低的数学代理模型，并将该模型用于优化。由于高精度分析模型计算一次所需时间的量级为小时，而代理模型计算一次所用时间的量级仅为秒甚至毫秒，因此与高精度分析模型的计算时间相比，构造代理模型以及基于代理模型优化的计算时间往往可以忽略不计。在最近10年里，许多公司都开始研究并推动了近似代理模型技术在设计和优化领域的应用，例如：engineous软件公司研发isight，vanderplaatsr&d公司研发的visualdoc，lms国际组织研发的optimus，phoenix公司研发的modelcenter和波音公司研发的designexplorer。

代理模型可以分为静态代理模型和动态代理模型。静态代理模型通过一次采取足够多的样本点构造代理模型，在优化过程中代理模型保持不变；而动态代理模型是序列采取样本点，在每次优化迭代过程中根据已知信息逐步改进和更新代理模型，直至优化收敛。与静态代理模型相比，动态代理模型在优化效率和结果精度方面更具有优势。

克里金(kriging)模型是常用的代理模型方法之一，能够较好地折中了计算效率和近似精度，从而获得广泛的使用。jones等提出了高效全局优化方法(efficientglobaloptimization，ego)。它在设计空间中获取少量样本点来构造kriging模型，并在误差较大处增加样本点更新代理模型直到找到近似的全局最优解。ego仅仅需要较少的模型调用次数就能找到较好的解，从而大大提高了优化的效率。尽管如此，ego还存在一定的不足——不能够处理约束优化问题。本发明将自适应罚函数方法与ego结合，提出了一种针对高耗时约束的飞行器快速近似优化方法，从而能够高效解决约束优化问题。

为了更好的说明本发明的技术方案，下面对应用到的物理规划理论相关数学基础进行具体介绍：

kriging模型：

kriging模型是1951年由南非地质学者daniekrige提出的一种估计方差最小的无偏最优估计模型，由全局模型和局部偏差叠加而成。从统计意义上而言，kriging模型是从变量相关性和变异性出发，在有限区域内对区域化变量的取值进行无偏、最优估计的一种方法；从插值角度讲，kriging模型是对空间分布的数据求线性最优、无偏内插估计的一种方法。

kriging模型的数学表达式可表示为

f(x)＝g(x)+z(x)(1)

其中，g(x)未关于x的函数，是设计空间范围内的全局近似模型，z(x)是均值为零、方差为σ²、协方差不为零的随机过程，z(x)提供在全局近似模型基础上的局部偏差。

通常，g(x)可取常值β，则式(1)变为

y(x)＝β+z(x)(2)

z(x)的协方差矩阵可表示为

cov[z(xi),z(xj)]＝σ²r[r(xi,xj)](3)

式中，r为相关矩阵，r为相关函数，i＝1,2,....,ns，j＝1,2,....,ns(ns为样本点数)。r为对称矩阵，且其对角线元素为1。

通常，r取高斯相关函数可表示为

其中，nv是设计变量个数,θk为未知相关参数向量，为了降低复杂度，通常取常值θ，这样，式(4)可表示为：

引入另一相关向量r(x)为

根据统计学知识，kriging模型可表示为

f(x)＝β+r^t(x)r^-1(y-gβ)(7)

式中，β是未知参数，σ²和r都是θ的函数，y是由采样点响应值组成的ns维列向量，当g(x)取常值时，g是元素全为1的n维列向量，即

g＝[1,1,...,1]^t(8)

β和σ²的最小二乘估计可由式(9)求出。

相关参数θ可通过极大似然估计(maximumlikelihoodestimation)，求解式(10)中的一维优化问题得到

高效全局优化方法：

高效全局优化方法(efficientglobaloptimization，ego)是一种贝叶斯全局优化算法，其核心是在误差较大处增加样本点更新代理模型直到找到近似的全局最优解。ego优化的基本流程为：利用实验设计方法生成初始样本点并构造kriging代理模型，以ei函数值为目标并将最大点作为校正点，采用高精度模型对校正点进行评估，将校正点及其响应值加入样本点数据库，重复此过程直至收敛。其中kriging模型的数学表达式为

f(x)＝s(x)+z(x)(11)

式中：s(x)为x的函数，是设计空间范围内的全局近似模型；z(x)是均值为零、方差为、协方差不为零的高斯随机过程，提供全局近似模型基础上的局部偏差。ei函数定义为

式中：φ为标准正态分布函数；φ为标准正态分布的概率密度函数；fmin为目标函数最小值；与s分别为样本点x处代理模型的预测值与标准差。式(12)的第一项将当前最小目标函数值与预测值的差乘以预测值小于当前最小目标函数值的概率，当预测值小于当前最小目标函数值时，第一项会变大，从而使搜索重点集中在当前最小目标函数值附近强化局部搜索能力。式(12)的第二项是预测值标准差与概率密度函数的积，当预测值与当前最小目标函数接近时，第二项会变大，从而使搜索重点集中在预测精度较差区域强化全局搜索能力。因此采用ei函数确定校正点的方法同时考虑了局部和全局的搜索能力，兼顾了收敛速度和全局搜索能力。

技术实现要素：

针对使用传统的全局优化算法在对高精度分析模型优化过程中计算费时，以及标准的高效全局优化方法不能求解约束问题的缺陷。本发明公开的一种针对高耗时约束的飞行器快速近似优化方法要解决的技术问题是，提高针对高耗时约束的飞行器快速近似优化方法在全局优化算法中对高精度分析模型优化过程的计算效率，并能够解决标准的高效全局优化方法不能求解约束问题的缺陷。

本发明公开的一种针对高耗时约束的飞行器快速近似优化方法，采用基于kriging模型的代理模型，同时采用自适应罚函数方法处理约束条件，其设计原理为：本发明在处理带约束的工程设计问题优化过程中，通过在每次迭代过程中根据已知样本点信息、样本点对应的约束以及每次更新的罚系数来构造代理模型，通过对代理模型进行优化并新增一个样本点，随后调用真实分析模型获取对应的函数值及约束值，若在优化所得样本点处满足约束则减小罚系数从而使得下一次优化迭代侧重于改善最优性，若在优化所得样本点处不满足约束则增大罚系数使得下一次优化迭代更倾向于提高可行性，如此进行反复迭代直至满足预设的最大模型调用次数，最终输出优化结果，完成针对高耗时约束的飞行器快速近似优化。

本发明公开的一种针对高耗时约束的飞行器快速近似优化方法，包括如下步骤：

步骤1，定义初始输入参数。

根据实际问题，定义分析模型、优化变量设计空间、约束个数以及初始样本点个数。

步骤2，生成初始样本点并计算模型响应值，并保存到设计样本点数据库中。

根据步骤1中定义的变量设计空间和初始样本点个数，采用采样方法生成相应的样本点，通过调用真实模型计算出每一个初始样本点所对应的真实模型响应值及约束值，并将上述的样本点对应的真实模型响应值和约束值保存到设计样本点数据库中。

步骤2中所述的采样方法优选拉丁超方采样方法能够使样本点均匀的分布在设计空间中，从而提高构造模型的精度。

步骤3，计算归一化因子。

根据步骤2中得到的样本点及其对应的真实模型响应值和约束值计算目标函数的归一化因子favg以及约束的归一化因子cavg，计算表达式如下：

favg＝fsum/npoints(13)

cavg＝csum/npoints(14)

其中npoints为当前样本点的个数，fsum为当前所有样本点目标函数值绝对值的和，csum为当前所有样本点约束值绝对值的和。

步骤4，计算伪目标函数值。

将步骤2得到的设计样本点数据库中所有样本点及其所对应的真实模型响应值和约束值提取出来，对约束进行归一化处理，其计算表达式如下：

cnorm(x)＝c(x)/cavg(15)

其中c(x)为样本点x处的约束值，cnorm(x)为归一化之后的约束值。在完成对约束的处理之后，对目标函数进行处理，其表达式如下：

fnorm(x)＝f(x)/favg(16)

其中f(x)为样本点x处的目标函数值，fnorm(x)为归一化之后的目标函数值。在完成目标函数值和约束值的归一化之后，构造包含约束值信息的伪目标函数值。当约束值c(x)小于或等于零时，目标函数值保持不变；当约束值c(x)大于零时，对目标函数值处理如下：

minf(x)＝fnorm(x)+λmax(cnorm(x),0)(17)

其中，λ为罚系数，初值为1，f(x)为根据归一化目标函数值和约束值构造的伪目标函数值。对所有样本点进行如上所示的归一化处理。

步骤5，构造ei函数并进行优化得到样本点xmin。

提取步骤2用到的所用样本点以及所构造的伪目标函数值，构造ei函数并以ei函数值最大为目标进行优化，ei函数的表达式如下：

式中：φ为标准正态分布函数；φ为标准正态分布的概率密度函数；fmin为当前所有伪目标函数值的最小值；与s分别为样本点x处代理模型的预测值与标准差。对ei函数进行优化后得到样本点xmin。

步骤6，计算样本点xmin的函数值及约束值。

调用真实分析模型计算xmin处的目标函数响应值fmin以及约束值cmin，将样本点及其对应的目标函数响应值和约束值保存到样本点数据库中。

步骤7，更新罚因子。

首先，将对cmin进行归一化处理，其表达式如下：

cnorm＝cmin/cavg(19)

其中cnorm为归一化之后的约束值，求出cnorm中的最大值即对应的最大约束违背度smax，如果最大约束违背度大于约束容差stol，则根据预先定义的增长因子μ1>1增大罚系数，从而使得下一次优化迭代更倾向于提高可行性；相反，如果最大约束违背度小于或等于容差，则根据初始定义的下降因子μ2≤1减小罚系数，从而使得下一次优化迭代侧重于改善最优性。另外，引入罚系数上下限约束(λmin,λmax)，一方面避免罚系数过小导致优化结果不可行；另一方面，在优化过程中避免罚系数过大所带来数值困难。自适应罚系数λ^(k+1)计算公式如下所示：

其中，增长因子μ1、下降因子μ2、上限约束λmax、下限约束λmin、约束容差stol为自适应罚系数λ^(k+1)计算需设定的系数。

作为优选，增长因子μ1取值为2；下降因子μ2取值为0.8；上限约束λmax取值为50；下限约束λmin取值为1；约束容差stol取值为0.001。

步骤8，生成伪目标函数值。

提取步骤6样本点数据库中的目标函数响应值和约束值，根据步骤7中得到的罚系数λ^(k+1)对目标函数值进行更新，更新公式如下：

minf(x)＝f(x)/favg+λmax(c(x)/cavg,0)(21)

式中，f(x)为样本点x处的伪目标函数值，f(x)为样本点x处的真实目标函数值，c(x)为样本点x处的约束值。

步骤9，在步骤8的基础上，令k＝k+1，转入步骤4进行下一次迭代直至达到最大迭代次数，最终输出优化结果，完成针对高耗时约束的飞行器快速近似优化。

还包括步骤10：将步骤1至步骤9所述的针对高耗时约束的飞行器快速近似优化方法应用于工程优化中，能够有效的降低优化计算成本，提高优化效率，有助于缩短工程优化设计的周期。

将步骤1至步骤9所述的针对高耗时约束的飞行器快速近似优化方法应用于工程优化中包括：标准工程算例方面，在满足约束条件下降低压力容器设计成本应用、在满足约束条件下减少压缩弹簧质量优化应用；卫星优化设计问题方面，在满足约束条件下降低对地观测卫星质量优化应用、在满足频率约束条件下降低卫星结构质量优化应用；导弹优化设计问题方面，在满足强度约束条件下增加导弹有效载荷优化应用、在满足燃料约束条件下增加导弹射程优化应用，本发明在飞行器系统优化工程设计中有较好的应用前景。

有益效果：

1、本发明公开的一种针对高耗时约束的飞行器快速近似优化方法，是一种动态代理模型方法，相比于静态代理模型方法，该方法在每一次的优化过程中新增一个样本点，从而不断更新构造代理模型，提高代理模型与真实分析模型的近似精度；相比于传统的高效全局优化方法，该方法引入了自适应罚函数方法来处理约束优化问题，在每一次优化过程中根据优化所得样本点不断更新罚系数并用于下一次优化，从而引导整个优化过程获得满足约束的最优解，即能够解决标准的高效全局优化方法不能求解约束问题的缺陷，拓展高效全局优化方法的适用范围。

2、本发明公开的一种针对高耗时约束的飞行器快速近似优化方法，能够克服传统的全局优化方法在工程优化问题时存在的计算耗时的缺点，并且与静态代理模型相比，本发明的优化方法可以通过调用较少次数分析模型而找到分析模型的全局最优点，能够有效的降低计算成本，提高优化效率，有助于缩短工程优化周期，在飞行器系统优化等工程设计中有较好的应用前景。

3、本发明公开的一种针对高耗时约束的飞行器快速近似优化方法，通过高效全局优化方法中的ei函数考虑局部近似精度与全局近似精度，同时具有处理约束的能力，具有全局搜索能力。

4、本发明公开的一种针对高耗时约束的飞行器快速近似优化方法，适用于各种带约束的复杂现代工程优化设计问题，对于许多工程设计和分析软件，如：飞行器气动分析软件cfd、飞行器结构分析软件nastran、飞行器结构中有限元分析软件fea等，在进行工程优化设计过程中，往往优化设计需要消耗几个小时甚至几天，而采用本发明对工程设计和分析软件中的真实模型进行近似，同时采用自适应罚函数处理约束条件，对这样的近似模型进行优化设计只需要几秒或几十秒，从而大大缩短工程优化设计周期，节约设计成本，显著提高设计效率，并且能够获得满足约束的最优设计方案。

附图说明

图1为针对高耗时约束的飞行器快速近似优化方法流程图；

图2为压力容器设计问题示意图；

图3为罚系数随函数值变化图；

图4为对地观测卫星设计结构矩阵图；

图5为对地观测卫星多学科优化问题优化结果对比箱线图，其中：图5(a)为目标函数值箱线图、图5(b)为目标和约束调用次数箱线图、图5(c)为最大约束违背度箱线图。

具体实施方式

为了更好的说明本发明的目的与优点，下面通过一个工程标准的测试算例与一个对地观测卫星多学科优化设计实例，结合附图对本发明做进一步说明，并通过与一次采样构造rbf代理模型技术结果比较，对本发明的综合性能进行验证分析。

实施例1：以在满足约束条件下降低压力容器设计成本应用为例。

压力容器设计问题的目的是在满足一定约束的前提下减少包括焊接、材料以及成型在内的成本，假设压力容器模型是工程设计中计算耗时的高精度分析模型，通过求解相应函数在设计空间中的最小值，验证c-ego的性能。本实施例的目的是提高优化设计效率，即减少求解该函数的次数。

给定的压力容器设计问题如下所示：

为了获得全局最优解，使用具有全局寻优能力的优化方法。在此使用本发明的针对高耗时约束的飞行器快速近似优化设计方法对该设计问题进行寻优。同时，为了比较该方法的效率，采用遗传算法(ga)以及处理约束的追峰采样方法(cimps)对该问题进行优化。

本实施例公开的一种针对高耗时约束的飞行器快速近似优化方法，具体实现步骤如下：

步骤1，定义初始输入参数。

根据压力容器设计问题采用m语言定义分析模型，优化变量设计空间下限为[1.0，0.625，25，25]、上限为[1.375,1.0,150,240]，约束个数为4个，初始样本点个数为20个。

步骤2，生成初始样本点并计算模型响应值。

根据步骤1中定义的变量设计空间和初始样本点个数，采用拉丁超方采样方法生成相应的样本点，通过调用真实模型计算出每一个初始样本点所对应的真实模型响应值及约束值，并将这些样本点对应的真实模型响应值和约束值保存到设计样本点数据库中。

步骤3，计算归一化因子。

根据步骤二中得到的样本点及其对应的真实模型响应值和约束值计算目标函数的归一化因子favg以及约束的归一化因子cavg，其计算结果如下：

favg＝fsum/npoints＝50850/20＝2542.5

cavg＝csum/npoints＝[14.54,6.05,14.52,2150]/20＝[0.727,0.303,0.726,107.5]

其中npoints为当前样本点的个数，fsum为当前所有样本点目标函数值绝对值的和，csum为当前所有样本点约束值绝对值的和。

步骤4，计算伪目标函数值。

将步骤2得到的设计样本点数据库中所有样本点及其所对应的真实模型响应值和约束值提取出来，对约束进行归一化处理，其计算表达式如下：

cnorm(x)＝c(x)/cavg(23)

其中c(x)为样本点x处的约束值，cnorm(x)为归一化之后的约束值。在完成对约束的处理之后，对目标函数进行处理，其表达式如下：

fnorm(x)＝f(x)/favg(24)

其中f(x)为样本点x处的目标函数值，fnorm(x)为归一化之后的目标函数值。在完成目标函数值和约束值的归一化之后，构造包含约束值信息的伪目标函数值。当约束值c(x)小于或等于零时，目标函数值保持不变；当约束值c(x)大于零时，对目标函数值处理如下：

minf(x)＝fnorm(x)+λmax(cnorm(x),0)(25)

其中，λ为罚系数，初值为1，f(x)为根据归一化目标函数值和约束值构造的伪目标函数值。对所有样本点进行如上所示的归一化处理。

步骤5，构造ei函数并进行优化。

提取步骤2用到的所用样本点以及所构造的伪目标函数值，构造ei函数并以ei函数值最大为目标进行优化，ei函数的表达式如下：

式中：φ为标准正态分布函数；φ为标准正态分布的概率密度函数；fmin为当前所有伪目标函数值的最小值；与s分别为样本点x处代理模型的预测值与标准差。对ei函数进行优化后得到样本点xmin。

步骤6，计算样本点xmin的函数值及约束值

调用真实分析模型计算xmin处的目标函数响应值fmin以及约束值cmin，将样本点及其对应的目标函数响应值和约束值保存到样本点数据库中。

步骤7，更新罚因子。

首先，将对cmin进行归一化处理，其表达式如下：

cnorm＝cmin/cavg(27)

其中cnorm为归一化之后的约束值，求出cnorm中的最大值即对应的最大约束违背度smax，如果最大约束违背度大于约束容差stol，则根据预先定义的增长因子μ1>1增大罚系数，从而使得下一次优化迭代更倾向于提高可行性；相反，如果最大约束违背度小于或等于容差，则根据初始定义的下降因子μ2≤1减小罚系数，从而使得下一次优化迭代侧重于改善最优性。另外，引入罚系数上下限约束(λmin,λmax)，一方面避免罚系数过小导致优化结果不可行；另一方面，在优化过程中避免罚系数过大所带来数值困难。自适应罚系数计算公式如下所示：

其中μ1取值为2；μ2取值为0.8；λmax取值为50；λmin取值为1；stol取值为0.001。在本实施例中，罚系数随目标函数值的变化不断发生变化，其变化过程图如图3所示。由该图可得：当c-ego获得可行解时，罚系数会随之减小，从而使下一次的迭代更加注重目标函数的改善；当c-ego获得不可行解时，罚系数会随之增大，从而使下一次优化限制在可行域内或者可行域附近。

步骤8，生成伪目标函数值。

提取步骤6样本点数据库中的目标函数响应值和约束值，根据步骤7中得到的罚系数对目标函数值进行更新，更新公式如下：

minf(x)＝f(x)/favg+λmax(c(x)/cavg,0)(29)

式中，f(x)为样本点x处的伪目标函数值，f(x)为样本点x处的真实目标函数值，c(x)为样本点x处的约束值。

步骤9，在步骤8的基础上，令k＝k+1，转入步骤4进行下一次迭代直至达到最大迭代次数100次，最终输出优化结果完成优化过程。

基于自适应罚函数的高效全局优化方法、遗传算法以及处理约束的追峰采样方法优化结果如下表所示。用目标函数调用次数(nfe)、约束条件调用次数(nce)以及模型调用总次数(ne)三个指标来衡量优化算法的效率。分别使用遗传算法(ga)、cimps和c-ego对优化问题进行连续10次优化。cimps使用算法工具包的缺省参数；ga使用matlab中ga函数的默认设置。表1中所列结果为三种算法10次优化所获得的最好结果，值得一提的是十次优化结果均满足约束。

表1压力容器设计问题优化结果对比

由表1结果可知，在结果最优性方面，ga最差，cimps最好，c-ego明显优于ga但略逊于cimps(差异小于7.9％)。在优化效率方面，虽然cimps仅需较少的目标函数调用次数，但cimps需要调用大量的约束条件来寻找可行解，导致该算法的优化效率较低。c-ego方法所需目标函数的调用次数与cimps相当，但是c-ego所需约束条件调用次数显著少于cimps。c-ego所需的模型调用总次数明显小于cimps，表明c-ego具有显著的效率优势。ga效率最低，无论模型调用次数还是约束调用次数都要远远多于另外两种算法。

通过上述优化过程可知，本实施例引入的自适应罚函数方法来处理约束优化问题，在每一次优化过程中根据优化所得样本点不断更新罚系数并用于下一次优化，从而引导整个优化过程获得满足约束的最优解，能够解决标准的高效全局优化方法不能求解约束问题的缺陷，拓展高效全局优化方法的适用范围。同时，本实施例的优化方法可以通过调用较少次数分析模型而找到分析模型的全局最优点，能够有效的降低计算成本，提高优化效率，有助于缩短工程优化周期，在飞行器系统优化等工程设计中有较好的应用前景。

实施例2：在满足约束条件下降低对地观测卫星质量优化应用为例。

对地观测卫星的设计目标是在保证分辨率和失真度的情况下，尽可能多的获取地面遥感图像和实现地球表面区域的覆盖，同时尽可能降低卫星成本。本文的多学科设计优化模型设计结构矩阵(designstructurematrix，dsm)如图4所示，其中横线代表学科输出，竖线代表学科输入，节点处符号表示学科间耦合参数。该多学科分析(multi-disciplinaryanalysis，mda)过程由轨道、控制、有效载荷、电源及结构等五个学科组成，同时包含了不同学科的耦合关系。针对该mdo问题，本实施例采用定点迭代法进行协调求解。

本实施例对地观测卫星设计指标为相机的对地观测性能最优，考虑地面分辨率和地面覆盖面积大小，设计方案可综合为加权形式的函数

f(x)＝m·(g/g0)+n·(ψ0/ψ)(30)

式中：m和n为加权系数，取m＝0.6，n＝0.4；g0和ψ0为无量纲化引入的固定值，取值分别为20m和12°。本文优化模型包含10个设计变量和14个约束，于是对地观测卫星多学科设计优化问题可描述为

式中：x为设计变量，包括轨道高度h、飞轮角动量hwh、太阳帆板面积asa、降交点时刻dt、本体高度hs、本体宽度hw、中心承力筒直径r、发动机推力fpro、蓄电池容量cs以及相机焦距f，x^lb和x^ub分别为设计变量的下边界和上边界，gi为各学科约束。

本实例所涉及的约束条件如表2所示。采用ga、cimps以及c-ego对卫星多学科设计优化问题进行连续10次优化。ga使用matlab中ga函数的默认设置；c-ego初始样本点为100个，迭代次数为200。研究中发现求解该问题时，cimps使用默认的收敛准则不能正常收敛，因此使用最大迭代次数作为cimps算法的收敛准则，当mda过程调用次数达到300次时，cimps优化终止并输出当前最好结果。采用定点迭代法进行多学科分析(multidisciplinaryanalysis,mda)确保学科间耦合状态变量的一致性。

本实施例公开的的针对高耗时约束的飞行器快速近似优化方法具体实现步骤如下：

步骤1，定义初始输入参数。

根据卫星多学科设计优化问题采用matlab的m语言定义分析模型，优化变量设计空间下限为[600,0.18,3.0,8.0,1.5,1.2,0.24,10,30,0.2]、上限为[800,4.0,20.0,11.5,5.0,2.475,0.825,200,100,1]，约束个数为14个，初始样本点个数为100个。

步骤2，生成初始样本点并计算模型响应值。

根据步骤一中定义的变量设计空间和初始样本点个数，采用拉丁超方采样方法生成相应的样本点，通过调用真实模型计算出每一个初始样本点所对应的真实模型响应值及约束值，并将这些样本点对应的真实模型响应值和约束值保存到设计样本点数据库中。

步骤3，计算归一化因子。

根据步骤2中得到的样本点及其对应的真实模型响应值和约束值计算目标函数的归一化因子favg以及约束的归一化因子cavg，其计算结果如下：

favg＝fsum/npoints＝95/100＝0.95

cavg＝csum/npoints＝[1.27e5,90.1,275,5.76e4,6.49e4,20.32,17.95,202.77,0.14,2.32,8.06e4,6.3

1e3,466,1.09]/100

＝[1.27e3,0.90,27.5,576,649,0.20,0.18,0.02,0.0014,0.023,806,63.1,46.6,0.01]

其中npoints为当前样本点的个数，fsum为当前所有样本点目标函数值绝对值的和，csum为当前所有样本点约束值绝对值的和。

步骤4，计算伪目标函数值。

将步骤2得到的设计样本点数据库中所有样本点及其所对应的真实模型响应值和约束值提取出来，对约束进行归一化处理，其计算表达式如下：

cnorm(x)＝c(x)/cavg(32)

其中c(x)为样本点x处的约束值，cnorm(x)为归一化之后的约束值。在完成对约束的处理之后，对目标函数进行处理，其表达式如下：

fnorm(x)＝f(x)/favg(33)

其中f(x)为样本点x处的目标函数值，fnorm(x)为归一化之后的目标函数值。在完成目标函数值和约束值的归一化之后，构造包含约束值信息的伪目标函数值。当约束值c(x)小于或等于零时，目标函数值保持不变；当约束值c(x)大于零时，对目标函数值处理如下：

minf(x)＝fnorm(x)+λmax(cnorm(x),0)(34)

其中，λ为罚系数，初值为1，f(x)为根据归一化目标函数值和约束值构造的伪目标函数值。对所有样本点进行如上所示的归一化处理。

步骤5，构造ei函数并进行优化。

提取步骤2用到的所用样本点以及所构造的伪目标函数值，构造ei函数并以ei函数值最大为目标进行优化，ei函数的表达式如下：

式中：φ为标准正态分布函数；φ为标准正态分布的概率密度函数；fmin为当前所有伪目标函数值的最小值；与s分别为样本点x处代理模型的预测值与标准差。对ei函数进行优化后得到样本点xmin。

步骤6，计算样本点xmin的函数值及约束值

调用真实分析模型计算xmin处的目标函数响应值fmin以及约束值cmin，将样本点及其对应的目标函数响应值和约束值保存到样本点数据库中。

步骤7，更新罚因子。

首先，将对cmin进行归一化处理，其表达式如下：

cnorm＝cmin/cavg(36)

其中cnorm为归一化之后的约束值，求出cnorm中的最大值即对应的最大约束违背度smax，如果最大约束违背度大于约束容差stol，则根据预先定义的增长因子μ1>1增大罚系数，从而使得下一次优化迭代更倾向于提高可行性；相反，如果最大约束违背度小于或等于容差，则根据初始定义的下降因子μ2≤1减小罚系数，从而使得下一次优化迭代侧重于改善最优性。另外，引入罚系数上下限约束(λmin,λmax)，一方面避免罚系数过小导致优化结果不可行；另一方面，在优化过程中避免罚系数过大所带来数值困难。自适应罚系数计算公式如下所示：

根据作者经验，本研究中μ1取值为2；μ2取值为0.8；λmax取值为50；λmin取值为1；stol取值为0.001。

步骤8，生成伪目标函数值。

提取步骤6样本点数据库中的目标函数响应值和约束值，根据步骤7中得到的罚系数对目标函数值进行更新，更新公式如下：

minf(x)＝f(x)/favg+λmax(c(x)/cavg,0)(38)

式中，f(x)为样本点x处的伪目标函数值，f(x)为样本点x处的真实目标函数值，c(x)为样本点x处的约束值。

步骤8，在步骤7的基础上，令k＝k+1，转入步骤4进行下一次迭代直至达到最大迭代次数200次，最终输出优化结果完成优化过程。

连续十次优化最好结果对应的约束值如表2所示，三种算法所得最好结果对应的设计变量值、目标函数值以及不同算法的平均模型调用次数如表3所示。图4中的箱线图给出了连续十次优化所有结果的对比。由表2可知，迭代终止时ga、cimps和c-ego的最好优化结果均满足约束。

表2对地观测卫星优化问题的约束条件及最优值

表3设计变量取值范围及优化结果对比

由表3和图5(a)可得，ga、cimps和c-ego得到的最优目标函数值较为接近，且在目标函数值的最优性上cimps最好、c-ego次之、ga最差，说明c-ego求解技术正确可行。结合图5(c)可知，尽管cimps能够获取较优的解，但其优化结果有违反约束的风险，而ga和c-ego的优化结果均在约束范围内。由图5(b)可得，在总的mda和约束调用次数上cimps>ga>c-ego，且cimps为c-ego的近2000多，ga为c-ego的20多倍，因而c-ego具有明显的效率优势。综合以上分析可得，尽管在最优性上c-ego稍差于cimps，但是相较于ga和cimps，c-ego在模型调用次数上具有较大的优势，且所得结果均能够满足约束，说明c-ego具有很强的工程实用性。

综合表2和表3的结果可知，初始设计方案不可行，未能满足卫星总重与固有频率约束。相比于初始设计方案，c-ego所得优化方案的中心承力筒直径增加83.3％，本体高度降低31.8％，本体宽度增加49.2％，从而提高了固有频率。此外，其它设计变量也有不同程度的改变，有效降低了整星质量，最终获取了可行的优化设计方案。

综上所得，c-ego的求解效率明显优于ga和cimps，且能够获得较优的可行解。因此，c-ego优化策略能够一定程度上提高卫星总体设计的效率，缩短设计周期，为卫星概念设计提供有力支撑。同样可以看出，本实施例考虑了局部近似精度与全局近似精度，同时具有处理约束的能力，具有全局搜索能力，并且能够获得满足约束的最优设计方案。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例，用于解释本发明，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。