Sensor-based Optimal Attitude Reorientation Control Scheme Based on Computational Programming Approach

Sensor-based optimal spacecraft attitude reorientation control by momentum exchange based on a computational programming approach is addressed in this study. The control problem of a rigid spacecraft actuated by more than three reaction wheels with an open time of maneuver is considered. The modified Rodrigues parameters for large principal rotations are applied to derive our kinematical model. The introduced algorithm can be realized by attitude sensors, such as rate gyros, with an appropriate arrangement. The cost function to be minimized is defined as a weighted performance index of the time of the maneuver and the integral of the squared sum of wheel-torque magnitudes. Instead of utilizing Pontryagin’s minimum principle, an iterative procedure is used to reformulate and solve the optimal reorientation control problem as a constrained nonlinear programming problem. To show the feasibility of the proposed method, numerical simulated results are included for illustration.


Introduction
Many current and future spacecraft have mission requirements of pointing, tracking, and multitarget acquisition within the physical limits of actuators. To achieve these tasks, sequences of revolutions about body principal axes are required to generate the required attitude reorientation maneuver. During the spacecraft attitude maneuver, the directional vectors of bright celestial objects, such as the Sun and the Earth, are measured by optical sensors (e.g., infrared sensors or star sensors). Then, the appropriate sensors communicate with an electronic module, which records the timing of triggers on the Sun and the Earth. The attitude can be determined by processing the recorded data. Many studies considering a rigid spacecraft with three attached control torques perpendicular to the three principal axes have been published, including studies on time-optimal attitude control, (1,2) energy-optimal control, (3) feedback control with input saturation, (4) finite-time output feedback control, (5) adaptive sliding mode control, (6)(7)(8)(9)(10) and path-planning control. (11,12) Because of their smooth operating modes, reaction wheel systems are usually utilized by momentum exchange to achieve and maintain the precise attitudes of a spacecraft. The control input signals of a reaction wheel system are the wheel torques. Many schemes of actuator fault detection (10) have been introduced to determine whether the control input signals are working precisely. Normally, three reaction wheels are used to control a spacecraft with the wheel axes aligned with the body principal axes. In practice, reaction wheel systems are assembled in a specific configuration with K wheels, where K > 3 provides redundancy in three-axis control. The K reaction wheels can be mounted in arbitrary orientations in the spacecraft frame, but at least three wheel directions have to be linearly independent for maneuver controllability.
Different aspects of spacecraft attitude control actuated by multiple reaction wheels have been addressed. Studies include spacecraft large-angle attitude control with a four-wheel pyramid configuration by the application of anti-windup control and intelligent integrators, (13) and the synthesis of a four-wheel control algorithm by applying the energy-shaping technique. (14) However, rapid and energy-efficient large-angle attitude maneuvers are important in many spacecraft missions. Consequently, research on attitude maneuvers of spacecraft actuated by multiple reaction wheels utilizing optimal control has also been reported. For example, Vadali and Junkins addressed the optimal large-angle reorientation of a rigid asymmetric spacecraft with multiple reaction wheels by using an integral of a weighted quadratic function for control as the cost function. (15) Carrington and Junkins developed a means of nonlinear optimal polynomial feedback control for the large-angle attitude maneuver of a spacecraft with four reaction wheels. (16) Wei and Lu introduced nonlinear feedback control logic for the near-minimum-time eigenaxis reorientation control of a rigid spacecraft actuated by four reaction wheels. (17) In Ref. 18, the problem of a constrained minimum-time maneuver by reaction wheels was solved by particle swarm optimization. In Ref. 19, kinematic steering attitude control by reaction wheels was proposed in the form of an inner-outer control loop.
In this study, the optimal large-angle reorientation attitude control (OLRAC) problem of a rigid spacecraft actuated by K reaction wheels between two different attitudes is considered. The rest-to-rest reoriented wheel-torque control inputs are determined to minimize a weighted performance index of the time of maneuver and the integral of the squared sum of wheel-torque magnitudes. Owing to the difficulty of applying Pontryagin's minimum principle (PMP), (20) the OLRAC problem of an asymmetric rigid spacecraft actuated by multiple reaction wheels with the free time of maneuver has not attracted much attention. Thus, in this paper, we extend the concept in Refs. 21 and 22 to solve the OLRAC problem of a rigid spacecraft with an open time of maneuver that is actuated by K > 3 reaction wheels mounted in arbitrary and linearly independent orientations.
We introduce a novel iterative procedure to reformulate and solve the OLRAC problem as a constrained nonlinear programming (NLP) problem. Then a method to generate initial feasible solutions of the NLP problem by using modified Rodrigues parameters (MRP) is also proposed. Since initial feasible solutions can be determined easily, the optimization of the NLP problem can be started from different points to find an optimal rest-to-rest reorientation maneuver between two different rotational attitudes under the constraints on wheel-torque inputs for a rigid spacecraft with K momentum wheels. Numerical simulations are performed to show the feasibility of the proposed method. The proposed algorithm can be realized by using attitude sensors, such as rate gyros, GPS receivers, and star sensors, with an appropriate arrangement. This paper is organized as follows. Section 2 reviews the spacecraft and reaction wheel dynamics and attitude representation by the MRP for large principal rotations. Then, the OLRAC of a rigid spacecraft between two rotational attitudes is formulated as an NLP problem in Sect. 3. In Sect. 4, our approach to generate initial feasible solutions of the NLP problem from different starting points is introduced. The procedure used to solve the OLRAC problem and the simulation results are shown in Sects. 5 and 6, respectively. Finally, conclusions are given in Sect. 7.

