CN104793613A - 一种航天器在日地系统不稳定平动点轨道间交会控制方法 - Google Patents

Abstract

本发明提供一种航天器在日地系统不稳定平动点轨道间交会控制方法包括：确定航天器在日地系统平动点轨道交会的初始轨道、目标轨道以及交会过程时间，建立受控航天器交会动力学模型；基于建立的受控航天器交会动力学模型，建立非线性最优控制问题；采用保辛数值法在有限长时间内求解非线性最优控制问题，得到协态变量；利用获得的协态变量更新当前子时间区间的控制输入，完成一个子时间区间的动力学仿真，利用导航方法获得当前航天器状态变量；进行子时间区间的递进，以当前航天器终端状态为下一个子时间区间的初始状态，依次重复步骤302至400，直至完成航天器在不稳定平动点轨道间交会任务。本发明能保证高精度的航天器交会的客观性和性能。

Description

一种航天器在日地系统不稳定平动点轨道间交会控制方法

技术领域

本发明涉及航天器深空探测技术领域，尤其涉及一种航天器在日地系统不稳定平动点轨道间交会控制方法。

背景技术

深空探测特别是平动点附近的空间探测是近年来航天方向的重要研究任务。图1为日地系统圆形限制性三体模型示意图。参照图1，由于圆形限制性三体问题的三个共线不稳定平动点重力与离心力平衡，这些空间位置是放置诸如空间站、空间望远镜等空间设备的理想场所。图2为航天器在日地系统L2平动点附近南北Halo间交会图。参照图2，对于空间站的一个典型且重要的问题就是在平动点附近Halo轨道族之间，进行航天器自动且高精度的交会控制任务。近地轨道附近的航天器交会任务已经有了很多成功的实际案例和成果，但是日地系统不稳定平动点Halo轨道族间还有进行过真正的航天器交会任务，也还没有进行过深入的研究。根据圆形限制性三体问题模型，平动点轨道相比于开普勒轨道有着不同的轨道周期，且相对速度也不同。在开普勒轨道上所做的工作不能应用于平动点轨道。为了实现平动点附近航天器的自动且高精度的交会，在航天器交会的控制问题中必须考虑圆形限制性三体问题轨道的非线性特性。进一步由于共线平动点周期轨道具有内在不稳定性，如果航天器相对与精确轨道的微小偏移也会导致一段时间后航天器相对于理想轨道的较大偏差。另外，来自天体的引力、导航误差、测量误差，动力装置的饱和界限以及推力的非线性特性等等都会不可避免的出现在航天器在平动点轨道的交会过程中。因此如果在上述情况下使用开环控制的话，自主及高精度的航天器交会的客观性和性能就会有所降低。另一方面，闭环反馈控制可以使用反馈信息抵消外部干扰和模型误差的影响。很明显，在上述情况下闭环反馈控制相比于开环控制有很大的优势。因此，亟需建立一种考虑控制系统非线性的闭环反馈控制策略和求解方法，为日地系统不稳定平动点附近Halo轨道族之间航天器高精度自主交会任务奠定基础。

发明内容

本发明主要解决航天器在日地系统不稳定平动点轨道族间交会的开环控制易受外部干扰和模型误差影响的问题，提出一种航天器在日地系统不稳定平动点轨道间交会控制方法，采用在线求解开环非线性最优控制，并不断更新控制输入构成闭环控制，可以使用反馈信息抵消外部干扰和模型误差的影响，保证高精度的航天器交会的客观性和性能。

本发明提供了一种航天器在日地系统不稳定平动点轨道间交会控制方法，包括以下步骤：

步骤100，确定航天器在日地系统平动点轨道交会的初始轨道、目标轨道以及交会过程时间，并建立受控航天器交会动力学模型；

步骤200，综合考虑航天器在交会过程中的能量消耗以及交会精度，并基于建立的受控航天器交会动力学模型，建立非线性最优控制问题；

步骤300，采用保辛数值方法在有限长时间内求解非线性最优控制问题，得到协态变量，包括以下步骤：

步骤301，将交会过程时间[t,t+T]等分为N个子时间区间，每一个子时间区间的时间步长为η＝T/N，其中，每一个子时间区间表示为t₀＝t,t₁＝t+η,...,t_j＝t+jη,...,t_N＝t_f；

步骤302，在第j个子时间区间中，获得开环最优控制系统的协态变量，其中，j从0开始到N结束；

步骤400，利用获得的协态变量更新当前子时间区间的控制输入，将得到控制输入作用到航天器交会动力学模型，完成一个子时间区间的动力学仿真，并通过导航方法获得当前航天器状态变量；

