Maximum-likelihood estimation optimizer for constrained, time-optimal satellite reorientation

Acta Astronautica 103 (2014) 185–192

Contents lists available at ScienceDirect

Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro

Maximum-likelihood estimation optimizer for constrained, time-optimal satellite reorientation Robert G. Melton n 229 Hammond Bldg., The Pennsylvania State University, University Park, Pennsylvania 16802, USA

a r t i c l e i n f o

abstract

Article history: Received 23 December 2013 Received in revised form 22 May 2014 Accepted 21 June 2014 Available online 7 July 2014

The Covariance Matrix Adaptation-Evolutionary Strategy (CMA-ES) method provides a high-quality estimate of the control solution for an unconstrained satellite reorientation problem, and rapid, useful guesses needed for high-fidelity methods that can solve timeoptimal reorientation problems with multiple path constraints. The CMA-ES algorithm offers two significant advantages over heuristic methods such as Particle Swarm or Bacteria Foraging Optimisation: it builds an approximation to the covariance matrix for the cost function, and uses that to determine a direction of maximum likelihood for the search, reducing the chance of stagnation; and it achieves second-order, quasi-Newton convergence behaviour. & 2014 IAA. Published by Elsevier Ltd. All rights reserved.

Keywords: Time-optimal reorientation Covariance matrix adaptation Optimal attitude control

1. Introduction This paper is one in a series that examines methods for determining approximate solutions to the constrained, time-optimal satellite reorientation problem. It evaluates the performance of a statistics-based method for use either as the first stage in a hybrid method or as a standalone approximate solver for the unconstrained problem. 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 sunline to avoid excessive solar heating. For scientific missions, observational instruments frequently must be kept

n

Tel.: þ 1 814 865 1185; fax: þ 1 814 865 7092. E-mail address: [email protected]

http://dx.doi.org/10.1016/j.actaastro.2014.06.032 0094-5765/& 2014 IAA. Published by Elsevier Ltd. All rights reserved.

beyond a specified minimum angular distance from highintensity 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] considered time-optimal slews for a rigid spacecraft whose mass distribution is spherically symmetric, and which has equal control-torque authority for all three axes. Despite the symmetry of the system, the intuitively obvious time-optimal solution is not a λ-rotation (rotation through a specified angle about an axis fixed in both the body and an external reference frame) about the eigenaxis. 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, forgetting the constraints imposed by Euler's equations of rigidbody motion. Bilimoria and Wie found that the timeoptimal solution includes precessional motion to achieve a lower time (approximately 10% less) than that obtained with an eigenaxis maneuver. Further, they found 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

186

R.G. Melton / Acta Astronautica 103 (2014) 185–192

in terms of a fictitious λ-rotation, with associated angle θ): for values of θ less than 72 degrees, the control history was found to contain seven switches between directions of the control torque components; larger values of θ require only five switches. Several subsequent papers 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 optimisation [5]. Bai and Junkins [6] reconsidered the original problem (spherically symmetric mass distribution, three equal control torque limits) and found that at least two locally optimum solutions exist for reorientations of less than 72 degrees (one of which requires only six switches) if the controls are independently limited. Further, they proved that if the total control vector is constrained to have a maximum magnitude (i.e., with the orthogonal control components not necessarily independent), then the time-optimal solution is the eigenaxis maneuver. Early work on constrained maneuvers focused upon generating feasible solutions [7,8] but did not attempt to find optimal solutions. The advent of pseudospectral methods and their application to various trajectory optimisation problems has made it possible to study problems with challenging constraints [9]. These methods are not always fully automatic and may require an intelligent initial guess. Hybrid methods, using heuristic algorithms such as Particle Swarm Optimisation (PSO) [10] and Bacteria Foraging Optimisation (BFO) [11] to generate initial guesses for the states and controls, have been shown to significantly reduce overall computation time [12]. The Swift spacecraft1 offers a useful model for the problem: the satellite 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 the these instruments, the slewing motion is constrained to prevent the telescopes' common axis (referred to as the sensor axis) from entering established “keep-out” zones, defined as cones with central axes pointing to the Sun, Earth and Moon, with specified half-angles of 47 degrees, 33 degrees and 23 degrees respectively. An example is shown in Fig. 1, which was generated using the DIDO software [13]. Initially the sensor axis is oriented out of the plane of the page toward the reader; its final orientation is on the opposite side of the Sun cone. The time-optimal solution carries the sensor axis through a narrow gap of 0.1 (not clearly visible) between the Sun and Moon cones. It should be noted that the obvious alternative path, around the lower side of the Sun cone, would require more time because of the specified final orientations of all three body axes.

