CN108519958B - 一种解析构造航天器追逃界栅和判断捕获逃逸区域的方法 - Google Patents

一种解析构造航天器追逃界栅和判断捕获逃逸区域的方法 Download PDF

Info

Publication number
CN108519958B
CN108519958B CN201810109434.3A CN201810109434A CN108519958B CN 108519958 B CN108519958 B CN 108519958B CN 201810109434 A CN201810109434 A CN 201810109434A CN 108519958 B CN108519958 B CN 108519958B
Authority
CN
China
Prior art keywords
spacecraft
equation
escape
boundary
differential
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
CN201810109434.3A
Other languages
English (en)
Other versions
CN108519958A (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 CN201810109434.3A priority Critical patent/CN108519958B/zh
Publication of CN108519958A publication Critical patent/CN108519958A/zh
Application granted granted Critical
Publication of CN108519958B publication Critical patent/CN108519958B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13Differential equations
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Operations Research (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Navigation (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)
  • Slot Machines And Peripheral Devices (AREA)

Abstract

本发明公开了一种解析构造航天器追逃界栅和判断捕获逃逸区域的方法,通过用微分对策方法构造初始模型,求解初始模型,得到界栅解析解从而解析构造航天器追逃界栅。本发明还通过对界栅进行快速计算和分析,明确判断航天器处在追逃博弈中的捕获区还是逃逸区,为空间威胁预警和博弈路径规划提供参考,设计方法正确合理,计算过程快速有效,对实际任务的适用性好。

Description

一种解析构造航天器追逃界栅和判断捕获逃逸区域的方法
技术领域
本发明涉及航天器追逃博弈路径规划方法,尤其涉及一种航天器追逃界栅的构造及捕获逃逸区域判断方法。
背景技术
空间交会规避机动是空间安全的一项关键技术。随着空间资源的日益紧张,在未来可能的空间对抗中,航天器将面临被追踪航天器捕获或杀伤的威胁。当对抗两航天器均能自主决策和机动时,该对抗便成为一个空间连续动态博弈的问题。
微分对策是一种用微分方法研究连续时间无限动态博弈的理论,非常适合对这类问题进行建模和分析。在应用该理论定性分析追逃的结局能否实现时,涉及到对捕获区、逃逸区以及二者的分界面——界栅的构造求解。而航天器追逃界栅和捕获逃逸区域分析方法则是将微分对策理论应用于航天器追逃博弈问题里的成果。对微分对策及界栅概念的详细介绍可以参看李登峰所著的《微分对策及其应用》一书。
构造求解界栅的方法通常可分为数值方法和解析方法。而目前对于航天器追逃界栅的求解基本上采用的都是数值方法。其原因之一是对策模型复杂:对抗在三维空间进行,同时进行追逃博弈的两航天器运动需要满足轨道动力学规律;二是对策求解困难:基于微分对策理论求解航天器轨道追逃的最优控制策略时,通常需要求解复杂的微分方程组。
但数值方法相比解析方法不便于快速简便地分析和计算,更不利于直接反映追逃博弈的运动规律。此外,采用数值方法求得的界栅应用起来也有一定局限性,对捕获区和逃逸区的分析更多止于两种区域的划分,不适合针对特定对策状态判断航天器当前所处区域是捕获区还是逃逸区。
发明内容
为解决上述问题,本发明特提出一种用解析方法构造求解航天器追逃界栅和判断捕获逃逸区域的方法。该种解析方法能够快速计算求解追逃界栅,直观地反映追逃博弈的运动规律,进而明确判别航天器在追逃博弈中处于捕获区还是逃逸区,其结果可为空间威胁预警和博弈路径规划提供有效的参考。
本发明的技术方案和实施步骤如下:
一种解析构造航天器追逃界栅的方法,包括以下步骤:
S1,用微分对策方法构造初始模型;
进一步的,所述S1包括以下步骤:
S101,基于C-W方程构造微分对策的Hamilton函数,
1)建立LVLH(Local Vertical Local Horizontal)坐标系,构造基于C-W方程的运动状态方程;
进一步的,以追逃两航天器附近的一颗虚拟航天器作为参考航天器建立当地轨道坐标系,该坐标系的原点位于参考航天器的质心o,ox轴沿参考航天器的径向,oz轴沿参考航天器轨道面的法向,oy轴沿参考航天器运动的轨迹切向,并与ox、oz轴构成右手坐标系;使用下式所示的C-W方程描述逃逸航天器或追踪航天器相对虚拟航天器的运动:
Figure BDA0001568749090000021
式中,ω为虚拟航天器做圆周运动的角速度,ax、ay、az为逃逸或追踪航天器分别在径向、迹向和法向推力加速度分量,若设追踪航天器和逃逸航天器的相对运动状态分别为
Figure BDA0001568749090000022
(本文中的下标P和E分别表示追踪航天器和逃逸航天器)且均满足式给出的C-W方程;令
Figure BDA0001568749090000023
为航天器追逃微分对策的状态变量,aE.x、aE.y、aE.z、aP.x、aP.y、aP.z分别为两航天器三个方向的推力加速度分量,则由式(1)可得式:
Figure BDA0001568749090000024
设两航天器均采用连续推力控制,最大加速度大小分别为Tp、TE,推力加速度方向用α和β表示,分别为航天器推力加速度方向的偏航角(Ti在xoy平面投影和x轴夹角)和俯仰角(Ti和xoy平面夹角),αi∈[0,2π]、
Figure BDA0001568749090000025
(下标i=P or E),则可得微分对策的运动状态方程:
Figure BDA0001568749090000026
2)构造微分对策的Hamilton函数;
将式(3)写成状态空间的表达式:
Figure BDA0001568749090000027
式中A、B为状态方程的系数矩阵,U=TE-Tp
构造微分对策的Hamilton函数:
Figure BDA0001568749090000028
式中λ=[λ123456]T为协态向量;
S102,求解微分对策的协态方程和最优控制方程;
1)求解微分对策的协态方程;
协态方程可以写为:
Figure BDA0001568749090000031
其解析表达式为:
λ(t)=Φλ(t,t0)λ(t0) (7)
式中状态转移矩阵Φλ(t,t0)满足:
Figure BDA0001568749090000037
设tf为对策结束时刻,记τ=tf-t为剩余捕获时间,可得:
Figure BDA0001568749090000032
假设追踪航天器的捕获半径为R0,则该微分对策的终端目标集边界为:
Figure BDA0001568749090000033
则终端时刻的横截条件为:
Figure BDA0001568749090000034
式中μ为乘子变量,为实数,在该边界下,优选的0.5;由上式可得协态量的末值条件:
Figure BDA0001568749090000035
式中
Figure BDA0001568749090000036
为目标集边界上一点相对坐标系的偏航角与俯仰角;
(2)求解微分对策的最优控制方程;
最优控制方程即为该微分对策的鞍点,所谓鞍点,就是泛函的极值点(除掉最大值点及最小值点),当微分对策的Hamilton函数为追逃双方控制量的连续可导函数,鞍点可用公式表达为:
Figure BDA0001568749090000041
将Hamilton函数代入上式得最优控制方程如下:
Figure BDA0001568749090000042
S2,求解初始模型,得到界栅解析解;
S201,据协态方程和末值条件,求解协态量和最优控制量的解析式;
进一步的,求解协态量的解析表达式为:
Figure BDA0001568749090000043
将协态量的解析式代入并在目标集边界τ=0附近进行泰勒展开,综合考虑计算的工作量和精度的,优选的,保留至二阶项,可得追逃两航天器最优控制量与剩余捕获时间τ的解析表达式:
Figure BDA0001568749090000044
S202,将最优控制量的解析式代回C-W方程积分,求解航天器追逃界栅解析表达式;
具体的,将最优控制量的解析式代回式进行积分求解,可得界栅关于剩余捕获时间τ的解析表达式:
Figure BDA0001568749090000051
式中各系数如下:
Figure BDA0001568749090000052
其中:
Figure BDA0001568749090000053
特别的,当两航天器在共面圆轨道追逃博弈时,界栅解析解可表达为:
Figure BDA0001568749090000054
式中各系数如下:
Figure BDA0001568749090000061
其中:
Figure BDA0001568749090000062
至此,完成界栅的构造求解。
本发明同时公布了一种基于界栅解析解判断航天器捕获逃逸区域方法,具体方案如下:
S1,用微分对策方法构造初始模型;
S2,求解初始模型,得到界栅解析解;
进一步的,给定两航天器博弈情景中的参数,初始化界栅的解析解模型,即将参数带入下式:
Figure BDA0001568749090000063
Figure BDA0001568749090000064
Figure BDA0001568749090000071
S3,用流程方法判断航天器所处区域;
进一步的,所述S3包括以下步骤:
S301:给定对策的当前状态
Figure BDA0001568749090000072
和最长博弈时间τS,令θ=0°、
Figure BDA0001568749090000073
τ=τS
S302:由式根据x=xS、y=yS、z=zS求解出另外三个参数R、S5和S4
S303:由式根据上步求解得到的一组参数,计算得到
Figure BDA0001568749090000074
Figure BDA0001568749090000075
S304:保持其余参数不变,对每一个τ∈[0,τS],分别计算出对应的
Figure BDA0001568749090000076
Figure BDA0001568749090000077
并从中找到使:
Figure BDA0001568749090000078
的τ,若存在这样的τ且此时R>0,则进入S306,否则进入S307;
S305:返回S301,逐步增大θ和
Figure BDA0001568749090000079
并重复S302~S304,直至两个方位角达到360°;
S306:对得到的解,若存在多组满足式的解,取其中R最小的那组解;
S307:根据得到的解对捕获逃逸区域进行判断:如果R>R0,则对策当前状态处于逃逸区内;如果R<R0,则对策当前状态处于捕获区内;如果R=R0,则对策当前状态恰好处于界栅上;若不存在这样的一组解,则对策当前状态处于逃逸区内;
至此,完成了航天器处于捕获逃逸区域的判断。
本发明一种解析构造航天器追逃界栅和判断捕获逃逸区域的方法的优点如下:
1.本发明采用解析方法构造求解航天器的追逃界栅,并基于解析法求得的界栅明确判断航天器当前状态所在区域是捕获区还是逃逸区。经验证,这种解析方法在计算结果上精度较高,与数值方法结果吻合度较好,却更为简便快捷,更好地反映博弈运动规律。这种对航天器所在区域的判断方法,相比于对捕获逃逸区域的划分方法,结果更加明确、实用。
2.特别地,航天器的追逃博弈更可能是在轨实时进行,由于星载计算机的计算能力相对较弱,因此解析计算追逃界栅和判断捕获逃逸区域的方法,更适合星载计算机快速准确地预警和求解博弈机动策略,对于实际工程任务具有更好的适用性。
附图说明
图1为本发明方法中构造求解界栅解析解的技术方案图。
图2为本发明方法中判断捕获逃逸区域的技术方案图。
图3为本发明实施例一的精度对比图。
图4为本发明实施例二的精度对比图。
图5为本发明实施例三的精度对比图。
图6为本发明实施例四的精度对比图。
图7为本发明实施例五的基本流程图。
图8为本发明实施例五的追逃验证图。
具体实施方式
实施例一:
如图1所示,本实施例一种解析构造航天器追逃界栅方法的步骤包括:
S1,用微分对策方法构造初始模型;
S101,基于C-W方程构造微分对策的Hamilton函数,
1)建立LVLH(Local Vertical Local Horizontal)坐标系,构造基于C-W方程的运动状态方程;
以追逃两航天器附近的一颗虚拟航天器作为参考航天器建立当地轨道坐标系,该坐标系的原点位于参考航天器的质心o,ox轴沿参考航天器的径向,oz轴沿参考航天器轨道面的法向,oy轴沿参考航天器运动的轨迹切向,并与ox、oz轴构成右手坐标系;使用下式所示的C-W方程描述逃逸航天器或追踪航天器相对虚拟航天器的运动:
Figure BDA0001568749090000081
式中,ω为虚拟航天器做圆周运动的角速度,ax、ay、az为逃逸或追踪航天器分别在径向、迹向和法向推力加速度分量,若设追踪航天器和逃逸航天器的相对运动状态分别为
Figure BDA0001568749090000082
且均满足式给出的C-W方程;令
Figure BDA0001568749090000083
为航天器追逃微分对策的状态变量,aE.x、aE.y、aE.z、aP.x、aP.y、aP.z分别为两航天器三个方向的推力加速度分量,则由式(1)可得:
Figure BDA0001568749090000084
设两航天器均采用连续推力控制,最大加速度大小分别为Tp、TE,推力加速度方向用α和β表示,分别为航天器推力加速度方向的偏航角(Ti在xoy平面投影和x轴夹角)和俯仰角(Ti和xoy平面夹角),αi∈[0,2π]、
Figure BDA0001568749090000085
(下标i=P or E),则可得微分对策的运动状态方程:
Figure BDA0001568749090000086
2)构造微分对策的Hamilton函数;
将式(3)写成状态空间的表达式:
Figure BDA0001568749090000087
式中A、B为状态方程的系数矩阵,U=TE-Tp
构造微分对策的Hamilton函数:
Figure BDA0001568749090000091
式中λ=[λ123456]T为协态向量;
S102,求解微分对策的协态方程和最优控制方程;
1)求解微分对策的协态方程;
协态方程可以写为:
Figure BDA0001568749090000092
其解析表达式为:
λ(t)=Φλ(t,t0)λ(t0) (7)
式中状态转移矩阵Φλ(t,t0)满足:
Figure BDA0001568749090000093
设tf为对策结束时刻,记τ=tf-t为剩余捕获时间,可得:
Figure BDA0001568749090000094
假设追踪航天器的捕获半径为R0,则该微分对策的终端目标集边界为:
Figure BDA0001568749090000095
则终端时刻的横截条件为:
Figure BDA0001568749090000096
式中μ=0.5;由上式可得协态量的末值条件:
Figure BDA0001568749090000101
式中
Figure BDA0001568749090000102
为目标集边界上一点相对坐标系的偏航角与俯仰角;
(2)求解微分对策的最优控制方程;
鞍点可用公式表达为:
Figure BDA0001568749090000103
将Hamilton函数代入上式得最优控制方程如下:
Figure BDA0001568749090000104
S2,求解初始模型,得到界栅解析解;
S201,据协态方程和末值条件,求解协态量和最优控制量的解析式;
本实施例中,求解协态量的解析表达式为:
Figure BDA0001568749090000105
将协态量的解析式代入并在目标集边界τ=0附近进行泰勒展开,保留至二阶项,可得追逃两航天器最优控制量与剩余捕获时间τ的解析表达式:
Figure BDA0001568749090000111
S202,将最优控制量的解析式代回C-W方程积分,求解航天器追逃界栅解析表达式:
Figure BDA0001568749090000112
式中各系数如下:
Figure BDA0001568749090000113
其中:
Figure BDA0001568749090000114
特别的,当两航天器在共面圆轨道追逃博弈时,界栅解析解可表达为:
Figure BDA0001568749090000121
式中各系数如下:
Figure BDA0001568749090000122
其中:
Figure BDA0001568749090000123
输入情景参数:
Figure BDA0001568749090000124
得到界栅分量随剩余捕获时间变化曲线,见图2中解析求解结果。
为验证解析求解结果计算的可靠性与精度,本实施例中也将解析方法求解得到的界栅结果与数值方法进行比较,见图3中数值求解结果。从图中可以看出,解析公式求解的结果和数值积分得到的结果吻合得较好,这表明本发明所提出的界栅求解方法精度较高。另外,从图中还可以看出,解析求解方法在方向上得到的位置和速度结果的精度比x方向上的要高,解析求解结果中精度最低的是方向x上的速度
Figure BDA0001568749090000125
随着剩余捕获时间τ的增大,
Figure BDA0001568749090000126
的计算结果发散较快。
实施例二:
本实施例与实施例一的步骤基本相同,其主要不同点为:本实施例中,捕获点方位角不再是给定值,而是作为自变量。因而通过计算可以得到不同捕获点方位角
Figure BDA0001568749090000131
下微分对策在τ=1200s时刻的界栅。
再比较数值积分的结果和解析求解的结果,见图4。从图中可以看出,方位角
Figure BDA0001568749090000132
仅对界栅y方向上位置的解析求解精度有一定影响,当
Figure BDA0001568749090000133
Figure BDA0001568749090000134
附近时误差明显,而在其它界栅状态量的求解上,方位角影响较小。
实施例三:
本实施例与实施例二的步骤基本相同,其主要不同点为:本实施例中,自变量是捕获点相对速度
Figure BDA0001568749090000135
图5中分别给出了不同捕获点相对速度
Figure BDA0001568749090000136
下计算得到的对策在τ时刻的界栅,并比较了数值积分结果和解析求解结果。从图中可以看出,对不同的
Figure BDA0001568749090000137
解析计算的结果与数值积分得到的结果均吻合较好。同时从图中还可以看出,在给定的范围内,界栅在x和y两个方向上的位置均随着
Figure BDA0001568749090000138
的增大而减小,而在两个方向上的速度则随之增大。
实施例四:
本实施例与实施例二的步骤基本相同,其主要不同点为:本实施例中,自变量是捕获半径R0
图6中分别给出了不同捕获半径下计算得到的对策在τ时刻的界栅,并比较了数值积分结果和解析求解结果。从图中可以看出,捕获半径R0对界栅的求解影响较小,且解析求解带来的误差随R0的变化也较小。
实施例五:
本实施例中,基于解析求解的界栅判断捕获逃逸区域方法的步骤包括:
S1,用微分对策方法构造初始模型,此步骤同实施例一;
S2,求解初始模型,得到界栅解析解:
给定两航天器博弈情景中的参数,初始化界栅的解析解模型,即将参数带入下式:
Figure BDA0001568749090000139
Figure BDA0001568749090000141
Figure BDA0001568749090000142
输入情景参数:
Figure BDA0001568749090000143
然后初始化界栅的解析解;
S3,用流程方法判断航天器所处区域,具体步骤的流程图见图7:
S301:给定对策的当前状态
Figure BDA0001568749090000144
和最长博弈时间τS,令θ=0°、
Figure BDA0001568749090000145
τ=τS
S302:由式根据x=xS、y=yS、z=zS求解出另外三个参数R、S5和S4
S303:由式根据上步求解得到的一组参数,计算得到
Figure BDA0001568749090000146
Figure BDA0001568749090000147
S304:保持其余参数不变,对每一个τ∈[0,τS],分别计算出对应的
Figure BDA0001568749090000148
Figure BDA0001568749090000149
并从中找到使:
Figure BDA00015687490900001410
的τ,若存在这样的τ且此时R>0,则进入S206,否则进入S207;
S305:返回S301,逐步增大θ和
Figure BDA00015687490900001411
并重复S302~S304,直至两个方位角达到360°;
S306:对得到的解,若存在多组满足式的解,取其中R最小的那组解;
S307:根据得到的解对捕获逃逸区域进行判断:如果R>R0,则对策当前状态处于逃逸区内;如果R<R0,则对策当前状态处于捕获区内;如果R=R0,则对策当前状态恰好处于界栅上;若不存在这样的一组解,则对策当前状态处于逃逸区内;
至此,完成了航天器处于捕获逃逸区域的判断。
下表给出了共面轨道上四种不同初始对策状态下捕获逃逸区域的计算结果,相对轨迹仿真结果见图8。
表1航天器追逃博弈捕获逃逸区域分析
Figure BDA0001568749090000151
从表中数据可以看出,第1、3种情况下对策状态处于捕获区内,而另外两种情况下的对策状态处于逃逸区内,其中第4种情况中在给定时间范围内不存在最小微分对策的解。