步骤500，进行子时间区间的递进，以当前航天器终端状态为下一个子时间区间的初始状态，依次重复步骤302至步骤400，直至完成航天器在不稳定平动点轨道间交会任务。

进一步的，建立受控航天器交会动力学模型，包括：

建立在日地系统两天体质心上的参考坐标系，参考坐标系与两天体以相同的轨道速度旋转，参考坐标系的x轴从系统质心处向地球延伸，z轴沿着系统的角动量方向延伸，y轴根据右手法则确定；

受控航天器交会动力学模型的动力学方程表示为：

$\left\{\begin{array}{c}\stackrel{\·\·}{x}-2\stackrel{\·}{y}-x=-\frac{(1-\mathrm{\μ})(x+1+1/\mathrm{\γ})}{{\mathrm{\γ}}^3{d}_1^3}-\frac{\mathrm{\μ}(x+1)}{{\mathrm{\γ}}^3{d}_2^3}+\frac{1-\mathrm{\μ}+\mathrm{\γ}}{\mathrm{\γ}}+u_x\\ \stackrel{\·\·}{y}+2\stackrel{\·}{x}-y=-\frac{(1-\mathrm{\μ})y}{{\mathrm{\γ}}^3{d}_1^3}-\frac{\mathrm{\μy}}{{\mathrm{\γ}}^3{d}_2^3}+u_y\\ \stackrel{\·\·}{z}=-\frac{(1-\mathrm{\μ})z}{{\mathrm{\γ}}^3{d}_1^3}-\frac{\mathrm{\μz}}{{\mathrm{\γ}}^3{d}_2^3}+u_z\end{array}\right.---\left(1\right)$

其中， $d_1=\sqrt{{(x+1+1/\mathrm{\γ})}^2+y^2+z^2},d_2=\sqrt{{(x+1)}^2+y^2+z^2},$ 方程中的点号代表旋转框架内的时间微分，μ是地球的质量与地球太阳质量之和的比值，γ表示地球和平动点L₂之间的距离，符号u_x、u_y和u_z分别表示受控航天器的控制输入，表示坐标x对时间求导。

进一步的，步骤200中建立的非线性最优控制问题的性能指标表示为：

$\begin{array}{c}J=\frac{1}{2}{[M_{f}x(t+T)-\mathrm{\ψ}]}^T{S}_f[M_{f}x(t+T)-\mathrm{\ψ}]\\ +\frac{1}{2}{\∫}_t^{t+T}[{(x-x_d)}^{T}Q(x-x_d)+{(u-u_d)}^{T}G(u-u_d)]\mathrm{d\τ}\end{array}---\left(4\right)$

其中，x∈R^n×1表示状态变量，u∈R^m×1表示控制变量，x_d∈R^n×1表示目标状态向量，u_d∈R^m×1表示目标控制输入，t表示时间，τ∈[t,t+T]表示预测未来系统状态的时间变量，ψ∈R^m×1表示线性混合终端状态M_fx(t+T)的目标值，Q∈R_n×n表示半正定矩阵，G∈R^m×m表示正定矩阵，S_f表示半正定终端权矩阵，M_f表示给定的矩阵。

进一步的，获得开环最优控制系统的协态变量，包括：

在第j个子区间中，用插入了r个等间距点的r-1阶Lagrange多项式来代替状态变量x(τ)，而用插入了s个等间距点的s-1阶Lagrange多项式来代替协态变量λ(τ)，状态变量x(τ)和协态变量λ(τ)表示为：

$x\left(\mathrm{\τ}\right)=(M\⊗I){\stackrel{\‾}{x}}_j---\left(11\right)$

$\mathrm{\λ}\left(\mathrm{\τ}\right)=N_1{\mathrm{\λ}}_{j-1}+(\underset{\‾}{N}\⊗I){\stackrel{\‾}{\mathrm{\λ}}}_j+N_n{\mathrm{\λ}}_j---\left(12\right)$

其中，符号代表Kronecker积，I代表n×n阶单位矩阵，向量λ_j-1与λ_j分别代表第j个子时间区间左端与右端的协态变量，

${\stackrel{\‾}{x}}_j={\{{\left({\stackrel{\‾}{x}}_j^1\right)}^T,{\left({\stackrel{\‾}{x}}_j^2\right)}^T,...,{\left({\stackrel{\‾}{x}}_j^r\right)}^T\}}^T,{\stackrel{\‾}{\mathrm{\λ}}}_j={\{{\left({\stackrel{\‾}{\mathrm{\λ}}}_j^2\right)}^T,{\left({\stackrel{\‾}{\mathrm{\λ}}}_j^3\right)}^T,...,{\left({\stackrel{\‾}{\mathrm{\λ}}}_j^{s-1}\right)}^T\}}^T;$

引入对偶变量变分原理如下：

$\mathrm{\δ}\stackrel{\‾}{S}={\∫}_t^{t+T}{\left(\mathrm{\δx}\right)}^T(-\stackrel{\·}{\mathrm{\λ}}-\frac{\∂H}{\∂x})\mathrm{d\τ}+{\∫}_t^{t+T}{\left(\mathrm{\δ\λ}\right)}^T(\stackrel{\·}{x}-\frac{\∂H}{\∂\mathrm{\λ}})\mathrm{d\τ}+{\mathrm{\λ}}^T\mathrm{\δx}{|}_t^{t+T}=0---\left(13\right)$

将离散的状态变量(11)和协态变量(12)代入上述对偶变量变分原理，获得如下线性方程

AZ＝B  (14)

其中，Z＝{λ(t₀)^T,λ(t₁)^T,…,λ(t_j)^T,…,λ(t_N)^T}^T；

$B={\{{(x_0-{\mathrm{\ζ}}_1^1)}^T,-{({\mathrm{\ζ}}_2^1+{\mathrm{\ζ}}_1^2)}^T,...,-{({\mathrm{\ζ}}_2^j+{\mathrm{\ζ}}_1^{j+1})}^T,...,-{({\mathrm{\ζ}}_2^{N-1}+{\mathrm{\ζ}}_1^N)}^T,-M_f^T{S}_f\mathrm{\ψ}-M_f^T{S}_f{M}_f{\mathrm{\ζ}}_2^N\}}^T;$

利用公式(14)得到开环最优控制系统的协态变量。

进一步的，通过以下公式利用获得的协态变量更新当前子时间区间的控制输入：

u＝u_d-G^-1(B(τ))^Tλ  (15)

其中，u∈R^m×1表示控制变量，u_d∈R^m×1表示目标控制输入，G∈R^m×m表示能量加权矩阵，是正定矩阵，τ表示时间，B(τ)表示受控非线性动力学方程对控制变量的偏导数矩阵，λ表示协态变量。

本发明提供的一种航天器在日地系统不稳定平动点轨道间交会控制方法，采用在线求解开环非线性最优控制，并不断更新控制输入构成闭环控制，通过步骤200建立航天器轨道交会的非线性最优控制问题，并进一步基于步骤300所提出的保辛数值方法求解非线性最优控制问题，最终可以在线快速解决航天器闭环交会的实时计算问题，得到最优航天器交会轨迹和推力，并满足实时性需求。本方法能够克服航天器在交会过程中由于导航误差、执行机构误差以及模型不确定性等，克服非线性开环控制由于上述因素导致航天器不能正常交会的局限性。

附图说明

图1为日地系统圆形限制性三体模型示意图；

图2为航天器在日地系统L2平动点附近南北Halo间交会图；

图3为本发明实施例提供的航天器在日地系统不稳定平动点轨道间交会控制方法的实现流程图。

具体实施方式

为使本发明解决的技术问题、采用的技术方案和达到的技术效果更加清楚，下面结合附图和实施例对本发明作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释本发明，而非对本发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与本发明相关的部分而非全部内容。

图3为本发明实施例提供的航天器在日地系统不稳定平动点轨道间交会控制方法的实现流程图。如图3所示，本发明实施例提供的航天器在日地系统不稳定平动点轨道间交会控制方法包括：

步骤100，确定航天器在日地系统平动点轨道交会的初始轨道、目标轨道以及交会过程时间，并建立受控航天器交会动力学模型。

具体选取的参数为：地球的质量与地球太阳质量之和的比值μ＝3.040360×10^-6，地球和平动点L₂之间的距离γ＝1.50767856×10⁶km，航天器的质量m＝1000kg。针对一个具体航天器轨道交会任务，可以选择初始轨道的z方向幅值为5×10⁵km，交会轨道的z方向幅值为7×10⁵km，交会过程时间选为100天，上述数值可以根据实际不同航天器交会情况修改。

具体的过程为：建立在日地系统两天体质心上的参考坐标系，参考坐标系与两天体以相同的轨道速度旋转，参考坐标系的x轴从系统质心处向地球延伸，z轴沿着系统的角动量方向延伸，y轴根据右手法则确定。

最终描述受控航天器交会动力学模型的动力学方程可以写成如下的无量纲形式：

$\left\{\begin{array}{c}\stackrel{\·\·}{x}-2\stackrel{\·}{y}-x=-\frac{(1-\mathrm{\μ})(x+1+1/\mathrm{\γ})}{{\mathrm{\γ}}^3{d}_1^3}-\frac{\mathrm{\μ}(x+1)}{{\mathrm{\γ}}^3{d}_2^3}+\frac{1-\mathrm{\μ}+\mathrm{\γ}}{\mathrm{\γ}}+u_x\\ \stackrel{\·\·}{y}+2\stackrel{\·}{x}-y=-\frac{(1-\mathrm{\μ})y}{{\mathrm{\γ}}^3{d}_1^3}-\frac{\mathrm{\μy}}{{\mathrm{\γ}}^3{d}_2^3}+u_y\\ \stackrel{\·\·}{z}=-\frac{(1-\mathrm{\μ})z}{{\mathrm{\γ}}^3{d}_1^3}-\frac{\mathrm{\μz}}{{\mathrm{\γ}}^3{d}_2^3}+u_z\end{array}\right.---\left(1\right)$

其中， $d_1=\sqrt{{(x+1+1/\mathrm{\γ})}^2+y^2+z^2},d_2=\sqrt{{(x+1)}^2+y^2+z^2},$ 方程中的点号代表旋转框架内的时间微分，μ是地球的质量与地球太阳质量之和的比值，γ表示地球和平动点L₂之间的距离，符号u_x、u_y和u_z分别表示受控航天器在x,y和_z方向的控制输入，表示坐标x对时间求导。

令x₁＝x，x₂＝y，x₃＝z，和将受控二阶动力学方程(1)式写成一阶状态空间形式如下：

$\left\{\begin{array}{c}{\stackrel{\·}{x}}_1=x_4\\ {\stackrel{\·}{x}}_2=x_5\\ {\stackrel{\·}{x}}_3=x_6\\ {\stackrel{\·}{x}}_4={2x}_5+x_1-\frac{(1-\mathrm{\μ})(x_1+1+1/\mathrm{\γ})}{{\mathrm{\γ}}^3{d}_1^3}-\frac{\mathrm{\μ}(x_1+1)}{{\mathrm{\γ}}^3{d}_2^3}+\frac{1-\mathrm{\μ}+\mathrm{\γ}}{\mathrm{\γ}}+u_x\\ {\stackrel{\·}{x}}_5=-2x_4+x_2-\frac{(1-\mathrm{\μ})x_2}{{\mathrm{\γ}}^3{d}_1^3}-\frac{\mathrm{\μ}{x}_2}{{\mathrm{\γ}}^3{d}_2^3}+u_y\\ {\stackrel{\·}{x}}_6=-\frac{(1-\mathrm{\μ})x_3}{{\mathrm{\γ}}^3{d}_1^3}-\frac{\mathrm{\μ}{x}_3}{{\mathrm{\γ}}^3{d}_2^3}+u_z\end{array}\right.---\left(2\right)$

进一步令状态变量x＝[x₁,x₂,x₃,x₄,x₅,x₆]^T，控制变量u＝[u_x,u_y,u_z]^T，则受控动力学方程(2)可以简化为

$\stackrel{\·}{x}=f(x\left(t\right),u\left(t\right))---\left(3\right)$

在本发明实施例中状态变量为航天器的位置和速度，控制变量代表的是航天器的推力加速度。

步骤200，综合考虑航天器在交会过程中的能量消耗以及交会精度，并基于建立的受控航天器交会动力学模型，建立非线性最优控制问题。

具体过程为：通过以下公式将航天器交会过程建立结合能量消耗最少和偏差最小的性能指标，采用加权性能指标方式建立非线性最优控制问题：

$\begin{array}{c}J=\frac{1}{2}{[M_{f}x(t+T)-\mathrm{\ψ}]}^T{S}_f[M_{f}x(t+T)-\mathrm{\ψ}]\\ +\frac{1}{2}{\∫}_t^{t+T}[{(x-x_d)}^{T}Q(x-x_d)+{(u-u_d)}^{T}G(u-u_d)]\mathrm{d\τ}\end{array}---\left(4\right)$

其中，x∈R^n×1表示状态变量，u∈R^m×1表示控制变量，x_d∈R^n×1表示目标状态向量，u_d∈R^m×1表示目标控制输入，t表示时间，τ∈[t,t+T]表示预测未来系统状态的时间变量，ψ∈R^m×1表示线性混合终端状态M_fx(t+T)的目标值，Q∈R^n×n表示半正定矩阵，G∈R^m×m表示正定矩阵，S_f表示半正定终端权矩阵，M_f表示给定的关于目标状态矩阵，根据最终设定的目标状态确定，例如可以选择单位矩阵。滚动时域控制问题的目标为寻找合适的控制输入u(t)使得性能指标J的值极小。在本发明中性能指标J综合了能量消耗最小和偏差最小。

采用拟线性化方法求解方程(3)和(4)，在k+1次迭代时将方程(3)和(4)在第k次迭代处进行泰勒级数展开，得到：

${\stackrel{\·}{x}}^{(k+1)}\left(\mathrm{\τ}\right)=A^{\left(k\right)}\left(\mathrm{\τ}\right)x^{(k+1)}\left(\mathrm{\τ}\right)+B^{\left(k\right)}\left(\mathrm{\τ}\right)u^{(k+1)}\left(\mathrm{\τ}\right)+w^{\left(k\right)}\left(\mathrm{\τ}\right),x^{(k+1)}(\mathrm{\τ}=t)=x\left(t\right)---\left(5\right)$

$\begin{array}{c}{J}^{(k+1)}=\frac{1}{2}{[M_f{x}^{(k+1)}(t+T)-\mathrm{\ψ}]}^T{S}_f[M_f{x}^{(k+1)}(t+T)-\mathrm{\ψ}]\\ +\frac{1}{2}{\∫}_t^{t+T}[{(x^{(k+1)}-x_d)}^{T}Q(x^{(k+1)}-x_d)+{(u^{(k+1)}-u_d)}^{T}G(u^{(k+1)}-u_d)]\mathrm{d\τ}\end{array}---\left(6\right)$

其中， $A^{\left(k\right)}\left(\mathrm{\τ}\right)=\frac{\∂f(x\left(\mathrm{\τ}\right),u\left(\mathrm{\τ}\right))}{\∂x}{|}_{x^{\left(k\right)}\left(\mathrm{\τ}\right),u^{\left(k\right)}\left(\mathrm{\τ}\right)},B^{\left(k\right)}\left(\mathrm{\τ}\right)=\frac{\∂f(x\left(\mathrm{\τ}\right),u\left(\mathrm{\τ}\right))}{\∂u}{|}_{x^{\left(k\right)}\left(\mathrm{\τ}\right),u^{\left(k\right)}\left(\mathrm{\τ}\right)},$

w^(k)(τ)＝f(x^(k)(τ),u^(k)(τ))-A^(k)(τ)x^(k)(τ)-B^(k)(τ)u^(k)(τ)；

采用变分法得到如下哈密顿两点边值问题：

$\left\{\begin{array}{c}{\stackrel{\·}{x}}^{(k+1)}=\frac{{\∂H}^{(k+1)}}{{\∂\mathrm{\λ}}^{(k+1)}}=A^{\left(k\right)}\left(\mathrm{\τ}\right)x^{(k+1)}-B^{\left(k\right)}\left(\mathrm{\τ}\right)G^{-1}{\left(B^{\left(k\right)}\left(\mathrm{\τ}\right)\right)}^T{\mathrm{\λ}}^{(k+1)}+B^{\left(k\right)}\left(\mathrm{\τ}\right)u_d+w^{\left(k\right)}\left(\mathrm{\τ}\right)\\ {\stackrel{\·}{\mathrm{\λ}}}^{(k+1)}=-\frac{{\∂H}^{(k+1)}}{\∂x^{(k+1)}}=-Q{x}^{(k+1)}-{\left(A^{\left(k\right)}\left(\mathrm{\τ}\right)\right)}^T{\mathrm{\λ}}^{(k+1)}+{\mathrm{Qx}}_d\end{array}\right.---\left(7\right)$

边界条件为：

x^(k+1)(τ＝t)＝x(t)  (8)

${\mathrm{\λ}}^{(k+1)}(\mathrm{\τ}=t+T)=M_f^T{S}_f[M_f{x}^{(k+1)}(t+T)-\mathrm{\ψ}]---\left(9\right)$

因此，方程(7)(8)(9)将航天器在不稳定平动点轨道族间交会的非线性最优控制问题通过迭代方式转化成序列线性最优控制问题。通过求解上述序列迭代线性最优控制问题，就可以得到在性能指标J最小的航天器轨道交汇所需能量消耗和最优交会轨迹。以下需要研究闭环反馈控制策略的构造以及相应的数值求解方法。

步骤300，采用保辛数值方法在有限长时间内求解非线性最优控制问题，得到协态变量。

步骤301，将交会过程时间[t,t+T]等分为N个子时间区间，每一个子时间区间的时间步长为η＝T/N，其中，每一个子时间区间表示为：

t₀＝t,t₁＝t+η,...,t_j＝t+jη,...,t_N＝t_f  (10)

步骤302，在第j个子时间区间中，获得开环最优控制系统的协态变量。

其中，j从0开始到N结束。具体过程为：在第j个子区间中，用插入了r个等间距点的r-1阶Lagrange多项式来近似代替状态变量x(τ)，而用插入了s个等间距点的s-1阶Lagrange多项式来近似代替协态变量λ(τ)，即

$x\left(\mathrm{\τ}\right)=(M\⊗I){\stackrel{\‾}{x}}_j---\left(11\right)$

$\mathrm{\λ}\left(\mathrm{\τ}\right)=N_1{\mathrm{\λ}}_{j-1}+(\underset{\‾}{N}\⊗I){\stackrel{\‾}{\mathrm{\λ}}}_j+N_n{\mathrm{\λ}}_j---\left(12\right)$

其中，符号代表Kronecker积，I代表n×n阶单位矩阵，向量λ_j-1与λ_j分别代表第j个子时间区间左端与右端的协态变量， ${\stackrel{\‾}{x}}_j={\{{\left({\stackrel{\‾}{x}}_j^1\right)}^T,{\left({\stackrel{\‾}{x}}_j^2\right)}^T,...,{\left({\stackrel{\‾}{x}}_j^r\right)}^T\}}^T,{\stackrel{\‾}{\mathrm{\λ}}}_j={\{{\left({\stackrel{\‾}{\mathrm{\λ}}}_j^2\right)}^T,{\left({\stackrel{\‾}{\mathrm{\λ}}}_j^3\right)}^T,...,{\left({\stackrel{\‾}{\mathrm{\λ}}}_j^{s-1}\right)}^T\}}^T.$ 其中，协态变量是在最优控制计算过程中引入的一个变量，其作用是为了求出控制推力变量。

引入对偶变量变分原理如下：

$\mathrm{\δ}\stackrel{\‾}{S}={\∫}_t^{t+T}{\left(\mathrm{\δx}\right)}^T(-\stackrel{\·}{\mathrm{\λ}}-\frac{\∂H}{\∂x})\mathrm{d\τ}+{\∫}_t^{t+T}{\left(\mathrm{\δ\λ}\right)}^T(\stackrel{\·}{x}-\frac{\∂H}{\∂\mathrm{\λ}})\mathrm{d\τ}+{\mathrm{\λ}}^T\mathrm{\δx}{|}_t^{t+T}=0---\left(13\right)$

将离散的状态变量(11)和协态变量(12)代入上述对偶变量变分原理，获得如下线性方程

AZ＝B  (14)

其中，Z＝{λ(t₀)^T,λ(t₁)^T,…,λ(t_j)^T,…,λ(t_N)^T}^T；

$B={\{{(x_0-{\mathrm{\ζ}}_1^1)}^T,-{({\mathrm{\ζ}}_2^1+{\mathrm{\ζ}}_1^2)}^T,...,-{({\mathrm{\ζ}}_2^j+{\mathrm{\ζ}}_1^{j+1})}^T,...,-{({\mathrm{\ζ}}_2^{N-1}+{\mathrm{\ζ}}_1^N)}^T,-M_f^T{S}_f\mathrm{\ψ}-M_f^T{S}_f{M}_f{\mathrm{\ζ}}_2^N\}}^T;$

利用公式(14)得到开环最优控制系统的协态变量。下一步需要根据协态变量求解控制输入。

步骤400，利用获得的协态变量更新当前子时间区间的控制输入，将得到的控制输入作用到航天器交会动力学模型，完成一个时间区间的动力学仿真，并通过导航方法获得当前航天器状态变量。

具体的过程为：可以将得到控制输入通过航天器的发动机等执行机构实施到航天器交会动力学模型。通过以下公式进行更新控制输入：

u＝u_d-G^-1(B(τ))^Tλ  (15)

其中，其中，u∈R^m×1表示控制变量，u_d∈R^m×1表示目标控制输入，G∈R^m×m表示正定矩阵，τ表示时间，B(τ)表示动力学方程对控制变量的偏导数矩阵，λ表示协态变量，将当前控制作用到受控的圆形限制性三体模型进行下一个时间步递进。

本发明通过导航方法获得当前航天器状态变量，例如，可以采用卡尔曼滤波或者扩展卡尔曼滤波等航天器导航算法确定当前航天器状态。

步骤500，进行子时间区间的递进，以当前航天器终端状态为下一个子时间区间的初始状态，依次重复步骤302至步骤400，直至完成航天器在不稳定平动点轨道间交会任务。

迭代的过程实际上是时间的递推，通过求解当前时间子区间(时间步)的状态，可以得到航天器当前的位置和速度，将当前的终端状态作为下一时间步的初始状态，可以得到下一时间步的状态，不断重复计算直至完成航天器在不稳定平动点轨道间交会任务。以上过程之后，可以得到航天器交会过程中的最优推力和最优轨迹。在每一个时间步中，保辛方法的作用是求解协态变量，协态变量的作用是求解控制变量，得到控制变量，再和导航方法得到的状态变量，将当前时间步的终端状态变量和控制变量作为下一时间步的初始变量。

本实施例提供的航天器在日地系统不稳定平动点轨道间交会控制方法，采用在线求解开环非线性最优控制，并不断更新控制输入构成闭环控制，通过步骤200建立航天器轨道交会的非线性最优控制问题，并进一步基于步骤300所提出的保辛数值方法求解非线性最优控制问题，最终可以在线快速解决航天器闭环交会的实时计算问题，得到最优航天器交会轨迹和推力，并满足实时性需求。本方法能够克服航天器在交会过程中由于导航误差、执行机构误差以及模型不确定性等，克服非线性开环控制由于上述因素导致航天器不能正常交会的局限性。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (5)

1.一种航天器在日地系统不稳定平动点轨道间交会控制方法，其特征在于，所述方法包括以下步骤：

步骤100，确定航天器在日地系统平动点轨道交会的初始轨道、目标轨道以及交会过程时间，并建立受控航天器交会动力学模型；

步骤200，综合考虑航天器在交会过程中的能量消耗以及交会精度，并基于建立的受控航天器交会动力学模型，建立非线性最优控制问题；

步骤300，采用保辛数值方法在有限长时间内求解非线性最优控制问题，得到协态变量，包括以下步骤：

步骤301，将交会过程时间[t,t+T]等分为N个子时间区间，每一个子时间区间的时间步长为η＝T/N，其中，每一个子时间区间表示为t₀＝t,t₁＝t+η,...,t_j＝t+jη,...,t_N＝t_f；

步骤302，在第j个子时间区间中，获得开环最优控制系统的协态变量，其中，j从0开始到N结束；

步骤400，利用获得的协态变量更新当前子时间区间的控制输入，将得到控制输入作用到航天器交会动力学模型，完成一个子时间区间的动力学仿真，并通过导航方法获得当前航天器状态变量；

步骤500，进行子时间区间的递进，以当前航天器终端状态为下一个子时间区间的初始状态，依次重复步骤302至步骤400，直至完成航天器在不稳定平动点轨道间交会任务。

2.根据权利要求1所述的航天器在日地系统不稳定平动点轨道间交会控制方法，其特征在于，建立受控航天器交会动力学模型，包括：

建立在日地系统两天体质心上的参考坐标系，参考坐标系与两天体以相同的轨道速度旋转，参考坐标系的x轴从系统质心处向地球延伸，z轴沿着系统的角动量方向延伸，y轴根据右手法则确定；

受控航天器交会动力学模型的动力学方程表示为：

$\left\{\begin{array}{c}\stackrel{..}{x}-2\stackrel{.}{y}-x=-\frac{(1-\mathrm{\μ})(x+1+1/\mathrm{\γ})}{{\mathrm{\γ}}^3{d}_1^3}-\frac{\mathrm{\μ}(x+1)}{{\mathrm{\γ}}^3{d}_2^3}+\frac{1-\mathrm{\μ}+\mathrm{\γ}}{\mathrm{\γ}}+u_x\\ \stackrel{..}{y}+2\stackrel{.}{x}-y=-\frac{(1-\mathrm{\μ})y}{{\mathrm{\γ}}^3{d}_1^3}-\frac{\mathrm{\μy}}{{\mathrm{\γ}}^3{d}_2^3}+u_y\\ \stackrel{..}{z}=-\frac{(1-\mathrm{\μ})z}{{\mathrm{\γ}}^3{d}_1^3}-\frac{\mathrm{\μz}}{{\mathrm{\γ}}^3{d}_2^3}+u_z\end{array}\right.---\left(1\right)$

其中， $d_1=\sqrt{{(x+1+1/\mathrm{\γ})}^2+y^2+z^2},d_2=\sqrt{{(x+1)}^2+y^2+z^2},$ 方程中的点号代表旋转框架内的时间微分，μ是地球的质量与地球太阳质量之和的比值，γ表示地球和平动点L₂之间的距离，符号u_x、u_y和u_z分别表示受控航天器的控制输入，表示坐标x对时间求导。

3.根据权利要求1所述的航天器在日地系统不稳定平动点轨道间交会控制方法，其特征在于，步骤200中建立的非线性最优控制问题的性能指标表示为：

$\begin{array}{c}J=\frac{1}{2}{[M_{f}x(t+T)-\mathrm{\ψ}]}^T{S}_f[M_{f}x(t+T)-\mathrm{\ψ}]\\ +\frac{1}{2}{\∫}_t^{t+T}[{(x-x_d)}^{T}Q(x-x_d)+{(u-u_d)}^{T}G(u-u_d)]\mathrm{d\τ}\end{array}---\left(4\right)$

其中，x∈R^n×1表示状态变量，u∈R^m×1表示控制变量，x_d∈R^n×1表示目标状态向量，u_d∈R^m×1表示目标控制输入，t表示时间，τ∈[t,t+T]表示预测未来系统状态的时间变量，ψ∈R^m×1表示线性混合终端状态M_fx(t+T)的目标值，Q∈R_n×n表示半正定矩阵，G∈R^m×m表示正定矩阵，S_f表示半正定终端权矩阵，M_f表示给定的矩阵。

4.根据权利要求1所述的航天器在日地系统不稳定平动点轨道间交会控制方法，其特征在于，获得开环最优控制系统的协态变量，包括：

在第j个子区间中，用插入了r个等间距点的r-1阶Lagrange多项式来代替状态变量x(τ)，而用插入了s个等间距点的s-1阶Lagrange多项式来代替协态变量λ(τ)，状态变量x(τ)和协态变量λ(τ)表示为：

$x\left(\mathrm{\τ}\right)=(M\⊗I){\stackrel{\‾}{x}}_j---\left(11\right)$

$\mathrm{\λ}\left(\mathrm{\τ}\right)=N_1{\mathrm{\λ}}_{j-1}+(\underset{\‾}{N}\⊗I){\stackrel{\‾}{\mathrm{\λ}}}_j+N_n{\mathrm{\λ}}_j---\left(12\right)$

其中，符号代表Kronecker积，I代表n×n阶单位矩阵，向量λ_j-1与λ_j分别代表第j个子时间区间左端与右端的协态变量， ${\stackrel{\‾}{x}}_j={\{{\left({\stackrel{\‾}{x}}_j^1\right)}^T,{\left({\stackrel{\‾}{x}}_j^2\right)}^T,\·\·\·,{\left({\stackrel{\‾}{x}}_j^r\right)}^T\}}^T,{\stackrel{\‾}{\mathrm{\λ}}}_j={\{{\left({\stackrel{\‾}{\mathrm{\λ}}}_j^2\right)}^T,{\left({\stackrel{\‾}{\mathrm{\λ}}}_j^3\right)}^T,\·\·\·,{\left({\stackrel{\‾}{\mathrm{\λ}}}_j^{s-1}\right)}^T\}}^T;$

引入对偶变量变分原理如下：

$\mathrm{\δ}\stackrel{\‾}{S}={\∫}_t^{t+T}{\left(\mathrm{\δx}\right)}^T(-\stackrel{.}{\mathrm{\λ}}-\frac{\∂H}{\∂x})\mathrm{d\τ}+{\∫}_t^{t+T}{\left(\mathrm{\δ\λ}\right)}^T(\stackrel{.}{x}-\frac{\∂H}{\∂\mathrm{\λ}})\mathrm{d\τ}+{\mathrm{\λ}}^T\mathrm{\δx}{|}_t^{t+T}=0---\left(13\right)$

将离散的状态变量(11)和协态变量(12)代入上述对偶变量变分原理，获得如下线性方程

AZ＝B                  (14)

其中，Z＝{λ(t₀)^T,λ(t₁)^T,…,λ(t_j)^T,…,λ(t_N)^T}^T；

$B={\{{(x_0-{\mathrm{\ζ}}_1^1)}^T,-{({\mathrm{\ζ}}_2^1+{\mathrm{\ζ}}_1^2)}^T,\·\·\·,-{({\mathrm{\ζ}}_2^j+{\mathrm{\ζ}}_1^{j+1})}^T,\·\·\·,-{({\mathrm{\ζ}}_2^{N-1}+{\mathrm{\ζ}}_1^N)}^T,-M_f^T{S}_f\mathrm{\ψ}-M_f^T{S}_f{M}_f{\mathrm{\ζ}}_2^N\}}^T;$

利用公式(14)得到开环最优控制系统的协态变量。

5.根据权利要求1所述的航天器在日地系统不稳定平动点轨道间交会控制方法，其特征在于，通过以下公式利用获得的协态变量更新当前子时间区间的控制输入：

u＝u_d-G^-1(B(τ))^Tλ                   (15)

其中，u∈R^m×1表示控制变量，u_d∈R^m×1表示目标控制输入，G∈R^m×m表示能量加权矩阵，是正定矩阵，τ表示时间，B(τ)表示受控非线性动力学方程对控制变量的偏导数矩阵，λ表示协态变量。