空间站作为长期在轨运行平台，其多舱段的复杂结构需要长期组装才能完成。从1987年质子1号发射到1996年完成组装，“和平号”用了9年时间；从1998年“曙光号”升空到2006完成装配，国际空间站用了近8年的时间。为了保证能量供应，在其组装早期阶段一般采用惯性系定向飞行模式。

目前研究姿态控制/角动量管理(Attitude Control/Momentum Management，ACMM)的成果较多，但大部分是在轨道系下进行的。Wie等^[1]利用重力梯度力矩和陀螺耦合力矩来平衡姿态，假设惯量积为零且姿态偏差为小量，将俯仰轴与滚动/偏航轴解耦，并且引入滤波变量采用线性二次型(Linear Quadratic Regulators，LQR)算法求解控制器，但该算法在姿态角偏差较大时并不适用。Parlos和Sunkel^[2-4]针对大角度偏差的情况采用自适应控制算法，但该算法所设计出的控制器较为复杂，实现起来有一定难度。

Harduvel^[5]采用三轴耦合的模型，引入滤波状态量抑制常值、1倍和2倍于轨道频率的扰动对系统的影响，在轨道系下设计了线性化的ACMM控制器，该控制器采用传统的LQR算法，即通过选取权重矩阵Q和代价矩阵R解出反馈增益矩阵，这也是目前国际空间站采用的ACMM算法。该算法的不足之处在于：所得到的三轴耦合状态空间方程达到了24维，要想选出合适的Q矩阵使得系统有较好的暂态和稳态性能并不是一件容易的事情，且线性化系统不适用于大角度偏差的情况。为了适应系统大角度偏差，朱孟萍等^[6-7]分别采用基于参数辨识的自适应算法和非线性最优控制算法对轨道系下ACMM进行了仿真验证，取得了较好的效果，但其采用θ-D^[8]算法求解状态相关黎卡提方程得到的是次优解。

在极点配置理论方面，针对LQR算法在求解高维状态方程时遇到Q矩阵选取困难的问题，Shieh等^[9-10]提出了带极点配置的LQR算法，该算法能够将极点配置到与虚轴夹角为±45°或±30°的扇形范围内，同时为了保证闭环系统具有一定的鲁棒性，该算法可以在虚轴左侧指定一条平行于虚轴的直线，然后将极点移动到该直线的左侧^[11-12]。这种设计方法的优越性在于：极点配置并不是通过指定每一个极点位置来实现的，而是指定一个扇形区域，将闭环系统配置到该区域内，实现起来较为简单。Sunkel等^[13-14]建立了轨道系下ACMM模型，同时利用该算法求解出控制器，但在建立模型时并没有抑制扰动对姿态的影响，导致其姿态角的周期性波动较大。

对于惯性系下的ACMM，张军等^[15]采用与文献[5]相似的方法验证了惯性系下的角动量管理的可行性，但该算法也存在选取Q矩阵困难的问题，同时没有验证系统的可控性，在俯仰轴方向只要存在常值扰动即有可能导致系统失控，即在俯仰轴方向系统没有鲁棒性。

本文对惯性系下ACMM问题进行分析，考虑到惯性系下气动力矩绕俯仰轴一周基本不积累，故将其作为外部扰动处理，仅采用引力梯度力矩来平衡姿态。由于角动量反馈控制对挠性附件具有较强的鲁棒性^[16]，并且柔性构件的固有频率一般在0.1Hz以上，远大于轨道频率，因此在控制器设计时忽略太阳帆板挠性模态的影响。同时分析了在俯仰轴方向角动量管理的不可行性，将俯仰轴与滚动/偏航轴解耦，将扰动项引入方程以抑制其对姿态的影响。采用带极点配置的LQR算法^{[9-10, 17]}求解出反馈增益矩阵，该算法根据给出的代价矩阵R、稳定裕度以及阻尼角, 通过迭代算法将系统闭环极点配置到指定的区域。

1 惯性系ACMM模型