Spacecraft and Reaction Wheel Dynamics and Attitude Representation
The large-angle maneuver of a rigid spacecraft actuated by a set of inertial reaction wheels is considered in this section. The reaction wheel assembly is assumed to comprise K reaction wheels mounted with arbitrary orientations to provide redundancy for three-axis control. Usually, K > 3 wheels allows the failure of at most K − 3 reaction wheels. The nonlinear threeaxis rotational dynamics of a rigid spacecraft can be represented as (13,23) , where J B and J R are the inertia matrices of the spacecraft and reaction wheels, respectively. h and ω B = [ω x ω y ω z ] T are the angular momentum vector of the total system and the angular velocity vector of the body in the body coordinate frame, respectively. ω R = [ω 1 ω 2 ... ω K ] T is the wheel speed vector in the wheel coordinate frame. u = [u 1 (t) u 2 (t) ... u K (t)] T is the wheeltorque input vector and C is the 3 × K coordinate transformation matrix that describes the orientations of all reaction wheels in the spacecraft frame. At least three column vectors of C must be linearly independent or (CC T ) −1 must exist for independent three-axis control. The vector Cu in Eq. (1) represents the three principal axis components of the spacecraft control torque inputs. It is clear that matrix C is reduced to a third-order identity matrix for a rigid spacecraft equipped with three reaction wheels that are mounted coaxially with the origin of the principal axes at the spacecraft's center of mass. For a spacecraft-command control input vector u com , the wheel-torque input vector is u = C + u com , where C + = C T (CC T ) −1 is the pseudo inverse of transformation matrix C.
On the other hand, the attitude representation for a rigid spacecraft is described by the MRP. (24) The physical interpretation of the MRP arises from Euler's theorem. It states that the general displacement of a rigid body with one point fixed is uniquely determined by a principal unit vector e = [e 1 e 2 e 3 ] T and a principal angle of rotation Φ. The MRP vector σ [σ 1 σ 2 σ 3 ] T is related to the principal axis and principal angle by σ = etan(Φ/4), which is well defined for all principal axis rotations in the range 0° ≤ Φ ≤ 360°. The body angular velocity ω B and the kinematics of rigid body motion are related by the celebrated formula (24) ( ) , Here, S(•) is the skew-symmetric matrix operator and I 3×3 is the third-order identity matrix.

