一种空间碎片的绳系拖曳最优离轨方法与流程

本发明涉及一种空间碎片的绳系拖曳最优离轨方法，尤其涉及航天器考虑姿态运动的转移轨道最优设计方法，属于航天器姿态和轨道动力学与控制领域。

背景技术：

在绳系卫星拖曳空间碎片离轨研究中，合理设计轨道控制策略，从而在平稳移除碎片的同时使燃料消耗尽可能少是碎片清除任务的研究重点。文献(wang,b.h.,meng,z.j.,huang,p.f."attitudecontroloftowedspacedebrisusingonlytether."actaastronautica138(2007):1-16)针对碎片离轨过程，仅研究了考虑系绳对拖船作用力的拖船轨道运动模型，并采用霍曼转移将空间碎片移除。文献(linskens,h.t.k.,andmooij,e."tetherdynamicsanalysisandguidanceandcontroldesignforactivespace-debrisremoval."journalofguidance,control,anddynamics(2016):1232-1243)设计了三个子系统组合完成碎片清除任务，分别是开环的降轨控制保证系统轨道运动正确，闭环的碎片拖船相对运动控制保证系统的构型正确，以及闭环的系统姿态控制保证系统的指向正确。该系统设计完备，但并未指出其在系统振荡，碎片旋转等干扰情况下推力控制律的可靠性和实用性。针对末时刻甩摆控制，文献(aslanov,v."swingprinciplefordeploymentofatether-assistedreturnmissionofare-entrycapsule."actaastronautica(2016):154-158)提出了控制系绳原长释放速率，从而在特定点剪断系绳，使碎片获取最小初速度降轨进入大气层的控制方法。

技术实现要素：

本发明公开的一种空间碎片的绳系拖曳最优离轨方法要解决的技术问题是：在保证绳系拖曳系统大范围最优轨道转移的同时，使绳系拖曳系统的姿态保持稳定，平稳地将空间碎片移除，同时减少碎片移除过程的燃料消耗。

本发明的目的通过以下技术方案实现。

本发明公开的一种空间碎片的绳系拖曳最优离轨方法，通过拉格朗日方法建立绳系拖曳系统的轨道平面内姿轨耦合动力学模型；采用轨道转移优化求解方法求解转移轨道最优燃料问题，在保证绳系拖曳系统大范围最优轨道转移的同时，使绳系拖曳系统的姿态保持稳定。在绳系拖曳末时刻对碎片实施甩摆释放控制，从而使碎片获取最大初速度，飞向更远轨道，以减少碎片移除过程的燃料消耗，完成碎片的移除任务。本发明能够在保证绳系拖曳系统大范围最优轨道转移的同时，使绳系拖曳系统的姿态保持稳定，平稳地将空间碎片进行移除，此外末时刻的甩摆控制也使碎片飞向更远轨道，减少碎片移除过程的燃料消耗。

本发明公开的一种空间碎片的绳系拖曳最优离轨方法，包括如下步骤：

步骤一：通过拉格朗日方法建立绳系拖曳系统的轨道平面内姿轨耦合动力学模型；

绳系拖曳系统由拖船航天器，空间碎片和系绳三部分组成，绳系拖曳系统简称系统，拖船航天器简称拖船，空间碎片简称碎片。在动力学建模中，将拖船和碎片均视为质点，系绳看作无质量的弹簧阻尼模型。定义惯性坐标系oexeye，绳系拖曳系统的轨道坐标系ooxoyo，其中，惯性系原点oe位于地球质心，xe和ye在轨道平面内相互垂直。xo由惯性系质心指向绳系拖曳系统的质心oo，yo在轨道平面内，指向绳系拖曳系统的速度方向。惯性系原点oe到oo的距离记为r，与xe的夹角记为θ(逆时针为正)，轨道角速度记为ω0，系绳的长度记为l，与yo轴的夹角记为α(顺时针为正)，拖船航天器的推力幅值为f，该推力与yo轴的夹角记为γ(顺时针为正)。拖船航天器和碎片的质量分别记为m1和m2，二者质量的总和记为m，地球引力常数记为μ，k＝es/l0是绳的弹性系数，其中e表示杨氏模量，s表示系绳的横截面积，l0表示系绳原长；c＝ds/l0是黏性阻尼系数，d表示黏性系数。选取系统的广义坐标r，θ，l和α。采用拉格朗日方法建立系统的动力学模型如公式(1)-(4)所示：

