CN113110559B - 一种小天体表面弹跳运动最优控制方法 - Google Patents

一种小天体表面弹跳运动最优控制方法 Download PDF

Info

Publication number
CN113110559B
CN113110559B CN202110522720.4A CN202110522720A CN113110559B CN 113110559 B CN113110559 B CN 113110559B CN 202110522720 A CN202110522720 A CN 202110522720A CN 113110559 B CN113110559 B CN 113110559B
Authority
CN
China
Prior art keywords
celestial body
small celestial
thrust
model
optimization
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.)
Expired - Fee Related
Application number
CN202110522720.4A
Other languages
English (en)
Other versions
CN113110559A (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.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
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 Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN202110522720.4A priority Critical patent/CN113110559B/zh
Publication of CN113110559A publication Critical patent/CN113110559A/zh
Application granted granted Critical
Publication of CN113110559B publication Critical patent/CN113110559B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05D—SYSTEMS FOR CONTROLLING OR REGULATING NON-ELECTRIC VARIABLES
    • G05D1/00—Control of position, course, altitude or attitude of land, water, air or space vehicles, e.g. using automatic pilots
    • G05D1/08—Control of attitude, i.e. control of roll, pitch, or yaw
    • G05D1/0808—Control of attitude, i.e. control of roll, pitch, or yaw specially adapted for aircraft
    • G05D1/0816—Control of attitude, i.e. control of roll, pitch, or yaw specially adapted for aircraft to ensure stability
    • G05D1/0833—Control of attitude, i.e. control of roll, pitch, or yaw specially adapted for aircraft to ensure stability using limited authority control
    • 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/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Software Systems (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Medical Informatics (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Data Mining & Analysis (AREA)
  • Automation & Control Theory (AREA)
  • Computing Systems (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Mathematical Physics (AREA)
  • Computer Security & Cryptography (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Feedback Control In General (AREA)

Abstract

本发明公开的一种小天体表面弹跳运动最优控制方法，属于深空探测技术领域。本发明首先，对不同位置的探测器表面弹跳运动轨迹进行优化求解，获得样本数据；接着对预测模型参数优化选取，通过对样本数据进行训练，进一步建立探测器最优控制推力的预测模型；由于该方法避免了复杂的建模和大规模的运算过程，从统计角度直接建立控制量与初值之间的映射关系，能够迅速且相对精确地得到推力值，因此求解效率得到提升。利用最优控制推力的预测计算结果能够快速得到探测器运动的优化轨迹。

Description

一种小天体表面弹跳运动最优控制方法
技术领域
本发明涉及一种小天体表面弹跳运动最优控制方法，属于深空探测技术领域。
背景技术
小天体探测有很多种形式，而对小天体进行表面探测可以帮助人类更加深入地了解小天体特性。由于小天体附近引力微弱，表面不能提供足够的摩擦力，传统的轮动式探测器无法在小天体表面有效移动，因此小天体表面探测通常采用弹跳移动方式。而探测器是通过控制推力进行运动的，探测器燃料储备有限，降低燃料消耗对探测任务的开展具有重要的意义。为了实现弹跳表面探测方案，计算出最优控制推力，学者们都做出了很多贡献，然而现有的小天体表面弹跳运动控制方法，专注于复杂的动力学建模，无法避免计算量大的问题，或者简化了模型从而降低了精度。如果想要获取小天体表面弹跳运动最优控制推力，仍需要花费大量的时间。
机器学习中通过挖掘数据特征给出模型输入输出关系的方法称为监督学习。监督学习是一种常用的机器学习算法，从训练集学习输入输出之间的映射关系。确定输入输出映射关系的方法通常有两种，一种为参数化回归，另一种为贝叶斯回归。贝叶斯回归定义了一个函数分布，赋予每一种可能的函数一个先验概率。通过假定的噪声分布可以得到训练集的似然函数。通过使似然函数最大化来得到输入输出的映射关系，无需找出一个具体的函数关系。本文所采用的方法建立于贝叶斯回归方法，能够利用已知的一部分观测值对输入输出模型进行训练学习，找到模型输入输出的映射关系，从而在给定新的输入值时，能够迅速的给出对应的输出结果。该方法不需要逐次寻优计算过程，避免了大规模的运算过程，从统计角度直接建立输入输出之间的映射关系，能够迅速且相对精确地得到预测值，因此求解效率得到提升。
发明内容
本发明的目的是为了解决小天体表面弹跳运动控制方法存在的计算量大、效率低的问题，提供一种小天体表面弹跳运动最优控制方法，该方法从数据统计规律挖掘的全新角度出发，挖掘探测器运动初始状态与控制推力之间映射关系的规律，将机器学习的思想引入到小天体表面弹跳运动最优控制推力的求解过程，能够显著提高计算效率。此外，应用本发明的小天体表面弹跳运动最优控制方法的计算结果，还能够进行探测器在小天体表面弹跳运动的轨迹规划，快速得到运动的最优轨迹。
本发明的目的是通过下述技术方案实现的。
本发明公开的一种小天体表面弹跳运动最优控制方法，从数据统计规律挖掘的全新角度出发，将机器学习的思想应用到小天体表面的弹跳运动上，寻找探测器运动初始状态与控制推力之间映射关系的规律；将机器学习的算法引入到小天体表面弹跳最优控制推力的求解过程中；首先，对不同位置的探测器表面弹跳运动轨迹进行优化求解，获得样本数据；接着对模型参数优化选取，通过对样本数据进行训练，进一步建立探测器最优控制推力的预测模型；由于该方法避免了复杂的建模和大规模的运算过程，从统计角度直接建立控制量与初值之间的映射关系，能够迅速且相对精确地得到推力值，因此求解效率得到提升。利用最优控制推力的预测计算结果能够快速得到探测器运动的优化轨迹。
一种小天体表面弹跳运动最优控制方法，包括如下步骤：
步骤1：建立小天体表面动力学模型。
考虑小天体的转动，采用如下动力学模型：
Figure BDA0003064692090000021
其中，r表示探测器相对于小天体质心的位置矢量，ω是旋转角速度矢量，v表示速度矢量，m表示探测器质量，T表示发动机推力矢量，g表示引力加速度矢量，Isp表示比冲，ge表示地球引力加速度常数。
通过牛顿恢复系数来描述碰撞前后法向状态之间的关系，求得碰撞后法向，切向速度的方法为：
牛顿恢复系数表示为：
e(v)＝-v+/v- (2)
其中v+和v-分别表示碰撞后和碰撞前物体的相对法向速度。
探测器会因为与小天体表面碰撞而发生轻微形变，因此定义恢复系数e。其中vn0和vn1是碰撞前后的法向速度，碰撞后法向速度是
vn1＝-evn0 (3)
在切向方向，通过引入瞬间摩擦来描述碰撞过程切向关系的变化，碰撞后的切向速度表示为
vt1＝vt0-μ(1-e)vn0 (4)
其中μ是小天体表面摩擦系数，vt0和vt1是碰撞前后的切向速度
步骤2：设计用于得到小天体表面弹跳运动最优控制推力的优化模型，实现小天体表面弹跳运动轨迹的优化求解。
轨迹优化它至少包含动力学模型和各种约束条件即初末状态约束及路径约束，在给定动力学模型的情况下，寻找满足各种约束条件的可行解。由于小天体表面探测器携带的燃料有限，降低燃耗对探测任务的开展具有重要的意义，因此以燃耗最优为目标函数进行优化求解。探测器在表面弹跳运动过程中需要满足动力学约束、初始和末端状态约束、路径约束以及推力幅值约束条件，将控制变量T作为优化变量，进行性能指标的计算。
步骤2具体实现方法为：
步骤2.1：设计小天体表面弹跳运动最优控制推力求解模型，通过优化方法得到小天体表面弹跳运动最优控制推力。
步骤2.1.1：建立小天体表面弹跳运动最优控制推力模型目标函数J。
由于小天体表面探测器携带的燃料有限，降低燃耗对探测任务的开展具有重要的意义，因此以燃耗最优为目标函数进行优化设计。
Figure BDA0003064692090000031
其中，tf表示飞行时间。
步骤2.1.2：确立天体表面弹跳运动最优控制推力模型过程约束。
在弹跳过程中使每次弹跳的路径约束为
Figure BDA0003064692090000041
Xi和Xf表示探测器的初末状态，各矩阵系数定义为：
Figure BDA0003064692090000042
Figure BDA0003064692090000043
r0和rf表示探测器初始位置和目标位置，
Figure BDA0003064692090000044
和ε表示探测器与地面的角度。
推力幅值约束：
||T||≤Tm (7)
Tm表示发动机最大推力。
步骤2.1.3：确立天体表面弹跳运动最优控制推力模型初始和终端约束。
满足的初始状态约束为：
r(t0)＝r0，v(t0)＝v0，m(t0)＝m0 (8)
终端状态约束：
r(tf)＝rf，v(tf)＝vf (9)
选定下一次弹跳的位置，根据碰撞动力学可以求出末端速度vf。v0和m0表示探测器初始速度和初始质量，vf表示末端速度，t0表示初始时间。
步骤2.2：通过优化方法实现小天体表面弹跳最优控制推力的计算。
将发动机推力矢量T作为优化变量，通过优化方法对初末两个时刻的控制推力进行优化，通过插值得到优化的控制推力曲线，再利用步骤1所建立的小天体表面动力学模型进行目标函数J的计算，实现小天体表面弹跳最优控制推力的计算。
步骤2.2所述的优化方法包括遗传算法、凸优化、高斯伪谱法等。凸优化问题的局部最优解就是全局最优解，从而保证优化结果的质量，步骤2.2所述的优化方法采用基于凸优化的优化求解器。
步骤2.3：将式(1)的问题转化为凸优化问题进行求解，实现小天体表面弹跳运动轨迹的优化求解。
Figure BDA0003064692090000051
Tx、Ty和Tz分别表示三轴的推力，Γ是松弛变量。
转化为凸优化问题进行求解，首先引入松弛变量，接着是更改变量，最终实现凸化。其中选择代表推力大小的松弛变量Γ对||T||进行替换，松弛后的约束条件为：
将式(1)线性化得到：
Figure BDA0003064692090000052
式(7)的不等式约束可以写成：
Figure BDA0003064692090000053
其中，γ0(t)＝ln(m0-Tmt/Ispge)，优化性能指标写成：
Figure BDA0003064692090000054
通过以上线性化处理，消除了式(1)中探测器质量m带来的非线性问题。
将时间区间离散化，分成N份，线性化后的动力学方程可以写成
Figure BDA0003064692090000055
的形式。
Figure BDA0003064692090000056
Figure BDA0003064692090000057
引力加速度矢量g是关于探测器相对于小天体质心的位置矢量r的非线性函数，这使得离散后的问题不是一个标准的二阶锥规划问题。采用逐次求解方法，以反复近似方程中的非线性动力学。令第(k-1)个位置矢量r(k-1)代替第k个位置矢量rk，将得到的新位置，作为下一次求解的参考值，反复迭代直到收敛。本文设定的收敛误差为0.5m,可以通过改变迭代次数来减小误差，将主导项放在A会减少迭代次数。
Figure BDA0003064692090000061
通过松弛约束，动力学线性化，引力加速度通过逐次求解来处理，对轨迹优化问题进行了转化，最终变成了以下形式：
Figure BDA0003064692090000062
实现小天体表面弹跳运动轨迹的优化求解。
步骤3：设计采用机器学习的小天体表面弹跳运动最优控制推力预测计算模型，实现基于数据统计规律挖掘的小天体表面弹跳运动最优轨迹计算。
利用机器学习的方法对行小天体表面弹跳最优控制推力预测模型进行设计，建立探测器初始状态和推力之间的输入输出模型。
首先在小天体表面一定范围内随机取点作为学习样本，取这些点的坐标作为训练模型的输入向量，训练模型的输出为初末位置处的控制推力，推力由文章第二部分的方法计算得出。xi和yi分别是模型的输入和输出，具体表示如下：
xi＝[r0,rt] 为模型的输入
yi＝Ti 为模型的输出
接着选择适当的均值函数和核函数，设计预测模型。选用零均值函数和平方指数协方差函数，表达式为：
Figure BDA0003064692090000063
加入噪声后，k(x,x')为：
Figure BDA0003064692090000071
结合联合分布，可以得到预测数据和已知训练数据的联合正态分布：
Figure BDA0003064692090000076
从而计算出预测数据的均值和方差，分布的均值作为y*的估计值即：
Figure BDA0003064692090000072
训练集中的学习样本输入xi＝[r0i,rti],i＝1,2,…,n,学习样本输出yi＝Ti,i＝1,2…n。检验集中的样本输入为x*＝[r0*,rt*]，输出为y*＝T*。于是,根据上式就能够得到y*的预测值。
通过训练优化模型超参数，即似然函数值最大时，就把模型超参数作为参数估计的输出；似然函数表示为：
Figure BDA0003064692090000073
在求log(p(y|X,θ))的最大值时，对其求导并用最优化方法求得概率最大时的θ值，以使得观测的结果无限接近己知样本；令θ为所有模型超参数的集合，共扼梯度法通过求取训练样本对数似然函数的极大值来获得最优超参数
Figure BDA0003064692090000074
具体如下：
Figure BDA0003064692090000075
最后，利用优化求解出的最优超参数，确定小天体表面弹跳运动推力预测计算模型。使在给定探测任务的某些初始参数时，能够利用预测模型快速、准确的对弹跳探测器最优控制推力进行预测计算。
还包括步骤4：利用步骤3所述的采用机器学习方法的小天体表面弹跳运动最优控制推力预测计算，实现对不同探测位置探测器最优控制推力进行预测计算。
还包括步骤5：利用步骤4预测计算结果对小天体表面弹跳运动的最优轨迹进行计算，进而解决小天体表面探测的相关工程问题。
有益效果
1.本发明公开的一种小天体表面弹跳运动最优控制方法，从数据统计规律挖掘的全新角度出发，建立探测器最优控制推力的预测模型；由于该方法避免了复杂的建模和大规模的运算过程，从统计角度直接建立控制量与初值之间的映射关系，能够迅速且相对精确地得到推力值，因此求解效率得到提升。利用本发明公开的方法计算最优控制推力，平均相对误差在5％以内，计算时间都小于1s，而传统优化方法至少需要几百秒以上，由此本发明所公开的方法在保证准确性的同时，效率大大提高。
2.本发明公开的小天体表面弹跳运动最优控制方法，采用凸优化算法作为优化求解，凸优化问题的局部最优解就是全局最优解，从而保证优化结果的质量。
3.应用本发明公开的小天体表面弹跳运动最优控制方法，能够对小天体表面弹跳运动最优控制推力预测计算，实现对不同探测位置探测器最优控制推力进行高效预测计算，利用计算结果可以对轨迹进行规划，完成最优轨迹的快速计算，从而解决小天体表面探测的相关工程问题。
附图说明
图1为本发明公开的小天体表面弹跳运动最优控制方法流程图；
图2为初始推力预测相对误差随训练数据大小变化曲线；
图3为末端推力预测相对误差随训练数据大小变化曲线。
图4为不同训练样本下检验样本的误差分布图，其中图a为200个训练样本时误差分布；图b为400个训练样本时误差分布；图c为800个训练样本时误差分布；图d为1000个训练样本时误差分布。
具体实施方式
为了更好的说明本发明的目的和优点，下面结合实例对发明内容做进一步说明。
本实例针对小天体表面弹跳运动最优控制方法，首先利用凸优化对小天体表面弹跳运动进行优化求解，以提供高质量的学习样本；其次选取合适的样本相关性描述参数，建立基于机器学习的小天体表面弹跳运动最优推力预测模型；最后验证本文所提出的预测模型的有效性，具体流程图如图1所示。
本实例公开的一种小天体表面弹跳运动最优控制方法，具体实施方式如下：
步骤1：建立小天体表面动力学模型。
考虑小天体的转动，采用如下动力学模型：
Figure BDA0003064692090000091
Figure BDA0003064692090000092
Figure BDA0003064692090000093
其中，r表示探测器相对于小天体质心的位置矢量，ω是旋转角速度矢量，v表示速度矢量，m表示探测器质量，T表示发动机推力矢量，g表示引力加速度矢量，Isp表示比冲，ge表示地球引力加速度常数。ω是0.00016559rad/s，发动机比冲Isp为300s，ge是9.8N/kg。
通过牛顿恢复系数来描述碰撞前后法向状态之间的关系，求得碰撞后法向，切向速度的方法为：
牛顿恢复系数表示为：
e(v)＝-v+/v-
其中v+和v-分别表示碰撞后和碰撞前物体的相对法向速度。
探测器会因为与小天体表面碰撞而发生轻微形变，因此定义恢复系数e。其中vn0和vn1是碰撞前后的法向速度，碰撞后法向速度是
vn1＝-evn0
在切向方向，通过引入瞬间摩擦来描述碰撞过程切向关系的变化，碰撞后的切向速度表示为
vt1＝vt0-μ(1-e)vn0
其中μ是小天体表面摩擦系数，vt0和vt1是碰撞前后的切向速度
其中e＝0.8，μ＝0.8，vn0和vn1是碰撞前后的法向速度。
步骤2：设计用于得到小天体表面弹跳运动最优控制推力的优化模型，实现小天体表面弹跳运动轨迹的优化求解。
轨迹优化它至少包含动力学模型和各种约束条件即初末状态约束及路径约束，在给定动力学模型的情况下，寻找满足各种约束条件的可行解。由于小天体表面探测器携带的燃料有限，降低燃耗对探测任务的开展具有重要的意义，因此以燃耗最优为目标函数进行优化求解。探测器在表面弹跳运动过程中需要满足动力学约束、初始和末端状态约束、路径约束以及推力幅值约束条件，将控制变量T作为优化变量，进行性能指标的计算。
步骤2.1：设计小天体表面弹跳运动最优控制推力求解模型，通过优化方法得到小天体表面弹跳运动最优控制推力。
步骤2.1.1：建立小天体表面弹跳运动最优控制推力模型目标函数J。
由于小天体表面探测器携带的燃料有限，降低燃耗对探测任务的开展具有重要的意义，因此以燃耗最优为目标函数进行优化设计。
Figure BDA0003064692090000101
其中飞行时间tf＝200s。
步骤2.1.2：确立天体表面弹跳运动最优控制推力模型过程约束。
在弹跳过程中使每次弹跳的路径约束为
||CXi-D1r0||+e1Xi-D2r0≤0
||CXf-D1rf||+e2Xf-D2rf≤0
Xi和Xf表示探测器的初末状态，各矩阵系数定义为：
Figure BDA0003064692090000102
Figure BDA0003064692090000103
推力幅值约束：
||T||≤Tm
其中
Figure BDA0003064692090000111
度，ε＝45度，发动机最大推力Tm为200N
步骤2.1.3：确立天体表面弹跳运动最优控制推力模型初始和终端约束。
满足的初始状态约束为：
r(t0)＝r0，v(t0)＝v0，m(t0)＝m0
终端状态约束：
r(tf)＝rf，v(tf)＝vf
选定下一次弹跳的位置，根据碰撞动力学可以末端速度vf。
其中探测器的初始位置r0＝[695.3 435.8 5484.1]m，探测器目标位置rf＝[795.2425.75481.0]m，末端速度vf＝[0.4383 1.4617 -0.1085]，探测器初始质量m0＝50kg。
步骤2.2：通过优化方法实现小天体表面弹跳最优控制推力的计算。
将发动机推力矢量T作为优化变量，通过优化方法对初末两个时刻的控制推力进行优化，通过插值得到优化的控制推力曲线，再利用步骤1所建立的小天体表面动力学模型进行目标函数J的计算，实现小天体表面弹跳最优控制推力的计算。
步骤2.2所述的优化方法包括遗传算法、凸优化、高斯伪谱法等。凸优化问题的局部最优解就是全局最优解，从而保证优化结果的质量，步骤2.2所述的优化方法采用基于凸优化的优化求解器。
步骤2.3：将式(1)的问题转化为凸优化问题进行求解，实现小天体表面弹跳运动轨迹的优化求解。
Figure BDA0003064692090000112
Tx、Ty和Tz分别表示三轴的推力，Γ是松弛变量。
转化为凸优化问题进行求解，首先引入松弛变量，接着是更改变量，最终实现凸化。其中选择代表推力大小的松弛变量Γ对||T||进行替换，松弛后的约束条件为：
将式(1)线性化得到：
Figure BDA0003064692090000113
式(7)的不等式约束可以写成：
Figure BDA0003064692090000121
其中，γ0(t)＝ln(m0-Tmt/Ispge)，优化性能指标写成：
Figure BDA0003064692090000122
通过以上线性化处理，消除了式(1)中探测器质量m带来的非线性问题。
将时间区间离散化，分成N份，线性化后的动力学方程可以写成
Figure BDA0003064692090000123
的形式。
Figure BDA0003064692090000124
Figure BDA0003064692090000125
引力加速度矢量g是关于探测器相对于小天体质心的位置矢量r的非线性函数，这使得离散后的问题不是一个标准的二阶锥规划问题。采用逐次求解方法，以反复近似方程中的非线性动力学。令第(k-1)个位置矢量r(k-1)代替第k个位置矢量rk，将得到的新位置，作为下一次求解的参考值，反复迭代直到收敛。本文设定的收敛误差为0.5m,可以通过改变迭代次数来减小误差，将主导项放在A会减少迭代次数。
Figure BDA0003064692090000126
通过松弛约束，动力学线性化，引力加速度通过逐次求解来处理，对轨迹优化问题进行了转化，最终变成了以下形式：
Figure BDA0003064692090000131
利用上述凸优化方法进行了4000余组轨迹优化仿真，最终得到了3000多组满足条件的最优解，为后续工作提供训练样本。
。步骤3：设计采用机器学习的小天体表面弹跳运动最优控制推力预测计算模型，实现基于数据统计规律挖掘的小天体表面弹跳运动最优轨迹计算。
利用机器学习的方法对行小天体表面弹跳最优控制推力预测模型进行设计，建立探测器初始状态和推力之间的输入输出模型。
首先在小天体表面一定范围内随机取点作为学习样本，取这些点的坐标作为训练模型的输入向量，由步骤2求出的初末位置处的发动机推力矢量T作为训练模型的输出。xi和yi分别是模型的输入和输出，具体表示如下：
xi＝[r0,rt] 为模型的输入
yi＝Ti 为模型的输出
接着选择适当的均值函数和核函数，设计预测模型。选用零均值函数和平方指数协方差函数，表达式为：
m(x)＝0
Figure BDA0003064692090000132
加入噪声后，k(x,x')为：
Figure BDA0003064692090000133
结合联合分布，可以得到预测数据和已知训练数据的联合正态分布：
Figure BDA0003064692090000141
从而计算出预测数据的均值和方差，分布的均值作为y*的估计值即：
Figure BDA0003064692090000142
训练集中的学习样本输入xi＝[r0i,rti],i＝1,2,…,n,学习样本输出yi＝Ti,i＝1,2…n。检验集中的样本输入为x*＝[r0*,rt*]，输出为y*＝T*。于是,根据上式就能够得到y*的预测值。
通过训练优化模型超参数，即似然函数值最大时，就把模型超参数作为参数估计的输出；似然函数表示为：
Figure BDA0003064692090000143
在求log(p(y|X,θ))的最大值时，对其求导并用最优化方法求得概率最大时的θ值，以使得观测的结果无限接近己知样本；令θ为所有模型超参数的集合，共扼梯度法通过求取训练样本对数似然函数的极大值来获得最优超参数
Figure BDA0003064692090000144
具体如下：
Figure BDA0003064692090000145
最后，利用优化求解出的最优超参数，确定小天体表面弹跳运动推力预测计算模型。使在给定探测任务的某些初始参数时，能够利用预测模型快速、准确的对弹跳探测器最优控制推力进行预测计算。
还包括步骤4：利用步骤3所述的采用机器学习方法的小天体表面弹跳运动最优控制推力预测计算，实现对不同探测位置探测器最优控制推力进行预测计算；
还包括步骤5：利用步骤4预测计算结果对小天体表面弹跳运动的最优轨迹进行计算，进而解决小天体表面探测的相关工程问题；
实施例1
为了检验模型的精度，对预测计算模型进行计算误差分析。随机选取了3000组数据点，从优化算法得到的数据集中抽取800组数据作为检验集，在其余的数据中分别选取200、400、800、1000、1200组数据作为训练集，对模型进行训练。初始推力预测相对误差随训练数据大小变化曲线如图2，末端推力预测相对误差随训练数据大小变化曲线如图3，得到不同训练样本下检验样本的误差分布图如图4。
由图2和图3可以看出，当训练数据达到800组左右时，随着训练数据的增加，预测模型的性能改善并不大，初始推力预测模型误差的平均相对误差稳定在6％附近，末端预测模型误差的平均相对误差稳定在5％附近。
图4横坐标是相对误差，纵坐标是检验点的个数，可以看出推力计算结果误差都主要集中在10％以内，并且主要分布在0％～6％之间。误差小于4％的数量占比在70％以上，而误差小于5％的数量占比达到了80％以上。由此分析表明，本发明的小天体表面弹跳运动最优控制方法精度较高。
实施例2
步骤一到四同实施例1
还包括步骤5：利用步骤4预测计算结果对小天体表面弹跳运动的最优轨迹进行计算，进而解决小天体表面探测的相关工程问题；
为了更充分的检验本发明提出方法的性能，对求得最优轨迹的运算时间进行了分析。对初末时刻施加控制推力，因此对模型训练两次，预测所用的时间如表1所示。并对凸优化方法所用时间和本发明提出的快速预测方法所用时间进行对比，结果如表2所示。
表1求解最优推力所用时间
Figure BDA0003064692090000151
表2求得优化轨迹所用时间对比
Figure BDA0003064692090000161
从表1可以看出机器学习的方法，随着训练样本的增加，计算时间仍然极短，求得推力所用时间都在1s左右。从表2得出随着运动距离的增大，凸优化算法求得优化轨迹所用时间显著变长，且所用时间都大于200s，本发明提出的小天体表面弹跳运动最优控制方法，求出最优轨迹所用时间都在1s量级以内，且随着运动距离的增加，时间变化不大。
以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种小天体表面弹跳运动最优控制方法，其特征在于：包括如下步骤：
步骤1：建立小天体表面动力学模型；
考虑小天体的转动，采用如下动力学模型：
Figure FDA0003484087390000011
其中，r表示探测器相对于小天体质心的位置矢量，ω是旋转角速度矢量，v表示速度矢量，m表示探测器质量，T表示发动机推力矢量，g表示引力加速度矢量，Isp表示比冲，ge表示地球引力加速度常数；
通过牛顿恢复系数来描述碰撞前后法向状态之间的关系；通过引入瞬间摩擦来描述碰撞过程切向关系的变化；则碰撞过程就简化成了碰撞前后速度的变化关系；
步骤2：设计用于得到小天体表面弹跳运动最优控制推力的优化模型，实现小天体表面弹跳运动轨迹的优化求解；
轨迹优化至少包含动力学模型和各种约束条件即初末状态约束及路径约束，在给定动力学模型的情况下，寻找满足各种约束条件的可行解；由于小天体表面探测器携带的燃料有限，降低燃耗对探测任务的开展具有重要的意义，因此以燃耗最优为目标函数进行优化求解；探测器在表面弹跳运动过程中需要满足动力学约束、初始和末端状态约束、路径约束以及推力幅值约束条件，将发动机推力矢量T作为优化变量，进行性能指标的计算；
步骤2具体实现方法为：
步骤2.1：建立小天体表面弹跳运动最优控制推力求解模型，通过优化方法得到小天体表面弹跳运动最优控制推力；
步骤2.1.1：建立小天体表面弹跳运动最优控制推力模型目标函数J；
以燃耗最优为目标函数进行优化；
Figure FDA0003484087390000021
其中，tf表示飞行时间；
步骤2.1.2：确立天体表面弹跳运动最优控制推力模型过程约束；
在弹跳过程中使每次弹跳的路径约束为
||CXi-D1r0||+e1Xi-D2r0≤0 (6)
||CXf-D1rf||+e2Xf-D2rf≤0
Xi和Xf表示探测器的初末状态，各矩阵系数定义为：
Figure FDA0003484087390000022
D2＝[0 0 1]
Figure FDA0003484087390000023
e2＝[0 0 -tanε 0 0 0 0]
r0和rf表示探测器初始位置和目标位置，
Figure FDA0003484087390000024
和ε表示探测器与地面的角度；
推力幅值约束：
||T||≤Tm (7)
Tm表示发动机最大推力；
步骤2.1.3：确立天体表面弹跳运动最优控制推力模型初始和终端约束；
满足的初始状态约束为：
r(t0)＝r0，v(t0)＝v0，m(t0)＝m0 (8)
终端状态约束：
r(tf)＝rf，v(tf)＝vf (9)
选定下一次弹跳的位置，根据碰撞动力学可以求出末端速度vf；v0和m0表示探测器初始速度和初始质量，vf表示末端速度，t0表示初始时间；
步骤2.2：通过优化方法实现小天体表面弹跳最优控制推力的计算；
将发动机推力矢量T作为优化变量，通过优化方法对初末两个时刻的控制推力进行优化，通过插值得到优化的控制推力曲线，再利用步骤1所建立的小天体表面动力学模型进行目标函数J的计算，实现小天体表面弹跳最优控制推力的计算；
步骤2.2所述的优化方法包括遗传算法、凸优化、高斯伪谱法；凸优化问题的局部最优解就是全局最优解，从而保证优化结果的质量，步骤2.2所述的优化方法采用基于凸优化的优化求解器；
步骤2.3：将式(1)的问题转化为凸优化问题进行求解，实现小天体表面弹跳运动轨迹的优化求解；
Figure FDA0003484087390000031
Tx、Ty和Tz分别表示三轴的推力，Γ是松弛变量；
转化为凸优化问题进行求解，首先引入松弛变量，接着是更改变量，最终实现凸化；其中选择代表推力大小的松弛变量Γ对||T||进行替换，松弛后的约束条件为：
将式(1)线性化得到：
Figure FDA0003484087390000032
式(7)的不等式约束可以写成：
Figure FDA0003484087390000033
其中，γ0(t)＝ln(m0-Tmt/Ispge)，优化性能指标写成：
Figure FDA0003484087390000034
通过以上线性化处理，消除了式(1)中探测器质量m带来的非线性问题；
将时间区间离散化，分成N份，线性化后的动力学方程可以写成
Figure FDA0003484087390000035
的形式；
Figure FDA0003484087390000036
X＝[rT,vT,γ]T
Figure FDA0003484087390000037
引力加速度矢量g是关于探测器相对于小天体质心的位置矢量r的非线性函数，这使得离散后的问题不是一个标准的二阶锥规划问题；采用逐次求解方法，以反复近似方程中的非线性动力学；令第(k-1)个位置矢量r(k-1)代替第k个位置矢量rk，将得到的新位置，作为下一次求解的参考值，反复迭代直到收敛；本文设定的收敛误差为0.5m,可以通过改变迭代次数来减小误差，将主导项放在A会减少迭代次数；
Figure FDA0003484087390000041
通过松弛约束，动力学线性化，引力加速度通过逐次求解来处理，对轨迹优化问题进行了转化，最终变成了以下形式，实现小天体表面弹跳运动轨迹的优化求解；
Figure FDA0003484087390000042
Figure FDA0003484087390000043
步骤3：设计采用机器学习的小天体表面弹跳运动最优控制推力预测计算模型，实现基于数据统计规律挖掘的小天体表面弹跳运动最优轨迹计算；
利用机器学习的方法对行小天体表面弹跳最优控制推力预测模型进行设计，建立探测器初始状态和推力之间的输入输出模型；
首先在小天体表面一定范围内随机取点作为学习样本，取这些点的坐标作为训练模型的输入向量，由步骤2求出的初末位置处的发动机推力矢量T作为训练模型的输出； xi和yi分别是模型的输入和输出，具体表示如下：
xi＝[r0,rt]为模型的输入
yi＝Ti为模型的输出
接着选择适当的均值函数和核函数，设计预测模型；选用零均值函数和平方指数协方差函数，表达式为：
m(x)＝0
Figure FDA0003484087390000051
加入噪声后，k(x,x')为：
Figure FDA0003484087390000052
结合联合分布，可以得到预测数据和已知训练数据的联合正态分布：
Figure FDA0003484087390000053
从而计算出预测数据的均值和方差，分布的均值作为y*的估计值即：
Figure FDA0003484087390000054
训练集中的学习样本输入xi＝[r0i,rti],i＝1,2,…,n,学习样本输出yi＝Ti,i＝1,2…n；检验集中的样本输入为x*＝[r0*,rt*]，输出为y*＝T*；于是,根据上式就能够得到y*的预测值；
通过训练优化模型超参数，即似然函数值最大时，就把模型超参数作为参数估计的输出；似然函数表示为：
Figure FDA0003484087390000055
在求log(p(y|X,θ))的最大值时，对其求导并用最优化方法求得概率最大时的θ值，以使得观测的结果无限接近己知样本；令θ为所有模型超参数的集合，共扼梯度法通过求取训练样本对数似然函数的极大值来获得最优超参数
Figure FDA0003484087390000056
具体如下：
Figure FDA0003484087390000057
最后，利用优化求解出的最优超参数，确定小天体表面弹跳运动推力预测计算模型；使在给定探测任务的某些初始参数时，能够利用预测模型快速、准确的对弹跳探测器最优控制推力进行预测计算。
2.如权利要求1所述的一种小天体表面弹跳运动最优控制方法，其特征在于：还包括步骤4：利用步骤3所述的采用机器学习方法的小天体表面弹跳运动最优控制推力预测计算，实现对不同探测位置探测器最优控制推力进行预测计算。
3.如权利要求2所述的一种小天体表面弹跳运动最优控制方法，其特征在于：还包括步骤4：还包括步骤5：利用步骤4预测计算结果对小天体表面弹跳运动的最优轨迹进行计算，进而解决小天体表面探测的相关工程问题。
4.如权利要求1所述的一种小天体表面弹跳运动最优控制方法，其特征在于：步骤一所述通过Newton恢复系数来描述碰撞前后法向状态之间的关系，求得碰撞后法向，切向速度的方法为：
牛顿恢复系数表示为：
Figure FDA0003484087390000061
其中v+和v-分别表示碰撞后和碰撞前物体的相对法向速度；
碰撞后法向速度是
vn1＝-evn0
探测器会因为与小天体表面碰撞而发生轻微形变，因此定义恢复系数e；其中vn0和vn1是碰撞前后的法向速度，在切向方向，碰撞后的速度表示为
vt1＝vt0-μ(1-e)vn0
其中μ是小天体表面摩擦系数，vt0和vt1是碰撞前后的切向速度。
CN202110522720.4A 2021-05-13 2021-05-13 一种小天体表面弹跳运动最优控制方法 Expired - Fee Related CN113110559B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110522720.4A CN113110559B (zh) 2021-05-13 2021-05-13 一种小天体表面弹跳运动最优控制方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110522720.4A CN113110559B (zh) 2021-05-13 2021-05-13 一种小天体表面弹跳运动最优控制方法

Publications (2)

Publication Number Publication Date
CN113110559A CN113110559A (zh) 2021-07-13
CN113110559B true CN113110559B (zh) 2022-03-18

Family

ID=76723135

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110522720.4A Expired - Fee Related CN113110559B (zh) 2021-05-13 2021-05-13 一种小天体表面弹跳运动最优控制方法

Country Status (1)

Country Link
CN (1) CN113110559B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114371734B (zh) * 2022-01-07 2023-07-04 中国人民解放军63921部队 基于高斯伪谱法的轨迹优化方法、电子设备和存储介质

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105739537A (zh) * 2016-03-29 2016-07-06 北京理工大学 一种小天体表面附着运动主动控制方法
CN106155076A (zh) * 2016-08-23 2016-11-23 华南理工大学 一种多旋翼无人飞行器的稳定飞行控制方法
CN106778012A (zh) * 2016-12-29 2017-05-31 北京理工大学 一种小天体附着探测下降轨迹优化方法
RU2654071C1 (ru) * 2017-07-24 2018-05-16 Акционерное общество "Научно-производственное объединение "Центральный научно-исследовательский институт технологии машиностроения" (АО "НПО "ЦНИИТМАШ") Способ прогнозирования ресурсоспособности стали для корпусов реакторов типа ввэр
CN108388135A (zh) * 2018-03-30 2018-08-10 上海交通大学 一种基于凸优化的火星着陆轨迹优化控制方法
CN109145490A (zh) * 2018-09-10 2019-01-04 北京理工大学 基于数据特征挖掘的行星进入可达集最优子集计算方法
CN109669471A (zh) * 2018-12-17 2019-04-23 北京理工大学 小天体悬停姿轨耦合自抗扰控制方法
CN110244344A (zh) * 2019-06-05 2019-09-17 中南大学 一种基于深度学习的tbm超前地质预报方法
CN112269390A (zh) * 2020-10-15 2021-01-26 北京理工大学 考虑弹跳的小天体表面定点附着轨迹规划方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7177785B2 (en) * 2003-08-22 2007-02-13 Honeywell International, Inc. Systems and methods for improved aircraft performance predictions
CN110736470B (zh) * 2019-11-06 2021-04-20 北京理工大学 一种不规则形状小天体附近连续推力轨道混合搜索方法

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105739537A (zh) * 2016-03-29 2016-07-06 北京理工大学 一种小天体表面附着运动主动控制方法
CN106155076A (zh) * 2016-08-23 2016-11-23 华南理工大学 一种多旋翼无人飞行器的稳定飞行控制方法
CN106778012A (zh) * 2016-12-29 2017-05-31 北京理工大学 一种小天体附着探测下降轨迹优化方法
RU2654071C1 (ru) * 2017-07-24 2018-05-16 Акционерное общество "Научно-производственное объединение "Центральный научно-исследовательский институт технологии машиностроения" (АО "НПО "ЦНИИТМАШ") Способ прогнозирования ресурсоспособности стали для корпусов реакторов типа ввэр
CN108388135A (zh) * 2018-03-30 2018-08-10 上海交通大学 一种基于凸优化的火星着陆轨迹优化控制方法
CN109145490A (zh) * 2018-09-10 2019-01-04 北京理工大学 基于数据特征挖掘的行星进入可达集最优子集计算方法
CN109669471A (zh) * 2018-12-17 2019-04-23 北京理工大学 小天体悬停姿轨耦合自抗扰控制方法
CN110244344A (zh) * 2019-06-05 2019-09-17 中南大学 一种基于深度学习的tbm超前地质预报方法
CN112269390A (zh) * 2020-10-15 2021-01-26 北京理工大学 考虑弹跳的小天体表面定点附着轨迹规划方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Predictor-Corrector Strategy Based Energy Suboptimal Obstacle Avoidance for Landing on Small Bodies;Jiateng Long;《2018 AIAA Guidance, Navigation, and Control Conference》;20180112;全文 *

Also Published As

Publication number Publication date
CN113110559A (zh) 2021-07-13

Similar Documents

Publication Publication Date Title
Kang et al. Deep convolutional identifier for dynamic modeling and adaptive control of unmanned helicopter
Kim et al. Estimating the non-linear dynamics of free-flying objects
CN105739537B (zh) 一种小天体表面附着运动主动控制方法
Gaudet et al. Adaptive pinpoint and fuel efficient mars landing using reinforcement learning
Lin et al. A gated recurrent unit-based particle filter for unmanned underwater vehicle state estimation
CN108959182B (zh) 基于高斯过程回归的小天体引力场建模方法
CN113110559B (zh) 一种小天体表面弹跳运动最优控制方法
CN115437406A (zh) 基于强化学习算法的飞行器再入跟踪制导方法
CN109543284B (zh) 基于Kriging空间插值的火星大气进入段最优制导方法
Zou et al. CNN based adaptive Kalman filter in high-dynamic condition for low-cost navigation system on highspeed UAV
CN112896560B (zh) 小天体表面安全弹跳移动轨迹规划方法
Zheng et al. Constrained trajectory optimization with flexible final time for autonomous vehicles
CN112560343B (zh) 基于深度神经网络与打靶算法的J2摄动Lambert问题求解方法
CN110231619B (zh) 基于恩克法的雷达交接时刻预报方法及装置
Hao et al. A novel framework for reliability assessment of payload fairing separation considering multi-source uncertainties and multiple failure modes
CN104715133B (zh) 一种待辨识对象的运动学参数在轨辨识方法和装置
He et al. Active compliance control of a position-controlled industrial robot for simulating space operations
US20220245310A1 (en) System and method for modeling rigid body motions with contacts and collisions
CN116068894A (zh) 基于双层强化学习的火箭回收制导方法
Yuqi et al. Time-varying parameters estimation with adaptive neural network EKF for missile-dual control system
CN109709805B (zh) 一种考虑不确定性因素的航天器鲁棒交会轨迹设计方法
El-Fakdi et al. Autonomous underwater vehicle control using reinforcement learning policy search methods
CN109583007B (zh) 一种火星进入飞行状态不确定性量化方法
CN112161626A (zh) 一种基于航路跟踪映射网络的高可飞性航路规划方法
Menon et al. Worst-case analysis of flight control laws for re-entry vehicles

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20220318