OLRAC Problem
The OLRAC problem of a rigid spacecraft actuated by K reaction wheels between two attitudes is to determine the wheel-torque inputs that will drive the rigid spacecraft system from an initial attitude to a desired final attitude, where the initial and final body angular velocities are assumed to be zero. The performance index to be minimized is defined as a weighted performance index of the time of maneuver and the integral of the squared sum of wheel-torque magnitudes. From the dynamics in Eqs. (1)-(5), the OLRAC problem can be formulated as follows: PROBLEM I: For a rigid spacecraft with the dynamics in Eqs. (1)-(5), assuming that the initial attitude and angular velocity vectors are given as where [σ 1_initial σ 2_initial σ 3_initial ] T and [σ 1_ final σ 2_ final σ 3_ final ] T represent the initial and the desired final attitude vectors of the spacecraft, respectively. In this problem, note that the time of maneuver t f is treated as a free variable and will be determined by an optimization process. The constant ρ is a weighting factor reflecting the relative importance of the control inputs with respect to the time of maneuver.
Problem I is clearly difficult to solve because of the nonlinear and coupled relation of the rigid spacecraft system. To cope with the difficulty, Problem I will be formulated and solved in the discrete-time domain by numerical methods. The first step is to divide the interval where N is the number of control steps. (21,22) That is, If the angular acceleration vector of the rigid spacecraft B  ω is assumed to be constant for each sub-interval, one obtains where ω B denotes ω B (i • Δt). Substituting Eq. (1) into Eq. (14) gives , one can obtain the iterative forms of the angular spinning rate vector of the reaction wheels as ( ) Assuming the time derivative of the MRP vector to be constant for each sub-interval, one also obtains Using the formula for the body angular velocity vector in Eq. (4) and the MRP vector in Eq. (5), one can represent Eq. (17) as 1 1 respectively. From the above equations, one finds that ω B (N) is a function of the angular velocity vectors ω B (0), ω B (1), ..., ω B (N − 1), the spinning rate vectors of reaction wheels ω R (0), ω R (1), ..., ω R (N − 1), the input vectors u(0), u(1), ..., u(N − 1), and the sampling period Δt. This means that In a similar way, one can also obtain Using Eqs. (19) and (21), Problem I is now transformed to a standard constrained NLP problem as follows: PROBLEM II: Given the initial attitudes in Eqs.
subject to where ω B (N) and σ(N) are respectively defined in Eqs. (19) and (21). Although the OLRAC problem of a rigid spacecraft actuated by K reaction wheels can be formulated as Problem II, there are still several difficulties to be overcome. One difficulty is the choice of the number of control steps N. A larger value of N will clearly give more freedom for the input variables. However, this also means a greater computational time for Problem II. For a linear system without constraints on the input variables, it has been shown that the initial value of N must be greater than the dimension of the state variables. (21,22) Although no similar criteria must be satisfied for nonlinear systems, an integer larger than the dimension of the state variables will be chosen as the initial value of N in this paper.

Initial Feasible Solutions to Form Different Starting Points
After formulating the OLRAC problem of a rigid spacecraft actuated by K reaction wheels as an NLP problem, there are only two ways to guarantee the global minimum of this NLP problem. The first one is to specify that the problem is convex (minimizing a convex objective function over a convex feasible region). If the problem is not convex, then the second way is to choose a sufficient number of different starting points so that all local minima can be obtained. (25) In the convex case, a local minimum is also a global minimum. In the nonconvex case, the global minimum is determined by evaluating the objective function at each local minimum. Since Problem II is clearly not a convex problem, only the second way can be used to find its global minimum.
In Problem II, an initial feasible solution is a set of u(0), u(1), ..., u(N − 1), and Δt satisfying the constraints in Eqs. (23)- (26). In this section, a systematic approach will be proposed to generate initial feasible solutions of Problem II. If initial feasible solutions can be found easily, then the optimization process of Problem II can be started from different points to find all local minima.
The first step of this approach is to find a maneuvering trajectory of the rigid spacecraft that satisfies the boundary conditions in Eqs. (6)-(8) and Eqs. (23)- (25), irrespective of the constraint in Eq. (26). If this trajectory satisfies the dynamics constraint in Eq. (26), then an initial feasible solution is found. Otherwise, this trajectory can be adjusted to a feasible one by reducing the velocities and accelerations but tracking the same positions. Such adjustment is similar to the situation of flying an airplane along a flight path. If a high flight speed cannot be achieved, then the airplane should fly at a lower speed with the same course.
To generate a maneuver trajectory that satisfies the boundary conditions in Eqs.
where the values of λ i , i = 2, ..., N − 1, are randomly selected from the interval [0, 1]. It is obvious that σ(0) and σ(N) are set to [σ 1_initial σ 2_initial σ 3_initial ] T and [σ 1_ final σ 2_ final σ 3_ final ] T in this procedure, respectively. Thus, the constraints in Eqs. (6) and (24) 7) and (25) are satisfied. Therefore, the trajectory determined by the above steps will satisfy the boundary conditions in Eqs. (6) and (7) and Eqs. (23)- (25). However, one still needs to check that the torque constraints in Eq. (26) are satisfied. Therefore, we discuss how to generate the corresponding input torques and how to adjust the input torques to make them feasible.
After determining the angular velocity vectors ω B (0), ω B (1), ..., ω B (N) and angular acceleration vectors  vectors u(0), u(1), ..., u(N − 1) can be computed sequentially provided that the value of ω R (0) given by Eq. (8) is known. Firstly, the initial total angular momentum vector h(0) can be obtained by substituting Eq. (8) into Eq. (2). The spacecraft command control input vector u com (0) is calculated as Thus, the wheel-torque input vector u(0) is determined as where C + = C T (CC T ) −1 is the generalized inverse matrix of the transformation matrix C.
Then h(0), and u(0) into Eq. (16). By applying the procedure repeatedly from i = 0 to i = N − 1, u(0), u(1), ..., u(N − 1) are determined sequentially. A flowchart to illustrate the generation of u(0), u(1), ..., u(N − 1) is shown in Fig. 1. In the above procedure, it is very likely that the input constraints in Eq. (26) will be violated. Therefore, a simple and effective method will be proposed for their adjustment. The basic idea is to increase the value of Δt to reduce the maneuvering velocities and accelerations of the rigid spacecraft while maintaining the same position. Since the velocities and accelerations are decreased, the required input torques will also be decreased such that the constraints in Eq. (26) are met. To provide a more detailed illustration, we will discuss how to generate and adjust a trajectory in the simulation examples. Because the maneuvering trajectory in Eq. (27) can be randomly generated, a sufficient number of different initial feasible solutions of Problem II can be induced so that all local minima are obtained. Thus, the global minimum is determined by evaluating the objective function at each local minimum.

