近年来，关于平动点任务的研究受到越来越多的关注。本文讨论的平动点属于三体问题范畴，是天体力学领域的基本问题之一。Euler和Lagrange分别在1765年和1772年研究限制性三体问题(Restricted Three Body Problem, RTBP)时，发现RTBP存在5个平衡点，称之为平动点或拉格朗日点^[1]，其中，位于两个主天体连线上的3个平动点称为共线平动点。由于平动点独特的动力学特性，将航天器部署到平动点上时，航天器在引力作用下，与主天体的相对位置保持不变。然而在真实空间环境下，平动点位置无法被精确确定，而相比于平动点，其附近的周期轨道更适用于航天任务。因此实际任务中通常将航天器部署在平动点附近的轨道上。

平动点的动力学内涵十分丰富^[2]，在其附近存在着丰富的轨道类型，包括周期轨道(如Lyapunov轨道、Halo轨道等)、拟周期轨道(如Lissajous轨道、拟Halo轨道等)以及不变流形等^[3-5]。这些轨道为一些特殊的深空探测任务以及空间低能转移提供了可能。例如，日地L₂点因其独特的空间位置和物理环境，为空间望远镜等深空观测卫星提供了良好的部署条件^[6]；日地L₁点位于日地连线之间，是部署太阳活动监测卫星的绝佳地点^[7]；地月L₂点位于地月连线上月球的背面，选择合适的轨道尺寸能够使部署在上面的航天器有效规避月掩，是为地球与月球背面提供中继通信的理想场所^[8]；此外，平动点在太阳系中广泛存在，这些平动点附近的不变流形相互连接，构成了一系列蜿蜒错综的管道，称为行星际高速公路(InterPlanetary Superhighway, IPS)，利用IPS，可以实现行星际的低能转移^[9]。

平动点轨道在实际任务中的应用起步相对较晚。1978年，美国航空航天局(NASA)成功发射ISEE-3探测器，是平动点轨道第1次在实际任务中得到应用。以ISEE-3的成功实践为起点，拉开了利用平动点进行深空探测的序幕^[10]。在接下来的几十年里，以NASA和欧洲航天局(ESA)为代表的航天机构，先后成功实施了SOHO、ACE、MAP、Genesis、ARTEMIS等多次平动点任务^[11]，在太阳活动监测、深空观测等领域取得了丰硕的成果。中国先后进行了嫦娥二号、嫦娥五号T1以及“鹊桥”卫星3次平动点任务，其中于2018年成功发射的“鹊桥”通信中继卫星实现了人类历史上首次月背中继通信，为嫦娥四号进行人类首次月背软着陆探测提供了有力支撑。

虽然共线平动点附近存在着丰富的轨道类型，但由于其本身是不稳定的，运行在其附近轨道上的航天器若不施加控制，将会很快偏离标称轨道。因此，设计平动点任务时，进行平动点轨道维持研究十分必要。针对平动点轨道维持问题，许多学者进行了研究，并取得了一些成果。Simo等^[12]利用Floquet模态理论针对拟Halo轨道和Halo轨道的维持问题进行了研究；Breakwell等^[13]最早利用线性二次型调节器(LQR)方法实现了对Halo轨道的维持控制；Howell和Pernicka^[14]提出了利用靶点法维持平动点轨道的策略，并在ARTEMIS任务中得到了应用；Lian等^[15]利用离散时间滑模控制实现了对真实星历下地月系平动点轨道的维持控制；Nazari等^[16]分别利用时变LQR、退步控制和反馈线性化实现了地月L₁附近Halo轨道的维持控制，并对3种方法的仿真结果进行了对比分析；徐明和徐世杰^[17]设计了一种线性周期控制策略实现了Halo轨道的维持控制。

基于特征模型的自适应控制方法是由吴宏鑫院士在20世纪80年代初提出的一种从实际应用角度出发的控制方法。经过30多年的发展，不仅在理论上取得了重要进展，同时还成功应用于神舟飞船再入返回控制、交会对接等重大工程中^[18]。所谓特征模型，即结合对象动力学特征、环境特征和控制性能要求而不是仅以对象精确动力学分析所建立的模型，其具有以下特点^[19]：①在同样控制作用下，对象特征模型和实际对象在输出上是等价的，在稳定情况下，输出是相等的；②特征模型的形式和阶次除考虑对象特征外，主要取决于控制性能要求；③特征模型建立的形式比原对象的动力学方程简单，易于控制器设计，工程实现容易、方便；④特征模型与高阶系统的降阶模型不同，它是把高阶模型有关信息都压缩到几个特征参量中，并不丢失信息，一般情况下用慢时变差分方程描述。

由于三体问题相对于二体问题更为复杂，而基于特征模型的控制方法具有将复杂问题简单化的优势，因此本文基于特征模型理论，研究了地月系L₂点附近Halo轨道的维持控制问题。本文结构安排如下：第1节介绍两类动力学模型及标称轨道的获取；第2节介绍黄金分割控制器和PD控制器的设计过程；第3节介绍仿真结果及结果分析；第4节总结全文，给出结论。

1 动力学模型与标称轨道 1.1 圆型限制性三体模型

