CN114200843B - 喉栓式固体姿轨控推力快速预示方法 - Google Patents

喉栓式固体姿轨控推力快速预示方法 Download PDF

Info

Publication number
CN114200843B
CN114200843B CN202210137215.2A CN202210137215A CN114200843B CN 114200843 B CN114200843 B CN 114200843B CN 202210137215 A CN202210137215 A CN 202210137215A CN 114200843 B CN114200843 B CN 114200843B
Authority
CN
China
Prior art keywords
training sample
orbit control
throat
solid attitude
type solid
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.)
Active
Application number
CN202210137215.2A
Other languages
English (en)
Other versions
CN114200843A (zh
Inventor
李国盛
武泽平
杨家伟
文谦
王政涛
王鹏宇
张为华
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202210137215.2A priority Critical patent/CN114200843B/zh
Publication of CN114200843A publication Critical patent/CN114200843A/zh
Application granted granted Critical
Publication of CN114200843B publication Critical patent/CN114200843B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/042—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Computation (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

一种喉栓式固体姿轨控推力快速预示方法，确定考虑沉积效应后的喉栓式固体姿轨控推力物理模型、设计变量、约束条件以及设计空间；通过实验设计方法在设计空间生成训练样本点的输入变量，计算各训练样本点处物理模型的真实输出响应；采用混沌多项式增广径向基方法建立喉栓式固体姿轨控推力的增广径向基函数近似模型；将近似模型中多个形状参数确定问题转化为各向异性缩放系数优化问题；采用K‑折交叉验证算法对各向异性缩放系数优化问题进行求解，选择交叉验证误差最小的各向异性缩放系数并将其带回所述近似模型，得到喉栓式固体姿轨控推力快速预示模型。基于喉栓式固体姿轨控推力快速预示模型即可进行喉栓式固体姿轨控推力的快速预示。

Description

喉栓式固体姿轨控推力快速预示方法
技术领域
本发明属于固体姿轨控发动机推力预示技术领域，特别涉及一种喉栓式固体姿轨控推力快速预示方法。
背景技术
固体姿轨控发动机喷管由于沉积效应影响，各喷管实际等效喉部面积与理想等效喉部面积不同，造成实际推力与预期理想推力出现偏差。快速精确预示出考虑沉积效应后的固体姿轨控发动机喷管推力，对实现发动机推力的精准调控具有重大意义。
目前国内外就固体姿轨控发动机推力预示所开展的研究较少。目前常用的基于神经网络的固体姿轨控发动机推力预示方法，包括：
（1）采用CFD对不同喉栓位移及沉积条件下的发动机喷管进行仿真，得出其实际推力。在不同燃烧室压强下对每个喉栓位移点以及耦合不同沉积情况都进行数值模拟分析。
（2）将神经网络技术用于固体姿轨控发动机推力预示，构造固体姿轨控发动机推力预示模型。在拥有一定规模的数据量后，采用神经网络技术构建推力预示的代理模型以降低计算资源消耗。随着实验数据的增加，发动机实际推力预示的精度会不断提高。常用的神经网络主要有误差反向传播（BP）神经网络和径向基（RBF）神经网络。
目前常用的基于神经网络的固体姿轨控发动机推力预示方法其缺点在于：
（1）误差反向传播（BP）神经网络是全局逼近网络，收敛速度较慢且易陷入局部极值；
（2）径向基神经网络是局部逼近网络，学习速度较快且可避免局部极小问题，但对于线性项预示效果不好。采用径向基神经网络进行系统近似建模，需要对径向基神经网络的形状参数的进行合理选择。
（3）交叉验证是估算形状参数的一种有效方法，现有的交叉验证方法每折的样本信息不均匀会导致构建的子代理模型不能很好地描述模型全局信息，导致训练的结果具有局部性。
因此，寻找一种预示精度高、预示速度快的固体姿轨控发动机推力预示方法具有重要的现实意义。
发明内容
针对固体姿轨控发动机喷管由于沉积效应影响，各喷管实际等效喉部面积与理想等效喉部面积不同，造成实际推力与预期理想推力出现偏差的问题，本发明提出了喉栓式固体姿轨控推力快速预示方法，其解决了现有基于神经网络技术固体姿轨控发动机推力预示方法所存在的预示精度低、预示速度慢以及预示信息具有局部性等问题。
为实现上述技术目的，本发明提出的技术方案为：
一方面，本发明提供一种喉栓式固体姿轨控推力快速预示方法，包括：
确定考虑沉积效应后的喉栓式固体姿轨控推力物理模型、设计变量、约束条件以及设计空间；
通过实验设计方法在设计空间生成训练样本点的输入变量，计算各训练样本点处物理模型的真实输出响应；
融合径向基函数以及混沌多项式建立喉栓式固体姿轨控推力的增广径向基函数近似模型；
对所述近似模型中径向基函数的形状参数进行归一化表达，将多个形状参数确定问题转化为各向异性缩放系数优化问题；
将所述训练样本点自适应分割成K层训练样本子集，分割后的各层训练样本子集中的训练样本点均匀分布于设计空间中，采用K-折交叉验证算法对各向异性缩放系数优化问题进行求解，选择交叉验证误差最小的各向异性缩放系数并将其带回所述近似模型，得到喉栓式固体姿轨控推力快速预示模型；
对于待测设计变量组成的输入向量，基于所述喉栓式固体姿轨控推力快速预示模型即可进行喉栓式固体姿轨控推力的快速预示。
在一实施例中，以固体姿轨控发动机轨控喷管作为物理模型，设计变量包括4个，分别是喷管沉积量、喉栓沉积量、燃烧室压强和喉栓位移。采用优化拉丁超立方实验设计方法生成多组均匀分布于设计空间的4维训练样本点。调用所述考虑沉积效应后的喉栓式固体姿轨控推力物理模型，基于CFD完成数值模拟，对固体姿轨控发动机在不同喉栓位移、喉栓沉积量、喷管沉积量、压强条件下进行仿真，得到各训练样本点处物理模型的真实输出响应，即得到发动机推力精细化计算结果。
另一方面，本发明提供一种计算机设备，包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时实现以下步骤：
确定考虑沉积效应后的喉栓式固体姿轨控推力物理模型、设计变量、约束条件以及设计空间；
通过实验设计方法在设计空间生成训练样本点的输入变量，计算各训练样本点处物理模型的真实输出响应；
融合径向基函数以及混沌多项式建立喉栓式固体姿轨控推力的增广径向基函数近似模型；
对所述近似模型中径向基函数的形状参数进行归一化表达，将多个形状参数确定问题转化为各向异性缩放系数优化问题；
将所述训练样本点自适应分割成K层训练样本子集，分割后的各层训练样本子集中的训练样本点均匀分布于设计空间中，采用K-折交叉验证算法对各向异性缩放系数优化问题进行求解，选择交叉验证误差最小的各向异性缩放系数并将其带回所述近似模型，得到喉栓式固体姿轨控推力代理模型；
对于任意设计变量组成的输入向量，基于所述喉栓式固体姿轨控推力代理模型即可进行喉栓式固体姿轨控推力的快速预示。
再一方面，本发明还提供一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现以下步骤：
确定考虑沉积效应后的喉栓式固体姿轨控推力物理模型、设计变量、约束条件以及设计空间；
通过实验设计方法在设计空间生成训练样本点的输入变量，计算各训练样本点处物理模型的真实输出响应；
融合径向基函数以及混沌多项式建立喉栓式固体姿轨控推力的增广径向基函数近似模型；
对所述近似模型中径向基函数的形状参数进行归一化表达，将多个形状参数确定问题转化为各向异性缩放系数优化问题；
将所述训练样本点自适应分割成K层训练样本子集，分割后的各层训练样本子集中的训练样本点均匀分布于设计空间中，采用K-折交叉验证算法对各向异性缩放系数优化问题进行求解，选择交叉验证误差最小的各向异性缩放系数并将其带回所述近似模型，得到喉栓式固体姿轨控推力代理模型；
对于任意设计变量组成的输入向量，基于所述喉栓式固体姿轨控推力代理模型即可进行喉栓式固体姿轨控推力的快速预示。
与现有技术相比，本发明的优点在于：
本发明提出增广径向基模型，针对增广径向基函数模型中参数确定困难的问题展开研究，提出基于样本局部密度的形状参数表征方法，将多个形状参数确定转化为缩放系数确定的单变量优化问题，实现形状参数优化与训练样本数量解耦，降低了径向基函数形状参数确定的计算复杂度。
本发明提出样本点自适应分割方法，使交叉验证中的辅助模型尽可能多地捕获真实模型在设计空间中的特征信息。
针对缩放系数合理确定，本发明提出通用快速K-折交叉验证算法，提升了增广径向基函数近似建模效率和精度。
本发明采用智能优化算法针对缩放系数进行选择，得到交叉验证误差最小的缩放系数，构建固体姿轨控发动机推力快速预示模型。
本发明大大降低了固体姿轨控推力预示所需要的仿真次数，节约了时间和精力，大幅度提高了计算效率；同时本发明建模精度高、计算速度快，大幅度降低了增广径向基代理模型构建的计算量。
本发明已经在多种不同的固体姿轨控发动机设计中得到了应用，与现有方法相比，人工参与过程少、设计速度快，自动化程度高，设计结果性能可靠。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
图1是本发明一实施例的流程图；
图2是本发明一实施例中计算网格划分与边界设置示意图，其中喷管的出、入口边界条件均由具体情况设定，图中对称轴为对称边界，喉栓壁和喷管壁均采用绝热无滑移壁面边界条件。
图3是本发明一实施例中推力快速预示相对误差图；
图4为本发明一实施例的结构示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚明白，下面将以附图及详细叙述来清楚说明本发明所揭示内容的精神，任何所属技术领域技术人员在了解本发明内容的实施例后，当可由本发明内容所教示的技术，加以改变及修饰，其并不脱离本发明内容的精神与范围。本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。
针对固体姿轨控发动机喷管由于沉积效应影响，各喷管实际等效喉部面积与理想等效喉部面积不同，造成实际推力与预期理想推力出现偏差的问题，本发明提出喉栓式固体姿轨控推力快速预示方法，该方法克服了现有基于神经网络技术固体姿轨控发动机推力预示方法所存在的预示精度低、预示速度慢以及预示信息具有局部性等缺陷，有效解决了已知沉积规律的固体姿轨控发动机推力快速精确预示问题。
参照图1，本发明一实施例中，提供一种喉栓式固体姿轨控推力快速预示方法，包括：
(S1) 确定考虑沉积效应后的喉栓式固体姿轨控推力物理模型、设计变量、约束条件以及设计空间；
(S2) 通过实验设计方法在设计空间生成训练样本点的输入变量，计算各训练样本点处物理模型的真实输出响应；
(S3) 融合径向基函数以及混沌多项式建立喉栓式固体姿轨控推力的增广径向基函数近似模型；
(S4) 对所述近似模型中径向基函数的形状参数进行归一化表达，将多个形状参数确定问题转化为各向异性缩放系数优化问题；
(S5) 将所述训练样本点自适应分割成K层训练样本子集，分割后的各层训练样本子集中的训练样本点均匀分布于设计空间中，采用K-折交叉验证算法对各向异性缩放系数优化问题进行求解，选择交叉验证误差最小的各向异性缩放系数并将其带回所述近似模型，得到喉栓式固体姿轨控推力快速预示模型；
(S6) 对于待测设计变量组成的输入向量，基于所述喉栓式固体姿轨控推力快速预示模型即可进行喉栓式固体姿轨控推力的快速预示。
喉栓式固体姿轨控发动机轨控喷管推力理论计算公式如下：
Figure 372558DEST_PATH_IMAGE001
其中，F为推力，C F 为推力系数，P c 为燃烧室压强，A t 为喷管等效喉部面积。
考虑沉积效应后，等效喉部面积与各喷管沉积物在喷管喉部以及喉栓头部的沉积量有关。在本发明一实施例中，所构建的考虑沉积效应后的喉栓式固体姿轨控推力物理模型为：
Figure 2
其中，F是推力，P c 是燃烧室压强，h d-Nozzle 喷管沉积量，h d-Pintle 喉栓沉积量，d Pintle 喉栓位移；
因此，将设喷管沉积量、喉栓沉积量、燃烧室压强和喉栓位移作为本发明的4个设计变量。以喉栓式固体姿轨控发动机轨控喷管物理模型，基于CFD完成数值模拟，对喉栓式固体姿轨控发动机在不同喉栓位移、喉栓沉积量、喷管沉积量、压强条件下进行仿真，得到发动机推力精细化计算结果。
本发明一实施例中，根据喷管沉积量、喉栓沉积厚度、燃烧室压强和喉栓位移范围，采用优化拉丁超立方实验设计方法在设计空间生成n个训练样本点的输入向量，并在各训练样本点的输入向量处调用所构建的考虑沉积效应后的喉栓式固体姿轨控推力物理模型，计算各训练样本的输出响应（即各训练样本点对应的真实响应值），生成训练样本集。其中输入向量由喷管沉积量、喉栓沉积量、燃烧室压强和喉栓位移四个设计变量构成。
本发明一实施例中，步骤(S3)中对于d维p次混沌多项式
Figure 663042DEST_PATH_IMAGE003
，其中，d为设计变量的数目，如前面实施例中设计变量有4个，则其d=4，p为给定数值如人为设定p=10，混沌多项式
Figure 895178DEST_PATH_IMAGE003
中待定的混沌多项式展开系数个数为m：
Figure 913950DEST_PATH_IMAGE004
，
由上式可知，完整的d维p次混沌多项式
Figure 795318DEST_PATH_IMAGE005
中的展开项个数随着设计变量个数和展开式阶数的增大而显著增加，但是混沌多项式模型具有稀疏性，可以将对输出响应影响较小的展开项系数视为0，只选择那些对输出响应具有较大影响的展开项构造稀疏正交多项混沌展开。
本发明一实施例中，步骤(S3)包括：
（S3.1)构建稀疏优化问题对d维p次混动多项式进行稀疏
Figure 495421DEST_PATH_IMAGE005
，选取混动多项式
Figure 145845DEST_PATH_IMAGE005
中多个对输出响应影响较大的项并求解其对应的混沌多项式展开系数：
Figure 69939DEST_PATH_IMAGE006
式中，
Figure 940068DEST_PATH_IMAGE007
为稀疏前的混沌多项式中待定的混沌多项式展开系数向量，
Figure 975020DEST_PATH_IMAGE008
表示
Figure 948792DEST_PATH_IMAGE009
中非零元素的个数，
Figure 43787DEST_PATH_IMAGE010
为稀疏后的混沌多项式展开系数向量；，
Figure 899748DEST_PATH_IMAGE011
为
Figure 941653DEST_PATH_IMAGE012
的转置；x i 、y i 分别为第i个训练样本点所对应的输入向量与输出响应，
Figure 799625DEST_PATH_IMAGE013
代表将x i 带入混沌多项式
Figure 65522DEST_PATH_IMAGE014
中。
通过优化算法（包括但不限于贪心算法或凸松弛算法）求解稀疏优化问题，得到稀疏后的混沌多项式展开系数向量
Figure 408778DEST_PATH_IMAGE015
，稀疏后的混沌多项式展开系数向量
Figure 51112DEST_PATH_IMAGE015
中的各个混沌多项式展开系数所对应的项即最终选取得到的混动多项式中对输出响应影响较大的项，最终选取得到的混动多项式中对输出响应影响较大的项组成稀疏后的混动多项式展开
Figure 733897DEST_PATH_IMAGE016
。本发明另一实施例中，采用贪心算法正交匹配追踪求解稀疏优化问题。
(3.2)，构建喉栓式固体姿轨控推力的增广径向基函数近似模型为：
Figure 639536DEST_PATH_IMAGE017
其中x为由喷管沉积量、喉栓沉积量、燃烧室压强和喉栓位移四个设计变量构成的输入向量，
Figure 971554DEST_PATH_IMAGE018
为近似模型在输入向量x处的预测响应值，
Figure 355262DEST_PATH_IMAGE019
为第i个训练样本点对应的径向基函数，记
Figure 689291DEST_PATH_IMAGE020
，ω i 为第i个训练样本点的径向基函数的权重系数，
Figure 296990DEST_PATH_IMAGE021
即稀疏后的混动多项式展开，
Figure 818101DEST_PATH_IMAGE022
是稀疏后的混沌多项式展开系数向量
Figure 802238DEST_PATH_IMAGE023
的第j个分量，
Figure 223730DEST_PATH_IMAGE024
为稀疏后的混动多项式展开中的项数，
Figure 205592DEST_PATH_IMAGE025
为稀疏后的混动多项式展开
Figure 745158DEST_PATH_IMAGE026
的第j项基函数。
本发明一实施例中各训练样本点对应的径向基函数采用Gauss基函数，其表达式为：
Figure 267406DEST_PATH_IMAGE027
式中：c i 为第i个训练样本点对应的Gauss基函数的形状参数，决定了训练样本点对周围空间的影响衰减程度。形状参数c越小，基函数衰减越快，其相应样本点对周围空间的影响范围越小，反之则越大。
Figure 576028DEST_PATH_IMAGE028
喉栓式固体姿轨控推力的增广径向基函数近似模型的插值条件，如下：
Figure 791108DEST_PATH_IMAGE029
构造喉栓式固体姿轨控推力的增广径向基函数近似模型的正交条件，如下：
Figure 257118DEST_PATH_IMAGE030
其中，P为
Figure 848637DEST_PATH_IMAGE031
的转置，为
Figure 11765DEST_PATH_IMAGE032
阶矩阵。ω代表所有训练样本点的径向基函数的权重系数构成的矩阵
Figure 335430DEST_PATH_IMAGE033
。
联立插值条件和正交条件，可以通过求解以下等式来计算ω：
Figure 849588DEST_PATH_IMAGE034
其中
Figure 713639DEST_PATH_IMAGE035
为径向基函数矩阵，y为输出结果向量。
对于喉栓式固体姿轨控推力的增广径向基函数近似模型中各径向基函数（Gauss基函数）的形状参数的确定问题，本发明一实施例中提出了基于样本局部密度的形状参数表征方法，将多个形状参数确定转化为缩放系数确定的单变量优化问题，实现形状参数优化与训练样本数量解耦，降低了径向基函数形状参数确定的计算复杂度。
基于Gauss基函数随形状参数的衰减特点，为充分利用样本点信息，需要增大稀疏区域样本点对周围空间的影响范围，减小密集区域样本点对周围空间的影响范围。本发明一实施例的（S4）中，令训练样本点影响范围与训练样本点局部密度成反比，即：
Figure 495387DEST_PATH_IMAGE036
式中：
Figure 52271DEST_PATH_IMAGE037
，
Figure 991408DEST_PATH_IMAGE038
分别为训练样本点x i 和x j 的影响体积。
令所有训练样本点的影响体积之和为1，求得各个训练样本点对应的Gauss函数的形状参数的基准核宽度，为：
Figure 659150DEST_PATH_IMAGE039
考虑到不同设计变量的非线性不同，需要对每一维度设计变量分配一个单独的形状参数。将训练样本点对应的Gauss函数的形状参数归一化表达为基准核宽度和各向异性缩放系数的乘积，有：
Figure 531291DEST_PATH_IMAGE040
其中，
Figure 259075DEST_PATH_IMAGE041
表示第i个训练样本点在第j维设计变量上的基准核宽度，α j 和c i 分别表示第i个训练样本点对应的各向异性缩放系数及形状参数。
为合理确定各向异性缩放系数α j ，建立各向异性缩放系数优化问题：
Figure 747825DEST_PATH_IMAGE042
式中：α为各向异性缩放系数向量；α j 为不同维度设计变量对应的各向异性缩放系数，按经验取值范围为[0.01，4]，R(.)为各向异性缩放系数向量α对应的交叉验证误差，用来评估所构建的近似模型精度。
上述实施例中，针对增广径向基函数近似模型中形状参数确定困难的问题，提出基于样本局部密度的形状参数表征方法，将多个形状参数确定转化为缩放系数确定的单变量优化问题，实现形状参数优化与训练样本数量解耦，降低了径向基函数形状参数确定的计算复杂度。
本发明一实施例的（S5）中，将所有训练样本进行切割，得到K个训练样本子集。为了使切割后的各层样本尽可能多地捕获真实模型在设计空间中的特征信息，要求各层训练样本子集中的训练样本点都具有较高的空间均匀性。
为了实现各层训练样本子集中的训练样本点都具有较高的空间均匀性，本发明一实施例采用的方法是：将划分到第k层训练样本子集中的第i个训练样本点定义为
Figure 466862DEST_PATH_IMAGE043
，将改善各层训练样本子集D k 中训练样本点空间均匀性的优化问题描述为：
Figure 459089DEST_PATH_IMAGE044
其中，k=1,2,...,K，
Figure 561037DEST_PATH_IMAGE045
，P为矩阵
Figure 5925DEST_PATH_IMAGE046
最后一列，即每个训练样本点编号；
Figure 281049DEST_PATH_IMAGE047
代表各层训练样本子集的空间均匀性准则，准则值越小，均匀性越好。采用优化算法（包括但不限于解决旅行商问题的排列优化算法）求解上述改善各层训练样本子集D k 中训练样本点空间均匀性的优化问题，得到各层训练样本子集D k 。这样使交叉验证过程中的辅助模型尽可能多地捕获真实模型在设计空间中的特征信息。
然后，用交叉验证方法对所构建的近似模型的预示结果精度进行验证，选择交叉验证误差最小的各向异性缩放系数并将其带回所述近似模型，得到喉栓式固体姿轨控推力快速预示模型。交叉验证方法可以运用本领域的常规算法，在此不再赘述。
本发明针对已知沉积规律的固体姿轨控发动机进行推力快速精确预示研究。采用实验设计方法生成在空间中均匀分布的训练样本点集，通过Fluent计算得出对应真实推力（即真实输出响应）。本发明进一步提出融合局部密度和快速交叉验证的增广径向基函数近似建模方法，对训练样本子集中的训练样本点进行分割，采用部分训练样本点作为训练集进行模型训练，剩余训练样本点作为测试集用交叉验证方法对所构建的近似模型的预示结果精度进行验证，选择交叉验证误差最小的各向异性缩放系数并将其带回所述近似模型，得到最终的喉栓式固体姿轨控推力快速预示模型。
在本发明另一实施例中，将所有训练样本点分割为K层训练样本子集，各层训练样本子集D k 中有
Figure 596624DEST_PATH_IMAGE048
行数据。K折交叉验证法(K-fold CV)误差
Figure 899167DEST_PATH_IMAGE049
为近似模型预测输出响应值与真实输出响应值之差，即：
Figure 96930DEST_PATH_IMAGE050
式中：
Figure 644586DEST_PATH_IMAGE051
为针对
Figure 876984DEST_PATH_IMAGE052
个训练样本点构建的近似模型在
Figure 851893DEST_PATH_IMAGE053
个测试样本点处的预测输出响应值；
Figure 536953DEST_PATH_IMAGE054
为测试样本点的真实输出响应值。
在K-fold CV中每折都需要训练一个新近似模型，即对一个
Figure 888299DEST_PATH_IMAGE055
阶矩阵进行求逆，整个过程需求解
Figure 414352DEST_PATH_IMAGE056
个
Figure 825742DEST_PATH_IMAGE057
阶矩阵的逆。
Figure 732518DEST_PATH_IMAGE058
和
Figure 887556DEST_PATH_IMAGE059
越大、
Figure 828967DEST_PATH_IMAGE060
越小，计算速度越慢。逆矩阵的求解占据了交叉验证过程大部分的计算时间，消耗了大量的计算资源。下面通过矩阵的快速求逆以及交叉验证的快速求解实现缩放系数的快速确定。
令L o 为维度为o的方阵，其中o=m+n。将L o 按照下式进行分块，
Figure 348941DEST_PATH_IMAGE061
通过左乘变行，右乘变列的原则对原系数矩阵L o 进行行列交换，将应作为测试子集的第i组训练样本子集与第1组训练样本子集换位后，得到新的系数矩阵
Figure 743013DEST_PATH_IMAGE062
。I(i,1)定义为第i组训练样本子集与第1组训练样本子集交换后得到的矩阵，L o 与
Figure 934698DEST_PATH_IMAGE062
关系如下式所示：
Figure 730616DEST_PATH_IMAGE063
Figure 483808DEST_PATH_IMAGE064
系数矩阵L o 中的
Figure 630756DEST_PATH_IMAGE065
被划分为k组后，设当前测试子集内有r个样本，将第i组训练样本子集与第1组训练样本子集交换位置，矩阵中的每个元素代表一个r×r的小矩阵，例β i,i 表示数据集
Figure 862017DEST_PATH_IMAGE065
中第i组数据。
对
Figure 715703DEST_PATH_IMAGE066
求逆，得
Figure 639797DEST_PATH_IMAGE067
由
Figure 775506DEST_PATH_IMAGE068
，可以得到下式：
Figure 13720DEST_PATH_IMAGE069
将
Figure 518651DEST_PATH_IMAGE070
分为4块，即
Figure 613646DEST_PATH_IMAGE071
，留出原矩阵的最前r行和最前r列。其中，A是r阶方阵B是r×(o-r)阶矩阵，C是(o-r)×r阶矩阵，D是(o-r)阶方阵，则
Figure 735185DEST_PATH_IMAGE072
由于
Figure 308249DEST_PATH_IMAGE073
由L o 行列变换所得，因此A、B和C包含测试样本的数据，D即为剩余样本训练近似模型的系数矩阵
Figure 635063DEST_PATH_IMAGE074
。
通过矩阵的求逆公式可知：
Figure 900959DEST_PATH_IMAGE075
采用上式即可直接计算出
Figure 978637DEST_PATH_IMAGE076
，运用r阶矩阵求逆替代了(o-r)阶矩阵直接求逆的过程，由于r＜＜o，因此转化后的矩阵求逆计算量会大幅降低。
K折交叉验证(K-fold CV)的误差
Figure 620971DEST_PATH_IMAGE077
可以表示为：
Figure 366073DEST_PATH_IMAGE078
式中：Y pre,r 为近似模型在r个测试样本处的预测值；Y r 为测试样本真实输出；ω r为测试样本点对应的r个权系数。
本发明另一实施例中，物理建模的边界设置如图2所示，其中P_in为喷管的入口边界条件，采用压力入口，总压为输入条件，总温3200K；Axis为对称边界；Pintle_wall和Nozzle_wall均采用绝热无滑移壁面边界条件，其型面根据输入条件可自动进行变化；P_out为喷管的出口边界条件，采用压力出口，压强101325Pa，环境温度288.15K。
采用200组仿真数据作为训练样本参数，使用本发明所提出的融合局部密度和快速交叉验证的增广径向基方法，构建固体姿轨控发动机推力预示模型。在选择缩放系数时，采用5-折交叉验证，每折中应有40个数据。在单次交叉验证过程中，160个样本点为测试数据，40个样本点为验证数据。在推力预示模型基础上完成20组阀门任意喉栓位移及不同沉积量下的推力快速精确预示，再对比仿真数据，验证预测精度。
采用随机生成的另外20组样本作为测试样本，进行数值模拟得到其实际推力，对比实际推力与模型预测推力偏差。以测试样本归一后的均方根误差(Normalized RootMean Square Error, NRMSE)作为模型预测精度评价指标。
Figure 802870DEST_PATH_IMAGE079
其中，n表示样本点个数，y(x i )表示样本点对应实际推力，
Figure 367844DEST_PATH_IMAGE080
表示模型预测推力。
采用20组数据对固体姿轨控发动机推力预示结果进行验证，各组数据相对误差如表1所示，NRMSE为0.0029。单次CFD计算平均耗时为0.5小时，采用预示模型仅需0.2秒，验证了本发明的计算能快速性。
表1 推力预示相对误差表
Figure 987437DEST_PATH_IMAGE081
由表1和图3可知，本方法可以很好地对固体姿轨控发动机推力进行快速预示，其相对误差均在2.5%以下，充分验证了本发明的有效性。本发明大大降低了固体姿轨控推力预示所需要的仿真次数，节约了时间和精力，大幅度提高了计算效率；本发明建模精度高、计算速度快，大幅度降低了增广径向基代理模型构建的计算量。
本发明已经在多种不同的固体姿轨控发动机设计中得到了应用，与现有方法相比，人工参与过程少、设计速度快，自动化程度高，设计结果性能可靠。
在本实施例中，提供了一种计算机设备，该计算机设备可以是服务器，其内部结构图可以如图4所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口和数据库。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统、计算机程序和数据库。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的数据库用于存储样本数据。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现上述实施例中喉栓式固体姿轨控推力快速预示方法的步骤。
在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述实施例中喉栓式固体姿轨控推力快速预示方法的步骤。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器（ROM）、可编程ROM（PROM）、电可编程ROM（EPROM）、电可擦除可编程ROM（EEPROM）或闪存。易失性存储器可包括随机存取存储器（RAM）或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM（SRAM）、动态RAM（DRAM）、同步DRAM（SDRAM）、双数据率SDRAM（DDRSDRAM）、增强型SDRAM（ESDRAM）、同步链路（Synchlink） DRAM（SLDRAM）、存储器总线（Rambus）直接RAM（RDRAM）、直接存储器总线动态RAM（DRDRAM）、以及存储器总线动态RAM（RDRAM）等。
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。

Claims (10)

1.喉栓式固体姿轨控推力快速预示方法，其特征在于，包括：
确定考虑沉积效应后的喉栓式固体姿轨控推力物理模型、设计变量、约束条件以及设计空间；
通过实验设计方法在设计空间生成训练样本点的输入变量，计算各训练样本点处物理模型的真实输出响应；
融合径向基函数以及混沌多项式建立喉栓式固体姿轨控推力的增广径向基函数近似模型；
对所述近似模型中径向基函数的形状参数进行归一化表达，将多个形状参数确定问题转化为各向异性缩放系数优化问题；
将所述训练样本点自适应分割成K层训练样本子集，分割后的各层训练样本子集中的训练样本点均匀分布于设计空间中，采用K-折交叉验证算法对各向异性缩放系数优化问题进行求解，选择交叉验证误差最小的各向异性缩放系数并将其带回所述近似模型，得到喉栓式固体姿轨控推力快速预示模型；
对于待测设计变量组成的输入向量，基于所述喉栓式固体姿轨控推力快速预示模型即可进行喉栓式固体姿轨控推力的快速预示。
2.根据权利要求1所述的喉栓式固体姿轨控推力快速预示方法，其特征在于，设计变量包括4个，分别是喷管沉积量、喉栓沉积量、燃烧室压强和喉栓位移。
3.根据权利要求2所述的喉栓式固体姿轨控推力快速预示方法，其特征在于，确定考虑沉积效应后的喉栓式固体姿轨控推力物理模型为：
Figure 4
其中，F是推力，P c 是燃烧室压强，h d-Nozzle 喷管沉积量，h d-Pintle 喉栓沉积量，d Pintle 喉栓位移；
基于CFD完成数值模拟，对固体姿轨控发动机在不同喉栓位移、喉栓沉积量、喷管沉积量、压强条件下进行仿真，得到各训练样本点处物理模型的真实输出响应。
4.根据权利要求1、2或3所述的喉栓式固体姿轨控推力快速预示方法，其特征在于，对于d维p次混沌多项式
Figure 994957DEST_PATH_IMAGE002
，其中，d为设计变量的数目，p为给定数值，混沌多项式
Figure 756719DEST_PATH_IMAGE002
中待定的混沌多项式展开系数个数为m，
Figure 824032DEST_PATH_IMAGE003
，构建稀疏优化问题对d维p次混动多项式进行稀疏，选取混动多项式中多个对输出响应影响较大的项以及对应的混沌多项式展开系数：
Figure 107246DEST_PATH_IMAGE004
式中，
Figure 664129DEST_PATH_IMAGE005
为稀疏前的混沌多项式中待定的混沌多项式展开系数向量，
Figure 400004DEST_PATH_IMAGE006
表示a中非零元素的个数，
Figure 67746DEST_PATH_IMAGE007
为稀疏后的混沌多项式展开系数向量，
Figure 907263DEST_PATH_IMAGE008
为
Figure 369469DEST_PATH_IMAGE009
的转置；x i 、y i 分别为第i个训练样本点所对应的输入向量与输出响应，
Figure 592640DEST_PATH_IMAGE010
代表将x i 带入混沌多项式
Figure 64072DEST_PATH_IMAGE011
中；
通过优化算法求解稀疏优化问题，得到稀疏后的混沌多项式展开系数向量
Figure 259561DEST_PATH_IMAGE012
，稀疏后的混沌多项式展开系数向量
Figure 423827DEST_PATH_IMAGE012
中的各个混沌多项式展开系数所对应的项即最终选取得到的混动多项式中对输出响应影响较大的项，最终选取得到的混动多项式中对输出响应影响较大的项组成稀疏后的混动多项式展开
Figure 134294DEST_PATH_IMAGE013
。
5.根据权利要求4所述的喉栓式固体姿轨控推力快速预示方法，其特征在于，由贪心算法或凸松弛算法求解稀疏优化问题。
6.根据权利要求4所述的喉栓式固体姿轨控推力快速预示方法，其特征在于，喉栓式固体姿轨控推力的增广径向基函数近似模型为：
Figure 582986DEST_PATH_IMAGE014
其中x为由喷管沉积量、喉栓沉积量、燃烧室压强和喉栓位移四个设计变量构成的输入向量，
Figure 695298DEST_PATH_IMAGE015
为近似模型在输入向量x处的预测响应值，
Figure 764886DEST_PATH_IMAGE016
为第i个训练样本点对应的径向基函数，ω i 为第i个训练样本点的径向基函数的权重系数，
Figure 962649DEST_PATH_IMAGE017
即稀疏后的混动多项式展开，
Figure 713567DEST_PATH_IMAGE018
是稀疏后的混沌多项式展开系数向量
Figure 680386DEST_PATH_IMAGE019
的第j个分量，
Figure 419410DEST_PATH_IMAGE020
为稀疏后的混动多项式展开中的项数，
Figure 104469DEST_PATH_IMAGE021
为稀疏后的混动多项式展开
Figure 455816DEST_PATH_IMAGE022
的第j项基函数。
7.根据权利要求6所述的喉栓式固体姿轨控推力快速预示方法，其特征在于，根据近似模型的插值条件和正交条件，计算由各个训练样本点的基函数的权重系数组成的系数矩阵。
8.根据权利要求6所述的喉栓式固体姿轨控推力快速预示方法，其特征在于，各训练样本点对应的径向基函数采用Gauss基函数；
令训练样本点影响范围与训练样本点局部密度成反比，所有训练样本点的影响体积之和为1，求得各个训练样本点对应的Gauss函数的形状参数的基准核宽度；
将训练样本点对应的Gauss函数的形状参数归一化表达为基准核宽度和各向异性缩放系数的乘积。
9.根据权利要求1或2或3或5或6或7或8所述的喉栓式固体姿轨控推力快速预示方法，其特征在于，将划分到第k层训练样本子集中的第i个训练样本点定义为
Figure 480404DEST_PATH_IMAGE023
，将改善各层训练样本子集D k 中训练样本点空间均匀性的优化问题描述为：
Figure 891793DEST_PATH_IMAGE024
其中，
Figure 798570DEST_PATH_IMAGE025
，P为矩阵
Figure 953607DEST_PATH_IMAGE026
最后一列，即每个训练样本点编号；
Figure 599746DEST_PATH_IMAGE027
代表各层训练样本子集的空间均匀性准则，准则值越小，均匀性越好，求解上述优化问题得到各层训练样本子集。
10.根据权利要求9所述的喉栓式固体姿轨控推力快速预示方法，其特征在于，采用解决旅行商问题的排列优化算法求解改善各层训练样本子集D k 中训练样本点空间均匀性的优化问题。
CN202210137215.2A 2022-02-15 2022-02-15 喉栓式固体姿轨控推力快速预示方法 Active CN114200843B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210137215.2A CN114200843B (zh) 2022-02-15 2022-02-15 喉栓式固体姿轨控推力快速预示方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210137215.2A CN114200843B (zh) 2022-02-15 2022-02-15 喉栓式固体姿轨控推力快速预示方法

Publications (2)

Publication Number Publication Date
CN114200843A CN114200843A (zh) 2022-03-18
CN114200843B true CN114200843B (zh) 2022-04-19

Family

ID=80658998

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210137215.2A Active CN114200843B (zh) 2022-02-15 2022-02-15 喉栓式固体姿轨控推力快速预示方法

Country Status (1)

Country Link
CN (1) CN114200843B (zh)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20190057614A (ko) * 2017-11-20 2019-05-29 주식회사 한화 고체 추진기관의 가변성 추력성능 시험장치
CN111783251A (zh) * 2020-07-16 2020-10-16 中国人民解放军国防科技大学 一种固体火箭发动机总体参数设计方法
CN112464387A (zh) * 2021-01-26 2021-03-09 中国人民解放军国防科技大学 一种喉栓式变推力固体发动机推力匹配设计方法
CN112528441A (zh) * 2021-02-18 2021-03-19 中国人民解放军国防科技大学 喉栓式变推力发动机总体参数设计方法、装置和设备
CN113255145A (zh) * 2021-06-03 2021-08-13 中国人民解放军国防科技大学 姿轨控发动机喉栓运动时序实验设计方法、装置和设备
CN113255235A (zh) * 2021-06-28 2021-08-13 中国人民解放军国防科技大学 飞行器复杂结构近似建模方法、装置、设备和介质
CN113496062A (zh) * 2021-07-26 2021-10-12 中国人民解放军国防科技大学 压强平稳约束下姿轨控发动机喉栓运动时序设计方法
CN113779704A (zh) * 2021-09-28 2021-12-10 中国人民解放军国防科技大学 喉栓式变推力固体火箭发动机喷管型面多目标优化方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20190057614A (ko) * 2017-11-20 2019-05-29 주식회사 한화 고체 추진기관의 가변성 추력성능 시험장치
CN111783251A (zh) * 2020-07-16 2020-10-16 中国人民解放军国防科技大学 一种固体火箭发动机总体参数设计方法
CN112464387A (zh) * 2021-01-26 2021-03-09 中国人民解放军国防科技大学 一种喉栓式变推力固体发动机推力匹配设计方法
CN112528441A (zh) * 2021-02-18 2021-03-19 中国人民解放军国防科技大学 喉栓式变推力发动机总体参数设计方法、装置和设备
CN113255145A (zh) * 2021-06-03 2021-08-13 中国人民解放军国防科技大学 姿轨控发动机喉栓运动时序实验设计方法、装置和设备
CN113255235A (zh) * 2021-06-28 2021-08-13 中国人民解放军国防科技大学 飞行器复杂结构近似建模方法、装置、设备和介质
CN113496062A (zh) * 2021-07-26 2021-10-12 中国人民解放军国防科技大学 压强平稳约束下姿轨控发动机喉栓运动时序设计方法
CN113779704A (zh) * 2021-09-28 2021-12-10 中国人民解放军国防科技大学 喉栓式变推力固体火箭发动机喷管型面多目标优化方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
喉栓式推力可调固体火箭发动机动态响应特性数值分析;李娟 等;《固体火箭技术》;20090228;第32卷(第1期);第48-52页 *
基于响应面法的喉栓式喷管型面优化设计;成沉 等;《固体火箭技术》;20140215;第37卷(第01期);第47-51页 *

Also Published As

Publication number Publication date
CN114200843A (zh) 2022-03-18

Similar Documents

Publication Publication Date Title
CN112434448B (zh) 基于多点加点的代理模型约束优化方法和装置
CN108304679A (zh) 一种自适应可靠性分析方法
CN112733997B (zh) 基于woa-lstm-mc的水文时间序列预测优化方法
CN111783209B (zh) 一种学习函数与kriging模型结合的自适应结构可靠性分析方法
CN113806991B (zh) 一种发动机燃烧噪声优化预测方法、装置及存储介质
CN115221801B (zh) 基于动态近似建模的飞行器不确定性传播分析方法和装置
CN113241132B (zh) 一种材料辐照损伤预测机器学习模型的建立方法和系统
Cheng et al. A multi-fidelity surrogate modeling method based on variance-weighted sum for the fusion of multiple non-hierarchical low-fidelity data
CN111444649B (zh) 基于强度折减法的边坡系统可靠度分析方法
CN114048544B (zh) 飞行器结构响应特性预示方法、装置、设备和介质
CN114200843B (zh) 喉栓式固体姿轨控推力快速预示方法
CN111339487A (zh) 基于径向基函数rbf的边坡系统失效概率计算方法
CN112906295B (zh) 固体火箭发动机燃速预示方法、装置和设备
CN113255082B (zh) 变推力固体发动机喉栓喷管型面一体化优化设计方法
CN114036805A (zh) 三维稳态热传导地温场的正演方法、装置、设备及介质
CN113312712B (zh) 一种面向飞行器优化设计的递归排列演化实验设计方法
CN109540089B (zh) 一种基于贝叶斯-克里金模型的桥面高程拟合方法
CN114934848B (zh) 一种面向柴油机燃烧性能优化控制的模糊神经网络建模方法
CN114036806B (zh) 基于热导率各向异性介质的三维地温场数值模拟方法
Agarwal et al. Spectral state space models
Guan et al. A data‐driven time‐step determination approach for dynamic simulation of heat‐electric coupled system
Zhou et al. An on-line multi-fidelity metamodel assisted multi-objective genetic algorithm
CN111651904B (zh) 基于机器学习的流场数据计算方法和装置
CN114595577A (zh) 基于克里金代理模型的分阶段高效约束优化方法和装置
Praneeth et al. Uncertainty quantification of tabulated supercritical thermodynamics for compressible Navier-Stokes solvers

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
GR01 Patent grant
GR01 Patent grant