Flight Trajectory Optimization Using Genetic Algorithm Combined

Document Sample
Flight Trajectory Optimization Using Genetic Algorithm Combined Powered By Docstoc
					                                                                                   e-journal: Information Technology
                                                                                     for Economics & Management
    ITEM, Volume 1, No 1, 2001, Paper 10                                                    ISSN 1643-8949

    Flight Trajectory Optimization Using Genetic Algorithm Combined
                          with Gradient Method
                 Nobuhiro YOKOYAMA: School of Aeronautics and Astronautics,
               the University of Tokyo 7-3-1 Hongo, Bunkyo-ku, Tokyo 113, JAPAN
                               email: tt06348@mail.ecc.u-tokyo.ac.jp
         Shinji SUZUKI: School of Aeronautics and Astronautics, the University of Tokyo,
                                   tshinji@mail.ecc.u-tokyo.ac.jp


Abstaract                                                        a nonlinear optimization problem by descritizing the
                                                                 time functions such as control or state variables using a
This paper considers the numerical method to solve tra-          set of discrete variables. Sophisticated algorithms with
jectory optimization problems where efficient flight tra-           gradient information for nonlinear optimization prob-
jectories are searched to minimize a performance index           lems make it easier to deal with complicated equal or
with specified constrained. While a wide variety of nu-           unequal constrains. However, the gradient algorithms
merical methods have been developed, the direct appli-           require the delicate selection of initial guess of the opti-
cation of a nonlinear programming method is frequently           mal solution. This paper tries to combine non-gradient
used in the recent applications, which transforms the            techniques with the gradient methods to solve the com-
original problem into a nonlinear optimization problem.          plicated trajectory optimization methods without accu-
Although these methods can be applied to complicated             rate initial solutions.
applications, it is very difficult to find reasonable initial             As a non-gradient optimization algorithm, a Ge-
solutions since practical problems have high nonlinear-          netic Algorithms (GA) is employed, which uses search
lity. This paper tries to apply the Genetic Algorithm            procedures based on the mechanics of natural genet-
(GA) method. Sine the GA approach is one of the ran-             ics, which combines a Darwinian survival-of-the-fittest
dom search methods, the initial solutions are selected           strategy to eliminate unfit characteristics and uses ran-
at random. As shown in numerical examples, the GA                dom information exchange(3) . While the original GA
approach has slow convergence characteristics and the            codes continuous variables into a string of bit code, the
meaningless fluctuation in the solutions. Therefore, the          continuous variables can be dealt with by an advanced
gradient approach is utilized to refine the solutions. As         continuous parameter GA(4) which defines a chromo-
the gradient method, the authors uses the BDH method             some as an array of floating point numbers. Since
which approximates the state and control variables us-           the GA approach is categorized as random search al-
ing linear interpolation and the collocation method is           gorithm, initial solutions are selected at random. How-
used to satisfy the dynamic equations. The obtained              ever, the convergence of the search process is usually
Hessian matrix can be approximated in a Block Diag-              very slow.
onal Form which is suitable for efficient computation.
The proposed approach is applied to the optimal accent                  This paper utilizes the continuos parameter GA
trajectory optimization of a spaceplane. The spaceplane          in order to select the initial solutions for the BDH
is a future space transportation vehicle and its concep-         method(5) which is one of the gradient methods pro-
tual design is investigated.                                     posed by us. The characteristics of the GA approach to
                                                                 the trajectory optimization are examined for a simple
                                                                 problem with an analytic solution and the limitation
1    Introduction                                                of this approach will be revealed. Additionally, as a
                                                                 practical complicated example, an ascent trajectory op-
