A Design Method for Redundant Microthruster Layout With Allocation Error

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

As a high-precision measurement platform, the drag-free satellite cannot install moving parts to meet the hardware design requirements for being a rigid body. Therefore, the microthruster is generally used to be the main executive mechanism of drag-free satellites in space. Thrust and torque produced by microthrusters should compensate for the nonconservative force and torque accurately so as to keep relatively static between the satellite and the internal test mass without contraction. The failure or damage of the thruster could occur in the long-term operation, which may lead to the instability of the control system and even the failure of the mission.

In order to avoid the above problem, one method is the master/replica actuator configuration. In 2015, the European Space Agency launched a satellite Laser Interferometer Space Antenna (LISA) Pathfinder (LPF) [1, 2]. LPF is equipped with the main and backup cold gas microthrusters as the execution system to operate safely and reliably in the halo orbit of the Lagrange L1 point in the Sun–Earth system. And the key technologies of drag-free control in space gravitational wave detection have been verified. Apart from the main backup method, the failure of the thruster can also be figured out by the redundant thruster configuration. The redundant thruster layout means that the thrusters used in the satellite are more than the degrees of freedom required for the control mission to make the system have a certain fault tolerance. In this research, Crawford [3] determined the relationship between the number of thrusters in the redundant configuration and the degrees of freedom of the satellite. For the LISA program [4, 5] in the future, the redundant microthruster configuration is proposed to complete the mission of drag-free control. In addition, Taiji-1 [6], a drag-free satellite launched by the Chinese Academy of Sciences in 2019, installed four radiofrequency ion thrusters and four Hall electric thrusters, which were used in single-degree-of-freedom drag-free control.

Most high-performance satellites could be equipped with redundant actuators in the control system to complete the control task. Xia et al. [7] extracted the design parameters of the GEO satellite thruster layout and optimized the process through programming to reduce the design iteration cycle and improve design efficiency. Gazzino et al. [8] introduced a control scheme to realize the satellite equipped with electric thrusters to minimize fuel consumption while maintaining position using four electric thrusters. Li et al. [9] designed a fault detection algorithm and developed two efficient station-keeping strategies in full and failure modes. The eight thrusters are used in normal mode, and the specified four thrusters are used in fault mode. However, apart from station-keeping, attitude control is also necessary for the system to be stable. Weiss, Kalabić, and Cairano [10] proposed a model predictive control strategy, which used six electrically powered thrusters and three axially symmetric reaction flywheels as actuators for orbit control and attitude control, respectively. In addition, Wei et al. [11] proposed a novel electric thruster configuration design combined with the 4-DOF robotic arm, which can provide flexible thrust in space to reduce propellant. Besides, to solve the problem of thruster failure, Faure et al. [12] studied the optimal layout of spacecraft (SC) thrusters under the constraints of fault diagnosability and fault recoverability.

Thrusters are more suitable for satellite models close to rigid body dynamics than robot arms and reaction flywheels. And the redundant thruster layout entails fewer thrusters at a lower cost than the master/replica actuator configuration because they could change the center of mass of the satellite. However, it should be noted that the redundant thruster layout could lead to unique thrust instructions. Many scholars proposed diverse solutions, such as the direct allocation method [13], daisy chain control allocation [14], and generalized inverse allocation methods [15]. Johansen and Fossen [16] summarized common methods under the nonlinear and linear models of actuators under different application backgrounds. Although varieties of algorithms can be used to obtain the proper instructions, the SC equipped with redundant thrusters still produce allocation errors in the control process, which will affect the control stability of the whole system [17]. Buffington and Enns [14] conducted a study on the stability of control systems with control allocation. To reduce the allocation error, the method based on mathematical programming [16] sets the minimum error as the objective function of the control allocation. Besides, Yang et al. [18] presented an optimization strategy for high-precision attitude control using the solid propellant microthruster array. Doman and Oppenheimer [19] proposed a closed-loop feedback control method to improve the stability of the system with control allocation, considering the effects of the efficiency matrix error caused by the installation deviation. Besides, the microthruster applied to a drag-free satellite needs an extremely strict thrust resolution. To decrease the actuator output deviation, Ma et al. [20] developed the radiofrequency ion microthruster μRIT-1 by optimizing the key structural components, which could make thrust output continuously adjustable. Moreover, Yang [21] developed the cold gas microthruster with an extremely strict thrust resolution to decrease output deviation.

Many researches only consider the relation between the thruster number and the orbit control or the attitude control. But the orbit control and the attitude control must be performed in the drag-free control process. And the design of the thruster layout for the drag-free control needs to take many factors into consideration, including the thruster number, position, and orientation. Moreover, a high microgravity level of drag-free control leads to high requirements for designing the thruster layout. But optimization algorithms and hardware optimization are unsuitable for the design of thruster layout to decrease errors. It is complex and difficult to design a suitable layout to optimize the thruster layout, considering allocation errors. Therefore, it is necessary for high-precision drag-free control to do research on designing the redundant thruster layout.

To realize the high precision attitude and orbit control, the paper proposes a design method for redundant thruster layout to reduce the distribution error generated by various influence factors. First, the main components of allocation error are derived from theoretical analysis. And the thruster configuration factor (TCF) is determined to represent the position and angle of the thrusters involved in the allocation error. In addition, the design criterion of minimum required thrusters is derived to be applied to the satellite model with the orbit control and the attitude control. The criterion is put forward to balance the thruster number and fault tolerance. On this basis, an optimization model is developed to reduce allocation error and improve the availability of layout design. Finally, the simulation results of the gravitational wave detection satellite show that the design of redundant thrusters can reduce the allocation error and improve the fault tolerance of the control system. It demonstrates the feasibility of the proposed method. This work is expected to provide a reference for the redundant thruster configuration design of the drag-free satellite or other satellites in the future.

The paper is organized as follows: Section 2 shows mathematical descriptions of thrust allocation. Section 3 presents the deduction of the optimization principle and method. Section 4 provides the application of the layout design and simulation strategy. Section 5 discusses the results. Section 6 presents the conclusions.

2. Mathematical Model of Thrust Allocation

The disturbing force exerted on the drag-free satellite with redundant thrusters is measured by the high-precision inertial sensor. And the measured disturbance is taken as input into the control algorithm to obtain the expected control quantity $\mathit{F}$. Then, the thrust command of each thruster is calculated by the control allocation algorithm. The calculation formula is as follows: $$\begin{array}{}\text{(1)}& {\mathit{F}}_c=\mathit{A}\mathit{T}\end{array}$$

The expected control quantity ${\mathit{F}}_c$ is the generalized active control force. $\mathit{T}$ is the thrust column vector. $\mathit{T}={T_1\cdots T_i\cdots T_n}^{\text{T}}$, where $T_i$ represents the thrust generated by the $i$-th thruster. $\mathit{A}$ is the control efficiency matrix, which is determined by the direction, position, and number of thrusters.

$\mathit{A}$ can be divided into the thrust efficiency matrix $\mathit{V}$ and the moment efficiency matrix $\mathit{L}$. $$\begin{array}{}\text{(2)}& \mathit{A}=\begin{array}{c}\mathit{V}\\ \mathit{L}\end{array}={\mathit{a}}_1\cdots {\mathit{a}}_i\cdots {\mathit{a}}_n\end{array}$$

The installation of one thruster is determined in the $xyz$-coordinate system, as shown below.

In Figure 1, $O$ is the origin of the satellite body coordinate system. $C$ is the center of mass of the satellite. $T$ is the installation position of the thruster. The $xyz$-coordinate system is parallel to the $x_b{y}_b{z}_b$-coordinate system. $\theta$ is the angle between the thruster direction and the $x$-$y$ plane. $\phi$ is the angle between the projection of the thruster in the $x$-$y$ plane and the $x$-axis.

[figure(s) omitted; refer to PDF]