定义轨道系f_o(O_ox_oy_oz_o)，原点在系统质心，z轴在轨道平面内指向地心，x轴正向指向空间站飞行方向，y轴按右手定则确定；惯性系f_i(O_ix_iy_iz_i)，其是一种中间坐标系，在t₀时刻与轨道系重合；本体系f_b(O_bx_by_bz_b)，姿态角为零时与惯性系重合。

忽略柔性构件，惯性系下空间站刚体部分姿态动力学方程为

式中：_sⁱ为惯性系下空间站本体角动量变化量；T_cⁱ为控制力矩陀螺(CMG)的输出力矩；T_gⁱ为引力梯度力矩；T_dⁱ为其他干扰力矩。

CMG动力学方程为

式中：_cⁱ为CMG角动量变化量。

引力梯度力矩在近圆轨道下为

式中：ω_o为轨道角速度；Rⁱ为惯性系下地心到卫星的单位矢量；Iⁱ为惯性系下本体转动惯量，即

由坐标系定义可知，在t₀时刻轨道系与惯性系重合，则在t时刻轨道系与惯性系之间的旋转角度即为轨道系绕y_o轴转到惯性系的角度υ_io=ω_o(t-t₀)，于是轨道系到惯性系的坐标转换矩阵可表示为

显然，地心到卫星的单位矢量在惯性系下的分量与轨道系下分量关系为

式中：R^o=00-1^T为轨道系下地心到卫星的单位矢量。

结合式(3)~式(5)，惯性系下引力梯度力矩为

采用3-1-2旋转顺序的欧拉角表示本体系到惯性系的旋转角度θ_ib=[φ υ ψ]^T，在小角度偏差时，从本体系到惯性系的转换矩阵近似为

定义本体系下航天器本体的惯量阵为

忽略高阶项，根据2阶张量之间的转换关系I^b=C_i^bIⁱC_bⁱ，得到本体系与惯性系下空间站转动惯量之间的近似关系为

将式(9)代入式(6)得到本体系下转动惯量表示的惯性系下引力梯度力矩表达式为

式中：

显然，引力梯度力矩被分为两部分：姿态相关项和姿态无关项。当本体系坐标轴与惯性系主轴重合时，惯量积为零，此时引力梯度力矩只有周期性的成分，而不会引起CMG角动量积累，因此惯性系下角动量管理就是通过姿态偏置来消除引力梯度力矩及其他扰动引起的常值积累。

惯性系下空间站角动量为

式中：ω_ibⁱ为本体系到惯性系的绝对角速度在惯性系下的分量，其与姿态角速度之间的关系为(_ib=)

对式(11)保留1阶项，得到线性化的姿态运动学方程为

将式(10)代入式(1), 再联合式(2)和式(13), 得到线性化的ACMM模型。将此模型写成状态方程的形式如下：

式中：E为3×3的单位矩阵。

2 ACMM状态方程分析

从引力梯度力矩表达式(10)可以看出，在姿态零偏置时，滚动/俯仰、俯仰/偏航方向的惯量积会导致CMG积累。因此在姿态平衡处，滚动/俯仰、偏航/俯仰轴的主惯量差产生控制力矩补偿滚动/俯仰、偏航/俯仰惯量积以及其他扰动带来的常值力矩。系统达到稳定姿态时，其姿态偏置可以通过引力梯度力矩及其他干扰力矩大致估算出：

式中：φ^*和ψ^*分别为滚动轴和偏航轴平衡姿态角；T_d^x为扰动力矩T_dⁱ在惯性系x轴方向的常值分量；T_d^z为扰动力矩T_dⁱ在惯性系z轴方向的常值分量。

下面对俯仰轴方向是否能进行角动量管理进行理论分析。空间站可以进行角动量管理的必要条件是需要存在额外力矩来平衡扰动。重新观察式(10)可以发现，在滚动/偏航方向上都存在与姿态相关的引力梯度力矩，通过姿态机动就可以调节引力梯度力矩来抵消扰动力矩对系统的影响，从而避免滚动/偏航方向的CMG角动量积累。但在俯仰轴方向上并不存在与姿态相关的引力梯度力矩，说明俯仰轴方向上没有额外力矩来平衡可能存在的常值扰动，即俯仰轴方向扰动力矩产生的角动量只能通过CMG吸收。需要强调的是，这并不是近似线性化产生的结果，在用惯性系下转动惯量表示引力梯度力矩的式(6)中，并没有进行近似处理，但可以看出俯仰轴方向也只存在2倍于轨道频率的扰动力矩。

