Time-Optimal Attitude Scheduling of a Spacecraft Equipped with Reaction Wheels

The time-optimal control problem of a spacecraft equipped with reaction wheels has been studied, in which the spacecraft is constrained to sequentially assume a set of attitudes, whose order is not specified. This attitude scheduling problem has been solved as a multiphase mixed-integer optimal control problem in which binary functions have been introduced to model the choice of the optimal sequence of target attitudes and to enforce the constraint of adopting once and only once each attitude. Given the dynamic model of the spacecraft, the initial and final attitudes, and a set of target attitudes, solving this problem consists in finding the control inputs, the sequence of attitudes with the corresponding passage times, and the resulting trajectory of the spacecraft that minimize the time of the maneuver. The multiphase mixed-integer optimal control problem has been converted into a mixed-integer nonlinear programming problem first making the unknown passage times through the target attitudes part of the state, then introducing binary variables to discretize the binary functions, and finally applying a fifth-degree Gauss-Lobatto direct collocation method to tackle the dynamic constraints. The resulting problem has been solved using a nonlinear programming-based branch-and-bound algorithm.


Introduction
This paper focuses on time-optimal attitude scheduling of spacecraft.This problem entails planning a sequence of slew maneuvers that join a set of attitudes of the spacecraft.It is important to point out that the order of these target attitudes to be assumed by the spacecraft is not specified and must be determined.Additionally, some constraints on the angular velocity of the spacecraft at the target attitudes must be satisfied.It is assumed that the spacecraft is rigid, asymmetric, and equipped with reaction wheels, which are subject to bound constraints on their torque and angular velocity.Since the reaction wheels generate torque on the spacecraft by accelerating, the latter assumption implies that they cannot operate when they reach their maximum angular velocity.Moreover, the spacecraft is assumed to be subject to bound constraints on its angular velocities.This kind of problem arises, for instance, in observation scheduling of Earth observation satellites and space telescopes with momentum actuators [1].
The problem of minimum time reorientation of a spacecraft between two attitudes has been studied extensively.See for example [2], in which the spacecraft is assumed to be equipped with reaction wheels, and references therein.In the rest of this section, only previous works on optimal control of spacecraft among several attitudes will be reviewed.
In [3], the Pontryagin's maximum principle has been extended to the multipoint case and applied to the problem of attitude scheduling of a satellite equipped with reaction wheels.In this work, the satellite is constrained to assume a set of attitudes whose order is specified.On the contrary, the passage times through the attitudes of the set are not specified and must be determined.A weighted combination of energy expended and elapsed time has been used as optimality criterion.In the problem formulation, bounds on the control torque have been taken into account.However, both the angular velocity of the satellite and the angular velocity of the reaction wheels are unconstrained.
In [4], the Legendre pseudospectral method has been employed to plan minimum-time slew maneuvers among several attitudes for image acquisition and minimum-time transition maneuvers for scanning operations of a space telescope equipped with reaction wheels.In both cases, the order of the target attitudes is specified.The formulation of the optimal control problem accounts for all of the relevant spacecraft actuator, sensor, and operational constraints, which include the requirement that the proper reaction wheel bias momentum be maintained at the beginning and end of each maneuver, and bounds on the reaction wheel torques, on the reaction wheel momentum, and on spacecraft body rates.The ability of this method to reliably generate flight-implementable shortest-time maneuvers has been demonstrated by means of flight implementation of the technique over multiple operational scenarios.
In [5], a new approach for designing shortest-time maneuvers of a control moment gyroscope Earth observation satellite has been proposed.In this work, the boresight axis of the satellite is required to move among several presorted collection regions and traverse through the center of each of them.A feedback scheme is employed in which, to eliminate the negative effects of unpredictable interactions between the open-loop and the feedback control laws, a Riemann Stieltjes optimal control approach with a tychastic constraint has been used.This constraint ensures that the variations in the control moment gyroscope singularity index caused by the feedback control can be appropriately managed.Experimental results obtained from tests performed on a momentum control system testbed have been reported to show the effectiveness of the method on a real attitude control system.
In this paper, the time-optimal attitude scheduling of a spacecraft is formulated as a multiphase mixed-integer optimal control problem (MIOCP) which can be stated as follows: given the dynamic model of a rigid spacecraft and a set of target attitudes, find the control inputs that steer in minimum time the spacecraft from an initial attitude to a final one adopting all the target attitudes of the set and satisfying some constraints on the angular velocity of the spacecraft at the target attitudes.
In this formulation of the time-optimal spacecraft attitude scheduling problem, the order in which the target attitudes must be assumed by the spacecraft is not specified and must be determined together with the times at which they are assumed.The fact that the order of the target attitudes is not specified characterizes this formulation with respect to previous ones and makes the resulting problem very difficult to solve.
Phases have been considered to model the rotation of the spacecraft between two target attitudes and to impose constraints at these attitudes.These constraints depend on the specific application.In general, they are constraints on the angular velocity of the spacecraft.Binary functions have been introduced to model the choice of the optimal sequence of target attitudes and to enforce the constraint of adopting once and only once each attitude.This problem can be regarded as an extension of the classical traveling salesman problem to dynamic systems, in which the problem consists in finding the inputs of a controlled dynamic system such that the resulting state trajectory visits all states of a given finite set in minimum time.
In [6], a numerical method for MIOCPs has been developed.This algorithm is based on direct multiple shooting method, convexification and relaxation of the original problem, adaptive refinement of the underlying discretization grid, and on both deterministic and heuristic integer methods.In [7], an algorithm for mixed-integer nonlinear model-predictive control has been presented.It is a combination of the direct multiple shooting method, a reformulation based on partial outer convexification and relaxation of the integer controls, a rounding scheme, and a realtime iteration scheme.
The main contribution of this paper is the solution of a spacecraft attitude scheduling problem as a multiphase MIOCP.To the best knowledge of the authors, the present paper is the first to employ and test this technique on this problem.In this approach, the multiphase MIOCP has been converted into a mixed-integer nonlinear programming (MINLP) problem first making the unknown passage times through the target attitudes part of the state, then introducing binary variables to discretize the binary functions, and finally applying a fifth-degree Gauss-Lobatto direct collocation method [8] to tackle the dynamic constraints.High-degree collocation permits the number of variables of the problem to be reduced for a given numerical precision, thus reducing the computational cost needed to solve the resulting MINLP.An important feature of the algorithm employed is a practical MINLP solution technique that combines the power of a stateof-the-art nonlinear programming (NLP) solver with a branch-and-bound (BB) strategy.Note that, even if the order in which the target attitudes are visited is fixed, the problem remains an optimal control problem which is very difficult to solve to global optimality.Our approach is based on the assumption that state-of-the-art direct numerical methods for optimal control problems [9] are able to compute locally optimal solutions that are good approximations of the optimal trajectory.The emphasis of this paper is on the algorithm used to find the discrete part of the solution, that is, the optimal sequence of target attitudes.In the numerical results presented in this paper, a fifth-degree Gauss-Lobatto direct collocation method [8] has been used, which permits the number of the variables of the problem to be reduced for a given numerical precision and, as a consequence, the computation times to be shortened.However, any numerical method for optimal control could be employed.Based on the above assumption, a BB algorithm has been used, which aims at finding the sequence of target attitudes for which the local optimal solution found by the NLP solver is best.In particular, the MINLP solver BONMIN [10] has been used, which includes a BB algorithm based on the NLP solver IPOPT [11].Several numerical experiments have been carried out over different instances of the problem with an increasing number of target attitudes.
The paper is organized as follows.In Section 2, the dynamic model of a rigid asymmetric spacecraft equipped 2 International Journal of Aerospace Engineering with reaction wheels is described.In Section 3, the control properties of this system are reviewed and in Section 4, the spacecraft attitude scheduling problem is stated.In Section 5, the general mathematical formulation of a multiphase MIOCP is given and in Section 6, it is particularized to model the spacecraft attitude scheduling problem.In Section 7, the methodology used for its resolution is outlined and in Section 8, the results of the application of the proposed method to several instances of the spacecraft attitude scheduling problem are reported.Finally, in Section 9, the conclusions are given.

