基于强化学习算法的飞行器再入跟踪制导方法
1.本发明属于飞行器制导控制技术领域，涉及一种飞行器再入跟踪制导方法，可用于火箭回收。
背景技术：
[0002][0003]
作为飞行器总体与控制系统的核心技术，轨迹与制导设计对于飞行器来说至关重要，且二者逐渐呈现出融合趋势。轨迹优化技术的研究早期来自协调总体参数设计、最大化资源利用率、降低任务运营成本等顶层设计需求。传统制导设计的主要任务是通过离线设计固定的参考轨迹并在线跟踪，确保完成飞行任务。由于传统的跟踪制导方案是建立在简化运动模型基础上，因而存在的对模型依赖性较大和适应性较差的问题，当建模偏差、外部扰动较大时，制导性能会严重恶化。例如“”号再入过程的标称轨迹制导方案。该方案的飞行器再入飞行环境非常复杂，其要历经外太空、稀薄大气层以及稠密大气层，由于在此过程中机械能迅速转化为热能，机体会剧烈升温，同时面临着巨大动压及过载的严苛考验。且由于再入过程呈现出多约束、强耦合、快时变、强非线性、强不确定性、大包络、长航时的特点，为制导控制系统设计带来了极大挑战，极大地制约了行器应有的飞行性能。
技术实现要素：
[0004]
本发明的目的在于针对上述现有技术存在的缺陷，提出一种基于强化学习算法的飞行器再入跟踪制导方法，以根据实时飞行状态快速生成攻角增量指令，通过跟踪期望的高度曲线，有效增强飞行器的任务执行能力，减小对飞行器模型的依赖、提高适应性及制导性能。
[0005]
为实现上述目的，本发明的技术方案包括如下：
[0006]
(1)将飞行器再入段轨迹优化描述为由数学模型、边界条件、容许控制、性能指标、过程约束构成的连续最优控制问题p0；
[0007]
(2)对p0进行更换形式、松弛变量、软化约束、逐次线性化方法的凸化处理，得到序列凸最优控制问题p1，采用伪谱法对该p1进行离散参数化处理，得到序列二阶锥规划问题p2；
[0008]
(3)采用内点法求解序列二阶锥规划问题p2，得最优参考轨迹；
[0009]
(4)对最优参考轨迹进行采样得到参考轨迹训练数据集；
[0010]
(5)构建神经网络actor网络：
[0011]
(5a)建立由第一输入层、第一隐藏层、第一输出层依次连接组成的动作评估子网络 actor_eval，该第一输入层输入的状态变量为当前状态智能体的地心距、经度、纬度、航迹角、航向角这五个状态变量，该第一输出层输出信息为当前状态的攻角增量指令；
[0012]
(5b)建立由第二输入层、第二隐藏层、第二输出层依次连接组成的动作目标子网络actor_target，该第二输入层输入的状态变量为智能体下一个状态的地心距、经度、纬
度、航迹角、航向角这五个状态变量，该第二输出层输出信息为下一个状态的目标攻角增量指令；
[0013]
(5c)将动作评估子网络actor_eval网络和动作目标子网络actor_target网络并联，构成actor网络，用于接收经验回放池的状态信息，并输出攻角增量指令信息；
[0014]
(6)构建神经网络critic网络：
[0015]
(6a)建立由第三输入层、第三隐藏层、第三输出层依次连接组成的价值评估子网络 critic_eval，该第三输入层输入的变量为当前状态智能体的地心距、经度、纬度、航迹角、航向角这五个状态变量以及actor_eval输出的攻角增量指令，该第三输出层输出信息为智能体这一状态采取指令后产生的累计奖励；
[0016]
(6b)建立由第四输入层、第四隐藏层、第四输出层依次连接组成的价值目标子网络 critic_target，该第四输入层输入的变量为下一个状态智能体的地心距、经度、纬度、航迹角、航向角这五个状态变量以及actor_target输出的攻角增量指令，该第四输出层输出信息为智能体下一状态采取指令后产生的目标累计奖励；
[0017]
(6c)将价值评估子网络critic_eval网络和价值目标子网络critic_target网络并联，构成critic网络，用于接受经验回放池的状态信息和actor网络输出的攻角增量指令信息，并输出智能体采取指令后产生的累计奖励信息；
[0018]
(7)设计actor网络和critic网络的奖励函数r均为：
[0019]
r＝r1+r2，其中r1为目标奖励，r2为约束奖励；
[0020]
(8)利用(4)中所得到的参考轨迹训练数据集使用深度确定性策略梯度算法ddpg，对actor网络和critic网络进行并行离线训练，当智能体的累计奖励收敛于一个最大值时，得到训练好的制导网络；
[0021]
(9)在线获取飞行器再入段制导指令并实时跟踪：
[0022]
飞行器在飞行期间，机载计算机读取飞行器导航系统测得的实时飞行状态量，将实时飞行的状态量作为制导网络的输入进行前向传播，得到实时的制导指令；飞行器根据指令将产生攻角变化实时跟踪预定高度曲线。
[0023]
本发明与现有技术相比，具有如下优点：
[0024]
1)本发明由于构建了神经网络模型，并将大量优化计算转移到离线训练过程中，减小了在线运算量，提高了生成制导指令的速度，使得飞行器能够实时对环境做出反应。
[0025]
2)本发明由于设计了关于目标高度和攻角增量的奖励函数，提高了制导精度，使得飞行器能够适应复杂的环境并实现对高度的精细跟踪。
附图说明
[0026]
图1为本发明的实现流程图；
[0027]
图2为本发明中在构建飞行器最优轨迹中使用的原始非凸集；
[0028]
图3为本发明中在凸化飞行器最优轨迹中得到的凸容许控制集；
[0029]
图4为本发明中求解飞行器最优轨迹得到的再入段倾侧角的参考轨迹；
[0030]
图5为用本发明对气动参数正向固定拉偏的累计收益曲线图；
[0031]
图6为用本发明对气动参数正向固定拉偏的跟踪效果图；
[0032]
图7为用本发明对气动参数正向随机拉偏的跟踪效果图；
[0033]
图8为用本发明对气动参数正向随机拉偏时的奖励曲线图；
[0034]
图9为用本发明对气动参数反向固定拉偏的累计收益曲线图；
[0035]
图10为用本发明对气动参数反向固定拉偏的跟踪效果图；
[0036]
图11为用本发明对气动参数反向随机拉偏的跟踪效果图；
[0037]
图12为用本发明对气动参数反向随机拉偏的奖励曲线图。
具体实施方式
[0038]
以下结合附图对本发明具体实施例和效果作进一步详细说明。
[0039]
参照图1，本实例的实现步骤如下：
[0040]
步骤1，建立飞行器运动模型。
[0041]
飞行器质心运动方程是研究其运动特性的基础，本实例以无动力可重复使用运载器 rlv为研究对象，在半速度系下建立其再入质心运动方程。轨迹优化问题通常不考虑控制力作用，并且由于附加哥氏力量级较小，忽略其影响，rlv再入过程主要依靠空气动力和地球引力改变飞行状态，假设地球为理想圆球，在半速度坐标系下，建立rlv的状态微分方程f1如下：
[0042][0043]
其中，r为rlv无量纲地心距，为rlv无量纲速度，ω为地球自转角速率，r0为地球半径，γ为航迹角，ψ为航向角，σ为倾侧角，l＝0.5r0ρv2sc
l
/m为升力加速度，d＝0.5r0ρv2scd/m为阻力加速度；e为rlv的无量纲能量，cd与c
l
分别为阻力系数与升力系数，这两个系数均为攻角α和马赫数ma的函数，s为rlv参考面积，ρ＝ρ0exp(-(r-1)/h)表示大气密度，其是地心距r的函数，ρ0为海平面大气密度，h为无量纲标高。
[0044]
对于状态微分方程，本实施例选择状态量为s＝(r,θ,φ,γ,ψ)，控制量为倾侧角u＝σ。
[0045]
步骤2，建立飞行器约束条件和性能指标。
[0046]
2.1)建立飞行器的约束条件：
[0047]
2.1.1)建立过程约束k：
[0048]
由于无动力可重复使用运载器rlv在大气层内飞行时，机体与大气摩擦产生大量热量，考虑到rlv的结构强度和设备安全工作条件，出于飞行安全，对动压q、热流率和法向过载n进行约束，记为k：
[0049][0050]
其中，kq为热流率计算系数，g0为海平面重力加速度，q
max
、和n
max
分别为rlv 可承受的最大动压、最大热流率和最大过载；
[0051]
2.1.2)建立初始状态约束l和终端状态约束φ1：
[0052]
根据rlv的飞行任务可以确定rlv的初始状态信息和终端状态信息，进而对rlv的初始状态和终端状态施加约束：
[0053]
初始状态约束l：
[0054][0055]
其中，x表示飞行器的五个状态变量，表示飞行器的无量纲初状态，e0表示rlv 的初始无量纲能量；
[0056]
终端状态约束φ1：
[0057][0058]
其中，ef表示rlv的终端无量纲能量，分别表示rlv终端无量纲地心距、经度、纬度、航迹角和航向角的设定值；
[0059]
2.1.3)建立容许控制c：
[0060]
在rlv飞行过程中，为确保rlv不会进行大角度的翻转，还应当对控制量倾侧角加以限制，因此建立容许控制约束c来约束倾侧角的幅值；
[0061]
容许控制c：
[0062]
c：σ
min
≤|σ|≤σ
max
[0063]
其中，σ
min
和σ
max
分别表示倾侧角的下限和上限。
[0064]
2.2)建立性能指标j
[0065]
在飞行器再入的轨迹优化中，通常对飞行器的期望要求包含飞行器再入时间最短，终端速度最大，飞行航程最长，消耗能量最小，本实例中性能指标j为再入时间最短，即：
[0066][0067]
步骤3，构建rlv再入段轨迹优化模型的最优控制问题p0。
[0068]
根据上述状态微分方程f1、过程约束k、初始约束l、终端约束φ1、容许控制c和性能指标j，将rlv再入段的轨迹优化模型描述为包括这些参数的最优控制问题p0：
[0069][0070]
步骤4，将飞行器再入连续最优控制问题p0转化为序列凸最优控制问题p1。
[0071]
4.1)过程约束k的形式转换：
[0072]
将过程约束k重新描述为关于状态变量r的线性不等式形式，满足凸问题的要求：
[0073][0074]
上式等价于
[0075][0076]
其中，lq(e)、lq(e)、ln(e)分别表示飞行器的热流率、动压和过载关于无量纲地心距r 的函数，表示lq(e)、lq(e)、ln(e)三者中最大值；。
[0077]
4.2)控制约束的松弛：
[0078]
4.2.1)转换控制变量
[0079]
由于选择倾侧角σ为控制量时，在对微分方程组进行逐次线性化时易产生不良震荡，因此需要更换控制变量来消除这种不良震荡，详细说明如下：
[0080]
首先，根据选择倾侧角σ为控制量时，质心运动方程f1关于状态量x和控制量u＝σ都是非凸的性质，采用逐次线性化方法凸化方程f1，逐次求解过程中第k次的解包含状态变量解x
(k)
(e)和控制量解u
(k)
(e)，记为{x
(k)
(e)；u
(k)
(e)}；
[0081]
其次，根据第k次迭代得到的解，将第k+1次迭代将方程线性化为：
[0082]
x
′
＝f
x
(x
(k)
,u
(k)
,e)x+fu(x
(k)
,u
(k)
,e)u+b(x
(k)
,u
(k)
,e)
[0083]
其中，x
′
表示五个状态变量对无量纲能量e的导数，表示微分方程组f对x的雅可比矩阵、表示微分方程组f对u的雅可比矩阵，b(x
(k)
,u
(k)
,e)＝f(x
(k)
,u
(k)
,e)-f
x
(x
(k)
,u
(k)
,e)x
(k)-fu(x
(k)
,u
(k)
,e)u
(k)
表示余矢量。
[0084]
接着，根据上述系数矩阵f
x
、fu以及余矢量b与u
(k)
有关，在迭代过程中控制量u
(k)
易产生震荡的情况，即在迭代过程中通过控制量u
(k)
的振荡控制第k+1次迭代得到的控制量 u
(k+1)
，随着迭代不断进行，振荡可能会进一步放大。这些振荡不是问题的解所固有的，是数值求解过程不稳定的表现，为消除线性化动力学方程与u
(k)
的关系，引入如下控制变量替换原控制量：其中，u1表示倾侧角余弦值，u2表示倾侧角正弦值，将替换后的控制变量从质心运动方程中分离出来，依据新的控制变量对微分方程组进行线性化，其雅可比系数矩阵近似为0矩阵，可消除不良的震荡传递。
[0085]
4.2.2)建立新的控制变量约束
[0086]
控制变量转换后，控制约束表示为：cosσ
max
≤u1≤cosσ
min
，由于控制量u1和u2不是相互独立的，还需要满足如下等式：
[0087][0088]
该等式约束是一个二次约束，约束确定的原始容许控制集如图2所示，该原始容许控制集是单位圆上在u1＝ω
l
＝cosσ
max
和u1＝ωh＝cosσ
min
之间的两段圆弧，这两段圆弧构成的集合是非凸的，需要进行凸化处理，即对约束进行松弛处理，将原约束松弛为如下二阶锥不等式约束：
[0089][0090]
通过松弛约束后确定的集合扩大了原始容许控制集，使得原非凸集合变为凸集，如图 3阴影区域所示。
[0091]
4.3)性能指标的逐次线性化：
[0092]
4.3.1)将阻力加速度表达式代入性能指标，将性能指标表示为如下状态量的函数：
[0093]
[0094]
其中，，d＝0.5r0ρv2scd/m为阻力加速度，系数不含任何状态量，为rlv无量纲速度，m为飞行器质量，ω为地球自转角速率，r0为地球半径，e为rlv的无量纲能量，cd为阻力系数，其为攻角α和马赫数ma的函数，s为rlv 参考面积，ρ＝ρ0exp(-(r-1)/h)表示大气密度，其是地心距r的函数，ρ0为海平面大气密度，h为无量纲标高。
[0095]
4.3.2)线性化处理性能指标
[0096]
性能指标式中的大气密度ρ是唯一可优化的量，相关状态分量为地心距r。大气密度近似采用指数模型，该性能指标的被积函数是关于r的凸函数，但为最终得到二阶锥规划问题，需将其转化为关于r的线性函数，为此，在前一次迭代的解r
(k)
(e)附近，对该性能指标的被积函数进行一阶泰勒展开：
[0097][0098]
其中，ρ
(k)
表示利用r
(k)
(e)计算得到的大气密度，表示大气密度ρ对r的导数在r
(k)
(e)处的值，分别表示在前一次迭代的解r
(k)
(e)附近性能指标关于无量纲地心距r的一阶项和余项；
[0099]
4.3.3)引入正则化项确保约束不变
[0100]
在线性化处理原性能指标后，还需引入正则化项，确保步骤4.2)中松弛约束后产生积极效果，即在线性化后的性能指标中引入
[0101]
其中，ε
ψ
≠0是一个常量，为确保得到的解近似于原问题的解，需要将ε
ψ
的量级设置的足够小。
[0102]
4.4)终端约束的松弛：
[0103]
终端约束虽然均为线性等式约束，但在初期寻优过程中，由于经、纬度两个硬约束可能难以满足，为此，在性能指标中引入如下惩罚项替代这两个终端约束：
[0104][0105]
其中，d(θ(ef),φ(ef))表示惩罚项，c
θ
＞0和c
φ
＞0为给定常量，
[0106]
由于该惩罚项是非线性的，故需要再通过引入两个松弛变量和将其转化为一个线性函数并入性能指标中：且受到线性不等式约束：
[0107]
将性能指标更新为：
[0108]
[0109]
其中，c1、c0分别表示性能指标逐次线性化过程中关于无量纲地心距r的一阶项和余项，余项，分别表示rlv终端无量纲地心距、经度、纬度、航迹角和航向角的设定值。
[0110]
4.5)质心运动方程的逐次线性化：
[0111]
根据更换后的控制量，质心运动方程可表示为：
[0112][0113]
其中，f1(x,e)为原微分方程组中不含控制量的项系数矩阵，b(x,e)为原微分方程组中含有控制量的项系数矩阵；
[0114]
采用逐次线性化方法凸化质心运动方程，逐次求解过程中第k次的解包含状态变量解 x
(k)
(e)和控制量解u
(k)
(e)，记为{x
(k)
(e)；u
(k)
(e)}；
[0115]
根据第k次迭代得到的解，将第k+1次迭代将方程线性化为：
[0116]
x
′
＝f
x
(x
(k)
,u
(k)
,e)x+b(x
(k)
,e)u+b(x
(k)
,u
(k)
,e)
[0117]
其中，表示雅可比矩阵；
[0118]
b(x
(k
),u
(k)
,e)＝f1(x
(k)
,e)-f
x
(x
(k)
,u
(k)
,e)x
(k)
表示余矢量；
[0119]
表示 b(x,e)u对x的偏导数，表示f1(x,e)对x的偏导数；
[0120]
l/d＝c
l
/cd，其与v均可看作能量e的单变量函数，因此矩阵b(x,e)中与状态量x相关的元素只有cosγ，导致中存在一非零项，由于飞行器再入过程中，航迹角|γ|很小，即该非零项很小，因此线性化方程可简化为：
[0121]
x
′
＝f
x
(x
(k)
,e)x+b(x
(k)
,e)u+b(x
(k)
,e)；
[0122]
进一步，从质心运动方程中分离出量级很小的地球自转相关项，将质心运动方程更新为如下简洁形式：
[0123]
x
′
＝f1(x,e)+b(x,e)u＝f0(x,e)+f
ω
(x,e)+b(x,e)u
[0124]
其中，表示微分方程中不含控制量的项，f
ω
(x,e) 为地球自转角速率相关项；
[0125]
再进一步，将更新过的质心运动方程进行逐次线性化处理，并省略符号中的无量纲能量e，根据第k次迭代得到的解，将第k+1次迭代将方程线性化为：
[0126]
x
′
＝a(x
(k)
)x+b(x
(k)
)u+b(x
(k)
)
[0127]
其中，b(x
(k)
)＝f0(x
(k)
)-a(x
(k)
)x
(k)
+f
ω
(x
(k)
)表示余矢量，a(x
(k)
)和b(x
(k)
)分别表示f0(x
(k)
)对x和u的雅可比矩阵，f0(x
(k)
)表示f0(x,e)第k次迭代得到的状态微分方程， f
ω
(x
(k)
)表示f
ω
(x,e)第k次迭代得到的状态微分方程；
[0128]
为确保逐次线性化的局部有效性，施加信赖域约束：|x-x
(k)
|≤δ，其中，δ为给定五维常矢量。
[0129]
4.6)根据步骤4.1)到步骤4.5)中对约束和性能指标的凸化处理，将rlv再入段的轨迹优化模型重新描述为序列凸最优控制问题p1：
[0130][0131]
步骤5，将序列凸最优控制问题p1转化为序列二阶锥规划问题p2。
[0132]
根据gauss伪谱法的原理，将序列凸最优控制问题p1离散参数化为序列凸规划问题 p2，其实现步骤如下：
[0133]
5.1)设置配点和离散点：
[0134]
将能量自变量e转换为伪能量自变量ε∈[-1,1]：
[0135][0136]
其中，e0、ef分别表示飞行器的初始无量纲能量和终端无量纲能量；
[0137]
在ε∈(-1,1)内设置n个lg点(ε1,ε2,
…
,εm,
…
,εn)作为配点，每个配点的坐
标是n 阶legendre多项式的每一个根，其中m＝1,2,
…
,n表示每一个lg点；
[0138]
设置离散点，包括伪能量自变量的两端点ε0、εf和所有区间内部的lg点，即共计 n+2个离散点ε0,ε1,
…
,εn,
…
,ε
n+1
，其中n＝0,1,
…
,n+1，表示每一个离散点，ε
n+1
＝εf表示伪能量自变量的右端点。
[0139]
5.2)质心运动微分方程的离散：
[0140]
5.2.1)利用拉格朗日插值多项式表示伪能量自变量ε∈[-1,1)的质心运动微分方程约束：
[0141]
以n+1个拉格朗日插值多项式作为基函数近似每个lg点的状态变量其中i＝0,1,
…
,n表示每一个插值多项式；
[0142]
利用近似的状态变量x(εm)对伪能量自变量
ε
求导，即其中，表示基函数li(εm)在lg点处的导数；
[0143]
利用代替原质心运动微分方程的x
′
，将质心运动微分方程约束转化为lg配点处的代数方程约束：
[0144][0145]
其中，表示状态变量在第m个lg点处第k次迭代的解，xm表示状态变量在第m个 lg点处的解，um表示控制量在第m个lg点处的解；
[0146]
5.2.2)利用gauss积分估计伪能量自变量右端点质心运动微分方程约束：
[0147]
由于状态变量近似表达式对应的伪能量区间为[-1,1)，未包含终端状态变量x
n+1
，因此采用gauss积分来估计x
n+1
，即终端状态x
n+1
需满足的质心运动微分方程约束：
[0148][0149]
其中，wm表示lg积分权重。
[0150]
5.3)设置离散后的其他约束条件：
[0151]
5.3.1)信赖域约束：
[0152][0153]
其中，表示第n个离散点处第k次迭代的状态变量解；
[0154]
5.3.2)控制约束：
[0155][0156]
其中，(u1)m和(u2)m分别表示控制量u1和u2在第m个lg点的值，ω
l
和ωh分别代表倾侧角余弦值的最小值和最大值；
[0157]
5.3.3)过程约束：
[0158][0159]
其中，rm表示飞行器的无量纲地心距在第m个lg点的值，表示飞行器热流率、动压和过载在第m个lg点的最大值；
[0160]
5.3.4)边界条件：
[0161][0162][0163][0164][0165][0166][0167]
其中，表示飞行器状态变量设定的初始值，r
n+1
表示飞行器无量纲地心距的终端状态，θ
n+1
表示飞行器经度的终端状态，φ
n+1
表示飞行器纬度的终端状态，γ
n+1
表示飞行器航迹角的终端状态，ψ
n+1
表示飞行器航向角的终端状态，分别表示飞行器终端无量纲地心距、经度、纬度、航迹角和航向角，和分别表示经度和纬度的松弛变量；
[0168]
5.3.5)性能指标：
[0169]
性能指标函数通过gauss积分近似为：
[0170][0171]
其中，c
θ
、分别表示经度和纬度的惩罚系数，分别表示性能指标逐次线性化过程中关于无量纲地心距r的一阶项和余项，ε
ψ
表示航向角正则化项系数，ψm表示航向角在第m个lg点上的值。
[0172]
5.4)根据步骤5.1)到步骤5.3)中对各约束和性能指标的离散化处理，将rlv再入段的轨迹优化模型重新描述为序列凸最优控制问题p2：
[0173][0174]
其中，n＝0,1,2...n+1表示每一个离散点，m＝1,2,3...n表示每一个lg点。
[0175]
步骤6，对无动力可重复使用运载器rlv再入段序列二阶锥规划问题p2进行离线求解，其流程如下：
[0176]
6.1)令迭代计数k＝0，给定离散点处的初始状态量6.1)令迭代计数k＝0，给定离散点处的初始状态量表示第n个离散点处第k次迭代的状态变量解，n＝0,1,2...n+1；
[0177]
6.2)将迭代计数k加一，采用内点法求解问题p2，得到解其中，z
(k+1)
表示第k+1次迭代的解向量，
[0178]
表示第k+1次迭代的状态变量解，表示第n个离散点处第k+1次迭代的状态变量解，n＝0,1,2...n+1，
[0179]
表示第k+1次迭代的控制量解，表示第m个离散点处第k+1次迭代的控制量解，m＝1,2,3...n，
[0180]
和表示第k+1次迭代的松弛变量解；
[0181]
6.3)设一个五维列向量ε表示收敛判断条件，判断每一个离散点处的状态变量解是否满足的收敛条件：
[0182]
若不满足，则令x
(k)
＝x
(k+1)
，返回步骤6.2；
[0183]
否则，输出最优解计算结束。
[0184]
最终，根据最优解中的u
(k+1)
得到再入段控制量倾侧角的参考轨迹，如图4所示。
[0185]
步骤7，对步骤6所得最优参考轨迹进行采样得到参考轨迹训练数据集。
[0186]
步骤8，构建神经网络actor和神经网络critic。
[0187]
8.1)构建actor网络：
[0188]
8.1.1)建立由第一输入层、第一隐藏层、第一输出层依次连接组成的动作评估子网络 actor_eval，该第一输入层输入的状态变量为当前状态智能体的地心距、经度、纬度、
航迹角、航向角这五个状态变量，该第一输出层输出信息为当前状态的攻角增量指令；
[0189]
8.1.2)建立由第二输入层、第二隐藏层、第二输出层依次连接组成的动作目标子网络 actor_target，该第二输入层输入的状态变量为智能体下一个状态的地心距、经度、纬度、航迹角、航向角这五个状态变量，该第二输出层输出信息为下一个状态的目标攻角增量指令；
[0190]
8.1.3)将动作评估子网络actor_eval网络和动作目标子网络actor_target网络并联，构成actor网络，用于接收经验回放池的状态信息，并输出攻角增量指令信息；
[0191]
8.2)构建critic网络：
[0192]
8.2.1)建立由第三输入层、第三隐藏层、第三输出层依次连接组成的价值评估子网络 critic_eval，该第三输入层输入的变量为当前状态智能体的地心距、经度、纬度、航迹角、航向角这五个状态变量以及动作评估子网络actor_eval输出的攻角增量指令，该第三输出层输出信息为智能体这一状态采取指令后产生的累计奖励；
[0193]
8.2.2)建立由第四输入层、第四隐藏层、第四输出层依次连接组成的价值目标子网络 critic_target，该第四输入层输入的变量为下一个状态智能体的地心距、经度、纬度、航迹角、航向角这五个状态变量以及动作目标子网络actor_target输出的攻角增量指令，该第四输出层输出信息为智能体下一状态采取指令后产生的目标累计奖励；
[0194]
8.2.3)将价值评估子网络critic_eval网络和价值目标子网络critic_target网络并联，构成critic网络，用于接受经验回放池的状态信息和actor网络输出的攻角增量指令信息，并输出智能体采取指令后产生的累计奖励信息；
[0195]
8.3)将actor网络和critic网络中的所有隐藏层均设置为100个神经元，且隐藏层的激活函数均使用relu函数，输出层均使用tanh函数。
[0196]
步骤9，设置神经网络actor和神经网络critic的奖励函数。
[0197]
奖励函数包括目标奖励和约束奖励，其设置如下：
[0198]
9.1)根据任务目标，设置目标奖励r1：
[0199][0200]
其中，设h
next
表示执行完动作δα后下一时间步的高度，δα表示攻角增量指令，h0表示下一时间步对应标称高度；
[0201]
9.2)根据攻角增量大小，设置约束奖励r2：
[0202]
r2＝-10|δα|
[0203]
9.3)将目标奖励r1与约束奖励r2相加，得到神经网络actor和神经网络critic的奖励函数为：r＝r1+r2。
[0204]
步骤10，对神经网络actor和神经网络critic进行同步离线训练，获得训练好的制导网络。
[0205]
10.1)参数初始化：
[0206]
随机初始化actor_eval网络参数θ
μ
、actor_target网络参数θ
μ
′
；
[0207]
随机初始化critic_eval网络参数θq、critic_target网络参数θq′
；
[0208]
初始化训练过程的相关超参数：即设置目标网络代替更新频率参数为τ、目标网络
更新周期回合数t、经验回放池大小为m、每个时间步从经验回放池采样得到的转移过程总数为batch、奖励折扣率为λ；
[0209]
10.2)使用马尔科夫决策过程mdp模型描述智能体的单步状态转移过程：
[0210]
10.2.1)结合再入运动模型描述单步转移过程：
[0211]
根据步骤1中构建的飞行器再入运动模型，设智能体的环境状态为轨迹状态变量s＝[r,θ,φ,γ,ψ]，智能体的动作为攻角增量δα；
[0212]
在每次状态转移过程中，智能体首先接收到一个当前环境状态信息s
t
＝[r
t
,θ
t
,φ
t
,γ
t
,ψ
t
]，而后做出相应的动作a
t
＝δα
t
，模型在执行该动作后，根据步骤7得到的参考轨迹训练数据集和s
t
、a
t
，通过龙格库塔积分将当前状态s
t
转移为下一个状态s
t+1
，同时智能体将收到一个收益r
t+1
，其中，t表示智能体当前状态，t+1表示智能体执行动作后转移到的下一个状态；
[0213]
10.2.2)在智能体动作中添加随机噪声并抽取新动作替代原动作：
[0214]
为使智能体能够探索到更多的动作，具备更强的环境适应能力，在智能体所选择的动作a
t
＝δα
t
中添加随机噪声；
[0215]
设置参数v
ar
作为标准差，并以强化学习算法选择出来的原动作a
t
作为平均值，构造一个正态分布函数，然后从正态分布函数中随机抽取一个新的动作a
t
代替a
t
。
[0216]
10.3)将步骤10.2中的单步转移过程以(s
t
,a
t
,r
t+1
,s
t+1
)的形式存储进经验回放池，其中，s
t
表示当前时刻状态、a
t
表示当前时刻采取的攻角增量动作、r
t+1
表示执行动作之后得到的奖励值、s
t+1
表示下一时刻状态；
[0217]
10.4)从经验回放池中随机抽取一些样本输入进actor网络和critic网络中，通过性能指标的策略梯度最大化累计奖励来更新actor网络中的actor_eval子网络参数，通过最小化损失函数即当前累计奖励与目标累计奖励的均方误差来更新critic网络中的critic_eval子网络参数，每过t 回合根据τ和原网络参数进行加权更新actor网络中的actor_target 子网络参数和critic网络中的critic_target子网络参数，
[0218]
其中，n表示从经验回放池中采样的样本数量，si为当前状态量，ai为神经网络根据当前状态量选择出的攻角增量指令，yi＝r
i+1
+λq
′
(s
i+1
,μ
′
(s
i+1
|w
μ
′
)wq′
)，表示目标累计奖励，q(si,ai|wq)表示critic_eval网络的输出值，w表示网络的权重和阈值参数，w的上标表示此参数的归属网络，上标q表示critic_eval网络，上标q
′
表示critic_target网络，上标μ表示actor_eval网络，上标μ
′
表示actor_target网络，μ
′
(s
i+1
|w
μ
′
)表示actor_target 网络的输出值；q
′
(s
i+1
,μ
′
(s
i+1
|w
μ
′
)|wq′
)表示critic_target网络的输出值，actor_eval网络的输出值为μ(si|w
μ
)。
[0219]
10.5)重复步骤10.4)，当智能体的累计奖励收敛于一个最大值时，得到训练好的actor 网络和critic网络，将actor网络中的动作评估子网络actor_eval网络作为制导网络。
[0220]
步骤11，在线获取飞行器再入段制导指令。
[0221]
飞行器在飞行期间，机载计算机读取飞行器导航系统测得的实时飞行状态量，将实时飞行的状态量作为制导网络的输入进行前向传播，得到实时的制导指令，根据指令飞行器将产生攻角变化实时跟踪预定高度曲线。
[0222]
以下结合仿真实验对本实施例的技术效果进一步说明：
[0223]
一、仿真条件
[0224]
飞行器初始条件为：
[0225]
初始地心距初始经度初始纬度初始航迹角初始航向角初始速度
[0226]
飞行器终端条件为：
[0227]
终端地心距终端经度终端纬度终端航迹角终端航向角终端速度
[0228]
所有数值仿真均在cpu为i7-8700f(3.20ghz)的pc上进行，软件环境为matlab 和pycharm。
[0229]
二、仿真内容
[0230]
仿真实验1：
[0231]
在上述仿真条件下，将本发明中每个离散时间步的升力系数和阻力系数固定在正向拉偏20％，测试本发明的收敛性，结果如图5所示。从图5可见，本发明的累计奖励曲线在 200回合以内收敛到一个最大值，保证了在正向拉偏情况下的收敛性和数据的的准确性；
[0232]
仿真实验2：
[0233]
在上述仿真条件下，将本发明中每个离散时间步的升力系数和阻力系数固定在正向拉偏20％，测试本发明的制导网络在正向拉偏情况下对高度的跟踪效果，结果如图6所示，其中虚线为标称气动参数条件下求解计算的最优高度轨迹，实线为本发明的制导网络生成的高度曲线。根据图6可计算出本发明的制导网络在正向拉偏情况下的终端高度制导精度为0.52％。
[0234]
仿真实验3：
[0235]
由于rlv的飞行环境复杂多变，为模拟更加真实的飞行环境，在上述仿真条件下，将本发明中每个离散时间步的气动系数在[0,20％]区间内随机拉偏，且升力系数和阻力系数的拉偏量不同，测试本发明的制导网络在正向随机拉偏的复杂环境情况下对高度的跟踪效果，结果如图7所示，其中虚线为标称气动参数条件下求解计算的最优高度轨迹，实线为本发明的制导网络生成的高度曲线。根据图7可计算出本发明的制导网络在正向随机拉偏的复杂环境下的终端高度制导精度为0.1％。
[0236]
仿真实验4：
[0237]
在上述仿真条件下，将本发明中每个离散时间步的气动系数在[0,20％]区间内随机拉偏，且升力系数和阻力系数的拉偏量不同，测试本发明中飞行器在正向随机拉偏的复杂环境下追踪高度的难易程度，结果如图8所示。从图8可见，本发明的奖励曲线在-1000到-1150 之间波动，轨迹末端的奖励下降较快，表示本发明中在正向随机拉偏的复杂情况下飞
行器在轨迹末端的跟踪难度较大。
[0238]
仿真实验5：
[0239]
在上述仿真条件下，将本发明中每个离散时间步的升力系数和阻力系数固定在反向拉偏20％，测试本发明的收敛性，结果如图9所示。从图9可见，本发明的累计奖励曲线在 200回合以内收敛到一个最大值，保证了反向拉偏情况的收敛性和数据的准确性。
[0240]
仿真实验6：
[0241]
在上述仿真条件下，将本发明中每个离散时间步的升力系数和阻力系数固定在反向拉偏20％，测试本发明的制导网络在反向拉偏情况下对高度的跟踪效果，结果如图10所示，其中虚线为标称气动参数条件下求解计算的最优高度轨迹，实线为本发明的制导网络生成的高度曲线，根据图10可计算出本发明的制导网络在反向拉偏情况下的终端高度制导精度为2.9％。
[0242]
仿真实验7：
[0243]
在上述仿真条件下，将本发明中每个离散时间步的气动系数在[-20％,0]区间内随机拉偏，且升力系数和阻力系数的拉偏量不同，测试本发明的制导网络在反向随机拉偏的复杂环境情况下对高度的跟踪效果，结果如图11所示，其中虚线为标称气动参数条件下求解计算的最优高度轨迹，实线为本发明的制导网络生成的高度曲线，根据图11可计算出本发明的制导网络在反向随机拉偏的复杂情况下的终端高度制导精度为2.4％。
[0244]
仿真实验8：
[0245]
在上述仿真条件下，将本发明中每个离散时间步的气动系数在[-20％,0]区间内随机拉偏，且升力系数和阻力系数的拉偏量不同，测试本发明飞行器在反向随机拉偏的复杂环境下追踪高度的难易程度，结果如图12所示。从图12可见，本发明的奖励曲线在-1000到-1025 之间波动，轨迹末端的奖励略有上升，表示本发明中在反向随机拉偏的复杂情况下飞行器在轨迹末端的跟踪难度较小。
[0246]
从上述随机拉偏条件下的仿真实验可以看出，本发明所提出的方案对在飞行器处于复杂环境下时对高度的跟踪效果好、精度高，且训练得到的制导网络能够较好地适应气动参数偏差。
技术特征：
1.一种基于强化学习算法的飞行器再入跟踪制导方法，其特征在于，包括如下步骤：(1)将飞行器再入段轨迹优化描述为由数学模型、边界条件、容许控制、性能指标、过程约束构成的连续最优控制问题p0；(2)对p0进行更换形式、松弛变量、软化约束、逐次线性化方法的凸化处理，得到序列凸最优控制问题p1，采用伪谱法对该p1进行离散参数化处理，得到序列二阶锥规划问题p2；(3)采用内点法求解序列二阶锥规划问题p2，得最优参考轨迹；(4)对最优参考轨迹进行采样得到参考轨迹训练数据集；(5)构建神经网络actor网络：(5a)建立由第一输入层、第一隐藏层、第一输出层依次连接组成的动作评估子网络actor_eval，该第一输入层输入的状态变量为当前状态智能体的地心距、经度、纬度、航迹角、航向角这五个状态变量，该第一输出层输出信息为当前状态的攻角增量指令；(5b)建立由第二输入层、第二隐藏层、第二输出层依次连接组成的动作目标子网络actor_target，该第二输入层输入的状态变量为智能体下一个状态的地心距、经度、纬度、航迹角、航向角这五个状态变量，该第二输出层输出信息为下一个状态的目标攻角增量指令；(5c)将动作评估子网络actor_eval网络和动作目标子网络actor_target网络并联，构成actor网络，用于接收经验回放池的状态信息，并输出攻角增量指令信息；(6)构建神经网络critic网络：(6a)建立由第三输入层、第三隐藏层、第三输出层依次连接组成的价值评估子网络critic_eval，该第三输入层输入的变量为当前状态智能体的地心距、经度、纬度、航迹角、航向角这五个状态变量以及actor_eval输出的攻角增量指令，该第三输出层输出信息为智能体这一状态采取指令后产生的累计奖励；(6b)建立由第四输入层、第四隐藏层、第四输出层依次连接组成的价值目标子网络critic_target，该第四输入层输入的变量为下一个状态智能体的地心距、经度、纬度、航迹角、航向角这五个状态变量以及actor_target输出的攻角增量指令，该第四输出层输出信息为智能体下一状态采取指令后产生的目标累计奖励；(6c)将价值评估子网络critic_eval网络和价值目标子网络critic_target网络并联，构成critic网络，用于接受经验回放池的状态信息和actor网络输出的攻角增量指令信息，并输出智能体采取指令后产生的累计奖励信息；(7)设计actor网络和critic网络的奖励函数r均为：r＝r1+r2，其中r1为目标奖励，r2为约束奖励；(8)利用(4)中所得到的参考轨迹训练数据集使用深度确定性策略梯度算法ddpg，对actor网络和critic网络进行并行离线训练，当智能体的累计奖励收敛于一个最大值时，得到训练好的制导网络；(9)在线获取飞行器再入段制导指令并实时跟踪：飞行器在飞行期间，机载计算机读取飞行器导航系统测得的实时飞行状态量，将实时飞行的状态量作为制导网络的输入进行前向传播，得到实时的制导指令；飞行器根据指令将产生攻角变化，实时跟踪预定高度曲线。2.根据权利要求1所述的方法，其特征在于：步骤(1)中连续最优控制问题p0表示为：
p0:mins.t.其中，f1为状态微分方程，表示为：f1：r、v、ω、l、d分别为无量纲的地心距、速度、地球自转角速率和升、阻力加速度，j表示连续最优控制问题的性能指标，x表示飞行器的五个状态变量，和分别表示飞行器的初末无量纲状态，e表示飞行器所具有的无量纲能量，e0和e
f
分别表示飞行器的初末无量纲能量，θ、φ分别为经度和纬度，r0为地球半径，γ为航迹角，ψ为航向角，σ为倾侧角，σ
min
和σ
max
分别为倾侧角的下限和上限，q、n分别为飞行器热流率，动压和过载，q
max
、n
max
分别为飞行器能承受热流率、动压和过载的最大值，k
q
和ρ分别表示热流率计算系数和大气密度。3.据权利要求1所述的方法，其特征在于：步骤(2)中得到的序列凸最优控制问题p1，表示为：p1:find x(e),u(e),mins.t.x
′
＝a(x
(k)
)x+b(x
(k)
)u+b(x
(k)
)|x(e)-x
(k)
(e)|≤δ(e)|≤δω
l
≤u1≤ω
h
其中，u＝[u1,u2]＝[cosσ,sinσ]表示控制量，r、v分别为无量纲的地心距、速度，j表示连续最优控制问题的性能指标，x表示飞行器的五个状态变量，表示飞行器的无量纲初状态，e表示飞行器所具有的无量纲能量，e0和e
f
分别表示飞行器的初末无量纲能量，θ、φ分别为经度和纬度，γ为航迹角，ψ为航向角，σ为倾侧角，σ
min
和σ
max
分别为倾侧角的下限和上限，分别表示经度和纬度的松弛变量，上标的(k)表示变量处于第k次迭代中，k表示迭代计数，a(x
(k)
)、b(x
(k)
)分别表示状态量和控制量的系数矩阵，b(x
(k)
)表示余矢量矩阵，c
θ
、分别表示经度和纬度的惩罚系数，c1、c0分别表示性能指标逐次线性化过程中关于无量纲地心距r的一阶项和余项，ε
ψ
表示航向角正则化项系数，δ表示信赖域约束向量，表示凸化后的过程约束，ω
l
、ω
h
分别表示倾侧角上限和下限对应的余弦值，分别表示飞行器末状态的无量纲地心距、经度、纬度、航迹角和航向角。4.据权利要求1所述的方法，其特征在于：步骤(2)中用伪谱法对该p1进行离散参数化处理，得到序列二阶锥规划问题p2表示为：p2:find x
n
,u
n
,mins.t.s.t.s.t.s.t.ω
l
≤(u1)
m
≤ω
h
,,,,其中，u＝[u1,u2]＝[cosσ,sinσ]表示控制量，r、v分别为无量纲的地心距、速度，j表示连续最优控制问题的性能指标，x表示飞行器的五个状态变量，表示飞行器的无量纲初状态，e表示飞行器所具有的无量纲能量，e0和e
f
分别表示飞行器的初末无量纲能量，θ、φ分别为经度和纬度，γ为航迹角，ψ为航向角，σ为倾侧角，σ
min
和σ
max
分别为倾侧角的下限和上限，分别表示经度和纬度的松弛变量，上标的(k)表示变量处于第k次迭代中，k表示迭代计数，a(x
(k)
)、b(x
(k)
)分别表示状态量和控制量的系数矩阵，b(x
(k)
)表示余矢量矩阵，c
θ
、分别表示经度和纬度的惩罚系数，c1、c0分别表示性能指标逐次线性化过程中关于无量纲地
心距r的一阶项和余项，ε
ψ
表示航向角正则化项系数，δ表示信赖域约束向量，表示凸化后的过程约束，ω
l
、ω
h
分别表示倾侧角上限和下限对应的余弦值，分别表示飞行器末状态的无量纲地心距、经度、纬度、航迹角和航向角，n＝0,1,2...n+1表示每一个离散点，m＝1,2,3...n表示每一个lg点，n表示表示选取的lg点数，w表示每一个lg点对应的积分权重，d
i
表示基函数在lg点的导数。5.据权利要求1所述的方法，其特征在于：步骤(3)采用内点法求解序列二阶锥规划问题p2，实现如下：(3a)令迭代计数k＝0，给定离散点处的初始状态量(3a)令迭代计数k＝0，给定离散点处的初始状态量表示第n个离散点处第k次迭代的状态变量解，n＝0,1,2...n+1；(3b)将迭代计数k加一，采用内点法求解问题p2，得到解其中，z
(k+1)
表示第k+1次迭代的解向量，表示第k+1次迭代的状态变量解，表示第n个离散点处第k+1次迭代的状态变量解，n＝0,1,2...n+1，表示第k+1次迭代的控制量解，表示第m个离散点处第k+1次迭代的控制量解，m＝1,2,3...n，和表示第k+1次迭代的松弛变量解；(3c)设一个五维列向量ε表示收敛判断条件，判断每一个离散点处的状态变量解是否满足的收敛条件：若不满足，则令x
(k)
＝x
(k+1)
，返回(3b)；否则，输出最优解计算结束。6.据权利要求1所述的方法，其特征在于：actor网络和critic网络中的隐藏层均包含100个神经元，且隐藏层的激活函数均使用relu函数，输出层均使用tanh函数。7.据权利要求1所述的方法，其特征在于：步骤(7)中目标奖励与约束奖励，分别设置如下：根据任务目标，设置目标奖励r1：其中，设h
next
表示执行完动作δα后下一时间步的高度，δα表示当前状态的攻角增量指令，h0表示下一时间步对应标称高度；根据攻角增量大小，设置约束奖励r2＝-10|δα|。8.据权利要求1所述的方法，其特征在于：步骤(8)中使用ddpg算法对actor网络和critic网络进行并行离线训练，实现如下：
(8a)参数初始化：随机初始化actor_eval网络参数θ
μ
、actor_target网络参数θ
μ
′
；随机初始化critic_eval网络参数θ
q
、critic_target网络参数θ
q
′
；初始化训练过程的相关超参数：设置目标网络代替更新频率参数为τ、目标网络更新周期回合数t、经验回放池大小为m、每个时间步从经验回放池采样得到的转移过程总数为batch、奖励折扣率为λ；(8b)使用马尔科夫决策过程mdp模型描述智能体的单步状态转移过程，即智能体接收到当前环境状态信息s
t
，依据s
t
执行攻角增量α
t
动作之后转移至新环境状态s
t+1
，同时智能体接收到一个标量奖励r
t+1
；(8c)将(8b)中的单步转移过程以(s
t
,a
t
,r
t+1
,s
t+1
)的形式存储进经验回放池；其中，s
t
表示当前时刻状态、a
t
表示当前时刻采取的攻角增量动作、r
t+1
表示执行动作之后得到的奖励值、s
t+1
表示下一时刻状态；(8d)每次训练的过程中，都从经验回放池中随机抽取一些样本输入进actor网络和critic网络中，通过性能指标的策略梯度最大化累计奖励来更新actor网络中的actor_eval子网络参数，通过最小化当前累计奖励与目标累计奖励的均方误差来更新critic网络中的critic_eval子网络参数，每过t回合根据τ进行加权更新actor网络中的actor_target子网络参数和critic网络中的critic_target子网络参数；(8e)重复步骤(8d)，当智能体的累计奖励收敛于一个最大值时，得到训练好的制导网络。
技术总结
本发明公开了一种基于强化学习算法的飞行器再入跟踪制导方法，用于解决现有技术中对系统模型依赖性较大和适应性及制导性能差的问题。其实现方案为建立半速度坐标系下飞行器再入连续最优控制问题；将飞行器再入连续最优控制问题转化为序列凸最优控制问题；将序列凸最优控制问题转化为序列二阶锥规划问题；对序列二阶锥规划问题进行求解，获取飞行器最优倾侧角；对最优倾侧角轨迹采样得到参考轨迹训练数据集；构建神经网络和奖励函数；利用训练数据集对神经网络进行离线训练，直到累计奖励收敛到最大值，得到制导网络；通过制导网络在线获取飞行器再入制导指令，实时跟踪目标高度。本发明适应性强，制导精度高，可用于火箭回收。可用于火箭回收。可用于火箭回收。