Claims (8)

1.一种解析构造航天器追逃界栅的方法,其特征在于,包括如下步骤:
S1,用微分对策方法构造初始模型;
所述S1包括以下步骤:
S101,基于C-W方程构造微分对策的Hamilton函数;
S102,求解微分对策的协态方程和最优控制方程;
S2,求解初始模型,得到界栅解析解;
所述S2包括以下步骤:
S201,据协态方程和末值条件,求解协态量和最优控制量的解析式;
S202,将最优控制量的解析式代回C-W方程积分,求解航天器追逃界栅解析表达式;
其中,所述S1包括以下步骤:
S101,基于C-W方程构造微分对策的Hamilton函数;
1)建立LVLH坐标系,LVLH是Local Vertical Local Horizontal的简称,构造基于C-W方程的运动状态方程;
以追逃两航天器附近的一颗虚拟航天器作为参考航天器建立当地轨道坐标系,该坐标系的原点位于参考航天器的质心o,ox轴沿参考航天器的径向,oz轴沿参考航天器轨道面的法向,oy轴沿参考航天器运动的轨迹切向,并与ox、oz轴构成右手坐标系;使用下式所示的C-W方程描述逃逸航天器或追踪航天器相对虚拟航天器的运动:
Figure FDA0003390340880000011
式中,ω为虚拟航天器做圆周运动的角速度,ax、ay、az为逃逸或追踪航天器分别在径向、迹向和法向推力加速度分量,若设追踪航天器和逃逸航天器的相对运动状态分别为
Figure FDA0003390340880000012
下标P和E分别表示追踪航天器和逃逸航天器,且均满足式给出的C-W方程;令
Figure FDA0003390340880000013
为航天器追逃微分对策的状态变量,aE.x、aE.y、aE.z、aP.x、aP.y、aP.z分别为两航天器三个方向的推力加速度分量,则由式(1)可得式(2):
Figure FDA0003390340880000014
设两航天器均采用连续推力控制,最大加速度大小分别为Tp、TE,推力加速度方向用α和β表示,分别为航天器推力加速度方向的偏航角和俯仰角,偏航角为Ti在xoy平面投影和x轴夹角,俯仰角为Ti和xoy平面夹角,αi∈[0,2π]、
Figure FDA0003390340880000015
下标i=P or E,则可得微分对策的运动状态方程:
Figure FDA0003390340880000021
2)构造微分对策的Hamilton函数;
将式(3)写成状态空间的表达式:
Figure FDA0003390340880000022
式中A、B为状态方程的系数矩阵,U=TE-Tp
构造微分对策的Hamilton函数:
Figure FDA0003390340880000023
式中λ=[λ123456]T为协态向量;
S102,求解微分对策的协态方程和最优控制方程;
1)求解微分对策的协态方程;
协态方程可以写为:
Figure FDA0003390340880000024
其解析表达式为:
λ(t)=Φλ(t,t0)λ(t0) (7)
式中状态转移矩阵Φλ(t,t0)满足:
Figure FDA0003390340880000025
设tf为对策结束时刻,记τ=tf-t为剩余捕获时间,可得:
Figure FDA0003390340880000031
假设追踪航天器的捕获半径为R0,则该微分对策的终端目标集边界为:
Figure FDA0003390340880000032
则终端时刻的横截条件为:
Figure FDA0003390340880000033
式中μ为乘子变量,为实数,由上式可得协态量的末值条件:
Figure FDA0003390340880000034
式中
Figure FDA0003390340880000035
为目标集边界上一点相对坐标系的偏航角与俯仰角;
2)求解微分对策的最优控制方程;
最优控制方程即为该微分对策的鞍点,当微分对策的Hamilton函数为追逃双方控制量的连续可导函数,鞍点可用公式表达为:
Figure FDA0003390340880000036
将Hamilton函数代入上式得最优控制方程如下:
Figure FDA0003390340880000041
2.根据权利要求1所述的一种解析构造航天器追逃界栅的方法,其特征在于,所述S102第1)步中,μ为0.5。
3.根据权利要求2所述的一种解析构造航天器追逃界栅的方法,其特征在于,所述S2包括以下步骤:
S201,据协态方程和末值条件,求解协态量和最优控制量的解析式;
求解协态量的解析表达式为:
Figure FDA0003390340880000042
将协态量的解析式代入式(15),并在目标集边界τ=0附近进行泰勒展开,可得追逃两航天器最优控制量与剩余捕获时间τ的解析表达式:
Figure FDA0003390340880000043
S202,将最优控制量的解析式代回C-W方程积分,求解航天器追逃界栅解析表达式;
将最优控制量的解析式代回式进行积分求解,可得界栅关于剩余捕获时间τ的解析表达式:
Figure FDA0003390340880000051
式中各系数如下:
Figure FDA0003390340880000052
其中:
Figure FDA0003390340880000053
至此,完成界栅的构造求解。
4.根据权利要求3所述的一种解析构造航天器追逃界栅的方法,其特征在于,所述S202步中,当两航天器在共面圆轨道追逃博弈时,界栅解析解可表达为:
Figure FDA0003390340880000061
式中各系数如下:
Figure FDA0003390340880000062
其中:
Figure FDA0003390340880000063
至此,完成界栅的构造求解。
5.一种基于界栅解析解判断航天器捕获逃逸区域方法,其特征在于,包括以下步骤:
S1,用微分对策方法构造初始模型;
S2,求解初始模型,得到界栅解析解;
S3,用流程方法判断航天器所处区域;
所述S3包括以下步骤:
S301:给定对策的当前状态
Figure FDA0003390340880000064
和最长博弈时间τS,令θ=0°、
Figure FDA0003390340880000065
τ=τS
S302:由式根据x=xS、y=yS、z=zS求解出另外三个参数R、S5和S4
S303:由式根据上步求解得到的一组参数,计算得到
Figure FDA0003390340880000066
Figure FDA0003390340880000067
S304:保持其余参数不变,对每一个τ∈[0,τS],分别计算出对应的
Figure FDA0003390340880000068
Figure FDA0003390340880000069
并从中找到使:
Figure FDA00033903408800000610
的τ,若存在这样的τ且此时R>0,则进入S306,否则进入S307;
S305:返回S301,逐步增大θ和
Figure FDA0003390340880000071
并重复S302~S304,直至两个方位角达到360°;
S306:对得到的解,若存在多组满足式的解,取其中R最小的那组解;
S307:根据得到的解对捕获逃逸区域进行判断:如果R>R0,则对策当前状态处于逃逸区内;如果R<R0,则对策当前状态处于捕获区内;如果R=R0,则对策当前状态恰好处于界栅上;若不存在这样的一组解,则对策当前状态处于逃逸区内;
至此,完成了航天器处于捕获逃逸区域的判断。
6.根据权利要求5所述一种基于界栅解析解判断航天器捕获逃逸区域方法,其特征在于,所述S2步中:给定两航天器博弈情景中的参数,初始化界栅的解析解模型,即将参数带入下式:
Figure FDA0003390340880000072
Figure FDA0003390340880000073
Figure FDA0003390340880000074
7.根据权利要求5或6任一项所述一种基于界栅解析解判断航天器捕获逃逸区域方法,其特征在于:所述S1包括以下步骤:
S101,基于C-W方程构造微分对策的Hamilton函数,
1)建立LVLH坐标系,LVLH是Local Vertical Local Horizontal的简称,构造基于C-W方程的运动状态方程;
以追逃两航天器附近的一颗虚拟航天器作为参考航天器建立当地轨道坐标系,该坐标系的原点位于参考航天器的质心o,ox轴沿参考航天器的径向,oz轴沿参考航天器轨道面的法向,oy轴沿参考航天器运动的轨迹切向,并与ox、oz轴构成右手坐标系;使用下式所示的C-W方程描述逃逸航天器或追踪航天器相对虚拟航天器的运动:
Figure FDA0003390340880000081
式中,ω为虚拟航天器做圆周运动的角速度,ax、ay、az为逃逸或追踪航天器分别在径向、迹向和法向推力加速度分量,若设追踪航天器和逃逸航天器的相对运动状态分别为
Figure FDA0003390340880000082
下标P和E分别表示追踪航天器和逃逸航天器,且均满足式给出的C-W方程;令
Figure FDA0003390340880000083
为航天器追逃微分对策的状态变量,aE.x、aE.y、aE.z、aP.x、aP.y、aP.z分别为两航天器三个方向的推力加速度分量,则由式(1)可得式:
Figure FDA0003390340880000084
设两航天器均采用连续推力控制,最大加速度大小分别为Tp、TE,推力加速度方向用α和β表示,分别为航天器推力加速度方向的偏航角和俯仰角,偏航角为Ti在xoy平面投影和x轴夹角,俯仰角为Ti和xoy平面夹角,αi∈[0,2π]、
Figure FDA0003390340880000085
下标i=P or E,则可得微分对策的运动状态方程:
Figure FDA0003390340880000086
2)构造微分对策的Hamilton函数;
将式(3)写成状态空间的表达式:
Figure FDA0003390340880000087
式中A、B为状态方程的系数矩阵,U=TE-Tp
构造微分对策的Hamilton函数:
Figure FDA0003390340880000088
式中λ=[λ123456]T为协态向量;
S102,求解微分对策的协态方程和最优控制方程;
1)求解微分对策的协态方程;
协态方程可以写为:
Figure FDA0003390340880000091
其解析表达式为:
λ(t)=Φλ(t,t0)λ(t0) (7)
式中状态转移矩阵Φλ(t,t0)满足:
Figure FDA0003390340880000092
设tf为对策结束时刻,记τ=tf-t为剩余捕获时间,可得:
Figure FDA0003390340880000093
假设追踪航天器的捕获半径为R0,则该微分对策的终端目标集边界为:
Figure FDA0003390340880000094
则终端时刻的横截条件为:
Figure FDA0003390340880000095
式中μ为乘子变量,为实数;由上式可得协态量的末值条件:
Figure FDA0003390340880000101
式中
Figure FDA0003390340880000102
为目标集边界上一点相对坐标系的偏航角与俯仰角;
2)求解微分对策的最优控制方程;
最优控制方程即为该微分对策的鞍点,鞍点可用公式表达为:
Figure FDA0003390340880000103
将Hamilton函数代入上式得最优控制方程如下:
Figure FDA0003390340880000104
8.根据权利要求7所述一种基于界栅解析解判断航天器捕获逃逸区域方法,其特征在于:所述S2包括以下步骤:
S201,据协态方程和末值条件,求解协态量和最优控制量的解析式;
求解协态量的解析表达式为:
Figure FDA0003390340880000105
将协态量的解析式代入式并在目标集边界τ=0附近进行泰勒展开,可得追逃两航天器最优控制量与剩余捕获时间τ的解析表达式:
Figure FDA0003390340880000111
S202,将最优控制量的解析式代回C-W方程积分,求解航天器追逃界栅解析表达式;将最优控制量的解析式代回式进行积分求解,可得界栅关于剩余捕获时间τ的解析表达式:
Figure FDA0003390340880000112
式中各系数如下:
Figure FDA0003390340880000113
其中:
Figure FDA0003390340880000114
至此,完成界栅的构造求解。
CN201810109434.3A 2018-02-05 2018-02-05 一种解析构造航天器追逃界栅和判断捕获逃逸区域的方法 Active CN108519958B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810109434.3A CN108519958B (zh) 2018-02-05 2018-02-05 一种解析构造航天器追逃界栅和判断捕获逃逸区域的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810109434.3A CN108519958B (zh) 2018-02-05 2018-02-05 一种解析构造航天器追逃界栅和判断捕获逃逸区域的方法

