空间引力波探测卫星微推进器组布局鲁棒优化设计

王继河，于振宁

（1.中山大学 航空航天学院，深圳 518107；2.中国航发控制系统研究所，无锡 214063）

引 言

从20世纪80年代开始，国内外陆续提出了若干项空间引力波探测计划[1]。按轨道类型划分，主要可分为如“激光干涉测量空间天线”（Laser Interferometer Space Antenna，LISA）[2]、“太极计划”[3]等日心轨道的引力波探测计划，以及如OMEGA（Orbiting Medium Explorer for Gravitational Astronomy）[4]、gLisa（geosynchronous Laser Interferometer Space Antenna）/GEOGRAWI（GEOstationary GRAvitational Wave Interferometer）[5]和“天琴计划”[6]等地心轨道的探测计划。由于空间引力波探测任务卫星为保障其中的检验质量在特定方向沿测地线自由运行，故需对空间引力波探测卫星开展无拖曳控制，微推进器组布局配置好坏，将直接影响无拖曳控制性能，本文展开空间引力波探测卫星微推进器组布局优化设计研究。

推进器的配置一般指航天器上配备的所有推进器的数量、推力大小、安装位置和指向等[7]。推进器作为执行机构，其配置方案的好坏直接决定了航天器轨道和姿态控制的精度和燃料消耗等性能指标，所以推进器布局配置的可行性分析以及优化设计对卫星任务的顺利完成至关重要[8]。

目前对于推进器配置布局的研究大多集中在针对已知的配置条件判断其能否满足任意方向的控制命令[9-11]，以及关于冗余度的分析[12-13]。Crowford[9]给出了一些关于推进器配置冗余度的结论。张兵[10]基于Crowford的研究简化了配置矩阵的判别方法，将配置矩阵的可行性判断问题转化为了判断配置矩阵的非奇异变换阵是否可负的问题。Sanchez[11]给出了对于确定配置矩阵控制能力的分析方法，以及对于推进器故障情况的冗余度分析。Milad等[13]研究了考虑燃料消耗、姿态控制误差、脉宽等因素的推进器配置问题。王敏等[14]给出了考虑燃料消耗的推进器配置优化的参考依据。Fauré等[15]在考虑航天器故障诊断和故障容错控制约束下，研究了航天器推力器组最优布局配置问题。Ghobadi等[16]采用遗传优化算法提出了航天器推力器冗余度约束下的最小推力器组最优配置方法。

上述对推进器布局配置的研究都没有考虑到推进器安装的位置及方向的约束，由于空间引力波探测卫星的结构特殊性，其推进器通常都只能安装在太阳能电池板的背面，无法产生电池板背面方向上的推力，在进行推进器布局可行性分析时要考虑到该种约束条件。而且目前对推进器布局配置的研究并没有考虑配置矩阵本身的鲁棒性，卫星的推进器在实际工作时，实际产生的推力方向会与喷嘴方向有一定的误差，同时推进器喷嘴安装方向也会存在误差，故推进器配置方案对误差的鲁棒性也是关键的优化目标。

综上，为得到空间引力波探测卫星微推进器布局的鲁棒优化设计结果，本文首先考虑空间引力波探测卫星推进器布局配置的可行性，然后以推进器配置的鲁棒性为目标，对空间引力波探测卫星微推进器的布局进行鲁棒优化设计，最后，通过数值仿真验证所得到的微推进器组布局方案的鲁棒性。

1 问题描述

地心高轨空间引力波探测卫星受限于卫星的结构，其微推进器组通常只能安装在太阳电池板的背面，布局通常如图1所示，这里暂定共6个微推进器，每2个为1组均匀地分布在X-Y平面上，每组间相差120°，且每组两个推进器的夹角是相同的，定义为角度φ；每个微推进器的单位向量与X-Y平面的夹角也是相同的，定义为角度θ。

图1 微推进器组布局示意图Fig.1 Micro-thruster group layout schematic

在这样配置下，需要对微推进器组的可行性进行分析，保证在该种微推进器组配置下能覆盖用于卫星无拖曳控制的控制力和用于卫星姿态控制的控制力矩输出要求。

地心高轨空间引力波探测卫星电池板通常固定安装在卫星一侧，3颗相同的卫星在轨道上形成一个等边三角形构形，3颗卫星的相位角相差120°，其中一个卫星的轨道参数如表1所示。

表1 卫星轨道根数[17]Table 1 Orbital elements of satellite