Thrust efficiency matrix $\mathit{V}={\mathit{v}}_1\cdots {\mathit{v}}_i\cdots {\mathit{v}}_n$, where ${\mathit{v}}_{\mathit{i}}={v_{xi},v_{yi},v_{zi}}^T$ is the unit thrust direction vector of a single thruster relative to the $x_b{y}_b{z}_b$-coordinate system. ${\mathit{v}}_{\mathit{i}}$ satisfies ${v_{xi}}^2+{v_{yi}}^2+{v_{zi}}^2=1$, determined by the ${\theta }_i$ and ${\phi }_i$. It can be described as $$\begin{array}{}\text{(3)}& \begin{array}{c}{v}_{xi}=\cos {\theta }_i\cos {\phi }_i\hfill \\ v_{yi}=\cos {\theta }_i\sin {\phi }_i\hfill \\ v_{zi}=\sin {\theta }_i\hfill \end{array}\end{array}$$

There are $m$ installation positions of thrusters. In Figure 1, $\mathit{r}$ represents the position vector of the thruster in the $x_b{y}_b{z}_b$-coordinate system. ${\mathit{r}}_j={x_j,y_j,z_j}^T$ moment efficiency matrix $\mathit{L}={\mathit{r}}_1\times {\mathit{v}}_1\cdots {\mathit{r}}_j\times {\mathit{v}}_i\cdots {\mathit{r}}_m\times {\mathit{v}}_n$.

The efficiency matrix $\mathit{A}$ is composed of the thrust efficiency matrix $\mathit{V}$ and the thrust moment efficiency matrix $\mathit{L}$. The expression of $\mathit{A}$ is as follows: $$\begin{array}{}\text{(4)}& \mathit{A}=\begin{array}{c}\mathit{V}\\ \mathit{L}\end{array}=\begin{array}{ccccc}{\mathit{v}}_1& \cdots & {\mathit{v}}_i& \cdots & {\mathit{v}}_n\\ {\mathit{r}}_1\times {\mathit{v}}_1& \cdots & {\mathit{r}}_j\times {\mathit{v}}_i& \cdots & {\mathit{r}}_m\times {\mathit{v}}_n\end{array}\end{array}$$

The generalized active control force $\mathit{F}$ consists of $\mathit{f}$ and $\mathit{l}$, where $\mathit{f}=\mathit{V}\mathit{T}$ and $\mathit{l}=\mathit{L}\mathit{T}$. $\mathit{f}$ controls the three degrees of freedom of satellite translation motion. $\mathit{l}$ controls the three degrees of freedom of the satellite rotation. $$\begin{array}{}\text{(5)}& \mathit{A}\mathit{T}=\begin{array}{c}\mathit{V}\\ \mathit{L}\end{array}\begin{array}{c}{T}_1\\ ⋮\\ T_n\end{array}=\begin{array}{c}\mathit{f}\\ \mathit{l}\end{array}=\mathit{F}\end{array}$$

Accordingly, the control allocation is established as the following model: $$\begin{array}{}\text{(6)}& \begin{array}{c}{F}_c=AT\hfill \\ \underset{¯}{T}\le T_i\le \overline{T}\hfill \end{array}\end{array}$$

The disturbing force measured is taken as the desired control vector ${\mathit{F}}_{\mathit{c}}$. The optimal solution $\mathit{T}$ is calculated considering specific constraints subject to the thrust limits. $T$ represents the lower limit of thrust, and $\overline{\text{T}}$ represents the upper limit of thrust. Based on Formula (6), the allocation errors of thruster configuration could be analyzed.

3. The Design Method of Thruster Layout

The advantage of a redundant thruster layout is that it could effectively solve problems such as failure or damage of thrusters. However, there exist errors between the actual result $\mathit{F}$ and the desired control vector ${\mathit{F}}_{\mathit{c}}$, called the allocation error. To reduce allocation errors, the number and the installation orientation of thrusters are designed in this paper.

3.1. The Allocation Error

The allocation error refers to the difference between the actual outcome $\mathit{F}$ and the desired control vector ${\mathit{F}}_{\mathit{c}}$. Assuming that the thrust output error of thrusters is $\mathit{\delta }\mathit{T}$ and the installation error matrix is $\mathit{\delta }\mathit{A}$, the actual result $\mathit{F}$ can be expressed as $$\begin{array}{}\text{(7)}& \mathit{F}=\mathit{A}+\mathit{\delta }\mathit{A}\mathit{T}+\mathit{\delta }\mathit{T}\end{array}$$

Thus, the allocation error $\mathit{\Delta }$ is formulated as the following equation: $$\begin{array}{}\text{(8)}& \mathbf{\Delta }=\mathit{A}+\mathit{\delta }\mathit{A}\mathit{T}+\mathit{\delta }\mathit{T}-{\mathit{F}}_c\end{array}$$

With a suitable allocation algorithm, it can be considered reasonable that the optimal solution $\mathit{T}$ satisfies $\mathit{A}\mathit{T}={\mathit{F}}_{\mathit{c}}$. Consequently, the formulation of allocation error can be obtained. $$\begin{array}{}\text{(9)}& \mathbf{\Delta }=\mathit{A}\mathit{\delta }\mathit{T}+\mathit{\delta }\mathit{A}\mathit{T}+\mathit{\delta }\mathit{A}\mathit{\delta }\mathit{T}\end{array}$$

$\mathit{\delta }\mathit{T}={\delta T_1\cdots \delta T_i\cdots \delta T_n}^T$, where $\delta T_i$ represents the output error of the $i$-th thruster. $\mathit{\delta }\mathit{A}=\mathit{\delta }{\mathit{a}}_1,\cdots ,\mathit{\delta }{\mathit{a}}_i,\cdots ,\mathit{\delta }{\mathit{a}}_n$, where the $i$-th column vector $\mathit{\delta }{\mathit{a}}_i$ can be described as $$\begin{array}{}\text{(10)}& {\mathit{\delta }\mathit{a}}_i=\begin{array}{c}{\mathit{\delta }\mathit{v}}_i\\ {\mathit{r}}_j\times {\mathit{\delta }\mathit{v}}_i+{\mathit{\delta }\mathit{r}}_j\times {\mathit{v}}_i+{\mathit{\delta }\mathit{r}}_j\times {\mathit{\delta }\mathit{v}}_i\end{array}\end{array}$$

$\mathit{\delta }{v}_i$ and $\mathit{\delta }{\mathit{r}}_j$ can be written as follows: $$\begin{array}{}\text{(11)}& \begin{array}{c}{\mathit{\delta }\mathit{v}}_i=\begin{array}{c}\cos {\theta }_i+{\delta \theta }_i\cos {\phi }_i+{\delta \phi }_i-\cos {\theta }_i\cos {\phi }_i\\ \cos {\theta }_i+{\delta \theta }_i\sin {\phi }_i+{\delta \phi }_i-\cos {\theta }_i\sin {\phi }_i\\ \sin {\theta }_i+{\delta \theta }_i-\sin {\theta }_i\end{array}\hfill \\ {\mathit{\delta }\mathit{r}}_j={\delta x_j,\delta y_j,\delta z_j}^{\text{T}}\hfill \end{array}\end{array}$$

Ignoring the term composed of infinitesimals above the second order, Equation (11) can be simplified as follows: $$\begin{array}{}\text{(12)}& \begin{array}{c}{\mathit{\delta }\mathit{v}}_i=\begin{array}{c}-{\delta \phi }_i\sin {\phi }_i\cos {\theta }_i-{\delta \theta }_i\sin {\theta }_i\cos {\phi }_i\\ {\delta \phi }_i\cos {\phi }_i\cos {\theta }_i-{\delta \theta }_i\sin {\theta }_i\sin {\phi }_i\\ {\delta \theta }_i\cos {\theta }_i\end{array}\hfill \\ {\mathit{\delta }\mathit{r}}_j={\delta x_j,\delta y_j,\delta z_j}^{\text{T}}\hfill \end{array}\end{array}$$

Accordingly, $\mathit{\delta }\mathit{A}\mathit{\delta }\mathit{T}$ can be also ignored, and Equation (9) can be described as $$\begin{array}{}\text{(13)}& \mathbf{\Delta }=\mathit{A}\mathit{\delta }\mathit{T}+\mathit{\delta }\mathit{A}\mathit{T}\end{array}$$

3.2. TCF

Based on the above analysis, the TCF is proposed to adjust the position and angle of thrusters to decrease $\mathit{\Delta }$. The relation between TCF and $\mathit{\Delta }$ can be derived from the covariance matrix of both sides of Equation (13). $$\begin{array}{}\text{(14)}& \mathrm{cov}\mathbf{\Delta }=\mathrm{cov}\mathit{A}\mathit{\delta }\mathit{T}+\mathrm{cov}\mathit{\delta }\mathit{A}\mathit{T}+2\mathrm{cov}\mathit{\delta }\mathit{A}\mathit{T},\mathit{A}\mathit{\delta }\mathit{T}\end{array}$$