圆型限制性三体问题(Circular Restricted Three Body Problem, CRTBP)是最简化形式的三体问题，其可以描述为：一个质量无穷小的天体在两个大质量天体的引力作用下的运动，其中两个大质量天体(又称为主天体)围绕其公共质心做圆形周期运动。在航天应用领域，常涉及的主天体主要是地球-月球或者太阳-地球/月球。

通常在如图 1所示的会合坐标系下对CRTBP进行研究。图中:m₃为航天器; M₁和M₂为两个主天体(M₁≥M₂); 坐标原点O为它们的公共质心，Ox由M₁指向M₂，Oz指向M₁和M₂运动的角动量方向，xOy平面为M₁和M₂的运动平面，Oxyz构成右手坐标系；r、r₁、r₂分别为由O、M₁和M₂指向m₃的矢径。为了方便计算，对各物理量进行无量纲化处理。归一化后的质量、长度、时间单位分别为

$ \left\{ {\begin{array}{*{20}{l}} {\left[ M \right] = {M_1} + {M_2}}\\ {\left[ L \right] = {L_{12}}}\\ {\left[ T \right] = {{\left[ {\frac{{L_{12}^3}}{{G\left( {{M_1} + {M_2}} \right)}}} \right]}^{\frac{1}{2}}}} \end{array}} \right. $	（1）

式中：等式右边的各物理量均采用国际单位制，M₁和M₂分别为两个主天体的质量；L₁₂为主天体之间的距离；G为万有引力常量，则可得到归一化后的主天体M₁和M₂的质量分别为1－μ和μ，μ为质量参数，形式为

$ \mu = \frac{{{M_2}}}{{{M_1} + {M_2}}} $	（2）

主天体在会合坐标系中的坐标分别为(－μ, 0, 0)和(1－μ, 0, 0)。在地月系中，μ=0.012 150 585 61。

令[x y z]^T为归一化条件下航天器的位置坐标，可以得到CRTBP中m₃在会合坐标系中的动力学方程为

$ \left\{ {\begin{array}{*{20}{l}} {\ddot x - 2\dot y = \frac{{\partial \mathit{\Omega }}}{{\partial x}}}\\ {\ddot y + 2\dot x = \frac{{\partial \mathit{\Omega }}}{{\partial y}}}\\ {\ddot z = \frac{{\partial \mathit{\Omega }}}{{\partial z}}} \end{array}} \right. $	（3）

式中：Ω为有效势，其形式为

$ \mathit{\Omega } = \frac{1}{2}\left( {{x^2} + {y^2}} \right) + \frac{{1 - \mu }}{{{r_1}}} + \frac{\mu }{{{r_2}}} $	（4）

其中：r₁=||r₁||和r₂=||r₂||分别为航天器与两个主天体之间的距离：

$ \left\{ {\begin{array}{*{20}{l}} {{r_1} = {{\left[ {{{(x + \mu )}^2} + {y^2} + {z^2}} \right]}^{\frac{1}{2}}}}\\ {{r_2} = {{\left[ {{{(1 - x - \mu )}^2} + {y^2} + {z^2}} \right]}^{\frac{1}{2}}}} \end{array}} \right. $	（5）

1.2 双圆限制性四体模型

CRTBP是最理想情况下的简化模型，在实际情况中，航天器在平动点轨道上的运动会受到第四体引力摄动、系统偏心率、太阳光压等摄动因素的影响。对于地月三体系统，太阳引力摄动是影响航天器运动的主要摄动源，因此在CRTBP基础上，引入太阳引力摄动，建立双圆限制性四体模型(BiCircular Restricted Four Body Model, BCRFBM)能够更加精确地模拟真实力学环境。

双圆限制性四体坐标系如图 2所示，在BCRFBM假设下，太阳在地月运动平面内绕系统质心做匀速圆周运动。图中:M₄为太阳，坐标轴定义与会合坐标系相同；r、r₁、r₂、r₄分别为由O、M₁、M₂和M₄指向m₃的矢径；d₄为O指向太阳的矢径，θ₄为d₄与坐标系x轴的夹角。

在BCRFBM假设下，航天器m₃的动力学方程为

$ \left\{ {\begin{array}{*{20}{l}} {\ddot x - 2\dot y = \frac{{\partial {\mathit{\Omega }_4}}}{{\partial x}}}\\ {\ddot y + 2\dot x = \frac{{\partial {\mathit{\Omega }_1}}}{{\partial y}}}\\ {\ddot z = \frac{{\partial {\mathit{\Omega }_4}}}{{\partial z}}} \end{array}} \right. $	（6）

式中：

$ \left\{ \begin{array}{l} {\mathit{\Omega }_4} = \frac{1}{2}\left( {{x^2} + {y^2}} \right) + \frac{{1 - \mu }}{{{r_1}}} + \frac{\mu }{{{r_2}}} + \frac{{{\mu _4}}}{{{r_4}}} - \\ \;\;\;\;\;\;\frac{{{\mu _4}}}{{d_4^2}}\left( {x\cos {\theta _4} + y\sin {\theta _4}} \right)\\ {r_4} = {\left[ {{{\left( {x - {d_4}\cos {\theta _4}} \right)}^2} + {{\left( {y - {d_4}\sin {\theta _4}} \right)}^2} + {z^2}} \right]^{\frac{1}{2}}}\\ {\theta _4} = {\omega _4}t + {\theta _{40}} \end{array} \right. $

（7）

其中:r₄=||r₄||; μ₄=328 900.54为归一化单位下太阳的质量; ω₄=0.925 2为太阳相对系统质心的角速度; θ₄₀为初始时刻太阳相位角; d₄=388.811 4为归一化单位下太阳与系统质心的平均距离，其他归一化单位下的参数与CRTBP下相同。

1.3 标称Halo轨道

标称轨道是Halo轨道维持研究的基础。Halo轨道是共线平动点附近的一类周期轨道，是水平Lyapunov轨道分岔的结果。由于三体动力学系统具有强非线性特征，无法获得Halo轨道的精确解析解，因此标称轨道通常采用数值积分的方式求解。

数值积分的初值，通常通过修正由近似解析解提供的近似初值获得。Richardson^[20]在1980年基于Lindstedt-Poincaré方法推导了Halo轨道的三阶近似解析解，其表达式为

$ \left\{ \begin{array}{l} x = {a_{21}}A_x^2 + {a_{22}}A_z^2 - {A_x}\cos {\tau _1} + \left( {{a_{23}}A_x^2 - } \right.\\ \;\;\;\;\;\left. {{a_{24}}A_z^2} \right)\cos 2{\tau _1} + \left( {{a_{31}}A_x^3 - {a_{32}}{A_x}A_z^2} \right)\cos 3{\tau _1}\\ y = k{A_x}\sin {\tau _1} + \left( {{b_{21}}A_x^2 - {b_{22}}A_z^2} \right)\sin 2{\tau _1} + \\ \;\;\;\;\;\;\left( {{b_{31}}A_x^3 - {b_{32}}{A_x}A_z^2} \right)\sin 3{\tau _1}\\ z = {\delta _n}\left[ {{A_z}\cos {\tau _1} + {d_{21}}{A_x}{A_z}\left( {\cos 2{\tau _1} - 3} \right) + } \right.\\ \;\;\;\;\;\left. {\left( {{d_{32}}{A_z}A_x^2 - {d_{31}}A_z^3} \right)\cos 3{\tau _1}} \right] \end{array} \right. $

（8）

式中：τ₁为Halo轨道的线性化相位角，τ₁=τ₀+ωt，τ₀为初始相位角；A_x、A_z分别为Halo轨道在平动点旋转坐标系中的线性化振幅，且A_x和A_z之间存在限制性条件：

$ {f_1}A_x^2 + {f_2}A_z^2 + \Delta = 0 $	（9）

式(8)和式(9)中的其他各参数定义见文献[20]。

利用式(8)获得的积分初值只是近似初值，直接利用此初值代入动力学方程式(3)中进行数值积分，结果将很快发散。因此，还需对近似初值进行修正，以获得满足精度要求的积分初值。

对近似初值的修正可以采用微分修正的方法进行。由于Halo轨道具有关于xOz平面对称的几何性质，因此选取Halo轨道与xOz平面的一个交点作为轨道初值，记为

$ {\mathit{\boldsymbol{X}}_0} = {\left[ {\begin{array}{*{20}{l}} {{x_0}}&0&{{z_0}}&0&{{{\dot y}_0}}&0 \end{array}} \right]^{\rm{T}}} $

若初值为精确值，则积分半个轨道周期T/2时，轨道第一次穿越xOz平面，且有

$ {\mathit{\boldsymbol{X}}_{\rm{d}}} = {\left[ {\begin{array}{*{20}{l}} {{x_{\rm{d}}}}&0&{{z_{\rm{d}}}}&0&{{{\dot y}_{\rm{d}}}}&0 \end{array}} \right]^{\rm{T}}} $

但由于初值不准确，导致实际Halo轨道穿越xOz平面时，X_d值会偏离精确值。通过调整初值x₀、z₀和${{\dot y}_0}$，使得

$ {y_{\rm{d}}} = 0,{{\dot x}_{\rm{d}}} = 0,{{\dot z}_{\rm{d}}} = 0 $

可达到修正初值的效果。实际修正中，可固定z₀，简化修正过程，则微分修正量可通过如下形式获得^[21]：

$ \begin{gathered} \left[ {\begin{array}{*{20}{c}} {\delta {{\dot x}_{\text{d}}}} \\ {\delta {{\dot z}_{\text{d}}}} \end{array}} \right] = \left( {\left[ {\begin{array}{*{20}{c}} {{\mathit{\Phi }_{41}}}&{{\mathit{\Phi }_{45}}} \\ {{\mathit{\Phi }_{61}}}&{{\mathit{\Phi }_{65}}} \end{array}} \right] - } \right. \hfill \\ \frac{1}{{{{\dot y}_0}}}\left[ {\begin{array}{*{20}{c}} {{{\ddot x}_{\text{d}}}} \\ {{{\ddot z}_{\text{d}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{\mathit{\Phi }_{21}}}&{{\mathit{\Phi }_{25}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\delta }{x_0}} \\ {{\delta }{y_0}} \end{array}} \right] \hfill \\ \end{gathered} $

（10）

式中：$\delta {{\dot x}_{\rm{d}}}、\delta {{\mathit{\dot z}}_{\rm{d}}}$为穿越点处速度偏差值；δx₀和δy₀为初值的修正量；Φ_ij为CRTBP下系统状态转移矩阵Φ(t, t₀)中的元素。

设置Halo轨道z向振幅为A_z=0.016 6，依据微分修正方法获得了修正后的积分初值，然后将此初值代入动力学方程式(3)，数值积分可得如图 3所示的Halo轨道，其轨道周期在归一化单位下为T=3.412 2，相当于约14.738 6天，图中L_EM表示地球与月球之间的距离。

值得注意的是，由于Halo轨道本身是不稳定的，经过微分修正初值后获得的Halo轨道也不是完全闭合的轨道，经过多圈积分之后仍会发散。考虑到实际任务需求，仅依靠微分修正后的Halo轨道作为标称轨道显然无法完成较长时间周期内的轨道维持，因此还需对微分修正后的轨道做进一步的修正，获得一条持续时间较长的标称轨道。

为获得满足需求的标称轨道，本文采用靶点法对Halo轨道进行了进一步的修正。靶点法的基本思想是通过最小化目标函数J，求解使轨道维持在标称轨道附近所需的最小速度脉冲，其中目标函数J是机动速度和靶点状态相对标称轨道偏差的加权和，具有以下形式：

$ J = \Delta {\mathit{\boldsymbol{v}}^{\rm{T}}}\mathit{\boldsymbol{Q}}\Delta \mathit{\boldsymbol{v}} + \sum\limits_{k = 1}^n {\left( {\mathit{\boldsymbol{p}}_k^{\rm{T}}{\mathit{\boldsymbol{R}}_k}{\mathit{\boldsymbol{p}}_k} + \mathit{\boldsymbol{v}}_k^{\rm{T}}{\mathit{\boldsymbol{T}}_k}{\mathit{\boldsymbol{v}}_k}} \right)} $	（11）

式中：Q、R_k和T_k为权值矩阵；p_k和v_k为航天器在第k个靶点处相对标称轨道的位置和速度偏差，可利用状态转移矩阵Φ(t_k, t_c)获得

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{p}}_k}}\\ {{\mathit{\boldsymbol{v}}_k}} \end{array}} \right] = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{t_k},{t_c}} \right)\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{p}}_{\rm{c}}}}\\ {{\mathit{\boldsymbol{v}}_{\rm{c}}} + \Delta \mathit{\boldsymbol{v}}} \end{array}} \right] = \\ \;\;\;\;\;\;\;\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{A}}_{\rm{c}}}}&{{\mathit{\boldsymbol{B}}_{\rm{c}}}}\\ {{\mathit{\boldsymbol{C}}_{\rm{c}}}}&{{\mathit{\boldsymbol{D}}_{\rm{c}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{p}}_{\rm{c}}}}\\ {{\mathit{\boldsymbol{v}}_{\rm{c}}} + \Delta \mathit{\boldsymbol{v}}} \end{array}} \right] \end{array} $