其中，h(·)表示heaviside阶跃函数。由于广义坐标中，轨道参数r远远大于其他各参数的数值，为避免数值运算中数字的相对舍入导致计算错误，对轨道参数r进行归一化处理，令

其中，r0为系统初始时刻的轨道半径。将式(5)带入式(1)-(4)中，得到系统的动力学方程如式(6)-(9)所示。

如式(6)-(9)所示的系统的动力学方程即为建立的绳系拖曳系统的轨道平面内姿轨耦合动力学模型。

步骤二：采用轨道转移优化方法求解转移轨道最优燃料问题，寻找最优推力控制方法，从而获取绳系拖曳系统最优转移方案，在保证绳系拖曳系统大范围最优轨道转移的同时，使绳系拖曳系统的姿态保持稳定。

作为优选，采用轨道转移优化方法求解转移轨道最优燃料问题即采用直接法将转移轨道分成若干段，采用内点法根据节点的数值及其一阶梯度信息寻找绳系拖曳系统最优推力控制方法，从而获取最优转移方案。在保证绳系拖曳系统大范围最优轨道转移的同时，使绳系拖曳系统的姿态保持稳定。

作为进一步优选，步骤二具体实现方法为：

在拖船航天器配备有小推力发动机，推力的幅值保持不变、推力方向角在轨道平面内变化范围为[-π/2π/2]的情况下，建立绳系拖曳系统的轨道转移问题。该情况下，燃料的消耗正比于推力作用的时间，即拖船推力作用时间越短，所消耗的燃料越少，燃料最优问题转化为时间最优问题。所述时间最优问题的目标函数为：

j＝tf(10)

其中，tf表示末端时刻。动力学约束函数为式(6)-(9)，其中γ为控制变量。初始边界条件如下：

末端时刻边界条件如下：

变量最值约束如下：

其中，下标0和f分别表示各个状态量在初始时刻和末时刻的取值，下标max和min则分别表示对应变量的最大取值和最小取值。时间最优问题描述为在动力学方程约束(6)-(9)，边界条件约束(11)和(12)和取值范围约束(13)下，寻找γ(t)，使得目标函数j取值最小。采用直接法求解该小推力转移时间最优问题。首先，将推力作用于轨道切线方向，通过数值积分获得系统状态量，同切向推力和末端时间一起作为最优问题的猜测初值。将连续的转移轨道状态变化量离散成若干个采样点上的状态列向量，采用hermit-simpson差值方法对各个区段中点的状态量及其一阶导数进行插值计算，控制变量采用牛顿插值方法，从而建立区间两端点和中点的动力学约束，保证轨道数值的准确性。最后采用内点法根据节点的数值及其一阶梯度信息寻找绳系拖曳系统最优推力控制方法，从而获取最优系统转移方案，在保证绳系拖曳系统大范围最优轨道转移的同时，使绳系拖曳系统的姿态保持稳定。

步骤三：在绳系拖曳碎片的末时刻对碎片实施甩摆释放控制，从而使碎片获取最大初速度，飞向更远轨道，以减少碎片移除过程的燃料消耗，完成碎片的移除任务。

绳系拖曳系统在到达目标轨道附近后，通过控制系绳原长变化率，对碎片进行末时刻甩摆控制。此时由于绳系拖曳系统已到达目标轨道附近，系统不存在大范围的轨道运动，因此仅考虑系统的姿态变化过程，动力学模型如公式(14)所示：

