Get documents about "Flight Trajectory Optimization Using Genetic Algorithm Combined
Shared by: bfb53718
Categories
Tags
genetic algorithm, genetic algorithms, trajectory optimization, evolutionary algorithms, guidance law, evolutionary computation, aerospace engineering, international conference, multi-objective optimization, multiobjective optimization, fitness function, design variables, optimization problems, kalyanmoy deb, objective function
-
Stats
- views:
- 28
- posted:
- 2/12/2010
- language:
- English
- pages:
- 8
Document Sample
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
Related docs
Other docs by bfb53718