（12）

式中：p_c和v_c为当前时刻t_c的位置和速度偏差。将p_k和v_k代入式(11)中，令

$ \frac{{\partial J}}{{\partial \Delta \mathit{\boldsymbol{v}}}} = {\bf{0}} $	（13）

即可求得Δv_min。

根据上述原理，以微分修正所获得的Halo轨道第一圈为基础，设计三靶点的修正策略，可以得到图 4所示的标称轨道。修正过程中，修正间隔为0.01 T，持续20个周期，图 5给出了速度修正量相对于速度的比例，可以看出除个别情况外，修正量均不超过速度的0.04%，因此可以认为标称轨道是连续的。

2 控制器设计

令$\boldsymbol{X}=\left[\begin{array}{llllll}{x} & {y} & {z} & {\dot{x}} & {\dot{y}} & {\dot{z}}\end{array}\right]^{\mathrm{T}}$，动力学方程式(3)可写为

$ \mathit{\boldsymbol{\dot X}} = F\left( \mathit{\boldsymbol{X}} \right) + \mathit{\boldsymbol{Bu}} $	（14）

式中：

$ F\left( \mathit{\boldsymbol{X}} \right) = \left[ {\begin{array}{*{20}{c}} {\dot x}\\ {\dot y}\\ {\dot z}\\ {2\dot y + \frac{{\partial \mathit{\Omega }}}{{\partial x}}}\\ { - 2\dot x + \frac{{\partial \mathit{\Omega }}}{{\partial y}}}\\ {\frac{{\partial \mathit{\Omega }}}{{\partial z}}} \end{array}} \right] $