Model of the Spacecraft
The spacecraft is assumed to be actuated by three reaction wheels.For i = 1, 2, 3, assume that the ith reaction wheel is spinning about an axis b i , fixed with respect to the spacecraft, such that the center of mass of the ith wheel lies on it.Assume that a torque T u i is supplied to the ith wheel about the axis b i by a motor fixed with respect to the spacecraft.Consequently, an equal and opposite torque is exerted by the wheel on the spacecraft.Assume that b i is a principal axes of the spacecraft and that it coincides with the principal axis of the ith wheel about which it is symmetric.Axes x, y, and z of the spacecraft reference frame are chosen to be coincident with axes b i , i = 1, 2, 3, respectively.In this paper, the attitude of the spacecraft with respect to the world frame is represented by quaternions [12].In space applications, quaternions are in general organized as a vector q = q 1 , q 2 , q 3 , q 4 T in which the real part is the last element.This convention will be used to derive the dynamic model of the spacecraft.For convenience, in Section 8, together with the quaternion representation, the roll, pitch, and yaw (RPY) representation of attitudes, with angles ϕ, θ, ψ , will be given.
In the rest of this section, [2] will be followed.The attitude kinematic equation in terms of quaternion is where q = q 1 , q 2 , q 3 , q 4 T is the quaternion vector that represents the attitude of the spacecraft, ω = ω 1 , ω 2 , ω 3 T is the angular velocity vector of the spacecraft, and The dynamic model of the spacecraft can be expressed in the following form: where I s is the inertia momentum matrix of the spacecraft, which, since axes b i , i = 1, 2, 3, coincide with axes x, y, and z of the spacecraft reference frame, can be expressed as where I xx , I yy , and I zz are the principal moments of inertia of the spacecraft.I w is the inertia momentum matrix of the reaction wheels T is the angular velocity vector of the reaction wheels, T u = T u 1 , T u 2 , T u 3 T is the vector of torques at the reaction wheels, and T ex = T ex 1 , T ex 2 , T ex 3 T is the vector of external disturbance torques about axes b i , i = 1, 2, 3, respectively.The disturbance torques of T ex are assumed to be of small amount and will be neglected.Thus, the dynamic equation of the spacecraft becomes The dynamic equation of the reaction wheels is The state and control variables of the problem are respectively.
The maximum torque and angular momentum of the reaction wheels are actually control and state constraints, respectively.The constraint on the angular momentum can be transformed into a constraint on the angular velocity of the reaction wheels where Ω is the maximum angular speed of the reaction wheels.The constraints on the maximum torque at the reaction wheels can be expressed as where T u is the maximum torque at the reaction wheels.
The constraint on the angular velocity of the spacecraft can be expressed as where ω is the maximum angular speed of the spacecraft about each control axis.The four components of the quaternion vector must satisfy the following condition: 3 International Journal of Aerospace Engineering Moreover, all the components of quaternion vector are bounded in the interval [−1, 1].
Equations ( 1), (6), and ( 7) can be rewritten in the general form which can be regarded as a differential constraint.The performance index for minimum-time problems is defined as t F − t I , where t I is the initial time, which is usually known, and t F is the unknown final time of the maneuver to be determined.The performance index for minimum energy problems is usually defined as the squared norm of the control vector u.
For numerical reasons, it is useful to scale state and control variables of the optimal control problem as follows: With all previous assumptions, the state equations with scaled state and control variables take the following form: which can be expressed as In the performance index of the minimum-time optimal control problem, without loss of generality, it can be set t I = 0.The scaled state and control variables are subject to the following bound constraints: where