1

http://heasarc.nasa.gov/docs/swift/swiftsc.html

2. Problem statement The problem is formulated as a Mayer optimal control problem, with performance index J ¼ tf

ð1Þ

where tf is the final time to be minimised. Euler's equations of rigid-body motion describe the system dynamics

ω_ 1 ¼ ½M1  ω2 ω3 ðI3  I2 Þ=I1 ω_ 2 ¼ ½M2  ω3 ω1 ðI1  I3 Þ=I2 ω_ 3 ¼ ½M3  ω1 ω2 ðI2  I1 Þ=I3

ð2Þ

These 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 3 2 3 2 ϵ_ 1 ϵ4  ϵ3 ϵ2 ϵ1 2 ω1 3 6_ 7 6 ϵ4  ϵ1 ϵ2 7 76 6 ϵ 2 7 6 ϵ3 ω2 7 7 76 6 7¼6 ð3Þ 7 76 6 ϵ_ 3 7 6  ϵ2 ϵ ϵ ϵ 4 1 4 3 5 ω3 5 4 5 4 ϵ_ 4  ϵ1  ϵ2  ϵ3 ϵ4 0 The problem to be considered is a rest-to-rest maneuver, consisting of a rotation about the z-axis through an angle of 3π =4 radians with boundary conditions

ω1 ð0Þ ¼ ω2 ð0Þ ¼ ω3 ð0Þ ¼ 0 ϵ1 ð0Þ ¼ ϵ2 ð0Þ ¼ ϵ3 ð0Þ ¼ 0; ϵ4 ð0Þ ¼ 1 ω1 ðt f Þ ¼ ω2 ðt f Þ ¼ ω3 ðt f Þ ¼ 0     ϵ1 t f ¼ ϵ2 t f ¼ 0     3π ϵ3 t f ¼ sin 8     3π ϵ4 t f ¼ cos 8

ð4Þ

ð5Þ

In the numerical calculations, time is scaled by the factor pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi I=M max and the control torques are scaled by the factor Mmax. For the optimal control problem formulated using Eqs. (1)–(3), the Hamiltonian is _ 1 þ λω1 ω _ 2 þ λω3 ω _3 H ¼ λω1 ω þ λϵ1 ϵ_ 1 þ λϵ2 ϵ_ 2 þ λϵ3 ϵ_ 3 þ λϵ4 ϵ_ 4

ð6Þ

For spacecraft of the type being modelled here (Swift, or other astronomical missions), one or more sensors is 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 αx from each of several high-intensity light sources. Denoting the directions to these sources as σ^ x , where the subscript x can be S (Sun), E (Earth), or M (Moon), the so-called keep-out constraints are then written as C x ¼ σ^  σ^ x  cos ðαx Þ r 0

ð7Þ

Without loss of generality, σ^ is assumed to lie along the bodyfixed x-axis and its orientation with respect to the inertial frame is then determined [14]

σ 1 ¼ 1  2ðϵ21 þ ϵ23 Þ σ 2 ¼ 2ðϵ1 ϵ2 þ ϵ3 ϵ4 Þ

R.G. Melton / Acta Astronautica 103 (2014) 185–192

187

Fig. 1. Constrained time-optimal path of the sensor axis.

σ 3 ¼ 2ðϵ3 ϵ1  ϵ2 ϵ4 Þ

