2.793

2018影响因子

(CJCR)

  • 中文核心
  • EI
  • 中国科技核心
  • Scopus
  • CSCD
  • 英国科学文摘

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于改进粒子群算法的飞行器协同轨迹规划

周宏宇 王小刚 单永志 赵亚丽 崔乃刚

周宏宇, 王小刚, 单永志, 赵亚丽, 崔乃刚. 基于改进粒子群算法的飞行器协同轨迹规划. 自动化学报, 2020, 46(x): 1−7. doi: 10.16383/j.aas.c190865
引用本文: 周宏宇, 王小刚, 单永志, 赵亚丽, 崔乃刚. 基于改进粒子群算法的飞行器协同轨迹规划. 自动化学报, 2020, 46(x): 1−7. doi: 10.16383/j.aas.c190865
Zhou Hong-Yu, Wang Xiao-Gang, Shan Yong-Zhi, Zhao Ya-Li, Cui Nai-Gang. Synergistic path planning for multiple vehicles based on an improved particle swarm optimization method. Acta Automatica Sinica, 2020, 46(x): 1−7. doi: 10.16383/j.aas.c190865
Citation: Zhou Hong-Yu, Wang Xiao-Gang, Shan Yong-Zhi, Zhao Ya-Li, Cui Nai-Gang. Synergistic path planning for multiple vehicles based on an improved particle swarm optimization method. Acta Automatica Sinica, 2020, 46(x): 1−7. doi: 10.16383/j.aas.c190865

基于改进粒子群算法的飞行器协同轨迹规划


DOI: 10.16383/j.aas.c190865
详细信息
    作者简介:

    哈尔滨工业大学讲师. 主要研究方向为飞行器轨迹优化.E-mail: sd4574462@foxmail.com

    哈尔滨工业大学副教授. 主要研究方向为飞行器制导、导航与控制. 本文通信作者.E-mail: wangxiaogang@hit.edu.com

    哈尔滨建成集团研究员. 主要研究方向为飞行器总体设计.E-mail: yongzhi0451@sina.com

    北京航天晨信科技有限公司研究员. 主要研究方向为优化设计与系统工程.E-mail: zhaoyali_aerocim@163.com

    哈尔滨工业大学教授. 主要研究方向为飞行器制导、导航与控制.E-mail: cui_naigang@163.com

  • 基金项目:  中国博士后科学基金(2019M661290)资助

Synergistic Path Planning for Multiple Vehicles Based on an Improved Particle Swarm Optimization Method

More Information
  • Fund Project:  Supported by National Postdoctoral Science Foundation of P. R. China (2019M661290)
计量
  • 文章访问数:  23
  • HTML全文浏览量:  8
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-12-20
  • 修回日期:  2020-03-16

基于改进粒子群算法的飞行器协同轨迹规划

doi: 10.16383/j.aas.c190865
    基金项目:  中国博士后科学基金(2019M661290)资助
    作者简介:

    哈尔滨工业大学讲师. 主要研究方向为飞行器轨迹优化.E-mail: sd4574462@foxmail.com

    哈尔滨工业大学副教授. 主要研究方向为飞行器制导、导航与控制. 本文通信作者.E-mail: wangxiaogang@hit.edu.com

    哈尔滨建成集团研究员. 主要研究方向为飞行器总体设计.E-mail: yongzhi0451@sina.com

    北京航天晨信科技有限公司研究员. 主要研究方向为优化设计与系统工程.E-mail: zhaoyali_aerocim@163.com

    哈尔滨工业大学教授. 主要研究方向为飞行器制导、导航与控制.E-mail: cui_naigang@163.com