图2为地心高轨空间引力波探测卫星的轨道示意图，可以看出，卫星电池板光照条件随地球绕太阳公转而不断变化，故地心高轨空间引力波探测卫星在1年的任务周期内以两种模式运行，分别为科学探测模式和非科学模式。当太阳矢量与轨道面夹角大于δ时，地心高轨空间引力波探测卫星处在星间激光链路指向对准保持状态，卫星电池板法向与轨道面垂直，卫星开展空间引力波探测任务，该段定义为科学探测模式，如图2中的A、C段所示；当太阳矢量与轨道面夹角小于δ时，卫星如果继续处在科学模式，其光照条件会变差进而不能满足空间引力波探测卫星供电需求，故为了满足该阶段的卫星供电需求，空间引力波探测卫星需进入非科学模式，即对日跟踪模式，控制卫星体贴太阳电池板的一面跟踪太阳矢量，具体如图2中的B、D段所示。

图2 地心高轨空间引力波探测卫星的轨道示意图Fig.2 Schematic diagram of orbit of geocentric gravitational wave detection satellite

空间引力波探测卫星微推进器配置角度θ和φ的设计要满足科学探测和非科学探测两种模式的任务需求，且考虑到推进器的安装角度以及推力方向可能存在误差，故本文对推进器布局的设计在满足任务可行性的同时，还对其配置方案进行鲁棒优化设计。

2 空间引力波探测卫星微推进器组布局可行性分析

考虑图1中所示的空间引力波探测卫星的推进器布局情况，针对空间引力波探测卫星的科学模式与非科学模式进行微推进器布局配置角度θ和φ的可行范围分析。

空间引力波探测卫星处在科学模式，微推进器的作用主要是用来抵消太阳光压力来保证引力波探测卫星无拖曳控制顺利进行。通过对运行在如表1所示的轨道上卫星所受的太阳光压力进行分析，可得如图3所示的，在2030年3个月内卫星运行在科学模式抵消太阳光压的目标控制力的变化曲线。

图3 空间引力波探测卫星科学模式下抵消太阳光压的目标力的变化曲线Fig.3 Schematic diagram of orbit of geocentric gravitational wave detection mission

从图3可知，在科学模式下，用于抵消太阳光压摄动力的微推力分量主要为卫星本体系的Z轴正向分量（卫星体贴太阳电池阵一面的法向），而在X轴和Y轴分量呈现周期振荡。微推进器组的布局配置角度要保证微推进器组能够产生覆盖太阳光压力的期望微推控制力的方向和幅值，即微推进器组能产生的推力范围要能够覆盖这些期望的控制力。

对于每一次指令力或者力矩的控制分配，采用线性规划的方法，优化目标函数为6个微推进器产生的力的大小之和最小，可通过Matlab中的linprog线性规划求解器求解。推力分配线性规划问题可描述为

式（1）是对指令力的控制分配问题描述，式（2）是对指令力矩的控制分配问题描述。这里的控制方式为三自由度的控制。式（1）和式（2）中，u为指令力矩，T为指令力，F= [F1，F2，F3，F4，F5，F6]T为6个微推进器各自的推力组成的推力向量；每个微推进器能产生的最大推力为100 μN；矩阵A∈R3×6，为推力的配置矩阵；矩阵B∈R3×6，为力矩的配置矩阵。矩阵A和B的计算表达式为

其中：d为每个推进器在本体坐标系中的位置矢量；e为每个推进器在本体坐标系中产生推力的方向矢量，取图1中所示的内圆半径为1 m。

将图3中的抵消光压的所有期望力作为输入，代入式（1）中，通过线性规划求解，以6个微推进器产生力的幅值之和最小为优化目标，每个推进器推力范围限制在0～100 μN，计算每次控制指令输入得到控制力分配结果。

定义微推进器布局的角度θ的范围为0°～90°，φ的范围为0°～180°。这两个角度在取值范围内排列组合，可得到不同的推进器配置矩阵，在每个配置矩阵下，都可通过线性规划得到控制分配后的每个推进器的推力。空间引力波探测卫星处在科学模式下，对整个取值空间内的推进器配置角度进行遍历，可得到在每个推进器配置角度下的单个微推进器产生的最大推力。在实际空间引力波探测任务中，微推进器受到各种噪声的影响，故要给推进器的推力留有裕度，尽量要让控制分配得到的每个推进器的推力远离最大推力[11]。也就是说，所有微推进器产生过的最大推力要尽量小。