ð8Þ

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. This assumption will be justified later in the numerical results. Inclusion of the path constraints, Eq. (7), requires augmenting the Hamiltonian with terms of the form μx C x , with the multiplier μx Z 0 if the constraint is active and μx ¼ 0 otherwise. While this problem could, in principle, be solved using an indirect method, the path constraints, Eq. (7), create especially difficult challenges. A direct approach, such as a pseudospectral method, offers the potential advantage of easier problem formulation; however, for this problem, an optimal solution, or even a feasible solution, is not automatically guaranteed. For example, for the sensor path shown in Fig. 1, the DIDO software required approximately 12 h to generate the high-fidelity approximate solution; the states and control torques are shown in Fig. 2, and the Hamiltonian and costates in Fig. 3. The near-constant value of the Hamiltonian gives confidence that the optimal solution has been determined. To understand the challenge posed by the path constraints, one must recognise that the direct methods all employ an implicit integration of the governing system equations. Such solutions require initial estimates (sometimes these are just outright guesses) of the states and controls at the discrete times (nodes) for which the overall solution is computed. Lacking any user-provided initial estimate of the states and controls at these nodes, the logical automated estimate for the states is simply a linear interpolation between the states at the initial and final times; the control is also assumed to be a linear function between the minimum and maximum specified values, mapped over the specified range of times. Such a linear relationship violates various dynamic and kinematic

constraints, viz. Eqs. (2) and (3). Consequently, the direct-method solver must then expend some computational effort to find a feasible solution. As determined in previous work [9,12], this is not always successful or can take considerable time (on the order of 1–12 h). A good initial guess for the states and controls can be obtained via explicit numerical integration of the equations of motion, using a nonlinear optimizer to determine the maneuver time and the control-torque values at the nodes. Such hybrid approaches result in lower overall computation time, but further improvement is sought, initially for research purposes but ultimately with the intent of onboard implementation. 3. Summary of the covariance matrix method Hansen et al. [15,16] present the Covariance Matrix Adaptation-Evolution Strategy (CMA-ES) in which the covariance matrix C of the objective (fitness) function J, in effect, replaces the inverse Hessian matrix of J. The method generates successive approximations to C via a relatively small number of evaluations of J; new search directions are generated in the directions of the principal axes of C, maximising the likelihood of finding improved solution points. In these regards, the method resembles quasi-Newton, variable-metric optimizers [17,18]. The complete algorithm for the basic method appears in the Appendix; a brief summary is presented here. The objective function is modelled as a quadratic surface 1 J ¼ xT Hx 2

ð9Þ

where x is the vector of unknowns. Taking search steps of the form xi ¼ x þ αC  1=2 δx

ð10Þ

where x is the mean value of x, α is the step size, and δxi is normally distributed, is consistent with the initial

188

R.G. Melton / Acta Astronautica 103 (2014) 185–192

Fig. 2. State and torque histories (from DIDO) for the sensor-axis path shown in Fig. 1. Hamiltonian −0.96 −0.97

H

−0.98 −0.99 −1 −1.01 −1.02

0

0.5

1

1.5

2

2.5

3

3.5

time Costates 4

λω

costates

2

λω

0

λω

−2

λε

−4

λε λε

−6 −8

λε

0

0.5

1

1.5

2

2.5

3

3.5

time

Fig. 3. Hamiltonian and costate histories (from DIDO) for the sensor-axis path shown in Fig. 1.

assignment of H  1 ¼ C ¼ I (the identity matrix) in Eq. (9). In each generation, i ¼ 1; …; λ samples of the search space are taken and sorted by their corresponding J values. The best half of these form a subset that is used to update C. Also the new step size α for the next generation is determined using information about the changes in x during the present generation. The update process employs an eigendecomposition of C, generating new sample points along the principal axes of C, with distances from the present mean of x determined by the corresponding eigenvalues of C. Hansen2 notes that the method

2

https://www.lri.fr/  hansen/cmatutorial.pdf

is less prone to stagnation, as compared with heuristic methods such as PSO and BFO; in the time-optimal reorientation cases examined thus far, stagnation has not been observed, and the convergence rate is very good when all terms in the cost function are quadratic.