摘要: 考虑气动、轨迹、约束、指标间的耦合关系, 以多高超声速飞行器同时到达为目标建立了协同规划模型; 设计了一种自动满足终端约束的全新滑翔飞行剖面, 减少了规划算法需要处理的约束数量; 推导了滑翔段高精度解析解, 实现了过程约束和性能指标的快速求解; 提出了一种改进粒子群优化(PSO)算法, 借助强化学习方法构建协同需求与惯性权重间的动态映射网络, 提高了在线规划效率. 最后通过数学仿真验证了方法的正确性和有效性.

English Abstract

周宏宇, 王小刚, 单永志, 赵亚丽, 崔乃刚. 基于改进粒子群算法的飞行器协同轨迹规划. 自动化学报, 2020, 46(x): 1−7. doi: 10.16383/j.aas.c190865
引用本文: 周宏宇, 王小刚, 单永志, 赵亚丽, 崔乃刚. 基于改进粒子群算法的飞行器协同轨迹规划. 自动化学报, 2020, 46(x): 1−7. doi: 10.16383/j.aas.c190865
Zhou Hong-Yu, Wang Xiao-Gang, Shan Yong-Zhi, Zhao Ya-Li, Cui Nai-Gang. Synergistic path planning for multiple vehicles based on an improved particle swarm optimization method. Acta Automatica Sinica, 2020, 46(x): 1−7. doi: 10.16383/j.aas.c190865
Citation: Zhou Hong-Yu, Wang Xiao-Gang, Shan Yong-Zhi, Zhao Ya-Li, Cui Nai-Gang. Synergistic path planning for multiple vehicles based on an improved particle swarm optimization method. Acta Automatica Sinica, 2020, 46(x): 1−7. doi: 10.16383/j.aas.c190865
  • 现代防空反导体系已具备全天域、多层次、多维度防御能力; 在此背景下, 传统弹道弹和巡航弹已无法满足复杂对抗态势下的作战需求. 因此人们开始研究具有多空域机动能力的临近空间高超声速滑翔飞行器. 美国于2003年便提出了HTV-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]分别用遗传算法和PSO算法实现了无人机编队控制, Xu[13]提出了一种用于故障条件下多巡航弹协同控制的时变容错算法.

    现有文献有效解决了滑翔段轨迹设计问题, 但协同规划方面的研究多以巡航弹/无人机为研究对象. 针对上述问题, 从提高在线规划效率出发, 提出一种自动满足终端约束的滑翔段飞行剖面. 推导滑翔段解析解, 实现了飞行剖面的快速重构. 提出一种基于惯性权重智能选择的改进PSO算法, 用于完成在线协同规划. 最后通过数学仿真验证了方法的正确性和有效性.

    • 高超声速飞行器滑翔段运动方程如下:

      $$\left\{ \begin{array}{l} \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{array} \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{\rm{ = }}\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}{\rm{ = }}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{\rm{ = }}0.156,{\beta _1}{\rm{ = }}2.6 \times {10^{ - 4}},{\beta _2} = 1.5 \times {10^{ - 6}}$ .

      定义从当前位置到终端位置间的航程为剩余航程, 记为 ${R_L}$ , 则有:

      $$\left\{ \begin{array}{l} {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{array} \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{array}{l} {{{c}}_{{\rm{end}}}} = {\left[ {{h_f} - {h_{{\rm{ter}}}},{\gamma _f} - {\gamma _{{\rm{ter}}}},{R_{Lf}},{\zeta _f}} \right]^{\rm{T}}} = {\bf{0}} \\ {{{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{array} \right.$$ (6)

      式中: ${{{c}}_{{\rm{end}}}}$ 为终端约束, ${{{c}}_{{\rm{path}}}}$ 为过程约束. 下标“ $f$ ”表示实际终端状态, 下标“ter”表示期望终端状态. 下标“max”表示该过程变量允许的最大值, ${N_m}$ 为法向过载的绝对值.

    • 设协同飞行器数量为 ${N_{{\rm{hpy}}}}.$ 协同规划应应考虑各个飞行器的总体性能以及各自飞行约束, 记为 ${{c}}_{{\rm{end}}}^{\left( 1 \right)}, \cdots ,{{c}}_{{\rm{end}}}^{\left( {{N_{hyp}}} \right)},{{c}}_{{\rm{path}}}^{\left( 1 \right)}, \cdots ,{{c}}_{{\rm{path}}}^{\left( {{N_{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{hpy}}}}$$ (7)
    • 设飞行高度为剩余航程的函数:

      $$h = f\left( {{R_L}} \right) = - \frac{1}{{\bar \beta }}\ln F\left( {{R_L}} \right) = - \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'\left( {{R_L}} \right)}}{{\bar \beta F\left( {{R_L}} \right)}} = - \frac{{\tan \gamma }}{{\cos \zeta }} $$ (9)
      $$\frac{{{{\rm{d}}^2}h}}{{{\rm{d}}R_L^2}} = \frac{{{{\left[ {F'\left( {{R_L}} \right)} \right]}^2} - F''\left( {{R_L}} \right)F\left( {{R_L}} \right)}}{{\bar \beta {F^2}\left( {{R_L}} \right)}} = - \frac{{{\rm{d}}\left( {\frac{{\tan \gamma }}{{\cos \zeta }}} \right)}}{{{\rm{d}}{R_L}}}$$ (10)

      式中: $F'({R_L})\!{\rm{ = }}\!\displaystyle\sum\nolimits_{i = 1}^5 {i{a_i}R_L^{i - 1}},F''({R_L})\!{\rm{ = }}\!\displaystyle\sum\nolimits_{i = 2}^5 i\left( {i \!-\! 1} \right)$ ${a_i}R_L^{i - 2}$ .

      定义 ${K_h} = {{\rm{d}}^2}h{\rm{/d}}R_L^2.$ 认为 ${\dot \psi _w} = 0,\dot \zeta = - \dot \psi; $ 由式(11)可得升力为:

      $$L = \dfrac{{{K_h}{{\cos }^2}\gamma {{\cos }^2}\zeta + \dfrac{g}{{{V^2}}} - \dfrac{1}{r} + \dfrac{{{K_b}}}{r}}}{{(\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})$ 需建立六个方程. 由于初始和终端高度已知, 由式(8)可得:

      $$\left\{ \begin{array}{l} {{\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{array} \right.$$ (12)

      式中: 下标“0”表示滑翔段初始运动状态.

      由于初始和终端飞行路径角已知, 认为 $\cos \zeta \approx 1,$ 由式(9)可得:

      $$\left\{ \begin{array}{l} \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{array} \right.$$ (13)

      ${h_1},{h_2}$ 分别为剩余航程2/3处(记为 ${R_L}_1)$ 和1/3处(记为 ${R_L}_2)$ 的高度值; 由式(8)可得:

      $$\left\{ \begin{array}{l} {{\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{array} \right.$$ (14)

      定义剖面设计变量 ${{{u}}_h} = {[{h_{{\kern 1pt} 1}},{h_2}]^{\rm{T}}}$ , 则式(7)对应的优化变量分别为 ${{u}}_{{\rm{syn}}}^{\left( 1 \right)} = [h_1^{\left( 1 \right)}, \cdots ,h_1^{({N_{{\rm{hpy}}}})}]$ ${{u}}_{{\rm{syn}}}^{\left( 2 \right)} = [h_2^{\left( 1 \right)}, \cdots ,h_2^{({N_{{\rm{hpy}}}})}]$ , 共 $2{N_{{\rm{hyp}}}}$ 个.

      最后, 设计倾侧角变化规律如下:

      $$\sigma = \zeta {\sigma _{\max }}/{\zeta _{\max }}$$ (15)

      式中: ${\zeta _{\max }}$ 为给定正数.

    • 滑翔段高度的解析解即式(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}{\rm{ = }}\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{\rm{ = }}0$ ${R_L}{\rm{ = }}{R_L}_0,$ $t{\rm{ = }}{t_{{\rm{go}}}}$ ${R_{Lf}}{\rm{ = }}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)
    • PSO算法用 ${M_{{\rm{pso}}}}$ 个粒子寻找 ${D_{{\rm{pso}}}}$ 个参数. 设粒子 $i$ 在空间中的位置和速度分别为 ${{{x}}_i} \!\!=\!\! {[{x_{i,1}}, \cdots ,{x_{i,{D_{{\rm{pso}}}}}}]^{\rm{T}}},$ ${{{v}}_i} = {[{v_{i,1}}, \cdots ,{v_{i,{D_{{\rm{pso}}}}}}]^{\rm{T}}},$ 其历史最佳位置为 ${{{p}}_i} = [{p_{i,1}},$ $\cdots ,{p_{i,{D_{{\rm{pso}}}}}}]^{\rm{T}};$ 整个群体的历史最佳位置为 ${{{p}}_g} = [{p_{g,1}},$ $\cdots ,{p_{g,{D_{{\rm{pso}}}}}}]^{\rm{T}},$ 相应性能指标为 ${g_{{\rm{best}}}}.$ PSO算法如下[14-15]:

      $$v_{i,j}^{\left( {\ell + 1} \right)}{\rm{ = }}wv_{i,j}^{\left( \ell \right)}{\rm{ + }}{c_1}{r_1}\left( {{p_{i,j}} - x_{i,j}^{\left( \ell \right)}} \right){\rm{ + }}{c_2}{r_2}\left( {{p_{g,j}} - x_{i,j}^{\left( \ell \right)}} \right)$$ (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}{\rm{ = }}{c_2}{\rm{ = }}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)避免了对粒子速度进行初始化及更新, 算法结构更简洁, 利于提高计算效率.

      调整权重是控制粒子寻优过程的常用方法. 借助强化学习方法, 通过快速智能选择惯性权重实现用较少的粒子数和迭代次数找出满足协同规划参数( ${{u}}_{{\rm{syn}}}^{\left( 1 \right)}$ ${{u}}_{{\rm{syn}}}^{\left( 2 \right)}$ ).

    • 设惯性权重为进化代数的分段函数:

      $$w\left( \ell \right) = \left\{ {\begin{array}{*{20}{l}} {{w_0} + \dfrac{{{w_1} - {w_0}}}{{{\ell _{\max }}{\rm{/3}}}}\ell }&{0 \le \ell < \dfrac{{{\ell _{\max }}}}{3}}\\ {{w_1} + \dfrac{{{w_2} - {w_1}}}{{{\ell _{\max }}{\rm{/3}}}}\ell }&{\dfrac{{{\ell _{\max }}}}{3} \le \ell < \dfrac{{2{\ell _{\max }}}}{3}}\\ {{w_2} + \dfrac{{{w_3} - {w_2}}}{{{\ell _{\max }}{\rm{/3}}}}\ell }&{\dfrac{{2{\ell _{\max }}}}{3} \le \ell < {\ell _{\max }}} \end{array}} \right.$$ (30)

      式中: ${w_0} \sim {w_3}$ 为未知量, ${\ell _{\max }}$ 为最大进化代数.

      定义自变量 ${{{u}}_w} = {[{w_0},{w_{{\kern 1pt} 1}},{w_2},{w_3}]^{\rm{T}}}$ 以及动作 ${{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}$ 为预置系数.

      设状态值 ${{s}} = {{\rm{[}}{p_{g,1}}{\rm{,}} \cdots {\rm{,}}{\kern 1pt} {p_{g,{D_{{\rm{pso}}}}}}{\rm{]}}^{\rm{T}}}$ 以及观测值 ${{{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}}}}$ 来表征状态 ${{s}}$ 和动作 ${{d}}$ 的优劣, 其中优劣性由奖励值 $r = {g_{{\rm{best}}}}$ 表示. 定义策略网络 ${P_{{\rm{net}}}}$ 用以表征动作 ${{d}}$ 的更新; 由式(35)可知, ${P_{{\rm{net}}}}$ 即为 ${{{u}}_w}$ 的求解器. 根据上述定义, 记 ${Q_{{\rm{net}}}}$ ${P_{{\rm{net}}}}$ 的值分别为 ${\theta _Q}$ ${{{\theta }}_P},$ 提出最优惯性权重训练方法如下:

      1) 随机初始化 ${\theta _Q}$ ${{{\theta }}_P},$ 并令 ${\theta '_Q} = {\theta _Q},{{{\theta}} '_P} = {{{\theta }}_P};$ $k = 1,$ 随机初始化 ${{{u}}_w}$ ;

      2) 根据 ${{{u}}_w}$ 执行PSO, 得到状态 ${{{s}}_k}$ ;

      3) 随机初始化扰动量 $\Delta {{{d}}_k},$ 得到动作 ${{{d}}_k}{\rm{ = }}{P_{{\rm{net}}}}({{{\theta}} _P},$ ${{{s}}_k}) + \Delta {{{d}}_k};$ 执行动作 ${{{d}}_k}$ 并更新 ${{{u}}_w},$ 运行PSO得到奖励 ${r_k}$ 和状态 ${{{s}}_k};$

      4) 将 $({{{s}}_k},{{{d}}_k},{r_k},{{{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} \cdot {\theta '_Q}\left( {{{{s}}_i},{{{d}}_i},{{{\theta}} _P}} \right)$$ (32)

      6) 以损失函数最小为目标更新 ${\theta _Q}$ :

      $$\min \;\;\;\frac{1}{{{N_R}}}\sum\limits_{i = 1}^{{N_R}} {{{\left( {{y_i} - {\theta _Q}\left( {{{{s}}_i},{{{d}}_i},{{{\theta}} _P}} \right)} \right)}^2}} \Rightarrow {\theta _Q}$$ (33)

      7) 通过梯度求解更新 ${{{\theta}} _P} = {{{\theta}} _P} + \Delta {{{\theta}} _P}$ :

      $$\Delta {{{\theta}} _P} = \sum\limits_{i = 1}^{{N_R}} {\left( {\frac{{\partial {\theta _Q}\left( {{{{d}}_i},{{{\theta}} _P},{{{s}}_i}} \right)}}{{\partial {{{d}}_i}}} \cdot \frac{{\partial {{{\theta}} _P}}}{{\partial {{{s}}_i}}}} \right)} $$ (34)

      8) 给定系数 ${\tau _\theta },$ ${\theta '_Q} = {\tau _\theta }{\theta _{^Q}} + (1 - {\tau _\theta }){\theta '_Q},{{{\theta}} '_P} =$ $ {\tau _\theta }{{{\theta}} '_P} + (1 - {\tau _\theta }){{{\theta}} '_P}$ .

    • 设从当前位置到起点的航程为 ${R_{{\rm{cov}}}},$ 更新周期为 $\Delta {R_{{\rm{cov}}}}.$ 图1所示, 在线协同规划具体步骤如下:

      图  1  在线规划算法流程图

      Figure 1.  The flowchart of the planning algorithm

      1) 协同飞行剖面重构

      计算观测值 ${{{o}}_t}$ 并和上一步获得的 ${{u}}_{{\rm{syn}}}^{\left( 1 \right)}$ ${{u}}_{{\rm{syn}}}^{\left( 2 \right)}$ (如果是第一步则随机初始化 ${{u}}_{{\rm{syn}}}^{\left( 1 \right)}$ ${{u}}_{{\rm{syn}}}^{\left( 2 \right)}$ )共同输入训练好的策略函数网络, 输出惯性权重; 然后用PSO计算 ${{u}}_{{\rm{syn}}}^{\left( 1 \right)}$ ${{u}}_{{\rm{syn}}}^{\left( 2 \right)}$ 的值, 使指标函数式(7)达到最小;

      2) 个体飞行剖面重构

      2a) 计算 ${R_{L1}},{R_{L2}}$ 及其对应高度值 ${h_1},{h_2};$

      2b) 求解系数 ${a_0} \sim {a_5}$ , 更新飞行剖面.

      3) 计算攻角和倾侧角, 更新运动状态;

      4) 计算航程 ${R_{{\rm{cov}}}};$ 根据 ${R_{{\rm{cov}}}},$ 每五个周期 $(5\Delta {R_{{\rm{cov}}}})$ 执行一次协同剖面重构, 每一个周期 $(\Delta {R_{{\rm{cov}}}})$ 执行一次个体剖面重构; 未到达更新周期时, 各飞行器只执行步骤3.

      循环步骤 $1\sim 4 $ , 直至飞行器到达目标点上空. 步骤2用于控制终端高度和飞行路径角, 但未考虑过程约束; 步骤1既要保证时间协同, 又要通过优化求解来满足过程约束, 同时还负责修正解析解与实际值间的偏差.

    • 参考相关文献[16-17], 设滑翔段过程约束为 ${q_{\max }}= $ ${\rm{ 60\;kPa}},{N_{\max }}{\rm{ = 3}}{\rm{.0}},{\alpha _{\max }}{\rm{ = 25^\circ }},{\sigma _{\max }}{\rm{ = 50}}^\circ . $ 考虑对四架总体参数不同的高超声速飞行器进行协同规划; 参考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 {\rm{ = }}{\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

      图  2  滑翔段飞行剖面

      Figure 2.  Flight profiles of different vehicles

      图  3  飞行路径角随时间变化情况

      Figure 3.  Time histories of the flight path angle

      图  4  过程约束随时间变化情况

      Figure 4.  Time histories of the path constraints

      图  5  攻角和倾侧角随时间变化情况

      Figure 5.  Histories of angle-of-attack and the bank angle

      若不进行协同规划, 各飞行剖面设计参数均取 ${h_1}$ =42 km, ${h_2}$ =38 km, 则各飞行器的到达时刻为694 s, 703 s, 733 s和744 s, 说明协同规划方法正确有效.

      考虑初始状态下第一次计算协同参数( ${{u}}_{{\rm{syn}}}^{\left( 1 \right)}$ ${{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σ值
      1 初始速度(m/s) 30
      2 初始飞行路径角(°) 0.2
      3 初始高度(m) 500
      4 初始航向偏差(°) 0.2
      5 初始剩余航程(km) 50
      6 气动系数(%) 10
      7 大气密度(%) 10

      以飞行器1为例, 给出打靶仿真结果如图6 $\sim $ 7所示; 结果表明各种约束条件依然能够得到满足. 在高度和飞行路径角的解析解推导过程中未为用到近似假设条件, 因此可忽略它们和真实值间的误差. 对于速度和剩余时间的解析解, 相应误差会随着时间而减小, 并且剩余时间是根据实时运动状态估计的, 因此不存在累计误差.

      图  6  攻角随时间变化情况(飞行器1)

      Figure 6.  Histories of the angle-of-attack (vehicle-1)

      图  7  动压随时间变化情况(飞行器1)

      Figure 7.  Histories of the dynamic pressure (vehicle-1)

    • 建立了多高超声速滑翔飞行器协同规划问题数学模型. 提出了一种新的滑翔飞行剖面, 实现了终端约束的自动满足, 并推导了滑翔段解析运动方程. 提出基于强化学习方法的改进PSO算法, 提高了在线协同规划效率. 仿真结果表明, 在初始状态和总体参数各异的情况下, 多弹时间协同偏差不超过2 s, 并且能够克服干扰因素的影响.

WeChat 关注分享

返回顶部

目录

    /

    返回文章
    返回