Control Properties of the Spacecraft
In [13] the necessary and sufficient conditions for the controllability of a rigid body in the case of one, two, and three independent control torques are provided.If the spacecraft is controlled by three independent torques, it is completely controllable.Consider an inertially symmetric spacecraft, that is, a spacecraft whose principal moments of inertia are equal.Assume that the control torques for each axis are bounded.The so-called Euler's principal axis or eigenaxis maneuver has often been considered as the fastest rotational maneuver, since it is the shortest angular path between two orientations.Whether the eigenaxis maneuver is optimal or not depends on the definition of the set of admissible control torques.For an inertially symmetric rigid spacecraft, it has been shown in [14] that, when the total magnitude of the control torque is constrained, the eigenaxis maneuver is indeed the time optimal maneuver.It has been shown in [15] that, when the three orthogonal components of the control torque are independently constrained, the associated optimal reorientation time is lower than in the eigenaxis maneuver case, as the nutational components can provide higher control torques about the rotation axis.An example of motion of the body frame of a spacecraft equipped with reaction wheels in a minimum-time restto-rest counterclockwise rotation of π/3 (rad) about the z axis is given in Figure 1.It can be observed from this figure, in which the traces of the x-, y-, and z-axes on the unit sphere are represented in yellow, blue, and purple, respectively, that the deviation from the eigenaxis rotation is large.This deviation does not appear so counterintuitive if one considers the solution to the classical brachistochrone problem in which the shortest path is not the time-optimal path.