4. Constraints While linear penalty functions serve well for handling constraints in conjunction with the heuristic methods [12], the CMA-ES method performs much better with penalties that are quadratic in the amount of the constraint violation.

R.G. Melton / Acta Astronautica 103 (2014) 185–192

For equality constraints of the form g i ðxÞ ¼ 0

ð11Þ

the penalty function is Gi ðxÞ ¼ keq;i ½maxf0; j g i ðxÞj  Δeq;i g2

hj ðxÞ r 0

ð13Þ

the penalty function is H j ðxÞ ¼ kineq;j ½maxf0; hj ðxÞ  Δineq;j g2

ð14Þ

where the keq;i and kineq;j are weights whose values are selected to make the penalty terms have the same order of magnitude as the principal part of the cost function; the user-specified tolerances Δeq;i and Δineq;j result in zero penalties if j g i j o Δeq;i and hj o Δineq;j . An additional consideration for CMA-ES is that simply bounding the control-torque values results in poor performance of the optimizer. For example, if the value of a particular torque exceeds the allowable bound at one of the discrete nodes, simply setting the value to the numerical bound disrupts the covariance-update process. Therefore, simple limits for the torques must be replaced by quadratic penalty terms of the form in Eq. (14). Further, because CMA-ES is a quasi-Newton method and therefore relies on the assumption of a well-defined Hessian, the principal part of the cost function cannot be linear, as in 2 Eq. (1); replacing it with tf is sufficient to give good convergence. The resulting modified cost function is then N eq

N ineq

i¼1

j¼1

J 0 ¼ t 2f þ ∑ Gi þ ∑ H j

ð15Þ

where Neq and Nineq are the number of equality and inequality constraints, respectively.

5. Modelling the control torques In this work, two models for the control torques have been considered: linear combinations of Chebyshev polynomials and cubic splines. In the case of Chebyshev polynomials, each torque is represented as NT  1

M i t ¼ ∑ cj T j t  j¼0

c0 2

ð16Þ

where the Tj are Chebyshev polynomials, and the cj are coefficients determined from the torque values calculated using any of the three heuristic algorithms (i.e., the decision variables in this formulation consist of the final time tf, and the NT torque values for each control). For each particle (or solution vector), the time variable t in the domain [0, tf] is converted to τ in the domain [  1,1], which is the domain of the Chebyshev functions, via the linear relation

τ¼

2t  ðt f þt 0 Þ tf  t0

and the Chebyshev coefficients are calculated using 2 3 1 1 1 k pjk  pjk  N 2 NT 6 27 2 ¼ 2 ∑T f cos 2 cj ¼ ∑ f 4 cos 5 cos NT NT NT NT k ¼ 1 NT k ¼ 1 k

ð12Þ

For inequality constraints of the form

ð17Þ

189

ð18Þ where the function f ðÞ represents a continuous torque function and the f k are torque values calculated by the optimizer at discrete times tk. Chebyshev polynomials are used here because of the property that all the Tj have extrema at 71 on the domain  1 r τ r1, making the representation in Eq. (16) readily compatible with the torque limits. Eqs. (16)–(18) are employed in the explicit integration to compute the control torques. Each Tj value is calculated using the Clenshaw recursion algorithm [18]. Another choice for representing the control torques is cubic splines. Splines have been used successfully in conjunction with PSOs for several types of trajectory optimisation and robot maneuver planning problems [19], where the authors used B-splines and included the spline-degree as one of the unknown parameters to be determined during the optimisation. Of particular interest is their solution for the robotic-arm problem that has a bang-bang control structure. 6. Results For comparison purposes, first consider the unconstrained problem of moving the sensor axis to a position 135 degree from its original direction; the final orientation of the body z-axis must be the same as its original orientation. This maneuver corresponds to the initial and final sensor-axis orientations shown in Fig. 1, but with the constraint cones removed. Using DIDO with 40 nodes (any fewer nodes gives insufficient resolution to capture the switching structure), the resulting states and control torques are as shown in Fig. 4. This calculation required approximately 70 min.3 Here it is important to verify that the assumption of a fast reorientation is reasonable. The results in Fig. 4 show a reorientation time of approximately 2.9 dimensionless time units. For a satellite with mass moment of inertia I ¼8000 kg m2 for each axis, and maximum control torque axis, the required physical time M max ¼ 100 N m forpeach ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi to reorient is 2:9  I=Mmax  27 s. This demonstrates that the reorientation can occur relatively fast and that the directions to the high-intensity sources (Sun, Earth, Moon) would remain essentially constant even for a satellite in a low orbit. For the same unconstrained case, using Chebyshev polynomials to model the control torques (number of torque values per control axis must equal the number of Chebyshev polynomials), and using the CMA-ES method, the resulting approximate optimal solution is that shown in Fig. 5. Note that the m3 torque (about the z-axis) is roughly the same as that in Fig. 4, although the two orthogonal torques bear little resemblance to their 3 All computation times refer to CPU times for Matlab code executed on an Intel(R) Core(TM) i3 2.27 GHz processor.