对系绳原长采用如公式(15)所示的控制方法：

其中，λ为比例系数。即根据系统的摆动角α对系绳的释放速度进行控制，使得系统的摆动角α不断增大。在不考虑系统轨道运动的情况下，碎片在单摆的最低点时具有最大速度，在单摆的最高点时速度为零。因此在碎片位于单摆最低点且速度方向与轨道切向速度方向一致时，将系绳剪断，绳系切断条件如式(16)所示。

此时，碎片将单摆的初速度和轨道速度叠加，从而获取最大初速度，飞向更远轨道，完成碎片的移除任务。

有益效果：

1、本发明公开的一种空间碎片的绳系拖曳最优离轨方法，建立轨道平面内绳系拖曳系统的姿轨耦合动力学模型，其中，将碎片和拖船视为质点，系绳看作无质量的弹簧阻尼，同时将系统的轨道运动纳入动力学模型中。通过对轨道参数进行归一化处理，使仿真计算的误差更小，通过考虑轨道和姿态的耦合作用使该模型更加精确，能真实反映绳系拖曳系统在实际应用中的运动状态。

2、本发明公开的一种空间碎片的绳系拖曳最优离轨方法，采用直接法，将转移轨道分成若干区间，然后采用内点法根据节点的数值及其一阶梯度信息寻找最优推力控制方法，从而获取最优转移方案，使系统完成大范围的最优轨道转移。

3、本发明公开的一种空间碎片的绳系拖曳最优离轨方法，在碎片移除的末时刻，采用控制系绳原长释放速率的方式不断增大系统的摆动角，然后在碎片具有最大摆动速度且方向与轨道运动方向一致时将系绳剪断，从而使碎片获取最大初速度飞向更远轨道，完成碎片的移除任务。

附图说明

图1为本发明中绳系拖曳系统的示意图；

图2为实施例的绳系拖曳系统的转移轨道；

图3为实施例的绳系拖曳系统的状态量变化情况；

图4为实施例的最优推力方向角；

图5为实施例的绳系拖曳系统的转移轨道；

图6为实施例的绳系拖曳系统的轨道变量；

图7为实施例的绳系拖曳系统的姿态变量；

图8为实施例的绳系拖曳系统系绳原长(算例a)；

图9为实施例的绳系拖曳系统姿态角(算例a)；

图10为实施例的碎片轨道的近地点高度(算例a)；

图11为实施例的绳系拖曳系统系绳原长(算例b)；

图12为实施例的绳系拖曳系统姿态角(算例b)；

图13为实施例的碎片轨道的近地点高度(算例b)。

具体实施方式

为了更好地说明本发明的目的和优点，下面结合附图对本发明的具体实施方式作进一步详细的说明。

本实施例公开的一种空间碎片的绳系拖曳最优离轨方法，为验证该方法，选取运行在半径为42164km的圆轨道上的绳系拖曳系统作为主要研究对象。其系统参数如下表所示。

表1绳系拖曳系统的结构参数.

步骤一：通过拉格朗日方法建立绳系拖曳系统的轨道平面内姿轨耦合动力学模型。

拉格朗日方法如公式(17)所示：

其中，t是系统动能，q表示系统的广义坐标，在本实施例中，vg表示系统的重力势能，ve代表弹性势能。广义坐标q对应的广义外力q的具体形式如式(18)所示：

其中，fd是广义弹性力，f是拖船推力矢量在惯性系下的分量，r1是拖船在惯性系中的位置矢量。系统的动能具体形式如式(19)所示：

其中，mi和vi分别代表第i个体的质量和绝对线速度，i＝1表示拖船，i＝2表示碎片。重力势能形式如式(20)所示：

系绳的弹性势能表示为式(21)：

