This study deals with the relative motion of two satellites in formation flying. One satellite, called leader, is assumed to be orbiting in a circular orbit, and the other satellite, called follower, is orbiting in the proximity of the leader’s orbit and is controlled to keep a desired relative orbit with respect to the leader. The relative motion of the follower satellite is described with respect to a leader-fixed frame, in which three orthogonal axes are defined: \( x\) -axis directs to a radial direction from the Earth, \( z\) -axis corresponds with the direction of the orbital angular momentum vector of the leader, and \( y\) -axis completes the right-handed coordinates (Fig. 1). For simplicity, no disturbances on the satellites, such as the effect of the Earth oblateness and solar radiation pressure, are assumed. The equations of in-plane motion of the follower, called Hill– Clohessy– Wiltshire (HCW) equations ([15]), are described as

where the state variable vector \( \boldsymbol{x}\) consists of the follower position and velocity in the leader-fixed frame and \( \Omega\) is the orbital rate of the leader satellite. The control input \( \boldsymbol{u}\) is assumed to be able to generate arbitrary accelerations. This paper considers the relative in-plane motion alone because the cross-track motion along the \( z\) -axis is a simple harmonic motion and is decoupled from the in-plane motion. The analytical solution of Eq. (1) is described as follows ([16]).

where \( c_{t}:=\cos\left[\Omega\left(t-t_{0}\right)\right]\) , \( s_{t}:=\sin\left[\Omega\left(t-t_{0}\right)\right]\) , and \( t_{f}\) and \( t_{0}\) are a terminated and initial time of a reconfiguration maneuver, respectively. The matrix \( \Phi\left(t,t_{0}\right)\) is the state transition matrix connecting the states \( \boldsymbol{x}\left(t_{0}\right)\) and \( \boldsymbol{x}\left(t\right)\) with no inputs, i.e., \( \boldsymbol{x}\left(t_{f}\right)=\Phi\left(t,\, t_{0}\right)\boldsymbol{x}\left(t_{0}\right)\) . The state transition matrix consists of a fundamental matrix \( \Phi_{xy}\) as

where

Pan/Zoom

The analytical solution of HCW equations without external forces is rewritten as

where

Equations (4) and (5) show that the follower satellite relatively flies around the leader when \( b=0\) . The relative orbit is ellipse with the ratio \( \left|x\left(t\right)\right|/\left|y\left(t\right)\right|=0.5\) and its center position is written as \( \left(x,\, y\right)=\left(d,\,0\right)\) . This paper considers a formation reconfiguration and assumes that the follower initially fly around the leader satellite, i.e., \( d=b=0\) .

In formation flying, input directional constraints should be considered according to mission requirements. For instance, an input directional constraint with respect to an inertial frame stems from mission requirements for astronomical observation or space solar power system. The limitation of the attitude change due to the power generation with fixed solar arrays restricts the thrust direction with respect to the inertial frame. The follower satellite thus needs to control a relative orbit while keeping the thrust direction along a specified direction in the inertial frame.

For redundancy, many thrusters are equipped on a satellite so that they can generate translational forces and rotational torques in arbitrary directions. This indicates that translational forces in the \( x\) -\( y\) plane can be generated by properly distributing the magnitude of each thrust. Thus this study considers thruster number and relative motion of satellites and thruster forces in the \( x\) -\( y\) plane without loss of generality.

Since the inertial frame can be set arbitrarily, without loss of generality, this paper assumes that the desired direction is the positive direction of the \( X\) -axis in the inertial frame as described in Fig. 2(a). In other words, by defining the thrust directional angle with respect to the inertial frame as

the input directional constraint denotes that the angle \( \delta_{i}\) should be kept zero throughout the maneuver. In Eq. (7), \( \theta\) is the true anomaly of the leader satellite.

In addition, this study investigates a formation control under the input directional constraint with respect to the leader satellite. Such constraint arises in an autonomous rendezvous and on-orbit servicing mission for avoiding plume impingement and/or contamination from the follower’ s thrusts to the leader. The thrust direction of the follower with respect to the leader is limited, and the follower must control the relative orbit to a desired orbit while keeping the thrust vector not to orient into the leader as shown in Fig. 2(b).

The angle between the position vector of the follower \( \boldsymbol{x}_{{\rm {pos}}}=\left[x,y\right]^{T}\) and the acceleration vector \( \boldsymbol{u}\) is defined as

Since the thruster plume is oriented into \( -\boldsymbol{u}\) , the forbidden region of the thruster plume is described as \( \left|\delta_{l}\right|\leq\bar{\delta}\) , where \( \bar{\delta}\) is an upper bound angle. The input directional constraint with respect to the leader satellite requires that the angle \( \delta_{l}\) should avoid being within the angle \( \bar{\delta}\) throughout the maneuver.

Pan/Zoom