图4为在每个推进器配置角度下，空间引力波探测卫星在科学模式期间内单个微推进器产生过的最大推力图。从图4可知，满足科学模式任务控制需求的推进器配置角度θ的最大值是41°，角度φ的可行范围广，且在不同的配置下，微推进器产生过的最大推力随着角度θ的增大而减小，并且在小角度内的变化率很大，而在30°到41°之间的变化率小且是平坦的；同时最大推力都在10 μN左右，能够给推力最大幅值留有足够的裕度，所以对于配置角度θ来说，应该取为30°～41°。具体取多少最佳，还需考虑推力器配置角度的鲁棒性问题，会在后面给出鲁棒性最好的配置角度取值。

图4 科学模式下抵消光压力时单个推进器产生过的最大推力随配置角度的变化情况Fig.4 Maximum force generated by a single thruster to compensate solar radiation pressure force with satellite in scientific mode

为满足空间引力波探测卫星非科学模式期间的姿态控制力矩需求，将卫星对日跟踪的所有指令力矩作为输入，代入式（2）中，通过线性规划求解器，以6个微推进器产生的力的幅值之和最小为目标，每个推进器推力范围限制在0～100 μN，计算每次指令力矩输入得到控制分配结果。针对空间引力波探测卫星非科学模式，对整个取值空间内的推进器配置角度进行遍历，可得到在每个推进器配置角度下的单个微推进器产生的最大推力如图5所示。

图5 非科学模式下卫星对日跟踪时单个推进器产生的最大推力随配置角度的变化情况Fig.5 Maximum force generated by a single thruster to meet control torque required for sun pointing mode

由图5可知，能满足非科学模式任务需求的配置角度θ的取值范围很广，由于选择安装角度θ和φ应该在同时满足科学和非科学模式任务需求的重合区域，故可依据满足科学模式控制需求的推力器配置角度范围来确定最终的推力器配置角度范围。

已经确定了角度θ的取值范围应该在30°～41°，而对于角度φ来说，其对推进器产生的最大推力的变化影响相对角度θ来说非常小。但是从图4可以发现，随着角度θ的增加，角度φ的取值范围逐渐缩小，由图4可知，角度φ的取值为60°时，微推进器组能产生的力的范围最广。

3 微推进器组布局的鲁棒优化设计

当推进器的安装角度存在误差，或者实际推力方向与推进器喷嘴方向有偏差时，会导致推进器配置矩阵有一定偏差，此时按照控制指令求解得到的推力与指令要求的推力存在误差。对于空间引力波探测卫星，其在科学模式的无拖曳控制精度要求非常高，为增强推进器配置对误差的鲁棒性，这里在满足任务可行性的基础上对微推进器布局配置角度进行了鲁棒优化设计。

3.1 微推进器控制分配的误差估计

引入一个数学定义，即矩阵的条件数，评价推进器布局的鲁棒性。对于矩阵A来说，其条件数cond(A)计算公式为

其中：σmax是矩阵A的最大奇异值；σmin是矩阵A的最小奇异值。

对于方程组AF=T，在数值计算中，矩阵A的条件数越大，解F对于A和T的误差的敏感性越强，最终计算得到的解F的误差越大，精度越低。所以，矩阵的条件数是评价微推进器安装矩阵鲁棒性的适宜指标。通过矩阵的条件数能够计算出误差的范围，从而衡量一个矩阵对误差的敏感性。条件数对线性方程组求解的误差估计参考由以下定理给出[18]。

对于线性方程组AF=T，有误差δA、δF和δT满足(A+δA)(F+δF) =T+δT，则有

当只考虑δT时，可以计算得到解F相对误差的下界以及上界，如式（7）所示

从式（6）和式（7）可以看出，当A和T的相对误差固定时，随着矩阵的条件数cond(A)的增大，解F的相对误差的最大值会随之增大，所以解的稳定性就会降低，导致鲁棒性差。当条件数cond(A)越小时，解F的相对误差的最小值越大，最大值越小。这说明条件数越小，则解的误差会越小，且越可控。从式（5）可以得到，矩阵的条件数最小值为1，当条件数cond(A)取到最小值1时，‖δF‖/‖F‖=‖δT‖/‖T‖，此时解F的误差是完全可控的。