$\mathit{T}$ is independent of $\mathit{\delta }\mathit{T}$, and we can obtain $$\begin{array}{}\text{(15)}& \mathrm{cov}\mathit{\delta }\mathit{A}\mathit{T},\mathit{A}\mathit{\delta }\mathit{T}=0\end{array}$$

Thus, Equation (14) can be written as follows: $$\begin{array}{}\text{(16)}& \mathrm{cov}\mathbf{\Delta }=\mathit{A}\mathrm{cov}\mathit{\delta }\mathit{T}{\mathit{A}}^{\text{T}}+\mathit{\delta }\mathit{A}\mathrm{cov}\mathit{T}\mathit{\delta }{\mathit{A}}^{\text{T}}\end{array}$$

Equation (12) shows that $\mathit{\delta }{\mathit{r}}_j\times \mathit{\delta }{\mathit{v}}_i$ is composed of infinitesimals above the second order. It can be ignored. As a result, $\mathit{\delta }{\mathit{a}}_i$ of $\mathit{\delta }\mathit{A}$ can be divided into $\mathit{\delta }{\mathit{a}}_{1i}$ and $\mathit{\delta }{\mathit{a}}_{2i}$. $$\begin{array}{}\text{(17)}& {\mathit{\delta }\mathit{a}}_i={\mathit{\delta }\mathit{a}}_{1i}+{\mathit{\delta }\mathit{a}}_{2i}=\begin{array}{c}{\mathit{\delta }\mathit{v}}_i\\ {\mathit{r}}_j\times {\mathit{\delta }\mathit{v}}_i\end{array}+\begin{array}{c}\mathit{0}\\ {\mathit{\delta }\mathit{r}}_j\times {\mathit{v}}_i\end{array}\end{array}$$

In the above equation, $\mathit{\delta }{\mathit{a}}_{1i}$ is the error caused by the installation angle and $\mathit{\delta }{\mathit{a}}_{2i}$ is the error caused by the installation position.

Given that the mean and the standard deviation of the output error $\delta T_i$ of each thruster are zero and ${\sigma }_{\delta T}$, respectively, we can obtain $$\begin{array}{}\text{(18)}& \mathrm{var}\mathit{\delta }\mathit{{\rm T}}={\sigma }_{\delta {\rm T}}^2\mathit{I}\end{array}$$

Consequently, Equation (16) can be expressed as $$\begin{array}{}\text{(19)}& \mathrm{cov}\mathbf{\Delta }={\sigma }_{\delta {\rm T}}^2{\mathit{A}\mathit{A}}^{\text{T}}+{\mathit{\delta }\mathit{A}}_1+{\mathit{\delta }\mathit{A}}_2\mathrm{cov}\mathit{T}{{\mathit{\delta }\mathit{A}}_1+{\mathit{\delta }\mathit{A}}_2}^T\end{array}$$

According to the definition of covariance, $$\begin{array}{}\text{(20)}& \mathrm{cov}\mathbf{\Delta }=\begin{array}{ccc}\mathrm{cov}{\mathbf{\Delta }}_1,{\mathbf{\Delta }}_1& \cdots & \mathrm{cov}{\mathbf{\Delta }}_1,{\mathbf{\Delta }}_D\\ ⋮& \ddots & ⋮\\ \mathrm{cov}{\mathbf{\Delta }}_D,{\mathbf{\Delta }}_1& \cdots & \mathrm{cov}{\mathbf{\Delta }}_D,{\mathbf{\Delta }}_D\end{array}\end{array}$$

Diagonal elements cov(${\mathit{\Delta }}_i$, ${\mathit{\Delta }}_i$) represent the covariance of the $i$-th control vector, $i=1,2,\cdots D$. Based on the definition of covariance and variance, $$\begin{array}{}\text{(21)}& \mathrm{cov}{\mathbf{\Delta }}_i,{\mathbf{\Delta }}_i=\mathrm{var}{\mathbf{\Delta }}_i=E{\mathbf{\Delta }}_i^2-E{{\mathbf{\Delta }}_i}^2\end{array}$$

Because the mean of the thrust output error $\delta T_i$ is zero, $E{\mathit{\Delta }}_i=0$. Equation (21) can be written as the following expression: $$\begin{array}{}\text{(22)}& \mathrm{cov}{\mathbf{\Delta }}_i,{\mathbf{\Delta }}_i=E{\mathbf{\Delta }}_i^2\end{array}$$

Combined with Equation (22), the expression of the diagonal elements of Equation (20) can be obtained. $$\begin{array}{}\text{(23)}& \text{tr}\mathrm{cov}\mathbf{\Delta }=\sum _{i=1}^{D}E{\mathbf{\Delta }}_i^2\end{array}$$

In the second item on the right of Equation (19), $\mathrm{cov}\mathit{T}$ represents the balance degree of the thrust distribution result, which is related to the selection of algorithms, constraint conditions, and control instructions. The thrust distribution result is computed with the same instruction in the same algorithm. We can obtain $$\begin{array}{}\text{(24)}& \mathrm{cov}\mathit{T}=\mathit{0}\end{array}$$

In this paper, we define the TCF described as follows: $$\begin{array}{}\text{(25)}& \text{TCF}=\text{tr}{\mathit{A}\mathit{A}}^T=\sum _{i=1}^n{{\mathit{a}}_i}^2\end{array}$$

On the basics of Equations (23) and (24), the traces of both sides of Equation (19) can be written as $$\begin{array}{}\text{(26)}& \text{tr}\mathrm{cov}\mathbf{\Delta }=\sum _{i=1}^{D}E{\mathbf{\Delta }}_i^2=\text{tr}{\sigma }_{\delta {\rm T}}^2{\mathit{A}\mathit{A}}^T={\sigma }_{\delta {\rm T}}^2\text{TCF}\end{array}$$

It shows that $\mathit{\Delta }$ is positively correlated with TCF. Therefore, TCF could become a design critical standard measuring the effects of the thruster configuration on $\mathit{\Delta }$.

3.3. TCF Optimization Principle

In order to decrease $\mathit{\Delta }$, the two key criterions about the thruster layout are derived from TCF. And it is effective and feasible to evaluate the installation of thrusters in the design of the thruster layout.

3.3.1. Design Principle of Thruster Number

A redundant thruster layout does not mean that the more thrusters, the better stability, while it can provide fault tolerance. Reducing unnecessary thrusters reasonably is conducive to lowering $\mathit{\Delta }$. The proof is as follows:

If one thruster is added to a redundant thruster layout with $n$ thrusters, $\text{TC}{\text{F}}_{n+1}$ can be described as follows: $$\begin{array}{}\text{(27)}& {\text{TCF}}_{n+1}=\text{tr}\mathit{A}+{\mathit{a}}_{n+1}{\mathit{A}+{\mathit{a}}_{n+1}}^T=\sum _{i=1}^{n+1}{{\mathit{a}}_i}^2\end{array}$$

According to Equation (4), $$\begin{array}{}\text{(28)}& {\mathit{a}}_i=\begin{array}{c}{\mathit{v}}_i\\ {\mathit{r}}_j\times {\mathit{v}}_i\end{array}\end{array}$$

Due to ${v_{xi}}^2+{v_{yi}}^2+{v_{zi}}^2=1$, $$\begin{array}{}\text{(29)}& {{{\text{a}}_n}_{+1}}^2\ge {{{\text{v}}_n}_{+1}}^2=1\end{array}$$

Corresponding to Equation (27), $\text{TC}{\text{F}}_{n+1}$ of $n+1$ thrusters can be rewritten as the following equation: $$\begin{array}{}\text{(30)}& {\text{TCF}}_{n+1}={{\mathit{a}}_{n+1}}^2+\sum _{i=1}^n{{\mathit{a}}_i}^2\end{array}$$

On the basics of Equation (29), $$\begin{array}{}\text{(31)}& {\text{TCF}}_{n+1}=\sum _{i=1}^{n+1}{a_i}^2\ge 1+\sum _{i=1}^n{a_i}^2>{\text{TCF}}_n\end{array}$$