The Spacecraft Attitude Scheduling Problem
The optimal control problem studied in this paper can be stated as follows: find the control inputs that steer a spacecraft, whose dynamic model is given by ( 16), from an initial attitude q I to a final attitude q F , through a set of n p target attitudes Q = q 1 , … , q n p , while minimizing the time of the maneuver.More specifically, the spacecraft is constrained to adopt once and only once all the n p attitudes of the set Q while rotating from q I to q F .For the sake of generality, operational constraints at the target attitudes are introduced in the problem formulation.More specifically, the angular velocity of the spacecraft about the z-axis at attitudes q 1 , … , q n p is constrained to be zero.However, the introduction in the problem formulation of other constraints to fulfil specific task requirements is straightforward.
It is important to point out that, since both the order and the times at which the orientations of Q are assumed are not specified, they must be determined.
This problem has been solved as a multiphase MIOCP. 4 International Journal of Aerospace Engineering

General Formulation of the Multiphase MIOCP
In this section, following [6], the general mathematical formulation of the multiphase MIOCP will be given.For the sake of brevity, specific assumptions about the properties of the functions involved in the definition have not been reported.A multiphase MIOCP is a constrained optimization problem of the form: subject to for k = 0, … , n p .In this definition, the variable t ∈ t I , t F represents time, where in general, the initial time t I is fixed, whereas the final time t F may be fixed or free.w t : t I , t F → 0, 1 n w is a binary time-dependent control function and w k t is the value of w t in phase k, that is, for t ∈ t k , t k+1 .If there is a discontinuity in at least one component of the binary control function w t at time t, the function w t switches and time t is called the switching time.Thus, times t k , k = 1, … , n p are the n p switching times between phases, such that As a consequence, for t ∈ t k , t k+1 , the system is in phase k.Γ Ψ is the set of admissible binary control functions w t , which is the set of piecewise constant functions with jumps only at time values contained in the set Ψ.The set Ψ is the feasible switching set, that is, the set of time values at which a discontinuity in the binary control function vector w t may occur.We can assume that Ψ ∈ t I , t F .The vector v ∈ 0, 1 n v is a time-independent vector of binary variables.The vector s ∈ R n s contains all the time-independent variables of the problem, including t F for problems with free final time.
x represents the state variables for phase k.They include both the differential variables y k t : t k , t k+1 → R n k y , that is, the state variables with time derivative, and the algebraic variables , which represents the control variables for phase k.
c represents path constraints.Equation (22e) represents the interior point inequality constraints with r ineq ∈ R n r ineq and t I < t 1 ≤ t 2 ≤ ⋯ ≤ t n r ineq < t F .Equation (22f) represents the interior point equality constraints with r eq ∈ R n r eq and t I < t 1 ≤ t 2 ≤ ⋯ ≤ t n r eq < t F .Interior point state equality and inequality constraints impose point-wise constraints on the state variables in a finite number of grid points.In (22g), the tr function represents the transition and u k t k+1 , and on v and s.
x t I ∈ R n 0 x represents the vector of initial conditions given at the initial time t I and x t F ∈ R n n p x represents the vector of final conditions given at the final time t F .The dimensions n x k , n u k , n w k are not necessarily identical for each phase k.
A comprehensive discussion about optimality conditions and solution methods for the multiphase MIOCP stated above can be found in [6].