Trajectory optimization problems have been widely con-           timization problem of a spaceplane is computed by the
sidered in aerospace engineering fields(1) . The problem          presented method which uses the GA approach com-
is formulated as the calculus of variations where an ob-         bined with the gradient method. The space plane is a
jective function such as a flight time or required fuel is        future space transportation vehicle which takes off hor-
minimized while satisfying initial/final conditions and           izontally, reaches to a space station with a single stage,
path constrains for a flying vehicle. Although some               and lands horizontally(6) . In order to perform this mis-
techniques such as the dynamic programming of Bell-              sion, the vehicle changes several types of engine, e.g., an
man and the maximum principle of Pontryagin have                 air breathing engine and a rocket engine according to
been developed and used, the recent development of               its flight phase. The efficient accent trajectory design is
computer is promoting new algorithms(2) . Some of                strongly desired to maximize payload or minimize con-
the new methods transform the original problem into              sumed fuel.

                                                             1
2     Trajectory Optimization Prob-                                   The initial Np population of the independent vari-
                                                                ables is prepared for the continuous parameter GA at
      lem                                                       random. The generation is updated from the following
Defining the state variables x(t) and the control vari-          process:
ables u(t), the dynamic system is defined by the state           1) Selection of parents: m + 2 parents are initially se-
equations as                                                        lected at random.
                   d
                     x(t) = f (x, u, t)            (1)          2) Generation of children: the gravity point of the first
                  dt
The state and control variables are optimized to mini-              m + 1 parents is determined as xp . The deviation
mize the following performance index                                vectors of each parents are defined as
                                                                                                m+1
                         J[x, u, t]                      (2)                                  1 X
                                                                                    xp =            xi                     (10)
                                                                                             m+1 1
while satisfying the initial and final conditions defined
as                                                                                      d i = xi − xp                      (11)
            ψ I (x(tI )) = 0, ψ F (x(tF )) = 0      (3)             The normalized orthogonal base vectors to the de-
where tI and tF are initial and final time respectively.             viation vector di , (i = 1, ..., m) are prepared as
In many cases, during a flight path some constrains are                                   e1 , ..., en−m                    (12)
introduced as
                                                                    where n is the number of the individual variables
               G(x, u) = 0, H(x, u) ≤ 0                  (4)        p.
                                                                    Using the deviation vectors and the base vectors,
3     GA Approach                                                   A set of children is generated in the following equa-
                                                                    tion;
The trajectory optimization problem stated in the pre-                                      m               n−m
                                                                                            X               X
ceding section can be converted to a nonlinear optimiza-                     xc = xp ±            wi di ±         vi Dei   (13)
tion problem. The time interval from tI to tF is divided                                    i=1             i=1
into N elements to transform the control variables u(t)
as a function of time to a set of discrete variables. By de-        where D is the distance of the orthogonal compo-
noting the nodal values of time by ti (i = 1, 2, ...N + 1),         nent of the vector dm+2 to the d1 , ..., dm vectors,
control variables are represented as a set of discrete val-         and wi and vi are the random variables with nor-
ues ui = u(ti ). If the initial states x(tI ) are specified in       mal deviates. Note thet the number of chirdren is
the initial conditions, and if the final time tF is assumed          specified as 2 × Nc . The above process is called
and the descritized control variables ui are assumed,               as UNDX-m method(7) which allows the uniform
the state variables can be calculated by an appropri-               search.
ate integration scheme. In this paper, the control vari-        3) Selection: the fitness function of each parent and
ables between nodal values of ui are interpolated using             generated child is calculated. The fitness function
a third order spline function, and the 4th order Runge-             can be generated by using the performance index
Kutta algorithm is used as the integration scheme. This             and the constraints as
approximation makes it possible to represent the state                                X j
variables at time node ti as a function of descretized                    F = J+          [w1 max[0, | Gj (p) | −²]2
control variables ui and the final time tF . Therefore,                                  j
the independent variables to be optimized are defined                                j
                                                                                  +w2   max[0, Hj (p)]2 ]                  (14)
as
                                                                             j         j
               pT = [uT , uT , ..., uT +1 , tF ]
                         1   2       N                    (5)       where w1 , and w2 are weighting values in the
                                                                    penalty method which represents the constraints,
It should be noted that unknown initial states variables
                                                                    and ² is a small number.