It can be proved that TCF increases with an increase in the number of thrusters. Besides, Equation (26) shows that $\mathit{\Delta }$ is positively related to TCF. Thus, it reveals that reducing unnecessary thrusters reasonably plays a critical role in lowering $\mathit{\Delta }$.

3.3.2. TCF Optimization of Thruster Configuration

Assuming that the number of thrusters is $n$, the control efficiency matrix $\mathit{A}$ is determined by the positions and angles of the thrusters. TCF can be set as the objective function to decrease $\mathit{\Delta }$. The corresponding optimization function of designing the thruster layout can be described as $$\begin{array}{}\text{(32)}& \begin{array}{c}J=\min \sum _{i=1}^n{{\mathit{a}}_i}^2\\ s.t.\hspace{1em}{\mathit{a}}_i\in \int \mathit{\Theta }\end{array}\end{array}$$where $\mathit{\Theta }$ represents the set of control vectors of thrusters applied on the satellite.

On the condition that the number of thrusts is constant, the position and angle of each thruster are limited by other hardware loads, which can be converted into the following constraints: $$\begin{array}{}\text{(33)}& \begin{array}{c}{\theta }_{i\min }\le {\theta }_i\le \le {\theta }_{i\max }\hfill \\ {\phi }_{i\min }\le {\phi }_i\le {\phi }_{i\max }\hfill \\ x_{j\min }\le x_j\le x_{j\max }\hfill \\ y_{j\min }\le y_j\le y_{j\max }\hfill \\ z_{j\min }\le z_j\le z_{j\max }\hfill \end{array}\end{array}$$

The above inequalities could construct simple linear matrix inequalities [22]. In order to improve the optimization efficiency of the model, it is useful to decrease the design parameters according to the actual situation, for instance, by minimizing the installation position and applying the paired thrusters. Hence, we can decrease the parameters of the constraints.

3.4. The Design Method of the Thruster Redundancy

The layout design of $n$ thrusters accommodates up to five $n$ parameters. These design parameters coupled with installation errors can significantly impact distribution outcomes and $\mathit{\Delta }$. Therefore, the minimum required thrusters [3] need to be figured out.

The redundant thruster layout needs a certain number of thrusters, which is determined by redundancy and degrees of freedom. The relation between the design value of $R$ redundancy and $n$ thrusters required to control the motion of the drag-free satellite with $D$ degrees of freedom is $$\begin{array}{}\text{(34)}& n\ge 2R+D+1\end{array}$$

The sufficient and necessary condition for Equation (34) is that each set of $D-1$ linearly independent column vectors of control efficiency matrix $\mathit{A}$ can form a $D-1$ dimensional hyperplane, which splits $D$ dimensional space into the same two parts where both there must be at least $R+1$ column vectors corresponding to thrusters. The inequation (34) means that apart from arbitrary $R$ failed thrusters, the remaining $n$-$R$ thrusters can still control satellite motion with $D$ degrees of freedom.

3.4.1. Calculation of the Redundancy

According to the above sufficient and necessary condition, the redundancy calculation method of the thruster layout can be deduced.

$\mathit{A}$ is a $D\times n$ matrix and $\mathit{I}$ is the $D$-order identity matrix, $\mathit{I}=\mathit{e}1\cdots \mathit{e}D$. $\tilde{\mathit{A}}$ is formed by a set of $D-1$ column vectors from $\mathit{A}$. The remaining vectors form the matrix $\widehat{\mathit{A}}$.

If $D-1$ column vectors are linearly independent, there must be a matrix of rank $D$ composed of ${\mathit{e}}_j$ in $\mathit{I}$ and $\tilde{\mathit{A}}$, $j\in 1,2,\cdots D$. The matrix must satisfy $$\begin{array}{}\text{(35)}& R{\tilde{\mathit{A}},{\mathit{e}}_j}^T=D\end{array}$$

And the normal vector of the $D-1$ dimensional hyperplane formed by $\tilde{\mathit{A}}$ is $\mathit{l}$. The relation between $\tilde{\mathit{A}}$ and $\mathit{l}$ is as follows: $$\begin{array}{}\text{(36)}& {\tilde{\mathit{A}}}^T\mathit{l}=\mathit{0}\end{array}$$

Thus, the calculation of $\mathit{l}$ is $$\begin{array}{}\text{(37)}& \mathit{l}={{\tilde{\mathit{A}}}^T;{{\mathit{e}}_j}^T}^{‐1}\begin{array}{c}0\\ ⋮\\ 0\\ 1\end{array}\end{array}$$

Then, all column vectors in matrix $\widehat{\mathit{A}}$ are compared with the normal vector $\mathit{l}$ in order to determine the partitioned space where they are located. $$\begin{array}{}\text{(38)}& \mathit{B}={\stackrel{\wedge }{\mathit{A}}}^T\mathit{l}\end{array}$$

$\mathit{B}$ is the set of cosines of the angle between each column vector of matrix $\widehat{\mathit{A}}$ and the normal vector $\mathit{l}$. The number of positive values in $\mathit{B}$ is equal to the number of column vectors in half space where the normal vector $\mathit{l}$ is located. The number of negative values in $\mathit{B}$ is equal to the number of column vectors in the other half space. The number of zero in $\mathit{B}$ is the number of column vectors located in the $D-1$ dimensional hyperplane formed by $\tilde{\mathit{A}}$.

$P$ is the number of positive values in $\mathit{B}$. $N$ is the number of negative values in $\mathit{B}$. $Z$ is the number of zero in $\mathit{B}$. And all the matrixes $\tilde{\mathit{A}}$ satisfying Equation (35) need to be found out so as to obtain $P$ and $N$ of each $\mathit{B}$ corresponding to each $\tilde{\mathit{A}}$.

$\mathit{P}$ is the set of $P$. $\mathit{N}$ is the set of $N$. $\mathit{Z}$ is the set of $Z$. The relationship among the actual redundancy $\underset{̱}{R}$, $\mathit{P}$, and $\mathit{N}$ can be expressed as follows: $$\begin{array}{}\text{(39)}& \underset{¯}{R}+1=\min \mathit{P},\mathit{N}\end{array}$$

3.4.2. Conditions for Minimum Redundancy

The above calculation method is performed on the basics of the fact that $D-1$ column vectors in $\tilde{\mathit{A}}$ must be linearly independent. A set of linearly dependent column vectors cannot form a $D-1$ dimensional hyperplane. And subsequent calculations can be avoided, which has no impact on the calculation of redundancy. Because the column vectors inapplicable to Equation (35) exist in $\widehat{\mathit{A}}$ of Equation (38) and the corresponding value of $\mathit{B}$ is zero, consequently, the redundancy is determined by the number of positive and negative in $\mathit{B}$.

$\underset{̱}{R}$ represents the actual redundancy, and $R$ represents the design redundancy. The calculation results of Equation (39) are analyzed as follows:

• If $\underset{̱}{R}<0$, the designed thruster layout cannot control the satellite motion with $D$ degrees of freedom. So, it needs to be redesigned.

• If $0<\underset{̱}{R}<R$, the layout does not meet the requirements for the design redundancy. So, it needs to be redesigned.

• If $\underset{̱}{R}\ge R$, the thruster layout meets the requirements of control and design redundancy.

Given the analysis of the relation between TCF and the number of thrusters, as few thrusters as possible should be used to complete the control target. It demands that the actual redundancy $\underset{̱}{R}$ should be equal to the design redundancy $R$.

Assuming that $\underset{̱}{R}=R$, the layout can use the minimum number of thrusters to reduce the error caused by the excessive number of thrusters. Hence, the thruster number is $$\begin{array}{}\text{(40)}& n=2R+D+1\end{array}$$

Based on the redundancy calculation, the conditions corresponding to $n\ge 2R+D+1$ can be written as follows: $$\begin{array}{}\text{(41)}& \begin{array}{c}\min \text{Z}\ge 0\hfill \\ \max \text{P},\text{N}\ge \min \text{P},\text{N}\hfill \end{array}\end{array}$$

If $n$ is the minimum number of thrusters with $R$ redundancy, namely, $n=2R+D+1$, the corresponding conditions are as follows: $$\begin{array}{}\text{(42)}& \begin{array}{c}\max \mathit{Z}=0\hfill \\ \max \mathit{P},\mathit{N}=\min \mathit{P},\mathit{N}\hfill \end{array}\end{array}$$

