Performance Analysis of Generating Function Approach for Optimal Reconfiguration of Formation Flying

  • cc icon
  • ABSTRACT

    The use of generating functions for solving optimal rendezvous problems has an advantage in the sense that it does not require one to guess and iterate the initial costate. This paper presents how to apply generating functions to analyze spacecraft optimal reconfiguration between projected circular orbits. The series-based solution obtained by using generating functions demonstrates excellent convergence and approximation to the nonlinear reference solution obtained from a numerical shooting method. These favorable properties are expected to hold for analyzing optimal formation reconfiguration under perturbations and non-circular reference orbits.


  • KEYWORD

    formation reconfiguration , generating function , optimal control

  • 1. INTRODUCTION

    For optimal formation reconfigurations of satellites in space, various analytical and numerical methods have been developed so far. Carter and Pardis analyzed for the optimal thrust which had the upper and lower limit (Carter & Pardis 1996). Palmer presented an analytic solution for the optimal reconfiguration problem based on the Fourier series expansion of thrust (Palmer 2006). Cho & Park derived a simple analytic solution for linear dynamics in the local vertical, local horizontal frame (Cho & Park 2009). Lee and Park derived approximated analytical solutions considering nonlinearities of the terrestrial gravity, eccentricities of a chief satellite, and J2 effects by perturbing the solution to the linear Hill equation (Lee & Park 2011). These analytical methods generally provide solutions based on the linear dynamics, which are limited to describe the actual dynamic environments precisely. While numerical methods can be employed to yield a precise solution based on nonlinear dynamics, it usually requires initial guess and iterative process for obtaining the initial costate vector accurately.

    Recently, Guibout & Scheeres formulated a formation reconfiguration problem near the L2 Lagrangian point of the Hill three-body dynamics as a Two Point Boundary Value Problem (TPBVP) of the Hamiltonian dynamics, and presented a new solution procedure without initial guess and iterative process of determining the initial momenta by using generating functions appearing in the Hamilton-Jabobi Theory (Guibout & Scheeres 2004). Since then, the generating function approach has been considered as an appropriate tool for reconfiguring formation flying satellites which usually needs to consider varying boundary conditions and time spans for a number of satellites (Guibout & Scheeres 2006). Moreover, Park et al. extended the generating function approach to an optimal control problem, and presented a solution by obtaining the initial costate vector for truncated nonlinear dynamics (Park et al. 2006).

    In this paper, Park et al.’s nonlinear optimal control method with generating functions is applied to the optimal reconfiguration problem on the Projected Circular Orbit (PCO) to analyze the performance of the generating function approach comprehensively. The analysis is conducted on the convergence of trajectories and performance index to those of a nonlinear reference solution obtained from a numerical shooting method. The performance index analyses also exhibit some insightful characteristics by varying PCO configuration parameters under nonlinear dynamics.

    2. METHODOLOGY

       2.1 Necessary conditions for optimal reconfiguration problems

    An optimal reconfiguration problem for satellites flying in formation can be formulated as an optimal rendezvous problem, and can be mathematically framed as follows: minimize the performance index

    image

    subject to nonlinear dynamic equations in affine form with initial and final boundary conditions

    image

    where xRn×1 and uRm×1 are the state and control vectors, respectively. The above optimal reconfiguration problem has the following first order necessary conditions by applying the Pontryagin Principle (Kirk 1970):

    image
    image
    image
    image

    λ is the costate vector, and the subscripts on the bottom of Eqs. (4) and (5) indicate partial derivatives. When Eq. (6) is substituted into Eqs. (3)-(5), the optimal reconfiguration problem is expressed as a TPBVP for a standard Hamiltonian System.

    image

    Once the initial costate vector λ0 is obtained, this TPBVP is transformed into an Initial Value Problem (IVP), and the solution for optimal reconfiguration is obtained through a simple forward integration. When the dynamic system is nonlinear, many numerical methods can be used for obtaining the initial costate vector. However, it requires initial guess and iterative procedure for determining the initial costate vector, which may cause computational burdens to design optimal reconfiguration trajectories for various boundary conditions and time spans of formation flying.

       2.2 Solution for optimal reconfiguration problems in terms of generating functions

    The nonlinear relative motion for optimal reconfiguration problems is expressed in the RSW coordinate system as shown in Fig. 1. The frame XYZ represents the Earth-Centered Inertial (ECI) frame with the center of the earth as the origin, where X directs to the vernal equinox, Y is in the 90° direction to the east from the vernal equinox on the equatorial plane, and Z completes the right hand coordinate system. In the relative RSW coordinate system with its origin located at the reference chief satellite, x is in the opposite direction to the center of the earth, z is perpendicular to the reference orbital plane, and y completes the right hand coordinate system. The distance from the center of the earth to the deputy satellite is expressed as R with the unit vectors

    image

    in the RSW coordinate system as follows.

    image

    Rref is the radius of the reference orbit, which is assumed to be circular in this paper. Eq. (8) is differentiated twice with respect to time to obtain the acceleration of the deputy:

    image

    ω is the angular rate of the reference orbit. Assuming that the deputy is only influenced by the central gravity of the

    earth, the relative dynamic motion of the deputy with respect to the chief in the RSW coordinate system is stated as

    image

    where R is the magnitude of R . Eq. (10) is normalized by Rref and the reciprocal of ω, and is expressed in the state space form as follows.

    image
    image
    image

    Here, [x1 x2 x3]T is the position vector in the RSW coordinate system, and [x4 x5 x6]T is the associated velocity vector.

    Now the optimal control method using the generating function is briefly summarized and applied to the optimal reconfiguration examples. The details regarding the generating function and the Hamiltonian system are given in Park et al. (2006) and Greenwood (1977). First, define another trivial Hamiltonian system for constant flows of initial state and costate vector as

    image

    where K represents the constant Hamiltonian of this system, and is set to be zero without loss of generality. Generating functions define canonical transformations between two different Hamiltonian systems, and thus provide the relationship between the Hamiltonian systems, Eqs. (7) and (14), as follows.

    image

    Theoretically, 2n kinds of generating functions can exist by the assignment of independent variables. Here two out of 4 principal kinds of generating functions are mainly used:

    image

    Each generating function satisfies the following relationships:

    image
    image

    When the associated Hamilton-Jacobi equation, Eq. (17), is solved for the F1 generating function, the initial costate vector is obtained by the relation λ0=-∂F1/∂x0 in Eq. (18). However, it must be noted that the F1 possesses a fundamental singularity, since the pair (X, X0), which is set to be independent variables, loses their independence when t=t0. In order to avoid this singularity, the Hamilton-Jacobi Equation is solved for the F2 generating function, which is then transformed into the F1 by the Legendre transformation which algebraically connects generating functions of different types.

    In order to solve the Hamilton-Jacobi equation, the Hamiltonian is first expanded as a Taylor series with respect to a nominal solution as follows:

    image

    Here, is the increment from the nominal solution. If the nominal solution is simply 0, Eq. (19) can be expressed without for simplicity:

    image

    Next, the F2 is expanded as a Taylor series with respect to the independent variables x and λ0 :

    image

    The maximum order of the F2-expansion is set to be the same as that of the Hamiltonian expansion. Since the control vector is obtained from partial derivatives of the generating function by Eqs. (6) and (18), the order of the optimal solution is k-1 when the order of the generating function is k. Obviously, higher k reflects higher degrees of nonlinearity. By substituting the relation of λ0=-∂F2/∂x into Eq. (20), the Hamilton-Jacobi equation can be arranged in the power series for x and λ0 as follows:

    image

    Balancing Eq. (22) for like kinds of variables results in the ordinary differential equation for the coefficients of Taylor expansion. It turns out that the initial condition F2(x, λ0 t = t0) = x(t0)Tλ0(t0) satisfies Eq. (18). The F2 can now be developed by propagating the initial value problem of the ordinary differential equations for the coefficients of Taylor expansion of generating functions until the desired final time. Then, the Legendre transformation converts the F2 into the F1 algebraically:

    image

    After introducing the boundary conditions for the reconfiguration problem, the initial costate vector and the performance index are expressed as follows (Park et al. 2006):

    image
    image

    The steps up to Eq. (24) show that the optimal control problem can be converted into the IVP, not the TPBVP, without any iterative procedure for determining the initial costate. Eqs. (24) and (25) state that it is necessary to obtain the associated generating function only once. That is, the initial costate vector and performance index can be obtained by simply introducing the given boundary conditions into Eqs. (24) and (25). These characteristics can be advantageous in solving for a number of boundary conditions and time spans for multiple satellites in formation.

    3. RESULT AND ANALYSIS

    In order to demonstrate the convergence of trajectories and performance index based on generating functions to a nonlinear reference solution, the linear, second and third order solutions are developed by using our approach. The expanded generating function and their coefficients are obtained by the Mathematica-based HJ package developed by Guibout, which implemented all the steps from Eq. (14) to Eq. (23) (Guibout & Scheeres 2004). The nominal solution for obtaining the generating function is defined as 0 which is a trivial equilibrium solution of the RSW coordinate system. The initial costate vector satisfying the boundary conditions is determined from the F1 generating function by Eq. (24), which allows one to integrate the Hamiltonian system with the given state vector to develop an optimal trajectory. The performance index is obtained by substituting the boundary conditions into Eq. (25).

       3.1 Error Analysis

    As a verification of the convergence of solution based on generating functions to a nonlinear reference solution, the satisfaction of the prescribed boundary conditions is tested through the error analysis. The linear, second and third order solutions based on generating functions are compared with the reference solution obtained from a numerical shooting method. For error analyses, the position and velocity vectors are considered separately, and are compared with prescribed terminal boundary conditions. In order to describe various formation reconfiguration examples, three cases are addressed as follows:

    Case 1. Errors according to the change of ρf

    Case 2. Errors according to the change of tf

    Case 3. Errors according to the change of α

    Fig. 2 represents the configuration parameters defining the boundary conditions on the PCO. ρ0 and ρf are the radius

    of the PCOs before and after reconfiguration, respectively. α is the phase angle of the initial PCO and △α is the increment in phase, and tf is the transfer time. Note that both △α and (normalized) tf are combined to represent total increment in phase. In all examples, the orbital period P and the radius of circular reference orbit Rref are given as

    image

    where μ is the gravitational constant of the earth.

    The position and velocity errors for Cases 1, 2 and 3 are presented in Figs. 3-8, respectively. The fixed configuration parameters for each case are given as follows:

    image
    image
    image

    In Figs. 3 and 4, while the position and velocity errors of the linear solution steadily increase as ρf increases, those of the second and third order solution are kept consistently low. The results of Case 2 and Case 3 show similar tendencies. Unlike the errors of the linear solution, those of the second and third order solution closely track the error of the nonlinear reference solution. Especially, the errors of the third order solution are quite close to those of the reference solution. As the errors of Cases 1-3 are displayed by changing different PCO configuration parameters, it is confirmed that the higher-order solution based on the generating function provides more accurate solutions than the linear solution for varying boundary conditions considered in the nonlinear optimal reconfiguration problem.

       3.2 Analysis of performance index

    In this section, the performance index based on generating functions is compared with that from a numerical shooting method. The performance index is displayed by varying the initial phase α and the operation time span tf , with other parameters fixed. While the performance indices of the linear to third order solutions are obtained from Eq. (25), that of the nonlinear reference solution is obtained from numerical integration. In addition, the performance index for the nonlinear dynamics shows some distinct characteristics from that for the linear Hill dynamic system. In the PCO reconfiguration problem based on the normalized linear Hill dynamic system, the phase angle α* which gives the minimum performance index can be obtained as follows (Yoo et al. 2007).

    image

    The period and radius of the reference orbit for all examples are defined in Eq. (26).

    Fig. 9 shows the performance index variation according to the initial phase α with tf = P. Other fixed parameters are given as follows:

    image

    As shown in Fig. 9, the performance index variation by the initial phase α from the second and third order solutions are close to that from the numerical shooting method. Especially,rmance index of the third order solution approximates the numerical performance index quite closely. Fig. 9 allows us to analyze α* in the nonlinear dynamic system. Whereas two α*’s exist symmetrically in case of the linear solution, which can be confirmed by Eq. (30), only one α* has the minimum performance index in case of the second and third order solutions as well as the nonlinear reference solution.

    Figs. 10-12 show the performance index variation by varying the initial phase α for a different transfer time tf, where the phase angle α is expressed in radians and the distance from the origin is the magnitude of the performance index of the associated phase. Fixed parameters in Figs. 10-12 are as shown in Eq. (31). The polar graph of the performance index shows that the optimal initial phase α* changes gradually as the transfer time tf changes. This change appears similarly in both the linear and nonlinear solutions. However, whereas two α*’s exist for the linear solution in the symmetric form, only one α* appears in the third order solution as well as the nonlinear reference solution. This result confirms the analysis in Fig. 9 that only one α* exists if nonlinearity prevails.

    When a number of satellites maneuver on the PCO, it is often considered that satellites should be distributed uniformly

    before and after reconfiguration for space interferometry. Figs. 13 and 14 show the sum of performance indices for three and six satellites that perform reconfigurations in equidistant intervals in terms of the initial phase α of the reference deputy satellite.

    Other constant parameters are set as follows.

    image

    Figs. 13 and 14 show that the sum of performance indices of the third order solution properly approximates that of the nonlinear reference solution obtained from the numerical shooting method. Fig. 13 shows the case of three satellites performing the reconfiguration, in which only the sum of performance indices of the linear solution has a fixed value regardless of α of the reference deputy satellite; those of the nonlinear solution show periodic variations. Even though the boundary conditions in Eq. (32) represent big differences in the initial and terminal PCO radius, and thus possibly incur high nonlinearities of the central gravity field of the earth, Fig. 14 shows that for six satellites, the sum of performance indices in both the linear and nonlinear solutions always shows a fixed value regardless of α of the reference satellite. This phenomena turns out to be consistent, as long as the number of satellite is more than 6.

    4. CONCLUSION

    The optimal control problem for the formation reconfiguration on the PCO has been solved using the generating function. The higher-order solution obtained by using the generating

    function converges to the nonlinear reference solution without any iterative procedure of obtaining the initial costate for the various boundary conditions in the nonlinear dynamic system. The semi-analytical performance index derived from the generating function closely approximates the performance index obtained from the direct numerical integration. It allows us to analyze performance index characteristics of the nonlinear formation reconfiguration, which shows distinct patterns from the linear reconfiguration. The analysis for the errors and performance indices demonstrate that the generating function approach is an appropriate method for nonlinear optimal reconfiguration of formation flying satellites. In the future, the generating function approach will be applied to more general case of perturbed elliptical reference orbit.

  • 1. Carter TE, Pardis CJ (1996) Optimal Power-Limited Rendezvous with Upper and Lower Bounds on Thrust [JGCD] Vol.19 P.1124-1133 google doi
  • 2. Cho HC, Park SY (2009) Analytic Solution for Fuel-Optimal Reconfiguration in Relative Motion [JOTA] Vol.141 P.495-521 google doi
  • 3. Greenwood DT 1977 Classical Dynamics P.187-271 google
  • 4. Guibout VM, Scheeres DJ (2004) Solving Relative Two-Point Boundary Value Problems: Spacecraft Formation Flight Transfers Application [JGCD] Vol.27 P.693-704 google doi
  • 5. Guibout VM, Scheeres DJ (2006) Spacecraft Formation Dynamics and Design [JGCD] Vol.29 P.121-133 google doi
  • 6. Kirk DE 1970 Optimal Control Theory P.184-325 google
  • 7. Lee S, Park SY (2011) Approximate Analytical Solutions to Optimal Reconfiguration Problems in Perturbed Satellite Relative Motion [JGCD] Vol.34 P.1097-1111 google doi
  • 8. Palmer P (2006) Optimal Relocation of Satellites Flying in Near-Circular-Orbit Formations [JGCD] Vol.29 P.519-526 google doi
  • 9. Park C, Guibout VM, Scheeres DJ (2006) Solving Optimal Continuous Thrust Rendezvous Problems with Generating Function [JGCD] Vol.29 P.321-331 google doi
  • 10. Yoo SM, Park SY, Choi KH 2007 A Fuel Balancing Method for Reconfiguration of Satellite Formation Flying [International Conference on Control, Automation and Systems 2007] google
  • [Fig. 1.] The Earth-Centered Inertial (ECI) frames and the RSW frames.
    The Earth-Centered Inertial (ECI) frames and the RSW frames.
  • [Fig. 2.] Configuration parameters of projected circular orbit.
    Configuration parameters of projected circular orbit.
  • [Fig. 3.] Case 1: Position error.
    Case 1: Position error.
  • [Fig. 4.] Case 1: Velocity error.
    Case 1: Velocity error.
  • [Fig. 5.] Case 2: Position error.
    Case 2: Position error.
  • [Fig. 6.] Case 2: Velocity error.
    Case 2: Velocity error.
  • [Fig. 7.] Case 3: Position error.
    Case 3: Position error.
  • [Fig. 8.] Case 3: Velocity error.
    Case 3: Velocity error.
  • [Fig. 9.] Cost Variation vs. Initial phase angle α (t=P).
    Cost Variation vs. Initial phase angle α (t=P).
  • [Fig. 10.] Cost Variation vs. Initial phase angle α and terminal time tf (Linear solution).
    Cost Variation vs. Initial phase angle α and terminal time tf (Linear solution).
  • [Fig. 11.] Cost Variation vs. Initial phase angle α and terminal time tf (third order solution).
    Cost Variation vs. Initial phase angle α and terminal time tf (third order solution).
  • [Fig. 12.] Cost Variation vs. Initial phase angle α and terminal time tf (Nonlinear reference solution).
    Cost Variation vs. Initial phase angle α and terminal time tf (Nonlinear reference solution).
  • [Fig. 13.] Total cost Variation vs. Initial phase angle α of reference spacecraft (3 Spacecrafts).
    Total cost Variation vs. Initial phase angle α of reference spacecraft (3 Spacecrafts).
  • [Fig. 14.] Total cost Variation vs. Initial phase angle α of reference spacecraft (6 Spacecrafts).
    Total cost Variation vs. Initial phase angle α of reference spacecraft (6 Spacecrafts).