The energy optimal controller in this paper is based on the method proposed by \cite{Palmer:2006tt}. The controller in (\cite{Palmer:2006tt}) uses the transformation of boundary conditions to facilitate the derivation. Such transformation causes a singularity that disables setting a maneuver terminated time \( t_{f}\) as the multiple of the orbital period, i.e., setting \( t_{f}=2\pi N/\Omega\) \( (N=1,2,\dots)\) is impossible. This singularity can be avoided by formulating the optimal controller without the transformation (\cite{Cho:2012ik}). The advantage of the singularity avoidance is that setting \( t_{f}=2\pi N/\Omega\) simplifies the calculations for satisfying the input directional constraints, as shown in later.

An energy optimal controller is designed for minimizing the following cost function:

In Eq. (9) and the following discussions, the initial time is set to \( t_{0}=0\) for simplicity. The control input is described with the Fourier series for expressing arbitrary inputs as

where \( a_{j0,}\, a_{jn},\, b_{jn}\,\left(j=x,y\right)\) are the Fourier coefficients. By substituting Eqs. (10) and (11) into Eq. (9) and using the Parseval’s theorem (\cite{Boas:2006ty}), the cost function in Eq. (9) is rewritten as follows.

The relative orbit control must satisfy a desirable state \( \boldsymbol{x}_{d}\) at the terminated time \( t_{f}\) . From the analytical solution of the HCW equations in Eq. (1), the follower state at the terminated time is written as

where \( \boldsymbol{x}_{h}\) is the homogeneous solution of the HCW equations. The boundary condition vector \( \boldsymbol{Q}\) is defined as follows.

The integrand in Eq. (14) is calculated as

The substitution of Eqs. (10), (11), and (15) into Eq. (14) yields

The definitions of the variables in Eq. (16) are summarized in Appendix.

The optimal control problem is finding the control inputs \( \boldsymbol{u}\) that minimizes the cost function in Eq. (9) and satisfies Eq. (13) that is calculated with given boundary states and terminated time. Furthermore, since the control inputs are expressed with the Fourier series as shown in Eqs. (10) and (11), the optimal control problem is reduced to finding the Fourier coefficients.

The Hamiltonian for minimizing the cost function (9) is described as

where \( \boldsymbol{\lambda}=\left[\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}\right]^{T}\) denotes the Lagrange multipliers. Substituting Eqs. (10) and (11) into Eq. (17) and the differentiating the Hamiltonian with respect to each Fourier coefficient yield the following equations.

From Eqs. (19)– (24), the Fourier coefficients are expressed with the Lagrange multipliers as follows.

Substituting these Fourier coefficients into Eq. (14) and using the Parseval’s theorem, we obtain the following equation.

This algebraic equation can be analytically solved for the Lagrange multipliers \( \boldsymbol{\lambda}\) with arbitrary initial and desired states and terminated time (see Appendix for the components of the matrix \( A\) ).

By substituting the multipliers back into Eqs. (25)– (30), we obtain the analytical forms of the Fourier coefficients, which are equivalent to the optimal control inputs:

or alternatively,

where

Note that, although the optimal controller in Eqs. (33) and (34) has the same form shown in (\cite{Palmer:2006tt}), the derived controller in the current paper has no singularity, that is, the terminated time can be set to \( t_{f}=2\pi N/\Omega\) . Such terminated time can simplify further calculations for satisfying the input directional constraints as shown in the following subsections.

Since the desired direction, the \( X\) -axis, is the axis in the inertial frame, the \( X\) -axis expressed in the leader-fixed frame is not constant but monotonically varying due to the orbital motion of the leader satellite as \( \dot{\psi}_{d}=-\Omega\) . If the input direction with respect to the leader-fixed frame also varies monotonically, the constraint on the input direction is expected to be satisfied, which is the control technique used in ([9]). In fact, the input direction can be designed for changing in one direction, because the derived controller in Eqs. (32) and (33) forms an ellipse when \( \lambda_{1}=0\) and furthermore an origin-centered ellipse when \( \lambda_{1}=\lambda_{4}=0\) . The drawback is that, the time derivative of the input angle is not constant, whereas that of the desired direction is constant. Such difference causes the oscillation of the input direction around the desired one ([9]). To solve this problem, an additional input that orients the input direction at the same rate as the orbital rate is used in this paper.

The boundary condition for \( \lambda_{1}=0\) has been previously derived by [9] as follows.

Starting and terminating the reconfiguration maneuver at these boundary states realize \( \lambda_{1}=0\) , resulting in an elliptic input direction as:

Under the assumption that \( \lambda_{4}\approx0\) , the input angle in the leader-fixed frame and its derivative are described as

Since \( \dot{\psi}_{d}=-\Omega\) , the input angle does not correspond with the desired direction throughout the maneuver even if the initial directions coincide. It is noted that, although the multiplier \( \lambda_{4}\) is not exactly zero, the larger \( N\) is, the smaller the multiplier \( \lambda_{4}\) becomes as seen in the 4th row of Eq. (63). That is, setting long terminated time \( t_{f}\) \( \left(=2N\pi/\Omega\right)\) makes the assumption \( \lambda_{4}\approx0\) reasonable.

