Spaceplane Trajectory Optimization with Vehicle Size Analysis

Copyright © IF AC Automatic Control in Aerospace, Seoul, Korea, 1998

SPACEPIANE TRAJECTORY OPTIMIZATION WITH VEmCLE SIZE ANALYSIS

Takeshi Tsuchiya, ShiIUi Suzuki

Department ofAeronautics and Astronautics, School ofEngineering, the University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan

Abstract: A numerical method to solve simultaneous design and control optimization problems is considered. Firstly, the simultaneous optimization problem in this paper is transformed into a new formulated nonlinear programming problem. Considering using a sequential quadratic programming method, the diagonal form of the Hessian matrix for the Lagrange function is introduced in a new formulation named Block Diagonal Hessian (BDH) method. Finally, the present method is applied to simultaneous spaceplane design and trajectory optimization problems for feasibility study of a future space transportation system. Through several calculations, the proposed BDH method is confirmed to be effective for simultaneous optimization problems. Copyright © 1998 IFAC Keywords: Numerical methods, Optimization problems, Nonlinear programming, Optimal trajectory, Space vehicles

flight trajectory optimization with dynamic variables can be formulated as optimal control problems which have been solved by using numerous solution approaches; e.g., calculus of variations, maximum principle, and dynamic programming (Bryson, 1969). Recent investigations (Goh and Teo, 1988; Suzuki and Yoshizawa, 1994; Hargraves and Paris; 1987) have revealed the advantage of using a nonlinear programming for optimal control problems. One reason for this is that computer capacities has been increased to deal with a large number of variables, and that efficient computation algorithms have been developed to deal with various constraints, such as differential equations and inequality constraints.

1. INTRODUCTION In a mechanical design process, both design parameters and control parameters have to be simultaneously optimized to obtain an extremely high performance. Simultaneous design and control optimization, however, is difficult, because both are in different technical fields. In a conventional aircraft design process, a configuration and geometrical parameters of an aircraft are optimized for assumed flight patterns, and, after an aircraft has been designed, flight paths or trajectories are optimized.

Numerical optimization techniques have been developed for both design and trajectory problems. Design optimization problems with static variables are solved by using a nonlinear programming categorized in a mathematical programming. Recent many papers report that design parameters, such as wing area, aspect ratio, total weight or length, and so on have been successfully optimized in an aircraft design field (Asbley, 1976; Grossman, 1988; Dovi, 1990). On the other hand, the

In this paper, a new nonlinear programming formulation is proposed to solve simultaneous design and control problems, and the numerical method is applied to spaceplane design optimization. Spaceplane is a future space transportation system from the ground to orbit. While many concepts have been proposed in the world (Mankins, 1998), this paper considers a single stage to

439

orbit (SSTO) reusable plane which takes off and lands horizontally and uses several engines operating in different environments. Both the design parameters and its flight trajectory are optirnized at the same time.

2. OPTIMIZATION TECHNIQUE