$ \mathit{\boldsymbol{u}} = {\left[ {\begin{array}{*{20}{l}} {{u_1}}&{{u_2}}&{{u_3}} \end{array}} \right]^{\rm{T}}} $

$ \mathit{\boldsymbol{B}} = {\left[ {\begin{array}{*{20}{l}} {{{\bf{0}}_{3 \times 3}}}&{{\mathit{\boldsymbol{I}}_{3 \times 3}}} \end{array}} \right]^{\rm{T}}} $

其中：u称为控制加速度。

2.1 速度跟踪的多变量黄金分割控制器

根据文献[19]，对多输入多输出系统建立特征模型，需输入输出同维。由于式(14)中控制量为3维，而为实现Halo轨道维持需要对全状态(6维)进行跟踪，因此无法对全状态建立特征模型，故考虑仅对速度跟踪特征建模。

在采样周期Δt满足一定条件下，速度跟踪的特征模型可用如下输出解耦型二阶差分方程组描述：

$ \begin{array}{l} {\mathit{\boldsymbol{X}}_1}\left( {k + 1} \right) = {\mathit{\boldsymbol{F}}_1}\left( k \right){\mathit{\boldsymbol{X}}_1}\left( k \right) + {\mathit{\boldsymbol{F}}_2}\left( k \right){\mathit{\boldsymbol{X}}_1}\left( {k - 1} \right) + \\ {\mathit{\boldsymbol{G}}_0}\left( k \right){\mathit{\boldsymbol{U}}_1}\left( k \right) + {\mathit{\boldsymbol{G}}_1}\left( k \right){\mathit{\boldsymbol{U}}_1}\left( {k - 1} \right) \end{array} $