Publications (2)

Publication Number Publication Date
CN108519958A CN108519958A (zh) 2018-09-11
CN108519958B true CN108519958B (zh) 2022-02-08

Family

ID=63432789

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810109434.3A Active CN108519958B (zh) 2018-02-05 2018-02-05 一种解析构造航天器追逃界栅和判断捕获逃逸区域的方法

Country Status (1)

Country Link
CN (1) CN108519958B (zh)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110673486B (zh) * 2019-10-22 2021-03-09 北京航空航天大学 一种基于动态博弈理论的多航天器追逃控制方法
CN111679592B (zh) * 2020-06-22 2023-04-07 中国人民解放军国防科技大学 一种航天器追逃博弈闭环半实物仿真系统及方法
CN112001120B (zh) * 2020-08-24 2022-03-01 哈尔滨工业大学 一种基于强化学习的航天器对多拦截器自主规避机动方法
CN113552872B (zh) * 2021-01-28 2024-05-14 北京理工大学 追捕机器人不同速度下的追逃控制方法
CN113435000B (zh) * 2021-04-30 2023-10-31 北京理工大学 基于几何异构2对1博弈问题的界栅构建及战况判断方法
CN113190033B (zh) * 2021-05-20 2022-05-10 北京理工大学 一种航天器飞行博弈中可交会的快速判别方法
CN113282096B (zh) * 2021-06-04 2023-06-09 中国人民解放军战略支援部队航天工程大学 地球静止轨道博弈航天器相对位置非线性误差的控制方法
CN113282097B (zh) * 2021-06-04 2022-07-29 中国人民解放军战略支援部队航天工程大学 一种geo博弈航天器相对位置非球形摄动误差的控制方法
CN113960926B (zh) * 2021-10-18 2024-04-16 北京理工大学 一种气动捕获制导参数边界的自适应调节方法
CN115729112B (zh) * 2023-01-10 2023-04-07 中国人民解放军国防科技大学 一种航天器追踪-逃逸-防御三方博弈的制导方法及系统
CN116449714B (zh) * 2023-04-20 2024-01-23 四川大学 一种多航天器追捕博弈轨道控制方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102902274A (zh) * 2012-08-08 2013-01-30 空军工程大学航空航天工程学院 一种自适应加权微分对策制导方法
CN105574261A (zh) * 2015-12-15 2016-05-11 北京理工大学 一种月球借力约束的地月平动点转移轨道设计方法
CN106647287A (zh) * 2017-02-20 2017-05-10 南京航空航天大学 一种基于自适应动态规划的输入受限微分对策制导方法
CN106697333A (zh) * 2017-01-12 2017-05-24 北京理工大学 一种航天器轨道控制策略的鲁棒性分析方法
CN107329935A (zh) * 2017-07-11 2017-11-07 沈阳航空航天大学 一种空中态势帕累托攻防策略的求解方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105302156B (zh) * 2015-12-03 2018-01-30 上海新跃仪表厂 一种地面验证系统及追踪航天器的轨迹规划方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102902274A (zh) * 2012-08-08 2013-01-30 空军工程大学航空航天工程学院 一种自适应加权微分对策制导方法
CN105574261A (zh) * 2015-12-15 2016-05-11 北京理工大学 一种月球借力约束的地月平动点转移轨道设计方法
CN106697333A (zh) * 2017-01-12 2017-05-24 北京理工大学 一种航天器轨道控制策略的鲁棒性分析方法
CN106647287A (zh) * 2017-02-20 2017-05-10 南京航空航天大学 一种基于自适应动态规划的输入受限微分对策制导方法
CN107329935A (zh) * 2017-07-11 2017-11-07 沈阳航空航天大学 一种空中态势帕累托攻防策略的求解方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Complex Differential Games of Pursuit-Evasion Type with State Constraints, Part 2: Numerical Computation of Optimal Open-Loop Strategies;M.H.BREITNER 等;《JOURNAL OF OPTIMIZATION THEORY AND APPLICATIONS》;19931231;第78卷(第3期);第443-462页 *
空间飞行器在视线坐标系中的追逃界栅;张秋华 等;《航天控制》;20071231;第25卷(第1期);第26-29页 *
航天器追逃微分对策界栅解析构造方法;祝海 等;《第十届全国多体动力学与控制暨第五届全国航天动力学与控制学术会议》;20171231;第167页 *