Procedure to Solve the OLRAC Problem
The procedure to solve the OLRAC problem for a rigid spacecraft between two different rotational attitudes can be summarized as follows.
Algorithm A: (Solution of Problem II) Step 1: Choose a value of Δt, an integer N, and an integer n feasible.
Step 3: Formulate the OLRAC problem as the NLP Problem II with the chosen value of N.
Step 5: Apply the procedure described in Sect. 4 to generate an initial feasible solution of Problem II.
Step 6: Use any NLP algorithm to obtain a local minimum of Problem II based on the initial feasible solution (starting point) specified in Step 5.
Step 7: If i ≤ n feasible , then go to Step 5. Otherwise continue.
Step 8: Choose the smallest local minimum among those found in Step 6.
Step 9: N • Δt is the time of maneuver. End.
In the above algorithm, one finds that (n feasible + 1) different starting points are generated for each value of N and at most (n feasible + 1) local minima can be obtained. Therefore, it is obvious that one can choose a large value of n feasible to generate a sufficient number of starting points. However, this will also result in a long computation time. Therefore, the choice of n feasible involves a trade-off between the computation time and the number of different starting points. As a compromise, when performing the simulation examples in Sect. 6, n feasible is chosen as 20.
A smaller value of Δt will result in higher discretization accuracy. Therefore, an upper limit on the discretization sampling period, denoted by Δt limit , can be chosen by considering the required accuracy of the discretization between the continuous system and the discrete system. If the optimal sampling period Δt(N) corresponding to a given value of Δt is greater than Δt limit , then the value of N needs to be adjusted. Since the choice of Δt limit and the adjustment of N are not the main issues of this paper, the value of N will be fixed in the following simulation examples. For a detailed discussion about the choice of Δt limit and the adjustment of N, see Refs. (21) and (22).

