CN108984864A - 基于高斯过程方差估计的高维复杂飞行器模型表征方法 - Google Patents
基于高斯过程方差估计的高维复杂飞行器模型表征方法 Download PDFInfo
- Publication number
- CN108984864A CN108984864A CN201810683292.1A CN201810683292A CN108984864A CN 108984864 A CN108984864 A CN 108984864A CN 201810683292 A CN201810683292 A CN 201810683292A CN 108984864 A CN108984864 A CN 108984864A
- Authority
- CN
- China
- Prior art keywords
- order
- constituent element
- point
- construction
- dimension
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 91
- 230000008569 process Effects 0.000 title claims abstract description 39
- 238000011156 evaluation Methods 0.000 title claims abstract description 21
- 239000000470 constituent Substances 0.000 claims abstract description 63
- 238000010276 construction Methods 0.000 claims abstract description 61
- 238000013461 design Methods 0.000 claims abstract description 45
- 238000005457 optimization Methods 0.000 claims abstract description 36
- 238000005070 sampling Methods 0.000 claims description 30
- 238000001514 detection method Methods 0.000 claims description 11
- 230000004044 response Effects 0.000 claims description 11
- 238000004458 analytical method Methods 0.000 claims description 10
- 238000012360 testing method Methods 0.000 claims description 9
- 238000005520 cutting process Methods 0.000 claims description 6
- 208000035051 Malignant migrating focal seizures of infancy Diseases 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000012217 deletion Methods 0.000 claims description 4
- 230000037430 deletion Effects 0.000 claims description 4
- 208000012054 malignant migrating partial seizures of infancy Diseases 0.000 claims description 4
- 238000010606 normalization Methods 0.000 claims description 4
- 239000012530 fluid Substances 0.000 claims description 3
- 239000000203 mixture Substances 0.000 claims description 3
- 230000001174 ascending effect Effects 0.000 claims description 2
- 238000012512 characterization method Methods 0.000 claims description 2
- 239000012141 concentrate Substances 0.000 claims description 2
- 230000008878 coupling Effects 0.000 claims description 2
- 238000010168 coupling process Methods 0.000 claims description 2
- 238000005859 coupling reaction Methods 0.000 claims description 2
- 230000011218 segmentation Effects 0.000 claims description 2
- 238000012795 verification Methods 0.000 claims 1
- 239000003795 chemical substances by application Substances 0.000 description 33
- 230000008901 benefit Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 238000004088 simulation Methods 0.000 description 6
- 238000004422 calculation algorithm Methods 0.000 description 4
- 238000004891 communication Methods 0.000 description 4
- 230000000052 comparative effect Effects 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 3
- 238000013178 mathematical model Methods 0.000 description 3
- 238000005314 correlation function Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000009826 distribution Methods 0.000 description 2
- 230000002068 genetic effect Effects 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 238000005192 partition Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- CNKHSLKYRMDDNQ-UHFFFAOYSA-N halofenozide Chemical compound C=1C=CC=CC=1C(=O)N(C(C)(C)C)NC(=O)C1=CC=C(Cl)C=C1 CNKHSLKYRMDDNQ-UHFFFAOYSA-N 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 210000001503 joint Anatomy 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000012916 structural analysis Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- 229910052724 xenon Inorganic materials 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/15—Vehicle, aircraft or watercraft design
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- General Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Computer Hardware Design (AREA)
- Automation & Control Theory (AREA)
- Aviation & Aerospace Engineering (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开的基于高斯过程方差估计的高维复杂飞行器模型表征方法，属于高维复杂飞行器工程优化技术领域。本发明利用Kriging逐个近似Cut‑HDMR的各组元函数，在各组元函数Kriging代理模型构造过程中，利用已构造的Kriging模型提供的高斯过程预测方差引导新增样本点的采集，加速组元Kriging代理模型构造的收敛，实现对高维复杂飞行器系统的近似建模表征。本发明通过Cut‑HDMR的全局近似框架，将原始高维复杂问题分解为多个低维问题求和的形式，完成各低维组元函数的近似建模后，加和即获得高维复杂问题的全局近似模型。本发明能有效缓解高维复杂飞行器系统设计优化面临的计算复杂性问题，对于高维复杂问题的全局近似建模、高效近似优化具有重要意义。
Description
技术领域
本发明涉及一种基于高斯过程方差估计信息的高维复杂飞行器系统模型表征方法，属于工程优化技术领域。
背景技术
随着建模仿真技术的发展和计算机技术的进步，为了充分挖掘产品设计潜力、提高设计质量与可信度，高精度数值仿真模型在工程优化中得到了越来越广泛的应用。然而高精度分析模型在提高分析精度和可信度的同时，也带来了计算耗时的困难。高维复杂飞行器系统设计优化问题往往涉及多个相互耦合的学科(如气动学科、结构学科等)，属于多学科耦合工程系统，涉及高维设计变量与大量约束条件，直接使用高精度仿真模型用于高维复杂飞行器系统设计优化的计算代价非常巨大，甚至难以接受。为了减轻基于数值仿真模型的高维复杂飞行器系统设计优化问题的计算代价，代理模型方法得到了广泛的关注与应用。代理模型方法通过调用少量计算耗时的飞行器数值仿真模型获得训练样本点，利用数学手段对数值仿真模型构造精度满足要求的近似模型，以极大的提高飞行器系统的设计优化效率。
传统的代理模型方法能很好地解决中低维近似优化问题，但当系统的维度升高时，构造代理模型的计算成本呈指数级增长，即所谓的“维度灾难”。高维模型表征方法(High Dimensional Model Representation,HDMR)是一种专门用于高维近似建模表征的全局代理模型方法(Global Metamodeling)，其旨在整个设计空间内构造一个全局近似精度较高的代理模型，然后基于该全局代理模型在设计空间内寻优。该方法能直接揭示变量之间的关联并具有完整的数学表达式，拥有较广阔的应用前景。
Cut-HDMR是HDMR的一种主要变形形式，其在输入空间中选取一个参考点作为切割中心，以各维度的方向向量为基，过切割中心的直线、平面以及超平面称为中心基，通过计算中心基上的各阶组元函数，叠加求和得到最终的Cut-HDMR模型。且对于一般的问题，Cut-HDMR仅扩展至二阶组元函数。可将传统代理模型方法与Cut-HDMR结合，对Cut-HDMR各一阶、二阶组元函数进行近似。
在传统代理模型方法中，Kriging，又称高斯过程回归(Gaussian ProcessRegression，GPR)，其近似能力最强且能提供任意未知设计点处的预测方差以表征各点的近似误差。为充分发挥Cut-HDMR的全局近似优势，本发明将Kriging与Cut-HDMR结合，即利用Kriging近似各低阶组元函数。同时在构造组元Kriging代理模型时，采用一种基于Kriging模型高斯过程方差估计信息的自适应序列采样方法，加速组元Kriging构造的收敛，以进一步提高全局近似精度与近似建模效率。
为了更好地说明本发明的技术方案，下面对涉及到的Cut-HDMR方法和Kriging模型的数学基础进行简要介绍。
Cut-HDMR方法：
HDMR有多种展开方式，其中Cut-HDMR理论基础相对完善且易于理解。其在输入空间中选取一个参考点x0＝(c1,c2,...,cd)作为切割中心，以各维度的方向向量为基，过切割中心的直线、平面以及超平面称为中心基，通过计算中心基上的各阶组元函数，叠加求和表示f(x)。Cut-HDMR各阶组元项的表达式如下所示：
式中，f0是原模型在切割中心的响应值。表示点(c1,c2,...,ci-1,xi,ci+1,...,cd)，即向量中除第i维坐标xi外，其余维度的坐标值均与切割中心对应的坐标值相同。依此类推，可知更高阶项的计算方式。每个一阶组元项fi(xi)均根据相应的中心基在xi的直线方向求得；每个二阶组元项fij(xi,xj)根据相应的中心基，在xi和xj组成的二维平面方向上求得。
Kriging模型：
Kriging(KRG)模型是一种针对空间分布数据的无偏最优估计插值模型，由全局模型和局部偏差叠加而成，如下式所示
其中，g(x)为多项式全局近似模型，反映近似对象在设计空间内的总体变化趋势，当近似对象数值特征未知时，可取常数μ。局部偏差项Z(x)是均值为零、方差为σ2、协方差非零的随机过程。KRG的近似能力主要由局部偏差项Z(x)决定。Z(x)的协方差矩阵可表示为
Cov[Z(xi),Z(xj)]＝σ2R[R(xi,xj)] (6)
其中，R为高斯相关函数，R为对称相关矩阵。
任意点x的相关矢量r(x)为
从而，式(5)中KRG数学模型可改写为
值得强调的是，KRG可通过式(10)预估任意点x处的方差s2(x)，可引导新样本点的选取，以加速代理模型构造的收敛。
发明内容
传统的代理模型方法在处理高维复杂飞行器问题时仍存在缺陷，即构造代理模型的计算成本随问题维度呈指数级增长。针对该缺陷，本发明公开的基于高斯过程方差估计的高维复杂飞行器模型表征方法(Kriging assisted High Dimensional ModelRepresentation method with Variance based Sampling,KRG-HDMR-VS)要解决的技术问题为：实现对高维复杂飞行器系统的高效高精度近似建模，缓解高维复杂飞行器系统设计优化面临的计算复杂性问题。
本发明的目的是通过下述技术方案实现的。
本发明公开的基于高斯过程方差估计的高维复杂飞行器模型表征方法，将Kriging用于Cut-HDMR展开式中各组元项的近似建模，并充分利用Kriging提供的预测方差信息引导组元项代理模型构造样本点的序列采集。即利用Kriging逐个近似Cut-HDMR的各组元函数，在各组元函数Kriging代理模型构造过程中，利用已构造的Kriging模型提供的高斯过程预测方差引导新增样本点的采集，加速组元Kriging代理模型构造的收敛，实现对高维复杂飞行器系统的近似建模表征。本发明能够克服传统代理模型方法在处理高维复杂问题时面临的维度灾难问题，通过Cut-HDMR的全局近似框架，将原始高维复杂问题分解为多个低维问题求和的形式，完成各低维组元函数的近似建模后，加和即获得高维复杂问题的全局近似模型，能有效缓解高维复杂飞行器系统设计优化面临的计算复杂性问题，对于高维复杂问题的全局近似建模、高效近似优化具有重要意义。
本发明公开的基于高斯过程方差估计的高维复杂飞行器模型表征方法，包括如下步骤：
步骤A：将各维设计变量归一化到[-1,1]区间，后续均利用归一化后的设计变量值构造各组元Kriging代理模型。
步骤B：选取设计空间的中心点作为切割中心点，并调用真实函数f(x)得到响应值f0。
步骤C：逐维进行一阶组元项Kriging代理模型构造，在对应维切割线上采点，直至构造的Kriging代理模型满足收敛条件。先使用设计变量的边界值作为构造样本点构造初始组元Kriging代理模型，分别调用真实函数得到响应值和然后用中心点x0检测fi(xi)的线性，若构造的初始组元Kriging代理模型通过中心点，则fi(xi)视作线性项，构造终止；否则，将中心点也加入到构造样本点集中，并采用最大均方差(Maximum Mean Square Error,MMSE)采点策略进行序列有偏采样直至组元Kriging代理模型构造收敛。一阶项收敛准则为相对误差准则及最大构造样本点规模准则，当或已有一阶项构造样本点数目大于给定阈值N1时，构造收敛。
最大均方差MMSE采点策略的实现方法为：在构造各组元Kriging代理模型过程中，通过式(1)预估各未知设计点处的高斯过程估计方差s2(xi)。
进而在当前高斯过程估计方差最大处新增样本点，即
式中：为新增样本点，s2(xi)为各未知设计点处的高斯过程估计方差。
步骤D：在完成所有一阶组元项构造后，判断一阶KRG-HDMR-VS模型是否满足精度要求，若足够精确，则二阶组元项影响可忽略不计，构造过程终止；否则继续二阶以及更高阶组元项构造。从每个一阶组元项构造样本点集中选取一个样本组成一个新的检测点xvalidate＝[x1,x2,L,xi,L,xd]T，在该点调用真实函数得到响应值f(xvalidate)，与对应的一阶KRG-HDMR-VS模型预测值相比较，如果则不存在更高阶项，构造过程终止；否则进入下一步，继续二阶组元项的构造。
步骤E：逐个判断各维变量组成的切割平面xi-xj上的二阶组元项存在与否，若存在，则利用对应切割线上的一阶构造样本点集在二维平面上形成的网格点阵Pij作为二阶备选点集，如式(3)所示。
然后从Pij中选取一点进一步检测的近似精度，如果通过该点，说明xi和xj没有耦合关系，并进行后续二阶组元项的构造；否则继续从Pij采点，构造该二阶组元项直至收敛。构造该二阶组元项具体实现方法包括如下步骤：
步骤E1：分别选取各维变量的下界值作为第一个检测样本点并从备选点集Pij中删除该点。调用真实函数得到响应值f(xdetect)，如果说明xi和xj没有耦合关系；否则进入下一步，构造该二阶组元项Kriging代理模型。
步骤E2：使用c-means聚类方法将备选点集Pij分为K个聚类Ci(i＝1,2,L,K)。针对二阶组元项在二维平面xi-xj构造的特点，将聚类数取为K＝4。依次从各聚类中选取一个点记作一代，即每一代从归一化后的xi-xj平面的四个象限中分别取一个点，以提升构造样本点的空间均布性。在各聚类中取点时，根据MPEI采样准则选取预测误差指标最大的点，并从相应聚类中将该点删除。
步骤E3：利用第一个检测样本点xdetect和已有的一阶构造样本点和及中心点生成初始二阶构造样本点集Tij，以构造二阶组元Kriging代理模型同时利用最大预测误差指标(Maximum Prediction Error Index,MPEI)采样准则重新获取点作为新的检测样本点，并从相应聚类中将其删除。重复该过程直至二阶组元项构造收敛。二阶组元项收敛准则为相对误差准则及最大采样代数准则，当或采样代数大于给定阈值G1时，构造收敛。
利用最大预测误差指标采样准则重新获取点作为新的检测样本点，具体实现方法为：在进行二阶组元项构造时已经完成所有一阶组元项的构造，一阶组元Kriging的预测方差也能用于引导二阶构造样本点的采集。二阶组元项方差如下式：
式中Vi为一阶方差项，表征了输入变量xi对总方差V(f(x))的影响。基于已构造的一阶组元Kriging代理模型其值可用预测方差s2(xi)表示。二阶方差项表征了输入变量xi和xj对总方差V(f(x))的耦合影响。则式(14)改写为：
式中G(xi,xj)为未知函数；(s2(xi)+s2(xj))代表一阶预测总方差，记作当较大时，二阶预测方差Vij较小，即当一阶预测总方差占比高时，相应二阶方差项占比小。因而能被直接用于引导二阶序列有偏采样。优先选取较大的点，以降低二阶组元项不确定性进而提升预测精度。
当计算的s2(xi)和s2(xj)量级不一致时，难以体现较小预测方差一项对的贡献。确定二阶备选点的采样顺序主要根据的排序结果，而不关心其具体取值。为此采用预测误差指标(Prediction Error Index，PEI)为每个二阶备选点标记一个PEI，进而消除一阶预测方差量级不一致对排序结果的影响。预测误差指标PEI计算方式如式(6)所示。
PEI(xi,xj)＝Scorei(s2(xi))+Scorej(s2(xj)) (16)
式中Scorei和Scorej分别表示一阶预测方差的分数。其取值根据和在各点处的预测方差值s2(xi)和s2(xj)分别升序排列，然后对应赋整数值{1,2,L,n}作为一阶预测方差的分数。显然高分数值Scorei代表了在此处的预测误差较大，反之则越小。与一阶组元项采用的MMSE方法类似，此处优先取出预测误差指标最大的样本点，表示为：
式中，为获取的二阶检测样本点。
步骤F：构造完所有二阶组元项后，输出二阶KRG-HDMR-VS近似模型即实现基于高斯过程方差估计的高维复杂飞行器模型表征。
步骤G：校验步骤F输出的二阶KRG-HDMR-VS近似模型的全局近似精度。采用分割样本点验证，即额外采集一定数量测试样本点，计算全局精度评价指标R2，RAAE和RMAE。
步骤H：当步骤F获取的近似模型满足步骤G校验精度要求时，将输出二阶KRG-HDMR-VS近似模型应用于高维复杂飞行器系统的近似优化领域，解决高维复杂飞行器系统领域相关工程问题。
所述高维复杂飞行器系统领域相关工程问题包括：含有高维设计变量的基于有限元分析的飞行器结构优化设计、基于计算流体力学分析的飞行器气动优化设计，能有效缓解高维复杂飞行器系统设计优化面临的计算复杂性问题，进而解决高维复杂飞行器系统领域相关工程问题。
有益效果：
1、本发明公开的基于高斯过程方差估计的高维复杂飞行器模型表征方法，通过将Kriging与典型的高维模型表征方法展开式Cut-HDMR相结合，充分利用Kriging对低维问题良好的近似能力及Cut-HDMR的全局近似框架，使用Kriging提供的高斯过程预测方差加速组元Kriging代理模型构造的收敛，该方法是一种新型的全局代理模型方法，对于高维复杂问题的全局近似建模、高效近似优化具有重要意义。
2、本发明公开的基于高斯过程方差估计的高维复杂飞行器模型表征方法，通过Cut-HDMR的全局近似框架，将原始高维复杂问题分解为多个低维问题求和的形式，完成各低维组元函数的近似建模后，加和即可获得高维复杂问题的全局近似模型，适用于需调用复杂耗时的数值仿真模型的飞行器系统设计优化问题，如含有高维设计变量的基于有限元分析的飞行器结构优化设计、基于计算流体力学分析的飞行器气动优化设计等，能有效缓解高维复杂飞行器系统设计优化面临的计算复杂性问题。
附图说明
图1为基于高斯过程方差估计的高维模型表征方法(KRG-HDMR-VS)流程图；
图2为2维问题的一阶项f1(x1)MMSE采点示意图，其中：(a)为第一次采点，(b)为第二次采点；
图3为2维问题的一阶项f2(x2)MMSE采点示意图，其中：(a)为第一次采点，(b)为第二次采点；
图4为2维问题的二阶构造样本点示意图，其中：(a)为二阶备选点集，(b)为二阶备选点聚类；
图5为优化前后翼型几何外形及压力系数分布对比，其中：(a)为几何外形，(b)为压力系数分布；
图6为卫星结构示意图，其中：1—顶板、2—通信舱南板、3—通信舱东板、4—服务舱东板、5—服务舱南板、6—中板、7—氙气瓶、8—对接环、9—底板、10—服务舱西板、11—服务舱北板、12—通信舱西板、13—服务舱隔板、14—中心承力筒、15—通信舱隔板、16—通信舱北板
图7为整星有限元模型；
图8为有限元模型的模态振型结果。
具体实施方式
为了更好地说明本发明的目的和优点，下面将结合附图和具体实例对本发明做进一步的说明，并通过与国际同类算法基于改进径向基函数的高维模型表征方法(AERBF-HDMR)进行对比，对本发明的综合性能进行验证分析。
实施例1：
下面采用一个二维问题，说明具体实施过程。
选取的二维算例为：
一阶项给定的收敛准则为ε1＝0.001、一阶项最大构造样本点规模N1＝4。二阶项给定的收敛准则为ε2＝0.001、最大采样代数G1＝1。
本发明公开的基于高斯过程方差估计的高维复杂飞行器模型表征方法，具体实现步骤如下：
步骤A-B：选择中心点作为切割中心，因各维设计变量均归一化到区间[-1,1]，则中心点为坐标系原点。不失一般性，一般在中心点基础上对各维变量加上一个较小的随机扰动，则x0＝[0.0496,0.0926]，f0＝f(x0)＝-9.4356。
步骤C：构造一阶项f1(x1)，首先选用边界值xiL＝-1，xiR＝1，分别调用真实函数得到响应值 利用两个边界点先构造初始组元Kriging代理模型，然后用中心点x0＝0.0496检测fi(xi)的线性。此时 显然fi(xi)是非线性项，采用MMSE采点策略进行序列有偏采样直至收敛。将中心点x0加入到构造样本点集得到{-1,1,0.0496}，重新构造组元Kriging代理模型，得到各未知设计点处的预测方差如图2(a)所示。根据MMSE准则在最大预测方差处新增点，即得到 将该点也加入到构造样本点集得到{-1,1,0.0496,-0.4949}，重新构造组元Kriging代理模型，得到预测方差如图2(b)所示，继续利用MMSE准则得到新增点得到 将加入到构造样本点集得到{-1,1,0.0496,-0.4949,0.6566}，重新构造组元Kriging代理模型，此时一阶构造样本点数目达到给定阈值N1，构造收敛，停止f1(x1)构造。
继续构造一阶项f2(x2)，与f1(x1)构造过程类似，利用MMSE采点的示意图如图3所示，获取的组元构造样本点集为{-1,1,0.0926,-0.4747,0.6364}，对应的组元响应值依次为{181.0475,31.1223,0,64.2117,-1.1227}。
步骤D：构造完所有一阶组元项后，判断一阶KRG-HDMR-VS模型是否满足精度要求，取每维变量下界值形成检测点xvalidate＝[-1,-1]，得到f(xvalidate)＝190，显然存在更高阶项，进入下一步，继续二阶项构造。
步骤E：判断x1与x2是否存在耦合关系，得到新的检测样本点xdetect＝[-1,-1]，得到fij(xdetect)＝-114.6837，f(xdetect)＝190，显然该二阶项存在。利用一阶样本点生成二阶备选点阵Pij如图4(a)所示。使用c-means聚类方法将Pij分为4个聚类如图4(b)所示，不同颜色标记的点分属不同聚类。依次从各个聚类中利用MPEI采样准则新增样本点，直至达到给定的最大采样代数G1，构造收敛，停止构造f12(x1,x2)。所采集的二阶样本点依次为{(-1,-1)；(1,-0.4747)；(1,0.6364)；(-0.4949,0.6364)}
步骤F：构造完二阶项f12(x1,x2)，即可输出二阶KRG-HDMR-PS近似模型所用构造样本点数量共为13。
步骤G：校验输出的二阶KRG-HDMR-VS近似模型采用分割样本点验证，额外采集500个测试样本点进行精度校验，计算全局近似精度评价指标，得到R2＝1.0000，RAAE＝0.0001和RMAE＝0.0003。
为了更好的说明本发明的目的和优势，选取标准数值算例和工程算例对本发明公开的基于高斯过程方差估计的高维模型表征方法(KRG-HDMR-VS)进行测试，并与国际同类算法AERBF-HDMR进行性能对比。
(一)标准测试算例
标准测试算例的表达式如表1所示。
同样采用R2，RAAE和RMAE三项评价指标进行精度校验，R2越接近1、RAAE和RMAE越接近0，表明全局近似精度越高。选择分割样本点验证作为校验方法，随机采集1000个样本点作为测试样本点。采用KRG-HDMR-VS和AERBF-HDMR对上述标准测试函数连续近似10次的结果统计均值如表2所示。近似成本即所需的模型调用次数(Number of Evaluation,NoE)为表2中最后一列所示，表征近似方法的近似效率。
显然除算例2外，KRG-HDMR-VS的精度均高于AERBF-HDMR。同时在近似成本上，KRG-HDMR-VS仅在算例3上略高于AERBF-HDMR。对于算例3，KRG-HDMR-VS一阶构造样本点数多于AERBF-HDMR是因为KRG-HDMR一般至少需要4个样本点以构造一阶组元项Kriging代理模型。
表1 测试函数表达式
表2 KRG-HDMR-VS与AERBF-HDMR近似建模结果对比
(二)工程算例
以一个翼型气动优化设计和卫星结构优化设计为例，介绍KRG-HDMR-VS方法在高维复杂飞行器设计优化中的应用。
翼型气动优化设计：
选取NACA64A816翼型作为初始2维翼型，通过形函数扰动法对翼型进行参数化建模，上下翼面分别选取5个权系数作为设计变量，即xui,xli(i＝1,2,3,4,5)，共10个设计变量。翼型气动优化问题的数学模型如下：
式中，-CL/D为负的升阻比，tmax表示翼型的最大厚度，为基准翼型的最大厚度，x是设计变量，xlb和xup分别是设计空间的下界和上界，各维设计变量的取值不超过±0.005。建立翼型参数化模型后，需调用计算耗时的计算流体力学(CFD)模型获取翼型升阻比。
针对马赫数为0.64，攻角为2°时的飞行工况，利用Kriging和KRG-HDMR-VS分别进行10次翼型的气动优化，最终结果如表3所示。
表3 KRG-HDMR-VS近似建模结果
由表3可知，在近似成本相同时，KRG-HDMR-VS的近似精度明显高于Kriging代理模型。
表4 翼型优化前后结果对比
近似建模结束后，使用遗传算法基于KRG-HDMR-VS模型进行优化，相较于直接基于真实的气动分析模型进行优化，显然基于KRG-HDMR-VS模型的优化能节约大量的计算成本。翼型优化前后结果对比如表4所示，与基准翼型相比，翼型升阻比提高了23.5％，同时翼型最大厚度满足约束。将优化后翼型代入CFD分析模型以校验KRG-HDMR-VS模型的精度，发现KRG-HDMR-VS的预测值与CFD值非常接近，进一步验证了KRG-HDMR-VS方法的准确性和有效性。图5为优化前后翼型几何外形及压力系数分布对比。
卫星结构优化设计：
所采用的案例卫星结构示意图如图6所示，其星体呈正六面体，采用中心承力筒构型。整理的结构学科设计变量共16个，分别为各舱板及承力筒的芯子厚度和铺层厚度。优化目标为整星质量最小，约束为X、Y、Z方向的一阶基频约束。卫星结构优化问题的数学模型如式所示。
式中，msatellite为整星质量，fX、fY分别为绕X、Y轴的横向一阶基频，fZ为沿Z轴的纵向一阶基频。该结构分析模型需通过Patran建立整星有限元模型并调用Nastran进行结构分析，所建立的整星有限元模型如图7所示。
采用KRG-HDMR-VS对目标函数msatellite进行近似建模，然而该优化问题存在三个复杂约束，因而利用KRG-HDMR-VS分别针对目标函数msatellite和三个复杂约束fX、fY、fZ构造近似模型。近似完成后采集500个随机样本点作为测试样本点，近似结果如表5所示。为对比KRG-HDMR-VS和Kriging的近似能力，采用相同数量的样本点构造目标函数msatellite的Kriging代理模型，Kriging的设置仍为零阶多项式回归模型和高斯相关函数，与KRG-HDMR-VS组元Kriging代理模型的设置保持一致，其近似结果如表5所示。显然，在相同构造样本点规模及相同Kriging设置下，KRG-HDMR-VS的近似能力显著优于Kriging。
表5 KRG-HDMR-VS与Kriging近似结果对比
采用遗传算法基于目标函数和约束条件的KRG-HDMR-VS代理模型进行优化，优化结果如表6所示。在满足三个基频约束的前提下，整星质量降低了44.8kg。用有限元模型校验最优解处的响应值，其与KRG-HDMR-VS的预测值非常接近，进一步验证了KRG-HDMR-VS方法的有效性。有限元模型得到的模态振型如图8所示。
表6 卫星优化前后结果对比
通过上述两类优化测试问题对比，本发明公开的一种基于高斯过程方差估计信息的高维模型表征方法(KRG-HDMR-VS)能充分利用Kriging对低维问题的近似优势以及Kriging提供的高斯过程预测方差信息，提升其对高维复杂问题的全局近似精度与近似效率。数值测试算例结果表明，相较于国际同类算法中性能最优的AERBF-HDMR，本发明公布的KRG-HDMR-VS方法在近似精度与近似效率上均具有明显优势。翼型气动优化及卫星结构优化算例结果表明，KRG-HDMR-VS能有效地实现对高维复杂工程问题的全局近似，从而实现近似优化，提升优化效率。
以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例，用于解释本发明，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (4)
1.基于高斯过程方差估计的高维复杂飞行器模型表征方法，其特征在于：包括如下步骤，
步骤A：将各维设计变量归一化到[-1,1]区间，后续均利用归一化后的设计变量值构造各组元Kriging代理模型；
步骤B：选取设计空间的中心点作为切割中心点，并调用真实函数f(x)得到响应值f0；
步骤C：逐维进行一阶组元项Kriging代理模型构造，在对应维切割线上采点，直至构造的Kriging代理模型满足收敛条件；先使用设计变量的边界值作为构造样本点构造初始组元Kriging代理模型，分别调用真实函数得到响应值和然后用中心点x0检测fi(xi)的线性，若构造的初始组元Kriging代理模型通过中心点，则fi(xi)视作线性项，构造终止；否则，将中心点也加入到构造样本点集中，并采用最大均方差(Maximum Mean Square Error,MMSE)采点策略进行序列有偏采样直至组元Kriging代理模型构造收敛；一阶项收敛准则为相对误差准则及最大构造样本点规模准则，当或已有一阶项构造样本点数目大于给定阈值N1时，构造收敛；
最大均方差MMSE采点策略的实现方法为：在构造各组元Kriging代理模型过程中，通过式(1)预估各未知设计点处的高斯过程估计方差s2(xi)；
进而在当前高斯过程估计方差最大处新增样本点，即
式中：为新增样本点，s2(xi)为各未知设计点处的高斯过程估计方差；
步骤D：在完成所有一阶组元项构造后，判断一阶KRG-HDMR-VS模型是否满足精度要求，若足够精确，则二阶组元项影响可忽略不计，构造过程终止；否则继续二阶以及更高阶组元项构造；从每个一阶组元项构造样本点集中选取一个样本组成一个新的检测点xvalidate＝[x1,x2,L,xi,L,xd]T，在该点调用真实函数得到响应值f(xvalidate)，与对应的一阶KRG-HDMR-VS模型预测值相比较，如果则不存在更高阶项，构造过程终止；否则进入下一步，继续二阶组元项的构造；
步骤E：逐个判断各维变量组成的切割平面xi-xj上的二阶组元项存在与否，若存在，则利用对应切割线上的一阶构造样本点集在二维平面上形成的网格点阵Pij作为二阶备选点集，如式(3)所示；
然后从Pij中选取一点进一步检测的近似精度，如果通过该点，说明xi和xj没有耦合关系，并进行后续二阶组元项的构造；否则继续从Pij采点，构造该二阶组元项直至收敛；
步骤F：构造完所有二阶组元项后，输出二阶KRG-HDMR-VS近似模型即实现基于高斯过程方差估计的高维复杂飞行器模型表征。
2.如权利要求1所述的基于高斯过程方差估计的高维复杂飞行器模型表征方法，其特征在于：还包括如下步骤，
步骤G：校验步骤F输出的二阶KRG-HDMR-VS近似模型的全局近似精度；采用分割样本点验证，即额外采集一定数量测试样本点，计算全局精度评价指标R2，RAAE和RMAE；
步骤H：当步骤F获取的近似模型满足步骤G校验精度要求时，将输出二阶KRG-HDMR-VS近似模型应用于高维复杂飞行器系统的近似优化领域，解决高维复杂飞行器系统领域相关工程问题。
3.如权利要求1或2所述的基于高斯过程方差估计的高维复杂飞行器模型表征方法，其特征在于：构造该二阶组元项具体实现方法包括如下步骤，
步骤E1：分别选取各维变量的下界值作为第一个检测样本点并从备选点集Pij中删除该点；调用真实函数得到响应值f(xdetect)，如果说明xi和xj没有耦合关系；否则进入下一步，构造该二阶组元项Kriging代理模型；
步骤E2：使用c-means聚类方法将备选点集Pij分为K个聚类Ci(i＝1,2,L,K)；针对二阶组元项在二维平面xi-xj构造的特点，将聚类数取为K＝4；依次从各聚类中选取一个点记作一代，即每一代从归一化后的xi-xj平面的四个象限中分别取一个点，以提升构造样本点的空间均布性；在各聚类中取点时，根据MPEI采样准则选取预测误差指标最大的点，并从相应聚类中将该点删除；
步骤E3：利用第一个检测样本点xdetect和已有的一阶构造样本点和及中心点生成初始二阶构造样本点集Tij，以构造二阶组元Kriging代理模型同时利用最大预测误差指标(Maximum Prediction Error Index,MPEI)采样准则重新获取点作为新的检测样本点，并从相应聚类中将其删除；重复该过程直至二阶组元项构造收敛；二阶组元项收敛准则为相对误差准则及最大采样代数准则，当或采样代数大于给定阈值G1时，构造收敛；
利用最大预测误差指标采样准则重新获取点作为新的检测样本点，具体实现方法为：在进行二阶组元项构造时已经完成所有一阶组元项的构造，一阶组元Kriging的预测方差也能用于引导二阶构造样本点的采集；二阶组元项方差如下式：
式中Vi为一阶方差项，表征了输入变量xi对总方差V(f(x))的影响；基于已构造的一阶组元Kriging代理模型其值可用预测方差s2(xi)表示；二阶方差项表征了输入变量xi和xj对总方差V(f(x))的耦合影响；则式(4)改写为：
式中G(xi,xj)为未知函数；(s2(xi)+s2(xj))代表一阶预测总方差，记作当较大时，二阶预测方差Vij较小，即当一阶预测总方差占比高时，相应二阶方差项占比小；因而能被直接用于引导二阶序列有偏采样；优先选取较大的点，以降低二阶组元项不确定性进而提升预测精度；
当计算的s2(xi)和s2(xj)量级不一致时，难以体现较小预测方差一项对的贡献；确定二阶备选点的采样顺序主要根据的排序结果，而不关心其具体取值；为此采用预测误差指标(Prediction Error Index，PEI)为每个二阶备选点标记一个PEI，进而消除一阶预测方差量级不一致对排序结果的影响；预测误差指标PEI计算方式如式(6)所示；
PEI(xi,xj)＝Scorei(s2(xi))+Scorej(s2(xj)) (6)
式中Scorei和Scorej分别表示一阶预测方差的分数；其取值根据和在各点处的预测方差值s2(xi)和s2(xj)分别升序排列，然后对应赋整数值{1,2,L,n}作为一阶预测方差的分数；显然高分数值Scorei代表了在此处的预测误差较大，反之则越小；与一阶组元项采用的MMSE方法类似，此处优先取出预测误差指标最大的样本点，表示为：
式中，为获取的二阶检测样本点。
4.如权利要求1或2所述的基于高斯过程方差估计的高维复杂飞行器模型表征方法，其特征在于：所述高维复杂飞行器系统领域相关工程问题包括，含有高维设计变量的基于有限元分析的飞行器结构优化设计、基于计算流体力学分析的飞行器气动优化设计，能有效缓解高维复杂飞行器系统设计优化面临的计算复杂性问题，进而解决高维复杂飞行器系统领域相关工程问题。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810683292.1A CN108984864A (zh) | 2018-06-28 | 2018-06-28 | 基于高斯过程方差估计的高维复杂飞行器模型表征方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810683292.1A CN108984864A (zh) | 2018-06-28 | 2018-06-28 | 基于高斯过程方差估计的高维复杂飞行器模型表征方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN108984864A true CN108984864A (zh) | 2018-12-11 |
Family
ID=64538655
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810683292.1A Pending CN108984864A (zh) | 2018-06-28 | 2018-06-28 | 基于高斯过程方差估计的高维复杂飞行器模型表征方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108984864A (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113657045A (zh) * | 2021-08-10 | 2021-11-16 | 北京理工大学 | 基于多层协同高斯过程的复杂飞行器模型降阶表征方法 |
| CN113935256A (zh) * | 2021-09-27 | 2022-01-14 | 北京理工大学 | 基于误差修正的高维复杂飞行器系统降阶表征方法 |
| CN114169192A (zh) * | 2021-11-23 | 2022-03-11 | 北京航空航天大学 | 基于粗糙集理论的热力耦合系统不确定性分析方法 |
| CN116341097A (zh) * | 2022-02-07 | 2023-06-27 | 西北工业大学 | 一种基于新型高维代理模型的跨音速机翼优化设计方法 |
-
2018
- 2018-06-28 CN CN201810683292.1A patent/CN108984864A/zh active Pending
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113657045A (zh) * | 2021-08-10 | 2021-11-16 | 北京理工大学 | 基于多层协同高斯过程的复杂飞行器模型降阶表征方法 |
| CN113657045B (zh) * | 2021-08-10 | 2024-03-08 | 北京理工大学 | 基于多层协同高斯过程的复杂飞行器模型降阶表征方法 |
| CN113935256A (zh) * | 2021-09-27 | 2022-01-14 | 北京理工大学 | 基于误差修正的高维复杂飞行器系统降阶表征方法 |
| CN114169192A (zh) * | 2021-11-23 | 2022-03-11 | 北京航空航天大学 | 基于粗糙集理论的热力耦合系统不确定性分析方法 |
| CN116341097A (zh) * | 2022-02-07 | 2023-06-27 | 西北工业大学 | 一种基于新型高维代理模型的跨音速机翼优化设计方法 |
| CN116341097B (zh) * | 2022-02-07 | 2023-09-12 | 西北工业大学 | 一种基于新型高维代理模型的跨音速机翼优化设计方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109685152B (zh) | 一种基于dc-spp-yolo的图像目标检测方法 | |
| CN108984864A (zh) | 基于高斯过程方差估计的高维复杂飞行器模型表征方法 | |
| CN110472846B (zh) | 核电厂热工水力安全分析最佳估算加不确定性方法 | |
| CN109559338B (zh) | 一种基于加权主成分分析法及m估计的三维点云配准方法 | |
| Praveen et al. | Low cost PSO using metamodels and inexact pre-evaluation: Application to aerodynamic shape design | |
| CN112100835B (zh) | 一种适用于复杂流动的高效高精度翼型绕流数值模拟方法 | |
| CN108376408A (zh) | 一种基于曲率特征的三维点云数据快速加权配准方法 | |
| CN113836634B (zh) | 一种大差异性气动数据的深度神经网络建模方法 | |
| CN105718634A (zh) | 一种基于非概率区间分析模型的翼型鲁棒优化设计方法 | |
| CN113777931B (zh) | 结冰翼型气动模型构造方法、装置、设备及介质 | |
| CN114119922A (zh) | 基于点云三维重构的复杂轮廓缺陷校验方法 | |
| CN115495965A (zh) | 复杂航空结构在混合不确定性下时变可靠性的分析方法 | |
| CN111815681A (zh) | 基于深度学习和判别式模型训练的目标跟踪方法及存储器 | |
| CN106407601B (zh) | 一种基于数据挖掘技术的气动特性数据处理方法 | |
| CN115345064A (zh) | 一种融入物理约束的集成深度神经网络气动力建模方法 | |
| CN111027250A (zh) | 一种基于网格变形技术的异形曲面加筋壳建模方法 | |
| CN113657029A (zh) | 异源数据驱动的飞行器高效近似优化方法 | |
| CN115206455B (zh) | 基于深度神经网络的稀土元素组分含量预测方法及系统 | |
| CN104699901A (zh) | 基于差量采样解的GappyPOD翼型反设计方法 | |
| CN116933386A (zh) | 一种基于mcok代理模型的飞行器气动数据融合方法 | |
| Chen et al. | An airfoil mesh quality criterion using deep neural networks | |
| CN109033678A (zh) | 一种基于虚拟样本生成的飞行器近似优化设计方法 | |
| CN110083946A (zh) | 一种基于无约束优化模型的多状态模型修正的方法 | |
| CN115470675A (zh) | 一种重载压力容器裂纹扩展路径的预测方法及系统 | |
| CN114357878A (zh) | 一种水陆两栖飞机典型截面的着水面载荷预测方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20181211 |