（15）

式中：

$ {\mathit{\boldsymbol{X}}_1}\left( k \right) = {\left[ {\begin{array}{*{20}{c}} {{x_1}\left( k \right)}&{{x_2}\left( k \right)}&{{x_3}\left( k \right)} \end{array}} \right]^{\rm{T}}} $

$ {\mathit{\boldsymbol{U}}_1}\left( k \right) = {\left[ {\begin{array}{*{20}{l}} {{u_{11}}\left( k \right)}&{{u_{12}}\left( k \right)}&{{u_{13}}\left( k \right)} \end{array}} \right]^{\rm{T}}} $

$ {\mathit{\boldsymbol{F}}_1}\left( k \right) = \left[ {\begin{array}{*{20}{c}} {{f_{11}}\left( k \right)}&0&0\\ 0&{{f_{12}}\left( k \right)}&0\\ 0&0&{{f_{13}}\left( k \right)} \end{array}} \right] $

$ {\mathit{\boldsymbol{F}}_2}\left( k \right) = \left[ {\begin{array}{*{20}{c}} {{f_{21}}\left( k \right)}&0&0\\ 0&{{f_{22}}\left( k \right)}&0\\ 0&0&{{f_{23}}\left( k \right)} \end{array}} \right] $

$ {\mathit{\boldsymbol{G}}_0}\left( k \right) = \left[ {\begin{array}{*{20}{c}} {{g_{0,11}}\left( k \right)}&{{g_{0,12}}\left( k \right)}&{{g_{0,13}}\left( k \right)}\\ {{g_{0,21}}\left( k \right)}&{{g_{0,22}}\left( k \right)}&{{g_{0,23}}\left( k \right)}\\ {{g_{0,31}}\left( k \right)}&{{g_{0,32}}\left( k \right)}&{{g_{0,33}}\left( k \right)} \end{array}} \right] $

$ {\mathit{\boldsymbol{G}}_1}\left( k \right) = \left[ {\begin{array}{*{20}{c}} {{g_{1,11}}\left( k \right)}&{{g_{1,12}}\left( k \right)}&{{g_{1,13}}\left( k \right)}\\ {{g_{1,21}}\left( k \right)}&{{g_{1,22}}\left( k \right)}&{{g_{1,23}}\left( k \right)}\\ {{g_{1,31}}\left( k \right)}&{{g_{1,32}}\left( k \right)}&{{g_{1,33}}\left( k \right)} \end{array}} \right] $

且有当Δt→0时

$ {f_{1j}}\left( k \right) \to 2\;\;\;j = 1,2,3 $

$ {f_{2j}}\left( k \right) \to - 1\;\;\;\;j = 1,2,3 $

$ {g_{0,jh}}\left( k \right) \to 0\;\;\;\;j,h = 1,2,3 $

$ {g_{1,jh}}\left( k \right) \to 0\;\;\;\;j,h = 1,2,3 $

$ \sum C = \sum\limits_{l = 1}^2 {{f_{ij}}} \left( k \right) + \sum\limits_{h = 1}^3 {\left[ {{g_{0,jh}}\left( k \right) + {g_{1,jh}}\left( k \right)} \right]} = 1 $

式(15)中的特征参量F₁(k)、F₂(k)、G₀(k)和G₁(k)由参数估计得到。尽管是多变量，对于式(15)的特征模型仍可按每路独立进行参数估计^[19]。每一回路的特征方程可写为

$ \begin{array}{l} {x_j}\left( k \right) = {f_{1j}}\left( k \right){x_j}(k - 1) + {f_{2j}}(k){x_j}(k - 2) + \\ \;\;\;\;\;\;\;\;\;\sum\limits_{h = 1}^3 {{g_{0,jh}}} (k){u_{1h}}(k - 1) + \\ \;\;\;\;\;\;\;\;\;\sum\limits_{h = 1}^3 {{g_{1,jh}}} (k){u_{1h}}(k - 2) = \phi _j^{\rm{T}}(k - 1){\mathit{\boldsymbol{\widehat \theta }}_j}(k) \end{array} $	（16）

式中：

$ \begin{array}{l} {\phi _j}(k - 1) = \left[ {{x_j}(k - 1)\quad {x_j}(k - 2)\quad {u_{11}}(k - 1)} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{u_{12}}(k - 1)\quad {u_{13}}(k - 1)\quad {u_{11}}(k - 2)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left. {{u_{12}}(k - 2)\quad {u_{12}}(k - 2)} \right]^{\rm{T}}} \end{array} $

