# Numerical Analysis of Constrained, Time-Optimal Satellite Reorientation

Authors: Robert G. Melton
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 1. Introduction

The problem of reorienting a spacecraft in minimum time, often through large angles (so-calledslew 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.Figure 1
Trajectory of sensor axis between Sun (yellow) and Moon (gray) cones.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(2.1)J=tf,
where tf is the final time to be minimized. Euler’s equations of rigid-body motion describe the system dynamics(2.2)ω̇1=[M1-ω2ω3(I3-I2)]I1ω̇2=[M2-ω3ω1(I1-I3)]I2ω̇3=[M3-ω1ω2(I2-I1)]I3.
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 εi and the angular velocity components ωj given by(2.3)[ε̇1ε̇2ε̇3ε̇4]=12[ε4-ε3ε2ε1ε3ε4-ε1ε2-ε2ε1ε4ε3-ε1-ε2-ε3ε4][ω1ω2ω30].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(2.4a)ω1(0)=ω2(0)=ω3(0)=0ε1(0)=ε2(0)=ε3(0)=0,ε4(0)=1,
and two possible sets of final conditions
(2.4b)ω1(tf)=ω2(tf)=ω3(tf)=0ε1(tf)=ε1f,ε2(tf)=ε2f,ε3(tf)=ε3f,ε4(tf)=ε4f,
for which the final orientation at the final time is completely specified, or
(2.4c)ω1(tf)=ω2(tf)=ω3(tf)=0f1(ε1,ε2,ε3,ε4)=0⋮fn(ε1,ε2,ε3,ε4)=0,for which only some aspect of the final orientation is specified (this is discussed further in Section 2.1).For theunconstrained optimal control problem formulated using (2.1)–(2.3), the Hamiltonian is(2.5)H=λω1ω̇1+λω2ω̇2+λω3ω̇3+λε1ε̇1+λε2ε̇2+λε3ε̇3+λε4ε̇4.
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 αA from each of several high-intensity astronomical sources. Denoting the directions to these sources as σ̂A, where the subscript A can be S (Sun), E (Earth), or M (Moon), the so-called keep-out constraints are then written as follows:(2.6)αA≤cos⁡-1(σ̂⋅σ̂A).
Without loss of generality, σ̂ is assumed to lie along the body-fixed x-axis and its orientation with respect to the inertial frame is then determined from (A.6), with A as the inertial frame and B as the body-fixed frame(2.7)σ1=1-2(ε12+ε32)σ2=2(ε1ε2+ε3ε4)σ3=2(ε3ε1-ε2ε4).
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:(2.8)C=cos⁡-1⁡(σ̂⋅σ̂A)-αA,
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 H̃ is formed by adjoining C̈ to H with Lagrange multiplier μ,(2.9)H̃=H+μC̈,
with the conditions(2.10)μ≥0ifC̈=0μ=0ifC̈>0.
Further, the so-called tangency conditions [11] must also be applied (i.e., for any part of the trajectory on the boundary, not only C, 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.11)σ1(tf)=σ1,f=1-2[ε12(tf)+ε32(tf)]σ2(tf)=σ2,f=2[ε1(tf)ε2(tf)+ε3(tf)ε4(tf)]σ3(tf)=σ3,f=2[ε3(tf)ε1(tf)-ε2(tf)ε4(tf)].


### 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 σ̂A (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 Mmax, the maximum torque along λ̂ has magnitude(2.12)|M⇀λ|=|Mmax⁡λ̂λmax⁡|,
where λmax=max(λ1,λ2,λ3). Equation (2.12) thus maximizes M⇀λ while obeying the limits on the individual control torques M1, M2, and M3. 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(2.13)tf,λ=2Iλθ|Mλ|,
where Iλ 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.Figure 2
Constrained rotation about the keep-out cone axis.


## 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.11)σ1(tf)=σ1,f=1-2[ε12(tf)+ε32(tf)]σ2(tf)=σ2,f=2[ε1(tf)ε2(tf)+ε3(tf)ε4(tf)]σ3(tf)=σ3,f=2[ε3(tf)ε1(tf)-ε2(tf)ε4(tf)].


## 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 σ̂A (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 Mmax, the maximum torque along λ̂ has magnitude(2.12)|M⇀λ|=|Mmax⁡λ̂λmax⁡|,
where λmax=max(λ1,λ2,λ3). Equation (2.12) thus maximizes M⇀λ while obeying the limits on the individual control torques M1, M2, and M3. 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(2.13)tf,λ=2Iλθ|Mλ|,
where Iλ 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.Figure 2
Constrained rotation about the keep-out cone axis.


## 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ω̅i, I̅i, and M̅i, respectively; the corresponding nondimensional quantities are defined as(3.1)ξ=I̅max⁡M̅max⁡ωi=ξω̅iIi=I̅iI̅max⁡Mi=M̅iM̅max⁡,
where ξ is the time unit, and I̅max and M̅max 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 (Figure3) 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 tf=1.9480 (nondimensional units). For comparison, if the rotation axis is constrained to lie along the constraint cone’s axis, the final time would be tf,λ=2.1078.Dynamic response and controls for the case BA1.
(a)(b)(c)Dynamic response and controls [9] for the motion shown in Figure 1.
(a)(b)(c)Figure5 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.Hamiltonian and costates for the case BA1.
(a)(b)


### 3.2. Case BA2

For this forced boundary arc, the constraint cone has half-angleαA = 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 σ̂f is calculated via (A.1), with a⇀=σ̂i=[1,0,0], b⇀=σ̂f, θ=-70 deg., and λ̂=[cos(αA),sin(αA),0].For this case, the solution uses 100 nodes and has a final time oftf=1.3020. For comparison, if the rotation axis is constrained to lie along the constraint cone’s axis, the final time would be tf,λ=2.0966. 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.Dynamic response and controls for the case BA2.
(a)(b)(c)Hamiltonian and costates for the case BA2.
(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 (M1 and M2) exhibit bang-bang behavior whereas M3 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 t=0.025 (corresponding to 22 pseudospectral nodes after t=0) 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].Dynamic response and controls for the case BP1.
(a)(b)(c)Hamiltonian and costates for the case BP1.
(a)(b)Figure 10
Distance from the sensor axis to the constraint boundary (case BP1).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 tf=1.2967.Dynamic response and controls for the case BP2.
(a)(b)(c)Hamiltonian and costates for the case BP2.
(a)(b)Figure 13
Distance from the sensor axis to the constraint boundary (case BP2).


## 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 (Figure3) 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 tf=1.9480 (nondimensional units). For comparison, if the rotation axis is constrained to lie along the constraint cone’s axis, the final time would be tf,λ=2.1078.Dynamic response and controls for the case BA1.
(a)(b)(c)Dynamic response and controls [9] for the motion shown in Figure 1.
(a)(b)(c)Figure5 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.Hamiltonian and costates for the case BA1.
(a)(b)


## 3.2. Case BA2

For this forced boundary arc, the constraint cone has half-angleαA = 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 σ̂f is calculated via (A.1), with a⇀=σ̂i=[1,0,0], b⇀=σ̂f, θ=-70 deg., and λ̂=[cos(αA),sin(αA),0].For this case, the solution uses 100 nodes and has a final time oftf=1.3020. For comparison, if the rotation axis is constrained to lie along the constraint cone’s axis, the final time would be tf,λ=2.0966. 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.Dynamic response and controls for the case BA2.
(a)(b)(c)Hamiltonian and costates for the case BA2.
(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 (M1 and M2) exhibit bang-bang behavior whereas M3 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 t=0.025 (corresponding to 22 pseudospectral nodes after t=0) 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].Dynamic response and controls for the case BP1.
(a)(b)(c)Hamiltonian and costates for the case BP1.
(a)(b)Figure 10
Distance from the sensor axis to the constraint boundary (case BP1).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 tf=1.2967.Dynamic response and controls for the case BP2.
(a)(b)(c)Hamiltonian and costates for the case BP2.
(a)(b)Figure 13
Distance from the sensor axis to the constraint boundary (case BP2).


## 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 componentMi as a sum of N Chebyshev polynomials of the form [18](4.1)Mi(t)≈[∑k=0N-1ckTk(t)]-12c0,
where Tk is the Chebyshev polynomial of degree k, and the coefficients ck 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 Figure14(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) Approximate solution for 1D slewing maneuver, generated by 50 iterations of a particle swarm optimizer. (b) Optimal control solution for the same maneuver, generated by DIDO (running in accurate mode). Total required cpu time for solution: 148 sec using only DIDO, with no initial guess; 76 sec for the combined PSO-DIDO solution.
(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.



---

*Source: 769376-2012-03-25.xml*