190

R.G. Melton / Acta Astronautica 103 (2014) 185–192

ωi

Angular Velocity 2

ω

0

ω2

−2

1

ω 0

0.5

1

1.5

2

2.5

3

3

time Euler Parameters

εi

2

ε1 ε

2

0 −2

ε3 0

0.5

1

1.5

2

2.5

ε

3

4

time

Mi

Control Torques 1

M

0

M2

−1

1

0

0.5

1

1.5

2

2.5

M3

3

time

Fig. 4. States and control torques for the unconstrained case, generated by DIDO. 1.5

m1 m2 m3

1

0.5

0

−0.5

−1

−1.5

0

0.5

1

1.5

2

2.5

3

3.5

time

Fig. 5. Torques for the unconstrained case, generated by CMA-ES, using Chebyshev polynomials.

counterparts in Fig. 4. Nevertheless, using this approximate solution (torques and states) as the initial guess in DIDO for the constrained case, yields results identical to those in (Figs. 1 and 2), but with significantly lower overall computation time. Note that some torque values exceed magnitudes of 1, a consequence of using a penalty method to constrain the torques. While the quality of the CMA-ES solution here would not be acceptable on its own, it does provide a physically feasible guess that satisfies the dynamic constraints (equations of motion), and nearly satisfies the kinematic constraints (final orientations of the body axes); this reduces the computational burden for the pseudospectral method. Table 1 compares the overall computation times for the constrained problem, using the pseudospectral method (DIDO) alone and the hybrid combinations where the first stage is either PSO, BFO or CMA-ES. Each hybrid case uses 40 nodes for the first stage and 300 nodes for the second (DIDO) stage. It should be noted that the result for DIDO alone use 40 nodes, but the quality of the solution is poor, with the resulting Hamiltonian varying between  1.2 and  0.8. For each hybrid method, the first-stage algorithm

Table 1 Computation times for the constrained maneuver. Method

Computation time (sec)

Pseudospectral alone PSO-pseudospectral BFO-pseudospectral CMAES-pseudospectral

4189 196 þ 1344 ¼1540 167 þ 1034 ¼1101 81 þ 695 ¼776

was executed 10 times with different seed values for the random-number generator; computation times shown in Table 1 are the smallest values observed (the largest execution times were no more than 20% greater). Finally, the CMA-ES method was applied to the unconstrained case, using cubic splines to represent the control torques (40 nodes). The computation required approximately 58 min, resulting in the control torques shown in Fig. 6 and the corresponding angular velocities shown in Fig. 7. Comparison with Fig. 4 shows close agreement for the control torques and the angular velocities about all three axes. It should be noted that similar attempts to find

R.G. Melton / Acta Astronautica 103 (2014) 185–192 m1

2

stages of the search) may prove more effective in reaching an approximate solution more rapidly, leading to its possible use as a stand-alone solver for the constrained problem.

0 −2

0

0.5

1

1.5

2

2.5

3

time m 2

Appendix A. CMA-ES algorithm

2

