Synergistic Path Planning for Multiple Vehicles Based on an Improved Particle Swarm Optimization Method
-
摘要: 考虑气动、轨迹、约束、指标间的耦合关系, 以多高超声速飞行器同时到达为目标建立了协同规划模型; 设计了一种自动满足终端约束的全新滑翔飞行剖面, 减少了规划算法需要处理的约束数量; 推导了滑翔段高精度解析解, 实现了过程约束和性能指标的快速求解; 提出了一种改进粒子群优化(Particle swarm optimization, PSO)算法, 借助强化学习方法构建协同需求与惯性权重间的动态映射网络, 提高了在线规划效率. 最后通过数学仿真验证了方法的正确性和有效性.Abstract: This paper researches the synergistic flight for multiple hypersonic vehicles. The synergistic planning problem is formulated in view of the nonlinear coupling among aerodynamics, the performance index, and the path constraints. Then, the gliding profile, which naturally satisfies the terminal constraints and decreases the constraints, is proposed. Meanwhile, accurate solutions are deduced in the glide phase, so path constraints and the performance index can be quickly derived. An improved particle swarm optimization (PSO) method is developed by building the network between synergistic requirements and the optimal inertial weight in PSO based on a reinforcement learning method. Thus, the efficiency online computational efficiency can be largely improved. Numerical simulation results indicate the efficiency of the proposed method.
-
现代防空反导体系已具备全天域、多层次、多维度防御能力; 在此背景下, 传统弹道弹和巡航弹已无法满足复杂对抗态势下的作战需求. 因此人们开始研究具有多空域机动能力的临近空间高超声速滑翔飞行器. 美国于2003年便提出了HTV-2 (Hypersonic technology vehicle 2)[1-2]计划, 并于2017年完成了高超声速滑翔试验[3]. 俄罗斯于2015年和2017年分别完成了高超声速飞行器“Yu-71”和“先锋”的飞行试验[4]. 日本于2019年起加紧研制高速滑翔飞行器[5]; 同年, 印度也完成了高超技术验证机首飞[6].
基于出色的突防能力, 多高超声速飞行器协同作战模式[7]同样备受关注. 但同时规划多条滑翔段轨迹具有一定难度: 1)约束多, 如载荷、热等过程约束以及位置、高度等终端约束; 2)耦合强, 如气动−约束−轨迹−指标; 3)干扰多、速度快、航程远. 综上, 在线规划方法必须具备很强的场景适应性, 同时还要精确控制各飞行器间的相对运动关系.
许多学者针对高超声速飞行器滑翔段轨迹规划问题展开了研究. 孙长银等[8]探讨了高超声速飞行器领域的新方法和新问题, 提出了一种通用的大包络飞行动力学模型; Liao等[9]基于有限时间非线性观测器实现了滑翔段轨迹的在线规划. 在协同规划方面, 王芳等[10]提出了一种时间最优导弹编队方法, Vincent等[11]和Chen等[12]分别用遗传算法和粒子群优化(Particle swarm optimization, PSO)算法实现了无人机编队控制, Xu等[13]提出了一种用于故障条件下多巡航弹协同控制的时变容错算法.
现有方法有效解决了滑翔段轨迹设计问题, 但协同规划方面的研究多以巡航弹/无人机为研究对象. 针对上述问题, 从提高在线规划效率出发, 本文提出一种自动满足终端约束的滑翔段飞行剖面. 推导滑翔段解析解, 实现了飞行剖面的快速重构. 提出一种基于惯性权重智能选择的改进PSO算法, 用于完成在线协同规划. 最后通过数学仿真验证了方法的正确性和有效性.
1. 飞行器协同轨迹规划问题建模
1.1 滑翔段运动数学模型
高超声速飞行器滑翔段运动方程为
$$\left\{ \begin{aligned} & \dot V = - \dfrac{D}{m} - g\sin \gamma \\ & \dot \gamma = \dfrac{{L\cos \sigma }}{{mV}} - \left( {\dfrac{g}{V} - \frac{V}{r}} \right)\cos \gamma \\ & \dot \psi = \dfrac{{L\sin \sigma }}{{mV\cos \gamma }} + \dfrac{V}{r}\cos \gamma \sin \psi \tan \phi \\ &\dot r = V\sin \gamma \\ &\dot \theta = \dfrac{{V\cos \gamma \sin \psi }}{{r\cos \phi }} \\ & \dot \phi = \dfrac{{V\cos \gamma \cos \psi }}{r} \end{aligned} \right.$$ (1) 式中,
$V$ 为速度,$\gamma $ 为飞行路径角,$\psi $ 为航向角,$r$ 为飞行器质心到地心的距离,$\theta $ 为经度,$\phi $ 为纬度;$m$ 为质量,$\sigma $ 为倾侧角,$g$ 为地球引力加速度.$D = {c_D}q{S_c}$ 为气动阻力,$L = {c_L}q{S_c}$ 为气动升力, 其中${c_D}$ 为阻力系数,${c_L}$ 为升力系数,${S_c}$ 为特征面积,$q$ 为飞行动压.$q=\rho {V^2}/2,$ 其中$\rho $ 为大气密度, 采用指数函数计算$$\rho = {\rho _0}{{\rm{e}} ^{ - \beta h}}$$ (2) 式中,
$h = r - {R_e}$ 为飞行高度$({R_e}$ 为地球平均半径),$\beta = 1.4064 \times {10^{ - 4}}\;{{\rm{m}}^{ - 1}}$ 且${\rho _0}=1.225\;{\rm{kg/}}{{\rm{m}}^{\rm{3}}}$ .设阻力系数按如下规律变化:
$${c_D} = {c_{D0}}{V^{{\beta _1}}}{{\rm{e}} ^{ - {\beta _2}h}}$$ (3) 式中,
$\alpha $ 为飞行攻角,${c_{D0}},{\beta _1}$ 和${\beta _2}$ 为给定参数. 取${c_D}_0=0.156,{\beta _1}=2.6 \times {10^{ - 4}},{\beta _2} = 1.5 \times {10^{ - 6}}$ .定义从当前位置到终端位置间的航程为剩余航程, 记为
${R_L}$ , 则有$$\left\{ \begin{aligned} & {R_L} = {\mu _L}{R_e} \\ & {{\dot R}_L} = - V\cos \gamma \cos \zeta \\ & {\mu _L} = \arccos (\sin {\phi _T}\sin \phi + \cos {\phi _T}\cos \phi \cos \Delta \theta ) \end{aligned} \right.$$ (4) 式中,
$\zeta = {\psi _w} - \psi $ 为航向偏差,$\Delta \theta = {\theta _T} - \theta ,{\theta _T}$ 和${\phi _T}$ 为目标经纬度;${\psi _w}$ 为期望航向角$${\psi _w} = \arccos \left( {\dfrac{{\sin {\phi _T} - \sin \phi \cos {\mu _L}}}{{\cos \phi \sin {\mu _L}}}} \right)$$ (5) 设滑翔段约束条件为
$$\left\{ \begin{aligned} & {{\boldsymbol{c}}_{{\rm{end}}}} = {\left[ {{h_f} - {h_{{\rm{ter}}}},{\gamma _f} - {\gamma _{{\rm{ter}}}},{R_{Lf}},{\zeta _f}} \right]^{\rm{T}}} = {\bf{0}} \\ & {{\boldsymbol{c}}_{{\rm{path}}}} = {\left[ {\dfrac{{{N_m}}}{{{N_{\max }}}},\dfrac{q}{{{q_{\max }}}},\dfrac{{\left| \alpha \right|}}{{{\alpha _{\max }}}},\dfrac{{\left| \sigma \right|}}{{{\sigma _{\max }}}}} \right]^{\rm{T}}} < {\bf{1}} \end{aligned} \right.$$ (6) 式中,
${{\boldsymbol{c}}_{{\rm{end}}}}$ 为终端约束,${{\boldsymbol{c}}_{{\rm{path}}}}$ 为过程约束. 下标“$f$ ”表示实际终端状态, 下标“ter”表示期望终端状态. 下标“max”表示该过程变量允许的最大值,${N_m}$ 为法向过载的绝对值.1.2 多飞行器协同轨迹规划模型
设协同飞行器数量为
${N_{{\rm{hyp}}}}.$ 协同规划应考虑各个飞行器的总体性能以及各自飞行约束, 记为${\boldsymbol{c}}_{{\rm{end}}}^{\left( 1 \right)}, \cdots ,{\boldsymbol{c}}_{{\rm{end}}}^{\left( {{N_{{\rm{hyp}}}}} \right)},{\boldsymbol{c}}_{{\rm{path}}}^{\left( 1 \right)}, \cdots ,{\boldsymbol{c}}_{{\rm{path}}}^{\left( {{N_{{\rm{hyp}}}}} \right)}.$ 考虑到各飞行器的初始运动状态各异, 设其剩余飞行时间为$t_{{\rm{go}}}^{(i)},$ 协同规划指标为$${J_{{\rm{syn}}}} = \max \left( {t_{{\rm{go}}}^{\left( i \right)}} \right) - \min \left( {t_{{\rm{go}}}^{\left( i \right)}} \right),\;\;\;i = 1,2, \cdots ,{N_{{\rm{hyp}}}}$$ (7) 2. 滑翔段飞行剖面设计
设飞行高度为剩余航程的函数为
$$h = f( {{R_L}}) = - \frac{1}{{\bar \beta }}\ln F( {{R_L}} ) = - \frac{1}{{\bar \beta }}\sum\limits_{i = 0}^n {\left( {{a_i}R_L^i} \right)} $$ (8) 式中,
${a_i}$ 为未知系数,$\bar \beta = \beta + {\beta _2},F({R_L})$ 为高阶多项式. 结合式(1)和式(8), 可得$$ \frac{{{\rm{d}}h}}{{{\rm{d}}{R_L}}} = - \frac{{F'( {{R_L}} )}}{{\bar \beta F( {{R_L}} )}} = - \frac{{\tan \gamma }}{{\cos \zeta }} $$ (9) $$\frac{{{{\rm{d}}^2}h}}{{{\rm{d}}R_L^2}} = \frac{{{{\left[ {F'( {{R_L}} )} \right]}^2} - F''( {{R_L}} )F( {{R_L}} )}}{{\bar \beta {F^2}( {{R_L}} )}} = - \frac{{{\rm{d}}\left( {\frac{{\tan \gamma }}{{\cos \zeta }}} \right)}}{{{\rm{d}}{R_L}}}$$ (10) 式中,
$$ F'({R_L})=\displaystyle\sum\limits_{i = 1}^5 {i{a_i}R_L^{i - 1}},F''({R_L})=\displaystyle\sum\limits_{i = 2}^5 i\left( {i - 1} \right){a_i}R_L^{i - 2} $$ 定义
${K_h} = \dfrac{{\rm{d}}^2h}{{\rm{d}}R_L^2}$ , 本文认为${\dot \psi _w} = 0,\dot \zeta = - \dot \psi.$ 由式(11)可得升力为$$L = \dfrac{{{K_h}{{\cos }^2}\gamma {{\cos }^2}\zeta + \frac{g}{{{V^2}}} - \frac{1}{r} + \frac{{{K_b}}}{r}}}{\frac{\cos \sigma - \sin \zeta \sin \gamma \sin \sigma }{m{V^2}\cos \gamma }}$$ (11) 式中,
$ {K_b} = 2\sin \gamma \cos \gamma \sin \psi \tan \phi \sin \zeta$ .设
$n = 5,$ 则确定$F({R_L})$ 需建立6个方程. 由于初始和终端高度已知, 由式(8)可得$$\left\{ \begin{aligned} & {{\rm{e}} ^{ - \bar \beta {h_0}}} = \displaystyle\sum\limits_{i = 0}^5 {{a_i}R_{L0}^i} \\ & {{\rm{e}} ^{ - \bar \beta {h_f}}} = \displaystyle\sum\limits_{i = 0}^5 {{a_i}R_{Lf}^i} \end{aligned} \right.$$ (12) 式中, 下标“0”表示滑翔段初始运动状态.
由于初始和终端飞行路径角已知, 认为
$\cos \zeta \approx 1,$ 由式(9)可得$$\left\{ \begin{aligned} & \displaystyle\sum\limits_{i = 1}^5 {{a_i}\left( {R_{L0}^i - \dfrac{{iR_{L0}^{i - 1}}}{{\bar \beta \tan {\gamma _0}}}} \right) + {a_0} = 0} \\ & \displaystyle\sum\limits_{i = 1}^5 {{a_i}\left( {R_{Lf}^i - \dfrac{{iR_{Lf}^{i - 1}}}{{\bar \beta \tan {\gamma _f}}}} \right) + {a_0} = 0} \end{aligned} \right.$$ (13) 设
${h_1},{h_2}$ 分别为剩余航程2/3处(记为${R_L}_1)$ 和1/3处(记为${R_L}_2)$ 的高度值, 由式(8)可得:$$\left\{ \begin{aligned} & {{\rm{e}} ^{ - \bar \beta {h_1}}} = \displaystyle\sum\limits_{i = 0}^5 {\left( {{a_i}R_{L1}^i} \right)} \\ & {{\rm{e}} ^{ - \bar \beta {h_2}}} = \displaystyle\sum\limits_{i = 0}^5 {\left( {{a_i}R_{L2}^i} \right)} \end{aligned} \right.$$ (14) 定义剖面设计变量
${{\boldsymbol{u}}_h} = {[{h_{{\kern 1pt} 1}},{h_2}]^{\rm{T}}}$ , 则式(7)对应的优化变量分别为${\boldsymbol{u}}_{{\rm{syn}}}^{\left( 1 \right)} = [h_1^{\left( 1 \right)}, \cdots ,h_1^{({N_{{\rm{hyp}}}})}]$ 和${\boldsymbol{u}}_{{\rm{syn}}}^{\left( 2 \right)} = [h_2^{\left( 1 \right)}, \cdots ,h_2^{({N_{{\rm{hyp}}}})}]$ , 共$2{N_{{\rm{hyp}}}}$ 个.最后, 设计倾侧角变化规律为
$$\sigma =\frac{ \zeta {\sigma _{\max }}}{{\zeta _{\max }}}$$ (15) 式中,
${\zeta _{\max }}$ 为给定正数.3. 滑翔段高精度解析解推导
滑翔段高度的解析解即式(8). 考虑到滑翔段中航向偏差
$\zeta $ 和飞行路径角$\gamma $ 较小, 由式(1)和式(9)可得飞行路径角的解析解为$$\gamma = - f'\left( {{R_L}} \right) = \dfrac{{\displaystyle\sum\limits_{i = 1}^5 {\left( {{a_i}R_{L0}^{i - 1}} \right)} }}{{\bar \beta \displaystyle\sum\limits_{i = 0}^5 {\left( {{a_i}R_{L0}^i} \right)} }}$$ (16) 考虑到滑翔段中阻力通常远大于引力沿速度方向的分量, 根据式(1)可得
$$\frac{{{\rm{d}}V}}{{{\rm{d}}h}} = - \frac{{{c_{D0}}{S_c}{\rho _0}}}{{2m}}\frac{{{V^{n + 1}}{{\rm{e}} ^{ - \bar \beta h}}}}{\gamma }$$ (17) 定义
$ {\tau _D} = {c_{D0}}{S_c}{\rho _0}/(2m)$ , 因此$$ - {V^{ - \left( {{\beta _1} + 1} \right)}}{\rm{d}}V = {\tau _D}F\left( {{R_L}} \right){\rm{d}}{R_L}$$ (18) 对上式两端积分, 得
$${V^{ - {\beta _1}}} = {\tau _D}\int_0^{{R_L}} {F\left( {{R_L}} \right){\kern 1pt} } {\rm{d}}{R_L} + {C_V}$$ (19) 式中,
${C_V}$ 为积分常数. 定义$$y\left( {{R_L}} \right) = \int_0^{{R_L}} {F\left( {{R_L}} \right){\kern 1pt} } {\rm{d}}{R_L}=\sum\limits_{i = 0}^n {\frac{{{a_i}R_L^{i + 1}}}{{i + 1}}} $$ (20) 故
${C_V} = {V_0}^{ - {\beta _1}} - {\tau _D}y({R_{L0}}),$ 速度解析解为$$V = {\left[ {{\tau _D}y\left( {{R_L}} \right) + {C_V}} \right]^{ - \frac{1}{{{\beta _1}}}}}$$ (21) 根据式(21), 动压的解析表达式为
$$q = 0.5{\rho _0}{{\rm{e}} ^{ - \beta h}}{\left[ {{\tau _D}y\left( {{R_L}} \right) + {C_V}} \right]^{ - \frac{2}{{{\beta _1}}}}}$$ (22) 令
$\cos \zeta = \cos \gamma = 1,$ 则升力的解析解为$$L = m\frac{{{K_h}{V^2} + g}}{{\cos \sigma }}$$ (23) 由式(23), 可以解析估计出法向过载
$${N_m} \le \frac{{g + f''\left( {{R_L}} \right){{\left[ {{\tau _D}y\left( {{R_L}} \right) + {C_V}} \right]}^{ - \frac{2}{{{\beta _1}}}}}}}{{g\cos {\sigma _{\max }}}}$$ (24) 最后, 由式(4)可得
$$ - \frac{{{\rm{d}}t}}{{{\rm{d}}{R_L}}} = \frac{1}{{V\cos \gamma \cos \zeta }} = {\left[ {{\tau _D}y\left( {{R_L}} \right) + {C_V}} \right]^{\frac{1}{{{\beta _1}}}}}$$ (25) 定义积分常数
${C_t},$ 考虑到$t=0$ 时${R_L}={R_L}_0,$ 且$t={t_{{\rm{go}}}}$ 时${R_{Lf}}=0;$ 结合式(20), 可得$${t_{{\rm{go}}}} = {\tau _D}\left[ {\sum\limits_{i = 0}^n {\frac{{{a_i}R_{L0}^{i + 2}}}{{\left( {i + 1} \right)\left( {i + 2} \right)}} + } \frac{{{C_V}}}{{{\tau _D}}}{R_{L0}}} \right]$$ (26) 4. 滑翔段轨迹协同规划算法
4.1 基本PSO算法
PSO算法用
${M_{{\rm{pso}}}}$ 个粒子寻找${D_{{\rm{pso}}}}$ 个参数. 设粒子$i$ 在空间中的位置和速度分别为${{\boldsymbol{x}}_i} = {[{x_{i,1}}, \cdots ,{x_{i,{D_{{\rm{pso}}}}}}]^{\rm{T}}},$ ${{\boldsymbol{v}}_i} = {[{v_{i,1}}, \cdots ,{v_{i,{D_{{\rm{pso}}}}}}]^{\rm{T}}},$ 其历史最佳位置为${{\boldsymbol{p}}_i} = [{p_{i,1}},$ $\cdots ,{p_{i,{D_{{\rm{pso}}}}}}]^{\rm{T}};$ 整个群体的历史最佳位置为${{\boldsymbol{p}}_g} = [{p_{g,1}},$ $\cdots ,{p_{g,{D_{{\rm{pso}}}}}}]^{\rm{T}},$ 相应性能指标为${g_{{\rm{best}}}}.$ PSO算法为[14-15]$$\begin{aligned}[b]v_{i,j}^{\left( {\ell + 1} \right)}=\;&wv_{i,j}^{\left( \ell \right)}+{c_1}{r_1}\left( {{p_{i,j}} - x_{i,j}^{\left( \ell \right)}} \right)+\\ &{c_2}{r_2}\left( {{p_{g,j}} - x_{i,j}^{\left( \ell \right)}} \right)\end{aligned}$$ (27) $$x_{i,j}^{\left( {\ell + 1} \right)} = x_{i,j}^{\left( \ell \right)} + v_{i,j}^{\left( {\ell + 1} \right)}$$ (28) 式中,
$j = 1,2, \cdots ,{D_{{\rm{pso}}}},\ell $ 为进化代数,$w \in [0,1]$ 为惯性权重,${c_1}={c_2}=2,{r_1},{r_2} \in [0,1]$ 由每代随机产生. 令${\phi _1} = {c_1}{r_1} + {c_2}{r_2},{\phi _2} = {c_1}{r_1}{p_{i,j}} + {c_2}{r_2}{p_{g,j}},$ 则粒子进化公式可改写为$$x_{i,j}^{\left( {\ell + 2} \right)} + \left( {{\phi _1} - w - 1} \right)x_{i,j}^{\left( {\ell + 1} \right)} + wx_{i,j}^{\left( \ell \right)} = {\phi _2}$$ (29) 式(29)避免了对粒子速度进行初始化及更新, 算法结构更简洁, 利于提高计算效率.
调整权重是控制粒子寻优过程的常用方法. 借助强化学习方法, 通过快速智能选择惯性权重实现用较少的粒子数和迭代次数找出满足协同规划参数(
${\boldsymbol{u}}_{{\rm{syn}}}^{\left( 1 \right)}$ 和${\boldsymbol{u}}_{{\rm{syn}}}^{\left( 2 \right)}$ ).4.2 最优权重在线生成方法
设惯性权重为进化代数的分段函数如下:
$$w\left( \ell \right) = \left\{ {\begin{array}{*{20}{l}} &{{w_0} + \dfrac{{{w_1} - {w_0}}}{{\frac{\ell _{\max }}{\rm{3}}}}\ell },&{0 \le \ell < \dfrac{{{\ell _{\max }}}}{3}}\\ &{{w_1} + \dfrac{{{w_2} - {w_1}}}{{\frac{\ell _{\max }}{{3}}}}\ell },&{\dfrac{{{\ell _{\max }}}}{3} \le \ell < \dfrac{{2{\ell _{\max }}}}{3}}\\ &{{w_2} + \dfrac{{{w_3} - {w_2}}}{{\frac{\ell _{\max }}{{3}}}}\ell },&{\dfrac{{2{\ell _{\max }}}}{3} \le \ell < {\ell _{\max }}} \end{array}} \right.$$ (30) 式中,
${w_0} \sim {w_3}$ 为未知量,${\ell _{\max }}$ 为最大进化代数.定义自变量
${{\boldsymbol{u}}_w} = {[{w_0},{w_{{\kern 1pt} 1}},{w_2},{w_3}]^{\rm{T}}}$ 以及动作${\boldsymbol{d}} = $ ${[{d_1},{d_2},{d_3},{d_4}]^{\rm{T}}},$ 其中${d_1} \sim {d_4} \in \left[ {0,1} \right]$ 为连续值, 分别对应${w_0} \sim {w_3}$ 的更新:$$w_j^{} = \left( {2{d_j} - 1} \right){n_j}w_j^{}$$ (31) 式中,
$j = 1,2,3,4;{n_1} \sim {n_4}$ 为预置系数.设状态值
${\boldsymbol{s}} = {{\rm{[}}{p_{g,1}}{\rm{,}} \cdots {\rm{,}}{\kern 1pt} {p_{g,{D_{{\rm{pso}}}}}}{\rm{]}}^{\rm{T}}}$ 以及观测值${{\boldsymbol{o}}_t} = $ $ {{\rm{[}}t_{{\rm{go}}}^{{\rm{(a)}}},t_{{\rm{go}}}^{{\rm{(s)}}},t_{{\rm{go}}}^{{\rm{(l)}}}{\rm{]}}^{\rm{T}}},$ 其中$t_{{\rm{go}}}^{{\rm{(a)}}},t_{{\rm{go}}}^{{\rm{(s)}}},t_{{\rm{go}}}^{{\rm{(l)}}}$ 分别为$t_{{\rm{go}}}^{(i)}$ 的平均值、标准差以及最大/小值之差. 定义函数网络${Q_{{\rm{net}}}}$ 来表征状态${\boldsymbol{s}}$ 和动作${\boldsymbol{d}}$ 的优劣, 其中优劣性由奖励值$r = {g_{{\rm{best}}}}$ 表示. 定义策略网络${P_{{\rm{net}}}}$ 用以表征动作${\boldsymbol{d}}$ 的更新; 由式(35)可知,${P_{{\rm{net}}}}$ 即为${{\boldsymbol{u}}_w}$ 的求解器. 根据上述定义, 记${Q_{{\rm{net}}}}$ 和${P_{{\rm{net}}}}$ 的值分别为${{\boldsymbol{\theta}} _Q}$ 和${{\boldsymbol{\theta }}_P},$ 提出最优惯性权重训练方法如下:步骤 1. 随机初始化
${{\boldsymbol{\theta}} _Q}$ 和${{\boldsymbol{\theta }}_P},$ 并令${{\boldsymbol{\theta}} '_Q} = {{\boldsymbol{\theta}} _Q},{{\boldsymbol{\theta}} '_P} = {{\boldsymbol{\theta }}_P};$ 设$k = 1,$ 随机初始化${{\boldsymbol{u}}_w}$ ;步骤 2. 根据
${{\boldsymbol{u}}_w}$ 执行PSO, 得到状态${{\boldsymbol{s}}_k}$ ;步骤 3. 随机初始化扰动量
$\Delta {{\boldsymbol{d}}_k},$ 得到动作${{\boldsymbol{d}}_k}= {P_{{\rm{net}}}}({{\boldsymbol{\theta}} _P},$ ${{\boldsymbol{s}}_k}) + \Delta {{\boldsymbol{d}}_k};$ 执行动作${{\boldsymbol{d}}_k}$ , 并更新${{\boldsymbol{u}}_w},$ 运行PSO得到奖励${r_k}$ 和状态${{\boldsymbol{s}}_k};$ 步骤 4. 将
$({{\boldsymbol{s}}_k},{{\boldsymbol{d}}_k},{r_k},{{\boldsymbol{s}}_{k + 1}})$ 放入经验回放池, 令$ k =$ $ k + 1,$ 返回步骤2, 直至$k = 1\;000$ ;步骤 5. 从经验池中
${N_R}$ 次随机取样; 设定系数${k_Q},$ 计算第$i{\kern 1pt} {\kern 1pt} \;(i = 1,2, \cdots ,{N_R})$ 次取样结果$${y_i} = {r_i} + {k_Q} \times {{\boldsymbol{\theta}} '_Q}\left( {{{\boldsymbol{s}}_i},{{\boldsymbol{d}}_i},{{\boldsymbol{\theta}} _P}} \right)$$ (32) 步骤 6. 以损失函数最小为目标更新
${{\boldsymbol{\theta }}_Q}$ $$\min \frac{1}{{{N_R}}}\sum\limits_{i = 1}^{{N_R}} {{{\left( {{y_i} - {{\boldsymbol{\theta}} _Q}\left( {{{\boldsymbol{s}}_i},{{\boldsymbol{d}}_i},{{\boldsymbol{\theta}} _P}} \right)} \right)}^2}} \Rightarrow {{\boldsymbol{\theta}} _Q}$$ (33) 步骤 7. 通过梯度求解更新
${{\boldsymbol{\theta}} _P} = {{\boldsymbol{\theta}} _P} + \Delta {{\boldsymbol{\theta}} _P}$ $$\Delta {{\boldsymbol{\theta}} _P} = \sum\limits_{i = 1}^{{N_R}} {\left( {\frac{{\partial {\theta _Q}\left( {{{\boldsymbol{d}}_i},{{\boldsymbol{\theta}} _P},{{\boldsymbol{s}}_i}} \right)}}{{\partial {{\boldsymbol{d}}_i}}} \times \frac{{\partial {{\boldsymbol{\theta}} _P}}}{{\partial {{\boldsymbol{s}}_i}}}} \right)} $$ (34) 步骤 8. 给定系数
${\tau _\theta },$ 令${{\boldsymbol{\theta}} '_Q} = {\tau _\theta }{{\boldsymbol{\theta}} _{Q}} + (1 - {\tau _\theta }){{\boldsymbol{\theta}} '_Q}, {{\boldsymbol{\theta}} '_P} =$ $ {\tau _\theta }{{\boldsymbol{\theta}} '_P} + (1 - {\tau _\theta }){{\boldsymbol{\theta}} '_P}$ .4.3 在线协同轨迹规划方法
设从当前位置到起点的航程为
${R_{{\rm{cov}}}},$ 更新周期为$\Delta {R_{{\rm{cov}}}}.$ 如图1所示, 在线协同规划具体步骤如下:步骤 1. 协同飞行剖面重构
计算观测值
${{\boldsymbol{o}}_t}$ 并与上一步获得的${\boldsymbol{u}}_{{\rm{syn}}}^{\left( 1 \right)}$ 和${\boldsymbol{u}}_{{\rm{syn}}}^{\left( 2 \right)}$ (如果是首次执行则随机初始化${\boldsymbol{u}}_{{\rm{syn}}}^{\left( 1 \right)}$ 和${\boldsymbol{u}}_{{\rm{syn}}}^{\left( 2 \right)}$ )共同输入训练好的策略函数网络, 输出惯性权重; 然后用PSO计算${\boldsymbol{u}}_{{\rm{syn}}}^{\left( 1 \right)}$ 和${\boldsymbol{u}}_{{\rm{syn}}}^{\left( 2 \right)}$ 的值, 使指标函数式(7)达到最小.步骤 2. 个体飞行剖面重构
a) 计算
${R_{L1}},{R_{L2}}$ 及其对应高度值${h_1},{h_2};$ b) 求解系数
${a_0} \sim {a_5}$ , 更新飞行剖面.步骤 3. 计算攻角和倾侧角, 更新运动状态.
步骤 4. 计算航程
${R_{{\rm{cov}}}};$ 根据${R_{{\rm{cov}}}},$ 每5个周期$(5\Delta {R_{{\rm{cov}}}})$ 执行一次协同剖面重构, 每1个周期$(\Delta {R_{{\rm{cov}}}})$ 执行一次个体剖面重构; 未到达更新周期时, 各飞行器只执行步骤3.循环执行步骤
$1\sim 4 $ , 直至飞行器到达目标点上空. 步骤2用于控制终端高度和飞行路径角, 但未考虑过程约束; 步骤1既要保证时间协同, 又要通过优化求解来满足过程约束, 同时还负责修正解析解与实际值间的偏差.5. 仿真分析
参考文献[16-17], 设滑翔段过程约束为
${q_{\max }}=$ ${\rm{ 60\;kPa}},{N_{\max }}= 3{\rm{.0}},{\alpha _{\max }}= 25^\circ ,{\sigma _{\max }} = 50^\circ .$ 考虑对4架总体参数不同的高超声速飞行器进行协同规划; 参考X-33[18], 设其特征面积分别为2.0 m2, 1.6 m2, 1.3 m2和1.0 m2, 质量分别为950 kg, 900 kg, 850 kg和800 kg.设PSO算法的最大进化代数为15, 粒子数为20, 约束条件用罚函数法[19]处理, 当
${J_{{\rm{syn}}}}$ < 1或$\ell ={\ell _{\max }}$ 时停止进化.仿真条件见表1, 初始航向角由式(5)求出, 故初始航向偏差为零. 取
$\Delta {R_{{\rm{cov}}}}$ = 30 km, 仿真结果见图2 ~ 5. 仿真结果表明, 各轨迹均能满足终端约束条件和过程约束条件, 且攻角和倾侧角曲线十分平滑. 各飞行器到达目标的时刻分别为735 s, 734 s, 734 s和734 s.表 1 初始状态和终端约束Table 1 The initial states and the terminal constraints初始值 经纬度 (°) 高度 (km) 速度 (m/s) 飞行路径角 (°) 剩余航程 (km) 飞行器 1 E 90, N 45 56 5400 −2.0 2304 飞行器 2 E 65, N 30 54 5300 −2.0 2263 飞行器 3 E 40, N 35 52 5200 −1.0 2324 飞行器 4 E 70, N 70 50 5100 −1.0 2285 终端值 E 60, N 50 25 — −1.0 0 若不进行协同规划, 各飞行剖面设计参数均取
${h_1}$ = 42 km,${h_2}$ = 38 km, 则各飞行器的到达时刻为694 s, 703 s, 733 s和744 s, 说明协同规划方法正确有效.考虑初始状态下第1次计算协同参数
$({\boldsymbol{u}}_{{\rm{syn}}}^{\left( 1 \right)}$ 和${\boldsymbol{u}}_{{\rm{syn}}}^{\left( 2 \right)})$ , 分别用基本PSO和改进PSO求解; 基本PSO中的惯性权重取0.7, 改进PSO中惯性权重由策略网络生成; 其他参数相同. 由于PSO算法受随机数影响, 分别进行30次独立仿真. 用进化代数表征计算速度并进行比较; 表2所示的结果表明, 改进PSO在计算效率上显著高于基本PSO.表 2 基本PSO和改进PSO计算效率对比Table 2 Comparison of the computation efficiency进化代数 最大值 最小值 平均值 标准差 基本 PSO 20 8 14 3.17 改进 PSO 16 6 10 2.69 最后, 为验证协同规划方法的抗干扰性, 设置服从高斯分布的干扰因素(见表3), 并进行500次蒙特卡洛打靶仿真. 500次仿真中, 最晚与最早到达时刻之差的最大值为5.37 s, 平均值为1.57 s, 说明算法能够在干扰因素的影响下完成在线协同规划.
表 3 滑翔段干扰因素设置Table 3 Disturbances in the glide phase序号 干扰因素 $3\sigma $值 1 初始速度 (m/s) 30 2 初始飞行路径角 (°) 0.2 3 初始高度 (m) 500 4 初始航向偏差 (°) 0.2 5 初始剩余航程 (km) 50 6 气动系数 (%) 10 7 大气密度 (%) 10 以飞行器1为例, 给出打靶仿真结果如图6和图7所示. 结果表明, 各种约束条件依然能够得到满足. 在高度和飞行路径角的解析解推导过程中未用到近似假设条件, 因此可忽略它们和真实值间的误差. 对于速度和剩余时间的解析解, 相应误差会随着时间而减小, 并且剩余时间是根据实时运动状态估计的, 因此不存在累计误差.
6. 结束语
本文建立了多高超声速滑翔飞行器协同规划问题数学模型. 提出了一种新的滑翔飞行剖面, 实现了终端约束的自动满足, 并推导了滑翔段解析运动方程. 提出基于强化学习方法的改进PSO算法, 提高了在线协同规划效率. 仿真结果表明, 在初始状态和总体参数各异的情况下, 多弹时间协同偏差不超过2 s, 并且能够克服干扰因素的影响.
-
表 1 初始状态和终端约束
Table 1 The initial states and the terminal constraints
初始值 经纬度 (°) 高度 (km) 速度 (m/s) 飞行路径角 (°) 剩余航程 (km) 飞行器 1 E 90, N 45 56 5400 −2.0 2304 飞行器 2 E 65, N 30 54 5300 −2.0 2263 飞行器 3 E 40, N 35 52 5200 −1.0 2324 飞行器 4 E 70, N 70 50 5100 −1.0 2285 终端值 E 60, N 50 25 — −1.0 0 
下载: 导出CSV
表 2 基本PSO和改进PSO计算效率对比
Table 2 Comparison of the computation efficiency
进化代数 最大值 最小值 平均值 标准差 基本 PSO 20 8 14 3.17 改进 PSO 16 6 10 2.69
下载: 导出CSV
表 3 滑翔段干扰因素设置
Table 3 Disturbances in the glide phase
序号 干扰因素 $3\sigma $值 1 初始速度 (m/s) 30 2 初始飞行路径角 (°) 0.2 3 初始高度 (m) 500 4 初始航向偏差 (°) 0.2 5 初始剩余航程 (km) 50 6 气动系数 (%) 10 7 大气密度 (%) 10
下载: 导出CSV
-
[1] 王颖, 闵昌万, 刘秀明, 刘全军, 王官宇. 类HTV-2飞行器横侧向稳定设计研究. 宇航学报, 2017, 38(6): 583-589.Wang Ying, Min Chang-Wan, Liu Xiu-Ming, Liu Quan-Jun, Wang Guan-Yu. Study on lateral-directional stable design of HTV-2 like vehicle. Journal of Astronautics, 2017, 38(6): 583-589. [2] Li J, He M, Li X P, Zhang C F. Multi-physics modeling of electromagnetic wave-hypersonic vehicle interactions under high-power microwave illumination: 2D case. IEEE Transactions on Antennas and Propagation, 2018, 66(7): 3653-3664. doi: 10.1109/TAP.2018.2835300 [3] Borg M P, Kimmel R. Ground test of transition for Hifire-5b at flight-relevant attitudes. Journal of Spacecraft and Rockets, 2018, 55(6): 1-12. [4] 宋巍, 梁轶, 王艳, 袁成, 王竹溪. 2018年国外高超声速技术发展综述. 飞航导弹, 2019, 5: 7-12.Song Wei, Liang Tie, Wang Yan, Yuan Cheng, Wang Zhu- Xi. Summary of the development of hypersonic technology in 2018. Aerodynamic Missile Journal, 2019, 5: 7-12. [5] 杨金龙, 林旭斌. 日本高速/高超声速导弹计划分析. 飞航导弹, 2019, 409(1): 27-30.Yang Jin-Long, Lin Xu-Bin. Analysis on Japan high-speed/ hypersonic missile plane. Aerodynamic Missile Journal, 2019, 409(1): 27-30. [6] Qu Feng, Sun Di, Zuo Guang. A study of upwind schemes on the laminar hypersonic heating predictions for the reusable space vehicle. Acta Astronautica, 2018, 147: 412-420. doi: 10.1016/j.actaastro.2018.03.046 [7] Jia X, Wu S T, Wen Y M, Yao Z. A distributed decision method for missiles autonomous formation based on potential game. Journal of Systems Engineering and Electronics, 2019, 30(4): 738-748. doi: 10.21629/JSEE.2019.04.11 [8] 孙长银, 穆朝絮, 余瑶. 近空间高超声速飞行器控制的几个科学问题研究. 自动化学报, 2013, 39(11): 1901-1913. doi: 10.3724/SP.J.1004.2013.01901Sun Chang-Yin, Mu Zhao-Xu, Yu Yao. Some control problems for near space hypersonic vehicles. Acta Automatica Sinica, 2013, 39(11): 1901-1913. doi: 10.3724/SP.J.1004.2013.01901 [9] Liao Y X, Li H F, Bao W M. Three-dimensional diving guidance for hypersonic gliding vehicle via integrated design of FTNDO and AMSTSMC. IEEE Transactions on Industrial Electronics, 2018, 65(3): 2704-2715. doi: 10.1109/TIE.2017.2736499 [10] 王芳, 林涛, 张克, 崔乃刚. 多阶段高斯伪谱法在编队最优控制中的应用. 宇航学报, 2015, 36(11): 1262-1269. doi: 10.3873/j.issn.1000-1328.2015.11.007Wang Fang, Lin Tao, Zhang Ke, Cui Nai-Gang. Application of multi-phase Gauss pseudospectral method in optimal control for formation. Journal of Astronautics, 2015, 36(11): 1262-1269. doi: 10.3873/j.issn.1000-1328.2015.11.007 [11] Vincent R, Mohammed T, Gilles L. Fast genetic algorithm path planner for fixed-wing military UAV using GPU. IEEE Transactions on Aerospace and Electronic Systems, 2018, 54(5): 2105-2117. doi: 10.1109/TAES.2018.2807558 [12] Chen Q Y, Lu Y F, Jia G W, Li Y, Zhu B J, Lin J C. Path planning for UAVs formation reconfiguration based on Dubins trajectory. Journal of Central South University, 2018, 25(11): 2664-2676. doi: 10.1007/s11771-018-3944-z [13] Xu X G, Wei Z Y, Ren Z, Li S S. Time-varying fault-tolerant formation tracking based cooperative control and guidance for multiple cruise missile systems under actuator failures and directed topologies. Journal of Systems Engineering and Electronics, 2019, 30(3): 587-600. doi: 10.21629/JSEE.2019.03.16 [14] (王东风, 孟丽. 粒子群优化算法的性能分析和参数选择. 自动化学报, 2016, 42(10): 1552-1561.Wang Dong-Feng, Meng Li. Performance analysis and parameter selection of PSO algorithms. Acta Automatica Sinica, 2016, 42(10): 1552-1561. [15] 吕柏权, 张静静, 李占培, 刘廷章. 基于变换函数与填充函数的模糊粒子群优化算法. 自动化学报, 2018, 44(1): 74-86.Lv Bai-Quan, Zhang Jing-Jing, Li Zhan-Pei, Liu Yan-Zhang. Fuzzy particle swarm optimization based on filled function and transformation function. Acta Automatica Sinica, 2018, 44(1): 74-86. [16] Zeng L, Zhang H B, Zheng W. A three-dimensional predictor corrector entry guidance based on reduced-order motion equations. Aerospace Science and Technology, 2018, 73: 223−231 [17] Li M M, Hu J. An approach and landing guidance design for reusable launch vehicle based on adaptive predictor corrector technique. Aerospace Science and Technology, 2018, 75: 13-23. doi: 10.1016/j.ast.2017.12.037 [18] Bollino K, Ross M, Doman D. Optimal nonlinear feedback guidance for reentry vehicles. In: Proceedings of the 2016 AIAA Guidance, Navigation and Control Conference and Exhibit. Keystone, USA: AIAA Press, 2006. [19] Garg H. A hybrid PSO-GA algorithm for constrained optimization problems. Applied Mathematics and Computation, 2016, 274: 292-305. doi: 10.1016/j.amc.2015.11.001 期刊类型引用(12)
1. 孙瑞胜，陈洁卿，陆宇，刘宣廷. 突防弹道技术研究进展综述. 南京理工大学学报. 2025(01): 1-14+140 . 
百度学术2. 郭博，铁鸣，范文慧，李传旭. 基于滑模控制的高升阻比飞行器协同制导方法. 系统工程与电子技术. 2025(02): 580-590 . 百度学术3. 李磊，苏中，吴学佳，雷明，王一静. 粒子群优化噪声参数的行人导航零速修正算法. 传感技术学报. 2024(01): 42-49 . 百度学术4. 冯子鑫，薛文超，张冉，齐洪胜. 可回收火箭大气层内动力下降的多阶段鲁棒优化制导方法. 自动化学报. 2024(03): 505-517 . 本站查看5. 王凯威，尉静娴. 基于智能优化方法的工业机器人时间最优轨迹规划方法. 价值工程. 2024(18): 127-129 . 百度学术6. 张远龙，廖宇新，付钰雯，殷泽阳. 高超声速飞行器编队协同飞行轨迹规划方法. 飞控与探测. 2024(04): 26-32 . 百度学术7. 徐慧，蔡光斌，崔亚龙，侯明哲，姚二亮. 高超声速滑翔飞行器再入轨迹优化. 哈尔滨工业大学学报. 2023(04): 44-55 . 百度学术8. 薛光伟，辛万青，傅瑜. 能量优化分配再入轨迹快速规划方法. 弹道学报. 2023(02): 1-8 . 百度学术9. 崔正达，魏明英，李运迁. 基于速度预测的防空导弹中制导末段协同弹道规划方法. 系统工程与电子技术. 2023(09): 2912-2921 . 百度学术10. 薛光伟，辛万青，傅瑜. 升力式飞行器集结轨迹实时规划方法. 宇航学报. 2023(08): 1195-1202 . 百度学术11. Hui XU，Guangbin CAI，Chaoxu MU，Xin LI. Analytical reentry guidance framework based on swarm intelligence optimization and altitude-energy profile. Chinese Journal of Aeronautics. 2023(12): 336-348 . 必应学术12. 彭科科，刘国群. 基于解析解的小范围高频率机动飞行器制导律设计. 飞控与探测. 2022(06): 38-44 . 百度学术其他类型引用(12)
-
下载:
计量
- 文章访问数: 1833
- HTML全文浏览量: 547
- PDF下载量: 507
- 被引次数: 24
