- Email: [email protected]
Launch vehicle trajectory optimization
Acta Astronautica Vol. 15, No. 11, pp. 845-850, 1987 Printed in Great Britain. All rights reserved
0094-5767/87 $3.00 + 0.00 © 1987 Pergamon Journals Ltd
LAUNCH VEHICLE TRAJECTORY OPTIMIZATION# V. ADIMURTHY Vikram Sarabhai Space Centre, Trivandrum~595022, India (Received 18 December 1986)
A~traet--Three dimensional trajectory optimization problem is solved through a diagonalized multiplier method in which the multiplier update formula of Tapia and Han is employed over a BFGS Hessian update. The diagonalized multiplier method (DMM) is a single sequence minimization technique for nonlinear programming and avoids the sequence of minimization problems inherent in the classical augmented Lagrangian methods. To avoid violent initial multiplier corrections in DMM, a new modification is proposed in which we minimize the Kuhn Tucker vector norm along the Han-Tapia multiplier update direction. A software package PYOPT is developed and validated at the Indian Space Research Organisation (ISRO) employing the concepts of diagonalised multiplier method for constrained pitch and yaw steering optimization of launch vehicle trajectories. The system dynamic model for trajectory optimization allows for effective and routine operational studies. The optimality of the solutions of PYOPT is established through sufficiency conditions of optimality in nonlinear programming via saddle point optimality theorem and augmentability theorem. The robustness of the package is demonstrated. The present package, which is faster than programmes based on techniques like projected gradient methods and recursive quadratic programming, has become a standard tool at ISRO for mission design studies.
1. INTRODUCTION
Trajectory Optimization with detailed launch vehicle and mission constraints is a complex nonlinear problem, and efforts to solve this effectively using the state of the art continue to be made. The reason for considerable distance between theory and practice and for the sparseness of reported papers in this area depends on the difficulty in establishing good convergence in numerical iterative procedures. Brusch[1] addresses to the trajectory optimization for Atlas/Centaur launch vehicle employing a classical augmented Lagrangian method, in which the solution is sought by alternately minimizing the augmented function with respect to independent variables and adjusting the Lagrange multipliers to satisfy constraints. POST, a generalized Program to Optimize Stimulated Trajectories[2], employs an accelerated projected gradient method. These approaches involve considerable computer time for convergence. Well and Tandon[3] effectively use the recursive quadratic programming approach to three dimensional trajectory optimization. They find that this method seems to be a robust and relatively fast algorithm for solving parameterized optimal control problems. But they conclude that further modifications of the algorithm have to be done in order to reduce CPU-time such that actual mission design studies can be carried out effectively. However, none of these papers address to establishing, theoretically or otherwise, the actual optimality of the solutions. tPaper IAF-86-235 presented at the 37th Congress of the International Astronautical Federation, Innsbruck, Austria, 4-11 October 1986
In the present paper, we solve the threedimensional trajectory optimization problem through a diagonalized multiplier method in which multiplier update formula of Tapia[4] and Han[5] is employed over a BFGS Hessian update. Like recursive quadratic programming (RQP), the diagonalized multiplier method (DMM) is a single sequence minimization technique for nonlinear programming and avoids the sequence of minimization problems inherent in classical augmented Lagrangian methods. This leads to the substantially faster performance of D M M and RQP compared to classical methods, as found by Bertocchi et al.[6], who make a comparative study of the DMM, RQP and classical methods. It is also found that the performance of the diagonalized multiplier methods is better in several case studies. A software package PYOPT is developed and validated at the Indian Space Research Organization (ISRO) employing concepts of diagonalized multiplier method for constrained pitch and yaw trajectory optimization. To avoid violent initial multiplier corrections, a new algorthm is proposed in which we minimize the Kuhn-Tucker vector norm along the Han-Tapia multiplier update directions. The optimality of the solutions of PYOPT is established through sufficiency conditions of optimality in nonlinear programming via saddle point optimality theorem and augmentability theorem. The system dynamics aspects of the trajectory optimization are modelled to allow for effective and routine operational studies. This paper describes the system modelling, algorithmic aspects, and studies made to establish optimality and robustness. The present method is faster than existing software packages like POST and those based on recursive quadratic pro845
846
V. AD1MURTHY
gramming, and has become a standard tool at-ISRO for mission design studies. 2. MATHEMATICAL FORMULATION
2.1. General statement of the problem Let vector s denote the state (say, position and velocity) of the dynamical system under consideration (launch vehicle in motion). Let the dynamics of the system depend on a vector u of control variables (say, pitch, yaw and other design parameters), which can be specified at our will (assuming that specified pitch and yaw histories can be realised through the control system). The dynamics of the system can now be represented by ds ~=F(s,u,t),
to<~t<~ O.
(1)
We are interested in maximizing a performance index (payload, injection velocity, range) which depends, among other things, on s(tj). Since this depends on the choice of u through the dynamical equations of the system, the performance index also depends on u. The problem is to select the control vector u, such that the performance index is maximized, and at the same time various specified objectives and constraints are satisfied. The objectives can be the final orbit size and inclination and the like. The constraints can be on the vehicle loading, flight safety, etc. To solve this problem one should be equipped with (a) a capability to simulate the system dynamics, (b) a procedure to ensure that the specified constraints are satisfied and (c) a methodology to vary the control vector so as to maximize the performance index. In the scheme adopted in the present program PYOPT, an attempt has been made to model some of the constraints in the system dynamics itself, while the rest are left for optimization algorithm. The system dynamics is considered as a black-box that satisfies these constraints and gives out the performance index, once a set of control variables is specified. This scheme is equivalent to making a transformation on the dynamics in order to make a highly constrained problem a much less constrained one. The system-black-box concept leads to a direct functional relationship between the performance and the control vector. Further, it is assumed that the control variables have been parametrized, so that the problem is to find a finite number of control parameters instead of a continuous control function. With this parametrization, the performance index becomes a function of a finite number of control parameters. The optimization of this function is done using methods of constrained nonlinear optimization. The approach described above differs from that of other workers like Brusch[l] in one essential feature. Our system dynamics for optimization is not simply an integration of a system of basic differential
equations. Various in-flight constraints have been embedded in the system dynamics. This feature of the present study has necessitated in a careful modelling of the transformed system dyanamics and the flight profile.
2.2. Dynamical equations The basic coordinate system, in which all the computations in PYOPT are carried out, is a Geocentric Inertial System (x,y,z). In this frame, the x-y plane is the equatorial plane. The z-axis is pointed towa~'ds the north pole. x-z plane contains the launch meridian frozen at the instant of initial stage ignition (at t = 0). The equations of motion in this frame can be written as: __d2_--T ~x +
dt
~.~cos ~
m
m
ga
x dS'
_
dt
3 a4
T ~y + QS ~y cos ~
Vay
m
(
~aa}
m-
( Xce--rc~ QS¢,.C ~ × Cs~ 1+ Xc~ / - m " - - g 0Y- [ ~ r d~
-
=
dt
-
T 4: +
a4 _ 5zZ']] + 53Jz~-g(1 rZ ./3
(3)
~, cos
m
Va J
z
[-a 2
dx
3 dy
dz
d~=.~, h5=?, ~5=_~
(5)
The forces modelled in the above differential equations are (a) thrust (corrected to local pressure of the atmosphere), (b) aerodynamic normal force, (c) ideal control force to balance aerodynamic moments, (d) aerodynamic drag and (e) gravity force for an oblate earth (up to second harmonic only).
2.3. Vehicle attitude representation In PYOPT, the vehicle attitude information is basically recognised through the direction cosines ~x, ~y and ~: of the attitude vector ~ with respect to the basic coordinate system. The histories of ~ , ~y and Cz are calculated by different means depending on the mode in which the package is used, and also depending on what arc of the trajectory is being integrated. In the performance mode, the ~,., ~. and ~z histories come from the pitch and yaw steering program. Given a set of pitch and yaw attitudes at an
Launch trajectory optimization
847
INPUTS
VEHICLE Propulsion Data, Aerodynamic Data
CONSTRAINTS
Stage Weights, Atmosphere, Winds
j
Structural Loads {DO() Land Mass Boundary Heat Shield Ejection Turing Rates Coast Times Stage Separation Orbit Inclination Orbit Size Perigee Location
Present Estimate of Steering and Sequence Parameters
APOGEE STAGE "Variable or Constant Orientation "Targeting for Orbit Inclination, Horizontal Injection or for given Perigee Location "Injection Point "Performance Index
STAGES 2 & 3
STAGES 0 & I Vertical Rise Pitch until ~ m a x . Transition to Gravity Turn Gravity Turn Arc Pitch & Yaw Steering Land Mass Constraint Coast after Burnout
_1 ~l
Pitch & Yaw Steering Heat Shield Ejection Land Mass Constraint Coast Times
Il
To Optimization Routine
=
Fig. I. Trajectory profile in PYOPT optimization mode. instant, the direction cosines of the vehicle axis are uniquely determined. The initial pitch and yaw attitude angles are not a very convenient set of angles for the purpose of optimization of trajectory; particularly when the vehicle is outside atmosphere. In PYOPT, a localized set of angles 0# and O,~zis employed in optimization mode in order to facilitate giving good initial guess estimates for a starting trajectory and easily control the trajectory during active-land-mass-constraint arcs. 0# is a local pitch angle measured from the local vertical to the vehicle axis. OAzis the azimuthal angle of the vehicle axis, measured from local north. When 0# and OAz are known, the direction cosines of
orientation vector ~ can be calculated by setting up a suitable local orthogonal triad. The time history of the local pitch angle 0p~will be parametrically optimized during exoatmospheric phase. OAz is not directly optimized, but optimized through an incremental angle above the azimuth calculated locally to obtain required orbital inclination, subject to land-mass-clearance constraint and maximum turning rate constraint. When an option for planar trajectory optimization is activated, the thrust direction will always be chosen parallel to the original pitch plane. In such a case, 0# history is parametrically optimized, with OAzchosen locally to satisfy planar thrusting condition.
INPUTS
I
I eh~o~elnp~ts ]
Perturbation Levels In Mission Critical Parameters
J [ Stet!ring&
Sequence
Inputs
I
t
Performance Mode
VEHICLE DYNAMICS Three Dimensional Motion with Ideal Control for Steering Program
Steering Parameters
_
[
I
l
t I I
I I
1
Optimization Mode
CONSTRAINED OPTIMIZATION ALGORITHMS FOR Maximum Payload, Velocity etc.
Maximum
I
Performance & Sensitivity
!-
l
]
~ Trial
[
MOde
I
i [Optimiz~d:~teering| and | Perfor...... :e | ~ - ~ Sequence
Fig. 2. Summary of major options in PYOPT program. A.A. 15/11--D
]
hRun for'
I
I OUTPUTS
Vehicle Performance History
F
Perturbation Mode
][.... I] New
I
~ission & Vehicle Cons t ra ints
848
V. ADIMURTHY
2.4. Software design
The main features of the trajectory profile in the optimization mode of PYOPT are described in Fig. 1. The major options available in the software are summarised in Fig. 2. Complete details of trajectory model, auxiliary calculations and optimization procedures are given in the software document[7].
Consider the augmented Lagrangian: L(x, 2, w) = f ( x ) + 2rh(x) + w h r ( x ) h ( x )
(6)
Parameters: a>~0,
b~>0,
~0>0,
a+b>l,
01e(0,1),
{~,}~ ~-~, > 0,
0,.e(0,1),
{~,,}--,0, ~ ~i< :~ i=0
3. O P T I M I Z A T I O N M E T H O D O L O G Y
The main drawback of the conventional methods of penalty functions, Lagrange multiplier or augmented Lagrangian is that they require a sequence of unconstrained minimization problems to be solved. This would require a large amount of computer time. Also questions like up to what accuracy each unconstrained minimization has to be carried out would arise. One approach to overcome these problems related to 'iteration within iteration' is to device continuous penalty functions whose unconstrained minimum lies at the solution of the original constrained problem. Such functions are called exact penalty functions. The main computational problem arising in these methods is the presence of first-order derivatives in the augmented Lagrangian. Any gradient-dependent minimization method would thus require second-order derivatives of the Lagrangian function. This approach is not promising for functions, calculation of whose derivatives is computationally costly because of numerical differentiation. In addition, inaccuracies in the evaluation of second derivatives by numerical methods will add to further limitation of these methods. From a different direction, Han[5] and independently Tapia[4] have presented an algorithm which takes only one unconstrained minimization of the augmented Lagrangian function. Tapia calls this method a diagonalized multiplier method. Motivation for this terminology is that if in the matrix A = (i, i), i,j = 0 , 1, 2 . . . . we let ( i , j ) correspond to thejth iteration of a method on the ith unconstrained problem (equivalently the ith multiplier update), then the diagonalized multiplier method corresponds to traversing the diagonal of A, while the usual multiplier method corresponds to going along each row of A in succession. The success of the diagonalised approach is mainly due to the special role played by the update formula for Lagrange multiplier proposed by Han-Tapia. Bertocchi et al.[6] make extensive numerical experiments to evaluate the performance of the diagonalized multiplier method and find that the strong superiority of this method is borne out in terms of accuracy, reliability and efficiency. The optimization scheme used in PYOPT is based on the diagonalised multiplier method and can be described in the following steps: Problem: minimize f (x ), f (x ) : R"--* R 1, subject to h(x) = O, h ( x ) : R " ~ R m.
Initiation set j = 0, K = 1, start with x 0 for minimum point, 20 for Lagrange multipliers, w0 for quadratic multiplier, H 0 for inverse Hessian of L
Step 1: Let x~=xo-~toVL(xo,2o, Wo) where ct0 is chosen on the basis of linear search Step 2: Scale H 0 by multiplying it by the scaling parameter of Shanno and Phua[8], obtain H~ from H0 by BFGS update formula Step 3: Calculate the multiplier increment of Han[5] and Tapia[4] given by d2 = (VrhHVh) ~(h - V r h H V L ) and l e t 2 x = 2 K
/~ = Step 4:
(7)
~+fld2where v r ( f + ~,rH)Vh d2 (Vh d2 )rVh d2
If 11fl d2 II2 ~< ~K go to Step 5
(8)
(9)
otherwise, set wj+l = a w j + b , j = j + 1 and go to Step 5 Step 5:
xx+ l = x x - ~rHxV~L(x x, 2 x, v~(x))
(10)
Here ctK may be taken equal to 1 or determined by a simple linear search, the later is preferred Step 6: Check for convergence (llh Iland IIVL tl are sufficiently small) If converged STOP. If not, generate H x +~ by BFGS formula 1 +0 2
let ~ x + l = O 2 ~ x + ~ - 2 set K = K + I
7K~l
(11)
and go to Step 3.
In the above scheme we use the Han-Tapia multiplier update formula in conjuction with BFGS inverse-Hessian update to minimize the HestenesPowell augmented Lagrangian. An initial scaling of the inverse Hessian matrix, on the basis of Shanno and Phua and a modification of quadratic multiplier up-datement given in Polak and Tits[9] are embedded. It may be noted that in Han's original algorthm, the initial value of 2 = 2o is of no consequence. In Han's method initially the Lagrange multipliers are updated before any updatement of the H matrix, which is independent of 20. If H0 is not proper, this leads to a bad initial approximation of 21, and in the subsequent application of the algorithm the convergence can be badly affected. To avoid this, we
Launch trajectory optimization modify the algorithm by minimizing the K u h n Tucker vector norm along the Han-Tapia update direction.
849
Primal problem P: minimize f(x),
subject to
g(x) <~O, h(x) = 0
Dual problem D: 4. OPTIMALITY OF SOLUTIONS
maximize 0(/~, 2),
As described above, the program PYOPT can be used to optimize the payload of a satellite launch vehicle under various specified constraints. A typical study of injection into 900 km circular Sun-synchronous orbit is carried out with three dimensional steering requirement due to land-mass constraints. Here, the function f, to be optimized is the payload mass and the terminal requirements are that the apogee altitude is 900 km (condition ht) and that the perigee altitude is also 900 km (condition h2). Several other vehicle and mission constraints are included in system dynamic model itself. A sixparameter model for the optimization problem is chosen governing the pitch and yaw steering program. The optimality of the solutions obtained by PYOPT is studied using three approaches (a) check a general form of Kuhn Tucker sufficiency conditions for optimality (b) using saddle point optimality theorem and (c) using augmentability theorem of nonlinear programming. All the theorems used in this study are described in our detailed document[10] and a summary is presented in this paper.
4.1. Kuhn- Tucker sufficiency conditions A general form of Kuhn-Tucker sufficiency conditions for optimality is proved by Bazaraa and Shetty[ll], where sufficiency is established though relaxed form of convexity conditions on f and h. First it has been found that strict convexity at the solution point of PYOPT is not existing. Then the Bazaraa and Shetty theorem is applied for the relaxed quasiconvexity conditions. Direct checking of quasiconvexity is difficult, and a characterisation by boardered determinants is utilised. The outcome of these tests also happens to be indeterminate in establishing the optimality.
subject to
/~/> 0,
where 0=infq~(x,/~,2);
q~=f+#g+2h
Let X be a nonempty set in E,. Suppose there exist x* e X and' #*, 2" with p* > 0, such that
q~(x*,l~,2)<~ q6(X*,l~*,2*)<~ ~(X,l~*,2*)
(12)
for all x e X and all #,2 with /~ t>0. Then x*, #* and 2" solve the primal problem P and the dual problem D. The left-hand inequality of eqn. (12) simply implies that the constraints of problem P are satisfied. In PYOPT we do get h~ = 0 (the required perigee and apogee altitudes) to reasonable accuracy of one meter. The right-hand inequality of eqn (12) implies that has a minimum with respect to x when the right/~* and 2* are used. There is no active inequality in PYOPT model. The 2~ and 22 corresponding to h I and h 2 are seen to be 21 = -0.43215 and 22 = -0.16916. If 21H(hl)+ 22H(h2) is positive definite at the solution point, then from continuity arguments it can be shown that the right-hand inequality in eqn (12) holds. Eigenvalues of 2~H(hl)+ 22H(h2) are: 0.024331, 0.167591, 0.307968, 0.761827, 4.02227, 11.14908. This shows that ~b(x,/~*, 2*) is positive definite at x = x*. Hence from saddle-point sufficient optimality theorem, we conclude that the solution obtained in PYOPT is local optimum.
4.3. Optimality using augmentability theorem Consider the problem PE: minimize f(x),
subject to h(x) = O.
The augmentability theorem (see [12]) shows that if
4.2. Saddle point optimality theorem
h(x °) = 0 and x ° is an unconstrained minimum of L(x, 2, w) for some 2 then x ° solves PE. By stopping criteria at PYOPT solution point h(x)= 0. Also as
In this approach, the saddle point optimality theorem based on duality theory of nonlinear programming is utilized to check if the solution given by PYOPT is optimal.
we established that Kuhn-Tucker necessary conditions are satisfied, it follows that VL(x, 2, w)= 0 for any w. Then to show L(x, 2, w) has minimum at solution point, it is enough to show that Hessian of
Table 1. Initial guesses in PYOPT convergence study Route No.
1
1 2 3 4 5 6 7 8 9
30.0 80.0 50.0 60.0 70.0 60.0 30.0 80.0 40.0
Initial guess parameters 2 3 4 40.0 80.0 60.0 60.0 70.0 70.0 75.0 50.0 40.0
50.0 80.0 70.0 60.0 30.0 80.0 80.0 50.0 80.0
0.0 --25.0 10.0 5.0 30.0 10.0 -30.0 --50.0 --25.0
Initial orbit (km)
5
Initial guess payload (kg)
Apogee ht
Perigee ht
0.0 -25.0 10.0 5.0 --30.0 10.0 --30.0 +50.0 --25.0
500.00 2000.0 500.0 600.0 1300.0 1500.0 100.0 700.0 500.0
2436.286 316.500 4712.657 2395.667 1092.772 573.788 2689.358 1161.214 1469.864
--51.362 --3334.965 1276.773 1435.236 --3387.172 -630.983 1228.014 - 3649.794 --2810.920
850
V. ADIMURTHY Table 2. Converged values from different initial guesses Route No.
1
1 2 3 4 5 6 7 8 9
51.91 l 52.011 51.990 51.989 51.970 50.991 51.989 51.998 53.083
Converged parameter 2 8 4 70.269 70.211 70.222 70.223 70.223 70.752 70.228 70.206 69.543
69.982 69.975 69.679 69.980 69.975 70.088 69.978 69.989 69.938
11.695 I 1.705 11.669 11.690 11.700 11.663 I 1.680 11.659 11.622
L(x, 2, w) is positive definite. We show this for w = 0.0001. Eigenvalues of the Hessian are: 0. ! 68009, 0.250361, 0.621234, 2.112908, 4.034475, 11.340258. Hence the Hessian is positive definite. We conclude t h a t the solution o b t a i n e d in P Y O P T is a local optimum. 5. ROBUSTNESS OF PYOPT SOFTWARE The robustness of the solutions generated by P Y O P T software is established t h r o u g h r u n n i n g the p r o g r a m with a n u m b e r of widely varying initial guess points for optimization. A s u m m a r y o f the experiments is given in Tables 1 a n d 2. Table 1 summarizes the initial guesses made. The payload guesses varied from 100kg to 2000 kg. Yaw angle varied from + 50 deg to - 50 deg. The initial orbits so o b t a i n e d are extremely wild. Table 2 summarizes the situation after convergence has been achieved in P Y O P T t h r o u g h all the different routes. All the solution points are the same. This shows t h a t the optimization procedure in P Y O P T is robust. This also suggests t h a t the solution o b t a i n e d is possibly n o t only locally optimal, but also globally optimal. 6. CONCLUDING REMARKS
A launch vehicle trajectory optimization software is developed using diagonalized multiplier approach. A new modification, in which the K u h n - T u c k e r vector n o r m is minimized along the multiplier u p d a t e direction, avoids violent initial multiplier corrections. The optimality of the solutions is established, a n d robustness of the package demonstrated. This package is presently employed as a s t a n d a r d tool for mission design studies at the I n d i a n Space Research Organisation.
Acknowledgements--The author would like to record the contributions of K. V. Joy, M. Karunakaran and Jessy Antony in programming the PYOPT software package
_ 5 13.075 13.074 13.078 13.075 13.066 13.062 13.086 13.117 13.208
Converged payload (kg) 1036.94 1036.95 1036.94 1036.94 1036.75 1037.01 1036.95 1037.01 1036.82
Apogeeht (km) 900.001 900.013 900.008 900.005 900.142 900.000 900.032 899.951 899.995
Perigeeht (km) 899.999 899.891 900.001 900.000 900.106 899.998 899.966 899.702 899.993
and running several case studies. Technical discussions with T. S. Jayaraman during the review of this software are also gratefully acknowledged.
REFERENCES
1. R. G. Brusch, Trajectory optimization for the Atlas/ Centaur launch vehicle. J. Spacecraft Rockets 14, 550-555 (1977). 2. G. L. Brauer, D. E. Cornick and R. Stevenson. Capabilities and applications of the Program to Optimize Simulated Trajectories (POST). NASA CR-2770 (1977). 3. K. H. Well and S. R. Tandon, Rocket ascent trajectory optimization via recursive quadratic programming. J. Astron. Sci. 30, 101--116 (1982). 4. R. A. Tapia. Diagonalized multiplier methods and quasi-Newton methods for constrained optimization. J. Optimiz. Theory Appl. 22, 135-194 (1977). 5. S. P. Han, Dual variable metric algorithms for constrained optimization. SIAM J. Control Optimization 15, 546-565 (1977). 6. M. Bertocchi, E. Cavalli and E. Spedicato, Computational performance of diagonalized multiplier quasi-Newton methods for nonlinear optimization with equality constraints. In Numerical Optimization o/' Dynamical Systems pp. 247 268, (Edited by L. C. W. Dixon and C. P. Szego), North-Holland, Amsterdam, (1980). 7. V. Adimurthy, K. V. Joy and M. Karunakaran, PYOPT 81, a trajectory optimization software package. VSSCTR-13-231-81, Vikram Sarabhai Space Centre, Trivandrum (1981). 8. D. F. Shanno and K. H. Phua, Matrix conditioning and nonlinear optimization. Math. Programming 14, 149 160 (1978). 9. E. Polak and A. L. Tits, A globally convergent implementable multiplier method with automatic penalty limitation. Appl. Math. Optimiz. 6, 335-360 (1980). 10. V. Adimurthy, Optimality in PYOPT. SSTC-ARD-FEB (1)/82, Vikram Sarabhai Space Centre Trivandrum (1982). 11. M. S. Bazaraa and C. M. Shetty, The Fritz John and the Kuhn Tucker optimality conditions, Nonlinear Programming, Theory and Algorthms. pp. 147 150, John Wiley, New York. (1979). 12. M. Avriel, Penalty function methods. Nonlinear Programming, Analysis and Methods'. p. 401. Prentice-Hall, N.J. (1976).
0094-5767/87 $3.00 + 0.00 © 1987 Pergamon Journals Ltd
LAUNCH VEHICLE TRAJECTORY OPTIMIZATION# V. ADIMURTHY Vikram Sarabhai Space Centre, Trivandrum~595022, India (Received 18 December 1986)
A~traet--Three dimensional trajectory optimization problem is solved through a diagonalized multiplier method in which the multiplier update formula of Tapia and Han is employed over a BFGS Hessian update. The diagonalized multiplier method (DMM) is a single sequence minimization technique for nonlinear programming and avoids the sequence of minimization problems inherent in the classical augmented Lagrangian methods. To avoid violent initial multiplier corrections in DMM, a new modification is proposed in which we minimize the Kuhn Tucker vector norm along the Han-Tapia multiplier update direction. A software package PYOPT is developed and validated at the Indian Space Research Organisation (ISRO) employing the concepts of diagonalised multiplier method for constrained pitch and yaw steering optimization of launch vehicle trajectories. The system dynamic model for trajectory optimization allows for effective and routine operational studies. The optimality of the solutions of PYOPT is established through sufficiency conditions of optimality in nonlinear programming via saddle point optimality theorem and augmentability theorem. The robustness of the package is demonstrated. The present package, which is faster than programmes based on techniques like projected gradient methods and recursive quadratic programming, has become a standard tool at ISRO for mission design studies.
1. INTRODUCTION
Trajectory Optimization with detailed launch vehicle and mission constraints is a complex nonlinear problem, and efforts to solve this effectively using the state of the art continue to be made. The reason for considerable distance between theory and practice and for the sparseness of reported papers in this area depends on the difficulty in establishing good convergence in numerical iterative procedures. Brusch[1] addresses to the trajectory optimization for Atlas/Centaur launch vehicle employing a classical augmented Lagrangian method, in which the solution is sought by alternately minimizing the augmented function with respect to independent variables and adjusting the Lagrange multipliers to satisfy constraints. POST, a generalized Program to Optimize Stimulated Trajectories[2], employs an accelerated projected gradient method. These approaches involve considerable computer time for convergence. Well and Tandon[3] effectively use the recursive quadratic programming approach to three dimensional trajectory optimization. They find that this method seems to be a robust and relatively fast algorithm for solving parameterized optimal control problems. But they conclude that further modifications of the algorithm have to be done in order to reduce CPU-time such that actual mission design studies can be carried out effectively. However, none of these papers address to establishing, theoretically or otherwise, the actual optimality of the solutions. tPaper IAF-86-235 presented at the 37th Congress of the International Astronautical Federation, Innsbruck, Austria, 4-11 October 1986
In the present paper, we solve the threedimensional trajectory optimization problem through a diagonalized multiplier method in which multiplier update formula of Tapia[4] and Han[5] is employed over a BFGS Hessian update. Like recursive quadratic programming (RQP), the diagonalized multiplier method (DMM) is a single sequence minimization technique for nonlinear programming and avoids the sequence of minimization problems inherent in classical augmented Lagrangian methods. This leads to the substantially faster performance of D M M and RQP compared to classical methods, as found by Bertocchi et al.[6], who make a comparative study of the DMM, RQP and classical methods. It is also found that the performance of the diagonalized multiplier methods is better in several case studies. A software package PYOPT is developed and validated at the Indian Space Research Organization (ISRO) employing concepts of diagonalized multiplier method for constrained pitch and yaw trajectory optimization. To avoid violent initial multiplier corrections, a new algorthm is proposed in which we minimize the Kuhn-Tucker vector norm along the Han-Tapia multiplier update directions. The optimality of the solutions of PYOPT is established through sufficiency conditions of optimality in nonlinear programming via saddle point optimality theorem and augmentability theorem. The system dynamics aspects of the trajectory optimization are modelled to allow for effective and routine operational studies. This paper describes the system modelling, algorithmic aspects, and studies made to establish optimality and robustness. The present method is faster than existing software packages like POST and those based on recursive quadratic pro845
846
V. AD1MURTHY
gramming, and has become a standard tool at-ISRO for mission design studies. 2. MATHEMATICAL FORMULATION
2.1. General statement of the problem Let vector s denote the state (say, position and velocity) of the dynamical system under consideration (launch vehicle in motion). Let the dynamics of the system depend on a vector u of control variables (say, pitch, yaw and other design parameters), which can be specified at our will (assuming that specified pitch and yaw histories can be realised through the control system). The dynamics of the system can now be represented by ds ~=F(s,u,t),
to<~t<~ O.
(1)
We are interested in maximizing a performance index (payload, injection velocity, range) which depends, among other things, on s(tj). Since this depends on the choice of u through the dynamical equations of the system, the performance index also depends on u. The problem is to select the control vector u, such that the performance index is maximized, and at the same time various specified objectives and constraints are satisfied. The objectives can be the final orbit size and inclination and the like. The constraints can be on the vehicle loading, flight safety, etc. To solve this problem one should be equipped with (a) a capability to simulate the system dynamics, (b) a procedure to ensure that the specified constraints are satisfied and (c) a methodology to vary the control vector so as to maximize the performance index. In the scheme adopted in the present program PYOPT, an attempt has been made to model some of the constraints in the system dynamics itself, while the rest are left for optimization algorithm. The system dynamics is considered as a black-box that satisfies these constraints and gives out the performance index, once a set of control variables is specified. This scheme is equivalent to making a transformation on the dynamics in order to make a highly constrained problem a much less constrained one. The system-black-box concept leads to a direct functional relationship between the performance and the control vector. Further, it is assumed that the control variables have been parametrized, so that the problem is to find a finite number of control parameters instead of a continuous control function. With this parametrization, the performance index becomes a function of a finite number of control parameters. The optimization of this function is done using methods of constrained nonlinear optimization. The approach described above differs from that of other workers like Brusch[l] in one essential feature. Our system dynamics for optimization is not simply an integration of a system of basic differential
equations. Various in-flight constraints have been embedded in the system dynamics. This feature of the present study has necessitated in a careful modelling of the transformed system dyanamics and the flight profile.
2.2. Dynamical equations The basic coordinate system, in which all the computations in PYOPT are carried out, is a Geocentric Inertial System (x,y,z). In this frame, the x-y plane is the equatorial plane. The z-axis is pointed towa~'ds the north pole. x-z plane contains the launch meridian frozen at the instant of initial stage ignition (at t = 0). The equations of motion in this frame can be written as: __d2_--T ~x +
dt
~.~cos ~
m
m
ga
x dS'
_
dt
3 a4
T ~y + QS ~y cos ~
Vay
m
(
~aa}
m-
( Xce--rc~ QS¢,.C ~ × Cs~ 1+ Xc~ / - m " - - g 0Y- [ ~ r d~
-
=
dt
-
T 4: +
a4 _ 5zZ']] + 53Jz~-g(1 rZ ./3
(3)
~, cos
m
Va J
z
[-a 2
dx
3 dy
dz
d~=.~, h5=?, ~5=_~
(5)
The forces modelled in the above differential equations are (a) thrust (corrected to local pressure of the atmosphere), (b) aerodynamic normal force, (c) ideal control force to balance aerodynamic moments, (d) aerodynamic drag and (e) gravity force for an oblate earth (up to second harmonic only).
2.3. Vehicle attitude representation In PYOPT, the vehicle attitude information is basically recognised through the direction cosines ~x, ~y and ~: of the attitude vector ~ with respect to the basic coordinate system. The histories of ~ , ~y and Cz are calculated by different means depending on the mode in which the package is used, and also depending on what arc of the trajectory is being integrated. In the performance mode, the ~,., ~. and ~z histories come from the pitch and yaw steering program. Given a set of pitch and yaw attitudes at an
Launch trajectory optimization
847
INPUTS
VEHICLE Propulsion Data, Aerodynamic Data
CONSTRAINTS
Stage Weights, Atmosphere, Winds
j
Structural Loads {DO() Land Mass Boundary Heat Shield Ejection Turing Rates Coast Times Stage Separation Orbit Inclination Orbit Size Perigee Location
Present Estimate of Steering and Sequence Parameters
APOGEE STAGE "Variable or Constant Orientation "Targeting for Orbit Inclination, Horizontal Injection or for given Perigee Location "Injection Point "Performance Index
STAGES 2 & 3
STAGES 0 & I Vertical Rise Pitch until ~ m a x . Transition to Gravity Turn Gravity Turn Arc Pitch & Yaw Steering Land Mass Constraint Coast after Burnout
_1 ~l
Pitch & Yaw Steering Heat Shield Ejection Land Mass Constraint Coast Times
Il
To Optimization Routine
=
Fig. I. Trajectory profile in PYOPT optimization mode. instant, the direction cosines of the vehicle axis are uniquely determined. The initial pitch and yaw attitude angles are not a very convenient set of angles for the purpose of optimization of trajectory; particularly when the vehicle is outside atmosphere. In PYOPT, a localized set of angles 0# and O,~zis employed in optimization mode in order to facilitate giving good initial guess estimates for a starting trajectory and easily control the trajectory during active-land-mass-constraint arcs. 0# is a local pitch angle measured from the local vertical to the vehicle axis. OAzis the azimuthal angle of the vehicle axis, measured from local north. When 0# and OAz are known, the direction cosines of
orientation vector ~ can be calculated by setting up a suitable local orthogonal triad. The time history of the local pitch angle 0p~will be parametrically optimized during exoatmospheric phase. OAz is not directly optimized, but optimized through an incremental angle above the azimuth calculated locally to obtain required orbital inclination, subject to land-mass-clearance constraint and maximum turning rate constraint. When an option for planar trajectory optimization is activated, the thrust direction will always be chosen parallel to the original pitch plane. In such a case, 0# history is parametrically optimized, with OAzchosen locally to satisfy planar thrusting condition.
INPUTS
I
I eh~o~elnp~ts ]
Perturbation Levels In Mission Critical Parameters
J [ Stet!ring&
Sequence
Inputs
I
t
Performance Mode
VEHICLE DYNAMICS Three Dimensional Motion with Ideal Control for Steering Program
Steering Parameters
_
[
I
l
t I I
I I
1
Optimization Mode
CONSTRAINED OPTIMIZATION ALGORITHMS FOR Maximum Payload, Velocity etc.
Maximum
I
Performance & Sensitivity
!-
l
]
~ Trial
[
MOde
I
i [Optimiz~d:~teering| and | Perfor...... :e | ~ - ~ Sequence
Fig. 2. Summary of major options in PYOPT program. A.A. 15/11--D
]
hRun for'
I
I OUTPUTS
Vehicle Performance History
F
Perturbation Mode
][.... I] New
I
~ission & Vehicle Cons t ra ints
848
V. ADIMURTHY
2.4. Software design
The main features of the trajectory profile in the optimization mode of PYOPT are described in Fig. 1. The major options available in the software are summarised in Fig. 2. Complete details of trajectory model, auxiliary calculations and optimization procedures are given in the software document[7].
Consider the augmented Lagrangian: L(x, 2, w) = f ( x ) + 2rh(x) + w h r ( x ) h ( x )
(6)
Parameters: a>~0,
b~>0,
~0>0,
a+b>l,
01e(0,1),
{~,}~ ~-~, > 0,
0,.e(0,1),
{~,,}--,0, ~ ~i< :~ i=0
3. O P T I M I Z A T I O N M E T H O D O L O G Y
The main drawback of the conventional methods of penalty functions, Lagrange multiplier or augmented Lagrangian is that they require a sequence of unconstrained minimization problems to be solved. This would require a large amount of computer time. Also questions like up to what accuracy each unconstrained minimization has to be carried out would arise. One approach to overcome these problems related to 'iteration within iteration' is to device continuous penalty functions whose unconstrained minimum lies at the solution of the original constrained problem. Such functions are called exact penalty functions. The main computational problem arising in these methods is the presence of first-order derivatives in the augmented Lagrangian. Any gradient-dependent minimization method would thus require second-order derivatives of the Lagrangian function. This approach is not promising for functions, calculation of whose derivatives is computationally costly because of numerical differentiation. In addition, inaccuracies in the evaluation of second derivatives by numerical methods will add to further limitation of these methods. From a different direction, Han[5] and independently Tapia[4] have presented an algorithm which takes only one unconstrained minimization of the augmented Lagrangian function. Tapia calls this method a diagonalized multiplier method. Motivation for this terminology is that if in the matrix A = (i, i), i,j = 0 , 1, 2 . . . . we let ( i , j ) correspond to thejth iteration of a method on the ith unconstrained problem (equivalently the ith multiplier update), then the diagonalized multiplier method corresponds to traversing the diagonal of A, while the usual multiplier method corresponds to going along each row of A in succession. The success of the diagonalised approach is mainly due to the special role played by the update formula for Lagrange multiplier proposed by Han-Tapia. Bertocchi et al.[6] make extensive numerical experiments to evaluate the performance of the diagonalized multiplier method and find that the strong superiority of this method is borne out in terms of accuracy, reliability and efficiency. The optimization scheme used in PYOPT is based on the diagonalised multiplier method and can be described in the following steps: Problem: minimize f (x ), f (x ) : R"--* R 1, subject to h(x) = O, h ( x ) : R " ~ R m.
Initiation set j = 0, K = 1, start with x 0 for minimum point, 20 for Lagrange multipliers, w0 for quadratic multiplier, H 0 for inverse Hessian of L
Step 1: Let x~=xo-~toVL(xo,2o, Wo) where ct0 is chosen on the basis of linear search Step 2: Scale H 0 by multiplying it by the scaling parameter of Shanno and Phua[8], obtain H~ from H0 by BFGS update formula Step 3: Calculate the multiplier increment of Han[5] and Tapia[4] given by d2 = (VrhHVh) ~(h - V r h H V L ) and l e t 2 x = 2 K
/~ = Step 4:
(7)
~+fld2where v r ( f + ~,rH)Vh d2 (Vh d2 )rVh d2
If 11fl d2 II2 ~< ~K go to Step 5
(8)
(9)
otherwise, set wj+l = a w j + b , j = j + 1 and go to Step 5 Step 5:
xx+ l = x x - ~rHxV~L(x x, 2 x, v~(x))
(10)
Here ctK may be taken equal to 1 or determined by a simple linear search, the later is preferred Step 6: Check for convergence (llh Iland IIVL tl are sufficiently small) If converged STOP. If not, generate H x +~ by BFGS formula 1 +0 2
let ~ x + l = O 2 ~ x + ~ - 2 set K = K + I
7K~l
(11)
and go to Step 3.
In the above scheme we use the Han-Tapia multiplier update formula in conjuction with BFGS inverse-Hessian update to minimize the HestenesPowell augmented Lagrangian. An initial scaling of the inverse Hessian matrix, on the basis of Shanno and Phua and a modification of quadratic multiplier up-datement given in Polak and Tits[9] are embedded. It may be noted that in Han's original algorthm, the initial value of 2 = 2o is of no consequence. In Han's method initially the Lagrange multipliers are updated before any updatement of the H matrix, which is independent of 20. If H0 is not proper, this leads to a bad initial approximation of 21, and in the subsequent application of the algorithm the convergence can be badly affected. To avoid this, we
Launch trajectory optimization modify the algorithm by minimizing the K u h n Tucker vector norm along the Han-Tapia update direction.
849
Primal problem P: minimize f(x),
subject to
g(x) <~O, h(x) = 0
Dual problem D: 4. OPTIMALITY OF SOLUTIONS
maximize 0(/~, 2),
As described above, the program PYOPT can be used to optimize the payload of a satellite launch vehicle under various specified constraints. A typical study of injection into 900 km circular Sun-synchronous orbit is carried out with three dimensional steering requirement due to land-mass constraints. Here, the function f, to be optimized is the payload mass and the terminal requirements are that the apogee altitude is 900 km (condition ht) and that the perigee altitude is also 900 km (condition h2). Several other vehicle and mission constraints are included in system dynamic model itself. A sixparameter model for the optimization problem is chosen governing the pitch and yaw steering program. The optimality of the solutions obtained by PYOPT is studied using three approaches (a) check a general form of Kuhn Tucker sufficiency conditions for optimality (b) using saddle point optimality theorem and (c) using augmentability theorem of nonlinear programming. All the theorems used in this study are described in our detailed document[10] and a summary is presented in this paper.
4.1. Kuhn- Tucker sufficiency conditions A general form of Kuhn-Tucker sufficiency conditions for optimality is proved by Bazaraa and Shetty[ll], where sufficiency is established though relaxed form of convexity conditions on f and h. First it has been found that strict convexity at the solution point of PYOPT is not existing. Then the Bazaraa and Shetty theorem is applied for the relaxed quasiconvexity conditions. Direct checking of quasiconvexity is difficult, and a characterisation by boardered determinants is utilised. The outcome of these tests also happens to be indeterminate in establishing the optimality.
subject to
/~/> 0,
where 0=infq~(x,/~,2);
q~=f+#g+2h
Let X be a nonempty set in E,. Suppose there exist x* e X and' #*, 2" with p* > 0, such that
q~(x*,l~,2)<~ q6(X*,l~*,2*)<~ ~(X,l~*,2*)
(12)
for all x e X and all #,2 with /~ t>0. Then x*, #* and 2" solve the primal problem P and the dual problem D. The left-hand inequality of eqn. (12) simply implies that the constraints of problem P are satisfied. In PYOPT we do get h~ = 0 (the required perigee and apogee altitudes) to reasonable accuracy of one meter. The right-hand inequality of eqn (12) implies that has a minimum with respect to x when the right/~* and 2* are used. There is no active inequality in PYOPT model. The 2~ and 22 corresponding to h I and h 2 are seen to be 21 = -0.43215 and 22 = -0.16916. If 21H(hl)+ 22H(h2) is positive definite at the solution point, then from continuity arguments it can be shown that the right-hand inequality in eqn (12) holds. Eigenvalues of 2~H(hl)+ 22H(h2) are: 0.024331, 0.167591, 0.307968, 0.761827, 4.02227, 11.14908. This shows that ~b(x,/~*, 2*) is positive definite at x = x*. Hence from saddle-point sufficient optimality theorem, we conclude that the solution obtained in PYOPT is local optimum.
4.3. Optimality using augmentability theorem Consider the problem PE: minimize f(x),
subject to h(x) = O.
The augmentability theorem (see [12]) shows that if
4.2. Saddle point optimality theorem
h(x °) = 0 and x ° is an unconstrained minimum of L(x, 2, w) for some 2 then x ° solves PE. By stopping criteria at PYOPT solution point h(x)= 0. Also as
In this approach, the saddle point optimality theorem based on duality theory of nonlinear programming is utilized to check if the solution given by PYOPT is optimal.
we established that Kuhn-Tucker necessary conditions are satisfied, it follows that VL(x, 2, w)= 0 for any w. Then to show L(x, 2, w) has minimum at solution point, it is enough to show that Hessian of
Table 1. Initial guesses in PYOPT convergence study Route No.
1
1 2 3 4 5 6 7 8 9
30.0 80.0 50.0 60.0 70.0 60.0 30.0 80.0 40.0
Initial guess parameters 2 3 4 40.0 80.0 60.0 60.0 70.0 70.0 75.0 50.0 40.0
50.0 80.0 70.0 60.0 30.0 80.0 80.0 50.0 80.0
0.0 --25.0 10.0 5.0 30.0 10.0 -30.0 --50.0 --25.0
Initial orbit (km)
5
Initial guess payload (kg)
Apogee ht
Perigee ht
0.0 -25.0 10.0 5.0 --30.0 10.0 --30.0 +50.0 --25.0
500.00 2000.0 500.0 600.0 1300.0 1500.0 100.0 700.0 500.0
2436.286 316.500 4712.657 2395.667 1092.772 573.788 2689.358 1161.214 1469.864
--51.362 --3334.965 1276.773 1435.236 --3387.172 -630.983 1228.014 - 3649.794 --2810.920
850
V. ADIMURTHY Table 2. Converged values from different initial guesses Route No.
1
1 2 3 4 5 6 7 8 9
51.91 l 52.011 51.990 51.989 51.970 50.991 51.989 51.998 53.083
Converged parameter 2 8 4 70.269 70.211 70.222 70.223 70.223 70.752 70.228 70.206 69.543
69.982 69.975 69.679 69.980 69.975 70.088 69.978 69.989 69.938
11.695 I 1.705 11.669 11.690 11.700 11.663 I 1.680 11.659 11.622
L(x, 2, w) is positive definite. We show this for w = 0.0001. Eigenvalues of the Hessian are: 0. ! 68009, 0.250361, 0.621234, 2.112908, 4.034475, 11.340258. Hence the Hessian is positive definite. We conclude t h a t the solution o b t a i n e d in P Y O P T is a local optimum. 5. ROBUSTNESS OF PYOPT SOFTWARE The robustness of the solutions generated by P Y O P T software is established t h r o u g h r u n n i n g the p r o g r a m with a n u m b e r of widely varying initial guess points for optimization. A s u m m a r y o f the experiments is given in Tables 1 a n d 2. Table 1 summarizes the initial guesses made. The payload guesses varied from 100kg to 2000 kg. Yaw angle varied from + 50 deg to - 50 deg. The initial orbits so o b t a i n e d are extremely wild. Table 2 summarizes the situation after convergence has been achieved in P Y O P T t h r o u g h all the different routes. All the solution points are the same. This shows t h a t the optimization procedure in P Y O P T is robust. This also suggests t h a t the solution o b t a i n e d is possibly n o t only locally optimal, but also globally optimal. 6. CONCLUDING REMARKS
A launch vehicle trajectory optimization software is developed using diagonalized multiplier approach. A new modification, in which the K u h n - T u c k e r vector n o r m is minimized along the multiplier u p d a t e direction, avoids violent initial multiplier corrections. The optimality of the solutions is established, a n d robustness of the package demonstrated. This package is presently employed as a s t a n d a r d tool for mission design studies at the I n d i a n Space Research Organisation.
Acknowledgements--The author would like to record the contributions of K. V. Joy, M. Karunakaran and Jessy Antony in programming the PYOPT software package
_ 5 13.075 13.074 13.078 13.075 13.066 13.062 13.086 13.117 13.208
Converged payload (kg) 1036.94 1036.95 1036.94 1036.94 1036.75 1037.01 1036.95 1037.01 1036.82
Apogeeht (km) 900.001 900.013 900.008 900.005 900.142 900.000 900.032 899.951 899.995
Perigeeht (km) 899.999 899.891 900.001 900.000 900.106 899.998 899.966 899.702 899.993
and running several case studies. Technical discussions with T. S. Jayaraman during the review of this software are also gratefully acknowledged.
REFERENCES
1. R. G. Brusch, Trajectory optimization for the Atlas/ Centaur launch vehicle. J. Spacecraft Rockets 14, 550-555 (1977). 2. G. L. Brauer, D. E. Cornick and R. Stevenson. Capabilities and applications of the Program to Optimize Simulated Trajectories (POST). NASA CR-2770 (1977). 3. K. H. Well and S. R. Tandon, Rocket ascent trajectory optimization via recursive quadratic programming. J. Astron. Sci. 30, 101--116 (1982). 4. R. A. Tapia. Diagonalized multiplier methods and quasi-Newton methods for constrained optimization. J. Optimiz. Theory Appl. 22, 135-194 (1977). 5. S. P. Han, Dual variable metric algorithms for constrained optimization. SIAM J. Control Optimization 15, 546-565 (1977). 6. M. Bertocchi, E. Cavalli and E. Spedicato, Computational performance of diagonalized multiplier quasi-Newton methods for nonlinear optimization with equality constraints. In Numerical Optimization o/' Dynamical Systems pp. 247 268, (Edited by L. C. W. Dixon and C. P. Szego), North-Holland, Amsterdam, (1980). 7. V. Adimurthy, K. V. Joy and M. Karunakaran, PYOPT 81, a trajectory optimization software package. VSSCTR-13-231-81, Vikram Sarabhai Space Centre, Trivandrum (1981). 8. D. F. Shanno and K. H. Phua, Matrix conditioning and nonlinear optimization. Math. Programming 14, 149 160 (1978). 9. E. Polak and A. L. Tits, A globally convergent implementable multiplier method with automatic penalty limitation. Appl. Math. Optimiz. 6, 335-360 (1980). 10. V. Adimurthy, Optimality in PYOPT. SSTC-ARD-FEB (1)/82, Vikram Sarabhai Space Centre Trivandrum (1982). 11. M. S. Bazaraa and C. M. Shetty, The Fritz John and the Kuhn Tucker optimality conditions, Nonlinear Programming, Theory and Algorthms. pp. 147 150, John Wiley, New York. (1979). 12. M. Avriel, Penalty function methods. Nonlinear Programming, Analysis and Methods'. p. 401. Prentice-Hall, N.J. (1976).