下面从可控性的角度来分析俯仰轴不能角动量管理对系统的影响。显然系统状态方程式(14)中至少有3个极点“0”，利用PBH判据来判断极点“0”的可控性，可以发现无论其他参数如何，式(14)的可控性指数均小于8，这就说明至少存在1个极点“0”是不可控的。结合本节分析可知，在俯仰轴方向上存在不可控状态量，即姿态角和CMG动量不能同时约束，只要在俯仰轴方向上存在常值扰动，就会导致姿态角发散或者CMG角动量积累，换句话说在俯仰轴方向不具备鲁棒性。在实际工程应用中，这是不能接受的。

考虑到ACMM首先是要能够控制姿态，并保证一定的姿态精度。为了解决上述问题，将俯仰轴与滚动/偏航轴解耦，不再约束俯仰轴方向的CMG角动量，因此去除状态方程中俯仰轴方向的CMG角动量，这样在俯仰轴方向上的姿态控制就不再考虑CMG角动量是否积累，滚动/偏航轴方向则保持不变。同时考虑到大气扰动主要是1倍于轨道频率，而从式(10)可以看出引力梯度力矩主要为2倍于轨道频率，为了让周期性扰动通过CMG吸收，而常值扰动通过姿态偏置来抵消，根据内模原理将系统状态方程扩维，引入滤波状态量：f=[f₀^Tf1₁^Tf1₂^Tf2₁^Tf2₂^T]^T，f₀用于抑制常值扰动，f1₁、f1₂用于抑制1倍于轨道频率的扰动，f2₁、f2₂用于抑制2倍于轨道频率的扰动，把常值、1倍和2倍于轨道频率的扰动引入状态方程中，得到最终的滚动/偏航轴状态方程及控制方程如下：

式中：；E_2×2为2×2的单位矩阵；I_xz=diag{I_x, I_z}；h_c^txz为CMG角动量在滚动/偏航轴方向的分量组成的列向量；T_c^ixz为T_cⁱ在滚动/偏航方向的分量组成的列向量；K为2×16的增益矩阵；x^T为h_cⁱ、θ_ib、H_sⁱ和f在滚动/偏航轴方向的分量所组成的列向量，即x^T=。

相应的，在俯仰轴方向上状态方程及控制器为

式中：T_c^iy为T_cⁱ在俯仰轴方向上的分量；y=为θ_ib、H_sⁱ和f在俯仰轴上的分量组成的列向量；K_y为1×7的反馈增益矩阵。

3 带极点配置的LQR算法 3.1 极点配置算法基本理论

本文采用的极点配置算法是Shieh等^[9-10]提出的，该算法通过给定的阻尼角、稳定裕度以及系统的代价矩阵R，通过迭代自动将极点配置到指定的扇形区域，该扇形区域与虚轴的夹角可以为30°~90°，为实现该算法，现介绍相关定理和引理如下。

考虑线性时不变连续时间系统：

设LQR代价函数为

式中：A和B分别为状态转移矩阵和控制输入矩阵；权重矩阵Q和代价矩阵R分别为R^n×n的半正定矩阵和R^m×m的正定矩阵。

最小化式(21)性能指标的反馈控制律为

式中：K=-R^-1B^TP为m×n反馈增益；P为n×n非对称阵，满足黎卡提方程：

引理1^[9]     记(A, B)为等式(20)描述的系统的可控对，取h≥0表示系统的稳定裕度，则闭环系统A-BR^-1B^TP所有的特征值将位于-h垂线的左侧，其中矩阵P是下列方程的解：

引理2^[17]    记{λ_i^-}n_i=1^-表示矩阵A在虚轴上和纯左半平面的特征子集，{λ_i⁺}n_i=1⁺表示在纯右半平面上的特征子集，则有：

1) 最优闭环系统的特征值集合:

当γ=1且矩阵P满足Q=0时的黎卡提方程式(23)。