2.1 Simultaneous Optimization Problem The simultaneous problem in this study is defined as follows. The time tE[O,t r ] is a scalar, the state variables x(t) is an n -vector, the control variables u(t) is an m -vector, and design variables d is [vector. The values of X(t) , u(t) and d are optimized to minimize the following performance index;

J

= tP[X(lr ),u(t r ),d]

Fig. 1. Discretized state Xi !!!!

x(ti )

(8)

Ui

U(ti )

(9)

E

and the time interval of the element is given

(1)

(10)

The dynamics of the system are defined by the state equations, initial conditions and terminal conditions

x = f[x(t),U(l),d]

(2)

({Jo[ x(O), u(O),d] = 0

(3)

({Jr[x(tr ),U(lr ),d] = 0

(4)

Supposing that the states on each element are approximated by two linear functions, as shown in Fig. 1, the following equations must be satisfied.

(11)

In addition, equality and inequality constraints of the form

e[ x(t),u(t),d] = 0

(5)

S[X(l),U(l),d]s 0

(6)

where Ki is a constant defined as,

(12)

must be satisfied.

which is introduced to obtain a block diagonal form of a Hessian matrix described in the following subsection. The variables Xi (i = 0,1, ... ,N) composed of the discretized state and control variables and the time interval are defined as

2.2 Nonlinear Programming Formulation DCNLP (Direct Collocation with Nonlinear Programming) method (Hargraves and Paris, 1987) is widely used as a numerical optimization method to solve the problem (Eqs. (1}-(6)). Here, a new formulation is proposed to improve computing load and convergence characteristics in DCNLP method.

(13)

To transform the optimal control problem (Eqs. (1}-(6)) into a nonlinear programming problem, the time t from initial time 0 to final time lr is divided into N elements. The points connecting the elements are called nodes and ti (i = O,l, ... ,N) denotes the time on each node. Note that

And then these variables and the design variables are combined to provide the following variables X used in a nonlinear programming problem:

X

•

(X 0

T

, ••• ,

X

N

T

,

d)T ER

(N +lXn+m+l)+1

(12)

(7) The simultaneous problem described above is transformed into a nonlinear programming problem defined as,

On each node the state and control variables depending on time are discretized as

440

minimize:

(14)

subject to: where p(k) = X(k+l) _ X(k)

(15)

q(k)

(23)

= VTL(X(k+I) ,A(k+!), ",(hi» _ VT L(X(kl,A (hi), ",(hi»

(24)

Keeping this in mind, the Hessian matrix for Eqs. (14)(16) is derived as

Functions F, G and H defined in Eqs. (14)-{16) are called an objective function, equality constraints and inequality constraints, respectively.

V xx 2 L 00

(25)

2.3 Block Diagonal Hessian Method A large amount of literature has been devoted to the solution method of nonlinear programming problems. A sequential quadratic programming (SQP) method (Ibaragi and Fukushima, 1991) is said to be most useful at present time. It follows from this that the SQP method is applied to solve the nonlinear programming problem defined above.

What is important in Eq. (25) is that, if the Hessian matrix is divided into four parts, the part in the upper left of the matrix becomes a block diagonal form. However, it is not advisable to directly apply the BFGS formulation to Eq. (25) since the non-diagonal parts in Eq. (25) require a huge computer memory and worsen the convergence characteristics of the iteration process.

This method finds the optimal solution by solving the following quadratic programming problem in a sequential manner:

minimize: VF(X(k»e(k) +.!.e(k)TB(k)e(k)

2

(17)

=0

(18)

H(X(k»+VH(X(k»e(k) sO

(19)

subject to: G(X(k» + VG(X(k»e(k )

Consequently the following method is adopted. The matrix B approximating Eq. (26) is introduced as

B=

where X(k) is the vector X at kth iteration step, V denotes the gradient of a function associated with the variables X, and e(k) is the update vector ( = X(k+I) - X(k». Matrix B in Eq. (17) denotes the approximation of the Hessian matrix H for the following Lagrange function L: L(X ,A,,,,) .,. F(X) + ATG(X) + JI.TH(X)

(20)

H. V2L

(21)

(26)

BI ER(n+m+l)x(n+m+l) (i = 0,1, . .. ,N )

(27)

B N+I ERlxl

(28)

The total Hessian matrix is approximated to a block diagonal one, where Bi and BN+I represent the approximated matrix for VXiX/L and Vdd2L respectively. As a result, the each block matrix can be updated by BFGS formulation independently. The presented method keeps sparse characteristics in the approximated Hessian matrix, and it is possible to save the computer storage by omitting the null parts in the matrix.

where A and '" are called the Lagrange variables vectors. The reason why the approximated matrix B is substituted for the Hessian matrix H is that it is difficult to compute the Hessian matrix actually. The SQP formulation assumes the initial B to be a unit matrix and modifies the matrix B through iteration steps by using gradient information. The B:oyde~­ Fletcher-Goldfarb-Shanno (BFGS) formulatIOn IS usually used in the following way:

The formulation Eqs. (14)-{16) and the method described above is named block diagonal Hessian (BDH) method (fsuchiya and Suzuki, 1997).

441

3. SPACEPLANE DESIGN STUDY

In this section, the BDH method is applied to simultaneous sizing and trajectory optimization for a future space transportation vehicle; i.e., a spaceplane. While the design phase and flight control phase for spaceplane will be separately explained to understand the outline of the problem easily, the two phases are simultaneously optimized by the proposed method.

v

p : the air density r

3.1 Spaceplane Design Phase Fig. 3. Symbols

dh . -=vsmy dt

(29)

dv = (TATR + TSCR + T ROC ) cos a -D ~siny dt m --r-2-

(30)

dy dt

-=

(TATR

+ TSCR + T ROC ) sin a +L mv

(V

~) vr2

+ - - - cosy r

(31)

Fig. 2. NAL spaceplane model The NAL model shown in Fig. 2 was adopted as a original model. This aerodynamic characteristic was experimentally studied from low speed to hypersonic speed by Nomura, et al. (1988).

where L and D are the lift and drag respectively, T ATR , TSCR and T ROC are the thrust of each engine, /SPATR' /SPscR and /SPROC are the specific impulse, ~ is the 2 gravity constant and go (= f.J./ Ro ) is the gravity acceleration at the ground level. The lift L and drag D are represented by the lift and drag coefficient, dynamic pressure q, wing area Sw and fuselage surface area. The thrust of ATR and SCR is represented by Kato (1989), where TATR and TSCR are proportional to S ATR and SSCR respectively. It should be noted that the motion equations change discontinuously since the operating engines are switched according to the flight conditions.

In the design phase, the design variables are the takeoff weight mo, wing area Sw, body length I, intake area of air-turboramjet engine (ATR) S ATR ' intake area of scramjet engine (SCR) SSCR' thrust of rocket engine (ROC) T ROC' maximum dynamic pressure q max , and maximum load factor nmax (:s 4 [G]). Note that the total fuel volume composed of a liquid hydrogen (LHz) and a liquid oxygen (LOX) must be less than the tank volume defined by the body size.

3.2 Flight Control Phase Initial conditions at time t Supposing the above design, in the flight control phase, the spaceplane rises above 90 km and is accelerated by ATR (to Mach 6), SCR (switched from ATR and useable to Mach 12) and ROC (useable with ATR and SCR), and then, after the engine is cut-off, it zooms up to 400 km with no thrust in an elliptical orbit. Finally, it is put into a 400 km circular orbit at the apogee in the elliptical orbit.

=0

are specified as

h(O)

=

0

[km]

(33)

y(O)

=0

[deg]

(34)

m(O) =mo

(35)

and the physical requirements at a runway are introduced as

The state variables are the altitude h ( = r - Ro , see Fig. 3), velocity v, flight-path angle y and weight m. The control variable is defined as the angle of attack a. The motion equations (Kato, 1989) of a spaceplane are expressed as

(36)

m v2 -:s1500[m]

_0

(37)

2TATR

Terminal conditions at the engine cut-off time t = t f is

442

3.4 Nwnerical Results

expressed as h(tf) ~ 90

[km]

(38)

y(tf) ~ 0 [deg]

(39)

H.[h(t f ), v(tf )' y(tf

)]= 400 [km]

The numerical solutions are shown in Fig. 4--6 and Table 1, where takeoff weight is fixed to 100 ton. The characteristics of the original design are summarized in Table 1. The flight trajectory of the original is calculated to minimize the fuel weight. The true required weights estimated by WATTS is 145.8 ton, therefore the performance index is 1.458. This indicates that it is necessary to decrease the true required weight more than 30 %.

(40)

where Ha is the function which defines the apogee altitude by using the cut-off states. In addition, the following path constraints are defined, which mean altitude limit, dynamic pressure limit, angle limit of attack, and load factor limit composed of the aerodynamic force and the thrust force. h~O

1

[km]

(41)

2

(42)

q="2{JV sqmax a slO

The table and figures also show the simultaneous optimal solution with the takeoff weight 100 ton. The minimized performance index value is 1.185. This results is approximately 0.273 smaller than that of the original, which shows the required weight to be reduced decreases drastically. Table 1 and Fig. 4 indicate that the vehicle size of the simultaneous solution is smaller than the original. The optimized wing area and the intake area of ATR is much smaller, and the intake area of SCR becomes to be zero, which means that SCR is unnecessary in this calculation. It can be considered from the simultaneous solution that the wing area and ATR respectively decrease to the limit size in order to take off and fly the vehicle against the aerodynamic drag, and that the volume of l l i2 is reduced because SCR isn't used. Another calculation revealed that SCR is useless, even if I SPSCR becomes twice as large.

[deg]

(43)

n s n max

(44)

3.3 Optimization After the design and flight trajectory are specified, the true required weight m TO is defined as the sum of the consumed fuel weight and the true vehicle weight. The true vehicle weight is estimated from the spaceplane size by using W AATS (Glatt, 1974), which is a computer program for statistical weight analysis of advanced transportation system. The coefficients are modified in parts (Shirouzu, 1989) because the original data are old. If the true required weight m TO is more than the takeoff weight mo, the spaceplane is unrealizable. For that reason, the performance index to be minimized is defined as the ratio between the both weights.

J

Although SCR and ROC are used at the same time in the flight by the original design, as for the simultaneous solution ATR and ROC isn't used at the same time. This can be understood from that the use of both engines at the same time causes the increase in the weight of the structure supporting the engines. Finally, the change of the takeoff weight which has been fixed to 100 ton is investigated. Figure 7 shows the performance index values for the various takeoff weights mo up to 500 ton. The performance index value slightly decreases with the takeoff weight. This indicates that the increase of the takeoff weight cannot improve the performance index greatly, and the weight reduction is the essential key technology to realize the future space transportation system.

(45)

= Tnro/mo

Iabl~] ChaIact~risti~ Qfspa~plaD~

Characteristics Body length [m1 Wing area [m2] Wing span [m] Intake area of ATR [m2] Intake area of SCR [m2] Thrust ofROC [ton] Max. thrust [ton] Max. dynamic pressure [kPa] Max. load factor [G] Takeoff weight [ton] Weight on orbit [ton] True reguired weight [ton] Performance index

Original Design 50.94 337.3

23.84 14.00 14.00 70.00 321.9 100.0 2.982 100.0 24.39 145.8 1.458

Optimal Design 45.89 96.63 15.85 6.095 0.000 77.63 77.63 56.82 4.000 100.0 18.78 118.5 1.185

Fig. 4. Spaceplane of 100 ton takeoff weight

443

100 Original Shape90 Optimal Shape--1/ 80 / .l" 70 ,l/ E 60 ~ 50 Am Roe C1.l '"0 40 .~ /------_/ 30 ,,,--._,,,, ~ " 20 // 10 ._____-:.\m Roe 00 200 400 600 800 1000 Time t[sec]

=-

10.0

C1.l

'"0

'E' 6.0 ..>0:

u

c
~ 4.0

....

~ 2.0

\

\.r"

Original ShapeOptimal Shape---

<:

-2.0

0

200

400 600 Time t[sec]

800

1000

Fig. 6. Time history of angle of attack Multiobjective Optimization Methods. Journal of Aircraft,27, 1043-1049. . Glatt, G. R. (1974). WAATS-A Computer Program for Weights Analysis of Advanced Transportation Systems, NASA CR-2420. Goh, C. J. and K L. Teo (1988). Control Parametrization: a Unified Approach to Optimal Control Problems with General Constraints, Automatica, 24, 3-18. Grossman, B., G. J. Strauch, W. M. Eppard, Z. Guardal, and R. T. Haftka (1988). Combined AerodynamicStructural Design of a Sailplane Wing. Journal of Aircraft, 25, 855-860. Hargraves, C. R. and S. W. Paris (1987). Direct Trajectory Optimization Using Nonlinear Programming and Collocation. Journal of Guidance, Contro~ and Dynamics, 10, 338-342. Ibaragi, H. and M. Fukushima (1991). FORTRAN77 Optimization Programming. Iwanami Shoten, Tokyo (in Japanese). Kato, K (1989). Spaceplane. University of Tokyo Press, Tokyo (in Japanese). Mankins, J. C. (1998). Lower Costs for Highly Reusable Space Vehicles. Aerospace America, March. 26-42. Nomura, S., K Hozumi, I. Kawamoto and Y. Miyamoto (1988). Experimental Studies on Aerodynamic Characteristics of SSTO Vehicle at Subsonic to Hypersonic Speeds. The 16th International Symposium on Space Technology and Science, 15471554. Shirouzu, M. (1989). Weight Estimation and Sensitivity Analysis of Spaceplane (part 1) SSTO 1st Report, Technical Memorandum of NAL, TM-598 (in Japanese). Suzuki, S. and T . Yoshizawa (1994). Multiobjective Trajectory Optimization by Goal Programming with Fuzzy Decisions. Journal of Guidance, Control, and Dynamics, 17, 297-303. Tsuchiya, T. and S. Suzuki (1997). A Study on Simultaneous Designffrajectory Optimization of Spaceplane. International Symposium on Optimization and Innovative Design, Paper#158.

lSOr~--------~----------~

><

\

00 c

Fig. 5. Time history of altitude

C1.l

\

00 8 .0

;f

..., 1

\,

Original Shape

/

'"0

.:: 1 1:l c

c
§ 1 .8

cE 100 Fig. 7. Performance index values for takeoff weights

4. CONCLUSION The design and control process of a system is integrated in a simultaneous optimization formulation. A huge number of design and control variables must be dealt with to solve simultaneous optimization problems. This paper proposes the nonlinear programming formulation called BDH (Block Diagonal Hessian) method in which the Hessian matrix for the Lagrange function is the block diagonal form. The simultaneous optimization problem for spaceplane design and its flight trajectory is studied for the design of the future space transportation system. Through this study, it was confirmed that the proposed BDH method was effective for the complex simultaneous optimization problems.

REFERENCES Ashley, H. (1976). On Making Things the BestAeronautical Uses of Optimization. Journal of Aircraft, 19, 5-28. Bryson, A E., Jr. and Y. C. Ho (1969). Applied Optimal Control. Blaisdell Publishing Company, MA Dovi, A R. and G. A Wrenn (1990). Aircraft Design for Mission Performance Using Nonlinear

444