航天器控制指令发出后，推力系统需要确定各个推力器的输出以实现满足控制指令需求的力和力矩，这个过程就是指令分配[1-5]。普通任务中变轨与姿态机动等航天任务推力器工作并不频繁，对指令分配策略的实时性要求并不高。但在无拖曳等需要持续推力的任务中，指令分配实时性就格外重要[6]。无拖曳航天器目前主要应用于引力波探测、重力场高精度测量等科学任务。无拖曳技术是指将惯性检测和微推力控制相结合，在航天器上施加持续的推力，用于抵消大气及太阳光压等非保守力的作用，为重力梯度仪等重力测量设备提供一个近于纯引力作用下的飞行环境，满足有效载荷正常工作对平台的技术要求。高精度无拖曳技术，可以极大推进卫星轨道及姿态控制技术、卫星平台技术的发展，提高卫星精密定轨、高分辨率对地观测、提高地震前兆信息观测能力等，具有重大社会和经济效益。针对无拖曳航天器等需要提供连续推力的航天任务，必须采用运算简单，响应速度快，且满足精度要求的推力分配方法[7-10]。
伪逆法是对于给定的伪控制量来确定可行的控制输入，其实现简单，实时性较好，但由于其会产生执行机构不可达的输出解，因此在实际应用时受到限制。有学者利用线性规划等方法对伪逆法进行了二次优化，但其运算复杂，在轨实时计算应用此方法有一定的难度[11]。本文基于航天器配置结构，提出了从解决输出力为负值与超出门限值2个方面对伪逆法进行二次优化的方法。
1 伪逆法求取推力初值首先不加证明给出伪逆法求取初值的方法。伪逆法可描述为：
|
式中，v表示控制指令需求的力或力矩，D为推力配置效率矩阵，T为要求解的推力器输出。最优解为：
|
式中，D+=DT(DDT)-1为推力配置效率矩阵D的伪逆阵。
由于无拖曳航天器需要同时跟踪无拖曳轨道与姿态，因此令
|
式中，A、B分别为力矩效率矩阵和力效率矩阵。则推进系统的实际输出为
|
伪逆法的计算会出现2个问题:①存在推力为负的情况; ②超出执行机构约束饱和限制。这两种情况下推力输出均无法实现，因此，对伪逆法所求结果必须进行二次优化才能作为推力输出。
2 伪逆法输出为负值情况优化用伪逆法计算得到初始推力输出后，首先，判断推力正负，如果为负，则寻找与此推力器成对的推力器产生与之相反的推力对其进行替代，因此，在进行推力分配前，需要建立成对的推力器之间的关联关系，这是二次优化方法的关键所在。
本文中对伪逆法的二次优化是基于航天器推力配置结构开展，这里参考GOCE卫星的配置提出一种配置结构，假设无拖曳航天器配置为2个离子推力器和8个冷气推力器构成，配置方案如图 1所示。
|
| 图 1 推力器配置方位示意图 |
针对图 1所示的推力配置，如果利用伪逆法计算得到推力器1输出力为负号，则需要用其他的推力器来对其输出进行替代。设推力器1输出力为-T，则其会产生推力器4指向方向的力和X-方向的力矩，因此可以用推力器4输出T来代替推力器1输出-T，但是推力器4产生的是X+方向的力矩，因此，需要推力器2和推力器3共同完成消除推力器4产生的力矩及代替推力器1产生的力矩，推力器2和推力器3输出均为T。
推力器之间的关联关系如下：
推力器1输出-T⇔推力器4输出T+推力器2输出T+推力器3输出T，简写为
1(-T)⇔4(T)+2(T)+3(T)
同理，
2(-T)⇔3(T)+1(T)+4(T)
3(-T)⇔2(T)+1(T)+4(T)
4(-T)⇔1(T)+2(T)+3(T)
5(-T)⇔8(T)+6(T)+7(T)
6(-T)⇔7(T)+5(T)+8(T)
7(-T)⇔6(T)+5(T)+8(T)
8(-T)⇔5(T)+6(T)+7(T)
如果伪逆法计算的推力输出为负，则用与其相关的推力器对其进行替换，可以达到同样的控制效果，并且所有推力输出均为正值。
假设推力输出用来跟踪气动干扰力和气动力矩、重力梯度力矩及地磁力矩等干扰力矩。根据各个推力器之间的对应关系，对伪逆法计算得到的负值输出进行处理，跟踪效果如图 2所示，可以看出在不考虑推力器输出门限的情况下，通过推力器间的相互替代，推力系统可以实现输出为正情况下高精度的干扰跟踪。
|
| 图 2 干扰跟踪 |
处理后冷气推力器输出与处理前对比如图 3所示。由于离子推力器输出均为正值，这里主要处理冷气推力器组，实线表示处理后的推力输出，虚线表示用伪逆法计算得到的推力输出。可以看出，在二次优化前，伪逆法计算的初值存在输出为负的现象，在二次优化后，所有推力器输出均为正值。
|
| 图 3 冷气推力器输出 |
伪逆法各个冷气推力器消耗分布如图 4所示，可以看出，推力器3和推力器4消耗较多，主要是用来抵消离子推力器9产生的纵向推力。
|
| 图 4 伪逆法各个冷气推力器消耗分布 |
冷气推力器输出范围为0.1~0.7 mN，根据图 3冷气推力器的输出可知，推力均未超出输出上限范围，但下限有小于0.1 mN的情况。本节提出对推力器进行限幅的策略，针对推力器输出推力阈值，对推力输出进行限幅。当推力器需求大于0.7 mN时，输出为0.7 mN。当推力器需求小于0.1 mN时，设置门限值S，若推力输出T满足0 < T < S，则输出为零，若满足S < =T < 0.1，则输出为0.1。S根据输出力和力矩跟踪曲线的方差变化来确定。表 1和表 2列举了S取0.01~0.1时输出力和力矩跟踪曲线的方差变化。
| S | Fx | Fy | Fz |
| 0.01 | 0 | 0.002 | 0.002 |
| 0.02 | 0 | 0.001 | 0.002 |
| 0.03 | 0 | 8.864×10-4 | 0.002 |
| 0.04 | 0 | 6.003×10-4 | 0.002 |
| 0.05 | 0 | 4.504×10-4 | 0.002 |
| 0.06 | 0 | 6.106×10-4 | 0.002 |
| 0.07 | 0 | 6.106 ×10-4 | 0.002 |
| 0.08 | 0 | 0.001 | 0.003 |
| 0.09 | 0 | 0.002 | 0.003 |
| 0.10 | 0 | 0.002 | 0.003 |
| S | Mx | My | Mz |
| 0.01 | 2.024×10-4 | 0.020 | 0.009 |
| 0.02 | 1.580×10-4 | 0.017 | 0.008 |
| 0.03 | 1.097×10-4 | 0.012 | 0.006 |
| 0.04 | 1.070×10-4 | 0.011 | 0.005 |
| 0.05 | 1.213×10-4 | 0.011 | 0.004 |
| 0.06 | 1.529×10-4 | 0.014 | 0.004 |
| 0.07 | 1.529×10-4 | 0.014 | 0.004 |
| 0.08 | 1.556×10-4 | 0.017 | 0.010 |
| 0.09 | 1.805×10-4 | 0.019 | 0.015 |
| 0.10 | 1.805×10-4 | 0.019 | 0.015 |
当S在0.01~0.1之间取值时，Fx、Fy、Fz、Mx、My、Mz的方差变化如图 5所示。
|
| 图 5 跟踪方差 |
可以看出，Fx跟踪误差均为0，除了Fy与Mz方差在S为0.05时取得最小值，其余曲线均在S=0.04时取得极小值，并且当S=0.04时，Fy与Mz方差也处于较小水平。因此，选取门限值为S=0.04。跟踪效果如图 6所示。由图可以看出，推力输出能够跟踪干扰力和干扰力矩，在Y方向跟踪误差较大。
|
| 图 6 推力限幅优化跟踪效果 |
冷气推力工作情况如图 7所示，由图 7可以看出冷气推力器工作在0.1~0.7 mN范围内，输出均为正值。
|
| 图 7 冷气推力工作情况 |
本文利用伪逆法实现无拖曳航天器推力分配，满足运算简单，响应速度快的任务需求。针对伪逆法的2个问题，一是存在推力为负的情况，二是超出执行机构约束饱和限制，从两方面进行二次优化。输出为负的问题通过建立推力器之间的关联关系替代来解决；执行机构饱和限制问题通过确定上限和下限值限幅解决。仿真表明，二次优化的伪逆法可以解决无拖曳航天器推力分配问题，此方法也可应用于其他需要连续推力的航天器推力分配任务。
| [1] |
岳晓奎, 段辰璐. 推力器布设对控制误差的影响分析[J]. 宇航学报, 2011, 32(7): 1522-1526.
Yue Xiaokui, Duan Chenlu. Analysis of Effect of Thruster Configuration on Control Error[J]. Journal of Astronautics, 2011, 32(7): 1522-1526. (in Chinese) |
| [2] | Servidia P A, Pena R S. Spacecraft Thruster Control Allocation Problems[J]. IEEE Trans on Automatic Control, 2005, 50(2): 245-249. DOI:10.1109/TAC.2004.841923 |
| [3] | Bodson M. Evaluation of Optimization Methods for Control Allocation[J]. Journal of Guidance, Control, and Dynamics, 2002, 25(4): 703-711. DOI:10.2514/2.4937 |
| [4] | Jeb S Orr. High-Efficiency Thrust Vector Control Allocation[J]. Journal of Guidance, Control and Dynamics, 2014, 37(2): 374-382. DOI:10.2514/1.61644 |
| [5] | Johansen T A, Fossen T I, Berge S P. Constrained Nonlinear Control Allocation with Singularity Avoidance Using Sequential Quadratic Programming[J]. IEEE Trans on Control Systems Technology, 2004, 12(1): 211-216. DOI:10.1109/TCST.2003.821952 |
| [6] |
吉莉, 刘昆, 项军华. 内编队重力场测量卫星全推力姿轨一体化控制研究[J]. 中国科学:技术科学, 2012, 42(2): 220-229.
Ji Li, Liu Kun, Xiang Junhua. On All-Propulsion Design of Integrated Orbit and Attitude Control for Inner-Formation Gravity Field Measurement Satellite[J]. Sci China:Tech Sci, 2012, 42(2): 220-229. (in Chinese) |
| [7] | Mark van der Meijdea, Pailb R, Binghamc R, et al. GOCE Data, Models, and Applications:a Review[J]. International Journal of Applied Earth Observation and Geoinformation, 2015, 35: 4-15. DOI:10.1016/j.jag.2013.10.001 |
| [8] | García M A Matatoros. GOCE Flight Dynamics Operations from an Orbital Perspective[J]. Journal of Aerospace Engineering, Sciences and Applications, 2012, 4(2): 93-106. DOI:10.7446/jaesa |
| [9] | Lange B. The Control and Use of Drag-Free Satellites[D]. Stanford, Stanford University, 1964 |
| [10] | Canuto E. Drag-Free and Attitude Control for the GOCE Satellite[J]. Automatica, 2008, 44: 1766-1780. DOI:10.1016/j.automatica.2007.11.023 |
| [11] |
刘冰. 组合航天器控制分配方法研究[D]. 哈尔滨: 哈尔滨工业大学, 2008 Liu Bing. The Research on Control Allocation for Cmbined Flight Control[D]. Harbin, Harbin Institute of Technology, 2008(in Chinese) http://d.wanfangdata.com.cn/Thesis/D254226 |