2) 记(λ_i^-)n_i=1^-为特征值{λ_i^-}n_i=1^-对应的特征向量，那么黎卡提方程的最大解半正定P满足：Pζ_i^-=0, i=1, 2，…。而对应特征值集合{λ(A-BR^-1B^TP)}保持不变。另外的n⁺个特征值在纯左半平面上，其中γ>0.5。

引理3^[17]    对于给定矩阵A∈R^n×n，定义2个相关矩阵A^-∈R^n×n和A⁺∈R^n×n，并且有

式中：所有参数同引理2。从引理2还可知，闭环系统A-BR^-1B^TP所有特征值λ_i^-、-λ_i⁺都是当Q=0时的黎卡提方程式(23)的解，从而有

式中：tr(·)为矩阵的迹。

引理4^[17]    若(A, B)可控，则(-A², B)、(A³, B)也是可控的。

引理5^[9]     考虑在复平面左半平面内的一个扇形区域：其边界与负实轴夹角为±45°/±30°直线包络而成。假定系统矩阵A为渐近稳定，其所有特征值分布在左半平面内，可分成扇形区域内、外2个特征值子集。定义保角映射，则对应矩阵A扇形区域内的特征子集映射为矩阵左半平面的特征值子集，而扇形区域外的特征值子集映射在的右半平面上。

定理1^[10]     设稳定的系统矩阵A有特征值{λ_i^-}_i=1n^-位于引理5所示的扇形区域内。而{λ_i⁺}_i=1n⁺位于扇形区域外，考虑以下2个黎卡提方程：

则对于闭环系统为

其特征值除{λ_i^-}_i=1n^-外，至少还会有一个实特征值或一对共轭复特征值位于扇形区域内，对于常增益系数γ满足：

式中：。

3.2 空间站ACMM控制器求解流程

从引理1~引理5可以看出，通过迭代算法可以将系统闭环极点配置到±45°或±30°以内。通常闭环极点在±45°左右时系统具有较短的调节时间，超调量也不会很大。为了使控制器求解有更大的灵活性，现做如图 1和图 2所示改进，即先根据引理1将极点配置到所需稳定裕度的左边，再根据所需要的角度Ω将虚轴向右移动h₂，此时得到的极点将可以大致配置在30°～90°内。

具体算法实现简述如下：

1) 对于状态方程；X(0)=0，设置稳定裕度-h₁、阻尼角及代价矩阵R，求解等式：

解出半正定矩阵P₀，则当前闭环系统为A₁=A-BR^-1B^TP₀，此时闭环系统所有极点在-h₁左侧，记i=1。

2) 判断系统所需的阻尼角，若为30°~45°(见图 2)，转步骤3)；若为45°~90°(见图 1)，转步骤7)。

3) 对于给定的阻尼角Ω，将复平面的虚轴向右移动：

得到新的状态矩阵A₁=A₁-E·h₂。

4) 解方程：

解出半正定矩阵，判断0.5tr()是否为0，若是转步骤11)；若不是则继续。

5) 解方程：

解出半正定矩阵P_i，对于新的闭环系统=A_i-γ_iBR^-1B^TP_i，然后解不等式(36)得到常增益系数γ_i：

式中：。

6) i=i+1;转步骤4)。

7) 对于给定的阻尼角Ω，将复平面的虚轴向右移动：

得到新的状态矩阵A₁=A₁-E·h₂。

8) 解方程：

解出半正定矩阵_i，判断0.5tr(BR^-1B^T_i)是否为0，若是转步骤11)；若不是则继续。

9) 解方程：

解出半正定矩阵P_i，对于新的闭环系统=A_i-γBR^-1B^TP_i，常增益系数γ通过式(40)给出：

式中：。

10) i=i+1, 转步骤8)。

11) 算法完成计算，最终系统增益矩阵为

将u=-Kx代入到系统状态方程，则闭环系统所有的极点都会移动到指定的区域。

分别给出系统俯仰轴与滚动/偏航轴的性能指标，然后将状态方程代入到上述算法即可求解出相应的ACMM控制器。

4 仿真验证