Once the number $n$ of thrusters is equal to $2R+D+1$, the above equalities combined with the calculation of the redundancy (Equations (35)–(39) in Section 3.4.1) could be translated to the linear matrix inequality problem [22]. The first equation indicates that all combinations of any five control vectors corresponding to thrusters are linearly independent. It means that any of the combinations can compose a $D-1$ dimensional hyperplane. The second equation indicates that all elements of $\mathit{P}$ and $\mathit{N}$ are equal to $R+1$. The implication is that there are $R+1$ thrusters in each of the two equally divided spaces by a $D-1$ dimensional hyperplane.

3.5. Assumptions

According to all the above deduction and analysis, the corresponding assumptions about the design procedure of thruster configuration are provided as follows:

• Assumption 1 is to apply a precise and reasonable allocation algorithm. It ensures that $\mathit{A}\mathit{T}={\mathit{F}}_c$.

• Assumption 2 uses a small-angle approximation to ignore the second-order small quantity.

• Assumption 3 is that the mean of the output error of each thruster is zero and the standard deviation is ${\sigma }_{\delta T}$.

Hence, the constraints of Equation (32) include not only the hardware limits but also the above conditions corresponding to $n=2R+D+1$. Combined with the above three assumptions, we can get the following constraints: $$\begin{array}{}\text{(43)}& \begin{array}{c}{a}_i\in \int \Theta \hfill \\ \min Z=0\hfill \\ \max P,N=\min P,N\hfill \\ {\theta }_{i\min }\le {\theta }_i\le \le {\theta }_{i\max }\hfill \\ {\phi }_{i\min }\le {\phi }_i\le {\phi }_{i\max }\hfill \\ x_{j\min }\le x_j\le x_{j\max }\hfill \\ y_{j\min }\le y_j\le y_{j\max }\hfill \\ z_{j\min }\le z_j\le z_{j\max }\hfill \end{array}\end{array}$$

The first constraint means the application conditions of the thrusters, which ensure that each thruster avoids interference. The second and third equations in constraints are derived from conditions for minimum redundancy. The last five inequations could show that the installation position and angle are limited by the actual physical constraints.

4. Designs of Redundant Thruster Layout

According to the above method of designing redundant thruster layout, the corresponding whole simulation strategy is shown in Figure 2. And the drag-free satellite LPF [23] is taken as an example for designing the thruster layout in this paper. A series of calculations would be performed to verify the redundancy of the thruster layout.

[figure(s) omitted; refer to PDF]

4.1. The Model of Redundant Thruster Layout

The satellite model is shown in Figure 3.

[figure(s) omitted; refer to PDF]

In order to reduce design parameters, paired thrusters are applied to each position, and as few installation positions as possible are used. The designed redundancy is set to be zero or one. Seven layouts are designed on the satellite model, as shown in Figure 4.

[figure(s) omitted; refer to PDF]

Layouts a, b, and c are designed to compare the effect of thruster number on $\mathit{\Delta }$. Different from the first layout, the latter two layouts add two and three thrusters, respectively, to compare the TCF of the corresponding layout. In addition, layouts d and e are arranged to contrast the influence of the installation position of thrusters on $\mathit{\Delta }$ with layout b on the basic of the same thrusters. Besides, the difference between layout f and layout d is the installation angle, the same as the difference between layout c and layout g.

4.2. Results of Redundancy

According to the redundancy calculation method given in Section 3.4.1, the calculation process of the redundancy is shown in Figure 5.

[figure(s) omitted; refer to PDF]

The designed layout determines the efficiency matrix $\mathit{A}$, which can be calculated by Equations (2)–(4). Efficiency matrices corresponding to all designed layouts in the previous section can be calculated. But with more factors taken into consideration, designing a redundant thruster layout is more complex. Using the paired thrusters can decrease the installation positions. And the design criterion minimum required thrusters can limit the number of thrusters.

Each unit vector applied on layout a is shown in Tables 1, 2, 3, 4, 5, and 6.

Table 1

Unit vectors of each thruster of layout a.

No.	${\mathit{v}}_{\mathit{x}}$	${\mathit{v}}_{\mathit{y}}$	${\mathit{v}}_{\mathit{z}}$
T1	−0.078310	−0.864364	0.496732
T2	−0.787716	0.364364	0.496732
T3	0.787716	0.364364	0.496732
T4	0.078310	−0.864364	0.496732
T5	−0.709406	0.500000	0.496732
T6	0.709406	0.500000	0.496732

Table 2

The extra unit vectors of layout b compared to layout a.

No.	${\mathit{v}}_{\mathit{x}}$	${\mathit{v}}_{\mathit{y}}$	${\mathit{v}}_{\mathit{z}}$
T7	0.787716	−0.364364	0.496732
T8	0.078310	0.864364	0.496732

Table 3

The extra unit vectors of layout c compared to layout a.

No.	${\mathit{v}}_{\mathit{x}}$	${\mathit{v}}_{\mathit{y}}$	${\mathit{v}}_{\mathit{z}}$
T7	0.433013	0.250000	0.866025
T8	−0.433013	0.250000	0.8660254
T9	0	−0.500000	0.866025

Table 4

Unit vectors of layout d and layout e.

No.	${\mathit{v}}_{\mathit{x}}$	${\mathit{v}}_{\mathit{y}}$	${\mathit{v}}_{\mathit{z}}$
T1	−0.078310	−0.864364	0.496732
T2	−0.787716	0.364364	0.496732
T3	0.787716	−0.364364	0.496732
T4	0.078310	0.864364	0.496732
T5	0.787716	−0.364364	0.496732
T6	0.078310	0.864364	0.496732
T7	−0.787716	−0.364364	0.496732
T8	−0.078310	0.864364	0.496732

Table 5

Unit vector of layout f.

No.	${\mathit{v}}_{\mathit{x}}$	${\mathit{v}}_{\mathit{y}}$	${\mathit{v}}_{\mathit{z}}$
T1	0.669550	−0.602866	0.433884
T2	−0.187322	0.881281	0.433884
T3	0.187322	0.881281	0.433884
T4	−0.669550	−0.602866	0.433884
T5	−0.669550	0.602866	0.433884
T6	0.187322	−0.881281	0.433884
T7	−0.187322	−0.881281	0.433884
T8	0.669550	0.602866	0.433884

Table 6

The extra unit vectors of layout g compared to layout c.

No.	${\mathit{v}}_{\mathit{x}}$	${\mathit{v}}_{\mathit{y}}$	${\mathit{v}}_{\mathit{z}}$
T7	0.813798	0.469846	0.342020
T8	−0.813798	0.469846	0.342020
T9	0	-0.939693	0.342020

Then, the installation positions of the thrusters should be determined. The installation positions of layouts c and g are the same as layout a. The position vectors of the corresponding installation place to the satellite centroid coordinate system are ${\mathit{r}}_1=-1.05,-0.61,0.62$, ${\mathit{r}}_2=1.05,-0.61,0.62$, and ${\mathit{r}}_3=0,1.21,0.62$.

Different from layout a, the additional position vector of layout b is ${\mathit{r}}_4=1.05,0.61,0.62$.

Layout d is the same as layout f. And the position vectors of each installation placed in the satellite centroid coordinate system are ${\mathit{r}}_1=-1.05,-0.61,0.62$, ${\mathit{r}}_2=1.05,-0.61,0.62$, ${\mathit{r}}_3=1.05,0.61,0.62$, and ${\mathit{r}}_4=-1.05,0.61,0.62$.

The position vectors of Layout e are r₁= (-0.7, -1.21, 0.62), r₂= (0.7, -1.21, 0.62), r₃= (0.7, 1.21, 0.62) and r₄= (-0.7, 1.21, 0.62).

Throughout the use of Equation (4), the efficiency matrix $\mathit{A}$ corresponding to each thruster layout is deduced, which is input into the redundancy calculation. In addition, $D=6$. Thus, the calculation results are shown in Table 7.

Table 7

The checking calculation results of redundancy.

Layout	a	b	c	d	e	f	g
The calculated redundancy	−1	−1	−1	−1	−1	−1	−1
The actual redundancy	0	1	1	1	1	1	1