由式（6）和式（7）可知，配置矩阵的条件数越大，其进行线性方程组求解时数值计算的误差越大。所以在优化微推进器组的配置角度θ和φ时，不仅要考虑推进器组产生推力的覆盖范围和推力幅值裕度，还要充分考虑配置矩阵的鲁棒性。

3.2 微推进器布局配置的鲁棒优化仿真分析

3.1节中已经论证了推进器配置矩阵的鲁棒性可通过矩阵的条件数来评价，本节通过式（5）计算式（3）矩阵的条件数随安装角度θ和φ的变化情况，结果如图6所示。

图6 推进器配置矩阵的条件数随角度θ和φ变化情况Fig.6 Condition number of thruster configuration matrix varies with angle θ and φ changes

由于在第2节已经计算得到推进器的配置角度θ的可行范围在41°以内，所以本节只计算了θ在0°～41°矩阵条件数的变化情况。从图6可以看出，微推进器安装矩阵的条件数随角度θ的变化很敏感，而对角度φ的变化灵敏度很低，几乎不会受到角度φ的影响。在角度θ接近0°时，微推进器安装矩阵的条件数变化率很大，而在20°之后，条件数的变化率就会变得平坦，如图6（b）所示。微推进器组的布局在角度θ取35°时，配置矩阵的条件数最小，也就是说，推进器的分配误差对于配置矩阵误差的鲁棒性最强。

综上所述，同时考虑微推进器配置可行性、其可产生的推力覆盖范围、推进器推力幅值的裕度以及配置鲁棒性，微推进器配置的安装角度θ和φ分别确定为35°和60°。

接下来本节针对不同的推进器配置角度，验证分析配置矩阵的条件数对分配误差的影响。仿真场景为：空间引力波探测卫星处在科学模式下，卫星将需施加的控制指令力（用以抵消的太阳光压摄动力）分配给6个推进器。仿真参数为：微推配置角度φ为60°，配置角度θ则取不同的值时，观察其角度θ对分配误差的影响。抵消光压的指令力如图3所示，将3个月的目标力通过采样取为512组数据作为目标输入力。仿真参数如表2所示。

表2 微推进器布局的仿真参数设置Table 2 Parameter setting of micro thruster layout

在不同角度θ的配置下，设置误差δθ为+0.5°，如图7所示，为配置矩阵误差对最终输出力的分配误差‖δT‖的影响。

图7 不同的角度θ下分配误差的影响Fig.7 Distribution error effect with respect to different angles of θ

从图7可以看出，随着角度θ的变化和条件数的降低，对最终输出力的分配误差有很大的影响。θ取35°比θ取5°时对于分配误差的影响会相差一个量级的大小，由此充分证明了3.1节中条件数对于推进器配置鲁棒性的影响，以及选择推进器配置角度θ为35°的优势。

接下来，要验证推进器配置角度φ对误差的影响。仿真场景为：选择相同的推进器布局配置，分别对配置角度θ和φ设置相同的误差大小，比较其造成的最终输出力的分配误差。仿真场景和输入力的设置和上面相同。仿真参数如表3所示。

表3 微推进器布局的仿真参数设置Table 3 Simulation parameter setting of micro thruster layout

取两种配置角度，分别将配置角度θ和φ设置误差δθ和δφ，比较角度θ和φ对最终分配误差的影响，如图8所示。

图8 θ和φ在相同角度误差下的误差灵敏度对比Fig.8 Comparison of error sensitivity of θ and φ with the same angle error

从图8可以看出，角度φ对配置矩阵的误差灵敏度更低，所以对于推进器布局配置矩阵，提高鲁棒性、降低误差灵敏度是以角度θ为主。

综上，微推进器布局的配置角度θ= 35°，φ=60°时配置矩阵鲁棒性最强。

4 结 论

本文针对空间引力波探测卫星微推进器组的特殊布局，优化设计了可满足任务需求同时鲁棒性强的微推进器组布局配置角度。为保证微推进器对指令力或力矩的实现，对微推进器组的布局配置进行了可行性分析，给出了可满足任务要求的配置角度区间，并以矩阵的条件数评价布局的鲁棒性，计算得到了鲁棒性最强的布局配置角度，然后考虑推进器布局配置的鲁棒性，以配置矩阵的条件数作为评价标准，优化设计出对布局配置误差鲁棒性最强的推进器布局结果为θ取35°，φ取60°；仿真验证了文中所提出的空间引力波探测卫星微推进器组布局的鲁棒性。