can be incorporated into the independent variables.
       The performance index J, the final conditions at              The finest chromosome is selected. Another chro-
tF , and path constrains are defined as a function of the            mosome is chosen from the roulette selection which
independent variables p. The trajectory optimization                determines a fine chromosome at random according
problem can be transformed into the following nonlinear             to the specified probability based on the ranking.
optimization problem:                                           4) Determination of the new generation: According to
                                                                    the MGG (minimum generation gap) method(7) ,
                 minimize : J(p)                         (6)
                                                                    the selected two chromosome are exchanged with
                subject to : ψ F (p) = 0                 (7)        the two selected at random from the old parents.
                                G(p) = 0                 (8)        Hence the new parents are generated. If the finest
                                H(p) ≤ 0                 (9)        solution does not converge, return to the step 2).
4    Gradient Approach                                        where X k is the vector X at kth iteration step, ∇ indi-
                                                              cates the gradient of a function associated with the vari-
While the many gradient approach methods have been            ables X, and dk is the update vector (= X k+1 − X k ).
presented for the trajectory optimization problem, this       Vector G0 is composed from all the equal constraints.
paper uses the BDH (5) method that is one of the direct       Matrix B in Eq. (26) denotes the approximation of the
collocation methods. In the direct collocation method,        Hessian matrix H for the following Lagrange function
not only the control variables but also the state variables   L:
are decritized. The BDH is using linear interpolation                    L = J(X) + λT G0 (X) + µT H(X)            (29)
for this descritization. In the similar way as the GA
approach, the time is divided into N elements from the                                          H = ∇2 L                             (30)
initial time tI to the final time tF . The state and control   where λ and µ are called the Lagrange variables
variables are denoted by a set of nodal values as,            vectors. The SQP formulation assumes the initial
                 x1 , u1 , ...., xN +1 , uN +1        (15)    B to be a unit matrix and modifies the matrix B
                                                              through iteration steps by using gradient informa-
and the time interval of the element is given                 tion. The Broyden-Fletcher-Goldfarb-Shanno (BFGS)
                                                              formulation(8) is used to update the matrix B in the
                      ∆ti = ti+1 − ti                 (16)    following way:
      The variables X i (i = 1, 2, ..., N + 1) composed of
                                                                                           q k (q k )T   B k pk (pk )T B k
the discretized state and control variables and the time               B k+1 = B k +                   −                             (31)
interval are defined as                                                                     (q k )T pk
                                                                                                          (pk )T B k pk

                   X i = [xT , uT , ∆ti ]T
                           i    i                     (17)    where
                                                                                          pk = X k+1 − X k                           (32)
The independent variables X are obtained as,
                                                                   k       T         k+1        k+1        k+1        T      k   k   k
                                                                   q = ∇ L(X               ,λ  ) − ∇ L(X , λ , µ )
                                                                                                      ,µ
                  X = [X T , ..., X T +1 ]T
                         1          N                 (18)                                                         (33)
                                                                       The Hessian matrix for Eqs. (19)-(24) is derived
     Using this approximation, the trajectory opti-
                                                              as
mization problem is transformed into a nonlinear pro-
gramming problem defined as,                                                                                                 
                                                                                   ∇2 1 X1 L
                                                                                    X                             0
                                                                                                 ..                         
                      minimize J(X)                   (19)             H =                            .                            (34)
                                                                                      0                     ∇2 N +1 XN+1 L
                                                                                                             X
subject to
             G(X) = 0                                 (20)    Since the Hessian matrix has a diagonal form, each
                                 ∆ti                          block in the matrix B is updated by BFGS formula-
             xi + f (xi , ui )                                tion independently. The fundamental characteristics of
                                  2
                               ∆ti+1 Ki                       the BDH methods were examined in Reference (10).
             −[xi+1 + f (xi , ui )      ]=0           (21)
                                 2 Ki+1
             Ki+1 ∆ti − Ki ∆ti+1 = 0                  (22)    5        Numerical Examples
             ψ F (X) = 0                              (23)
             H(X) ≤ 0                                 (24)    5.1        Rocket Thrust Direction Control
                                                              To check the validity of the GA approach, a planar pow-