The difference between the calculated redundancy and the actual redundancy is caused by the virtual thruster produced by solar radiation pressure. Because the solar radiation pressure exerted on the drag-free satellite is a relatively stable force with respect to the SC, it can be regarded as a virtual thruster ${\mathit{T}}_v$, as shown in Figure 6.

[figure(s) omitted; refer to PDF]

The position vector of the virtual thruster ${\mathit{T}}_v$ is $\mathit{r}=0,0,-0.4\text{m}$. $\theta ={60}^{°}$. $\phi$ of ${\mathit{T}}_v$ in the SC body coordinate system varies with the satellite position in the orbit. The actual redundancy of each layout should be computed with the virtual thruster, which is reliable. Consequently, the actual redundancy of each thruster in Table 7 can fulfill the design requirement for six degrees of freedom.

5. Simulation

According to the simulation strategy and the control allocation model, a series of simulations of the above layouts are completed to evaluate the relationship between $\mathit{\Delta }$ and TCF. Besides, the propellant consumption of each layout is figured out, and the fault tolerance can be verified.

5.1. Simulation of the Orbit

The orbits of the Taiji Project [24] are taken for simulation test in this paper, as shown in Figure 7.

[figure(s) omitted; refer to PDF]

In the above figure, three SCs (SC1, SC2, and SC3) in an equilateral triangle formation [25] are behind the Earth at a certain angle $\gamma$ ($\gamma ={20}^{°}$). The arm length $l$ between each other is $3\times {10}^6$ km. The satellite formation is in the Earth’s orbit around the sun ($ɑ=1A.U.$). The angle $\varphi$ between the plane of the satellite formation and the ecliptic plane is 60° ($\varphi ={60}^{°}$). $i$ represents the orbital inclination of SC1, expressed as $$\begin{array}{}\text{(44)}& \tan i=\frac{l}{2a+l/\sqrt{3}}\end{array}$$

The semimajor axis $ɑ$ of the orbit satisfies the following equation: $$\begin{array}{}\text{(45)}& {a1+e}^2={a+\frac{l}{2}\sqrt{3}}^2+{\frac{l}{2}}^2\end{array}$$

Besides, the expression of $e$ is as follows: $$\begin{array}{}\text{(46)}& e=\sqrt{1+\frac{l}{\sqrt{3}a}+\frac{l^2}{3a^2}}-1\end{array}$$

$l=3\times {10}^6$ km. And the results of $i$ and $e$ are about 0.009968838702433 and 0.005839000853703, respectively.

The relation among the longitude of ascending node $\Omega$, the argument of pericenter $\omega$, and the mean anomaly $M$ is $$\begin{array}{}\text{(47)}& M+\omega +\Omega =2\text{k}\pi \end{array}$$

The results are shown in Table 8.

Table 8

The values of $\Omega$, $\omega$, and $M$.

Satellite	$\mathit{\Omega }$	$\mathit{\omega }$	$\mathit{M}$
SC1	270°	270°	180°
SC2	150°	270°	300°
SC3	30°	270°	60°

The formula of the eccentric anomaly $E$ is $$\begin{array}{}\text{(48)}& E-e\sin E=M\end{array}$$

And true anomaly ${\theta }_a$ can be calculated. $$\begin{array}{}\text{(49)}& \tan \frac{{\theta }_a}{2}=\sqrt{\frac{1+e}{1-e}}\tan \frac{E}{2}\end{array}$$

Hence, the corresponding values of ${\theta }_a$ are 180° 299.5636083225304° and 60.436391677469594°, respectively.

Based on the above results, the orbits are simulated through MATLAB software and the TU Delft Astrodynamics Toolbox (TUDAT). The mass of each satellite model is 1500 kg. The diameter of the sail is 4.2 m. The operation time is 2 years. The step size is 1000 s. The simulation orbits are shown in Figure 8.

[figure(s) omitted; refer to PDF]

5.2. Simulation of the Disturbing Force

The dynamics of the drag-free satellite around the sun is $$\begin{array}{}\text{(50)}& {\mathit{M}}_{sc}{\ddot{\mathit{q}}}_{sc}={\mathit{F}}_c+{\mathit{f}}_d\end{array}$$

${\mathit{q}}_{sc}$ is the generalized displacement of the satellite. ${\mathit{M}}_{sc}$ is the generalized mass of the satellite. ${\mathit{f}}_d$ is the generalized disturbing force applied to the satellite. ${\mathit{F}}_c$ is the desired control force. The SC needs to become a static and stable platform. ${\mathit{f}}_d$ and ${\mathit{F}}_c$ should satisfy the following equation: $$\begin{array}{}\text{(51)}& {\mathit{F}}_c+{\mathit{f}}_d=0\end{array}$$

5.2.1. Simulation of the Solar Radiation Pressure

The generalized disturbing force ${\mathit{f}}_d$ acting on the satellite mainly comes from solar radiation pressure [5]. Combined with the solar radiation pressure model [26], the perturbation force model ${\mathit{f}}_{srp}$ in this paper is $$\begin{array}{}\text{(52)}& {\mathit{f}}_{\mathit{s}\mathit{r}\mathit{p}}=\frac{{\text{W}}_0}{c}{\frac{1\text{A}.\text{U}.}{R_l}}^2{S}_{Z}2\frac{R_d}{3}+R_s\mathit{z}·\mathit{s}\mathit{z}+R_a+R_d\mathit{s}\max -\mathit{z}·\mathit{s},0\end{array}$$

The speed of light in vacuum is $c=299792458$ m/s. $1\text{A}.\text{U}.=1.495978707\times {10}^{11}$ m. The standard solar constant is $W_0=1361$ W/m² [27]. $R_l$ is the relative distance between the satellite and the sun. $S_z$ is the solar sail area. $R_a$ is the absorption coefficient. $R_d$ is the diffuse reflection coefficient. $R_s$ is the specular reflection coefficient. These three coefficients satisfy $R_a+R_d+R_s=1$. $\mathit{z}$ is the unit vector of the $z$-axis. $\mathit{s}$ is the unit direction vector of the satellite pointing towards the sun.

Three drag-free satellites construct an equilateral triangular formation. The calculation of the disturbing force acting on SC1 is the same as the other two. Assuming $R_s=0$ and $R_a=0.14$ [28], the calculation results of ${\mathit{f}}_{srp}$ in the ecliptic coordinate system are obtained, as shown in Figure 9.

[figure(s) omitted; refer to PDF]

5.2.2. Enter Instruction of Control Allocation

Results in Figure 9 show the axial components of solar radiation perturbation exerted on the satellite in the heliocentric ecliptic coordinate system. But Equation (6) is established in the satellite centroid coordinate system. Therefore, the coordinate system needs to be transformed by the coordinate transformation matrix. $$\begin{array}{}\text{(53)}& \mathit{f}={\mathit{C}}_{\mathit{n}/\mathit{i}}{\mathit{f}}_{\mathit{s}\mathit{r}\mathit{p}}=\begin{array}{c}{{\stackrel{\wedge }{\mathit{x}}}_i}^{\mathit{\tau }}\\ {{\stackrel{\wedge }{\mathit{y}}}_i}^{\mathit{\tau }}\\ {{\stackrel{\wedge }{\mathit{z}}}_i}^{\mathit{\tau }}\end{array}{\mathit{f}}_{\mathit{s}\mathit{r}\mathit{p}}\end{array}$$

The force $\mathit{f}$ in the generalized disturbing force ${\mathit{f}}_d$ is calculated, as shown in Figure 10.

[figure(s) omitted; refer to PDF]

The calculation of the torque $\mathit{l}$ is as follows: $$\begin{array}{}\text{(54)}& \mathit{l}={\mathit{r}}_{\mathit{c}}\times \mathit{f}\end{array}$$

${\mathit{r}}_{\mathit{c}}={0,0,-0.4\text{m}}^{Т}$. The torque $\mathit{l}$ in the generalized disturbing force ${\mathit{f}}_d$ is calculated, shown in Figure 11.

[figure(s) omitted; refer to PDF]

$\mathit{f}$ and $\mathit{l}$ compose ${\mathit{f}}_d$. According to Equation (51), ${\mathit{F}}_c$ can be figured out and used as the enter instruction to calculate the thrust distribution of each thruster at each time point.

