Abstract
Previous work on time-optimal satellite slewing maneuvers, with one satellite axis (sensor axis) required to obey multiple path constraints (exclusion from keep-out cones centered on high-intensity astronomical sources) reveals complex motions with no part of the trajectory touching the constraint boundaries (boundary points) or lying along a finite arc of the constraint boundary (boundary arcs). This paper examines four cases in which the sensor axis is either forced to follow a boundary arc, or has initial and final directions that lie on the constraint boundary. Numerical solutions, generated via a Legendre pseudospectral method, show that the forced boundary arcs are suboptimal. Precession created by the control torques, moving the sensor axis away from the constraint boundary, results in faster slewing maneuvers. A two-stage process is proposed for generating optimal solutions in less time, an important consideration for eventual onboard implementation.
1. Introduction
The problem of reorienting a spacecraft in minimum time, often through large angles (so-called slew maneuvers) and subject to various constraints, can take a number of forms. For example, the axis normal to the solar panels may be required to lie always within some specified minimum angular distance from the sun-line. In cases where the control authority is low and the slew requires a relatively long time, certain faces of the vehicle may benefit from being kept as far as possible from the sun-line to avoid excessive solar heating. For many scientific missions, observational instruments must be kept beyond a specified minimum angular distance from high-intensity light sources to prevent damage.
Before addressing the time-optimal, constrained problem, it is useful to review what is known about the unconstrained problem. In a seminal paper, Bilimoria and Wie [1] consider time-optimal slews for a rigid spacecraft whose mass distribution is spherically symmetric, and which has equal and independently limited control-torque authority for all three axes. Despite the symmetry of the system, the intuitively obvious time-optimal solution is not a λ-rotation about the eigenaxis (the appendix contains a discussion of λ-rotations and the eigenaxis of a direction cosine matrix). Indeed, the fallacy here is that one easily confuses the minimum-angle of rotation problem (i.e., the angle about the eigenaxis) with the minimum-time problem, ignoring the constraints imposed by Euler’s equations of rigid-body motion. Bilimoria and Wie find that the time-optimal solution includes precessional motion to achieve a lower time (approximately 10% less) than that obtained with an eigenaxis maneuver. Further, they determine that the control history is bang-bang, with a switching structure that changes depending upon the magnitude of the angular maneuver (referring here to the final orientation in terms of a fictitious λ-rotation, with associated angle ): for values of less than 72 degrees, the control history is found to contain seven switches between directions of the control torque components; larger values of require only five switches.
Several subsequent papers have revisited the unconstrained problem, including such modifications as axisymmetric mass distribution and only two-axis control [2], asymmetric mass distribution [3], small reorientation angles [4], and combined time and fuel optimization [5]. Recently, Bai and Junkins [6] have reconsidered the original problem (spherically symmetric mass distribution, three equal control torques) and find that at least two locally optimum solutions exist for reorientations of less than 72 deg. (one of which requires only six switches) if the controls are independently limited. Further, they prove that if the total control vector is constrained to have a maximum magnitude (i.e., with the orthogonal control components not independent), then the time-optimal solution is the eigenaxis maneuver.
Hablani [7] and Mengali and Quarta [8] consider constrained maneuvers, but focus upon generating feasible solutions without attempting to find optimal solutions. Melton [9] considers time-optimal, constrained slewing maneuvers for cases involving multiple path constraints. That work uses the Swift spacecraft [10] as an example of a vehicle that must be rapidly reoriented to align two telescopes at a desired astronomical target, namely, a gamma-ray burst. The satellite’s burst alert telescope (wide field of view) first detects the gamma-ray burst and the spacecraft then must reorient to allow the X-ray and UV/optical instruments to capture the rapidly fading afterglow of the event. To prevent damage to these instruments, the slewing motion is constrained to prevent the telescopes’ common axis from entering established “keep-out” zones, defined as cones with central axes pointing to the Sun, Earth, and Moon, with specified half-angles. A somewhat surprising result is that all of the cases studied yield trajectories of the sensor axis that neither travel along the boundary of the keep-out cone nor even touch it (so-called boundary arcs and boundary points). Figure 1 shows an example of this behavior, in which the sensor axis traverses a narrow gap (0.1 deg.) between the Sun and Moon cones, but does not intersect either cone.
This paper presents a preliminary study of boundary arcs and boundary points in this same problem. A full analytical solution is not possible; however, some insight to the problem can be gained by examining instances where the sensor axis is constrained to follow the constraint boundary, and those where the initial and final sensor axis directions lie exactly on the boundary. The paper also addresses the practical challenge of implementing onboard optimal control for this type of maneuver and proposes a means for reducing computation time for the reference trajectory.
2. Problem Statement
The problem is formulated as a Mayer optimal control problem, with performance index where is the final time to be minimized. Euler’s equations of rigid-body motion describe the system dynamics Note that the control torques are assumed to be independently bounded, and no assumption is made about the type of control actuator being used. Euler’s equations must be augmented with kinematic relationships in order to determine the orientation of the body over time. In this work, a formulation that uses Euler parameters is employed, with the relationship between the Euler parameters and the angular velocity components given by
The appendix discusses the relationship between the Euler parameters, λ-rotations, and direction cosine matrices. The problem to be considered is a rest-to-rest maneuver, with initial conditions and two possible sets of final conditions for which the final orientation at the final time is completely specified, or for which only some aspect of the final orientation is specified (this is discussed further in Section 2.1).
For the unconstrained optimal control problem formulated using (2.1)–(2.3), the Hamiltonian is For spacecraft of the type being modeled here (Swift, or other astronomical missions), one or more sensors are fixed to the spacecraft bus; these sensors all have the same central axis for their fields of view and this axis is designated here with the unit vector and referred to as the sensor axis. This axis must be kept at least a minimum angular distance from each of several high-intensity astronomical sources. Denoting the directions to these sources as , where the subscript can be S (Sun), E (Earth), or M (Moon), the so-called keep-out constraints are then written as follows: Without loss of generality, is assumed to lie along the body-fixed -axis and its orientation with respect to the inertial frame is then determined from (A.6), with as the inertial frame and as the body-fixed frame It is further assumed that the reorientation maneuver can occur quickly enough that the spacecraft’s orbital position remains essentially unchanged, and that therefore, the inertial directions to the high-intensity sources also remain constant during the slew maneuver.
Analytically, the path constraint given by (2.6) can be adjoined to the Hamiltonian by creating the following constraint function: then substituting for using (2.7). The result must then be differentiated twice with respect to time (and using (2.2) and (2.3)) in order to get a form in which the control torques appear [11]. Finally, a new Hamiltonian is formed by adjoining to with Lagrange multiplier , with the conditions Further, the so-called tangency conditions [11] must also be applied (i.e., for any part of the trajectory on the boundary, not only , but also its first and second time derivatives must be zero). The resulting form is analytically intractable for assessing whether necessary conditions are met along a boundary arc or at a boundary point; however, by using a direct method (Legendre pseudospectral), it is possible to obtain numerical solutions that meet the necessary conditions of optimality. Fleming and Ross [12] show that a pseudospectral method is effective for solving the unconstrained minimum-time reorientation problem.
2.1. Semifree Final Orientation
For some missions, the reorientation strategy may be altered if the final orientation is not completely specified. An example would be the need to reorient the sensor axis to a desired target direction in minimum time, with no constraints on the orientation of the other body-fixed axes at the final time. In practice, some subsequent rotation about the sensor axis might be required to optimize some other parameter (e.g., maximizing illumination of solar panels, or minimizing solar heating of sensitive components), but the principal reorientation maneuver could be achieved faster. The corresponding optimal control problem is the same as before, but with the final conditions on the Euler parameter values given in (2.4c) specified by
2.2. Constrained Rotation Axis
Consider now the suboptimal approach of a simple λ-rotation that carries the sensor axis from initial to final orientation along the constraint boundary, with the other body-fixed axes undergoing the same λ-rotation This amounts to a λ-rotation about the keep-out cone axis (Figure 2). Such a constrained rotation is a simple one-dimensional problem whose solution is bang-bang, with the maximum torque being applied along the direction of . Because the control torque components are each limited to a maximum magnitude , the maximum torque along has magnitude where . Equation (2.12) thus maximizes while obeying the limits on the individual control torques , , and . Assuming that a feasible rotation axis can be identified, then the rotation angle can be determined from (A.5), and the rotation time is given by where is the moment of inertia about the -axis. This provides a useful benchmark to which the optimal solution can be compared. As a practical matter, such a forced-axis rotation could be an acceptable suboptimal solution if the optimal path required too much computing time or resources to calculate.
3. Results
Four cases are considered: two of these (cases BA1 and BA2) constrain the sensor axis to move along the constraint boundary; the initial and final sensor axis directions also lie on the boundary. These cases serve as proxies for boundary arcs, at least in that they provide some indication of the slewing time and other qualitative properties of the motion. The other spacecraft axes are unconstrained at the final time. The other two cases (BP1 and BP2) have the sensor axis beginning and terminating on the constraint boundary (and of course, prohibited from entering the constraint cone); the other spacecraft axes are unconstrained at the final time.
The numerical results are generated using a Legendre pseudospectral method, implemented in the software package DIDO [13]. In all cases, the discrete solution for the control history produced by the pseudospectral method has been subsequently employed to numerically propagate the dynamics from the given initial conditions; the results give solutions that match the discrete solutions to within an absolute error of 10−16 (corresponding units) at each node. These cases all require significant computation time (as much as 72 hours on a computer with an Intel Core 2 2.0 GHz processor, with the number of pseudospectral nodes in the range of 100–250.
Note that the Hamiltonian and costate values are not evaluated directly during the pseudospectral solution process, but rather are reconstructed via the covector mapping principle [14] after the problem solution is found. In all cases here, the Hamiltonian is found to be reasonably constant given the relatively small number of nodes and level of accuracy specified for the nonlinear programming aspect of the calculations.
A system of nondimensional units is employed, partly to provide somewhat more general results, but chiefly because this system of units provides the kind of scaling needed for the pseudospectral method to perform well. In physical units, the angular velocities, moments of inertia, and control torques about the spacecraft’s principal axes are denoted as , , and , respectively; the corresponding nondimensional quantities are defined as where is the time unit, and and are the maximum values of the principal inertias and control torques, respectively. In all of the cases presented, the spacecraft is assumed to have a spherically symmetric mass distribution (i.e., all three principal moments of inertia are equal) and three-axis control capability, with equal and independently limited control authority about all three axes.
3.1. Case BA1
This is a forced boundary arc trajectory of the sensor axis . The constraint cone has half-angle of 45 deg. (corresponding approximately to the Sun keep-out cone for the Swift spacecraft) and the sensor axis must traverse an arc of −90 deg. about the cone’s central axis. This solution uses 151 nodes in the pseudospectral method (even a modest increase in the number of nodes (to 171) resulted in nonconvergence with the available computing resources).
A notable feature of the motion (Figure 3) is its qualitative similarity to the unconstrained time-optimal solution (shown in Figure 4). Referring to the angular velocity components, it is evident that the motion includes precession, with portions of the trajectory including fully saturated control torques. Some finite time intervals have intermediate torques, necessitated by the boundary arc constraint. This solution yields a final time of (nondimensional units). For comparison, if the rotation axis is constrained to lie along the constraint cone’s axis, the final time would be .
(a)
(b)
(c)
(a)
(b)
(c)
Figure 5 depicts the Hamiltonian and costate histories. It is evident that the Hamiltonian is fairly constant, giving some confidence that the optimal solution has been determined.
(a)
(b)
3.2. Case BA2
For this forced boundary arc, the constraint cone has half-angle = 23 deg. (corresponding to the Moon cone for Swift), and the sensor axis must traverse an angle of −70 deg. about the cone’s central axis. The final orientation of the sensor axis is calculated via (A.1), with , , deg., and .
For this case, the solution uses 100 nodes and has a final time of . For comparison, if the rotation axis is constrained to lie along the constraint cone’s axis, the final time would be . Figure 6 depicts the dynamic response and control torques; Figure 7 shows the Hamiltonian and costate histories. The motion and controls are qualitatively similar to those in case BA1. It should be noted that in both BA1 and BA2, the path of the sensor axis is verified to lie within 10−16 radians of the constraint boundary.
(a)
(b)
(c)
(a)
(b)
3.3. Case BP1
This case is identical to case BA1, except that only the initial and final directions of the sensor axis are constrained to lie on the constraint boundary, corresponding to two forced boundary points. Two of the control torques ( and ) exhibit bang-bang behavior whereas shows some chatter (even with 250 nodes employed), as seen in Figure 8. The somewhat larger variation in the Hamiltonian (Figure 9) occurs near the initial time; however, this happens at a finite time (corresponding to 22 pseudospectral nodes after ) after the sensor axis has moved away from the constraint boundary (Figure 10). Such behavior is therefore not the theoretically expected discontinuity in the Hamiltonian and costates at a point where the trajectory leaves the constraint boundary [11]; indeed, such discontinuities may not be observable in numerical solutions where the equality condition in (2.10) is unlikely to occur. The numerical solutions may be improved by using a Bellman chain, based upon a sequence of embedded optimal solutions each of which can use a relatively small number of nodes [15].
(a)
(b)
(c)
(a)
(b)
Nevertheless, the trend is clear: precession created by the control torques, works to reduce the final time to 1.9258, approximately 1% faster than the solution in BA1.
3.4. Case BP2
This case is identical to case BA2, except that only the initial and final directions of the sensor axis are constrained to lie on the constraint boundary. Unlike case BP1, the control torques display more intermediate behavior (Figure 11); solution accuracy is not as good here since convergence could not be obtained for more than 100 nodes in the 72 hours of computer time available. This is also evident in the slight variation in the Hamiltonian (Figure 12). As with case BP1, the only points where the sensor axis contacts the constraint boundary (see Figure 13) are at the initial and final times. The final time achieved is .
(a)
(b)
(c)
(a)
(b)
4. Practical Consideration
Regardless of whether boundary points or boundary arcs exist as part of an optimal trajectory, the numerical determination of the solution via a pseudospectral method frequently requires considerable computation time; indeed, the actual slewing maneuver of the spacecraft can be accomplished in seconds or minutes (depending upon the control authority) while the pseudospectral solver requires from 20 minutes to 72 hours to compute the solution. Experience shows that providing even a rudimentary estimate of the states and controls as an initial guess for the pseudospectral solver can reduce the computation time significantly. A two-stage process is proposed, wherein a random-process algorithm, such as a particle swarm optimizer (PSO) which can rapidly explore the solution space, provides the initial guess to the pseudospectral algorithm. The literature abounds with hybrid methods used in related control problems. For example, Ahmed et al. [16] successfully apply just a PSO to the problem of tuning a satellite’s attitude controller while Sentinella and Casalino [17] examine a hybrid evolutionary algorithm that comprises differential evolution, genetic algorithms, and a PSO applied to the problem of spacecraft trajectory optimization.
Initial efforts that employ a two-stage process for optimal slewing maneuvers show promising results. A PSO produces a low-quality approximation of the states, controls, and node times for the solution, followed by the pseudospectral solver, which takes the approximate solution as its initial starting point. An efficient method for generating the first-stage solution is to represent each control torque component as a sum of Chebyshev polynomials of the form [18] where is the Chebyshev polynomial of degree , and the coefficients are determined by the PSO using an explicit integration of the equations of motion ((2.2), (2.3)). Implementation of (4.1) makes use of the Clenshaw recursion relation [18] for rapid evaluation of the Chebyshev polynomials.
As an example, the control torque for a one-dimensional slewing maneuver (with no keep-out cones) is shown in Figure 14(a). The PSO runs for only 50 iterations (requiring 64 seconds of computation time), producing a crude approximation to the bang-bang control solution that is known to be optimal. That approximate solution is then employed by the DIDO pseudospectral solver to calculate the optimal solution (Figure 14(b)), requiring 12 seconds. Using this two-stage method, the total computation time (PSO plus pseudospectral solver) requires approximately half the time needed using only the pseudospectral method (148 seconds) with no initial guess for the solution. Future study should examine the utility of the two-stage method for the full three-dimensional, constrained problem.
(a)
(b)
5. Conclusion
This preliminary study indicates that, although boundary arcs and boundary points may exist in time-optimal spacecraft slewing maneuvers with path constraints, they are at best part of a suboptimal solution. The numerical calculations (completed via a Legendre pseudospectral method) show that even if the initial and final states are boundary points, the solution moves away from the constraint boundary, resulting in a lower final time than if the motion is forced to move along the boundary (a forced boundary arc). The necessary conditions lead to an unwieldy set of relations, making it impossible to determine analytically if boundary points or boundary arcs are excluded. Further examination of the problem using a Bellman chain to improve the numerical accuracy may provide additional insight. A two-stage method for generating optimal solutions in less time than that required by the pseudospectral method alone shows some promise, but further work is needed to determine its utility for the three-dimensional constrained problem.
Appendix
λ-Rotations and the Eigenaxis
Consider a dextral orthonormal basis set , fixed in reference frame . A copy of this set, denoted , is rotated in a right-handed sense, through an angle , about a unit vector , which has fixed orientation with respect to . The new orientation is denoted as reference frame and the rotation axis will have the same orientation with respect to as it has to . This representation of the orientation of with respect to is called a λ-rotation [19]. For two vectors and , fixed in frames and , respectively, if initially, and undergoes a λ-rotation, then the subsequent relationship between and is given by a result named after Olinde Rodrigues (Goldstein [20] claims that the relationship was known before Rodrigues, and that the vector form used here was first published by Gibbs [21]). This form is useful in describing rotations of a sensor axis about an axis fixed in both the inertial and spacecraft frames.
If one uses a direction cosine matrix to express the orientation of with respect to , then the constituent basis vectors are related by and the elements of are given by
Every direction cosine matrix has one unity-valued eigenvalue with corresponding eigenvector (in matrix form) , that is, and the elements of clearly constitute the components of the associated λ-rotation vector. The vector is commonly referred to as the eigenaxis for the rotation that generates from .
By taking dot products of both sides of (A.1) with and then with , and using appropriate vector identities, it can be shown that which permits the calculation of for given , , and . Another useful relationship gives the direction cosine matrix in terms of the Euler parameters
Nomenclature
| : | Constraint function |
| : | Hamiltonian without the path constraint |
| : | Hamiltonian with the path constraint |
| : | Principal moment of inertia (nondimensional) |
| : | Principal moment of inertia (in physical units) |
| : | Performance index |
| : | Control torque (nondimensional) |
| : | Control torque (in physical units) |
| : | Final time |
| : | Half-angle of the keep-out cone for object |
| : | Euler parameter |
| : | Rotation axis |
| : | Costate corresponding to state |
| : | Lagrange multiplier associated with constraint function |
| : | Rotation angle about the -axis |
| : | Time unit |
| : | Direction of the sensor axis |
| : | Direction of the central axis of the keep-out cone for object |
| : | Angular velocity component (nondimensional) |
| : | Angular velocity component (in physical units). |
Acknowledgment
This paper has been presented at the 6th International Workshop and Advanced School “Spaceflight Dynamics and Control,” Covilha, Portugal, March 28–30, 2011, http://www.aerospace.ubi.pt/workshop2011/.