$ \begin{array}{l} {{\mathit{\boldsymbol{\hat \theta }}}_j}\left( k \right) = \left[ {{{\hat f}_{1j}}\left( k \right)\quad {{\hat f}_{2j}}\left( k \right)\quad {{\hat g}_{0,j1}}\left( k \right)\quad {{\hat g}_{0,j2}}\left( k \right)} \right.\\ \;\;\;\;\;\;\;{\left. {\begin{array}{*{20}{c}} {{{\hat g}_{0,j3}}\left( k \right)}&{{{\hat g}_{1,j1}}\left( k \right)}&{{{\hat g}_{1,j2}}\left( k \right)}&{{{\hat g}_{1,j3}}\left( k \right)} \end{array}} \right]^{\rm{T}}} \end{array} $

参数估计采用如下最小二乘递推算法：

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{k}}_j}\left( k \right) = \frac{{{\mathit{\boldsymbol{p}}_j}\left( {k - 1} \right){\phi _j}\left( {k - 1} \right)}}{{\lambda + \phi _j^{\rm{T}}\left( {k - 1} \right){\mathit{\boldsymbol{p}}_j}\left( {k - 1} \right){\phi _j}\left( {k - 1} \right)}}\\ {{\mathit{\boldsymbol{\hat \theta }}}_j}\left( k \right) = {{\mathit{\boldsymbol{\hat \theta }}}_j}\left( {k - 1} \right) + \\ \;\;\;\;\;\;\;{\mathit{\boldsymbol{k}}_j}\left( k \right)\left[ {{x_j}\left( k \right) - \phi _j^{\rm{T}}\left( {k - 1} \right){{\mathit{\boldsymbol{\hat \theta }}}_j}\left( {k - 1} \right)} \right]\\ {\mathit{\boldsymbol{p}}_j}\left( k \right) = \frac{1}{\lambda }\left[ {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{k}}_j}\left( k \right)\phi _j^{\rm{T}}\left( {k - 1} \right)} \right]{\mathit{\boldsymbol{p}}_j}\left( {k - 1} \right) \end{array} \right. $

（17）

为了实现速度跟踪，设期望输出为X₁^r(k)，实际输出为X₁(k)，$\hat{\boldsymbol{F}}_{1}(k), \hat{\boldsymbol{F}}_{2}(k), \hat{\boldsymbol{G}}_{0}(k)$和${\mathit{\boldsymbol{\hat G}}_1}(k)$为参数估计值，则可设计如下黄金分割控制控制器：

$ \begin{array}{l} {\mathit{\boldsymbol{u}}_1}\left( k \right) = - {\left[ {{{\mathit{\boldsymbol{\hat G}}}_0}\left( k \right) + \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}} \right]^{ - 1}}\left[ {{l_1}{{\mathit{\boldsymbol{\hat F}}}_1}\left( k \right){{\mathit{\boldsymbol{\tilde Y}}}_1}\left( k \right) + {l_2} \cdot } \right.\\ \;\;\;\;\;\;\;\left. {{{\mathit{\boldsymbol{\hat F}}}_2}\left( k \right){{\mathit{\boldsymbol{\tilde Y}}}_1}\left( {k - 1} \right) + {{\mathit{\boldsymbol{\hat G}}}_1}\left( k \right){\mathit{\boldsymbol{u}}_1}\left( {k - 1} \right)} \right] \end{array} $

（18）

式中：$\tilde{\boldsymbol{Y}}_{1}(k)=\boldsymbol{X}_{1}^{\mathrm{r}}(k)-\boldsymbol{X}_{1}(k)$；l₁=0.382、l₂=0.618分别为黄金分割数，文献[19]证明，当参数选用黄金分割数时，由控制器(18)组成的闭环系统动态性能较好；Λ=diag(λ₁, λ₂, λ₃)为正常值矩阵，避免矩阵求逆过程出现奇异。

2.2 位置跟踪的PD控制器

为了实现Halo轨道的长期维持，除了跟踪速度量以外，还需要实现对位置量的跟踪。上述多变量黄金分割控制控制器设计中仅考虑的速度误差的反馈，因此还需设计控制器引入位置误差反馈。令X₂^r为期望输出，X₂为实际输出，设计如下简单的PD控制器实现位置跟踪：

$ {\mathit{\boldsymbol{u}}_2} = {\mathit{\boldsymbol{u}}_{\rm{p}}} + {\mathit{\boldsymbol{u}}_{\rm{d}}} $	（19）

式中：

$ {\mathit{\boldsymbol{u}}_{\rm{p}}} = {\mathit{\boldsymbol{K}}_{\rm{p}}}{{\mathit{\boldsymbol{\tilde Y}}}_2}\left( k \right) $	（20）

其中：${\mathit{\boldsymbol{\widetilde Y}}_2}(k) = \mathit{\boldsymbol{X}}_2^{\rm{r}} - {\mathit{\boldsymbol{X}}_2};{\mathit{\boldsymbol{K}}_{\rm{p}}} = {\mathop{\rm diag}\nolimits} \left({{k_{{\rm{p}}1}}, {k_{{\rm{p}}2}}, {k_{{\rm{p}}3}}} \right)$为常值矩阵。