0 −2

0

0.5

1

1.5

2

2.5

3

time m 2

3

0 −2

0

0.5

1

1.5

2

2.5

3

time

Fig. 6. Torques for the unconstrained case, generated by CMA-ES, using cubic splines.

ω

1

1.4

Set parameters n ¼dimension of the unknown vector x gmax ¼ maximum number of generations   λ λ ¼ 4 þ ⌊3 ln nc; μ ¼ 2 ! λ 1 w0  lni; i ¼ 1; …; μ wi ¼ ∑μ i w0 ; w0i ¼ ln þ 2 2 j ¼ 1 j ! rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi μeff þ 2 μeff  1  1 þ cα ; dα ¼ 1 þ 2 max 0; cα ¼ n þ μeff þ 5 nþ1 cc ¼

4þ μeff =n n þ 4 þ 2μeff =n 2 ðnþ 1:3Þ2 þ μeff

cμ ¼ min 1  c1 ; αμ

ω2 ω

1.2

3

!

μeff  2þ 1=μeff ; ðn þ 2Þ2 þ αμ μeff =2

αμ ¼ 2

Initialization Set evolution paths pα ¼ 0, pc ¼ 0, covariance matrix C ¼ I, and g ¼ 0 Choose distribution mean x A Rn and step size α A R

1 0.8 ωi

This algorithm is taken from Hansen's tutorial4 on CMA-ES, but with some modifications to the notation. The symbol Nðβ; ΓÞ indicates a normally distributed random variable with mean β and covariance Γ.

c1 ¼

1.6

while (g r g max ) and maxD=minD r 1014 ) do Sample new population of search points, for k ¼ 1; …; λ zk  Nð0; IÞ yk ¼ BDzk  Nð0; CÞ [The n columns of B consist of the eigenvectors of C, and D is diagonal, formed from the eigenvalues of C.] x ¼ x þ αyk  Nðx ; CÞ Selection and recombination5 μ μ 〈y〉w ¼ ∑i ¼ 1 wi yi:λ where ∑i ¼ 1 wi ¼ 1; wi 40

0.6 0.4 0.2 0 −0.2 −0.4

191

μ

0

0.5

1

1.5

2

2.5

3

time

Fig. 7. Angular velocities, generated by CMA-ES, with the splinetorque model.

the unconstrained solution using the cubic spline model and either the PSO or BFO method resulted in very poor performance, with seemingly random oscillations in the control torques, and repeated stagnation in the search progress.

x ¼ x þ α〈y〉w ¼ ∑i ¼ 1 wi xi:λ Step-size control pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 pα ¼ ð1  cα Þpα þ hα cc ð2  cc Þμeff C  2 〈y〉w

pα J α ¼ α  exp dcαα E JJNð0;IÞ J 1

Covariance matrix adaptation 8   J pα J 2 > > EJ N ð0; IÞ J o 1:4þ < 1; if qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n þ 1 hα ¼ 1 ð1  cα Þ2ðg þ 1Þ > > : 0; otherwise

δðhα Þ ¼ ð1 hα Þcc ð2 cc Þ pc ¼ ð1  cc Þpc þ hα

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cc ð2  cc Þμeff 〈y〉w

μ

C ¼ ð1  c1  cμ ÞðCÞ þ c1 ðpc pTc þ δðhα ÞC þ cμ ∑ wi yi:λ yTi:λ i¼1

end while