Simulation Results and Discussion
In this section, several simulations are used to verify the feasibility of the proposed method. During each simulation, the K reaction wheels used to actuate the rigid spacecraft are chosen to be the same. The inertia matrix of the spacecraft is assumed to be J Numerical examples are presented for three major kinds of K-wheel reaction control systems. For the case of K = 3, the transformation matrix is assumed to be the third-order identity matrix, which means that three reaction wheels are aligned with the body axes. For K = 4, the transformation matrix is given by In the case of K = 5, the transformation matrix of this wheel configuration is set to 1 0 0 1 2 1 2 In general, the wheel-torque inputs corresponding to the maneuvering trajectory found in Sect. 4 will violate the dynamic constraint in Eq. (26). For example, a plot of the set of infeasible wheel-torque inputs for K = 4, N = 35, and Δt = 16 s is shown in Fig. 2, in which the wheel torque u 3 exceeds the saturation limits in some control steps. A simple way to adjust this trajectory to be feasible is to increase Δt. The plot of a set of feasible wheel-torque inputs corresponding to K = 4, N = 35, and Δt = 20 s is shown in Fig. 3. After generating the initial feasible solutions, the OLRAC problem can be solved by applying Algorithm A. By selecting the weighting factor ρ = 0.01, Δt and N • Δt are found to be 0.6918 and 24.2124 s, respectively, for K = 4 and N = 35. In comparison, for the time of maneuver (t f = N • Δt), the values of 81.6625 and 42.4974 s were given in Refs. (22) and (26), respectively. This indicates the advantageousness of the present method.
For the case of ρ = 0.001 the simulation results of rigid spacecraft attitude maneuvers with K = 3, 4, and 5 reaction wheels are shown in Fig. 4. According to this figure, the rigid spacecraft accomplished the maneuver from the initial attitude to the final attitude in less than 30 s with different values of K. In addition, the greater the number of actuators, the faster the attitude maneuver is achieved. Figure 5 depicts plots of the wheel-torque inputs for K = 5. At least one of the wheel-torque inputs is always at its extreme value, which has been observed in many cases of time-optimal attitude control.    Fig. 6. The time of maneuver approaches infinity as ρ→∞. The minimal time-slewing control is achieved as ρ→0. Meanwhile, a larger integral of the squared sum of wheel-torque magnitudes is required to perform this minimal time-slewing maneuver. Increasing the value of ρ will decrease the control effort and increase the time of maneuver. For ρ ≤ 0.1, increasing the number of reaction wheels K will reduce the time of maneuver but require a larger control effort. On the other hand, the control effort is almost the same with different K for large values of the weighting factor (ρ > 0.1). It is concluded that the trade-off between the time of maneuver and the control effort for a spacecraft's maneuvers can be determined by adjusting the weighting factor ρ. Nearly minimum-time control with an appropriate control effort is obtained by choosing a suitable weighting factor of about ρ ≈ 0.1. This is consistent with the main results provided in Ref. 22. Figure 7 shows the rigid spacecraft angular velocity with respect to the normalized duration t/t f for K = 5 with three different weighting factors. By comparing the results in Fig. 7 with the previous results in Refs. 1, 6, 24, 26, and 27, it can be seen that these angular velocity trajectories with a smaller value of ρ have similar characteristics of the time-optimal solutions. By increasing the value of ρ, smoother spacecraft maneuvers are produced.
In Fig. 8, wheel-torque magnitudes are plotted with respect to the normalized duration t/t f for K = 5 and three different weighting factors ρ with N = 35 control steps. It has been shown that for the case of rest-to-rest slewing with ρ = 10, all the wheel-torque magnitudes are approximately linear functions of time and do not reach their saturation levels. By decreasing the weighting factor ρ, some of the wheel-torque magnitudes can be expected to reach their limits and contribute more to the rotation. By successively decreasing the weighting factor ρ, minimum-time control can be obtained, where some or all of the wheel-torque magnitudes are of the bang-bang type.

Conclusion
We presented a novel method to solve the OLRAC problem of a rigid asymmetric spacecraft actuated by K reaction wheels mounted in arbitrary orientations with respect to the spacecraft frame. The performance index to be minimized is chosen as a weighted performance index of the time of the maneuver and the integral of the squared sum of wheel-torque magnitudes. Different from the conventional shooting method, which utilizes PMP, this optimal control problem is first transformed into an NLP problem by an iterative procedure. The main advantage of the proposed method is that one does not need to solve a set of highly nonlinear differential equations. A systematic method to generate initial feasible solutions is also proposed. This allows the optimization process of the NLP problem to be started from many different points to search for all local minima. Thus, the global minimum can be obtained by evaluating the objective function at each local minimum.
The solutions obtained can be verified to satisfy the Kuhn-Tucker condition, (25,28) which is a criterion used to check a local minimum. If a solution does not satisfy the Kuhn-Tucker condition, then the solution is not a global minimum. However, even if a solution satisfies the Kuhn-Tucker condition, one still cannot determine whether the solution is globally optimal or not. Many phenomena in this paper have indicated that the solutions obtained should be globally optimal. However, more effort will be needed to obtain a mathematical proof.