where Ki is a constant defined as,
                                                              ered flight was investigated. A point mass m is guided
                    Ki = ∆ti /(tF − tI )              (25)    by controlling the direction angle β of a constant thrust
                                                              with magnitude ma on a planar inertial coordinate sys-
Note that the Eq. (21) indicates the continuity condi-        tem as shown in Fig.1.
tion of states at the mid point in time element ∆ti .               The state equations are
      In order to solve the nonlinear programming                                                        
problem, a sequential quadratic programming (SQP)(8)                         0 0 1 0             
                                                                                                      0    
                                                                                                            
method is employed. The SQP tries to solve the fol-                d        0 0 0 1                      
                                                                                                       0
                                                                      x=   0 0 0 0 
                                                                                           x +                    (35)
lowing quadratic programming problem in a sequential               dt                             a cos β 
                                                                                                           
manner:                                                                                                    
                                                                             0 0 0 0                a sin β
                               1
          minimize ∇J(X k )dk + (dk )T Bdk            (26)                                x = [x, y u, v]T                           (36)
                               2
subject to                                                    where x and y are the inertia coordinates and u and v
                                                              are the particle’s velocity.
             G0 (X k ) + ∇G0 (X k )dk        =   0    (27)          Let’s consider the transfer of the particle to a spec-
             H 0 (X k ) + ∇H 0 (X k )dk      ≤   0    (28)    ified height h with the zero vertical velocity in a fixed
given time T so as to maximize the final horizontal ve-
                                                                                Table 1: Numerical Data
locity. The initial and final conditions and performance
index for this problem are                                                  take off weight m0         300 [ton]
           ∙               ¸         ½    ¾                                 wing area Sw              336.8 [m2 ]
              0 1 0 0                   h                                   ATR intake area           17.96 [m2 ]
     ψF =                    x(T ) −        =0     (37)
              0 0 0 1                   0                                   SCR intake area           17.96 [m2 ]
                                                                            RE thrust TRE             233.6 [ton]
                        J = −v(T )                      (38)
                                                                            RE Isp                    450 [sec]
Dividing the time scale into N = 20 elements, the inde-
pendent variables to be optimized are represented as a
set of βi .
                                                                of attack α. Note that the engine holds the max thrust.
      In this calculation, the constant values are a = 2,       The state equations are defined as
h = 20, and T = 10, respectively. The following values
are selected in the GA method:                                       dh
                                                                           = v sin γ                                  (40)
                                                                     dt
              Np = 100, m = 4, Nc = 100                 (39)         dv      T cos α − D            µ
                                                                           =             + (ω2 r − 2 ) sin γ        (41)
                                                                     dt            m                r
Figure 2 indicates the transition of the finest perfor-               dγ      T sin α + L    v     µ    ω2r
mance index during the GA process. Figures 3,4 and                         =             +( − 2 +           ) cos γ
                                                                     dt          mv         r    vr      v
5 are the horizontal velocity, the vertical velocity, and
                                                                             +2ω                                    (42)
the trajectory, respectively. Those values are compared
                                                                    dm           T
with analytical solutions. The numerical solutions show                    = −                                      (43)
fluctuation around the analytical ones. This is remark-               dt        Isp g0
able in Fig. 6 showing the control variables, i.e., the
                                                                where ω and g0 (= µ/R2 ) are the angular velocity of
                                                                                           0
time history of the thrust direction. These results in-         the earth and the gravity constant at the earth surface.
dicate that the GA approach has a great advantage in            Additionally, D and L are the lift and the drag. T and
the selection of initial solutions since those are selected
                                                                Is p are the engine thrust and the specific impulse which