本文针对“T”字构型的空间站进行仿真验证，仿真中动力学模型采用带挠性附件的中心刚体模型，执行机构只采用CMG。考虑到CMG角动量测量较为精确，其测量误差不会影响到控制器的稳定性。另外，由于平衡姿态并非本体系与惯性系重合，不失一般性，取初始姿态角和姿态角速度为零，系统干扰力矩如表 1所示。本体系下转动惯量为

空间站的运行轨道为400km的近圆轨道，故轨道角速度近似取ω_o=0.0011rad/s。

从状态方程式(16)和式(18)可以看出，扩维后的方程在滚动/偏航方向上常值扰动通过姿态偏置吸收，避免了CMG角动量的积累，周期性的扰动通过CMG吸收；而在俯仰轴方向常值扰动、1倍和2倍于轨道频率的扰动均由CMG吸收。从式(15)估算出稳态时平衡姿态角为

在俯仰轴方向和滚动/偏航轴方向代价矩阵R均取单位矩阵。为了使系统有一定的鲁棒性，同时又有比较快的响应速度，仿真时在俯仰轴方向稳定裕度取0.2ω_o，阻尼角取45°。在滚动/偏航方向由于需要约束CMG角动量，如果阻尼角取太小，会使得反馈增益阵过大，从而导致暂态时CMG输出力矩较大，因此在滚动/偏航方向上稳定裕度取0.25ω_o，阻尼角取60°，用第3节算法求解出闭环极点如表 2所示，闭环反馈增益矩阵如表 3所示。为了更加直观地显示闭环极点的位置，现在复平面上表示出，如图 3和图 4所示。

从图 3中可以清楚看出，在俯仰轴方向，所有极点均被配置到了离虚轴0.2ω_o~6ω_o的地方，同时闭环极点都在±45°以内，说明在俯仰轴方向会有较短的调节时间和较小的超调量；从图 4可以看出，在滚动/偏航轴方向所有极点均在±60°内，保证了较小的超调量，且都在离虚轴0.25ω_o~4ω_o内，保证了较短的调节时间，也避免了出现暂态CMG输出力矩较大的情况。系统角动量管理仿真结果如图 5~图 8所示。

从图 5可以看出，3个轨道周期后系统基本达到稳定状态，俯仰轴基本没有偏差，周期性的波动幅值也很小；滚动轴最终稳定在-1.4°左右，偏航轴最终稳定在3.6°左右，这与式(42)估计结果相符。由于在状态方程中引入了滤波变量，抑制了1倍、2倍于轨道频率的扰动对姿态的影响，稳定时的姿态角仅有一些幅值较小的高频波动。

结合图 5~图 8可以看到，CMG的角动量周期性波动但基本不积累，抵消掉了周期性的扰动对姿态的影响。同时在滚动/偏航方向引力梯度力矩几乎没有周期性扰动，而只存在常值偏置，正是因为滚动/偏航方向的姿态偏置才抵消掉了滚动/偏航方向的常值扰动以及俯仰/偏航、俯仰/滚动轴之间的惯量积。俯仰方向的引力梯度力矩引起的角动量周期性波动比较大，这主要通过俯仰方向的CMG来吸收。对于其他干扰力矩来说，滚动/偏航方向的少量的常值项积累通过滚动/偏航方向的姿态机动来与重力梯度力矩相抵消，而俯仰方向几乎没有常值积累，但周期性扰动比较大。总体来说控制器的设计达到了目的。

5 结论

本文综合分析了惯性系下空间站利用引力梯度力矩平衡姿态的可行性，经过理论分析及仿真验证表明：

1) 提出了一种较为实用的惯性系下角动量管理控制器求解算法，在设计空间站构型时，为了保证平衡姿态角较小，需要增大俯仰轴与滚动/偏航轴方向转动惯量差，同时应尽量减少惯量积。

2) 基于极点配置的LQR算法具有较好的通用性，特别适合于高维的状态方程，利用计算机求解有较大优势。

3) 仿真结果说明，尽管俯仰轴方向没有额外力矩来平衡姿态，但在俯仰轴方向常值扰动较小CMG角动量积累较慢，并不需要经常卸载。

为了使本文算法更好地满足实际需求，还需要进一步探索存在不确定性参数时的ACMM问题。