其中，h(·)表示heaviside阶跃函数，当l≤l0时，该函数取值为零。系绳的广义阻尼力用式(22)表示：

将以上力与能量的表达式带入式(17)中，并对轨道半径r进行归一化处理，然后得到一组表示绳系拖曳系统轨道平面内姿轨耦合动力学模型，如式(6)-(9)所示。

步骤二：采用轨道转移优化方法求解转移轨道最优燃料问题，采用直接法将转移轨道分成若干段，采用内点法根据节点的数值及其一阶梯度信息寻找绳系拖曳系统最优推力控制方法，从而获取最优轨道转移方案，在保证绳系拖曳系统大范围最优轨道转移的同时，使绳系拖曳系统的姿态保持稳定。

考虑拖船航天器具有幅值为1n的连续推力，沿轨道系yo轴方向作用于整个绳系拖曳系统。广义坐标的仿真初值取为：

v＝0,θ＝0,f＝7.29×10^-5,l＝1000m,u＝0,α＝0,w＝0(23)

仿真时长为10⁶s。系统质心在拖船推力的作用下的运动轨迹如图2所示。绳系拖曳系统由geo轨道出发(图中蓝色空心圆)，在拖船推力持续作用下，系统轨道呈螺旋式上升。在整个仿真过程中，系统质心高度上升75436km，实现升高碎片轨道的目标。在升轨过程中，系绳长度和系统姿态变化情况如图3所示。由于拖船航天器受轨道推力的作用产生加速度，从而导致系绳发生伸长变形，在系绳阻尼的作用下，其变形量的幅值逐渐减小，最终稳定在某一数值附近，保证系统的稳定。绳系拖曳系统在重力梯度力矩的作用下产生微小的姿态角运动，姿态角不为零时，拖船推力对系统产生恢复力矩，使姿态角α趋于零值，表现为振荡稳定过程。整体上，系统姿态角α数值较小，系统整体上实现稳定离轨。

考虑拖船航天器配备有小推力发动机，推力的幅值保持不变，推力方向在轨道平面内变化范围为[-π/2π/2]，将碎片由geo轨道转移至坟墓轨道(高于geo轨道300km的圆轨道)的轨道转移最优问题。值得注意的是，拖船的推力不仅是系统进行轨道转移的能量来源，也是系绳保持时刻张紧的必要条件。由于连接拖船和碎片的系绳仅能产生拉力对碎片进行控制，为保证系统的稳定，避免碎片产生无规律的随机运动对系统产生扰动，拖船推力需时刻作用于系统，以保证系绳张紧，从而时刻产生拉力对碎片进行姿态和轨道的控制。考虑推力幅值不变，方向角在一定范围内可变的燃料最优问题。该情况下，燃料的消耗正比于推力作用的时间，即拖船推力作用时间越短，所消耗的燃料越少，燃料最优可等效为时间最优。该时间最优问题的目标函数为：

j＝tf(24)

其中，tf表示末端时刻。动力学约束函数为(6)-(9)，其中γ为控制变量。初始边界条件如下：

末端时刻边界条件如下：

变量最值约束如下：