为了便于调节，u_d采用逻辑微分，其形式为

$ \left\{ \begin{array}{l} {u_{{\rm{d}}j}} = {c_j}\left[ {{x_{2j}}(k) - {x_{2j}}(k - 1)} \right]\sqrt {\sum\limits_{i = 1}^{{N_j}} {\left\{ {{{\left[ {{{\tilde y}_{2j}}(k - i)} \right]}^2} + {{\left[ {{{\tilde y}_{2j}}(k - i) - {{\tilde y}_{2j}}(k - i - 1)} \right]}^2}} \right\}} } \\ {\mathit{\boldsymbol{u}}_{\rm{d}}} = {\rm diag} \left( {{u_{{\rm{d}}1}},{u_{{\rm{d}}2}},{u_{{\rm{d}}3}}} \right) \end{array} \right. $

（21）

式中：c_j、N_j为常数。

至此，可以得到最终的控制器为

$ \mathit{\boldsymbol{u}} = {\mathit{\boldsymbol{u}}_1} + {\mathit{\boldsymbol{u}}_2} $	（22）

由MIMO特征模型式(15)和黄金分割控制器式(18)组成的闭环系统的稳定性，一直是相关研究学者关注的问题，其中文献[22-23]均给出了闭环系统稳定的充分条件，但该充分条件仍有很大的局限性，因此基于MIMO特征模型的黄金分割控制器的稳定性证明仍有待进一步研究和完善。由于控制器(22)为离散形式，而原系统(14)为连续模型，因此由控制器(22)和原系统组成的闭环系统是一个典型的混杂系统，目前关于混杂系统稳定性证明尚无有效的工具，仍有待进一步研究。基于以上分析，本文中暂不涉及控制器稳定性的证明。

3 仿真结果

为了验证第2节中设计的控制器的有效性，选取了地月L₂点附近的Halo轨道进行了仿真，并与LQR方法下的仿真结果进行了对比。设置轨道z向振幅为A_z=0.016 6，相当于约6 391.5 km；周期为T=3.412 2，相当于14.738 6天；仿真总时长为20个周期，约300天。

3.1 CRTBP模型下的仿真结果

针对地月L₂点附近的Halo轨道维持控制问题，利用第2节中设计的控制器，首先在CRTBP模型下进行了仿真。考虑到实际任务中，航天器进入Halo轨道时，不可避免地会产生入轨误差，导致实际轨道偏离标称轨道，因此，仿真时设置了初始位置误差和速度误差，其大小为

$ \left\{ \begin{array}{l} \Delta {\mathit{\boldsymbol{X}}_0} = {\left[ {\begin{array}{*{20}{c}} {0.0001}&{0.0001}&{0.0001} \end{array}} \right]^{\rm{T}}}\\ \Delta {\mathit{\boldsymbol{V}}_0} = {\left[ {\begin{array}{*{20}{c}} {0.0001}&{0.0001}&{0.0001} \end{array}} \right]^{\rm{T}}} \end{array} \right. $

约相当于国际单位制下的：

$ \left\{ \begin{array}{l} \Delta {\mathit{\boldsymbol{X}}_0} = {\left[ {\begin{array}{*{20}{c}} {38.44\;{\rm{km}}}&{38.44\;{\rm{km}}}&{38.44\;{\rm{km}}} \end{array}} \right]^{\rm{T}}}\\ \Delta {\mathit{\boldsymbol{V}}_0} = {\left[ {\begin{array}{*{20}{c}} {0.{\rm{1}}\;{\rm{m/s}}}&{0.{\rm{1}}\;{\rm{m/s}}}&{0.{\rm{1}}\;{\rm{m/s}}} \end{array}} \right]^{\rm{T}}} \end{array} \right. $

设置系统采样周期为Δt=0.001，仿真结果如图 6所示。

位置误差曲线和速度误差曲线如图 7(a)和7(b)所示。由图 7(a)可以看出，采用黄金分割控制器和PD控制器的情况下，位置和速度误差均在很短时间内收敛到0；而采用LQR的情况下，速度误差虽然仍能很快收敛，但位置误差则表现出明显的震荡。

稳态后采用两种控制器情况下的位置和速度误差均值如表 1所示。由表中数据可以看出，相比于LQR方法，采用本文设计的黄金分割控制器和PD控制器进行轨道维持，位置误差均值和速度误差均值均有显著下降，由此可以说明本文设计的控制器相比于LQR方法具有更高的控制精度。

控制加速度曲线如图 7(c)所示。可以看出，采用黄金分割控制器和PD控制器情况下，控制加速度最大值出现在起始时刻，之后很快趋向于0；而采用LQR方法情况下，控制加速度则一直相对较小。对控制加速度进行时间积分可得到轨道维持过程中所需要的速度增量，如表 2所示。

从表 2中数据可以看出，20个周期的仿真过程中，采用本文设计的控制器，所需要的总速度增量为95.513 0 m/s，采用LQR方法所需要的总速度增量为31.670 8 m/s。由于仿真初值加入了入轨误差，起始阶段为了克服初始误差，需要比较大的控制加速度，同时，由于参数辨识未收敛，也会导致初始控制量较大，采用本文设计的控制器时，第一个周期所需要的速度增量为73.191 7 m/s，相比于LQR方法下的1.858 5 m/s，控制消耗较大。而达到稳态后，采用本文设计的控制器时，每个周期所需要平均速度增量则仅为1.322 4 m/s，相比于第1周期有显著降低；而采用LQR方法，稳态后平均每个周期所需的速度增量为1.569 9 m/s，相比于第1周期也有所下降，但仍高于采用本文设计的控制器时的所需的速度增量。

由此可知，在CRTBP模型下，相比于LQR方法，利用第2节所设计的控制器进行轨道维持，可实现精度较高的控制效果，且进入稳态后，控制消耗较小。另外，提高入轨精度，减小入轨时的位置和速度误差，可以显著降低起始阶段的控制消耗。

3.2 BCRFBM模型下的仿真结果

为验证控制器在太阳引力摄动影响下的有效性，本文基于BCRFBM模型，利用第2节设计的控制器以及LQR方法分别进行了仿真。仿真过程中，仍然考虑了入轨误差的影响，初始误差大小与式(23)相同，采样周期同3.1节，初始时刻太阳相位角取θ₄₀=0。仿真结果如图 8所示。

位置误差曲线、速度误差曲线以及控制加速度曲线如图 9所示。可以看出，采用第2节设计的控制器时，位置和速度误差同CRTBP模型下的仿真结果一样，均能在很短时间内收敛到0；而采用LQR方法时，位置误差和速度误差均呈现明显的震荡，不能收敛。

稳态后采用两种控制器情况下的位置和速度误差均值如表 3所示。由表中数据可以看出，采用第2节设计的黄金分割控制器和PD控制器进行轨道维持时，相比于初始误差，位置误差均值和速度误差均值仍能维持在较低的水平；而采用LQR方法时，位置和速度误差均值相比于初始误差均有较大幅度上升。由此可以再次说明本文设计的控制器在控制精度方面的优势。

从控制加速度曲线图可以看出，与CRTBP模型下的仿真结果类似，采用第2节设计的控制器时，由于存在初始误差，且参数辨识尚未收敛，控制加速度在起始时刻较大，之后很快趋向于0；而采用LQR方法，控制加速度一直比较平稳。

对控制加速度进行时间积分可得到轨道维持过程中所需要的速度增量，如表 4所示。从表中数据可以看出，20个周期的仿真过程中，采用第2节设计的控制器时，所需要的总速度增量为745.024 6 m/s；由于仿真初值加入了入轨误差，起始阶段为了克服初始误差，同时由于参数辨识未收敛，需要比较大的控制加速度，第1周期所需要的速度增量为130.037 4 m/s；达到稳态后，每个周期所需要平均速度增量则仅为33.166 3 m/s，相比于第1周期有显著降低。而采用LQR方法时，所需的总速度增量、第1周期所需的速度增量以及稳态后每个周期所需要的平均速度增量，分别为610.551 6、30.294 2和30.544 7 m/s，相比于使用本文设计的控制器情况下均有所下降。

与CRTBP模型下仿真结果对比可以看出，采用第2节设计的控制器时，在轨道维持过程中引入太阳引力摄动，所需要的速度增量提升了一个数量级；x轴和y轴的位置跟踪误差同样也提升了一个数量级，分析原因，主要在于太阳引力摄动主要作用于航天器x轴和y轴运动，对z轴运动则影响不大，但考虑到Halo轨道本身尺度较大，这样的位置误差仍是可以接受的；z轴的位置跟踪误差以及三轴的速度跟踪误差则与CRTBP下的仿真结果相差不大。同时，与CRTBP模型下的仿真结果类似，提高入轨精度，起始阶段的控制消耗有望进一步降低。而采用LQR方法时，虽然控制消耗有所下降，但其在消除跟踪误差方面表现很差，因此综合来看，在达到稳态后，本文设计的控制器仍具有优势。

从仿真位置和速度误差来看，本文所设计的控制方法是一种跟踪精度较高的策略。在Halo轨道的维持控制中，控制精度与控制消耗互相制约，更高的控制精度通常也意味着更高的控制消耗。此外，根据文献[24]中的研究结论，标称轨道的精确度与控制消耗密切相关。标称轨道越精确，维持控制所需要的控制消耗则越低。本文所采用的基准Halo轨道是在CRTBP模型下修正的结果，对于考虑太阳引力的BCRFBM模型下的维持控制来说，精确度不够。若能采用更精确的标称轨道，维持控制所需的总的速度增量有望进一步降低。

4 结论

本文首先基于Richardson三阶近似解析解、微分修正及打靶法获得了CRTBP模型下的地月L₂点附近的标称Halo轨道，然后设计了基于特征模型理论速度跟踪的黄金分割控制器以及位置跟踪的PD控制器，最后分别在CRTBP模型和BCRFBM下进行了地月L₂点Halo轨道维持控制仿真，并与LQR方法下的仿真结果进行对比，验证了控制器的有效性。通过分析仿真结果，可以得到以下结论：

1) 在CRTBP模型下，仿真结果跟踪精度较高，控制消耗很小。

2) 在BCRFBM模型下，位置跟踪精度有所下降，但相比于Halo轨道的尺度，仍然可以接受；但由于太阳引力作用以及标称轨道精度影响，控制消耗显著增大。

3) 提高入轨精度，可以显著降低起始阶段的控制消耗。