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 eﬃcient ﬂight 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 speciﬁed 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 diﬃcult to ﬁnd 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-ﬁttest dom search methods, the initial solutions are selected strategy to eliminate unﬁt 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 ﬂuctuation in the solutions. Therefore, the continuous variables can be dealt with by an advanced gradient approach is utilized to reﬁne the solutions. As continuous parameter GA(4) which deﬁnes a chromo- the gradient method, the authors uses the BDH method some as an array of ﬂoating 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 eﬃcient 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 ﬁelds(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 ﬂight time or required fuel is future space transportation vehicle which takes oﬀ hor- minimized while satisfying initial/ﬁnal conditions and izontally, reaches to a space station with a single stage, path constrains for a ﬂying 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 ﬂight phase. The eﬃcient 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 Deﬁning the state variables x(t) and the control vari- process: ables u(t), the dynamic system is deﬁned 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 ﬁrst 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 deﬁned as m+1 J[x, u, t] (2) 1 X xp = xi (10) m+1 1 while satisfying the initial and ﬁnal conditions deﬁned 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 ﬁnal time respectively. viation vector di , (i = 1, ..., m) are prepared as In many cases, during a ﬂight 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 speciﬁed in mal deviates. Note thet the number of chirdren is the initial conditions, and if the ﬁnal time tF is assumed speciﬁed 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 ﬁtness function of each parent and ables between nodal values of ui are interpolated using generated child is calculated. The ﬁtness 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 ﬁnal time tF . Therefore, j the independent variables to be optimized are deﬁned 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 ﬁnal conditions at The ﬁnest chromosome is selected. Another chro- tF , and path constrains are deﬁned as a function of the mosome is chosen from the roulette selection which independent variables p. The trajectory optimization determines a ﬁne chromosome at random according problem can be transformed into the following nonlinear to the speciﬁed 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 ﬁnest 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 ﬁnal 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 modiﬁes 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 deﬁned 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 deﬁned 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 deﬁned as, ered ﬂight 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) iﬁed height h with the zero vertical velocity in a ﬁxed given time T so as to maximize the ﬁnal horizontal ve- Table 1: Numerical Data locity. The initial and ﬁnal conditions and performance index for this problem are take oﬀ 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 deﬁned 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 ﬁnest 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) ﬂuctuation 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 speciﬁc 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 ﬂuctua- 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 speciﬁed to consider take-oﬀ 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 oﬀ 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 speciﬁed as lected according to the ﬂight environment. A typical Ms ≤ 12 (50) accent ﬂight 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- oﬀ are speciﬁed as gine RE sequentially. The RE is cut-oﬀ 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, ﬂight-path angle γ and where Ha is the function which deﬁnes the apogee alti- weight m. The control variable is deﬁned 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 Artiﬁcial The objective function to be minimized is the ﬁnal 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 ﬁnd 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 reﬁne the solutions. The ﬁnal 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 tiﬁcial 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 ﬁnement of the trajectory obtained by the BDH which Optimization Programming”, Iwanami Shoten, strictly satisﬁes 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 ﬂight 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 ﬁnd initial of Tokyo, Tokyo, (1989). solutions for the BDH method. While the GA approach can select initial solution at random, it is diﬃcult 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 speciﬁed. 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 ﬂight 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 ﬁtness 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

DOCUMENT INFO

Shared By:

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 |

OTHER DOCS BY bfb53718

How are you planning on using Docstoc?
BUSINESS
PERSONAL

By registering with docstoc.com you agree to our
privacy policy and
terms of service, and to receive content and offer notifications.

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.