To orient the input direction at the same rate as the desired direction, the following additional input \( \boldsymbol{u}_{{\rm {add}},i}\) is superimposed onto the nominal optimal input \( \boldsymbol{u}\) .

The resulting input is described as

This additional input can make the input direction circular, which has the same rate as the orbital rate of the leader. In fact, the time derivative of the input angle is calculated as \( \dot{\psi}=-\Omega\) using Eqs. (42) and (43).

The additional input affects the transition of the follower to the target orbit and causes an error. Due to this error, the relative orbit of the follower does not converge to a relative target orbit. Thus this undesirable error is compensated by setting a virtual target orbit. The effect of the additional input to the follower position can be calculated as follows.

Although the effect of the additional input is analytically described in Eq. (44), the convergence of the follower to the target relative orbit is not guaranteed due to the additional input. Thus the state variation due to the additional inputs should be reduced.

Setting the terminated time \( t_{f}=2\pi N/\Omega\) can simplify Eq. (44) as

From Eqs. (13) and (30), the Lagrange multipliers can be easily calculated using MAPLE or Mathematica as follows:

Then, the variable \( \Lambda\) is also obtained from Eq. (34) as

Substituting \( \Lambda\) into the first component of Eq. (45) yields the virtual target position along the \( x\) -axis as

where \( x_{d}\) is the target position along the \( x\) -axis. The formation reconfiguration from \( \boldsymbol{x}\left(0\right)\) to \( \boldsymbol{x}_{d}\) is achieved using the suboptimal input in Eqs. (42) and (43), in which the virtual target orbit is used.

The proposed control method under the constraint with respect to the inertial frame is summarized as follows. The terminated time of the maneuver \( t_{f}\) is determined as the multiple of the orbital period. Then, substituting \( t_{f}\) into Eqs. (36) and (37) provides the maneuver starting position and terminated one of the follower. The control input shown in Eqs. (42) and (43) is implemented, in which the Lagrange multipliers are calculated using the virtual target position along the \( x\) -axis in Eq. (48)

This subsection extends the control method to accommodate the input directional constraint with respect to the leader-fixed frame. Since the follower satellite is relatively orbiting around the leader, the forbidden region of the thrust direction monotonically varies. Thus the input direction that changes elliptically is exploited also to keep the input direction within the desired one in the leader-fixed frame.

The initial input direction must be in the admissible region to satisfy the input directional constraint throughout the maneuver. From Eqs. (32) and (33) and \( \lambda_{1}=0\) , the initial input angle is expressed with

The input angle \( \psi\left(0\right)\) is a function of the Lagrange multipliers, and furthermore the multipliers are the function of the terminated time \( t_{f}\) , making \( \psi\left(0\right)\) highly nonlinear. Nevertheless the initial input direction can be simplified by setting the terminated time as \( t_{f}=2\pi N/\Omega\) , which is the reason that the optimal controller without singularities is derived. Substituting \( t_{f}=2N\pi/\Omega\) into Eq. (36) yields

The initial input angle \( \psi\left(0\right)\) is then determined as

Therefore the initial thrust direction \( -\boldsymbol{u}\) is analytically obtained and is orthogonal to the follower position vector regardless of the initial and target orbit.

Both the follower trajectory and the input direction vary elliptically. However, the maneuver starts when the follower is on the \( x\) -axis as derived in Eq. (50), whereas the initial input angle is \( \pi/2\) as shown in Eq. (51). Thus the initial phase angle differs by 90 deg and this different phase causes an error. To eliminate this phase difference, an additional input is superimposed onto the nominal inputs described in Eqs. (32) and (33).

Since the nominal input direction becomes ellipse with the ratio \( \left|u_{x}\right|/\left|u_{y}\right|=1/2\) , the additional input making the ratio \( \left|u_{x}\right|/\left|u_{y}\right|=2\) is required for compensating the phase difference. The following additional input for the acceleration along the \( y\) -axis is added to the nominal input \( \boldsymbol{u}\) .

This additional input can make the input direction ellipse at the ratio \( \left|u_{x}\right|/\left|u_{y}\right|=2\) , resulting in the same phase as the formation reconfiguration of the follower.

In the similar manner to the previous subsection, the effect of the additional input to the follower position can be calculated as follows.

It is noted that the position variation \( \Delta\boldsymbol{x}\) in Eq. (53) is also simplified by setting the terminated time \( t_{f}=2\pi N/\Omega\) , whereas using the other terminated time results in a more complicated form. The virtual target orbit to compensate the variation in Eq. (53) can be analytically obtained as

The proposed control method under the constraint with respect to the leader-fixed frame is described as follows. The terminated time of the maneuver \( t_{f}\) is determined as the multiple of the orbital period. The maneuver starting and terminated positions are obtained by substituting \( t_{f}\) into Eqs. (36) and (37). In the maneuver, the control input \( \boldsymbol{u}+\boldsymbol{u}_{{\rm {add},}l}\) is implemented, in which the Lagrange multipliers are calculated using the virtual target position along the \( x\) -axis in Eq. (54).