CN114384803A - 考虑地影约束的小推力最优变轨方法 - Google Patents

考虑地影约束的小推力最优变轨方法 Download PDF

Info

Publication number
CN114384803A
CN114384803A CN202210029706.5A CN202210029706A CN114384803A CN 114384803 A CN114384803 A CN 114384803A CN 202210029706 A CN202210029706 A CN 202210029706A CN 114384803 A CN114384803 A CN 114384803A
Authority
CN
China
Prior art keywords
constraint
optimal
thrust
shadow
spacecraft
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.)
Pending
Application number
CN202210029706.5A
Other languages
English (en)
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 CN202210029706.5A priority Critical patent/CN114384803A/zh
Publication of CN114384803A publication Critical patent/CN114384803A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/042—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Computation (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明公开的考虑地影约束的小推力最优变轨方法，属于航天器轨道动力学与控制领域。本发明实现方法：建立基于角度差判断的地影约束模型，将对航天器是否位于地影区的判断转化为对角度差符号的判断，不受航天器相对于地球方位限制；将地影约束添加在动力学中，通过最优控制建立考虑地影约束的小推力轨道最优控制模型；设计平滑函数和平滑参数，对不连续点平滑处理并使地影区推力大小随平滑参数的变化而相应变化，使得以平滑参数作为同伦参数的同伦过程顺利进行；结合同伦法和间接法求解考虑地影约束的小推力轨道最优控制问题，获得满足地影约束的小推力轨道最优控制和最优轨迹，实现满足地影约束的小推力最优变轨。

Description

考虑地影约束的小推力最优变轨方法
技术领域
本发明涉及考虑地影约束的小推力最优变轨方法，尤其涉及航天器小推力轨道最优控制和地影约束处理方法，属于航天器轨道动力学与控制领域。
背景技术
小推力以其高比冲、节省燃料的优势在航天器轨道机动中具有重要的应用价值。在小推力最优变轨中，建立小推力轨道最优控制模型，结合间接法和同伦法进行求解是一种有效的手段，但在考虑地影约束时存在一定困难。文献(Woollands R,Taheri E.Optimallow-thrust gravity perturbed orbit transfers with shadow constraints[C]//The2019AAS/AIAA Astrodynamics Specialist Conference.2019.)基于距离差判断的地影约束模型，采用双曲正切函数将地影约束进行平滑处理，结合间接法和同伦法求解了小推力燃料最优轨道转移问题。方法理论研究无误，但在实际小推力变轨中，小推力为航天器提供的加速度很小，使得变轨时间很长，航天器进出地影区次数较多，导致动力学系统中出现大量不连续点。此时，该文献中的方法受限于其地影约束模型和平滑函数，从不考虑地影约束的解到考虑地影约束的解间的同伦过程难以进行，无法完成对考虑地影约束的小推力最优变轨问题的求解。为解决这一问题，必须使用一种新的方法来平滑处理地影约束，从而在航天器多次进出地影区的情况下，仍能够有效地将地影约束进行平滑，并使以平滑参数为同伦参数的同伦过程顺利进行，实现考虑地影约束的小推力最优变轨。因此建立一种新的地影约束模型，设计更为有效的平滑函数，顺利求解考虑地影约束的小推力最优变轨问题，实现满足地影约束的小推力最优变轨，对于小推力在航天器轨道机动中的应用具有重要意义。
发明内容
为了解决航天器多次进出地影区情况下小推力最优变轨问题难以求解的问题，本发明的目的是提供考虑地影约束的小推力最优变轨方法，该方法以角度差作为判断依据对地影约束进行建模，适用于全部轨道段，在此基础上，对不连续点进行平滑处理的同时使地影区推力大小相应变化，保证同伦过程顺利进行，完成对地影约束下小推力变轨的优化，实现满足地影约束的小推力最优变轨，包括时间最优和燃料最优。本发明具有求解精确、适用性广、可靠性高、实现简便等优点。
本发明的目的是通过以下技术方案实现的。
本发明公开的考虑地影约束的小推力最优变轨方法，在圆锥形或圆柱形地影模型的基础上，建立基于角度差判断的地影约束模型，将对航天器是否位于地影区的判断转化为对角度差符号的判断，不受航天器相对于地球的方位限制，适用于全部轨道段；将地影约束添加在动力学中，通过最优控制原理建立考虑地影约束的小推力轨道最优控制模型；设计平滑函数和平滑参数，对不连续点进行平滑处理并使地影区推力大小随平滑参数的变化而相应变化，使得以平滑参数作为同伦参数的同伦过程顺利进行；结合同伦法和间接法求解考虑地影约束的小推力轨道最优控制问题，从而获得满足地影约束的小推力轨道最优控制和最优轨迹；根据所获得的最优控制和最优轨迹，实现满足地影约束的小推力最优变轨。所述最优包括时间最优和燃料最优。
本发明公开的考虑地影约束的小推力最优变轨方法，包括如下步骤：
步骤一：在圆锥形或圆柱形地影模型的基础上，建立基于角度差判断的地影约束模型，将对航天器是否位于地影区的判断转化为对角度差符号的判断，不受航天器相对于地球的方位限制，因此适用于全部轨道段。
针对圆锥形地影模型建立基于角度差判断的地影约束模型，具体实现方法如下：
将地球和太阳均视为球体，此时地影区域为圆锥形投影。S表示太阳中心，E表示地球中心，RS表示太阳半径，RE表示地球半径，SE所在直线为地影区域的中心线，V在中心线上表示本影区的顶点，αp为半影区圆锥半角，αu为本影区圆锥半角。航天器运行在轨道上，r为其在地心惯性坐标系下的位置矢量，r为其地心距，即位置矢量r的模，rS为太阳中心在地心惯性坐标系下的位置矢量。当航天器位于本影区或半影区时，均认为推力器无法工作，即不区分本影区和半影区。在某一时刻，航天器可能位于地影区外，也可能位于地影区内。以地心为圆心、r为半径作圆，该圆与地影区边界相交于以中心线为对称轴的两点，D点为与航天器位于同侧的交点。βsc为航天器位置矢量r与地影区中轴线间的夹角，满足
βsc∈[0 180°] (1)
βp为DE与地影区中心线间的夹角，有
βp＝αp+β (2)
其中
Figure BDA0003465798020000021
为辅助角，因为r＞RE，所以β∈(0 90°)。
令rS＝||rS||为日地距离，
Figure BDA0003465798020000022
为太阳在地心惯性坐标系下的方向向量，则
Figure BDA0003465798020000023
Figure BDA0003465798020000024
并将βsc和αp分别表示为
Figure BDA0003465798020000031
Figure BDA0003465798020000032
令ρp表示角量βsc和βp间的差值，用于表征航天器与地影区的位置关系，即
Figure BDA0003465798020000033
至此，在任意时刻，对航天器是否位于地影区的判断转化为对该时刻ρp的符号判断，即
Figure BDA0003465798020000034
式(7)(8)不受航天器相对于地球的方位限制，因此适用于全部轨道段。
针对圆柱形地影模型建立基于角度差判断的地影约束模型，具体实现方法如下：
将太阳光视为平行光，此时地影区域为圆柱形投影。E表示地球中心，RE表示地球半径，
Figure BDA0003465798020000035
为太阳在地心惯性坐标系下的方向向量，其所在直线为地影区域的中心线。航天器运行在轨道上，r为其在地心惯性坐标系下的位置矢量，r为其地心距。在某一时刻，航天器可能位于地影区外，也可能位于地影区内。以地心为圆心、r为半径作圆，该圆与地影区边界相交于以中心线为对称轴的两点，D点为与航天器位于同侧的交点。βsc为航天器位置矢量r与地影区中轴线间的夹角，满足式(1)(5)。βp为DE与地影区中心线间的夹角，有
Figure BDA0003465798020000036
则角量βsc和βp间的差值ρp为
Figure BDA0003465798020000037
至此，在任意时刻，对航天器是否位于地影区的判断转化为对该时刻ρp的符号判断，即式(8)。式(10)(8)不受航天器相对于地球的方位限制，因此适用于全部轨道段。
对航天器受到的地影约束描述为：当航天器位于地影区时，推力器无法工作，即地影约束使得原容许控制集U改变，记为Up，但等效为容许控制集不变，而在原推力幅值比u上添加系数up，得到考虑地影约束后的推力幅值比ue，即
Figure BDA0003465798020000038
其中系数up为
Figure BDA0003465798020000041
因此，将ρp称为地影约束开关函数。针对圆锥形地影模型，式(7)(8)(11)(12)即表述了基于角度差判断的地影约束模型；针对圆柱形地影模型，式(10)(8)(11)(12)即表述基于角度差判断的地影约束模型。
步骤二：将地影约束添加在动力学中，通过最优控制原理建立考虑地影约束的小推力轨道最优控制模型。
采用改进春分点根数描述航天器的轨道运动，改进春分点根数与经典轨道根数间的关系为
Figure BDA0003465798020000042
其中a为轨道半长轴，e为轨道偏心率，i为轨道倾角，Ω为升交点赤经，ω为近地点幅角，ν为真近点角，p为半通径，L为真经度。令状态变量x为
x＝[p f g h k L]T (14)
m为航天器质量，再根据Gauss型摄动方程，航天器的动力学方程为
Figure BDA0003465798020000043
其中u∈[0 1]为推力幅值比，up为地影约束系数，Tmax为最大推力幅值，Isp为推力器比冲，g0为标准重力加速度，α为推力方向向量。Tmax和Isp在航天器飞行过程中保持不变。式(15)中的M、D分别为
Figure BDA0003465798020000051
Figure BDA0003465798020000052
其中μ为地球引力常数，辅助变量s2和w分别为
Figure BDA0003465798020000053
通过式(11)(12)(15)将地影约束对容许控制集的改变转移到了动力学中。飞行时间和燃料消耗是小推力轨道机动中常考虑的两个指标，考虑地影约束后时间最优和燃料最优的性能指标为
Figure BDA0003465798020000054
其中J表示性能指标，t0为初始时刻，tf为终端时刻。根据最优控制理论，动力学系统的哈密顿函数为
Figure BDA0003465798020000055
其中与状态量x关联的协态变量λx为
λx＝[λp λf λg λh λk λL]T (21)
λm是与质量m关联的协态变量，Lg为拉格朗日函数
Figure BDA0003465798020000061
令xi＝x(i),λxi＝λx(i),i＝1,2,…6，则协态变量λx和λm的微分方程为
Figure BDA0003465798020000062
其中
Figure BDA0003465798020000063
为
Figure BDA0003465798020000064
由式(7)(10)(12)可知，up＝up(t,r(t))，即up是关于时间t和航天器位置矢量r的函数，则up是关于状态量x的函数，且是不连续的，所以无法直接应用庞特里亚金极小值原理，需要用连续函数去拟合up。由式(12)知，up＝up(ρp)，因此up的不连续性是阶跃函数sign(ρp)带来的。将通过步骤三设计平滑函数
Figure BDA0003465798020000065
用以拟合up使动力学系统连续，为了方便，仍将
Figure BDA0003465798020000066
记为up，并在该步骤中将up视作拟合后的连续函数。up关于状态量x的偏导数将在步骤三中给出。
根据庞特里亚金极小值原理，最优控制使得哈密顿函数式(20)最小化，因为0≤u≤1，0≤up≤1，所以最优推力方向α*为
Figure BDA0003465798020000067
最优推力幅值比u*为
Figure BDA0003465798020000071
其中，不考虑出现奇异弧段的情况，ρ为燃料最优问题下的开关函数
Figure BDA0003465798020000072
由式(26)得燃料最优问题下的最优推力幅值比是bang-bang的，通过扩展对数同伦使最优推力幅值比连续，将燃料最优问题的优化指标改写为
Figure BDA0003465798020000073
其中ε∈(01]为燃料最优同伦参数，当ε趋近于0的时候，式(28)趋近于最优的燃料指标。燃料最优问题的Lg和
Figure BDA0003465798020000074
分别改写为
Figure BDA0003465798020000075
Figure BDA0003465798020000076
i＝1,2,3,4,5,6,燃料最优 (30)
根据庞特里亚金极小值原理，因为0≤up≤1，所以燃料最优问题的最优推力幅值比u*变为
Figure BDA0003465798020000077
因此，对燃料最优问题的求解转换成了对ε取一系列值(ε＝1,0.1,0.001,…)时对应问题的求解，当ε的取值趋近于0时，其对应问题的解作为燃料最优解。
将终端约束和横截条件整体表述为
Figure BDA0003465798020000078
其中
Figure BDA0003465798020000079
为向量函数，包含终端约束和横截条件。终端约束由具体问题给出，横截条件则由最优控制原理相应给出。z为未知量，包括协态变量初值，由具体问题给出。
至此，考虑地影约束的小推力轨道最优控制模型已建立。
步骤三：在步骤一建立的基于角度差判断的地影约束模型基础上，设计平滑函数和平滑参数对地影约束进行平滑处理，在平滑过程中使地影区推力大小随平滑参数的变化而相应变化，使得以平滑参数作为同伦参数的同伦过程顺利进行。
步骤二已建立考虑地影约束的小推力轨道最优控制模型，其中平滑函数up＝up(ρp)及其关于状态量x的偏导数将在本步骤中给出。
在对不连续点进行平滑的过程中，使地影区推力幅值比同时变化，设计如下平滑函数
Figure BDA0003465798020000081
其中εp为平滑参数，
Figure BDA0003465798020000082
和
Figure BDA0003465798020000083
分别为平滑参数εp的上界和下界，即
Figure BDA0003465798020000084
为系数up在地影区外(ρp＞0)随ρp变化的取值上界，
Figure BDA0003465798020000085
为系数up在地影区内(ρp＜0)随ρp变化的取值下界，
Figure BDA0003465798020000086
为
Figure BDA0003465798020000087
随εp变化的取值上界。随着εp从
Figure BDA0003465798020000088
变化到
Figure BDA0003465798020000089
up在地影区内的取值下界从
Figure BDA00034657980200000810
变化到0，也即在对不连续点进行平滑处理的过程中，使地影区的推力幅值比u从
Figure BDA00034657980200000811
逐步变为0，最终满足地影约束。
给出up关于状态量x的偏导数，有
Figure BDA00034657980200000812
且
Figure BDA00034657980200000813
因此只需要给出
Figure BDA00034657980200000814
和
Figure BDA00034657980200000815
即可。根据式(33)得
Figure BDA00034657980200000816
其中sech表示双曲正割函数。根据式(7)得
Figure BDA0003465798020000091
其中
Figure BDA0003465798020000092
Figure BDA0003465798020000093
其中
Figure BDA0003465798020000094
为航天器位置矢量关于改进春分点根数的偏导数。将位置矢量r用改进春分点根数表示为
Figure BDA0003465798020000095
则偏导数
Figure BDA0003465798020000096
为
Figure BDA0003465798020000097
其中
Figure BDA0003465798020000098
Figure BDA0003465798020000099
Figure BDA0003465798020000101
Figure BDA0003465798020000102
Figure BDA0003465798020000103
Figure BDA0003465798020000104
步骤四：将步骤三中设计的平滑函数及其关于状态量的偏导数代入步骤二中建立的最优控制模型，并以平滑参数为同伦参数，结合同伦法和间接法求解考虑地影约束的小推力轨道最优控制问题，得到满足地影约束的小推力轨道最优控制和最优轨迹。
根据步骤三中给出的平滑函数和平滑参数，将平滑参数εp作为地影约束同伦参数，结合同伦法和间接法求解地影约束同伦参数εp从
Figure BDA0003465798020000105
到
Figure BDA0003465798020000106
取一系列值对应的最优控制问题，当
Figure BDA0003465798020000107
时，若
Figure BDA0003465798020000108
小于等于0.001，则对应的解为考虑地影约束的小推力轨道最优控制问题的解，即得到满足地影约束的小推力轨道最优控制和最优轨迹。
对于时间最优问题，根据步骤二中建立的考虑地影约束的小推力轨道最优控制模型，首先采用间接法求解不考虑地影约束的时间最优解，然后将其作为初值对地影约束同伦参数εp进行同伦，直到获得
Figure BDA0003465798020000109
的解，该解为考虑地影约束的小推力轨道时间最优问题的解，即得到满足地影约束的小推力轨道时间最优控制和时间最优轨迹。
对于燃料最优问题，根据步骤二中建立的考虑地影约束的小推力轨道最优控制模型，首先不考虑地影约束，采用间接法求解燃料最优同伦参数ε取值从1到一个趋近于0的值的同伦过程，当ε趋近于0时，对应的解为燃料最优解，然后将其作为初值对地影约束同伦参数εp进行同伦，直到获得
Figure BDA00034657980200001010
的解，该解为考虑地影约束的小推力轨道燃料最优问题的解，即得到满足地影约束的小推力轨道燃料最优控制和燃料最优轨迹。
步骤五：根据步骤四中得到的满足地影约束的小推力轨道最优控制和最优轨迹，对航天器进行相应制导控制，进而实现满足地影约束的小推力时间最优或燃料最优变轨。
有益效果：
1、本发明公开的考虑地影约束的小推力最优变轨方法，在圆锥形或圆柱形地影模型的基础上，建立基于角度差判断的地影约束模型，将对航天器是否位于地影区的判断转化为对角度差符号的判断，不受航天器相对于地球的方位限制，因此适用于全部轨道段，有利于对地影约束的平滑处理，为平滑函数的设计及同伦过程的顺利进行提供基础。
2、本发明公开的考虑地影约束的小推力最优变轨方法，通过设计平滑函数和平滑参数对地影约束进行平滑处理，在平滑过程中使地影区推力大小随平滑参数的变化而相应变化，使得以平滑参数为同伦参数的同伦过程顺利进行，实现对考虑地影约束的小推力最优变轨问题的求解，进而实现满足地影约束的小推力时间最优或燃料最优的变轨。
附图说明
图1为本发明中在圆锥形地影模型的基础上建立的基于角度差判断的地影约束模型的示意图；
图2为本发明中在圆柱形地影模型的基础上建立的基于角度差判断的地影约束模型的示意图；
图3为本发明中设计的地影约束平滑函数的示意图；
图4为实施例中本发明公开的考虑地影约束的小推力最优变轨方法流程图；
图5为实施例中地影约束下小推力时间最优变轨的轨迹示意图；
图6为实施例中地影约束下小推力时间最优变轨的推力序列、开关函数及地影约束开关函数的时间历程图，其中：图6(a)为地影约束下小推力时间最优变轨的推力序列的时间历程图、图6(b)为地影约束下小推力时间最优变轨的开关函数的时间历程图、图6(c)为地影约束下小推力时间最优变轨的地影约束开关函数的时间历程图；
图7为实施例中地影约束下小推力时间最优变轨的协态变量的时间历程及其局部放大图，其中：图7(a)为地影约束下小推力时间最优变轨的协态变量λp和λf的时间历程及其局部放大图、图7(b)为地影约束下小推力时间最优变轨的协态变量λg和λh的时间历程及其局部放大图、图7(c)为地影约束下小推力时间最优变轨的协态变量λk和λL的时间历程及其局部放大图；
图8为实施例中地影约束下小推力时间最优变轨的哈密顿函数的时间历程及其局部放大图；
图9为实施例中地影约束下小推力时间最优变轨的半长轴、偏心率及轨道倾角的时间历程图，其中：图9(a)为地影约束下小推力时间最优变轨的半长轴a的时间历程图、图9(b)为地影约束下小推力时间最优变轨的偏心率e的时间历程图、图9(c)为地影约束下小推力时间最优变轨的轨道倾角i的时间历程图；
图10为实施例中地影约束下小推力燃料最优变轨的轨迹示意图；
图11为实施例中地影约束下小推力燃料最优变轨的推力序列、开关函数及地影约束开关函数的时间历程图，其中：图11(a)为地影约束下小推力燃料最优变轨的推力序列的时间历程图、图11(b)为地影约束下小推力燃料最优变轨的开关函数的时间历程图、图11(c)为地影约束下小推力燃料最优变轨的地影约束开关函数的时间历程图；
图12为实施例中地影约束下小推力燃料最优变轨的协态变量的时间历程及其局部放大图，其中：图12(a)为地影约束下小推力燃料最优变轨的协态变量λp和λf的时间历程及其局部放大图、图12(b)为地影约束下小推力燃料最优变轨的协态变量λg和λh的时间历程及其局部放大图、图12(c)为地影约束下小推力燃料最优变轨的协态变量λk和λL的时间历程及其局部放大图；
图13为实施例中地影约束下小推力燃料最优变轨的半长轴、偏心率及轨道倾角的时间历程图，其中：图13(a)为地影约束下小推力燃料最优变轨的半长轴a的时间历程图、图13(b)为地影约束下小推力燃料最优变轨的偏心率e的时间历程图、图13(c)为地影约束下小推力燃料最优变轨的轨道倾角i的时间历程图。
具体实施方式
为了更好地说明本发明的目的和优点，下面结合附图对本发明的具体实施方式作进一步详细的说明。
本实施例公开的考虑地影约束的小推力最优变轨方法，为验证该方法，首先，选取航天器从GTO变轨至GEO的任务为主要研究对象，变轨的初始时刻为2023年3月20日21时37分12秒，轨道参数如表1所示，航天器的系统参数如表2所示，物理参数如表3所示。
表1 GTO和GEO的经典轨道根数
Figure BDA0003465798020000121
表2航天器的系统参数
Figure BDA0003465798020000131
表3物理参数
Figure BDA0003465798020000132
步骤一：在圆锥形或圆柱形地影模型的基础上，建立基于角度差判断的地影约束模型，将对航天器是否位于地影区的判断转化为对角度差符号的判断，不受航天器相对于地球的方位限制，因此适用于全部轨道段。
选取圆锥形地影模型，根据图1、式(7)(8)(11)(12)，即可建立基于角度差判断的地影约束模型。
步骤二：将地影约束添加在动力学中，通过最优控制原理建立考虑地影约束的小推力轨道最优控制模型。
根据式(19)(28)设定性能指标，包括时间最优和燃料最优，以式(15)(23)为状态变量和协态变量的微分方程；根据式(25)(26)(31)给出最优控制；在此实施例中，终端约束和横截条件为
Figure BDA0003465798020000133
Figure BDA0003465798020000134
其中xf为目标状态。为构建两点边值问题，给出初始猜测的未知量
Figure BDA0003465798020000141
根据打靶法，将对考虑地影约束的小推力轨道最优控制问题的求解描述为寻找未知量(式(50))使得动力学系统在终端满足终端约束和横截条件(式(48)(49))。
步骤三：在步骤一建立的基于角度差判断的地影约束模型基础上，设计平滑函数和平滑参数对地影约束进行平滑处理，在平滑过程中使地影区推力大小随平滑参数的变化而相应变化，使得以平滑参数作为同伦参数的同伦过程顺利进行。
根据式(33)，令
Figure BDA0003465798020000142
则系数up随ρp变化的曲线如图3所示，根据式(8)，ρp＞0表示航天器位于地影区外，ρp＜0表示航天器位于地影区内，当平滑参数εp＝1时，up随ρp变化的取值始终为1，即等价于不考虑地影约束；当平滑参数εp＝0.5时，up的取值在ρp符号变化点附近是光滑的，并且在地影区内取值的下界为0.5而不是0；当平滑参数εp＝0.001时，up的取值随ρp的变化接近于阶跃函数，但up的取值在ρp符号变化点附近仍然是光滑的。
步骤四：将步骤三中设计的平滑函数及其关于状态量的偏导数代入步骤二中建立的最优控制模型，并以平滑参数为同伦参数，结合同伦法和间接法求解考虑地影约束的小推力轨道最优控制问题，得到满足地影约束的小推力轨道最优控制和最优轨迹。
将步骤三中设计的平滑函数up＝up(ρp)代入步骤二的模型中，并根据式(34)(35)(36)(37)给出up关于状态量x的偏导数，代入步骤二的模型中，使得模型完整。
结合同伦法和间接法求解考虑地影约束的小推力轨道优化问题。将平滑参数εp作为地影约束同伦参数，从1同伦到0.001，每一个εp∈[0.0011]都对应一个最优控制问题，其中平滑参数εp＝0.001所对应的最优控制问题的解即为对应的考虑地影约束的小推力最优控制问题的解。对于时间最优问题，首先求解εp＝1所对应的最优控制问题的解，然后逐渐减小εp，将上一个问题的解作为下一个问题的初值，直到求解出εp＝0.001所对应的最优控制问题的解。对于燃料最优问题，令εp＝1，首先求解燃料最优同伦参数ε＝1所对应的最优控制问题的解，然后逐渐减小ε，将上一个问题的解作为下一个问题的初值，直到求解出ε＝0.0001所对应的最优控制问题的解，作为燃料最优解，接下来以该燃料最优解为初值，将εp从1同伦到0.001，流程与时间最优问题中的类似，直到求解出εp＝0.001所对应的最优控制问题的解。
为了方便，在本实施例中，令
Figure BDA0003465798020000143
表示满足地影约束的最优推力序列。
步骤五：根据步骤四中得到的满足地影约束的小推力轨道最优控制和最优轨迹，对航天器进行相应制导控制，进而实现满足地影约束的小推力时间最优或燃料最优变轨。
本发明公开的考虑地影约束的小推力最优变轨方法流程图如图4所示。该实施例的结果如下：
(1)满足地影约束的小推力时间最优变轨：
在小推力时间最优变轨中，如表2所示，航天器的初始质量设置为1000kg，最大推力幅值为0.8N。图5为地影约束下小推力时间最优变轨的轨迹示意图，其中黑色细线表示推力段而黑色粗细表示航天器位于地影区，推力器不可开机，航天器处于滑行段。在图6中，图6(a)为地影约束下小推力时间最优变轨的推力序列的时间历程图、图6(b)为地影约束下小推力时间最优变轨的开关函数的时间历程图、图6(c)为地影约束下小推力时间最优变轨的地影约束开关函数的时间历程图，因航天器频繁进出地影区，推力序列在前半段出现了在地影区关机的bang-bang控制形式，关机部分对应于ρp＜0的部分。根据最优控制理论，地影约束将使得协态变量和哈密顿函数出现不连续点，反映在图像中即为跳变现象。在图7中，图7(a)左侧为地影约束下小推力时间最优变轨的协态变量λp和λf的时间历程，右侧为对应的局部放大图；图7(b)左侧为地影约束下小推力时间最优变轨的协态变量λg和λh的时间历程，右侧为对应的局部放大图；图7(c)左侧为地影约束下小推力时间最优变轨的协态变量λk和λL的时间历程，右侧为对应的局部放大图。在图8中，上方为地影约束下小推力时间最优变轨的哈密顿函数的时间历程，下方为对应的局部放大图像。图7和图8中的局部放大图像清晰的显示了协态变量和哈密顿函数在航天器进出地影区时刻的跳变现象，即不连续点。此外，图8还显示出，当航天器不再进出地影区后，哈密顿函数变为了0，符合无地影约束下小推力时间最优变轨的最优控制条件。在图9中，图9(a)为地影约束下小推力时间最优变轨的半长轴a的时间历程图、图9(b)为地影约束下小推力时间最优变轨的偏心率e的时间历程图、图9(c)为地影约束下小推力时间最优变轨的轨道倾角i的时间历程图，表明航天器由初始的GTO轨道变轨至GEO轨道。在该实例中，航天器经过地影区共计31次，给动力学系统带来了62个不连续点，表明本发明在处理航天器多次进出地影区的小推力最优变轨问题上的优势。
(2)满足地影约束的小推力燃料最优变轨：
在燃料最优问题中，如表2所示，航天器的初始质量设置为125kg，最大推力幅值为0.8N。图10为地影约束下小推力燃料最优变轨的轨迹示意图，其中灰色细线表示最优控制给出的滑行段，黑色细线表示最优控制给出的推力段，黑色粗线表示航天器在地影区的滑行段。在图11中，图11(a)为地影约束下小推力燃料最优变轨的推力序列的时间历程图，为bang-bang控制形式，切换时刻包括最优控制给出的开关机时刻和航天器进出地影区时刻；图11(b)为地影约束下小推力燃料最优变轨的开关函数的时间历程图；图11(a)为地影约束下小推力燃料最优变轨的推力序列的时间历程图。在图12中，图12(a)左侧为地影约束下小推力燃料最优变轨的协态变量λp和λf的时间历程，右侧为对应的局部放大图；图12(b)左侧为地影约束下小推力燃料最优变轨的协态变量λg和λh的时间历程，右侧为对应的局部放大图；图12(c)左侧为地影约束下小推力燃料最优变轨的协态变量λk和λL的时间历程，右侧为对应的局部放大图，可以清晰的看到协态变量在航天器进出地影区时刻的跳变现象，即不连续点。在图13中，图13(a)为地影约束下小推力燃料最优变轨的半长轴a的时间历程图、图13(b)为地影约束下小推力燃料最优变轨的偏心率e的时间历程图、图13(c)为地影约束下小推力燃料最优变轨的轨道倾角i的时间历程图，表明航天器由初始的GTO轨道变轨至GEO轨道。在该实例中，燃料最优本身的bang-bang控制形式加上地影区关机约束，使得推力序列共计40次开关机切换，进一步表明本发明在处理航天器多次进出地影区的小推力最优变轨问题上的优势。
以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.考虑地影约束的小推力最优变轨方法，其特征在于：包括如下步骤，
步骤一：在圆锥形或圆柱形地影模型的基础上，建立基于角度差判断的地影约束模型，将对航天器是否位于地影区的判断转化为对角度差符号的判断，不受航天器相对于地球的方位限制，因此适用于全部轨道段；
步骤二：将地影约束添加在动力学中，通过最优控制原理建立考虑地影约束的小推力轨道最优控制模型；
步骤三：在步骤一建立的基于角度差判断的地影约束模型基础上，设计平滑函数和平滑参数对地影约束进行平滑处理，在平滑过程中使地影区推力大小随平滑参数的变化而相应变化，使得以平滑参数作为同伦参数的同伦过程顺利进行；
步骤四：将步骤三中设计的平滑函数及其关于状态量的偏导数代入步骤二中建立的最优控制模型，并以平滑参数为同伦参数，结合同伦法和间接法求解考虑地影约束的小推力轨道最优控制问题，得到满足地影约束的小推力轨道最优控制和最优轨迹。
2.如权利要求1所述的考虑地影约束的小推力最优变轨方法，其特征在于：步骤五，根据步骤四中得到的满足地影约束的小推力轨道最优控制和最优轨迹，对航天器进行相应制导控制，进而实现满足地影约束的小推力时间最优或燃料最优变轨。
3.如权利要求1或2所述的考虑地影约束的小推力最优变轨方法，其特征在于：针对圆锥形地影模型建立基于角度差判断的地影约束模型，具体实现方法如下：
将地球和太阳均视为球体，此时地影区域为圆锥形投影；S表示太阳中心，E表示地球中心，RS表示太阳半径，RE表示地球半径，SE所在直线为地影区域的中心线，V在中心线上表示本影区的顶点，αp为半影区圆锥半角，αu为本影区圆锥半角；航天器运行在轨道上，r为其在地心惯性坐标系下的位置矢量，r为其地心距，即位置矢量r的模，rS为太阳中心在地心惯性坐标系下的位置矢量；当航天器位于本影区或半影区时，均认为推力器无法工作，即不区分本影区和半影区；在某一时刻，航天器可能位于地影区外，也可能位于地影区内；以地心为圆心、r为半径作圆，该圆与地影区边界相交于以中心线为对称轴的两点，D点为与航天器位于同侧的交点；βsc为航天器位置矢量r与地影区中轴线间的夹角，满足
βsc∈[0 180°] (1)
βp为DE与地影区中心线间的夹角，有
βp＝αp+β (2)
其中
Figure FDA0003465798010000011
为辅助角，因为r＞RE，所以β∈(0 90°)；
令rS＝||rS||为日地距离，
Figure FDA0003465798010000021
为太阳在地心惯性坐标系下的方向向量，则
Figure FDA0003465798010000022
Figure FDA0003465798010000023
并将βsc和αp分别表示为
Figure FDA0003465798010000024
Figure FDA0003465798010000025
令ρp表示角量βsc和βp间的差值，用于表征航天器与地影区的位置关系，即
Figure FDA0003465798010000026
至此，在任意时刻，对航天器是否位于地影区的判断转化为对该时刻ρp的符号判断，即
Figure FDA0003465798010000027
式(7)(8)不受航天器相对于地球的方位限制，因此适用于全部轨道段；
针对圆柱形地影模型建立基于角度差判断的地影约束模型，具体实现方法如下：
将太阳光视为平行光，此时地影区域为圆柱形投影；E表示地球中心，RE表示地球半径，
Figure FDA0003465798010000028
为太阳在地心惯性坐标系下的方向向量，其所在直线为地影区域的中心线；航天器运行在轨道上，r为其在地心惯性坐标系下的位置矢量，r为其地心距；在某一时刻，航天器可能位于地影区外，也可能位于地影区内；以地心为圆心、r为半径作圆，该圆与地影区边界相交于以中心线为对称轴的两点，D点为与航天器位于同侧的交点；βsc为航天器位置矢量r与地影区中轴线间的夹角，满足式(1)(5)；βp为DE与地影区中心线间的夹角，有
Figure FDA0003465798010000029
则角量βsc和βp间的差值ρp为
Figure FDA00034657980100000210
至此，在任意时刻，对航天器是否位于地影区的判断转化为对该时刻ρp的符号判断，即式(8)；式(10)(8)不受航天器相对于地球的方位限制，因此适用于全部轨道段；
对航天器受到的地影约束描述为：当航天器位于地影区时，推力器无法工作，即地影约束使得原容许控制集U改变，记为Up，但等效为容许控制集不变，而在原推力幅值比u上添加系数up，得到考虑地影约束后的推力幅值比ue，即
Figure FDA0003465798010000031
其中系数up为
Figure FDA0003465798010000032
因此，将ρp称为地影约束开关函数；针对圆锥形地影模型，式(7)(8)(11)(12)即表述了基于角度差判断的地影约束模型；针对圆柱形地影模型，式(10)(8)(11)(12)即表述基于角度差判断的地影约束模型。
4.如权利要求3所述的考虑地影约束的小推力最优变轨方法，其特征在于：步骤二实现方法为，
采用改进春分点根数描述航天器的轨道运动，改进春分点根数与经典轨道根数间的关系为
Figure FDA0003465798010000033
其中a为轨道半长轴，e为轨道偏心率，i为轨道倾角，Ω为升交点赤经，ω为近地点幅角，ν为真近点角，p为半通径，L为真经度；令状态变量x为
x＝[p f g h k L]T (14)
m为航天器质量，再根据Gauss型摄动方程，航天器的动力学方程为
Figure FDA0003465798010000034
其中u∈[0 1]为推力幅值比，up为地影约束系数，Tmax为最大推力幅值，Isp为推力器比冲，g0为标准重力加速度，α为推力方向向量；Tmax和Isp在航天器飞行过程中保持不变；式(15)中的M、D分别为
Figure FDA0003465798010000041
Figure FDA0003465798010000042
其中μ为地球引力常数，辅助变量s2和w分别为
Figure FDA0003465798010000043
通过式(11)(12)(15)将地影约束对容许控制集的改变转移到了动力学中；飞行时间和燃料消耗是小推力轨道机动中常考虑的两个指标，考虑地影约束后时间最优和燃料最优的性能指标为
Figure FDA0003465798010000044
其中J表示性能指标，t0为初始时刻，tf为终端时刻；根据最优控制理论，动力学系统的哈密顿函数为
Figure FDA0003465798010000045
其中与状态量x关联的协态变量λx为
λx＝[λp λf λg λh λk λL]T (21)
λm是与质量m关联的协态变量，Lg为拉格朗日函数
Figure FDA0003465798010000051
令xi＝x(i),λxi＝λx(i),i＝1,2,…6，则协态变量λx和λm的微分方程为
Figure FDA0003465798010000052
其中
Figure FDA0003465798010000053
为
Figure FDA0003465798010000054
由式(7)(10)(12)知，up＝up(t,r(t))，即up是关于时间t和航天器位置矢量r的函数，则up是关于状态量x的函数，且是不连续的，所以无法直接应用庞特里亚金极小值原理，需要用连续函数去拟合up；由式(12)知，up＝up(ρp)，因此up的不连续性是阶跃函数sign(ρp)带来的；将通过步骤三设计平滑函数
Figure FDA0003465798010000055
用以拟合up使动力学系统连续，为了方便，仍将
Figure FDA0003465798010000056
记为up，并在该步骤中将up视作拟合后的连续函数；up关于状态量x的偏导数将在步骤三中给出；
根据庞特里亚金极小值原理，最优控制使得哈密顿函数式(20)最小化，因为0≤u≤1，0≤up≤1，所以最优推力方向α*为
Figure FDA0003465798010000057
最优推力幅值比u*为
Figure FDA0003465798010000058
其中，不考虑出现奇异弧段的情况，ρ为燃料最优问题下的开关函数
Figure FDA0003465798010000061
由式(26)得燃料最优问题下的最优推力幅值比是bang-bang的，通过扩展对数同伦使最优推力幅值比连续，将燃料最优问题的优化指标改写为
Figure FDA0003465798010000062
其中ε∈(01]为燃料最优同伦参数，当ε趋近于0的时候，式(28)趋近于最优的燃料指标；燃料最优问题的Lg和
Figure FDA0003465798010000063
分别改写为
Figure FDA0003465798010000064
Figure FDA0003465798010000065
根据庞特里亚金极小值原理，因为0≤up≤1，所以燃料最优问题的最优推力幅值比u*变为
Figure FDA0003465798010000066
因此，对燃料最优问题的求解转换成了对ε取一系列值(ε＝1,0.1,0.001,)时对应问题的求解，当ε的取值趋近于0时，其对应问题的解作为燃料最优解；
将终端约束和横截条件整体表述为
Figure FDA0003465798010000067
其中
Figure FDA0003465798010000068
为向量函数，包含终端约束和横截条件；终端约束由具体问题给出，横截条件则由最优控制原理相应给出；z为未知量，包括协态变量初值，由具体问题给出；
至此，考虑地影约束的小推力轨道最优控制模型已建立。
5.如权利要求4所述的考虑地影约束的小推力最优变轨方法，其特征在于：步骤三实现方法为，
步骤二已建立考虑地影约束的小推力轨道最优控制模型，其中平滑函数up＝up(ρp)及其关于状态量x的偏导数将在本步骤中给出；
在对不连续点进行平滑的过程中，使地影区推力幅值比同时变化，设计如下平滑函数
Figure FDA0003465798010000071
其中εp为平滑参数，
Figure FDA0003465798010000072
和
Figure FDA0003465798010000073
分别为平滑参数εp的上界和下界，即
Figure FDA0003465798010000074
Figure FDA0003465798010000075
为系数up在地影区外(ρp＞0)随ρp变化的取值上界，
Figure FDA0003465798010000076
为系数up在地影区内(ρp＜0)随ρp变化的取值下界，
Figure FDA0003465798010000077
为
Figure FDA0003465798010000078
随εp变化的取值上界；随着εp从
Figure FDA0003465798010000079
变化到
Figure FDA00034657980100000710
up在地影区内的取值下界从
Figure FDA00034657980100000711
变化到0，也即在对不连续点进行平滑处理的过程中，使地影区的推力幅值比u从
Figure FDA00034657980100000712
逐步变为0，最终满足地影约束；
给出up关于状态量x的偏导数，有
Figure FDA00034657980100000713
且
Figure FDA00034657980100000714
因此只需要给出
Figure FDA00034657980100000715
和
Figure FDA00034657980100000716
即可；根据式(33)得
Figure FDA00034657980100000717
其中sech表示双曲正割函数；根据式(7)得
Figure FDA00034657980100000718
其中
Figure FDA0003465798010000081
Figure FDA0003465798010000082
其中
Figure FDA0003465798010000083
为航天器位置矢量关于改进春分点根数的偏导数；将位置矢量r用改进春分点根数表示为
Figure FDA0003465798010000084
则偏导数
Figure FDA0003465798010000085
为
Figure FDA0003465798010000086
其中
Figure FDA0003465798010000087
Figure FDA0003465798010000088
Figure FDA0003465798010000089
Figure FDA0003465798010000091
Figure FDA0003465798010000092
Figure FDA0003465798010000093
6.如权利要求5所述的考虑地影约束的小推力最优变轨方法，其特征在于：步骤四实现方法为，
根据步骤三中给出的平滑函数和平滑参数，将平滑参数εp作为地影约束同伦参数，结合同伦法和间接法求解地影约束同伦参数εp从
Figure FDA0003465798010000094
到
Figure FDA0003465798010000095
取一系列值对应的最优控制问题，当
Figure FDA0003465798010000096
时，若
Figure FDA0003465798010000097
小于等于预设值，则对应的解为考虑地影约束的小推力轨道最优控制问题的解，即得到满足地影约束的小推力轨道最优控制和最优轨迹；
对于时间最优问题，根据步骤二中建立的考虑地影约束的小推力轨道最优控制模型，首先采用间接法求解不考虑地影约束的时间最优解，然后将其作为初值对地影约束同伦参数εp进行同伦，直到获得
Figure FDA0003465798010000098
的解，该解为考虑地影约束的小推力轨道时间最优问题的解，即得到满足地影约束的小推力轨道时间最优控制和时间最优轨迹；
对于燃料最优问题，根据步骤二中建立的考虑地影约束的小推力轨道最优控制模型，首先不考虑地影约束，采用间接法求解燃料最优同伦参数ε取值从1到一个趋近于0的值的同伦过程，当ε趋近于0时，对应的解为燃料最优解，然后将其作为初值对地影约束同伦参数εp进行同伦，直到获得
Figure FDA0003465798010000099
的解，该解为考虑地影约束的小推力轨道燃料最优问题的解，即得到满足地影约束的小推力轨道燃料最优控制和燃料最优轨迹。
CN202210029706.5A 2022-01-12 2022-01-12 考虑地影约束的小推力最优变轨方法 Pending CN114384803A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210029706.5A CN114384803A (zh) 2022-01-12 2022-01-12 考虑地影约束的小推力最优变轨方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210029706.5A CN114384803A (zh) 2022-01-12 2022-01-12 考虑地影约束的小推力最优变轨方法

Publications (1)

Publication Number Publication Date
CN114384803A true CN114384803A (zh) 2022-04-22

Family

ID=81201127

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210029706.5A Pending CN114384803A (zh) 2022-01-12 2022-01-12 考虑地影约束的小推力最优变轨方法

Country Status (1)

Country Link
CN (1) CN114384803A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115343960A (zh) * 2022-10-19 2022-11-15 北京理工大学 一种地月系统中航天器光照阴影规避控制方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106379555A (zh) * 2016-09-05 2017-02-08 北京理工大学 一种考虑j2摄动的近地卫星有限推力最优变轨方法
CN110736470A (zh) * 2019-11-06 2020-01-31 北京理工大学 一种不规则形状小天体附近连续推力轨道混合搜索方法
CN112084581A (zh) * 2020-09-24 2020-12-15 中国人民解放军国防科技大学 一种航天器小推力摄动交会轨迹优化方法及系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106379555A (zh) * 2016-09-05 2017-02-08 北京理工大学 一种考虑j2摄动的近地卫星有限推力最优变轨方法
CN110736470A (zh) * 2019-11-06 2020-01-31 北京理工大学 一种不规则形状小天体附近连续推力轨道混合搜索方法
CN112084581A (zh) * 2020-09-24 2020-12-15 中国人民解放军国防科技大学 一种航天器小推力摄动交会轨迹优化方法及系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
杨盛庆 等: "连续小推力条件下星座轨道机动方法研究", 中国空间科学技术, vol. 40, no. 04, pages 69 - 77 *
潘迅 等: "基于同伦方法三体问题小推力推进转移轨道设计", 深空探测学报, vol. 4, no. 03, pages 270 - 275 *
王帅 等: "小推力地球卫星圆轨道同轨调相设计方法研究", 宇航学报, vol. 34, no. 01, pages 1 - 8 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115343960A (zh) * 2022-10-19 2022-11-15 北京理工大学 一种地月系统中航天器光照阴影规避控制方法

Similar Documents

Publication Publication Date Title
Schaub et al. Impulsive feedback control to establish specific mean orbit elements of spacecraft formations
CN105631095B (zh) 一种等间隔发射的多约束地月转移轨道簇搜索方法
CN109625323A (zh) 一种卫星化学推进变轨方法及系统
US5124925A (en) Method for controlling east/west motion of a geostationary satellite
CN109539903A (zh) 一种固体运载火箭椭圆转移轨道迭代制导控制方法
CN105353621A (zh) 一种地球静止轨道卫星电推力器故障模式推力分配方法
CN105607478A (zh) 地球静止轨道航天器电推进转移轨道控制方法
CN109240340B (zh) 一种基于拟周期轨道的洛伦兹力多星编队构型方法
CN112016187B (zh) 一种基于混合动力的近地小行星交会任务轨道优化方法
CN114169072A (zh) 一种复杂约束下航天器连续推力轨道转移控制及优化方法
CN111994304B (zh) 一种静止轨道卫星小推力长期位置保持方法
CN114384803A (zh) 考虑地影约束的小推力最优变轨方法
CN115373423A (zh) 一种用于商业卫星的编队捕获方法
CN104252548A (zh) 一种燃料最优的火星探测器入轨目标点设计方法
Kos et al. Altair descent and ascent reference trajectory design and initial dispersion analyses
CN110723315B (zh) 一种天体表面弹道式飞行探测的轨迹生成方法
Yan et al. Optimal design of satellite formation relative motion orbits using least-squares methods
CN112009727B (zh) 平动点轨道最优小推力转移分段设计方法
Golan et al. Minimum fuel lunar trajectories for a low-thrust power-limited spacecraft
Melbourne Three-dimensional optimum thrust trajectories for power limited propulsion systems
CN111208847B (zh) 带太阳帆板的倾斜轨道卫星对日最优固定偏航角确定方法
Yang Earth-moon trajectory optimization using solar electric propulsion
JP2527895B2 (ja) 衛星管制方法
Gordienko et al. Analysis of trajectories of spacecraft launching into high circular orbits of the Moon artificial satellite
Boudad et al. Energy and Phasing Considerations for Low-Energy Transfers from Cislunar to Heliocentric Space

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