at random but has a disadvantage in the convergence             are determined according to the Mach number M as
characteristics. It should be noted that the fluctua-
tion in the GA method can be reduced if the number                                    M ≤6        :   ATR             (44)
of initial population increases and the huge number of                           6 ≤ M ≤ Ms       :   SCR             (45)
generations is allowed.
       Figures 7 and 8 show the optimal solutions ob-                                  Ms ≤ M     :   RE              (46)
tained from the BDH approach in which the initial solu-                The aerodynamic characteristics of the spaceplane
tions are obtained from the GA approach. These results          are determined from Reference (6) and the ATR and
indicate the good accuracy of the BDH method.                   SCR engine characteristics are obtained from Reference
                                                                (10). The other data used are listed in the following
5.2    Space Plane Accent Trajectory Op-                        table.
       timization                                                      Initial conditions are specified to consider take-off
                                                                conditions at time tI = 0 as
As a practical and complicated trajectory optimization
problem, an accent trajectory of a spaceplane (Fig. 9)                    h(0) = 0 [km], v(0) = 150 [m/s]             (47)
is investigated. The space plane is planned as a fu-                      γ(0) = 0 [deg], m(0) = 300 [ton]            (48)
ture space transportation vehicle which takes off hori-                                          dγ
zontally, accents to the space station directly, and de-                                           (0) ≥ 0            (49)
                                                                                                dt
scends to the earth to land horizontally. In order to
perform this mission in a single stage, the plane is de-        The Mach number Ms which indicates the switching
signed to equip with several types of engines to be se-         timing from SCR engine to RE engine is specified as
lected according to the flight environment. A typical                                     Ms ≤ 12                      (50)
accent flight path is shown in Fig. 10 where the space
plane rises above 90 km by using an air-turboramjet             The terminal conditions at t = tF when the RE is cut-
engine ATR, a scramjet engine SCR, and a rocket en-             off are specified as
gine RE sequentially. The RE is cut-off and the vehicle
zooms up to 400 km in the elliptic orbit. At the apogee                                     h(tF ) ≥       90 [km]    (51)
in the elliptic orbit, the vehicle is put on a circular orbit                                 γ(tF ) ≥     0 [deg]    (52)
as shown in Fig.11.                                                       Ha [h(tF ), v(tF ), γ(tF )] =    400 [km]   (53)
      As shown in Fig. 12, the state variables are the
altitude h (= r −R0 ), velocity v, flight-path angle γ and       where Ha is the function which defines the apogee alti-
weight m. The control variable is defined as the angle           tude on the elliptic orbit. Several path constraints must
be considered in the following items:                         Reference

             altitude :    h ≥ 0 [km]                (54)     (1) Bryson, A. E., Jr., and Ho, Y. C., “Applied Opti-
                               1 2                                mal Control”, Blaisdell Publishing Company, MA,
    dynamic pressure :     q = ρv ≤ 100 [kPa]        (55)         (1969).
                               2
      angle ofattack :     α ≤ 20 [deg]              (56)
                                                              (2) Hargraves, C. R. and Paris, S. W., “Direct Trajec-
        heating rate :     Q˙ = 4.4575 × 10−8 √ρv 3.07
                                                                  tory Optimization Using Nonlinear Programming
                                           2                      and Collocation”, Journal of Guidance, Control,
                           ≤ 300 [kW/m ]               (57)
                           L cos α + D sin α                      and Dynamics, Vol. 10, No. 4 (1987), pp.338-342.
          load factor :                      ≤4        (58)
                                  mg0                         (3) Holland, J., “Adaptation and Natural and Artificial
     The objective function to be minimized is the final           System , Univ. of Michigan Press, Ann Arbor, MI,
mass at the circular orbit and is given as.                       1975.

                        ∆v[h(tF ), v(tF ), γ(tF )]            (4) Wright, A. “Genetic Algorithms for Real Parame-
     J = −m(tF )exp(−                              )   (59)       ter Optimization, in G.J.E. Rawlins (Ed.), Foun-
                                 g0 Isp
                                                                  dations of Genetic Algorithms, San Matero, CA:
where ∆v is the velocity change required to put the               Morgan Kaufmann, pp. 205-218.
veihcle on the circular orbit.
      As the same way in the previous example, the GA         (5) Tsuchiya, T., and Suzuki, S., “Spaceplane Tra-
approach is applied to find initial solutions for the BDH          jectory Optimization with Vehicle Size Analysis,
method. In the GA calculation, the following data are             14th IFAC Symposium on Automatic Control in
used:                                                             Aerospace, 1998, pp.444-449.
             Np = 500, m = 10, Nc = 100               (60)
                                                              (6) Nomura, S., Hozumi, K., Kawamoto, I., and
       Figures 13-16 show the comparison between the              Miyamoto, Y., “Experimental Studies on Aero-
GA results after 10,000th generation and the BDH re-              dynamic Characteristics of SSTO Vehicle at Sub-
sults in which the time scale is divided into N = 80              sonic to Hypersonic Speeds”, The 16th Interna-
elements. Note that the initial solutions in the BDH              tional Symposium on Space Technology and Sci-
are obtained from the GA method. Figure 13, 14, and               ence, (1988) pp.1547-1554.
15 are the time histories of the velocity, the altitude and
the control variables, respectively. It should be noted       (7) Kita, H., Ono, H., and Kobayashi, S, “A New Gen-
that the BDH can refine the solutions. The final mass               eration Alternation Method of Genetic Algorithms
increases from 75.93 ton to 80.34 ton and the terminal            and Its Assessment (in Japanese), ” Journal of Ar-
time tF is reduced from 1387.8 sec to 835.9 sec. Figure           tificial Inteligence, Vol. 12, No. 5 (1997), pp.734-
16 shows the v-h diagram which shows the relation the             744.
velocity and the altitude. This clearly indicates the re-
                                                              (8) Ibaragi, H. and Fukushima, M., “FORTRAN77
finement of the trajectory obtained by the BDH which
                                                                  Optimization Programming”, Iwanami Shoten,
strictly satisfies the constraints of the dynamic pressure
                                                                  Tokyo, (1991)
and the heating rate.
                                                              (9) Tsuchiya, T., and Suzuki, S., “Computational
                                                                  method of Optimal Control Problem Using Math-
6     Conclusion                                                  ematical Programming (2nd Report), Introduction
                                                                  of Block Hessian Method (in Japanese),” Journal of
The numerical method for the flight trajectory opti-
                                                                  Japan Society for Aeronautical and Space Sciences,
mization problem is investigated in this paper. The
                                                                  Vol. 46, No.536 (1998), pp.497-503.
new approach is combined the Genetic Algorithm (GA)
method and the BDH method which is one of the gra-            (10) Kato, K., “Spaceplane (in Japanese)”, University
dient approach. The GA approach is used to find initial            of Tokyo, Tokyo, (1989).
solutions for the BDH method. While the GA approach
can select initial solution at random, it is difficult to ob-
tain the smooth converged solution. On the other hand,
the BDH method has the good convergence character-
istics if the good initial solutions are specified. The
proposing approach utilizes the good characteristic in
each method. The simple example which has an ana-
lytic solution is used to check the validity of the method,
and the complicated example of the accent flight tra-
jectory optimization of a spaceplane is demonstrated to
show the applicability of the proposing approach.
          Figure 1: Coordinate system




                                                        Figure 5: Flight trajectory




Figure 2: History of fitness function vs generation




                                                     Figure 6: Time history of control




  Figure 3: Time history of horizontal velocity




                                                        Figure 7: Flight trajectory

   Figure 4: Time History of vertical velocity
   Figure 8: Time history of control         Figure 11: Orbit transfer of spaceplane




         Figure 9: Spaceplane                 Figure 12: State and control variables




Figure 10: Spaceplane accent trajectory   Figure 13: Time history of spaceplane altitude
    Figure 14: Time history of spaceplane velocity




    Figure 15: Time history of spaceplane control




Figure 16: H-V diagram of spaceplane accent trajectory