5.3. Modified Pseudoinverse Allocation

In this paper, the modified pseudoinverse method is used to work out the thrust distribution results $\mathit{T}$ with the constraint condition ($T\le T_i\le \overline{\text{T}}$).

First, the calculation of the initial result ${\mathit{T}}_p$ is $$\begin{array}{}\text{(55)}& {\mathit{T}}_p={\mathit{A}}^{†}{\mathit{F}}_c={\mathit{A}}^{†}\begin{array}{c}\mathit{f}\\ \mathit{l}\end{array}\end{array}$$

${\mathit{A}}^{†}$ is the pseudoinverse matrix of $\mathit{A}$. ${\mathit{A}}^{†}={{\mathit{A}}^T\mathit{A}}^{-1}{\mathit{A}}^T$.

Second, the null space matrix needs to be calculated. Although the pseudoinverse method has high efficiency, it cannot solve the problem of limited thrust. The null space matrix $\mathit{N}$ can be deduced from the following calculation: $$\begin{array}{}\text{(56)}& \mathit{\text{AN}}=\mathbf{0},\mathit{N}\in {\mathit{R}}^{n\times n-m}\end{array}$$

$\mathit{N}=\mathit{n}1\mathit{n}2\cdots \mathit{n}n-m$.

Finally, the new thrust instruction $\mathit{T}$ could be derived by using $\mathit{N}$. The thrust boundary value is assigned to elements of the initial thrust result ${\mathit{T}}_p$, which cannot meet the constraint condition. The result $\mathit{T}$ is expressed as the following calculation: $$\begin{array}{}\text{(57)}& \mathit{T}={\mathit{T}}_p+\mathit{N}\mathit{K}\end{array}$$where $\mathit{K}$ is the adjustment coefficient matrix.

5.4. Simulation Results and Analysis

Because of a larger thrust range and lower noise [21], the cold gas microthruster is chosen as the actuator. The limited range of $T_i$ in Equation (6) is set as [1 μN, 100 μN]. Different simulations of all layouts in normal mode are compared to find the relationship between the TCF of each layout and $\mathit{\Delta }$ of the counterpart. And the simulation in thruster failure mode is performed to verify the fault tolerance of the redundant thruster layout.

5.4.1. Simulation Results of All layouts

In order to simulate the thrust output error, a random value taken from 0 to 0.1 μN is added to each element of the distribution result $\mathit{T}$. The actual thrust results of all layouts are shown in Figure 12.

[figure(s) omitted; refer to PDF]

From top to bottom and left to right, the seven pictures represent the thrust results of layouts a, b, c, d, e, f, and g. All the figures show that the actual thrust of each thruster in each picture fluctuates periodically within the restricted range, which meets the constraint condition. The smooth curves indicate that the correspondent thrusters are barely assigned a constant value equal to the lower limit or upper limit of thrust. In addition, if the thrust is in a flat line, the output of the thruster must be equal to the limits of the thrust.

Throughout the use of the above results, we can obtain the actual control forces and torques, which constitute the actual control vector. And the allocation error $\mathit{\Delta }$ is the difference between the actual control vector and the desired control vector. The results of $\mathit{\Delta }$ are shown in Figure 13.

[figure(s) omitted; refer to PDF]

Two pictures at each line in Figure 13 represent $\mathit{\Delta }$ of the correspondent layout. The left picture shows the error $\Delta$ of each control force. The right picture shows the error $\Delta$ of each control torque. The results show that the allocation error in all pictures is at the level of 10⁻⁷. The difference in pictures in the third line implies that thrust redistribution is executed in the simulation of layout c. With the increase in thruster number, thrust redistribution could be more likely to occur in the control process. But given the possibility of thruster failure, we should keep the balance between thruster number and allocation error.

5.4.2. Analysis of Allocations

The main source of noise is allocation error. And the influence of TCF on allocation error is studied by frequency spectrum analysis to obtain the ability to restrain the error of each layout. The spectrum of $\mathit{\Delta }$ corresponding to each layout is shown in Figure 14.

[figure(s) omitted; refer to PDF]

The left picture in each row represents the control force noise, and the right picture in each row represents the control torque noise. Except for layout c, the results of other layouts show that the noise from 0.1 to 1000 Hz is between 16 and 17 orders of magnitude. The results of layouts c and g with the same thrusters reveal that thrust redistribution also influences the error, which belongs to the algorithm error. Subsequently, the variance of $\Delta$ is calculated to compare the data from different layouts intuitively. And the sum of $\mathrm{var}{\Delta }_i$ is worked out. The value of TCF corresponding to each layout is computed. The calculations are shown in Table 9.

Table 9

The variance of $\mathit{\Delta }$ of each layout.

Layout	a	b	c	d	e	f	g
$$\mathrm{var}{\mathit{\Delta }}_i/{10}^{-15}$$	1.88	2.40	16.17	2.09	2.09	1.62	3.00
	1.89	2.63	16.19	2.90	2.94	3.77	2.98
	1.24	1.63	11.12	1.62	1.66	1.27	1.53
	0.70	0.98	8.23	1.10	1.21	1.61	1.95
	0.69	0.89	8.17	0.75	1.07	1.10	1.94
	3.71	4.94	31.36	5.00	6.59	7.14	3.70
$$\text{Sum}/{10}^{-15}$$	10.1	13.5	91.2	13.5	15.6	16.5	15.1
TCF	12.1	16.2	20.7	16.2	18.6	19.9	18.1

It can be seen that TCF is positively associated with thruster number, derived from the data in Columns 1, 2, and 3 corresponding to layouts a, b, and c. It verifies the validity of the deduction in Section 3.3.1. And from the data in Columns 2, 4, 5, and 6 in the table, it can be found that the sum of $\mathrm{var}{\mathit{\Delta }}_i$ augments with the increase of TCF on the condition of the same thrusters. In addition, because $\mathrm{cov}{\mathit{\Delta }}_i,{\mathit{\Delta }}_i=\mathrm{var}{\mathit{\Delta }}_i$, the sum of $\mathrm{var}{\mathit{\Delta }}_i$ is equal to $\text{tr}\mathrm{cov}\mathit{\Delta }$. Thus, it can be found that TCF is proportional to $\text{tr}\mathrm{cov}\mathit{\Delta }$, which demonstrates the deduction of Equation (26). Moreover, the sum of $\mathrm{var}{\mathit{\Delta }}_i$ of layout c is the maximum among all results. It illustrates that thrust redistribution could lead to an increase in the sum. The limitation of the obtained results is that the relation between TCF and allocation error must be evaluated on the condition of the same number of thrusters. It can be found that the TCF of layout g is smaller than the TCF of layouts e and f, while the thrusters of layout g are more than the thrusters of layouts e and f. But if the extra thrusters are added to the designed layout, TCF can be compared between them. For example, the results of layouts a, b, c, and g show that the TCF of layout c error is the biggest, as well as the allocation.

5.4.3. Analysis of Fuel Consumption and Fault Tolerance

On the basics of the above thrust distribution results, the total impulse consumed in 2 years can be calculated through the time integral of thrust. Assuming that the specific impulse $I_{sp}$ of the cold gas microthruster is 60 s ($N_2$) and $g$ is 9.8 N/kg, the total mass corresponding to each layout is shown in Table 10.

Table 10

Propellant consumption of each thruster layout.

Layout	Total impulse /N·s	Mass/kg	TCF
a	9878.85	16.80	12.1
b	9885.18	16.81	16.2
c	7494.40	12.75	20.7
d	9885.06	16.81	16.2
e	9855.10	16.81	18.6
f	11313.31	19.24	19.9
g	10742.96	18.27	18.1

The above results show that the fuel consumption of redundant thrusters is related not only to the number of thrusters but also to thruster configuration. And it can be seen that layouts a, b, and e consume almost the same mass of fuel, while they have different values of TCF. Besides, four layouts with eight thrusters indicate that fuel consumption has some relationship with thruster configuration. A reasonable configuration could balance thruster number, allocation error, and fuel consumption.

A redundant thruster layout could improve the fault tolerance caused by the permission of $R$ failure thrusters. $R$ of layout a is 0. Hence, one of the redundant thrusters designed in other layouts is chosen to be invalid. The results are shown in Figure 15.