时间最优问题描述为在动力学方程约束(6)-(9)，边界条件约束(25)-(26)和取值范围约束(27)下，寻找γ(t)，使得式(24)取值最小。求解该小推力转移轨道最优问题的方法主要分为间接法和直接法。间接法将小推力轨道优化问题描述为最优控制问题，然后基于变分法推导最优解存在的一阶必要条件，从而将最优控制问题转化为两点边值问题，求解两点边值问题获得原问题的解。在求解过程中，协态变量初值未知，末值满足特定的非线性方程组，即打靶方程，对于协态变量初值的猜测较为困难和复杂。直接法则是将小推力轨道优化问题描述为一个最优控制问题，将连续无穷维的最优控制问题离散成有限维的非线性规划(nlp)问题，将离散点上的状态量和控制量等参数作为nlp问题的优化向量，将动力学约束方程通过一定插值方法转换为代数约束，同初末状态约束，边界约束等约束一起构成nlp问题的约束列向量，通过求解该非线性规划问题获得原问题的解。具体地，在本问题中将推力作用于轨道切线方向，通过数值积分获得系统任意时刻的状态量猜测值，同切向推力一起作为最优问题的猜测初值。将连续的转移轨道状态变化量离散成200个采样点上的状态列向量，采用hermit-simpson差值方法对各个区段中点的状态量及其一阶导数进行插值计算，控制变量中间点的数值采用牛顿插值方法，从而建立区间两端点和中点的动力学约束，保证轨道数值的准确性。最后调用matlab中的fmincon函数对该最优问题进行求解。采用的具体算法是内点法，它是一种求解非线性凸优化问题的算法，通过引进罚函数，其作用相当于在可行域的边界上设置障碍，使求解的迭代过程始终在可行域内部进行。由于这种罚函数使得迭代点保持在可行域内部，故称为内部罚函数或障碍函数。

考虑绳系拖曳系统系绳不可伸长的情况，将系统的姿态运动和轨道运动均纳入最优转移问题的计算中，采用动力学模型(式(28))进行nlp问题的计算。

获得的最小转移时间下的拖船推力方向角曲线如图4所示。整体上，推力方向角从47.98度出发，逐渐增大到90度(约束上限)，在该过程中，轨道高度不断上升，但同时轨道径向速度v也不断增大。由于末端状态为圆轨道，v取值为零，因此在一定时间后推力作用于相反方向，以消除速度v，从而满足末端状态约束条件。在15000s附近出现推力方向角的多次突变。事实上，当推力方向角γ与系统摆动角(又名姿态角)α的差值不为零时，即推力产生力矩作用于绳系拖曳系统时，会导致系统产生绕质心的姿态运动，同时由于geo轨道上重力梯度力矩十分微弱，使得由推力产生的力矩作为主要作用力矩为姿态角α提供角加速度，使系统发生旋转情况。但由于仿真约束中要求α角的变化区间为[-π/2π/2]，这使得推力方向角不断发生变化，以消除姿态角的偏移，保证其在区间约束范围内。将该推力方向角变化情况带入完整的系统动力学模型(6)-(9)中，对系统轨道转移过程进行数值积分，各状态量的变化情况如图5-7所示。图5中，绳系拖曳系统在拖船推力的作用下采用椭圆曲线进行轨道转移，其中初始猜测轨道是切向推力情况下系统的运动轨迹，最优转移轨道即求解所得的系统最小时间转移轨迹，最终转移轨道是采用推力方向角最优解重新对系统动力学方程积分所得的系统真实轨迹，由图可知，初始猜测轨道，最优转移轨道和最终转移轨道之间的差值较小。图6和7反映采用最优解积分所得的系统各参数变化情况。由于最优解中拖船推力方向角在一段时间后频繁发生突变，这使得在最终转移轨道中系绳长度，系统的姿态角均发生明显振荡，为系统带来不稳定的影响。总体来说，系统的姿态角和相对距离均在一定范围内变化，无明显单方向增大的趋势，这表明系统在转移过程中较为稳定。本实施例中最优轨道转移问题的转移精度及转移时间如下表所示：

表2绳系拖曳系统的状态变量(算例6b)