Also Published As

Publication number Publication date
CN108519958A (zh) 2018-09-11

Similar Documents

Publication Publication Date Title
CN108519958B (zh) 一种解析构造航天器追逃界栅和判断捕获逃逸区域的方法
Ender et al. Systems-of-systems analysis of ballistic missile defense architecture effectiveness through surrogate modeling and simulation
CN109238287B (zh) 一种航天器逃逸路径规划方法及系统
McLain et al. A decomposition strategy for optimal coordination of unmanned air vehicles
US6529821B2 (en) Route planner with area avoidance capability
Lu et al. Rapid generation of accurate entry landing footprints
CN108759839B (zh) 一种基于态势空间的无人飞行器路径规划方法
CN110017832B (zh) 一种基于Gauss解群优选的短弧初轨确定方法
CN111486851B (zh) 航天器近距离相对运动三维避障轨迹规划方法和装置
Jha et al. Cooperative guidance and collision avoidance for multiple pursuers
CN114812569B (zh) 一种追逃博弈机动航天器相对状态估计方法、装置和设备
Ben-Asher et al. New proportional navigation law for ground-to-air systems
Morgan et al. Minimum energy guidance for aerodynamically controlled missiles
CN114911167A (zh) 一种航天器有限时间追逃博弈控制的解析求解方法与系统
Wang et al. Equal-collision-probability-curve method for safe spacecraft close-range proximity maneuvers
Du et al. The object-oriented dynamic task assignment for unmanned surface vessels
Lee et al. Weapon target assignment problem with interference constraints
CN110231619B (zh) 基于恩克法的雷达交接时刻预报方法及装置
Ming et al. A real-time monocular visual SLAM based on the bundle adjustment with adaptive robust kernel
CN108897029B (zh) 非合作目标近距离相对导航视觉测量系统指标评估方法
Miller et al. Simultaneous tracking of multiple ground targets from a multirotor unmanned aerial vehicle
Wang et al. Force-based delay compensation for hardware-in-the-loop simulation divergence of 6-DOF space contact
Zhang et al. Evaluation method of UAV cluster simulation models with different granularity based on sensitivity analysis
Günaydin et al. Response Surface Based Performance Analysis of an Air-Defense Missile System Against Maneuvering Targets
Schuppe Modeling and simulation: a Department of Defense critical technology

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