[figure(s) omitted; refer to PDF]

Based on the calculation of the actual redundancy corresponding to each layout, all layouts but layout a can perform simulation experiments with one thruster failure. All the above pictures show the distribution results of each redundant thruster layout. The missing curve T7 in the first picture indicates that the Thruster 7 is out of use. The same goes for other pictures. Compared with the former simulation without failed thrusters, the latter results reveal that the redundant thruster layout has a certain fault tolerance, which could use other thrusters to compensate for the effect of the failed thruster on control stability. It demonstrates the feasibility of the designed layout with redundant thrusters.

6. Conclusion

This paper presents a design method for the redundant thruster layout to improve the control stability of the satellite and reduce the allocation error. In the proposed method, the TCF of redundant thruster layout is put forward and derived through theoretical analysis. The simulation results can verify that TCF can be set as the key criterion for evaluating the allocation error of the attitude and orbit control systems. And the geometric relationship between the thrusters should be designed to reduce the allocation error with a smaller TCF. In addition, the minimum number of thrusters satisfying the redundancy requirement is defined in the layout design. The calculation method for thruster redundancy is provided in this paper. The results of the layouts with one failed thruster verify the validity of the design criterion of minimum required thrusters. And the advantage of a redundant thruster layout is that it can effectively solve problems such as failure or damage of thrusters. Moreover, the simulation results can verify the feasibility and effectiveness of the proposed design method for attitude and orbit control and show the influence of thruster layout on fuel consumption. Future work will be to consider establishing the relation between fuel consumption and requests for allocation error to seek a more suitable design method for the redundant thruster layout for a wider practical application.

Funding

This work was supported by the Ministry of Science and Technology of China (grant number 2021YFC2202604).

References

[1] P. McNamara, G. Racca, "Introduction to LISA Pathfinder," ESA Science & Technology, vol. 1 no. 1, 2009.

[2] G. E. Morris, N. Dunbar, P. Bianco, C. Edwards, "Design of a cold-gas micropropulsion system for LISA Pathfinder," 49th AIAA/ASME/SAE/ASEE Joint Propulsion Conference, pp. AIAA 2013-2854, DOI: 10.2514/6.2013-3854, .

[3] B. S. Crawford, Operation and Design of Multi-Jet Spacecraft Control Systems, 1968.

[4] T. J. Sumner, D. N. A. Shaul, "The observation of gravitational waves from space using LISA," Modern Physics Letters A, vol. 19 no. 11, pp. 785-800, DOI: 10.1142/S0217732304013647, 2004.

[5] O. Jennrich, "LISA technology and instrumentation," Classical and Quantum Gravity, vol. 26 no. 15, pp. 153001-153032, DOI: 10.1088/0264-9381/26/15/153001, 2009.

[6] J. He, P. Liu, R. Gao, C. Xue, L. Ma, L. Duan, Q. Kang, "Research on the neutralization control of the RF ion micropropulsion system for the ‘Taiji-1’ satellite mission," Plasma Science and Technology, vol. 22 no. 9, article 094002,DOI: 10.1088/2058-6272/ab7dd5, 2020.

[7] Y. Xia, H. Jiang, R. Wu, X. Gao, Y. Wang, L. Yang, "Thruster layout design on GEO satellite," 2019 IEEE 3rd Information Technology, Networking, Electronic and Automation Control Conference (ITNEC), pp. 1676-1680, DOI: 10.1109/ITNEC.2019.8729047, .

[8] C. Gazzino, D. Arzelier, L. Cerri, D. Losa, C. Louembet, C. Pittet, "A three-step decomposition method for solving the minimum-fuel geostationary station keeping of satellites equipped with electric propulsion," Acta Astronautica, vol. 158, pp. 12-22, DOI: 10.1016/j.actaastro.2018.08.015, 2019.

[9] L. Li, J. Zhang, Y. Li, S. Zhao, "Geostationary station-keeping with electric propulsion in full and failure modes," Acta Astronautica, vol. 163, pp. 130-144, DOI: 10.1016/j.actaastro.2019.03.021, 2019.

[10] A. Weiss, U. V. Kalabić, S. Di Cairano, "Station keeping and momentum management of low-thrust satellites using MPC," Aerospace Science and Technology, vol. 76, pp. 229-241, DOI: 10.1016/j.ast.2018.02.014, 2018.

[11] Z. Wei, T. Long, R. Shi, X. Song, "Electric thruster configuration design optimization for geostationary satellites with robotic manipulators," Advances in Space Research, vol. 69 no. 10, pp. 3798-3813, DOI: 10.1016/j.asr.2022.02.048, 2022.

[12] M. Faure, D. Henry, J. Cieslak, P. Lachevre, F. Ankersen, "Optimization of spacecraft thrusters configuration under fault diagnosability and recoverability constraints," IEEE Transactions on Aerospace and Electronic Systems, vol. 59 no. 5,DOI: 10.1109/TAES.2023.3256976, 2023.

[13] W. C. Durham, "Constrained control allocation," Journal of Guidance, Control, and Dynamics, vol. 16 no. 4, pp. 717-725, DOI: 10.2514/3.21072, 1993.

[14] J. M. Buffington, D. F. Enns, "Lyapunov stability analysis of daisy chain control allocation," Journal of Guidance, Control, and Dynamics, vol. 19 no. 6, pp. 1226-1230, DOI: 10.2514/3.21776, 1996.

[15] D. Bindel, M. Schlotterer, S. Theil, "Thruster actuation algorithms for scientific space missions," 56th International Astronautical Congress of the International Astronautical Federation, the International Academy of Astronautics, and the International Institute of Space Law. AIAA 2005-C1.P.10, .

[16] T. A. Johansen, T. I. Fossen, "Control allocation—a survey," Automatica, vol. 49 no. 5, pp. 1087-1103, DOI: 10.1016/j.automatica.2013.01.035, 2013.

[17] J. Fan, J. Ma, Z. Zheng, M. Lv, "An overview of stability design methods for overactuated systems based on control allocation," Journal of System Simulation, vol. 22, 2010.

[18] B. Yang, J. Miao, Z. Fan, J. Long, X. Liu, "Allocation optimization strategy for high-precision control of picosatellites and nanosatellites," International Journal of Aerospace Engineering, vol. 2018,DOI: 10.31193/ssap.01.9787520126595, 2018.

[19] D. B. Doman, M. W. Oppenheimer, "Improving control allocation accuracy for nonlinear aircraft dynamics," AIAA Guidance, Navigation, and Control Conference and Exhibit, pp. 2002-4667, DOI: 10.2514/6.2002-4667, .

[20] L. Ma, J. He, C. Yang, L. Duan, Q. Kang, "Structure optimization of micro-newton class radio-frequency ion thruster," Journal of Propulsion Technology, vol. 42 no. 2, pp. 474-480, 2021.

[21] C. Yang, Study on Micro-Newton-Scale Cold Gas Thruster and Thrust Performance, 2022.

[22] C. Liu, X. Yue, K. Shi, Z. Sun, Spacecraft attitude control: a linear matrix inequality approach,DOI: 10.1016/B978-0-323-99005-9.00001-8, 2022.

[23] M. Armano, H. Audley, J. Baird, P. Binetruy, M. Born, D. Bortoluzzi, E. Castelli, A. Cavalleri, A. Cesarini, A. M. Cruise, K. Danzmann, "LISA Pathfinder micronewton cold gas thrusters: in-flight characterization," Physical Review D, vol. 12 no. 99, pp. 122003-122013, DOI: 10.1103/PhysRevD.99.122003, 2019.

[24] Y. Xia, Orbit Design and Optimization for the LISA Gravitational Wave Observatory, 2009.

[25] S. V. Dhurandhar, K. R. Nayak, S. Koshti, J.-Y. Vinet, "Fundamentals of the LISA stable flight formation," Classical and Quantum Gravity, vol. 22 no. 3, pp. 481-488, 2005.

[26] J. Wertz, Spacecraft Attitude Determination and Control, 1978.

[27] J. L. Wright, A. Kantrowitz, "Space sailing," Physics Today, vol. 45 no. 12, pp. 85-85, DOI: 10.1063/1.2809919, 1992.

[28] R. Gray, P. Kemp, S. Barraclough, LISA Pathfinder Solar Array Requirement Specification, 2006.