表格最左侧一列表示系统的各个状态量和末端时刻，第二列表示状态量在末端时刻的目标值(最优问题的末端边界条件)，第三、四和五列分别表示初始猜测的末端时刻值，最优解中的系统末端时刻值和使用最优解积分所得的系统轨道末端时刻值与第二列的差值，没有末端目标值的量则取其具体数值代替差值表示在表格中。从表格得知，采用切向解获得的初始猜测轨道在末端时刻具有较大的径向速度，不满足geo轨道向坟墓轨道(圆轨道)转移的目标，实际满足条件的转移则需在一定时间之后对不为零的径向速度进行消除，因此转移时间将大于切向推力的情况。最优解中的轨道半径r，径向速度v和轨道角速度f和目标值的差距在10^-9以内，与fmincon函数求解过程中容许误差的设置值有关，在一定的误差容许度下，认为计算所得的最优解满足末端边界条件。采用推力方向角最优解对完整的系统动力学模型进行数值积分，获得的系统各状态量末端时刻取值与目标值的差距比最优解的差距大，主要是由于完整模型考虑系统的姿态运动等因素，推力方向角发生反复突变时，系绳的变形和系统的姿态角均产生很大的扰动，进而影响系统沿最优轨道运动的精度。总体来说，轨道转移的相对误差在1.92％，精度仍然较高。

步骤三：在绳系拖曳的末时刻对碎片实施甩摆释放控制，从而使碎片获取最大初速度，飞向更远轨道，以减少碎片移除过程的燃料消耗，完成碎片的移除任务。

绳系拖曳系统在到达目标轨道附近后，采用式(14)所示的姿态动力学模型，采用式(15)的控制方法对碎片实施末端时刻甩摆控制。在不考虑系统轨道运动的情况下，碎片绕拖船航天器做单摆运动，摆动角为α，在单摆的最低点时碎片具有最大速度，因此可在碎片位于单摆最低点且速度方向与轨道切向速度方向一致，即满足式(16)时，将系绳剪断。此时碎片将单摆的初速度和轨道速度叠加，从而具有最大的初始速度。数值仿真初值设定如下表所示：

表3绳系拖曳系统的初始状态参数

其中，算例a考虑系绳初始张紧，碎片初始时刻位于拖船下方，姿态角为80度，角速度为零的情况，此时，单摆在重力作用下沿轨道系x轴进行最大摆动角为10度的往复摆动；算例b考虑系绳初始张紧，碎片初始时刻位于单摆运动的中间某点，单摆的摆动幅度大于目前的姿态角。算例a的仿真结果如图8-10所示。在式(15)的控制律作用下，系绳原长，绳系拖曳系统的姿态角变化幅度不断增大，在曲线中星号点位置释放碎片，此时碎片处于拖船的正下方，且与轨道速度方向一致，碎片可获得最大初始速度，抬升其轨道高度。在星号之后的下一个摆动周期里，绳系拖曳系统的摆动角将会超出[0180]deg范围，此时可能出现系绳松弛，碎片进入无规则运动状态，系统将不再稳定。因此系绳的剪断操作应当在摆动角不超出[0180]deg范围的前一个摆动周期进行。算例b的剪断时刻为图11-13中星号位置，即碎片位于单摆的最低点，且速度方向与轨道方向一致的时刻，剪断系绳，碎片可获得最大初始速度。

系绳控制律中λ取值不同时，系统的摆动角增速不同，表4表示同一绳系拖曳系统在λ取不同数值的控制律作用下，系绳长度变化范围，碎片在最大位置释放所需时间和系绳拉力等参数的取值情况。系绳原长为1000m，初始时刻姿态角α为80deg，取值为零。

表4绳系拖曳系统的系统参数.

由表中数据可知，λ越大，系绳的释放速度越大，系统的摆动角增长越快，系绳的拉力也越大，到达碎片释放条件的时间越短。但由于λ数值较大，导致系绳原长的变化速度和范围也相应增大，说明若需要尽快剪断系绳释放碎片，要求系绳原长余量较大。此外，系绳释放速度越大，对释放装置的要求也越高。

对不同系绳长度情况下，甩摆控制过程中各个参数进行对比，结果如下表所示：

表5绳系拖曳系统的系统参数.

由表中数据可知，系绳的原长越长，碎片释放时的轨道高度增量越大，这说明在碎片释放前，应尽可能将系绳完全展开，从而获得更大的轨道高度增量。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。