7. Conclusion The Covariance Matrix Adaptation-Evolutionary Strategy (CMA-ES) shows promise for use as the first stage in a hybrid method involving a second-stage pseudospectral method, when applied to the constrained, time-optimal reorientation maneuver problem. Significant overall computation time improvement was found when compared to other first-stage methods, such as Particle Swarm Optimisation (PSO) and Bacteria Foraging Optimisation (BFO). CMA-ES also performs well as a stand-alone optimizer for the unconstrained problem. Modified versions of CMA-ES (in which the eigendecomposition is removed in the early

References [1] K.D. Bilimoria, B. Wie, Time-optimal three-axis reorientation of a rigid spacecraft, J. Guid. Control Dyn. 16 (3) (1993) 446–452.

4 An excellent development of the method, complete with justifications for the particular choices of the various parameters, as well as Matlab code for the basic method, is available at https://www.lri.fr/  hansen/cmatutorial.pdf. 5 The notation i: λ indicates that i is selected from the best μ points (i.e., lowest values of J 0 ) in the population of size λ.

192

R.G. Melton / Acta Astronautica 103 (2014) 185–192

[2] H. Shen, P. Tsiotras, Time-optimal control of an axi-symmetric rigid spacecraft sing two controls, J. Guid. Control Dyn. 22 (5) (1999) 682–694. [3] R.J. Proulx, I.M. Ross, Time-optimal reorientation of asymmetric rigid bodies, Adv. Astronaut. Sci. 109 (part II) (2002) 1207–1227. [4] S.L. Scrivener, R.C. Thompson, Time-optimal reorientation of a rigid spacecraft using collocation and nonlinear programming, Adv. Astronaut. Sci. 85 (part III) (1993) 1905–1924. [5] S.-W. Liu, T. Singh, Fuel/time optimal control of spacecraft maneuvers, J. Guid. Control Dyn. 20 (2) (1997) 394–397. [6] X. Bai, J. Junkins, New results for time-optimal three-axis reorientation of a rigid spacecraft, J. Guid. Control Dyn. 32 (4) (2009) 1071–1076. [7] H. Hablani, Attitude commands avoiding bright objects and maintaining communication with ground station, J. Guid. Control Dyn. 22 (6) (1999) 759–767. [8] G. Mengali, A. Quarta, Spacecraft control with constrained fast reorientation and accurate pointing, Aeronaut. J. 108 (1080) (2004) 85–91. [9] R.G. Melton, Constrained time-optimal slewing maneuvers for rigid spacecraft, Adv. Astronaut. Sci. 135 (2010) 107–126. [10] J. Kennedy, R. Eberhart, Particle swarm opimization, in: Proceedings of IEEE International Conference on Neural Networks, IV, 1995, pp. 1942–1948. [11] H. Chen, Y. Zhu, K. Hu, Adaptive bacterial foraging optimization, Abstract and Applied Analysis, http://dx.doi.org/10.1155/2011/ 108269.

[12] R.G. Melton, Hybrid methods for determining time-optimal, constrained spacecraft reorientation maneuvers, Acta Astronaut. 94 (1) (2014) 294–301. http://dx.doi.org/10.1016/j.actaastro.2013.05.007. [13] I. Ross, A beginner's guide to DIDO (ver 7.3): A Matlab application package for solving optimal control problems, Document #TR-711, Elissar, LLC, 2007. [14] T.R. Kane, P.W. Likins, D.A. Levinson, Spacecraft Dynamics, McGrawHill, New York, 1983 (Chapter 1). [15] N.A. Hansen, A. Ostermeier, A. Gawelczyk, On the adaptation of arbitrary normal mutation distributions in evolutionary strategies: the generating set adaptation, in: Eshelman, Ed., Proceedings of the Sixth International Conference on Genetic Algorithms, Pittsburgh, PA, 1995, pp. 57–64. [16] N.A. Hansen, A. Ostermeier, Adapting arbitrary normal mutation distributions in evolution strategies: the covariance matrix adaptation, in: Proceedings of the 1996 IEEE International Conference on Evolutionary Computation, 1996, pp. 312–317. [17] D. Luenberger, Introduction to Linear and Nonlinear Programming, Addison-Wesley, Reading, MA, 1973, pp. 189–204. [18] W. Press, S. Teukolsky, W. Vetterling, B. Flannery, Numerical Recipes: The Art of Scientific Computing, 2nd ed. Cambridge University Press, Cambridge, England, 2005, pp. 196–199, 430–432. [19] P. Ghosh, B. Conway, Numerical trajectory optimization with swarm intelligence and dynamic assignment of solution structure, J. Guid. Control Dyn. 35 (4) (2012) 1178–1191.