Specification of the Multiphase MIOCP
In this section, the general mathematical formulation of the multiphase MIOCP given in Section 5 is particularized to model the spacecraft attitude scheduling problem stated in Section 4.
Phases are introduced in the formulation of the MIOCP only to reflect the rotation of the spacecraft between two different attitudes.Thus, each phase is identified by two attitudes of the spacecraft, the dynamic model of the spacecraft is assumed to be the same across the different phases, and the number of phases is assumed to be known.However, both the untimed sequence of phases and the switching times between them must be determined.
Let p I and p F be the intersection points of the z-axis of the spacecraft with the unit sphere at orientations q I and q F of the spacecraft.Similarly, let P = p 1 , … , p n p be the set of intersection points of the z-axis of the spacecraft with the unit sphere at the attitudes Q = q, … , q n p .Since n p is the number of points of the set P, n p + 1 phases can be identified during the rotation of the spacecraft from t I to t F , with n p switching times between them such that t with p 0 = p I and p n p +1 = p F and Q ′ = q 0 , q 1 , … , q n p , q n p +1 with q 0 = q I and q n p +1 = q F .The objective functional (22a) to be minimized is in this case The dynamic model of the spacecraft (15) takes the following scalar for Equation (26) correspond to the differential constraints (22b) of the problem, which, in this case, are independent of k.
Variable w t is, in this case, the n p + 1-dimensional binary control function that indicates which attitude of the set Q′ is taken by the spacecraft at time t k , k = 0, … , n p + 1.As said before, it is a vector function whose components are piecewise constant in t I , t F with jumps at times t k , k = 0, … , n p + 1, with a single nonzero component in each phase.More precisely, a phase starts when the orientation of spacecraft coincides with q I or when it takes an attitude of the set Q and ends when the spacecraft takes another attitude of the set Q or when it reaches the final attitude q F .Since function w t during phase k is denoted by w k t , then w k j t k = 1, that is, component j of w k t k at time t k equal to 1, means that the spacecraft takes the attitude q j at time t k .Let q k t be the orientation variable of the spacecraft during phase k.Thus, the constraints that impose taking the n p attitudes of the set Q can be expressed as follow: which means that, if w k j t k = 1, it will be q k t k = q j , that is, at time t k the spacecraft will assume the attitude q j and the z-axis will intersect the unit sphere at point p j .Since q 0 t 0 = q 0 = q I and q n p t n p +1 = q n p +1 = qF, the following conditions must be fulfilled: 6 International Journal of Aerospace Engineering Additional constraints to be taken into account are Constraints (29) mean that the spacecraft must assume a single attitude of the set Q at time t k , k = 1, … , n p , whereas constraints (30) mean that the spacecraft must assume each attitude q j , j = 1, … , n p , once.Equations ( 27), (29), and (30) are algebraic constraints of type (22c).Another algebraic constraint to be considered is (12), which is independent of k.Equations ( 19) and ( 20) are path constraints of type (22d), which are independent of k.
The angular velocity of the spacecraft about the z-axis at attitudes of the set Q is constrained to be zero.These are interior point constraints of type (22f) and can be expressed as In the present case, interior point constraints (22f) are given only for one component of the state vector, ω 3 t k , whereas boundary conditions (22j) are given for all the components of the state vector at t I and t F .Actual boundary constraints are specified in the description of each experiment.
Transition conditions (22g) are usually introduced to enforce continuity and therefore are referred to as linkage conditions.In the reformulation of the multiphase MIOCP as a MINLP presented in this paper, a special technique is used to avoid the necessity to introduce the linkage conditions explicitly.This technique will be described in the next section.

Numerical Solution of the Multiphase MIOCP
In this section, reformulation of the multiphase MIOCP as a MINLP problem and the numerical method used to solve it will be briefly outlined.
The multiphase MIOCP has been tackled first making the unknown switching times part of the state and then introducing a new independent variable with respect to which the switching times are fixed.In this reformulated problem, there is a linear relation between the new variable and time, but the slope of this linear relation changes in each interval between any two switches.These slopes, which are part of the solution to the multiphase MIOCP, are actually time-scaling factors that determine the optimal switching times.Thus, the reformulated problem takes the form of a MIOCP in which linkage conditions are not necessary.More details on this technique can be found in [16], Chapter 6.
The MIOCP has been converted into a MINLP applying a fifth-degree Gauss-Lobatto direct collocation method to transcribe the dynamic constraints [8].This transcription only involves the continuous part of the MIOCP.The values of the binary function at the switching times are taken as binary variables of the MIOCP.Thus, a MINLP is obtained which consists in minimizing a nonlinear function subject to nonlinear constraints in a space where some of the variables take values in the set {0, 1}, while the others take values in R.
Notice that fixing the binary variable w k j , where k = 0, … , n p and j = 0, … , n p + 1, is equivalent to fixing the sequence of target attitudes of the spacecraft and, if this is done, the MIOCP takes the form of a continuous optimal control problem.Therefore, the corresponding MINLP takes the form of a NLP problem.Thus, the method employed in this paper is based on the observation that, when the sequence of target attitudes is fixed, NLP algorithms, which compute locally optimal solution to transcribed continuous optimal control problems, find good trajectories in short computational times.
A simple algorithmic approach could therefore be to enumerate all possible values for w k j , solve the corresponding continuous optimal control problems, and pick the best solution.However, this approach is impractical because the number of problems to solve would grow exponentially with the number of target attitudes and would be too large even with a handful of them.To deal with these difficulties, the NLP BB framework BONMIN [10] based on the interior point NLP solver IPOPT [11] has been used.

Numerical Results
In this section, the results of the application of the numerical method described in Section 7 to several instances of the spacecraft observation scheduling problem stated in Section 4 will be reported.
In all the numerical experiments, the following initial and final attitude parameters of the spacecraft have been assumed:

32
which correspond to the following RPY angles: ϕ t I = θ t I = ψ t I = θ t F = ψ t F = ϕ t F = 0 rad .The initial and final angular velocities of both the spacecraft and the reaction wheels have been set to 0. This maneuver can be regarded as a rest-to-rest maneuver in which the spacecraft is required to rotate from the initial attitude to the final one adopting several target attitudes.At all the target attitudes, the angular velocity of the spacecraft about the z-axis has been constrained to be zero.As said before, although in our experiments only velocity constraints on the angular velocity of the body reference frame about the z-axis at the target attitudes have been considered, the introduction in the problem formulation of other constraints to fulfil specific task requirements is straightforward.This can be 7 International Journal of Aerospace Engineering done by adding equality and inequality constraints on the angular velocity of the spacecraft at the time at which it reaches a target attitude.
The target attitudes used in the experiments are listed in Table 1, parameterized using both unit quaternion and RPY angles.Three groups of target attitudes can be distinguished depending on the range of the angle ϕ.For the attitudes a-f, this range is 15 The number of intervals of the discretization of the fifthdegree Gauss-Lobatto direct collocation method has been 2 for each phase.From the results of the error analysis based on the number of intervals of the discretization reported in [8], it can be deduced that this value is a good trade-off between good numerical accuracy and reduced number of variables.Since in an experiment with n p target attitudes, the maneuver has to be subdivided in n p + 1 phases; a total number of 2 n p + 1 intervals are obtained.
To check the numerical precision of the results, a comparison has been carried out using the software DIDO (http://www.elissarglobal.com/),a commercial MATLAB (https://www.mathworks.com/)optimal control toolbox which implements the Legendre-Gauss-Lobatto pseudospectral method, a widely used optimal control technique in aerospace engineering [17].The results of the comparison will be discussed at the end of this section.
Thus, the total number of intervals used in the discretization has been 14, 26, and 38 in the experiments with 6, 12, and 18 target attitudes, respectively.In the figures, 5 points for each phase are represented, since for each phase, state and control variables are evaluated at the independent collocation points, namely, at the endpoints and midpoints of the grid.The numerical experiments have been carried out on a 4 GHz i7 processor with 16 GB RAM.
8.1.Experiment A with 6 Target Attitudes.In this experiment, the spacecraft is constrained to assume the 6 target attitudes of the set Q A = q a , q b , q c , q d , q e , q f while moving from the initial attitude, q I , to the final one, q F .The computation time to find a solution was 5.48 [s]; the minimum time to execute the maneuver is 592.49[s]; and the optimal sequence of target attitudes is q a , q f , q e , q d , q c , q b 33 The corresponding motion of the body frame of the spacecraft is represented in Figure 2. In this figure, the trace of intersection points of the z-axis of the spacecraft with the unit sphere at the attitudes of the set Q A .For the sake of clarity, the traces of the x-and y-axes are represented in grey.The state variables of the optimal solution are represented in Figures 3-5, and the control variables are represented in Figure 6.In these figures, phases are represented by vertical grey dashed lines located at the optimal switching times.The numerical values of the optimal switching times between phases are listed in Table 2.They correspond to the times at which the spacecraft assumes the attitudes of the set Q A .It can be seen from Figure 6 that the control is of bang-bang type and from Figure 5 that all the reaction wheels saturate during the maneuver.It can be observed from Figure 4 that the spacecraft does not reach its maximum angular velocity during the maneuver and from Figure 4(c) that the angular velocity of the spacecraft about the z-axis at the target attitudes is zero as required.It can be noticed from Figure 5(c) that the angular velocity of reaction wheel 3 at the target attitudes is zero.
Since, in this case, the set of target attitudes contains only 6 elements, the number of their permutations is not prohibitive for an exhaustive search of the best solution.The exhaustive check of the 720 possible sequences of target attitudes confirmed that the sequence q a , q f , q e , q d , q c , q b 35 gives rise to the global optimal solution.The second best solution is obtained for the sequence q a , q b , q c , q d , q e , q f , 36 which corresponds to a maneuver whose minimal duration is 628.46 [s] and the third best solution is obtained for the sequence q a , q c , q d , q f , q e , q b , 37 which corresponds to a maneuver whose minimal duration is 647.22 [s].
8.2.Experiment B with 12 Target Attitudes.In this experiment, the spacecraft is constrained to assume the 12 target attitudes of the set Q B = q a , q b , q c , q d , q e , q f , q g , q h , q i , q j , q k , q l 38 while moving from the initial attitude, q I , to the final one, q F .The computation time to find a solution was 23.27 [s].The minimum time to execute the maneuver is 744.03 [s].The optimal sequence of target attitudes is q a , q b , q h , q c , q d , q e , q f , q l , q k , q j , q i , q g 39 The corresponding motion of the body reference frame of the spacecraft is shown in Figure 7.In this figure, the trace of the z-axis on the unit sphere is represented in purple together with the set P B = p a , p b , p c , p d , p e , p f , p g , p h , p i , p j , p k , p l 40 of intersection points of the z-axis of the spacecraft with the unit sphere at the attitudes of the set Q B .For the sake of clarity, the traces of the x-and y-axes are represented in grey.The switching times between phases are listed in Table 3. Since, in this case, the set of target attitudes contains 12 elements, the number of their permutations is prohibitive for an exhaustive search of the best solution as done in experiment A with 6 target attitudes.8.3.Experiment C with 18 Target Attitudes.In this experiment, the spacecraft is constrained to assume the 18 target attitudes of the set Q C = q a , q b , q c , q d , q e , q f , q g , q h , q i , q j , q k , q l , q m , q n , q o , q p , q q , q r 41 while moving from the initial attitude, q I , to the final one, q F .The computation time to find a solution was 35.50 [s]; the minimum time to execute the maneuver is t F = 897 46 [s]; and the optimal sequence of target attitudes is q g , q m , q n , q o , q p , q k , q e , q f , q l , q r , q q , q d , q j , q i , q c , q h , q b , q a 42 11 International Journal of Aerospace Engineering of intersection points of the z-axis of the spacecraft with the unit sphere at the attitudes of the set Q C .For the sake of clarity, the traces of the x-and y-axes are represented in grey.The switching times between phases are listed in Table 4.
It is worth mentioning that in Figures 2, 7, and 8, only the intersection points between the z-axis of the body frame of the spacecraft with the unit sphere at the target attitudes have been represented.As a consequence, these points represent only the spatial position of the z-axis of the body frame, whereas the spatial position of the other two axes of the body frame is not represented.For this reason, the trace of this intersection point in these figures should not be considered as unnecessarily erratic.
For the sake of comparison, the minimum times to execute the maneuvers and the computation times to find the solutions in all the numerical experiments are reported in Table 5.In spite of the fact that the number of permutations of the set of target attitudes is 720 in experiment A and 12! = 479,001,600 in experiment B, the computation time in experiment B is only 4.24 times higher than the value obtained in experiment A. Similarly, although the number of permutation of the set of target attitudes in experiment C is 18! = 6.4023737057 × 10 15 , the computation time in experiment C is only 6.47 times higher than the value obtained in experiment A. It is worth observing that 18! is a huge number.Indeed, supposing that computing a minimum-time solution for a given sequence of target attitudes takes about 0.1 [s], an exhaustive search of the optimal sequence would take more than 20,301,793 years.This shows the effectiveness of the BB strategy used to find the optimal solution.
Although the number of target attitudes in experiment B is twice the number of target attitudes in experiment A, the obtained value of the minimum time to execute the maneuver in experiment B is only 1.26 times higher than the value obtained in experiment A. The same occurs in experiment C, in which the number of target attitudes is three times the number of target attitudes in experiment A, and the obtained value of the minimum time to execute the maneuver is only 1.51 times higher than the value obtained in experiment A. This is due to the fact that the optimal maneuver obtained in experiment A itself exhibits large deviations from the eigenaxis rotation when rotating from one target attitude to another.
For the cases considered herein, the angular velocities of the spacecraft about the three axes are limited by the bound    constraints on the angular velocity of the reaction wheels and therefore do not reach their maximum value.As said before, to check the numerical precision of the results, a comparison has been carried out using the software DIDO.Since this MATLAB optimal control toolbox does not support mixed-integer formulations of the optimal control problem, a test has been carried out over a minimum-time rest-to-rest maneuver of the spacecraft between two target attitudes, namely, q a and q b , given in Table 1.The attitude parameters q of the solution are shown in Figure 9.In this figure, the solution obtained with the fifth-degree Gauss-Lobatto direct collocation method is represented with continuous lines, whereas the solution obtained with the pseudospectral method is represented with dashed lines.The optimal maneuver time has been 57.0503[s] for the fifth-degree Gauss-Lobatto direct collocation method with 2 intervals and 57.58834 [s] for the pseudospectral method with 50 intervals.Assuming that the solution obtained with the pseudospectral method is reliable, we can conclude that, despite the significant difference in the number of intervals of the discretization, the obtained solutions are very similar.The main difference is in the computation times.Indeed, the computation time was 0.585451 [s] for the fifth-degree Gauss-Lobatto direct collocation method and 59.638876 [s] for the pseudospectral method.This confirms that using a high-degree direct collocation method, notable reductions of the computation time can be achieved without significantly reducing the accuracy of the solutions.

Conclusions
One of the main difficulties of the spacecraft attitude scheduling problem studied in this paper, in which the spacecraft is constrained to assume a set of attitudes, is that the order in which they must be adopted is not predefined and must be determined.This combinatorial problem has been formulated as a multiphase mixed-integer optimal control problem.
A key feature of the method used for its resolution is the use of a fifth-degree direct collocation method to tackle the dynamic constraints, which permits the number of variables of the problem to be reduced for a given numerical precision.
Another key feature is the use of a branch-and-bound strategy to deal with the combinatorial aspects of the problem.

Nomenclature b i :
Principal axes of the spacecraft, rotation axes of the reaction wheels, axes of the body reference frame q = q 1 , q 2 , q 3 , q Figure 9: Attitude parameters of a minimum-time rest-to-rest maneuver of the spacecraft between the target attitudes q a and q b given in Table 1.The solution obtained with the fifth-degree Gauss-Lobatto direct collocation method is represented with continuous lines whereas the solution obtained with the pseudospectral method is represented with dashed lines.

Figure 1 :
Figure 1: Motion of the body frame of the spacecraft in the solution of a minimum-time rest-to-rest counterclockwise rotation of π/3 (rad) about the z-axis.The traces of the x-, y-, and z-axes on the unit sphere are represented in yellow, blue, and purple, respectively.
-25 [deg]; for attitudes g-l, this range is 40-50 [deg]; and for the attitudes m-r, the range of the angle ϕ is 85-95 [deg].For all attitudes, the ψ angle takes the same set of values in the range 10-290 [deg] and the θ angle is zero.In this way, several numerical experiments have been designed with a different number of target attitudes and involving small, medium, and large angle rotations in ϕ or a combination of them.In all the numerical experiments, the following parameters have been used for the dynamic system: I xx = 6 63 Kg m 2 , I yy = 8 90 Kg m 2 , I zz = 9 63 Kg m 2 , I w 1 = I w 2 = I w 3 = 0 000169 Kg m 2 , Ω = 710 rad/s , T u = 0 0101 Nm , and ω = 0 04 rad/s .

Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :
Figure 2: Motion of the body frame of the spacecraft in experiment A with 6 target attitudes.The traces of the x-and y-axes on the unit sphere are represented in grey, whereas the trace of the z-axis is represented in purple.

Figure 7 :
Figure 7: Motion of the body frame of the spacecraft in experiment B with 12 target attitudes.The traces of the x-and y-axes on the unit sphere are represented in grey, whereas the trace of the z-axis is represented in purple.

Figure 8 :
Figure 8: Motion of the body frame of the spacecraft in experiment C with 18 target attitudes.The traces of the x-and y-axes on the unit sphere are represented in grey, whereas the trace of the z-axis is represented in purple.

Table 1 :
Target attitudes used in the experiments.
of the z-axis on the unit sphere is represented in purple together with the set P A = p a , p b , p c , p d , p e , p f 34

Table 2 :
Switching times expressed in [s] of the optimal solution of experiment A with 6 target attitudes.

Table 3 :
Switching times expressed in [s] of the optimal solution of experiment B with 12 target attitudes.

Table 4 :
Switching times expressed in [s] of the optimal solution of experiment C with 18 target attitudes.

Table 5 :
Minimum times and computation times in [s] obtained in the numerical experiments for different numbers of target attitudes.