STOP - A computer program for supersonic transport trajectory optimization

CONTENTS (Cont. ) Page Program Equations, Variables, and Constants APPENDIX B - Defined Equations of Motion Partial Derivatives 306 Auxiliary Printout Variables VAR Array CONST Array 316 Miscellaneous COMMON Internal Program Indicators 319 APPENDIX C -Program Control Logic vii STOP- A COMPUTER PROGRAM FOR SUPERSONIC TRANSPORT TRAJECTORY OPTIMIZATION By Lawrence H. Stein, Malcolm L. Matthews, and Joel W. Frenk Boeing Aerospace Group SUMMARY An IBM 7094 digital computer program using a steepest-ascent procedure has been developed for optimizing the flight path of supersonic transport aircraft from the start of climb through cruise and descent. This document describes the analytical development of the supersonic transport optimization program (STOP) and presents a manual for program users.

Program capability is summarized below.

The program is capable of optimizing the entire flight of an SST from a given 1.

low-speed, low-altitude condition at the start of the flight to a given low-speed, low-altitude condition at the end of the flight. The program will also optimize the climb, cruise, and descent phases of a mission separately.

2. The program incorporates a 3-D point-mass simulation for a vehicle moving in relation to a spherical, rotating earth. The inverse-square law for gravity is used. The 1962 U.S. standard atmosphere is included as a single subroutine.

3. The optimization is accomplished with an automatic step-size controller so that, in general, only one pass on the computer is required for a solution.

The program will generate a nominal trajectory for starting the iterative procedure. Automatic plotting capability is included. STOP is sufficiently versatile so that other vehicles besides the SST may have their flight paths optimized. These vehicles include: space boosters, ICBM’s, reentry vehicles, scramjets, and vehicles with air-augmented propulsion. Even problems com- pletely divorced from flight path optimization may be solved with a minimum of reprogramming.

4. Payoff functions, terminal constraints, and stopping parameters may be selected from the list of 40 state variables (defined by equations of motion and enroute placards) given in figure 1. The program will optimize any one of the state variables while simultaneously satisfying 14 terminal constraints (one of which is considered to be the stopping condition). Inequality con- straints, imposed by the user, are considered as terminal constraints.

Pitch angle, bank angle, wing sweep, and throttle setting are the control vari- 5.

ables for the program.

The user may select any subset of these variables for a given problem.

.

VARIABLE VARIABLE NC PLOT FORTRAN SYMBOL DESCRIPTION UNITS TYPE INDEX INDEX INDEX NAME STATE 121 1 X(Kl) W WEIGHT LB (EQUATION : 122 2 XlK2) H ALTITUDE Fr OF MOTION) 3 I23 3 xK3t Y RELATIVE FLIGHT PATH ANGLE DEG 4 124 4 X(K4) RELATIVE VELOCITY FPS 5 I25 5 XW) LATITUDE ANGLE DEG ; 6 126 6 XlK6) RELATIVE HEADlNG ANGLE DEG 7 127 7 X(K7) LONGIlUDE ANGLE !t" DEG a 128 a X(K8) TD DUMMY TIME SEC 9 9 X(K9) RNG PATH RANGE ALONG EARTH'S SURFACE N MI iii 10 10 XiKlD) AHI AERODYNAMIC HEATING INTEGRAL F-T-LB/FT2 11 11 X(Kl1) A V IDEAL RELATIVE DELTA VELOCITY FPS 12 z: X(K12) GL GRAVllY LOSS FPS E I3 X(Kl3) DL DRAG LOSS FPS 14 ;i 14 X(K14) TVL THRUST VECTORING LOSS FPS 135 15 X(K15) ER RELATIVE SPECIFIC ENERGY Fl ii 136 16 X(K161 137 17 X(K17) AVAILABLE FOR EXPANSION :i 138 la X(K18) NOT DEFINED 19 139 19 X(K19) 20 140 20 X(K2D) ‘I STATE 141 21 X(K21) 0" PITCH ANGLEPLACARD = Fit) DEG2SEC (ENROUTE 142 22 ii X(KiZ) q+ BANKANGLEPLACARD =F(t) DEG'SEC PLACARD) 23 143 23 X(K23) ?* THROTTLEPLACARD = F(H,M) SEC 24 144 24 X(K24) 11' WING SWEEP PLACARD = F(H,M) DEG'SEC 145 25 X(K25) (Y. ANGLE OFAllACK PLACARD - F1H.M) DEG2SEC iii 146 26 X(K26) NOT DEFINED 27 147 27 X(K27) HD' ALTITUDE RATE PLACARD = F(t) FPS ,PSF12 SEC 28 148 28 X(K28) Q DYNAMIC PRESSUREPLACARD = F(t) 149 29 XlK291 DYNAMIC PRESSURE PLACARD * F(M) (PSFI'SEC :i 150 30 X(K30) L Q"PlACARD = F(M) (PSF-DEGj2SEC I"R12SEC 31 31 X(K31) TEMT' TOTAL TEMPERATUREPLACARD - F(t) :: 32 32 X(K32) N' NORMAL LOAD FACTOR PLACARD = F(H,M) 33 X(K33) RPA" RESULTANTPHYSIOLOGICALACCEL. PLACARD = Fit1 ZSEC~I~SEC z 34 ALTITUDE PLACARD = F(M) FT2SEC ;'4 X(K34) H' 155 35 ~(~35) AP SONIC BOOMOVERPRESSURE PLACARD nF(A,B) (PSFI2SEC 156 36 MACHNUMBER PLACARD =F(H) SEC ;'6 X(K36t M" 157 37 xK37l NOT DEFINED ;'B Xw38l Y" GAMMA PLACARD =F(H,M) DEG2 SEC si ;; X(K39) NOTDEFINED * 40 160 40 X(K40) NOTDEFINED

L

-. -..-- - --_--. .- ~

FIGURE I: STATEVARIABLES

-.- . -. .---.--- _I... ._ 6. Aerodynamic and engine options are available for receiving the data in the forms commonly used for most types of vehicles.

7. Multistaged vehicles and those where external stores are jettisoned as a function of time may be optimized.

Inequality constraints may be imposed on parameters affecting geopolitics, 8.

passenger comfort, control limitations, structural loads, and engine opera- tion. The geopolitical constraints include sonic boom overpressure and maximum and minimum altitude.

9. The initial conditions can be varied by the program to obtain increased per- fo rmance . This option eliminates the need for a preliminary search to deter- mine the neighborhood of initial conditions for an optimal flight path.

INTRODUCTION The need to optimize transport aircraft flight paths is becoming apparent because (IL) the next generation will include supersonic transports for which the climb and acceleration phase consumes an appreciable part of the vehicle fuel, and (2) a large part of the flight is constrained by enroute placards.

The steepest-ascent method has been successfully used to optimize rocket- reentry-vehicle trajectories, and orbital transfers giving boost trajectories, There have been numerous attempts to optimize substantial performance gains.

flight paths of airbreathing interceptor-type aircraft with a relatively high thrust- to-weight ratio using the steepest-ascent technique (ref. 1). These methods have been successful to a degree, but, when applied to low-thrust-to-weight-ratio air- craft, have resulted in flight path instability problems that have made convergence difficult.

Recognizing that, in the past, optimization techniques have produced significant performance gains for many classes of vehicles, NASA initiated the present study The purpose of the study was development of a digital under RFP L-5347 (ref. 2).

computer program that would optimize the flight path of a supersonic transport with realistic operational constraints. A sufficient number of cases would be run to substantiate the program and define the most appropriate flight paths for selected SST configurations. A computer program general enough to optimize the SST would have the capability to optimize trajectories for many classes of vehicles including rockets, air-augmented rot kets, airbreather vehicles, and gliders.

The technique used for optimization is the steepest-ascent method given by Bryson (ref. l), which follows the direct approach of determining a maximum.

The optimum flight path for the present study is the solution to the nonlinear dif- ferential equations of motion that satisfy the imposed constraints and maximizes or minimizes one of the state variables.

The development of a computer program to satisfy the SST requirements and overcome problems associated with the flight path instabilities and complex engine characteristics required some changes in the usual methods used in optimization program s. Significant contributions to the program were made by NASA/LRC personnel. The use of pitch angle as a control variable was suggested by J. R.

Elliott as a technique to overcome the flight path instabilities. The formulation of jet engine data into a form compatible with the optimization procedure was suggested by W. E. Foss, Jr. This form of data uses thrust input as T/P and weight flow as W/P, where P is the ambient atmospheric pressure.

Contribu- tion to the form of the SST aerodynamic data representation was made by C. M.

Jackson, Jr.

SY!MBOLS A axial force (pounds) speed of sound (feet/second) a inertial acceleration vector (feet/second2) a indicial response at time t due to a unit step in f(t) at A (t - 7) time T (see Duham ells integral, equation 81) nozzle exit area (square inches) Ae aerodynamic heating integral (foot-pounds/feet2) AH1 component of acceleration normal to relative velocity vector an component of acceleration along the velocity vector “t component of acceleration completing the right-hand set aP (an, at, ap) = ii axial force coefficient, A/q S CA drag force coefficient, D/q S CD minimum drag force coefficient cDo lift force coefficient, L/q S CL (per degree) aCJJaa QJo!

lift force coefficient at Q = 0 CL,,0 lift force coefficient for minimum drag -!M normal force coefficient, N/q S CN acN/a@ (per degree) WY (per degree3) acN/i%t3 cNor3 thrust force coefficient, T/q (square feet) CT D drag force (pounds) drag constants (defined by equation 151) Do, Dl, D2 velocity loss caused by drag forces (feet/second) DL control variable perturbation magnitude (see equation 100) w2 df predicted change in a function f (see fig. 8) trial value for df df0 df* actual change in a function f (see fig. 8) i trial value for df* df,” I ‘_” - -___ .- combined changes in the constraints and initial state variables

dB

change in payoff, constraint, and stopping functions

dtD> CG, dfi

trial values for d@ and dT &no, dii, n x n matrix of partial derivatives, .afi/axj

F, F(t)

an n x 1 matrix of functions defining the time derivatives f, f (x(t)> u(t), of the state variables a value determining the degree of nonlinearity of the func- fNL tion f forcing function of Duhamel’s integral (see equation 81) f(t) n x m matrix of partial derivatives, 6fi/6uj G, G(t) local acceleration of gravity (feet/second)

g

forcing function of the convolution integral (see equation 99)

g(t)

gravitational acceleration at earth’s surface (feet/secondz)

go

velocity loss due to earth’s gravitational field (feet/second) GL GF function defined by equation 140 h altitude above earth’s surface (feet) reference altitude used for skin friction drag calculations (feet) href response of a system at time t, due to a unit control varia-

h (t 7)

ble impulse at time 7 (see convolution integral, equation 99) specific impulse, Tfi (pounds thrust/pounds fuel per second) IsP integrals defined by equations 108b, c , and d bP

Idw’ Iw

k step size coefficient (see fig. 8) fuel flow correction factor for rockets Km constant used in sonic boom calculations KR total impulse correction factor for rockets KT lift force (pounds); reference length for sonic boom calcu- L lations Mach number M m number of control variables N force normal to body axis (pounds) n number of state variables plus enroute placards number of terminal and enroute constraints P sonic boom overpressure (pounds/fee@) AP ambient pressure (pounds/feet2) p, p* G (pounds/feet2) pref sea-level ambient pressure (pounds/feet2) PSL dynamic pressure, l/2 p V2 (pounds/feet2) Y dynamic pressure times angle of attack (pounds-degree/feet2) YQf r position vector measured from earth’s center to vehicle (feet) desired final altitude for circular orbit (feet) ‘D RPA resultant physiological acceleration (feet/second2) radius of a spherical earth (feet) Ro aerodynamic reference area (square feet) S SL6 integrated payoff function sensitivities (see equation 134) u CD instantaneous payoff function sensitivities (see equation 130) su specific fuel consumption (pounds fuel per hour/pounds thrust) SFC actual. thrust component along body axis (pounds) ; static T temperature (degrees Rankine) ; final trajectory time (seconds) t independent trajectory variable, time (seconds) initial trajectory time (seconds) t0 dummy time (seconds) tD uncorrected thrust as input for airbreathers (pounds) Tref total corrected thrust for rockets (pounds) TT total temperature (degrees Rankine) Ttota1 T uncorrected vacuum thrust as input for rockets (pounds) vat velocity loss caused by thrust vectoring (feet/second) TVL u, u(t) an m x 1 matrix of control variables 7, v, inertial velocity vector (feet/second) circular satellite velocity (feet/second) vcs relative velocity vector (feet/second) TR function defined by equation 150 VF W weight (pounds); an m x m weighting matrix elements of the weighting matrix W (equation 137) wii fuel weight (pounds) Wf path range (feet) X coordinates of points on the forward shock signature (see x9 Y fig. 7) - - an n x 1 matrix of state variables and enroute placards x.3 x(t) lateral extent of forward shock signature (see fig. 7) yc angle of attack (degrees) angle of attack of the previous integration step (degrees) QO latitude of vehicle (degrees) B latitude at which the overpressure intersects the ground BG (degrees) (see fig. 7) latitude of a point on the forward shock signature (see fig. 7) Blocal dix-l Bm flight path angle (degrees); specific heat ratio Y Kronecker delta; indicates variance thrust cant angle (degrees) % throttling control variable r) pitch angle ; control variable, y -t a (degrees) sweepback control variable (degrees) A longitude of vehicle (degrees) h adjoint variables A, x(t) longitude at which the overpressure intersects the ground XG (see fig. 7) longitude of a point on the forward shock signature (see fig. 7) hlocal payoff function adjoint variables, measures sensitivity of XQ m xQy Q at time T, to state variable changes at time t constraint function adjoint variables, measures sensitivity

A@, x$(t)

of a constraint at time T, to state variable changes at time t stopping function adjoint variables, measures sensitivity of xgp An(t) the stopping function at time T, to state variable changes at time t see equations 98a and b Lagrange multiplier (see equation 112); first Lagrange multiplier (see equation 713) f a 1 x p matrix of Lagrange multipliers (see equations 108 v, vs.

and 114) ambient density (slugs/foot3) P dummy time used in Duhamel’s and the convolution integral payoff function; bank angle (degrees) Q heading angle, measured in degrees north of east @ a p x 1 constraint vector T sz stopping function angular velocity of earth (radians/second) w Subscripts 0 sea level or earth surface; initial; variable based on the angle of attack of the previous integration step and current state variables refers to relative coordinate system R I refers to inertial coordinate system maximum max required req vat vacuum ref reference m ambient condition Superscripts transpose of a matrix ( ) ( )’ inverse of a matrix ( ) ( j-l indicates nominal value ( )* - - ANALYTICAL DEVELOPMENT Equations of Motion Coordinate system. -The flight-path coordinate system is used to define the vehicle position and atitude with respect to a spherical, rotating earth. This sys- tem is sometimes regarded as the Watural” coordinate system because the state variables of flight-path angle and velocity are explicit coordinates. The basic state variable coordinates are longitude (X), latitude (fl) , altitude (h) , relative velocity (VR), relative flight-path angle (YR), and relative heading angle (#R) (fig. 2).

Control variables. -The total force acting on the vehicle has three distinct sources: (1) gravitational force as a result of mutual mass attraction between the vehicle and earth; (2) aerodynamic force resulting from the vehicle motion through the atmosphere; and (3) thrust force from the vehicle propulsion system. There is no way of controlling the gravitational force since it is a function only of the state of the system.

The aerodynamic forces are determined by the geometry of the vehicle and its attitude with respect to the free-stream air mass. For fixed-geometry aircraft, the aerodynamic force is dependent on the angle of attack (defined as the angle be- tween the vehicle longitudinal axis and the velocity vector). The angle of attack is used primarily to establish lift that is normally in a vertical plane. Out-of- plane maneuvers are made by banking the vehicle to direct the lift vector out of the vertical plane, thus permitting lateral translations. The bank angle, Q , is defined as the angle between the vertical plane containing the velocity vector and the vehicle plane of symmetry as viewed along the velocity vector (fig. 3).

A problem exists with angle of attack as a control variable for the point mass equations of motion when simulating the flight of a low thrust/weight, winged, air- breathing vehicle. The angle-of-attack control variable has been satisfactorily used for high thrust/weight vehicles, but the use of angle of attack for SST-type vehicles results in an unstable, oscillatory flight path during cruise which is not representative of the actual path. The use of a pitch angle (8) given by e=y+a!

(1) as a control variable resolves the problem and results in a stable, well behaved flight path in all cases checked. The analysis of this problem is given in appendix A. Note that if the vehicle is banked, 8 is not the pitch angle in the usual sense (that is, the angle between the longitudinal axis and horizontal plane), but instead is the algebraic sum of the flight-path angle and the angle of attack as given by equation 1. Therefore, for nonzero bank, 8 is not a physical angle.

Since some of the more advanced aircraft permit changing vehicle geometry by varying the wing sweep, the effect of this variable must be considered in !

--,--,e- ,.

.I : ,_ .I ,..;” i, ,,,/ --- ..:’ ,. .‘: , ~’ ::; J .. ‘.,P ‘_ I

VR

HORILONTAL

figure 2. BASIC COORDINATE SYSTEM

Figure 3. VEHICLE ORIENTATION

The wing sweep, h, is de- establishing the aerodynamic force on the vehicle.

fined as the angle between the lateral body axis and a reference line on the wing (e. g. , leading edge, quarter chord, and half chord).

The remaining control variable to be considered is the throttle setting, r) , for controlling the thrust. This parameter is defined as the ratio of thrust at a given state to the reference thrust at that same state. Thrust is considered to act along the vehicle body axis.

Therefore, to determine the aerodynamic and thrust forces required for con- trolling the vehicle, four control variables are required: 8, Q , h, and q.

Kinematic equations. - The basic state equations of motion are defined in this section. Details of the development are not given here because the methods are standard.

The kinematics of a point mass vehicle for a rotating, central-force field is given by where V is the inertial velocity vector vR is the relative velocity vector z is the earth’s rotation rate r is the position vector of the vehicle a is the inertial acceleration vector.

Substituting equation 2 into equation 3 gives diT, -+2wxVR+WX(wXr) ‘= dt Using the coordinate system of figure 2, equation 4 can be written in corn - ponents along the flight path, normal to the flight path in the plane containing the r and the v vector, and perpendicular to the plane containing r and v.

at = +R + a2 r ‘cos 6 (sin fl cos yR sin $R - cos fi sin yR) vR2 =VR~R-~COSYR-~UVRCOS~COS$R an - u2 r cos fl ,(cos p Cos yR + sin fl sin yR sin QR)

(5)

vR2 r tan fi cos2 yR cos #R ap = VR cos YR & f - 2 cd VR (sin yR sin $R cos 6 - sin fl cos yR) + a2 r sin p cos fl cos qR YR is the relative flight-path angle where qR is the relative heading measured north of east.

Other equations which follow from the procedures are vR cos yR sin #R j= r VR cos YR cos #R >;=

r cos /3 (6)

i-l = vR sin YR . RO x= c”s YR vR r X is the longitude angle where h is the altitude x is the path range.

The change in weight of the vehicle is given by the differential equation

&-

weight flow

(7)

where the weight flow is a combination of fuel flow and inerts.

and thrust forces are required to complete The gravitational, aerodynamic, the system. Other “auxiliary” state variables that are required will be defined in a later section.

Applied forces. -The forces acting on the vehicle which are not a result of , the vehicle kinematics are discussed in this section.

I i i I ----.._ I :, _,_ :.

. ‘, Gravitational forces: The gravitational forces are determined for the cen- tral force field by representing the gravitational acceleration by the inverse square law

(8)

Aerodynamic forces: The aerodynamic forces for aircraft are written for the wind axes system with lift (L) normal to the flight path and drag (D) parallel to the flight path. Side-slip is assumed to be zero, i.e. , all turns are assumed to be coordinated such that no side-slip results.

Drag Figure 4.

AERODYNAMIC FORCES- WIND AXES

If the bank angle is not zero, the lift vector is rotated about the velocity vector through the angle Q (bank angle).

The aerodynamic forces may be expressed as L=cLqs (9) D = CD q S (10) L is the lift force where D is the drag force CL is the lift coefficient CD is the drag coefficient y?m 2 q is the dynamic pressure = i p V2 or - M S is the reference area.

Aerodynamic data input options are available for accepting data in the most common forms. For vehicles with fixed wings, the data are considered as C L = f (h, M, q) (11) CJ-j = f 01, M, a) (12) if not represented in polar form. If available as drag polars, the data are repre- sented as CL = CL, 01 (13) aCD CD = CD0 + -

(14)

CL2

acL2

where = lift curve slope, f (M) per degree % = minimum drag coefficient, f (M) cDO aCD - = induced drag constant, f (M).

acL2 The aerodynamics for variable-geometry aircraft with wing sweep as a con- trol variable are formulated as CL = f (A, M, a) (15) CD = f (A, M, a) (16) or if the data are given as drag polars (17) CL = CLo,.-o + CL, cx aCD CD = CD0 + -

(‘L - ‘LMj2 (18)

8cD2 where = lift coefficient for Q = 0, f (A, M) CLQ=O = lift curve slope, f (A, M) per degree

%Y

CD0 = minimum drag coefficient, f (A, M) acD - = induced drag constant, f (A, M) acL2 = lift coefficient for minimum drag, f (A, M).

cLM The effect of altitude on drag (skin friction, AC,,) is included as a linear function of altitude, *CD = K fh - href) (19) A ‘Do where &CD/Ah is a function of the reference altitude.

The aerodynamic forces on a missile are quite often referenced to the body axis with the normal force perpendicular to the body axis and the axial force along the body axis.

N=cNqs (20) A = CA q f-3 (21) where N is the normal force A is the axial force CN is normal force coefficient CA is axial force coefficient.

FlGURE5: AERODYNAMIC FORCES - BODY AXES

The normal and axial force coefficients are given by CN = cNa Q! + CNa3 a3

(22)

=C Ao (not dependent on Q)

(23)

cA = f (M) per degree where WY CQ~ = f (IM) per degree3 = f (M).

cAO The forces in the body axis system are transformed to the wind axis by = cN cos o! - CA sin 01 (24) CL CD = CN sin a + CA cos a (25) for use in the program.

Thrust forces: The thrust forces on a vehicle are considered to act in a direction fixed relative to the body axis so that the resultant thrust along the body axis is T = TT COS 6, (26) where T is the thrust along body axis TT is the total thrust force 6, is the thrust cant angle.

(The angle 6, is used only for symmetric nozzle configurations where the normal thrust component from two opposing engines cancel.)

T,

Figure 6. THRUSTREPRESENTATION

The computation of thrust for airbreathing engines follows the form

(27)

TT = Tref q Tref is the reference thrust = F (cy , h, M, or V) where TJ is the throttling parameter.

The thrust may also be given by T

ref = cT q (28)

where CT is thrust coefficient. = F (dl , h, M, or V) A thrust option that has particular use for the jet engine is T ref

(29)

- = f (O!, h, M) P which was formulated to give a smooth, consistent set of data over the entire operating range of the engine. P is the atmospheric pressure (psf).

vacuum conditions. The The reference thrust for rocket engines is based on total thrust is given by T? = KT TV,, r) - P, A,/144 (30) where Tvac is the vacuum thrust (pounds), f (time) P, is atmospheric ambient pressure (psf) A, is the exit area (square inches), f (time) KT is the total impulse correction factor.

Fuel flow: The fuel flow for airbreathing engines can be specified in terms of fuel flow directly or in terms of thrust specific fuel consumption.

The fuel flow and specific fuel consumption can be expressed as tif, SFC = f (q, h, V) = f (o!, h, V) (32) = f 01, M, h) = f (a!, M, h) depending on the independent variables required for adequate representation.

The fuel flow for jet engines is also formulated as .

wf -= fth, M, q (33) P and is designed to give a consistent set of data from idle speeds through maximum augmentation. Again P is. the atmospheric pressure (psf).

For rocket engines, the fuel flow is obtained from the specific impulse by the relation (34) Km is fuel flow correction factor.

where Isp is specific impulse (seconds), f (time).

The options of thrust and fuel flow available will handle most types of engine data formulations.

Basic equations of motion. -The basic equations are summarized below with the effects of gravity, aerodynamics, and propulsion included.

T cos cc - D - g sin yR + a2 r cos B (cos fl sin yR SR = m - sin fl cos YR sin $R) (35) T sin O! + L) cos cp _ + 2w cos fl cos GR ‘OS YR +R = ( m vR n 0’: r + - cos p (cos p cos yR + sin @ sin YR sin $R) (36) VR k = VR sin yR (37) .

W = - tif - tiinert (38) $R = _ (T sin cy + L) sin ge vR - y tan p cos yR cos GR m VR cos yR - 2 w (sin fl - tan yR sin $R cos 8) W2 r sin /3 cos /I3 cos #R, (3% VR co.5 ~ VR cos yR sin @R (40) r VR cos YR cos GR (41) r cos /3 These equations are sufficient to define a vehicle trajectory for a spherical, rotating earth. The equation of motion for path range on the earth’s surface is RO j,=- r VR c”s YR (42) Auxiliary state variables. -Variables required for optimization, stopping or constraining a flight path, but which do not fall in the category of basic state The differential equations variables, are defined as auxiliary state variables.

The additional for these variables are integrated when desired by the user.

variables and their defining equations are given below.

The dummy time, introduced so that it may be used as a Dummy time: state variable, is defined by (43) The aerodynamic heating integral, a meas- Aerodynamic heating integral: ure of the heat encountered by a vehicle during flight, is defined by dAHI = $ p VR3 (44) dt where pis the atmospheric density.

Ideal relative AV: The ideal relative AV, the velocity change a vehicle would experience in the absence of atmospheric and gravitational effects, is defined by dAV = 2 + 2 1: cos /3 (cos /3 sin y - sinBcosy R sin 9,)

(45)

R dt where m is the vehicle mass.

Drag loss: Drag loss, the velocity decrement resulting from motion through the atmosphere, is defined by dDL D -=- (46) dt m Gravity loss, the velocity decrement resulting from motion Gravity loss: in the earth’s gravitational field, is defined by d GL - = g sin yR (47) dt Thrust vector loss: Thrust vector loss, the velocity decrement resulting from the thrust vector not aligned with the flight path velocity vector, is defined by dTVL (48) dt Inequality constraints. -During the course of a flight path, a function may be required to be less than or equal to and/or greater than or equal to a particular value or set of values. Inequality constraints arise from a number of different considerations: (1) geopolitical limitations; (2) passenger comfort; (3) control limits; (4) structural limits; and (5) engine limits.

Geopolitical: This class of inequality constraints is related to consideration of other than aircraft limitations. The maximum and minimum altitude may be imposed by the FAA for traffic control. Aircraft noise and sonic boom (ref. 7) are important in vicinity of cities, airports, etc. The sonic boom is of special consideration to SST flight paths because the boom signature and overpressure (AP) on the ground may design the flight path. The sonic boom overpressure on the ground, in the plane of the velocity vector, may be found from the following expression: 3/4 AP @I/L) = f (M, h, pm CLS/ZL~)

(49)

KR 8m1'4 'ref where P = 21 P ref SL poJ = 1.9 KR and f (M , h, &-, CL S/2L2) is a tabular function. The overpressure generated by an airplane at a latitude and longitude (p , A) will intersect the ground ahead of the airplane in the plane of the velocity vector at a latitude and longitude of

R (50)

m 0 h cos QR (51) ‘(2 = A + pm R. cos fl The shock pattern on the ground is considered to be the intersection of a cone (with the half apex angle equal to the complement of the Mach angle) with the ground plane (fig. 7). This intersection, or shock signature, is a hyper- bola with the equation h2 + y2 =x

(52)

where h is the vehicle altitude y is lateral distance from the flight path x is distance along the flight path from the vehicle to the point on the shock signature.

9 -/

’ ‘\

1 \

/

/

“I \

/

/’

(“LOCAL, pLOCAL)

/

/

/

/

/I(W)

Figure 7: SONIC BOOM GEOMETRY

The lateral extent of the shock wave is given as a tabular function = f 01, M) (53) yC The geographical location of a point along the shock signature is expressed in term s of longitude, latitude , heading angle, and altitude =fi+$ o si;GR KyT+ y cos GR]

B (54)

local m [ =A+ IL ‘OS ‘R {w - y sin GR (55) x1ocal Rocos P [ B m 1 where y can go from - yc to + yc to cover the entire shock signature.

The overpressure along the shock signature is related to that given in equation 49 (APy=6) by h AP(Y) = Apy=6 (56) J/m

[ I

The maximum allowable overpressure over the ground is given by = AP AP

(57)

max max ( local’ ’ local) are the points along the shock signature.

where hocal and ~local

The next step is to derive a function in terms of the instantaneous vehicle state that represents the integral of the overpressure violation along the shock The overpressure violation is given by signature from - yc to + yc.

AP -Ax’ ‘o S max where the AP, is given by equation 56.. A true integral of the violation along the signature is not practical, therefore an approximation is made by summing the overpressure violations for several selected points as follows.

d BOOM PF (AP - AP )

(5 9)

dt = S max i 6 = 0 i.fAps<Apmax 6 = 1 ifAps ?.APmax At present 7 points NP is the number of points selected for the approximation.

Equation 59 which defines the sonic are used spaced at an interval of yc/3.

boom penalty function, is integrated as one of the equations of motion along the path with a terminal constraint value of zero. The partial derivatives of equa- tion 59 with respect to the state and control variables are determined numerically (appendix B) .

The comfort of a passenger is significant to the air- Passenger comfort: lines who are trying to please the customer. To insure that the vehicle does not fly into areas of discomfort, placards may be imposed on normal load factor (n), pitch angle ( 6)) bank angle (cp), resultant physiological acceleration (RPA), and The altitude rate controls such items as rate of change of cabin altitude rate (h).

pressure. The magnitude of cabin pressure is determined by altitude and mechan- ical or compressor limitations.

Control : These placards result from physical or mechanical limitations on the system. The control variables limited are pitch angle, bank angle, throttle setting, and wing sweep.

Structural: The structural requirement on airframe requires that the state of the system be constrained so that dynamic pressure, dynamic pressure times angle of attack (Qa), stagnation temperature, and cabin pressure differential do not exceed prescribed limits. The stagnation temperature is given by the adia- batic equation y-l 2 T = T (1 + -y-M) total where y = specific heat ratio and T is the ambient temperature. Cabin pressure differential can be formulated as a placard on altitude rate.

Engine : The engine is protected from flight conditions that would produce undesirable results by bounding the altitude Mach number region in which the vehicle can operate. In addition, maximum and minimum throttle setting may be bounded by a control placard as a function of altitude and Mach number.

Several methods of handling inequality constraints have been devised, i. e. , integral method (ref. 3)) penalty function method (ref. 4)) and a method by Bryson, Denham, and Dreyfus (refs. 5 and 6). The method used during this study is based on the square of the violation. This method is similar to the integral method but produces a form more satisfactory for use with analytical partial derivatives. Consider a typical problem as shown in figure 8.

Let C be the time history of the function to be constrained, which violates the placard shown as P. The region of constraint violation is shaded. A meas- ure of the total constraint violation is the areas shaded, or T A= (C - P) 6 dt (61) s to

I --

where 6 = 0 when there is no violation 6 = 1 when there is a violation.

This formulaf;ion is not satisfactory for analytical partials since the partial de- rivative of dA/dt with respect to the state or control variables does not have a measure of the violation. A form which does have the required nature is A= (C - P)2 6 dt (62) where the partial derivative with respect to the state variables does have the proper form

ai 6C 6P

- = 2 (C - P) 6j - e

(63)

8X

( )

which is still dependent on the violation.

Each of the inequality constraints is formulated in the same manner and permits a minimum as well as a maximum limit. Each is treated in the same manner.

Inequality constraints can be imposed for the placards shown in figure 1.

PLACARD VIOLATED -P I t

T

I ‘0

Figure 8: INEQUALITY CONSTRAINT

Steepest-Ascent Technique The analysis of trajectories by the steepest-ascent method has been treated extentensively in the literature (refs. 1,. 3, and 4). The derivation given here follows that of Bryson and Denham, reference 1. Since the steepest-ascent method is basic to this study, the derivation is included here for completeness.

- The steepest-ascent solution is formulated to deter- Statement of problem.

so as to maximize a function (payoff mine the control variable history, u(t), function) @ = G+(T), T] (64) and meet the terminal constraints c, = +(T), T] = 0 (65) where cp = performance index to be optimized x = state variable matrix (order n x 1) T = time at the stopping condition $ = a matrix of terminal constraint functions (order P x 1).

The stopping time, T, is established by the condition St3 = s+(T), T] = 0 o-33 where a is referred to as the stopping condition.

The equations of motion defining the state of the system can be written in terms of first-order differential equations. The form of these equations is x(t) = T [x(t), u(t), t] (67) where i(t) is an (m x 1) matrix of control variables f is an (n x 1) matrix of functions defining the time derivatives of the state variables.

The initial values of the state variables, x(t,), are generally known for a given problem. If some of the initial values are not specified, they may be deter- mined along with the control matrix to optimize ge. If it is desired to optimize some function other than one of the normal state variables, it is necessary to introduce an additional state variable and another differential equation defining the state variable and solve it simultaneously with the required equations of motion.

The steepest-ascent process is started by _assuming a control variable time history, u(t). About the only requirement on u(t) is that the trajectory generated must reach the stopping condition.

I -

The nominal trajectory is generated using the assumed u(t), vehicle and engine characteristics, and the initial conditions. The equations of motion (equation 67) are integrated numerically to obtain the nominal state variable program until the For convenience, the nominal conditions will be stopping condition is satisfied.

designated by ( )*, i.e., i(t)nomiml = ;*lt) w9 -* XWnomiml = x @) . (69) The steepest-ascent procedure uses the concept of local linearization about a nominal flight path. For small perturbation in the control variables u(t) = ii*(t) + s;(t) (70) the state variables are perturbed x(t) = x*(t) + 6x(t) (71) Substituting equations 70 and 71 into equation 67, the linear differential equa- tion for 6; is obtained -$ 6x(t) = F 6; + G 6; (72) w’her e (order n x n) (73) G= (order n x m)

(74)

Again the ( )* indicates that the partial derivatives are evaluated along the nomi- nal flight path.

Adjoint equations. -The effects of the control variable perturbations on the payoff quantity, stopping condition, and constraint functions must now be determined.

Equation 72 is a linear equation that describes the small perturbations about the nominal path. To this is added a linear differential equation adjoint to equa- tion 72 dX - = - F’ X (75) dt where the X Is are a matrix of adjoint variables. The significance of the X’s and their role as influence functions for Q , T, and S2 will be demonstrated.

To show this, premultiply equation 72 by X ’ , equation 75 by 6 x’ (where ( ) ’ indicates the transpose of the matrix), and add the transpose of the second product to the first 1 do; dh’ x - + --g 6,; = A’ F 6;; + A’ G 6; - A’ F 6; (76) dt which can be written as (77) Integrating equation 77,

(78)

Examining this equation gives much information on the nature of the adjoint vari- The product G 6U gives the rate of change of the state variables due to ables.

changes of the control variables only: 6x = G 6u (79) therefore, equation 78 can be written T .

(ii’ 6ii)t=T = A’ 6; dt + (X’ 6x)t=to

(80)

s

t0 but this is exactly the form of the Duhamel integral, which is

t

z(t) = A (t - 7) y dr + A (t

- to) Wo) (81)

S

t0 where A (t - T) is the indicial response at time t due to a unit step in f(t) at time 7. The f(t) is the arbitrary forcing function; z(t) is the response of the system to the arbitrary forcing function. Making a comparison between the two equations for t = T, the X’s are the indicial response of the system due to a unit step in each element of the state variable vector, 6 2. The variables of interest for this work are the payoff function cp , the stopping condition Sz , and the terminal constraints T. A set of X’s can now be defined corresponding. to Q, a, and the vs. The adjoint variables are identified as a(8

(82)

A’Q= 23

6X

and X’Q isa (lxn) is a (p x n) matrix, is a (1 x n) matrix, X’ where A’ he X’s are influence fun&ions and may be considered as !f matrix.

(83) ‘Qi = 2 Similar interpretation can be or the change of Q due to a unit step change of Xi.

A set of boundary conditions can now be written for the made for X9 and X0.

X’s at t = T giving j = 1, . . . . . n (84) XQj(T) = OGi j=l,...., n (85) i = 1, . . . . . p Ai (T) = - % axj T j = 1, . . . . . n (86) 1f the payoff, constraint, and stopping functions are restricted to state variables, the X(T)‘s are either zeroes or ones.

Starting from the terminal values of the X’s, equation 75 may be integrated from T to to) uying martial de’rivatives which backwards along the flight path (i. e. , were stored during the forward trajectory to obtain X (4 A’$, and X 0 as functions of time.

84, 85, and 86 into equation 78 and Substituting the functions from equations noting that

(h’Q 6;)~ = 6@(T)

(87) (A’$ 6G)T = 6F(T) (88) 6i)T = sa(T) (89) (‘C’n the following are obtained ~Q(T) = IT X’Q G 6; dt + (X’Q 6$t=t, t0

(91)

T G 63 dt + (x’o 6;;,t=t

6n(T) = (92)

X1, These equations would give the changes in the payoff function, stopping conditions, and terminal constraints if the final time (T) at the stopping condition did not due to small perturbations, the stopping change. In general, this is not true and, To account for this change equations 90, 91, and time is at some time, T + AT.

92 become dQ = G 6; dt + (X’Q 6&, + Q(T) AT (93) X’Q T G &dt+ (A’ d+ = 6$+, + z(T) AT

(94)

Q t0 .

dC-2 = A’, G 6; dt + (X’a d$to + n(T) AT

(95)

This set of equations gives Note that these equations are now total differentials.

and Sz for the perturbed trajectory.

the changes in (8, F, The stopping condition does not change from one iteration to the next, there- fore dsI (T) = 0, giving the relation G 6; dt - - (96) b2l fiZ(T) the equation for Q and 7 from equations 93, 94, and 95 can now be written as dQ = G 6; dt + (X’Qsl 6&,

(97)

xlQS2 where

@a)

Pb)

are modified influence functions, which account for a change The xQ(j- and X9n in the time at the stopping condition.

The integral terms of equations 97 and 98 may be compared with the convolu- tion integral used for determining the response of a system to a continuous control variable forcing function by first obtaining the impulse response functions, i. e. , + L z(t) = h (t - T) g(r) d7 (99)

S

t0 where h (t - T) is the response of a system at time t due to a unit control vari- able impulse at time T. The function g(T) is the continuous control variable forcing function and z(t) is the response at some arbitrary time t to that forcing function.

If the arbitrary time is taken to be the final time T, and the forcing functions are the perturbed control variables 6;) the functions X ‘Qa G and x ‘Q 0 G are shown to be the responses of the final values of the payoff function and terminal constraints to a unit impulse of each of the elements of 6u at some time t.

Variational equations. -For steepest ascent, the control variable history, Si, that maximizes dQ given by equation 97 for a given dF and dSZ = 0 is de- sired. An additional requirement is that dU defined from (dU)2 = JT 6$(t) W s;(t) dt to be chosen to ensure that the perturbations, 6<, will be small enough for the linearizations leading to equations 84, 85, and 86 to be within reason. W is an (m x m) symmetric weighting matrix chosen to improve convergence during the steepest-ascent procedure. The automation of logic for the selection of (dU)2 and the weighting matrix are discussed in later sections.

The values of dv are chosen to bring the nominal solution closer to the de- sired terminal constraints, T-= 0. It is desirable, in many cases, to restrict the change of constraints, d# , allowed for each iteration.

The procedure now follows the calculus of variations using the method of Lagrange multipliers. A linear combination of equations 98 and 100 with equa- tion 97 gives T dQ = G 6u dt + A’,, o (t 1 6;; Go) X’Qi-2 s t0 T G 6; dt - X’en (to) 6; (to) dT - + v’ &&2 J [ “0 T &ii’ W 6u dt

(101)

I t0 or combined gives T -p 6&v 6udt dQ = XIasl G - v’ dqn G S[

to

(to) -v’ X1@ (to> (to)

+ v’ dF -c c1 (dU)2

(102)

where v is a (1 x p) matrix of constant Lagrange multipliers, and P is a con- stant . Both u and p can be chosen for convenience, since they were assumed arbitrary.

Taking the variation of equation 102 with respect to the control variable gives T - Cc 62,’ w 6; dt (103) G - P 6i’ W) S2u G - v’ XfQsl 6(dQ) = (X’QS2 S[ t0 and dU are considered to be constants. Since this is a where 6; (to), dF,, scalar equation s21? w 6; = 6U’ w ti2u (104) i.e., the transpose is equal to the matrix where W is a symmetric matrix.

Therefore the expression for the variance of dQ is reduced to T 2 Cr 6ii’ w tj2ii dt v’ X’qa G -

X’QQ G - (105)

&dQ) = J[ t0 The optimal dQ will occur where its variance is zero. Therefore, since b2u # 0, A’& - v$,G - 2p &iii’ w = 0

(106)

Solving for 6ii gives LW 6U = -’ G’ (XQn - $,V)

(107)

2P -1 Substituting the expression indicates the inverse of a matrix.

where the ( ) back into equation 98 and solving for V gives where dj=d5,- A' (108a) -1 (108b) XiaGW G’ $dt T -1 (108~) I+~= to xlllinGW G’ A&t s T -1 (108d) I A’ G’X dt lpnGW QQ = Qn

S

t0 In the computer program, reference is made to the ‘It1 matrix, which is defined as h’ GW-1 G’ h dt

( 109)

which may be represented as a partitioned where A is the n x (p + 1) matrix, matrix as

(110)

nxl nxp The I matrix can be partitioned as follows I I I QQ IQ@ I --- _----- I=

(111)

I I$Q Iw I [ symmetric, where

This is a symmetric matrix since IQ$ = I’$$,@and I @+ is

are defined by equation 108b, c, d.

lQQ’l$Qf and1$$ Using the expression for v (equation 108) and the equation for 6< (equation 107)) substituting into equation 100 and solving for p gives -l l/2 a-

(112)

The Lagrange multipliers referred to in STOP are of different form than those given by equations 108 and 112. For convenience, the “‘first Lagrange multiplier” in the program, /l,, is l/2 (d u)2’ - d6 ’ I$@ dfi p, +=* (113) IQQ - I&o ‘s$k I$Q [ and the array of Lagrange multipliers are denoted by

( 114)

v

%=-G =-%

The form of the perturbed control variable (equation 107) therefore can be written as ( 115 ) The complete expression for 6; can now be obtained by substituting the equations for ,u and v into equation 107, which gives

hi(t) = zkw-l G’ (XQn - X

@ Ii; ‘@Q) -1 ’ +w G XQn I;‘$ dB (116) This equation is the foundation for changing the control variable to improve the performance for each iteration. The payoff function change predicted for an iteration is obtained by substituting equation 116 for 6U (t) into equation 97 for d 0, giving l/2 dll) = f - d;’ I,; d/3 IQQ - I;icaI-;/,I+Q -I( + I;iQ I,; dP + X;b#,) 6% (to) (117) where the + sign is used to increase Q and the - sign to decrease Q. This pro- cedure will give the optimum trajectory as discussed under the statement of the problem. The method for perturbing the free initial condition 6ii (to) will be discussed in a later section. The e_ffects of varying the initial conditions enter the problem through its effect on d @ as well as 6ji (to).

Automatic Convergence -A major problem encountered in steepest-ascent cal- Iterative procedure.

culations is the choice of the step size, (dU)2. Too small a step results in regu- lar, but slow convergence; too large a step results in irregular progress and Therefore, to insure convergence in a minimum probable convergence failure.

amount of computer time, it is desirable to automatically select the largest step size possible, per iteration, consistent with the linearity requirements of the The selection technique used in this program is pat- steepest-ascent technique.

terned after a method developed by D. S. Hague (ref. 3).

A function fi (a payoff or constraint function) rarely is -nor does it have to The be -perfectly linear to be acceptable for steepest-ascent calculations.

degree of linearity or nonlinearity is, however, important and is defined as fNLi = Jdfi*iiai 1 018) where fNLi is the degree of nonlinearity representing the percentage error of the linear prediction to a parabolic approximation of the actual change of a pay- off variable or a constraint. dfi* is the actual change in fi and dfi is the linear prediction of the change in fi. Both the degree of linearity and choice of the func- tio,l used to measure linearity must be chosen with care and are discussed below.

On a perturbed trajectory the predicted payoff function change is given by Essentially the step size may be defined as a choice of (dU2) and d$. Con- sider the one parameter set of perturbations (dU)2 = k2 (dUo)2 WO) d$i = k d$oi (121) and d~oi are arbitrary nominal changes. It follows on substituting where duo2 equations 120 and 121 into equation 119 that W (k) = k dq,

022)

Here dcP, is the change in performance resulting from the nominal step, (dUo)2 which results from the choice of k = 1. Consider figure 9 where and d&i, the nominal step is denoted by the “trial” step and fi is any of the functions Cp or iJ.

Predicted linear change C .-o ;j bolic approximation 5 dfo i 2.0 .5 I .o paiameter k.

Stepsize - Fiqure 9. PARABOLIC APPROXIMATION The predicted linear variation with k is the straight line through the origin For a small enough step the actual nonlinear change and the point (l. 0, df,).

As the step size increases the two results will substantially agree with this line.

Suppose a parabola is fitted through the origin will, however, tend to diverge.

and the point (1.0, df,? having the same slope at the origin as the linear change.

The equation of this parabola is

df” = @fo* - dfo)k2 + dfok

023)

Provided a trial step is reasonably chosen, this expression will approximate the actual nonlinear change over a considerably greater region of step sizes than did the linear prediction. Essentially, the combination of a trial result and the linear prediction allow a second-order approximation to the variation of a function with step size.

Assuming for the moment that df, > dfo*, and df, (0)/k > 0, we can antici- pate from equation 123 that the greatest change in f will occur when df0 k M’

024)

2 (dfo” - dfo) Further, on substituting this result into equations 123 and 118, it follows that will be about 0.5. If the actual nonlinear at this point the nOIdineaI?ity, fNL, in variation were to be parabolic this value would be the desired nonlinearity.

Using 0.5 in conjunction practice, it has been found that this value is too high.

with the parabolic approximation leads to steps that are “over the hump. l’ Experi- 0.3 practically eliminates this ence has shown that using a nonlinearity of fNL = problem and leads to satisfactory steps.

Now that a satisfactory degree of nonlinearity has been determined, the sec- to use?” has still to be answered. This is equi- ond question, “which Of the fNL valent to deciding which of the functions of interest should be used to control the step size. One approach is to control with the worst-behaved function, i.e., the one whose nonlinearity is greatest. If this approach is taken, difficulties immedi- ately appear. For example, it may be that only one function is behaving badly, the remainder being extremely linear. In this case, small steps with resulting slow convergence, or failure to converge, may result. Failure to converge here, however, does not mean divergence, but that (dU)2 is so significantly reduced due to the nonlinearity of the control function that the changes in constraints and pay- This results in negligible or random off functions become lost in computer noise.

variation of the constraint and payoff variables from iteration to iteration. The conclusion is made, therefore, that the control must be made with the best be- haved function, that is, the one having the most linearity.

To ensure that the remaining functions do not wander too far, various tests are designed to maintain convergence consistent with the linearity requirements These tests, though based on experience to some of the steepest-ascent method.

extent, have proven to be a successful means of ensuring convergence.

The majority vote test consists of the examination of the Majority vote test: changes in the constraint and the payoff functions based upon the results of the previous trial. Only those functions that axe more than a specified tolerance If at least the same number of functional from their desired value are examined.

changes are in the proper direction as those with adverse travel, then the test is satisfied. Otherwise, the step size, (du)2, is reduced and a new trial is com- puted. The majority vote test, in addition to the adverse Q test discussed below, is also used to determine whether or not a valid step is acceptable, prior to going into a reverse integration of the adjoint equations of motion. A valid step is performed after a limited number of trials or a successful trial and consists of the forward integration of the equations of motion together with the partial deri- vative calculations.

Step size coefficients: The step size coefficients are based on the func- tional nonlinearities as shown in figure 9. Using the parabolic approximation, the step size that causes the best-behaved function to have a nonlinearity of 0.3 is computed. If the resulting step size coefficient, k , fails to satisfy the condi- tion 0.5 < k < 2.0, then a further trial at the upper or lower limit is undertaken.

Experience has indicated that these bounds are reasonable ones to impose on inter- polation or extrapolation using the parabolic approximation.

Adverse Q test: The adverse Q test follows the calculation of the step size coefficient. The test ensures that, if adverse travel occurs in the payoff function Q, it will not exceed a specified tolerance. When the adverse travel exceeds this tolerance, the problem is regarded as too nonlinear and (dU)2 is again reduced, followed by a new trial computation.

Successful completion of the above tests is necessary to perform a valid step.

Other logic, however, is incorporated to force a valid step under certain circum- stances. Such tests are discussed in the user’s manual.

Figure 10 is presented to illustrate the flow of the automatic convergence process.

Variable initial conditions. - The steepest-ascent technique determines, as a by-product of the optimization process, a set of influence coefficients that defines the perturbation of the payoff quantity per unit-step-change in the initial state vector. The problem is to determine the correct direction to perturb the initial conditions in order to increase performance on the next iteration. In addi- tion, logic must be incorporated to guarantee convergence of the initial conditions to their optimum values. It is usually desirable to allow the initial condition vector to be free only in a bounded region.

The performance predicted for an iteration, given by equation 117, provides the basis for selecting the optimum initial conditions. Consider the equation For a given problem solution, the payoff function has been optimized; that is, dQ (I’) = 0, the constraints have been met, d q = 0, and the initial conditions are at their optimum values so that 6x (to) = 0. Then, since (dU)2 # 0, the relation must hold that

- f&QI$$ I@Q = 0

026) 'QQ Now restricting this to the case where the constraints are met but the initial conditions have not yet been optimized, equation 125 can be reduced to dQ p) =[(& (to) - I;Q I;$ $0 (to))]~~ (to) if it is assumed (as is the case for the optimum) that equation 126 holds.

This expression allows the variation of the free initial conditions for each iteration but is exact only for the optimum path. The correct direction to perturb x (to) is the direction to improve performance. In obtaining optimum perform- ance, perturbation of each free initial condition is considered independently since equation 127 shows that the effects of the perturbations are uncoupled.

Convergence is guaranteed by bounding the perturbations. The perturbation size is controlled using the rule that if the perturbation is in the same direction Calculate G(t) Based on New (dlJ12 I

I

W ‘)

Calculate Step Size Coefficient

I I

NO Adverse < > Phi OK Yes JL Valid Step

I I

AUTOMATIC CONVERGENCE FLOW

Figure IO

for two successive iterations, then the value of the perturbation is doubled, or halved if the perturbation changes sign. The maximum initial condition change per iteration and the limiting values are specified for each free variable.

- A cursory glance at the steepest-ascent method Weighting functions matrix.

might lead to the conclusion that the choice of weighting matrices is not a signifi- cant problem, but in practice this is not always so. The control variables, for example, may have widely differing abilities to affect the final values of the pay- off and constraint functions. It is necessary, therefore, to differentiate between powerful and weak control variables, based on their ability to affect the optimi- zation functions, particularly in their ability to affect the payoff function. Even with a single control variable, this problem remains, since along some subarcs the control may be weak and on others strong.

Two types of weighting are implied, one resulting from the differences in sensitivity of the multiple control variables and the other resulting from differ- ences in sensitivity of a single control variable along a trajectory. There are two basic reasons for the use of weighting matrices. One is to speed up conver- It is fairly obvious that a well- gence and the other is to ensure convergence.

chosen weighting matrix can improve convergence even in a well-behaved problem that involves both a weak and a strong control variable. The steepest-ascent process can be speeded up by forcing larger perturbations in the weak control variable in the initial iterations, rather than waiting for the stronger one to con- verge before making significant changes in the weaker one. What is not so obvious is the fact that, without a reasonably chosen weighting matrix, it is pos- sible that the steepest-ascent method may fail to converge entirely. Generally, it is a relatively straightforward process to develop the control history of a single control variable until it lies within the neighborhood of the optimal solution, pro- vided the problem does not involve weak and powerful subarcs.

However, the control history can only be developed into the neighborhood of the optimum since there is a region about the optimal history in which the con- trol variable is ineffective by definition. If the problem involves terminal con- straints, however, it is only ineffective for those perturbations that leave the terminal constraints unchanged. For a perturbation that merely seeks to gain performance, it may well, and usually does, remain powerful. If now there are two control variables, one weak, the other strong, it may be that even though the strong variable is in the neighborhood of the optimal solution, it still dominates the perturbations that occur as it alternates about the ideal solution. In this case, the method may fail to converge the weak control variable at all, since the effect of the weaker control can be completely obscured by the noise level of the stronger variable. Further discussion on weighting functions is given in reference 3.

Control variable power: The control variable power is used as a basis for determining a weighting function that will ensure convergence with respect to the payoff function.

The change in the payoff function resulting from a perturbation in the control variable was given by equation 97

028)

For a unit impulse perturbation in the control variables &i, it was shown in a previous section that the response of the payoff function was 6 (dQ) = A;aG

029)

where again hIQo is a 1 x n matrix and G is an n x m matrix resulting in 6 (d@) as a 1 x m matrix giving the response of the payoff function to a unit impulse of each of the m control variables. The elements of 6 (dcp) are referred to as the instantaneous payoff function sensitivites and will be designated as G

030)

These quantities measure the power of a control variable with respect to the payoff function, provided no restrictions are imposed on the terminal constraint changes.

The control variable perturbations are closely related to the instantaneous payoff sensitivities. For the case with no terminal constraints, the control pertur- bations are -1 6ii=+W G’h (131) - or in terms of the instantaneous payoff sensitivities Thus the optimum (steepest-ascent) perturbation varies directly with the inverse weighting matrix and the instantaneous sensitivities. If terminal con- straints are considered in the problem, the control variable perturbations are (from equation 116 with fixed initial conditions ) _ l/2 - d$’ I;;dz,6 (W2 6ii = f W-l s; - G’ x &; I#Q c133) I +W-lG’A @I& d 3) The results above suggest an approach to the problem of false convergence.

The problem is due to small perturbations in the weak control variable. The inverse weighting matrix based on the control variable sensitivities can be used to amplify the effect of the weaker control variables. Using this type of weighting matrix, effectively the basis of optimization is changed from that perturbation having the greatest change in the payoff function to that perturbation having the greatest change in o with all control variables being equally important and must therefore be perturbed by a reasonable amount.

The possibility of failing to converge to the desired constraints is more re- mote than that of failing to converge the payoff function. The dominant control variables for the payoff function are very often the dominant control variables for the constraints and hence will continue to be perturbed until the constraints are achieved. The terminal constraints may often be met without the optimum con- trol history. The failure to meet the terminal constraints is immediately recog- nized whereas the optimum performance is verified only by starting from a con- siderably different nominal. Weighting matrices based on constraint sensitivities (G’ Ati,) and a mixed payoff and constraint sensitivity have been .discussed by Hague in reference 3 but are not considered here.

The overall control variable power may be determined by integrating the instantaneous payoff sensitivities where the constraints are again ignored The elements of this column matrix are referred to as the integrated payoff sensitivities . Both the integrated payoff sensitivities and the instantaneous pay- off sensitivities may be used as the basis for the weighting functions.

The instantaneous payoff sensitivities can be used to construct a weighting matrix to balance weak and strong subarcs for a given control variable. The integrated payoff sensitivities are used to obtain the weighting matrix to balance the power of the multiple control variables. The matrix is of the form .g SC W.. =- l+ m . . . . 035) i = 13 1, 2, J m-t-1 1 [

[ -1 11 1

sQ U i where the wii are the diagonal elements of a square matrix. Note that when all

the control variables are equally powerful, i.e., SQ =sQ

ui -1

Q -

=- 1+mS = I 036) m+l

[ 1 SQ

the identity matrix.

This weighting matrix will ensure that each control variable is perturbed to the same order. A time-varying weighting matrix based on the instantaneous payoff sensitivity can have the form CWii j-’ = [l + T$ax ] (137) u: where s is the largest value of s $ along the trajectory. This weighting umax matrix was not automated in STOP since it was felt that a satisfactory matrix could be input based on the knowledge of the problem by the user. This also allows the user to weight regions of the trajectory as his experience dictates.

Nominal Trajectory Generation The steepest-ascent method begins with a nominal flight generated by an initial choice for the control variable history. Essentially, the only requirements on the nomin.al flight path are that it meet the stopping condition and be in the neighborhood of the optimum. A reasonable flight path can sometimes be gener- ated open loop, that is, by inputing a table for the control variable history. In many cases, however, it is difficult and time consuming to determine a control table that will even meet the stopping condition. Therefore, capability is in- cluded to generate closed-loop trajectories for the nominal. The nominal guid- ance modes are shown in table I.

Table I: NOMINAL GUIDANCE OPTIONS Mode Guidance Type Control Variable 0 8 =f(t) 8 1 CL=0 2 cr=f(t) Y =f(t) 4 ^j=f(t) 8 5 h=f(V) 8 6 h =f (M) 8

7 p=0 8

8 8

y=O

9 v=o 8

10 iI= 0 $=o r) 12 i!I=O rl The basic method for determining the control variable in STOP is a table lookup of all the selected variables as functions of time, The function of the guidance mode is to override the, tabulated input value for one of the control Therefore a control table must always be input with values for each variables.

control variable. The values of the control variables in the table which will be A brief discussion of overriden may be a dummy value since it is not used.

several methods for obtaining angle of attack or throttling is given prior to dis- cussing each guidance mode. The three schemes are: (1) determine angle of @tack to produce a given y; (2) determine the angle of attack to give a required V; and (3) determine the throttling parameter to produce a given ?.

equation of motion for flight path angle is Angle of attack from known ? , -The given by sin CY+ L) cos Q -jS - (e - ?) cos YR R m vR a2 r + yy cos B (cos f3 cos y, + sin B sin y, sin G-R)

+ 2 ocos fi cos 2j.q 038)

or for simplicity T sin 01 + L) cos cp f GF 039) mV R where + sin fi sin y, sin $R) + 2 wcos fl cos +R If the angle of attack and bank angle are small so that CYE sinol and COSQ 4 1, equation 139 can be written as (T + La) CY + L,=o ;= + GF V-41) m vR For a given G, the angle of attack required is obtained by solving equation 141 for cy, - L,+) yreq -

( m VR

GF) a=

(142)

T + L, angle of attack for a required V is Angle of attack for a known G , -The based on equation 35 = T cos 01 - D - g sin yR m req

+ a2 1” cos fl (cos B sin yR - sin fl cos yR sin z,!&)

043) or for simplicity _ T cos o! - D + VF

Ti

444)

req - m where VF = - g sin yR + a2 r cos fl (cos fl sin yR - sin p cos yR sin # ) R 045) Making the assumption that cos CYGZ 1 - c~2/2 and expanding the drag in terms of angle of attack gives D = Do + DI ~11 + D2 a2 (146) where Do = zz 0 - %M)’ ’ ’ acD D1=2-

'L, ('L, = 0 - 'LM) q '

acL2

Equation 144 can be written as T - DO - D1 o!

- (; + 4 Q2 + VF c147) t req = m For a given f,,, the angle of attack can be obtained by solving the quadratic equation (:+D2)cu2+D10L+(iTreq-VF)m-p-Do)=0 (148) Equation 148 yields two values for CX. The larger value is used since it represents the usual case.

In some cases the ? required may be larger than the maximum V that can be attained. In this case the best that can be done is to fly at the angle of attack for maximum 0. This indicates that the input h - V or h - M curves exceed the capability of the vehicle through adjusting 9. The angle of attack for maximum V is obtained by setting -Dl - 2 (5 + D2) (Y o = 049) m and solving for CY, Throttling for a known c. -The equation for thrust is given in general as

pm Ae

cos 6

051)

Solving for r) gives Pa Ae T cos 6c +144 T)= 052) ELr Tvac where T is obtained by solving equation 144 @req - VF) m + D T=

V-53)

cos a Substituting into equation 152 gives .

(Vreq - VF) m + D P Ae +L cos CYcos tic 7’ (154) KT Tvac There are no approximations here since T is defined as a linear function of Q.

Guidance Modes. - Each of the nominal guidance modes, given in table I, is detailed below. Modes 3 through 10 calculate the 8 control variable from O=C!!+y.

Mode 0: This open-loop guidance table is used as it was input by the user.

Mode 1: The angle of attack is zero. The control variable 8 is obtained by

e=y

055) Mode 2: The program interprets the first control variable input in the control table as angle of attack for the nominal trajectory only. The control variable 8 is calculated as e = ainput + Y c156) Mode 3: This mode calculates 8 where the flight path angle is given as a function of time by a guidance table input, and uses the techniques of calculating the angle of attack required to fly a given 3 (equation 142). The required 3 is obtained by (Y - YI) Preq=?I- At 057) where fI is the slope of the input y with respect to time y is the state variable yI is the desired flight path angle from the input table bt is an increment of time to correct for the y error (currently bt = 1 set) Mode 4: This mode calculates 8 when the 9 is given as a function of time by a guidance table input. The angle of attack is calculated by equation 142 from the 9 in the guidance table.

Mode 5: This mode calculates 8 where the altitude is given as a function of velocity by a guidance table input. The required V is obtained from h c

058)

req = (a/dv) req The-desired h - V curve is flown by adjusting angle of attack to give the proper V. Inaccuracies in the numerical procedure may cause drift from the - V. This is compensated by adding a feedback term as desired h

(159)

(~)req = (~)input * ~ -2~ut’ The + sign is used for descent phases and the - sign is used for ascent along - V or h - M curves. The angle of attack is obtained from equation 148.

the h Mode 6: This mode calculates 8 where altitude is given as a function of Mach number by a guidance table input.

The slope dh/dM is obtained from the input data. The technique of mode 5 is used where

dh 1 dh

( 1 TiV ( )

= a Zii input input where a is the speed of sound along the flight path Mode 7: This mode calculates 8 which produces 7 = 0. Equation 142, for calculating angle of attack for a given 3, is used as in mode 4 but no input guid- ance table is required.

Mode 8: This mode calculates 8 which will give y = 0. The method described If the flight path angle is under mode 3 is used but no guidance tables are input.

not zero at the start of the phase, the feedback term causes the vehicle to pitch until y = 0.

Mode 9: This mode calculates 8 for which the flight path acceleration &) is zero. Equation 148 is used to calculate the angle of attack for which q = 0.

Mode 10: This mode calculates 8 for which the Mach number is constant.

The requirement of l% = 0 is obtained through a creq. Since V = M a which, since k = 0 for this mode, reduces to

=MEi

062)

req

The angle of attack is obtained from equation 148.

Mode 11: This mode c’alculates the throttling (q) for constant velocity. The throttling is calculated from equation 154 for 3,, = 0.

Mode 12: This mode calculates the throttling (7) for constant Mach number.

The throttling is calculated using equation 162 for Vreq and equation 154.

Use of guidance options -A nominal trajectory generation relies on the user The to input various modes that will take the vehicle to the stopping condition.

flight path is divided into various guidance phases , each controlled by one of the guidance modes described above. A phase ends upon reaching an input value of any of the 85 state, auxiliary, and control variables or time. The number of phases is limited to 12, of which only 5 can use the input tables reserved for guidance modes. The tables required for modes 0 and 2 are not included in the five reserved tables a An example statement of a guidance phase buildup is: Phase 1. Fly a prescribed h - V (mode 5) to an altitude of 30,000 feet Phase 2. Fly a 9 = -0.2 degree per second (mode 4) and decrease to a y of 0 degree Phase 3. Fly at constant velocity using throttling as the control variable (mode 11) to a range of 1000 nautical miles Phase 4. Fly a 9 = -0.2 degree per second (mode 4) decreasing to a,y of -1 degree Phase 5. Fly a prescribed h - V (mode 5) to a velocity of 500 feet per second.

The generation of nominal trajectories using modes 3 through 10 is limited to vehicles with aerodynamics given by (163)

064)

cD or = CL,, Q

065)

cL acD

(166)

CD +y-z CL2 = ‘Do ac This restriction is made because the and thrust is independent of angle of attack.

amount of computer time required to generate a nominal can be reduced by using the simpler representations given in equations 141 and 147.

Additional Options Several options are included in STOP which permit the user to optimize special problems with fewer inputs than required in the usual sense. These are referred to as the gamma tilt, maximum payload, and circular satellite options.

Each is discussed briefly below.

Gamma tilt. -The gamma tilt option is required for vertical takeoff (VTO) The tilt angle, rtilt, is input as a function of time in a table. Termi- boosters.

nation of the tilt maneuver is accomplished by a stage stopping time; consequently, the tilt maneuver must be counted as a stage. If flight path angle, y, is selected as a free initial condition, then the perturbation on y occurs at the end of tilt.

Circular satellite. -This option permits the usage of relative circular satellite velocity, V as a stopping condition or constraint. The option csp ’ eliminates the calculations by the user as given below V rD ucos fi- v cos $$ + (VCSI sin $g2y2 c167) CSR = CSI where V CSI = and Cc = 14.081718 x 1015 = desired final radius rD ----If weight is the payoff variable and also a free initial Maximum Payload.

condition, then the program will add payload and adjust the last-stage propellant A plot of payload weight at launch so that the burnout weight will be a maximum.

versus last-stage propellant illustrates the manner in which this option operates (figure 11) .

k .

-YE maximum payload maximum paylo -- -- ----

B

/ .k F Last-Stage Propellant Last-Stage Propellant

Figure II. MAXIMUM PAYLOAD OPTtON

Notice that when the fuel required for maximum payload is less than the maximum last-stage propellant, the true maximum payload is attainable. Cther- wise, the payload limit is restricted by the amount of fuel the last stage may hold.

Conclusions and Recommendations The computer program developed for the contract meets all requirements of the statement of work (ref. 1) and in addition has the capability to select the opti- mum initial conditions for a case with free initial conditions.

For multistage vehicles, the staging times must be input in the present program. It would be desirable, in a future modification, to have the optimal staging capability using the algebraic steepest-ascent method that allows staging on variables other than time. This option would permit optimization of coast and burn times to increase performance and give added flexibility to the program.

Multiple-vehicle capability is also desirable in a program of this type.

This would enable optimization of intercept problems involving fighter aircraft, rendez- vous between spacecraft and satellites, space rescue, etc.

PROGRAM USER’S MANUAL Program Assumptions and Limitations The assumptions and limitations that were made in the development of STOP are detailed in this section. The manner in which the assumptions affect the utility of the program are indicated. The effects of the limitations on the prog- ram are apparent.

-Simulation: The basic assumption for STOP is that the Assumptions.

vehicle is represented as a point mass lifting body. This precludes any effect of the vehicle short-period dynamics on the optimum flight path. The results of the program in optimizing the flight path for a vehicle having high pitch rates have not been thoroughly substantiated. For aircraft and most boosters, the pitch rates do not affect the optimal results appreciably.

Earth model: The earth is assumed to be a rotating sphere with an inverse- square gravity field. The effect of the differences between the actual and assumed earth representation on the optimum flight path are small for the type of vehicles for which the program is designed.

Aerodynamics and propulsion: The aerodynamics and propulsion data for the vehicle can be input in several ways, depending on the availability of the data or type of vehicle. The most used aerodynamic representation for aircraft is the drag polar and lift as a linear function of angle of attack.

Generation of nominal flight paths: The generation of nominal flight paths by STOP for modes 3 through 10 assumes the drag polar and the linear lift representation and that thrust is not dependent on angle of attack. This restric- tion is not considered serious because this is the most common aerodynamic and thrust representation for aircraft. A method using more general aerodynamics and engine data can be formulated but the calculations require iterative procedures and are very time consuming.

Control variables: The control variables selected for STOP are pitch angle, bank angle, throttling, and wing sweep. The pitch angle is defined as ,g = y + Q.

The bank angle is used to perform out-of-plane maneuvers and assumes coordinat- ed turns with no slideslip. Note that, for nonzero bank, 8 is not a physical angle.

Payoff, constraints, and stopping conditions : The steepest-ascent procedure as formulated for STOP allows only state variables to be used as the payoff, constraint, and stopping functions. Enroute placards are formulated as auxiliary state variables for which terminal constraints are applied.

Atmosphere: The earth atmosphere is input as a single subroutine, using a calling sequence to transmit the data. Any subroutines using the same calling sequence can be used for the atmosphere calculation.

Staging: The trajectory can be staged only on time. Weight jettisons can be performed only at staging.

Limitations. - State variables: The number of state variables is limited to 15. These may be selected from the 20 state variables defined by equations of motion and 20 state variables defined by enroute constraints. Altitude, flight path angle, and velocity must always be selected.

Control variables: The maximum number of control variables is four.

Pitch angle, bank angle, throttling, and wing sweep are incorporated.

Constraints: A maximum of 13 constraints plus a stopping condition can be imposed on a problem. Of the 13, the number of enroute constraints is limited only by the number of data tables allowed for input of maximum and minimum placards. Twenty tables are reserved for the placards, which allows a maxi- mum of 10 placards.

Stages: A maximum of 14 computer stages are permitted. The number of physical stages is not necessarily the same as the number of computer stages.

Data deck: A maximum of 1000 data cards are allowed for tables 1 through 30 plus 31 through 34 for any stage.

Equations of motion: The point mass equations of motion include the follow- ing limitations +=O if VR < 100 fps or during a tilt maneuver I I 4 = 0 if IVR\ < 100 fps or < 0.1 degrees or Igo * YRI < 0.1 degrees p-B 1 >; = 0 if I90 & /3 1 < 0.1 degrees.

Input Data Preparation The following section presents a brief listing of the purpose of each of the nine (9) card sets that make up the Input Data Deck. This is followed by detailed instructions for preparing the data deck for a problem solution. Comments and notes are included to give the user further information to aid the setup. Prior to sending the data through the computer, it should be checked for the most com- mon errors made in data preparation, as discussed under “Trouble Shooting. I’ Input Deck Summary - The Input Data Deck is composed of nine (9) card sets of which card sets 1 through 5 and 9 are required. Card sets 6, 7, and 8 may or may not be included in the Input Data Deck, depending on the problem to be solved.

Number Reference Card of Cards Pages Purpose Set 1 1 53 Title 2 1 53 Print and Storage frequencies, type of solution, number of iterations and computer stages.

3 l-23 54-63 These cards, by use of control (NC) numbers, (See also define the input options, output options, the state, auxiliary, and control variables as well figs. 12 and 13) as the use of a number of special-purpose op- tions. Included in this selection are the equa- tions of motion and the enroute placards.

64-65 4 3-5 Initial condition data including various physical constants, value of stopping condition, integra- tion error limits, stopping condition tolerance, DUSQ, and initial values of the various state variables.

65-66 5 8-16 Computer stage dependent parameters such as: initial weight, final stage time (trajectory time), aerodynamic reference area, nozzle cant angle, total impulse correction factor, weight flow correction factor, inert weight flow and mini- mum allowable integration step.

66-67 6 o-13 Constraint dependent parameters; desired values and final tolerances for the various terminal constraints and enroute placards.

o-7 Free initial conditions, inputs maximum and minimum value and initial increment per iter- ation of each free initial condition.

67-68 8 o-12 Nominal guidance data; defines the various op- tions used in the nominal trajectory generation for each stage.

Card Number Reference Pages Purpose Set of Cards All the tabulated input such as control function, 9 2- 58 68-76 weighting function, gamma tilt program, en- route placards, aerodynamic data (lift, drag, etc. ), and engine data.

Preparation - The following are detailed instructions for preparing the data deck for a problem solution: Card Set 1 Format A72 1 Card Title Card Card 1 It may contain The contents of this card will be printed at the start of the case.

any desired information in columns l-72 or left blank, but it must be the first card in the data case.

Card Set 2 1 Card Format 1415 Controls Card 1 Column Name Description IPRNTl Printout of forward trajectory occurs every IPRNTl l-5 integration steps.

6-10 IPRNTS Printout of influence coefficients occurs every IPRNTB integration steps during backward integration if NC(7) = 1.

11-15 ISTORl Storage of time, partials, and plot data occurs every (Recom- ISTORl integration steps along forward trajectory.

mended value = 1 for rockets, 2 for airbreathers.)

16-20 ISTORB Storage of impulse response functions occurs every ISTORB integration steps along backward trajectory.

(Recommended value = 1 for rockets, 2 for airbreathers. ) 21-25 NARBY Degree of problem solution indicator = 0 Integrate nominal trajectory only and quit = 1 Integrate nominal trajectory and backward trajectory only and quit = 2 Iterate (NARBY = 2 required for optimization) Column Name Description Maximum number of iterations 26-30 NITC Number of (computer) stages 31-35 MSTAGE (1s MSTAGE I: 14) Card Set 3 Format l4I5 1-23 Cards NC Controls Each card contains seven pairs of integers. The first integer is the index of the Every card NC to be set and the second integer is the desired value of the NC.

must contain seven pairs of integers, except the last card. The NC’s may be in- put in any order. The reading of the NC’s is terminated by a blank field in col- umns 61-65 of the last card. Thus, if the number of NC’s to be read is not a multiple of 7, the last card will contain one to six pairs of integers. Otherwise, The explanation of the NC’s is as follows. NC’s not the last card will be blank.

read in will be zero.

Description NC Index Name 1 Plot paper option NC (1) = 0 use 11 x 17 millimeter paper = 1 use 8-l/2 x 11 millimeter paper 2 Not used NC (2) Plot identification option NC(3) ORTHOMAT (NC(ll) < 0 ) = 0 Use symbols = 1 Use colors SC 4020 (NC(ll) > 0 ) = 0 9 x 9 vellum only ZZ 1 microfilm only = 2 9 x 9 vellum only = 3 9 x 9 vellum and microfilm plots 4 Number of constraints + 2 NC(4) (Count both terminal constraints and enroute placards) (2s NC(4) S 15) Description NC Index Name 5 Number of free initial conditions NC(5) OS NC(5)< 7 Note: If a tilt maneuver is selected (NC (14) = l), and if gamma is a free initial condition, then the free initial condition will be applied at the end of tilt (start of com- puter stage 2) 6 Maximum payload option NC (6) = 0 Ignore option = 1 Last-stage propellant will be limited to the maximum last-stage propellant (MLSP, input in CARD SET 4).

Weight must be a free initial condition for this option and NC(19) = 1 is required.

Note: If NC(6) = 0 and NC(19) = 1, final weight is maximized by minimizing last-stage fuel.

Influence coefficient print option NC(7) = 0 Do not print influence coefficients = 1 Print influence coefficients (Usual value = 0) 8 Punch option NC(g) = 0 Do not punch control and restart tables used for last iteration = 1 Punch control and restart tables used for last iteration (Usual value = 1) 1) If a guidance mode is selected (NC(13) > 0) Notes : for a nominal trajectory (NARBY=O), then if NC (8)=1, STOP will punch a control table (table 0), which may be used to generate the same trajectory open loop (NC(13)=0).

The control table punched by STOP is in 2) octal (Format 6 8 12).

Extremal option NC (9) = 1 maximize payoff variable = -1 minimize payoff variable = 0 solve boundary value problem NC Index Name Description 10 Circular satellite option NC(l0) = 0 Ignore option = -1 Stop on V at desired final altitude. (NC(18) = 4 is required or this option ES = N Nth constraint is V = Vcs This requires NC(N+19)=4 11 Plot option NC(ll) = 0 Do no plot = N Plot first, last, and every Nth iteration on SC 4020 plotter ZZ -N plot first, last, and every Nth iteration on ORTHO- MAT plotter Note: Only stored points are plotted.

Numerical partial check option 12 NC(12) = 0 Do not check partials = N Check partials every Nth stored integration step and printout analytical and numerical values (Usually NC(12) = 0) Nominal trajectory guidance mode option 13 NC(13) = 0 Open loop (control history from table 0) = N Closed loop, N phases (see card set 8 if NC(13) > 0) (Usually NC(13) = 0), l< N 512 Tilt maneuver option 14 NC(14) = 0 Do not use tilt maneuver = 1 Use tilt maneuver for stage 1 (table 2) (NC(14) = 1 for VT0 systems) and the tilt maneuver counts as a stage Skin friction option 15 NC(15) = 0 NoCDF = 1 Href, a CD/a H = f (M). This requires that table 23 be input Control table input option 16 NC(16) = 0 input decimal table = 1 input octal table Note: NC(16) = 1 is required to accept tables punched by the program NC Index Name Description 17 NC(17) Restart table option = 0 No restart table (use for nominal trajectory) = 1 Input restart table after table 0 Note: A restart table is obtained whenever a control variable table is generated by the program as a result of previous iterations. See punched card output discussion; 18 NC(18) State variable index for stopping condition (figure 12, See input instructions check list.

column 2).

19 NC(19) State variable index for payoff quantity (figure 12, column 2) 20 NC (20) State variable index for first constraint (figure 12, column 2) 21 NC (21) State variable index for second constraint (figure 12, column 2) 1 I I 1 1 1 I t 32 NC (32) State variable index for 13th constraints (figure 12, column 2) Note: Enroute placards are state variables and count as See figure 12 for list of state variables.

constraints.

Engine option for stage 1 33 NC(33) 34 NC(34) Engine option for stage 2 r r 1 1 46 NC (46) Engine option for stage 14 Note: The engine option for each stage is input as a two- digit number. The left digit specifies the type of thrust table as follows: 'ORTRAN jYMBOL DESCRIPTION VARIABLE VARIABLE NC PLOT JNITS NAME TYPE INDEX INDEX INDEX STATE 1 121 1 K(KU N WEIGHT LB (EQUATION 2 122 2 K(K2) i ALTITUDE FT 3 123 3 X(K3) OF-I ITION) RELATIVE FLIGHT PATH ANGLE Y DEG 124 4 X(K4) RELATIVE VELOCITY FPS Ii 5 125 5 X(K5) LATITUDE ANGLE DEG 6 126 6 X(K6) RELATIVE HEADING ANGLE DEG tR 127 7 XiK7) h LONGITUDEANGLE DEG B 128 8 XIK8) TD DUMMY TIME SEC 9 129 9 X(K9) PATH RANGEALONG EARTH'S SURFACE RNG N Ml l3Q 10 X(K10) 4HI AERODYNAMIC HEATING INTEGRAL FT-LBIFT2 11 11 X(K11) IDEAL RELATIVE DELTA VELOCITY AV FPS 12 z: X(K12) SL GRAVITY LOSS FPS l3 I33 :: X(K13) DL DRAG LOSS FPS 14 134 14 X(K14) rvL THRUST VECTORINGLOSS FPS 15 135 15 X(K15) ER RELATIVESPECIFIC ENERGY FT 16 136 16 X(K16) 17 137 17 X(K17) AVAILABLEFOR EXPANSION 138 18 X(K18) NOT DEFINED 19 139 19 X(K19) I 20 140 20 X(K20) 21 141 21 X(K2U STA PITCH ANGLEPLACARD= F(t) DEG2SEC 22 142 X(K22) 22 (ENROUTE BANK ANGLEPLACARD = F(t) DEG'SEC 143 23 X(K23) RDI THROTTLE PLACARD * F&(M) SEC 24 144 24 X(K24) WING SWEEP PLACARD = F(H,M) DEG'SEC 145 25 X(K25) ANGLE OF ATrACK PLACARD = F(H,M! DEG'SEC 26 146 26 X(K26) NOT DEFINED 27 147 27 X(K27) HD' ALTITUDERATE PLACARD= F(t) FPS 28 148 X(K28) iPSF1' SEC 28 Q" DYNAMIC PRESSUREPLACARD= F(t) X(K29) DYNAMIC PRESSUREPLACARD = F(M) 29 149 29 (PSF)*SEC 30 150 30 X(K30) El* Q'IPLACARD = F(M) (PSF-DEGI*SE 31 X(K3l.l TEMT' TOTAL TEMPERATURE PLACARD = F(t) I"R?iEC 31 I51 NORMAL LOAD FACTOR PLACARD = F(H.M) 32 152 32 X(K32) N+ SEC 33 153 33 X(K33) RPA' RESULTANT PHYSIOLOGICALACCEL. PLACARD = F(t) (FT/SEC212SEC 34 34 H' ALTITUDE PLACARD = F(M) w X(K34) FT2SEC A P* SONIC BOOMOVERPRESSURE PLACARD = F(A,BI 35 I.55 35 X(K35) (PSF+SEC 36 156 36 X(106) M' MACHNUMBER PLACARD = F(H) SEC 37 X(K37-l 157 37 NOT DEFINED 38 158 38 XIK38) SAMMA PLACARD = F(H,MI DEG'SEC 39 39 X(K39) NOT DEFINED 40 X(K40) NOTDEFINED

FIGURE12: STATEVARIABLES

Thrust Option Thrust Table Digit -- 0 Coast (no thrust table) 1 Vacuum thrust (lb), nozzle exit area (in2), vacuum specific impulse (set) = F (stage time). Table format 1 for 3 depen- dent variables 2 Vacuum thrust (lb), nozzle exit area (in2), vacuum specific impulse (set) = F (total time). Table format 1 for 3 depen- dent variables T = F (ar, H, V) T = F (a, H, M) T = F (CY, M, II) Table format 2 for one, two, or T/P.= F (ar, I-I, M) three independent variables CT = F ((u, H, V) CT = F (01, II, M) CT = F (a, M, H) Note: (1) P is atmospheric pressure for option 6 ; (2) For thrust options 3 through 9, the number of engines is input in columns 11 to 20 on the first data card in table 31. If this field is left blank, STOP will assume one engine. The thrust from table 31 (and fuel flow from table 32 if a WDOT option is selected) will then be multiplied by the number of engines.

The right digit specifies the type of fuel flow table as follows: Fuel Flow Option Digit Fuel Flow Table Coast (no fuel flow table) WDOT = F (H, M, q) Table format 2 for one, two, or 2 SFC = F (q, H, V) three independent variables SFC = F (a, H, V) Fuel Flow Table (Continued) Digit 4 WDOT = F (?j, H, V) WDOT/P = F (H, M, ?j) 6 SFC = F (Q, M, H) Table format 2 for one, two, or three independent variables 7 SFC = F ((u, M, H) 8 WDOT = F (7, M, H) 9 WDOT = F (a!, M, H) Note: P is atmospheric pressure for option 5 Description NC Index Name Aerodynamics option for stage 1 47 NC(47) Aerodynamics option for stage 2 48 NC(48)

I

I 1

Aerodynamics option for stage 14 6b NC{60) Note: The aerodynamics options available are given below and are identified by the value assigned to NC(47) - NC(60).

Aerodynamic Option Aerodynamic Table NC Value 0 Vacuum (no aerodynamic table) 1 CD = F (H, M, Q ) CL = F (H, M, a) 2 CD = F (A, M, (u) CL = F (A, M,ar) WD/ aCL2, 6CL/dOr = F (M) CD,, CA, aCN/aar, aCN/dCY3 = F (M) CLcrco= F (A, M) CDo, dCD/dCL2, CLmin 5 acL/d~!, drag Note: (1) Options 1 and 2 use table format 2 for one, two, or three independent variables ; (2) Options 1 and 2 require two tables -the first is a CD table; the second, a CL table; (3) Options 3 and 4 are table format 1 for three dependent variables; (4) Option 5 uses table format 3 for two independent and five dependent variables; (5) CD, in option 5 is CDmfn; (6) All angles and slopes are in degrees.

Description Name NC Index 61 Index of dependent variable for 1st plot NC(61) Index of independent variable for 1st plot 62 NC (62) Index of dependent variable for 2nd plot 63 NC (63) 64 NC (64) Index of independent variable for 2nd plot I I Index of dependent variable for 7th plot 73 NCi73) 74 Index of independent variable for 7th plot NC(74) Note: See figures 12 and 13 for plot index 75 Sonic boom placard option (Table 30 must be input) NC(75) = 0 Planar sonic boom = N Area sonic boom N is number of points to be used along the shock to deter- mine the placard violation. N must be odd and greater (Recommended value is 7. ) than one.

76 Not used NC(76) 77 NC(77) Pitch angle option = 0 Do not use pitch angle = 1 Use pitch angle 78 NC(78) Bank angle option = 0 Do not use bank angle = 1 Use bank angle 79 NC(79) Engine throttling option = 0 Do not use throttling = 1 Use throttling 80 Wing sweep option NC(80) = 0 Do not use wing sweep = 1 Use sweep AK (1) print option 81 NC(81) = 0 Do not print AK (1) = 1 Print AK (1) AK (2) print option 82 NC (82) = 0 Do not print AK (2) = 1 Print AK (2) I 1 I 1

1c PLOT ;YMBOL DESCRII'TION JNITS

VARIABLE VARIABLE FORTRAN NDEX TYPE INDEX INDEX NAME THRUST .B 4UXILIARY I ;I 11 AK(l) , 'RINTOUT 2 12 AK(2) LlFl .B I 13 AK(3) ; DRAG B DYNAMIC PRESSURE 'SF I Ii I4 AK(4) 3 - i I5 15 AK(5) H MACH NUMBER 5 I6 16 AK(6) XR GREATCIRCLE RANGE '4 Ml I I7 17 AK(7) NDOT WEIGHTFLOW -B/SEC :Ps 3 I8 I8 AK181 YI INERTIAL VELOCITY 3 19 19 AK(9l INERTIAL FLIGHT PATH ANGLE IEG 10 K1 50 AKilDl ;I INERTIAL HEADING ANGLE IEG ANGLEOF ATTACK IEG.

11 51 AK(11) AK1121 'R 12 ii 52 F;MP TOTALTIMPERATURE l3 53 AK(U) VA SPEEDOF SOUND 'PS 14 54 AK(14) P AMBIENT PRESSURE %F s P 15 35 55 AK(l5) AMBIENT DENSITY SLUGIFT' 16 96 56 AK061 TEM AMBIENT TEMPERATURE 'R 17 31 57 AK(17) VEWTONS WM WEIGHT(METRIC) METER 18 98 58 AK(18) HM ALTITUDE (METRIC) RESULTANT PHYSIOLOGICALACCELERATION FTlSEC 19 w 59 AKU9) RPA - 20 60 AK(20) NORMAL LOAD FACTOR 61 AK121) DYNAMIC PRESSURE TIMES ANGLEOF Al-TACK PSF-DEG 21 101 !a AK(22) c LOCAL GRAVITY FTnSEc2 22 1CQ RADIUS VECTOR Ff 23 103 ii AK1231 PSF 24 104 64 AK(24) SONIC BOOMOVERPRESURE Lb PSF AK(25) LIFTCOEFFICIENT 25 105 65 CL DRAGCOEFFICIENT 26 106 66 AK(26) CD LIFT/DRAGRATIO 27 107 67 AK(271 UD SPECIFIC FUELCONSUMPTION 68 AK(28) SFC 28 108 l/SEC AK1291 VM 29 109 69 RELATIVE VELOCITY(METRIC) M/SEC TM 30 110 70 AK(301 THRUST(METRIC) NEWTONS 31 111 71 AK(31) LM NEWTONS LIFT (METRIC) 32 112 72 AKUP) DRAG (METRIC) NEWTON T 73 AK(33) :M" 33 113 DYNbMIC PRESSURE(METRICI NWTSIM WDOTM 34 114 74 AK041 WEIGHTFLOW(METRIC) KGFlSEC 35 115 75 AK(35) VIM INERTIALVELOCITY(METRIC) M/SEC 76 AK06) 36 116 QaM DYNbMIC PRESSURETIMES ANGLEOFAllACKtMETRIC NWT-RAD/M* AK(371 37 117 77 iPM SONIC BOCMOVERPRESSURE (METRIC) NWISIM2 38 118 78 AK(38) SFCM SPECIFIC FUELCONSUMPTION (METRIC) KFGINWT-SEC 39 119 19 AK(39) NOT DEFINED > 40 120 80 AK1401 U(Il) PITCH ANGLE DEG CONT*ROL 1 17 81 U(l2) DEG 2 78 BANK ANGLE - 3 79 ii lJ(l3) ENGINETHROTTLE 4 80 84 U114) WING SWEEPANGLE DEG T TIME SEC INDEPENDEF 85

FIGURE 13 AUXILIARY, CONTROLAND INDEPENDENTVARIABLES

NC Index Name Description 120 AK (40) print option NC (120) = 0 Do not print AK (40) = 1 Print AK (40) 121 X (1) state variable option NC (121) = 0 Do not integrate X (1) = 1 Integrate X (1) NC (122) X (2) state variable option = 0 Do not integrate X (2) = 1 Integrate X (2) X (20) state variable option 140 NC (140) = 0 Do not integrate X (20) = 1 Integrate X (20) 141 NC(141) X (21) state variable option = - 1’ Observe minimum placard = 0 Ignore placard = 1 Observe maximum placard = 2 Observe both maximum and minimum placards 142 NC (142) X (22) state variable option = -1 Observe minimum placard = 0 Ignore placard = 1 Observe maximum placard = 2 Observe both maximum and minimum placards X (40) state variable option 160 NC(160) = -1 Observe minimum placard = 0 Ignore placard = 1 Observe maximum placard = 2 Observe both maximum and minimum placards Note: (1) kigures 12 and 13 list state variables and AK’s; (2) X (2), X (3), and X (4) must always be selected; (3) Two tables are reserved for each placard in the order selected, beginning with table 3; (4) Two tables must be input for a placard only if table format 2 is required and both max and min are $0 be observed.

Otherwise, one table is required; (5) For table format 2, the odd numbered table is reserved for the min placard, the even num- bered table for the max placard; (6) For table format 1, both limits are always input in the odd numbered table with min plac- ard first; (7) A maximum of 15 state variables may be used.

Card Set 4 3-5 Cards Format 7FlO Initial Condition Data Card 1 Column Name Description l- 10 XSTP Desired value of stopping condition (checked in last stage only) Earth gravitational constant times lo-l5 (Recommended 11-20 GR value = 14.081718) 21-30 RO Earth radius (feet) (Recommended value = 20,902,992) 31-40 OMEGA Earth rotational rate (rad/sec) times lo4 (Recommended value = 0.72921152) 41-50 GO Sea level gravity (ft/sec2) (Recommended value = 32.174) Card 2 T (may be left blank) l-10 DUSQ Initial value of 6ii’W 6ii dt Forward trajectory integration error limit for variable 11-20 FINNER step integration (Recommended value = 18 for rockets, 15 for airbreathers) Backward trajectory integration error limit for variable 21-30 BINNER step integration (Recommended value = 18 for rockets, 15 for airbreathers) Note: FINNER and BINNER are in some respect the number of binary bits of accuracy to be held during integration.

Eighteen yields four to five significant figures of accuracy.

Changing FINNER and BINNER by 3 amounted to changing accuracy by about 1 significant figure. Decreasing FINNER and BINNER is the direction of decreasing accuracy and decreasing computer time 31-40 EPSLN Maximum allowable tolerance in stopping condition.

(Usually EPSLN should be adjusted to give six to eight For example, if XSTP = significant digits for XSTP.

25,000, EPSLN should be about 0.01. If XSTP = 0, EPSLN should be about 0.000001) :ard 2 :01umn Name Description -- - X-50 MLSP Maximum last-stage propellant (A value required if NC(6) = 1) 11-60 Trajectory starting time (Usual value = 0) TSTART il- 70 L Reference length for sonic boom calculations lard 3 Description Zolumn Name i-10 Initial value of first state variable selected from figure 12 x0 (1) i l-20 Initial value of second state variable selected from figure 12 x0 (2) , f 1 , , , 31- 70 Initial value of seventh state variable selected from figure 12 x0 (7) (Additional cards are required if more than seven state variables are selected. ) Note: These values are input in ascending variable index order (as selected from figure 12) with no blank fields.

Card Set 5 Format 7FlO 8-16 Cards Stage-Dependent Parameters Each card subset contains the stage-dependent parameter values in fields of 10; Two i.e., columns 1 to 10, value for stage 1; 11 to 20, value for stage 2, etc.

cards are required in each subset if MSTAGE > 7.

Card Subset Parameters --~ Name 1 Initial weight (lb) of each stage WO (NSTAGE) Note: (1) A positive value will be used as an initial weight; (2) A zero value will use the end weight of the previous stage at the start weight of the stage; (3) A negative value will use the end weight of the previous stage minus the input value for the stage.

Stopping trajectory time for each stage 2 TSTP (NSTAGE) Note: The last stage does not stop on time (except dummy time if it is a state variable) unless the trajec- tory exceeds the input stopping time for the last stage, in which case the run will be aborted.

Card Set 5 Parameters Card Subset Name Aerodynamic reference area (ft2) for each stage SRE F (NSTAGE) DE LC (NSTAGE) Nozzle cant angle (deg) for each stage measured from body centerline XKT (NSTAGE) Total impulse correction factor (XKT) for each stage Note: These values may be used to correct total impulse or each stage since m ‘BO I (XKT) (TVAC) dt s TOT = TSTART (Usually use XKT’s = 1) 6 XKM (NSTAGE) Weight flow correction factor (XKM) for each stage Note: These values may be used to correct the final weight of each stage since W (XKM) (WDOT) S’“” [ = Winitial - final TSTART - WDoTinerts dt I (Usually use XKM’s = 1) 7 WDOTI (NSTAGE) Inert weight flow (lb/set) for each stage 8 DTMIN (NSTAGE) Minimum allowable integration step for each stage (Recommended Value = 0.1 second) Card Set 6 Constraint-Dependent Parameters Format 7FlO O-13 Cards Parameters Card 1 Desired value for first terminal constraint l-10 PSI 11-20 Desired final tolerance for first terminal constraint DPSI Card 2 l-10 Desired value for second terminal constraint PSI 11-20 Desired final tolerance for second terminal constraint DPSI Note: (1) Enroute placards are treated as terminal constraints and require a desired terminal value of zero; (2) The desired - Card Set 6 (Continued) value of the constraints and the tolerances must be ordered as selected by the NC(20) - NC(32) controls. For instance, if NC(20) = 5, then the first terminal constraint will be lati- tude as indicated in figure 12; (3) One card required for each constraint.

Card Set 7 -Free Initial Conditions Format 110, 3FlO O-7 Cards One card is input for each desired variable initial condition as follows: Description Column Name l-l.0 IX Index of state variable with free initial condition (See figure 12) (Right adjusted) 11-20 XOIMAX Maximum allowable value 21-30 XOIMIN Minimum allowable value 31-40 DELX Initial increment per iteration Note: (1) This card set will be missing if NC(5) = 0; (2) If YTILT is a free initial condition, then the maximum allow- able value must be 200.

Card Set 8 Format 215, 6FlO Nominal Guidance Data O-12 Cards Description Column Name l-5 MODE Guidance mode index 6-10 NSW Index of variable used to switch guidance modes (see figures 12 and 13, plot index). NSW is input positive if SWVAL is to be approached from below, or negative if it is to be approached from above 11-20 SWVAL Value of the variable on which guidance mode is switched 21-30 AMIN Minimum allowable angle of attack 31-40 AMAX Maximum allowable angle of attack 41- 50 TMIN Minimum allowable throttling 51-60 TMAX Maximum allowable throttling Note: (1) The guidance modes available (MODE) are as follows : 6’7 Card Set 8 (Continued) Guidance Mode Options Mode Descriptions Table Format 0 f(t), open loop e= 1 cY= 0, calculate 8 01= 2 f(t), calculate 8 3 Y= f(t), calculate 8 jb f(t), calculate 8 5 h= f(V), calculate 8 6 h= f(M), calculate 8 i/Z 7 0, calculate 8 Y= 8 0, calculate 8 += 9 0, calculate 8 & 10 0, calculate 8 += 11 0, calculate 7 R;I= 12 0, calculate rl (2) The maximum number of phases is limited to 12, of which only five may be phases requiring tables; (3) Modes 0 and 2 use the input control variable table (table 0) and therefore are not included in the five reserved guidance tables (Note 2); (4) TMIN and TMAX are input only for MODES 11 and 12; (5) If SWVAL occurs before XSTP on the last phase, the re- mainder of the trajectory will be flown open loop using table 0; (6) The trajectory will always stop on XSTP; (7) The guidance mode option is for the nominal trajectory only. If NARBY = 2, STOP will iterate from the nominal in the usual open loop mode; (8) This card set will be missing if NC(13) = 0; (9) Aerodynamic option 3 or 5 must be used when modes 3 through 10 are used; (10) For modes 5 and 6, MODE is input as a positive number for ascent or negative for descent. In either case, the H-M or H-V table must be monotonically increasing in the independent variable.

Card Set 9 Tables Format 215, A60/ (7FlO) 2- 58 Tables Every table must begin with a title card as follows: Column Description 1-5 Table number 6-10 Number of data cards in the table (do not count this card) ll- 70 Any desired information Card Set 9 (Continued) Input Data Tables Table Number Description 0 Control function table (U) Table format 1, NU dependent variables vs time where NU is the number of control variables selected. This table must always be present as the first table and will not be extrapolated.

Note : Restart table follows table 0 if NC(17) = 1 and is also assigned as table 0.

-1 1 Weighting function table (W diagonal elements). Table format 1, NU, dependent variables vs time. ’ This table must always be pres- ent as the second table and allows the user to aid in assigning the value of DUSQ along the trajectory and among the control variables.

A zero value will cause no change in a U, and a 1. will allow the full change determined by the gradient. As an example, it is a good idea to restrict pitch angle changes during a tilt maneuver and For instance, during the low-speed portion of a VT0 trajectory.

this table might be as follows for a VT0 liquid rocket tilting at 8 seconds after launch ( 8 control): 0. 0. 15. 0. 30. 1. 10000.

1.

2 Gamma tilt table Table format 1, 1 dependent variable vs time (Used only if NC(14) = 1) First table for 1st enroute placard (Format specified by card set 3 and figure 12) Second table for 1st enroute placard .

.

.

21 First table for 10th enroute placard Second table for 10th enroute placard ACD 23 l)H = F(M) Skin friction drag table (required if NC(15) = REF’ H-HmF Table format 1, 2 dependent variables vs Mach number First guidance table (Format specified by card set 8) Second guidance table : .

Table Number Description 28 Last guidance table Lateral cutoff distance (required if NC (75) >l), yC0 = F (M, H) Table Format 2 Sonic boom table (required if NC(104) or NC(135) = 1 APJ = F (H, M, DENOM) Table Format‘2 where : l/4 H 3/4 APJ = (AP) y iKRfi pref .

AP = Overpressure - psf = Sonic boom constant z 1.9 KR

B

DENOM = 2 CL+ 31 1st engine table for stage 1 Note: Engine options 10 and 20 require only one table 32 2nd engine table for stage 1 1st aerodynamics table for stage 1 Note: Aero options 3, 4, and 5 require only one table 2nd aerodynamics table for stage 1 31 1st engine table for stage 2 32 2nd engine table for stage 2 33 1st aerodynamics table for stage 2 34 2nd aerodynamics table for stage 2 (Repeat tables 31 through 34 for each stage) Note: (1) Some or all of tables 2 through 34 may be missing, depend- ing on the options selected from card set 3; (2) The first-stage.depend- ent table (31-34) of each stage must be preceded by a card with the word STG in columns l-3 (remaining columns must be blank); (3).The last- stage dependent table for the last stage must be followed by a card with the word END in columns l-3 (remaining columns must be blank); (4) The STG cards are required even for stages with no tables; (5) All tables except table 0 may be extrapolated; (6) The total number of data cards for tables 1 through 30 + 31 through 34 for any .&age is limited to 1000.

Card Set 9 (Continued) Three formats are available for data table input. The table lookup routines use linear interpolation and will extrapolate, although the extrapolation feature should be avoided if possible.

by 1, 2, 3, 4, or 5 Dependent Table Format 1 1 Independept X = independent variable Yi = dependent variables NDV = number of dependent variables, x, Y1, Y2, --- Y NDV’ x, Y1’ Y2’ --- YNDV’ x, --- 7 entries per card Figure 14 shows an example for two dependent variables Table Format 2 1, 2, or 3 Independent by 1 Dependent Cardset 1 -- Column Name Description l-10 NPL Number of planes I @nay = 1) 11-20 ENG Number of engines (used only for thrust options 3-9) Cardset 2 (All data for first plane) -- Card Subset 1 Description Column Name l-10 Xl Value of the first plane (Xl = 0, if only one plane) 11-20 NCV Number of curves on first plane pay = 1) Card Subset 2 (All data for first curve) Card 1 Column Name Description l-10 x2 Value of the first curve (X2 = 0, if only one curve) 11-20 NPT Number of data points on first curve (must be > 1) 21-30 X3 Value of the 1st point 31-40 Y Value of the dependent variable 41-50 x3 Value of the 2nd point 51-60 Y Value of the dependent variable Value of the 3rd point 61-70 X3 Card Set 9 (Continued) Card Subset 2 Card 2 Value of the dependent variable l-10 Y Value of the 4th point 11-20 x3 I 1 !

1 I 7 values per card until all points for the first curve are input Card Subset 3 (All data for the 2nd curve) Same format as card subset 2.

Card subsets continue until all curves on the first plane are input.

The next plane is input in the same manner, beginning with card set 2.

Note: (1) Each plane begins with a new card; -(2) Each curve begins with a new card; (3) All data in the table is input floating point (Every value must have a decimal point); (4) The table lookup routine uses linear interpolation and will extrapolate linearly in all directions.

Figure 15 shows an example for Table Format 2.

Table Format 3 2 Independent by 5 Dependent This format is used only for the aerodynamics option 5 as follows: Cardset 1 column Name Description l-10 Number of constant sweep curves (may be 1) NSWP Cardset 2 (All data for first sweep) Card Subset 1 Column Name Description l-10 Sl Value of first sweep curve 11-20 NM1 Number of Mach number points on first curve Card Subset 2 Card 1 Column Description Name l-10 Ml Value of first Mach number 11-20 CD0 C ardset 9 (Continued) Card Subset 2 Card 1 Description Column Name 21-30 BCD/BCL2 31-40 CL min 41- 50 6C L/8 OL 51-60 CL0 61-70 M2 Value of second Mach number Card 2 Description Column Name -- ~ l-10 CD0 11-20 6 CD/BCL2 (7 values per card until all data for first sweep or input) Repeat cardset 2 for each sweep angle.

Figure 16 shows an example for Table Format 3.

I I I I I I I I I I

t t

OV OV / /

0 0 .2 .2 .4 .4 .6 .6 .8 .8 MACH NUMBER MACH NUMBER 1111111111111111111111lll,l,lll,ll,l,,l,, NOTE: ADDITIONAL POINTS INCLUDED TO PREVENT TABLE LOG KUP EXTRAPOLATION

Figure14 EXAMPLEOFTABLE FORMAT1

-1 .o t ANGLE OF ATTACK ANGLE OF ATTACK 0. 4. !-510. I I- .!3 -110. -. 3 .3 .

.3 ! I I I I I I I I I I I I .s 6. -5 0. ;d. 5 -1’0. -. 5 Ilj-5!.(1IlIIII T I rrt I I I I r I I I i 5.’ .5 ‘. 3 .3 50.

I I I-JI I I I I I I I I 5.

.6 6. L5 0,. 7 -1 0,. -. 7 , , ,-, , , , , , , , , -. 5, 5. .5 ;o’ * .I 50 . 1 1 I.71 1 1 1 1 1 1 1 1 I 10 . I 2. -5ln n ml n. I I I I I I I Ii I II 1 1 1 1 r-li0.i i i i- . is i i i C ’ i i i i i i i i 51.

I I AI I I I hl I I I I Al I I I I I I I I NOTE: ADDITIONAL POINTS INCLUDED TO PREVENT TABLE LOOKUP EXTRAPOLATION

Figure 15 EXAMPLEOFTABLEFORMAT2

A-0 h-200 2- cDo b

I I

I

a 1

I A

.5 I .o 0 .5 1.0 MACH NUMBER MACH NUMBER

Figure 16 EXAMPLEOFTABLE FORMAT3

‘76 PROGRAM OPERATION Nomenclature - FORTRAN COMMON. - The table that follows is a list of the FORTRAN names, mathematical symbols, units, descriptions, and sub- routines where computed for all items in COMMON.

FORTRAN Subroutine Name Symbol Units Description Computed set Maximum step size that the AHS STEP1 RKVS integration package STEP2 (STEP1 and STEPB) will take for its next integration step.

(This may be decreased, but never increased for the step to be taken) AK( 40) Auxiliary printout variables AKSTP (see appendix B) STPl cl AK11 Angle of attack AKSTP deg AMAG(15) Array used by RKVS for the INITAL integration of the forward trajectory AP(40) A packed array containing only STPl those auxiliary variables (AK’s) to be printed at each time point on forward trajectories.

BINNER Input backward trajectory inte- INITAL gration error limit for Runge- Kutta variable step integration CAKll cos Q Cosine of the angle of attack AKSTP CD Drag coefficient TLD cD Cl1 cos 8 Cosine of the pitch angle INITAL (Set to 1. if NC(77) = 0) C@NTR c12 Cosine of the bank angle cos q INITAL (Set to 1. if NC(78) = 0) C@NTR CK3 Cosine of the flight path angle cos YR AKSTP CK5 cos J3 Cosine of the latitude (Set to INITAL 1. if NC(125) = 0) AKSTP CK6 Cosine of the heading angle INITAL cos +R (Set to 1. if NC(126) = 0) AKSTP Lift coefficient TLD CL CL C@EFK Step-size coefficient KCALC An array containing aerody- TLD C@F(2,5) coefficients (see Appendix B) IMTCG Array of constraint(s) and INITAL C@N(W performance name(s) FORTRAN Subroutines Name Symbol units Description Computed Array of program constants (see INITAL CaNST(30) appendix B) BLOCK Nondimensional constrain CPsI(13) INITAL change(s) asked for VALID DEG Radians to de-grees conversion BLOCK factor DELC(14) 6 Input nozzle cant angle for each INITA L deg C stage measured from body center- line. Assumes zero thrust compo- nent normal to body centerline.

SJ Increment added to u vector for DE LU(4) UCALC next trajectory Array of incremental changes in INITAL DELX( 7’) Go free initial conditions VARIC Input maximum allowable incre- INITAL DELXI( 7) ment by which a free initial con- dition can change per iteration.

Array of incremental changes in INITAL DELXS( 7) variable initial conditions in order VARIC of selected free initial conditions DENaiM Denominator of first Lagrange MATRIX multiplier An array of alphanumeric DESXUK(170) BLOCK characters containing all of the state, auxiliary (AK), and control variable names that can be selected for printout via the NC array. Each name is defined by two sequential words in this array.

A packed array of alphanumeric TITLES DESXUP(120) characters containing the printout headings for a given computer run.

DPHI Change in performance (with MATRIX respect to last accepted iteration) MATRIX DPHID Predicted change in performance due to DUSQ DPHIP Total predicted change in MATRIX performance FORTRAN Subroutines Name Symbol Units Description Computed MATRIX Predicted changed in per- DPH13 formance due to variable initial condition(s) .

DPSIM( 13) Input desired final tolerance(s) INITAL for terminal constraints.

Value of DUSQ at the start of MATRIX DPBBAR the first trial of any iteration.

Smallest value of DUSQ attained DPBHI MATRIX on any trial of an iteration that is greater than the current DUSQ.

Largest value of DUSQ attained MATRIX DP2L0 on any trial of an iteration that is less than the current DUSQ.

DP2NMl Value of DUSQ used on previous MATRIX trial or valid step.

Equal to dDPSBAR*DUSQ and set MATRIX DPBTMY only when too many trials for an iteration have been run.

Constant used by RKVS integration BLOCK DRHQ) package set Input minimum allowable integra- INITAL DTMIN( 14) tion step for each stage Defined by equation 100 INITAL DUSQ VALJD MATRIX MATRIX DUSQ required for DBETA DXNUM Input maximum allowable tolerance INITAL EPSLN in stopping conditions Array of predictor-corrector STEP1 ERR(15) errors FPR@G Array of derivatives for the F(15) equations of motion (see appendix PLAC B@Q)M B) MATRIX I Value of I FEFEI w w Input forward trajectory integra- INITAL FINNER tion error limit for Runge-Kutta variable step integration FORTRAN Subroutines Description Name Symbol Units Computed Values of F array from previous EXEC FSAVE (15) integration step Stopping condition indicator EXEC FSIGN VALID Constraint change factor. That INITAL FXTRAl part of the constraint error re- quested for the next iteration.

MATRIX DU2min Minimum allowable DU2 FXTRA2 Not used FXTRA3 FXTRA4 Saved value of time when LAMBDA DVALB DVALB was last called Input mass to weight conversion Go ft/sec2 INITAL factor Input earth gravitational constant GR P times lo-l5 INITAL GUIDN( 10,3) Array of guidance phase names BL@CK plus words used in printout of stopping condition and optirniza- tion parameter HITIM The last, or highest, time stored INITAL in the TIMEU array EXEC UCALC Current minimum step size EXEC HMIN for a given stage RKVS diagnostic print indicator IDIAG BL@CK ID#NE RKVS step completed indicator STEP1 STEP 2 IERR@R RKVS error message indicator BL@CK IN RKVS control parameter (see STEP1 EXEC ref. 8) STEP2 LAMBDA GUIDE PARTIAL Subroutines FORTRAN Description Computed Symbol Units Name Number of meaningful data words INITAL EXEC INC@R in the TIMEU array at a given time UCALC DUSQ ‘bounce test” indicator MATRIX INDBNL Number of consecutive accepted VALID INDSIC (13) iterations since a constraint has been inside its temporary tolerance band.

Too many trials indicator INDTMY INITAL MATRIX End point search mode indicator MATRIX INDZER Printout counter for forward and EXEC INPe, LAMBDA backward trajectory STPl STP2 Storage counter for forward and EXEC INST LAMBDA backward trajectory STPl STP2 INT#T Counter for the total of words EXEC UCALC read into the TIMEU array Packed array containing state TITLES IPF(10) variable indices of the placards Input variable used to control INITAL IPRNTl printout of forward tra jet tory Input variable used to control IPRNTB INITAL printout of influence coefficients occurs during backward inte- gration if NC (7) = 1 Input variable used to control INITAL IST@Ri storage of time and partials along forward trajectory Input variable used to control INITAL IST@R2 storage of adjoints along back- ward trajectory Current iteration counter INITAL ITC MATIX Number of iterations that have INITAL ITR8 been plotted EXEC FORTRAN Subroutines Description Name Symbol Units Computed ITW@27 Constant used by RKVS integra- BLOCK tion package. Computer word size dependent Input array of indices for the INITAL

1X(7)

free initial condition desired Packed array of the plot indices IXKU(G0) TITLES of the state variables, AK’s, and control variables selected IXTRAl Thrust option indicator INITC@, TLD IXTRAB Fuel flow option indicator INITC@, TLD IXTRA3 Aerodynamic option indicator INITCQi, TLD IXTRA4 Not used KDAT Logical tape number assigned to BLOCK the data tape. KDAT contains all input tables except table 0. KDAT is written in INITAL and LAMBDA and read in EXEC and CARDS KINP Logical tape number assigned to BLOCK the input tape, KINP, is read in INITAL KLAM Logical tape number assigned to BL@CK the adjoint tape. KL4M is written in STP2 and contains T, UULAM array, and U array of each point stored during adjoint integration.

KLAM is read in UCALC to build new control history.

K@U T Logical tape number assigned to BL@CK the output tape. K(bUT is written inINITAL, CARDS, STPl, STPB, STEPl, STEPB, MAT@UT, VALID, MATRJX, KCALC, VARIC, PRTIAL, AKSTP KPAR Logical tape number assigned to BLOCK the partial tape. KPAR contains NETA records of T, PFX, PFU, and U. KPAR is written in STPl and read in DVALB. Also, sub- routine UCALC uses KPAR for scratch storage.

Subroutines FORTRAN Description Computed Name Symbol Units KPLT Logical tape number assigned to BL@CK the plot tape. KPLT is written in STP 1 and contains the data to be plotted. KPLT is read by PLOTZ Logical tape number assigned to BLOCK KPUN the punch tape. KPUN is written in CARDS KSCR Logical -tape number assigned to BI@CK the scratch tape. KSCR is written in STPl and contains the control variable history of every accepted valid trajectory.

Logical tape number assigned to KTAN BLOCK the control history overlay tape.

KTAN is written in INITAL and UCALC and contains the part of the control history that would not fit in the TIMEU array. KTAN is read by EXEC.

KTC Index of weight in the IX array onlyINITAL if weight is a free initial condition.

Array of program control indicators LP) (see appendix B) Number of stored integration points STP2 LETA generated during a backward trajectory LAMBDA Printout line counter during the LINES EXEC LAMBDA forward and backward trajectory.

Subroutines STPl and STP2 will STPl STP2 page when LtNES exceeds 57.

LINETl Number of lines (including blank INITAL lines) printed at each printed time during the forward trajectory Number of lines (including blank LINET2 INITAL lines) printed at each printed time during the backward trajectory Total number of equations being LV INITAL integrated during backward trajectory Subroutines FORTRAN Units Symbol Description Name Computed Current closed-loop mode number IMTAL MAL GUIDE Largest value of J in the PFX(I , J) INITAL MI TC array for a given problem Desired maximum number of BLOCK M NDIAG diagnostics printed from RKVS per (computer) stage RKVS parameter that determines M@DZ EXEC!

the method of integration LAMBDA Input total number of (computer) INITAL MSTAGE &ages Two alphanumeric words (STG and BLOCK NAM(2) END) used as indicators during First word deter- reading of data.

mines when stage-dependent data is to be read. Second word deter- mines when all the data for a case has been read.

Input degree of problem solution NARBY INITAL indicator See input description NC(160) INITAL ND(112) Start and end location of all input INITAL stage dependent tables in Z array Starting locations of all input NDS(34) INITAL tables in Z array for a given stage EXEC NE@M Number of equations of motion TITLES being integrated, excluding enroute constraints NERR(15) Array used by RKVS for error con- INITAL trol during integration of forward trajectory NERRSl RKVS error message control BLOCK NERRSB RKVS error message control BLOCK RKVS error message control NERRS3 BL@CK NETA Number of stored integration points EXEC generated during a valid forward trajectory STPl FORTRAN Subroutines De scription Name Symbol Units Computed Number of dated points saved for EXEC NETAP(6) each trajectory to be plotted NETLST Number of integration steps taken IAMBDA for the previous valid forward trajectory Location in Z array that contains NEZ INITAL last data point of last nonstage dependent table Input maximum number of iterations INITAL NITC desired Number of selected auxiliary print- TITLES NK out variables (AK’s) Total number of points in the con- INITAL UCALC N@UP trol table at any time NP(40) Packed array of auxiliary printout TITLES variable indices Number of penalty functions to be TITLES NPF integrated INITAL NPHASE Current closed-loop phase number Number of plots to be made for a TITLES NPLTZ given data case Index of the state variable to be NPP INITAL used for the payoff variable Number of sets of partial deriva- NR;CC LAMBDA tives that will fit in allowable stor- age during backward integration at any one time.

NSE Abort error indicator MATRIX INITAL AKSTP EXEC STPl STP2 NST Index of state variable to be used INITAL as stopping condition NSTAGE Current (computer) stage INITAL Closed-loop guidance table counter. INITAL NTB GUIDE (Initialized to 23 and increm ented each time subroutine GUIDE re- quires a new guidance table) FORTRAN Subroutines Name Units Symbol Description Computed NTRY Current trial number (set to 1 at INITAL MATRIX the start of every iteration) VALID Nu Number of control variables TITLES selected NUMNC Total number of input nonzero TITLES plot indicating NC’s (NC(61) - NC(74)) NUPl Number of control variables INITAL .

selected plus 1 NV Subset of the equations being INITAL integrated in the predictor- corrector mode during the forward trajectory NV4M Location in the FB array preced- INITAL ing the first adjoint derivative NWPR Number of words per record INITAL written on KPAR tape on sub- routine STP 1 Total number of equations of motion TITLES being integrated during the forward trajectory NXUK The sum of the selected number of TITLES state variables, auxiliary printout variables, and control variables plus 1 NBUK Two times NXUK and indicates the TITLES total number of words needed to print the trajectory heading N4 Number of constraints plus the TITLES performance and stopping condition N4Ml N4-1 TITLES N4M2 N4-2 TITLES @MEGA w rad/sec Earth rotational rate INITAL PFU(15,4) Array of 82, /6UJ I = 1, NX, ANPRTL J = 1,NU F’WIX’I’JXAN Subroutines Name Units Symbol Description C om puted PFX( 15,7) Array of a’FI/aXJ I = 1, NX, ANPRTL J = 1, MITC PHIGRN Largest value in magnitude MATRIX VALID that performance parameter has attained on any previous accepted iteration Step-size coefficient with respect PHIK KCXLC to performance PHIMl Performance value for the last VALID accepted iteration PHINL Nonlinearity of performance KCALC variable PHINLD Desired nonlinearity of performance BLOCK PLTNAM(28) Array containing hollerith data TITLES used for plot titles Each plot title is defined by four words of this array PMPH 8 mach/c?altitude AKSTP PMPV amach/dvelocity AKSTP PQPH d dynamic pressure/aaltitude AKSTP PQPM adynamic pressure/amach AKSTP PQPV adynamic pressure/bvelocity AKSTP PSI( 13) Input desired value(s) for terminal INITAL constraint(s) PSIMl(13) Terminal constraint value(s) for VALID previous accepted iteration RAD Degrees to radians conversion BLOCK factor RADIC Absolute value of first Lagrange MATRIX multiplier ft Input earth radius INITAL Array used to transmit placard PLAC value and derivatives of the placard input table (see subroutine PLAC) FORTRAN Subroutines Name Symbol units Description Computed SAKll Sin Q! Sine of the angle of attack AKSTP Sin 8 Sine of the pitch angle (set to zero SIl C@NTR if NC(77) = 0) Sin tp Sine of the bank angle (set to zero S12 C@NTR if NC(78) = 0) Sine of the flight path angle AKSTP SK3 Sin yR Sin p Sine of the latitude (set to zero AKSTP SK5 if NC(125) = 0) Sine of the heading angle (set to SK6 Sin $R AKSTP zero if NC(126) = 0) SREF(14) S ft2 Input aerodynamic reference area INITAL for each stage Trajectory time measure from T t set EXEC LAMBDA zero at launch STEP1 STEP2 set Total time increment to be inte- EXEC LAMBDA TI grated for a stage The first 1000 words of the input TIMEU( 1000) INITAL EXEC control table. The remainder of UCALC the control table (if any) is written on KTAN in record sizes of NUPI words Alphanumeric case title TITL( 12) INITAL TOMAX set Terminal trajectory time -up- MATRIX dated every accepted iteration TSTP(14) set Input stopping time for each stage INITAL measured from time zero at launch Temporary tolerances (one for each VALID TToL(13) - constraint) U Array of control variables CQ)NTR GUIDE U(4) Ull Pitch angle (set to zero if C@NTR deg NC(77) = 0) Bank angle (set to zero if UI2 C@NTR cp deg NC(78) = 0) Throttling parameter (set to 1 if U13 C@NTR GUIDE rl NC(79) = 0) Sweep angle (set to zero if UI4 A C@NTR deg NC(80) = 0) Subroutines FORTRAN units Description Name Symbol Computed Array of variables (see appendixB) AKSTP VAR(50) TLD INITAL Automatic weighting matrix LAMBDA

mN(4)

wii elements Input weighting matrix elements DVALB VINP(4) lb/set Input inert weight flow for each INITAL WD@T( 14) *inert stage lb Input initial weight of each stage INITAL W&14) Array of state variables being INITCO STEP1 X(15) integrated EXEC INITAL Input weight flow correction factor XKM (14) INITAL for each stage XKT( 14) Input total impulse correction fac- INITAL tor for each stage lb XKl W Weight (set to initial weight if INITCQ) NC(121) = 0) XNAME(19) Alphanumeric adjoint headings BLOCK Floating point value equivalent XNC9 INITAL to NC(9) Numerator of first Lagrange XNUM MATRIX multiplier Input initial values Of state INITAL X0(15) variables X&M AX( 7) Input maximum allowable value of INITAL free initial condition Input minimum allowable value of INITAL XJ;?IMIN(7) free initial condition XRHJb Constant used by RKVS integration BLOCK package XSAVES(15) Terminal value of state variables VALID from previous valid trajectory XSTP Input desired value of stopping INITAL STPl condition XXLAM (15,14) Influence coefficients corrected DVALB A for stopping time variation FORTRAN Subroutines Name Units Symbol Description Computed Alphanumeric adjoint headings YNAME(14) BL@CK z (7000) General data storage array for INITAL EXEC non&age-dependent tables (per CARDS UCALC stage) during forward trajectory.

Used for storage of partial sets during backward integration (See program equivalencing) Output Description The program produces three types of output: the printout, punched cards, and magnetic tape.

Printed output. -The printout gives listings of the basic input data, controls for the program, trajectory data for both forward and backward integrations, con- vergence information, a control variable printout, and a trajectory summary. A partial derivative check is printed out during the nominal forward trajectory if selected.

Preliminary trajectory printout: Page 1 Program title Data case title Initial condition data Stage dependent parameters Constraint dependent parameters Optimization parameter Stopping parameter Free initial condition parameters Guidance phase dependent parameters Page 2 Control card (see card set 2 ) NC array Page 3 (Number of pages are dependent on number and size of tables ) The input tables are printed in the order they are input The nonstage dependent tables are printed first, followed by the tables for each stage.

Forward trajectory: The forward trajectory printout follows the tabular out- put. At the top of each page, the names of each variable selected for printout The appear. These headings are printed eight per line, beginning with TIME.

state, auxiliary, and control variable headings are printed next, beginning with the first state variable and ending with the last control variable. Any heading followed by an asterisk (9 indicates metric units for that variable. The print- out of the forward trajectory data occurs every IPRNTl integration steps, with two blank lines between each printed step. A line counter, assuming 57 lines per page, ensures that each page will begin with the headings described above. For data cases with more than one stage, the beginning of each stage starts a new page.

Partial derivative check printout: The forward trajectory printout format is altered if the numerical partial check option is being used. The partial deriv- ative printout will occur every NC (12) stored integration steps immediately following the forward trajectory printout of the same time point. The’format for The word TIME and the current time are printed on the partial check follows: the first line. The partials are printed eight per line with lPEl4.7 format.

The analytical partial derivatives of the first equation of motion with respect to the first MITC state variables are printed first, followed by the equivalent numerical partials. This is repeated for each of NX equations of motion in fol- lowing rows. The printout of the partial derivatives with respect to the state variables has the form of the array PFX (I, J) where I (the row) represents the equation of motion and J (the column) represents the derivative state variable.

The analytical partial derivatives of the NX equations of motion with respect to the first control variable are printed next, followed by the equivalent numerical partials. This is repeated for NU control variables. The printout of the partial derivatives of the equations of motion with respect to the control variables has the form of the array PFU (I, J) where I (the column) represents the equation of motion and J (the row) represents the derivative control variable.

The backward trajectory is printed only if Backward trajectory printout: NC (7 ) is nonzero. The backward trajectory begins at the final time of the valid step and integrates to the trajectory starting time, printing the influence coeffi- cients and instantaneous impulse responses every IPRNTB integration steps.

Both the influence coefficients and impulse responses are corrected for vari- ations in the stopping condition. At the top of each page is printed a block of headings identifying each variable. The work TIME and the current time are printed on the first line. The second line consists of the influence coefficients for MITC state variables on performance. Next are printed the influence coef- ficients for MITC state variables on the first constraint, followed by the influence coefficients for MITC state variables on the second constraint, etc., through all the constraints. Following the influence coefficient printout, the instantaneous impulse responses for the first control variable on performance, on the first con- straint, the second constraint, etc., are printed. The next row of print is the impulse response for the second control variable on performance, on the first constraint, the second constraint, etc., through all constraints. This is repeated for all the control variables selected. The backward trajectory printout uses the same line count and paging control as the forward trajectory printout.

Matrix printout (convergence information) : The convergence printout gives information used by the stepsize controller in making decisions concerning trial trajectories, valid steps, etc. The printout is given as notes and as arrays of values used by the program. The information aids the user in determining the The printout information and notes are given below. The nature of convergence.

underlined statements and statements preceding a colon are for clarification in describing regions of the printout and are not part of the printout.

Printout from convergence logic after a valid step or trial.

Iteration number Trial number Number of integration points stored during forward trajectory integration (NETA) Number of integration points stored during backward trajectory integration (LETA) Printout from convergence logic after a valid step only, Majority vote counter Messages printed under special conditions are: “Number of tries plus rejected valid steps for current iteration exceeds 10 -abort” 7YI’oo many functions with adverse travel” “Adverse phi too great” “DUSQ modified because of too many tries: “Reject valid step” Printout from convergence logic after a reverse integration only.

I matrix Minimum allowable DUSQ Denominator of first Lagrange multiplier Diagonal elements of automatic weighting matrix *Derivatives of performance with respect to variable initial condition (DPDX) *Elements of DPDX (SUM, XXLAM) Note: (I) The items with * are printed only when the variable initial condition option is selected, NC(5) # 0.

(2) XXLAM is the performance change resulting from a unit step in X0 (first term in the bracket in equation 127) and STJM is the performance change that would result if the control function were modified to bring the constraint vector back to the state prior to the application of X0 (second term in the bracket in equation 127).

Messages printed under special conditions are: “Convergence failed” “Converging but gradient too large for extremal” “The last trajectory appears to be an extremal” “Gradient of phi too negative to continue” “Overflow occurred in matrix inversion-abort” “Singular matrix-abort” Printout from convergence logic after a trial only.

Majority vote counter *Maximum permissible changes in constraints *Controlling function and the step-size coefficient *Function or functions reducing step-size coefficient and the new value of the step-size coefficient *Final step-size coefficient Note: The notes above with * are written only if the majority vote is 1 zero Messages printed under special conditions are: “Controlling with constraint within tolerance band which reduces -repeat-step-size-coefficient calculation” the step-size coefficient “Too many functions with adverse travel” “Adverse phi too great” “End point search and step-size coefficient greater than 2 - force valid step” “Penalty function violation - force valid step” “Performance improved and step-size coefficient equals 0.5, accept current trial and force valid step” “All constraints outside tolerance band improved but DUSQ reduced too much, accept current trial and force valid step” TfDUSQ was increasing and wants to decrease accept previous trial and force valid step” “DUSQ was decreasing and wants to increase accept current trial and force valid step” Printout from convergence logic after a reverse integration, trial, or a rejected valid step.

Step-size (DUSQ) required for A B

Numerator of first Lagrange multiplier Lagrange multipliers Accuracy check for step-size (DUSQ) (this value should agree with the value of DUSQ for next trajectory) Heading block of constraint and performance names Current end point values of constraints and performance Previous valid step end point values of constraints and performance Change in constraint and performance end point values Constraint travel indicator (0 if inside tolerance band, 1 if constraint change is toward desired constraint value, and -1 if constraint change is away from desired constraint value) Constraint tolerance indicators (number of consecutive iterations a constraint has been outside its tolerance band) Allowable nondimensional forward constraint change (in the direc- tion of the desired constraint value) Allowable nondimensional backward constraint change (in the direction away from the desired constraint value) Nondimensional constraint change asked for (percentage of constraint error to be eliminated by this iteration) Constraint tolerance bands Nonlinearities of constraints and performance Step-size coefficients of constraints and performance Maximum step-size coefficients of constraints (based on maximum permissible constraint travel) Predicted changes in constraint and performance end points (total) Predicted changes in constraint and performance end points due to step-size (DUSQ ) Predicted changes in constraint and performance end points due to changes in variable initial conditions Changes in variable initial conditions to be made on next trial or valid step (printed only if variable initial condition option is selected, NC@) # 0) Final step-size coefficient DUSQ for next trajectory Messages printed under special conditions: “DUSQ too small to meet DBETA” “DUSQ set equal to DUSQ required for DBETA” “End point search” “Scale DBETA to match DUSQ” “DUSQ set to minimum allowable value” “DUSQ too small for optimization -DUSQ reset” “G tries are too many -force valid step” “8 tries are too many -force valid step” “Next trajectory will be a trial” “Next trajectory will be a valid step” The control variable table and the Control variable table and restart table: restart table are printed following the last iteration. The restart table is not printed following the running of a nominal trajectory only.

Trajectory summary: The trajectory summary follows the printout of the control variable table and restart table if one or more iterations have been run.

The heading block that appears at the beginning of the trajectory summary is identical to the one which appears during the printout of the forward trajectory.

The format of the trajectory summary is also identical to the forward trajectory format. The summary contains the end points for each successful iteration.

Punched card output. - The program produces punched cards in the form of control variable and restart tables. These tables are used to restart a problem from the point where it was stopped.

Punched output is produced when either or both of the following options are selected: NC(8) >OandNARBY=2 1) 2) NC(8) >OandNC(13)>0 The control table is punched six words per card in octal format (6812) and the restart table is punched in 2413/6E12.5. Both tables are preceded by title cards in the format required for input to the program. A restart table is not punched following a nominal trajectory.

The restart table consists of previous valid step convergence information, used by the MATRIX subroutine. Included in this table are DUSQ, performance and constraint end points, constraint temporary tolerance bands, and the stopping time of the previous valid step.

Note: If punched cards are not produced due to machine error or program abort, the control variable table can be recovered from the magnetic tape (KSCR).

Restart table output description: Card 1 Name Column Description l-5 Restart table number (0) 6-10 Number of cards in restart table (title card not in- cluded in this count) 12-23 Alphanumeric word RESTART 24-47 Case title as input on title card of data case Format 2413 Card 2 Column Name Description l-3 INDTMY Forced valid step indicator = 0 valid step not forced = 1 valid step forced because too many,trials were run = 2 valid step forced because control logic is using the minimum allowable value of DUSQ 4-6 INDSIC (1) INDSIC (2) Number of consecutive iterations a constraint has , been outside its respective temporary tolerance band.

10-12 INDSIC (3) A value is present for each constraint (None of these I 1 I I values appear for a no-constraint problem. ) 1 t I I 40142 INDSIC(13) Card 3 Format 6E12.5 Column Name Description 1-12 DUSQ Value of DUSQ used for -last successful valid step 13-24 FXTRAl Amount of nondimensional constraint error asked for on the last valid step. (This value will be zero for a no-constraint problem. ) This is the maximum amount of constraint error that can be asked for by any one constraint.

25-36 DPBTMY (DPBBAR x DUSQ)1/2 where DPBBAR is the value of DUSQ used on the second-to-last su&essful valid step and DUSQ is the current value of DUSQ.

Note: DUSQ is set equal to DP2TMY only when the indicator INDTMY has been set equal to 1.

T@MAX Trajectory time for last successful valid step.

37-48 49-60 PHIMl Value of the performance function on the second-to- 1 I last successful valid step.

1 1 I 1 I 1 1 I (The following values occur only if the problem has one or more constraints. ) 61- 72 Terminal value of the first constraint on the second- PSIMl(l) to-last successful valid step.

Card 4 Format 6E12. 5 (If needed) Column Name Description l-12 The temporary tolerance band corresponding to the TT@L(l) first constraint (expressed as half of the total band width).

13-24 CPSI(1) The amount of nondimensional constraint error asked for on the last valid step.

Note: This value is normally equal to FXTRAl but under the conditions mentioned in appendix C, the value may be 0.1.

25-36 PSIMl(2) 37-48 Values for second constraint TTtiL(2) 49-60 CPSI(2) Special notes : The value of DUSQ input through the use of a

1)

restart table will override the value of DUSQ input in the initial condition data section of the input; The values of the temporary tolerances input 2) through the restart table will override those normally computed by the control logic using the final input constraint tolerances; 3) If it is desired to change the amount of constraint error asked for all that is necessary is to change the value of FXTRAl.

(FXTRAl must never be greater than 1. ) The individual values of CPSI are set equal to FXTRAl or . 1 (under special conditions men- tioned above).

The values of CPSI never need to be reset.

Magnetic tape output. - The program produces a magnetic tape when plotting (KPLT) and another magnetic tape when iterating (KSCR).

Plot tape (KPLT): The tape KPLT is produced when NC(11) # 0. This tape contains the data for the variables selected for plotting. For each valid trajectory that is plotted, NETA records are written on KPLT in binary format. The number of points stored when plotting (NETA) must be less than or equal to 1000. (The absolute value of NC(l1) defines which iterations will be plotted. ) The number of words written per record is two times the number of plots to be made, or the num- ber of nonzero elements in the NC array from NC(61) through NC( 74). The plot data is written on tape every ISTORl integration steps in the order selected by input.

Control table tape (KSCR): The tape KSCR contains the control variable history for the nominal and each accepted valid trajectory of a data case in binary format. This tape is generated when NC (8) > 0; and NARBY = 2 and/or NC (13) > 0. For each valid trajectory there are NETA records written on tape, where each record consists of TIME and NU control variables. This tape (KSCR) is used to recover the control table for any iteration, in the event a data case is interrupted and doesn’t punch cards.

The procedure used to recover a control table is described below.

A short FORTRAN IV program called ffTHETAff is included in the source decks for the STOP program. This program reads the binary tape (KSCR) on logical unit 2 and produces a listing and punched cards (OCTAL FORMAT 61d’12) of the control variable table selected by input. Input to this program is described below.

CARD NUMBER 1 (315, 4A6) CCLS 1-5, ITC, Number of control tables written on KSCR UP to, but not including, the desired table. (The nominal and every accepted valid step will have a control table written on tape.) Equal to ITC printed immedi- ately following desired valid step.

CCLS 6-10, NU, Number of control variables used.

CCLS 11-15, NTABLE, Number of copies of control table desired.

COLS 21-44, TITLE Any information desired on the title card of the punched control table.

CARD NUMBER 2 (145) COLS l-5, NETAS(1) Number of points stored during nominal 6-10, NETAS(2 ) Number of points stored during first valid trajec- tory.

11-15, NETAS(3) Number of points stored during second valid trajectory.

16-20, NETAS( 4) Number of points stored during third valid trajec- tory.

l l l l l l l NETAS( ITC ) Number of points stored during valid trajectory preceding the one from which restart is desired.

NETA Number of points stored during valid trajectory from which restart is desired.

Both the printed and punched tables produced contain one additional point at the end of the data to prevent extrapolation off the end of the table.

Sample Problem A sample problem is given which demonstrates the data setup and output for One iteration is given showing the .’ a simple, single-stage, air-launched rocket.

nominal trajectory, the control logic output, one trial, a valid step, an output control table and restart table, and the trajectory summary.

Statement of problem. -The flight path and control variable history are required to maximize the final weight of the rocket with the following initial con- ditions.

= 1500 pounds wO = 500 feet hO = 60 degrees yO = 600 fPs V The terminal constraints imposed are: = 2000 fps vf h = 1200 feet f The final velocity is selected as the stopping condition. The initial weight and flight path angles are free to change to improve performance. The initial conditions W, and y, above are the starting values.

The The aerodynamics and propulsion data are given in tables II and III.

reference area is 2 square feet.

TABLE II Aerodynamic Data cD M .2 .5 .075 .4 .5 .075 1.5 .4 .075 10 .5 TABLE III Propulsion Data Thrust 7000 pounds Specific impulse 280 seconds Nozzle exit area 144 square inches Sample input. -The statement of problem is transformed into data input cards through the use of the STOP input instructions given under “Input Data Preparation. ” The data input cards are listed below.

EXAMPLE ROCKET CHECK CASE 1 5 1 4 --+w-- 3 2 19 I3 .2 16 0 17 3 77 1 7 18 4 19 1 20 2 33 1; 47 81 1 82 1 87 1 84 1 85 1 87 1 171 1 122 1 123 124 1 91 1 32.174 2000. 14.08171812090Z992. Oa -100. 18. . 18. 0001 ~ .- 1500a 500e 609 600e 1500e 0.

-. -- l* O* 85 400.

1 THETA TABLE Oa 60m 100. 60e 1 1 WEIGHTING TABLE 0. 16 lOO* 1.

-. ~ STG 31 2 ENGINE TABLE 144m 280, lOO* 7000. 144.

OI 7000.

280a 33 2 AERO TABLE 00 a2 05 9075 1.5 l 4 -5 04 05 9075 10.

-The results of the problem for one iteration are given Sample output.

below. A complete description of the printed output is given in the section on output description.

--- .---_- ---._ .--.--..- __.._.. ____ SUPERSCNIC TR4NSPORT OPTIMIZATION PROCRPM (STDPI THETA CONTROL 7094 IRSVS 13 VERSICK Eg)rPLE ROCKET CFECK CI\SE ~-___-_____ ----- - -.-_- -- - -. . - ._ -. . ._ ..-_ _ __ INITIAL CONDITIGN DATA - ________ -_--___ -__--_ .--_ _ = .?.OOOODOOE 03 CR = 1.40E1718E 16 RO = 2.0902992E 07 OMEG4 = O.GUOOUUOE-39 GO = 3.2174WCE Gl XSTP 1.0000000E 02 FIhNER= 1 e~.O~u~C.E__~.L__.~~N.~~~~._l,IOI)’)C~?E_~l . ..FRSLN_ .~.1.-,O.~ou.~cE._C3. _.MLrP .=,C.C~j?C’O((?E-39 _ _, _ _ _ DUSO p.I--.- -.------.1.

X0( =,I= b.CCOP%DE 01 TSTART=-0.0000000E-39 L =-O.COCOOOOE-39 XOI I)= 1.5000000E 03 XIII 21= 5.ooo1)ucck 02 X0( 4)= 6.00COOOOE 02 X0( .___----_----._. ___ STAGE DEPENDENT PARAMETERS IrE ICHT TSTP SREF DELC XKT XK’! wnuf UTPIN STAGE --. . -- ._ .--.. 2;-~ob~~.-~--~.;o.~o~~ --‘---I;ob’Bab----- i.obbii’o--. ..-.- .

L;. l(‘,,j(),i ...

I, .cl~oc:~ 1 1500.00000 4CC.00000 CONSTRAINT DEPENDENT PARAPETERS -___-_.-- .-_. --__----- DPSIIIAX VARIABLE PSI - .-_-- - -..-__- ---- -----. .- .-_. -_. . _.

ALTI T%E l2ood.coocc 1CO.CcC00 hEIGHT HAXIYIZE -_---.

STOP ON VELOCITY R _--__ FREE INITIAL CONDITION PPRAMETERS --.-.---.-. ._ z Ip XMAX XM IN DELX VARIABLE _--- .---_-_-- --__- -...... -.. --.

IdEIGHT 2000.OCCCO lOoO.COCcO -rlsmGm0 GAMMA R 70.00000 50. cc000 2.00000 --- .-____.__----.- --..

GUIOANCE PHASE DEPENDENT PARACETERS PHASE. MCIDE- SkITCH VARIABLE SWITCH VfUE ALPHA MIN ALPHA YAX THROTTLE YIN THROTTLE MAX -__---_ -_-._-.- ..- - _.._ -__- _.. ,_..... __,_. _.______^_ ALPHA INCREASE TO VELOCITV R 1000.3n0c~0 -0.0000G -0.00070 -lJ.000P0 1 -0. ooocc --. -- -. ___--__-.-.__ .._- -.-. --- 2 ALPHA INCREASE TO TIME 400.00000 -o.ouoco -0.00000 -o.couc’l -C.C9OC% -Ts&q ._. .-__ __ --__--..-..---.- ._____-..- .-_____-.-..-_.-.._-.- .._ ____- .___._ --.-._-_ INTEGER TABLES _ __-..- _.___... -__-- -__--- _-. ___ --__-.-.----_- __I. _- 1----------- +fsrAcE = IPRNTI = 1 IPRNTZ = 5 ISTORI = I lSTCR2 = 1 NARdY = 2 NlTC = 1 ~cmIm&p---- ..__. - .-----.- - ---.---.- -.-.- .-- --- -------- ----.--..- - -...-.. .--..- NCI II= 0 NCI 2l= 0 NC1 3l= 0 NCt 4l= 3 NCI 5l= 2 NCI 61; 0 NC1 7)=-O NC t 9l= I NCI 9l= 1 WI lOI= 0 = NC1 111 0 NCl 12) 0 NCI 131= L_“CC.L1_4I.~~O__NCL.U~~UC_~_16r~O!C(._1Z)~=~._.~C( .ldlY +- NCL?_9)~...-~-~c!-~~~=-?~ .- NCI 21,: 0 NCI 22l= 0 NC1 23l= 0 hCI 24). 0 NCI 251= ‘3 YCI 2bl= 0 NCt 271= 0 NC{ 2YI= II NCI 29l= 0 NCI 3?1= c.

NCI 31l= 0 NCI 32l= 0 NC1 ?3l=lO NCI 34)= 0 WI 35l= 0 YCI ?bl= 0 NCl 37)= 0 NC1 3811 0 NCI 391. 0 MI 4Ul= F XI 4bl= 3 NCI 471: 3 WI 4ttJ= 0 4ci. 4Ql= 6 NCI 41l= 0 FICI.~Jr_~.NC1._43L~_C._hCL 44)= 0. NC{ 45l= 0 Y’cI 50): c .._. -.. -.

h;Cl 5il= C!

NCI 51l= 0 WI 5211 0 NCI 53l= C hCl 54). C NCI 55)= 0 NC1 56)= 0 NC( 57J= r) NC{ 581. 0 NC1 bOl= C NCI bll= 0 NCI 621= 0 NCI .t31= 0 &Cl C41= C NC( 65l= 0 htl tbl= 0 NC( bTI= 0 Xl hAI= 0 NCI tSl= 0 hCl 7n1= c = 0 hC( 72l= 0 N~_/4.l=-~_,NCl_75)~.0.._ NC{ 761=..?-‘JCI. 77!=_.l_NCI .7P!? O--.hCt .7?l.=.._C_.P;C!_Pc)!.~..!-___- NC-1 7311 0 -NC1 71) NCI 81l= 1 NCI 8215 1 NCI i?3l= I hC1 84)~ I NC1 .35l= 1 ;YCf !!bl= 0 NC1 87l= 1 NC( 89l= c htl 8$1= 0 VCI $01= c 1 NC( 92)~ 0 NC1 53)~ 0 hCI 54l= 0 MI 95)= 0 NCI 96)= 0 NCI 97J= C NCI 9Sj-5 0 WI 55)= I! NCflOrJl- is Ntl 91l= _ NCl~~).= 0 hC.~l~2,~~-_N~.!~.C,3l~O_- hC(l.C4); O~..NCllO51= !l..yC(lObl= 0 NCI 1371.= C .UCl 109l.= 0. hCIlC9)= ‘I- ‘!C.l.llq)= C __.__.. ___ 0 hCl112l= 0 KCllL?l= 0 hC( 114)= 0 NCll151= 0 NClllhl= 0 Nc11171= c Xl 1191= 0 kCI119~= n NCllZOl= lj NCfllll= NCllZll= 1 NCl122l= 1 hCl123l= 1 hC(124)= 1 NCl125l= 0 hCl1261= 0 NC(1271= C hCi12Jl= 0 NCI lZSl= 0 NCl13Ul= 0 YCl1371~~,~. NCl131l= 0 NCl132)= 0 NCl133!=~._h_C~~3_41= 0 NCl1351= 0 N!Cl 13R!:.-C-?!Cl 13,91.=..1’__.. NTl.~~?l.= 0 NFl.l!CI= 0 ..--_._. .__ -_ _--._ NCI 141)= 0 NC1142l= 0 NCt14?1= 0 NCl144)= 0 NCl145)= 0 YCI 14bJ= 0 NCI 147l= C wzilca)= c hC1145)= r, NC I 153 I= f-s 0 NCt1551= 0 !dCl15bl= 0 NCl157l= C W-ul15Al= (r NCl155l? C NCI lbOl= 0 NC(151l= 0 NCl152l= 0 hCl15?1= 0 NCl154l= -.____-.-. __ __.__ -_._-__--_---.---.--------.-----.----.

--- -- G -- cn NON-STAGE OEPENOENT TAELES TABLE 0. 1 CARD(S) THETA fA8L E 0.000000 60.000000 100 .cooooo bO.OOi3000 -0.000000 -0.000000 -o.ooo_oco TABLE 1. 1 CARCISI ME IGbTING TABLE 0.000000 1.cooooo 100.c00000 1.000000 -0.000000 -0.00G000 -c. 000000 STAGE DEPENDENT TA8LES SiAGE 1 TABLE 31. 2 CARDIS) ENGINE TABLE 0.000000 7000.0000c0 144.cooooo 280.060000 7000.000000 -0 .cncooo -u.oooooo 100: 0 0 oo”ocLJ”o” 00 0.000000 1”4ao+,“O”,’ . c 280.000000 -0.000000 -- -- TABLE 33. 2 C*RolSl AERC TAELE 0.000000 0.200000 o.5coooo 0.075000 1.500000 P.401)OOO o.5cooPo 0.075000 10.000000 0.400000 0.500000 0.075000 -0.000000 -c..ooooco --__----- -. -. ------.- -_--- _--__- .------ _--- __--- -- LIFT GAFlMPLR -mu.--VECUCLIY-R.- THRUST D!G(dC TIVF AEU;tiI--ULLll!2E ALPHA PITCH ANGLE 0 MACH NUMBER ;IE IGHr FLCW 1500.000 60.@30 229.249 mo..coP--~ - 6CO.C~Om ____.__C~2l_.~2.?-~-_0,P_nr-_-~- -.. --- _-.--- _____ _ ~o...ow.~~- . -_- 0.000 60.000 421.757 0.538 25.oco b 3&..5_69 4941.750 49.412 261.514 1487.5&C 7 6 6.7 ZIZ70 1 O.5ao 0.699 60.000 471.027 0.572 25.coc 4973 434 3eL8.30.!

1.L5U14&3.7%-- LLe..l!tL 5a.ftB3 - ___. -__ b92eO75 ---.-~..--_12_5_t!45 60.000 549. a01 0.622 25.CCO 1.517 2.375 1440.625 ia93,w 776.490 5 023-.. 4 237.181 4l&467 5 7_..aLp --.--- 677.911 0.700 25.000 2.332 60.000 4.062 1398.438 3Ce_LS?9- 306.9OC --.- 591.8&-.--.--- 000.706 0.819 25.COO 5.1,4?p....____l)~~-~~~___5109.458--_ 2.902 .

5.750 1356.25C 4461 .oa9 57.013 1035.342 5203.335 5oc.c92 782.826 _-.- 1116.151 0.942 25.000 2.987 bO.CW 6.250 1343.75c 56.574 1074.978 5232,699 C.COC 785.@42 -------ilI37.299 0.979 4s_oZ,Tj~L--_ 25.000 0.000 . . - ---. 56.574 - ---___ -.._. -- 6.583 1335.417 5204.788 56.262 1101.560 5252.571 O.O’lC P25.16A -_- 1235.379 1.005 25.COC 0.000 56.262 7.083 1322.917 5669.E62 -- 55.803 - -.^- - 1141.597 --._.- --_---.--.-______ 5282.aco O.OCb 887.150 1308.197 1.043 25.OOC o.coo 55.803 7.833 1304. Lb7 6393.E07 55.135 1201.988 5329.009 U.CO@ 983.961 ---- 1418.583 1.101 25.&C 0.000 5 5x35 a.958 i276.042 7538.446 54.174 1293.208 5400.005 o.ooc l137.t-‘28 ^ ---- - .-_-__---.- --.__I-._-..

d-ooo~.-.~.-. 54.174 1585.342 1.190 25.COG 10.646 1233.8fr4 9 3 85.487 52 .-S_l7 1431.095 5509.389 o.ooe 13sl.lc2 --- 1833.245 1.325 25.COO 0.000 52.817 13.177 1170.573 12440.834 50.943 1639.415 5676.950 o.coc 1747.49r 2ia4.373 .____--_-.--___- 1.536 2s.coo -.- 0.000 Sd;94~-.----‘--- 15.708 1107.292 15829.469 49.234 1954.644 5044.557 o.ocr 2CJs 2 b 2505.158 1.760 25.000 0.000 49.234 17.345 1066.387 la198.317 48.207 ~-~zoco,o_o_o 5951.141 0.000 215u.912 -__ ---- -_---__ 2693.640 1.915 25.000 0.000 48.207 ___-- __- __--.-. ____ - liERATION 0 TRIAL 0 NETA= 26 LETA= 0 ____--_ __.. - ..__ - _._._.. -_._---. -.-~ - ___ - - .--- -- -.~-. --- I MArRlX 0.33_~iP_b_8E-O1 -0.517~~~g8_!E~cl -0.51741087E 01 0.29956934E 04 IrlNlWM ALLOUABLE DUSO 0.12824789E 03 I PHI PSI*1 PSI PSl_lNyERSE+I qS1 P4.1. il,89365292Ey02 _ ____.- . -_-.. ..-- 0.24814069F-01 DENO?llNATOR OF FIRST LAGRANGE PULTIPLIEP AUTOYATIC UEIGHTING YATRIX ELEHEhTS 0.10000000E q1 0.7051853rE O!, VARIC CHECKO~IT PRexr - DPDXISM~XXLAM 0.72270353E 00 -0.ltEl8235E-01 ----__-- - .-._ -___-__---- __- -..

~-_-.- 0.260563156 00 -0.33655244E-01 0.22689191E 00 VARIC CHECKOUT PRINTOUT - DPOX,SL~.XYLAf4 DUSO REOUIRED FOR ORETA 0.58396113E 03 NUMERATOR OF FIRST LAGRANGE P~~LTIPLIEH -0.38396113E 03 ___ __-.---- __ _ - ---- - ._.. ---- DUSO TOO SMALL TO MEET DRETA SCALE DEETA TO MATCH DUSO OUSQ TOO SMALL FOR OPTIMILATIGN-CLSO RESET -.- DUSC REOIJIRED FOR OBETA CI.~OOCOCOCE 03 0.21WOOUOE 03 NUMERATOR OF FIRST LAGRANGE PULTIPLIER LAGRANGE MULTIPLIERS C.91994244E 02 -0.99493583E-01 --.___- -.- - .-.._ - --~ - -___ ACCURACY CHECK FOR DUSP 0.40955999E--03 AI TII CURRENT END PCINT V4LUES “t----r%%520 18198.317139 .

0.000000 o.ooooco PREVIOUS VALlO STEP END PCINT VALUES CHANGE IN END POINTS o.oocJooo 1066.386520 --.

CONSTRAINT TRAVEL INiTnATOR IS 1 23159 CONSTRAINT CROSS-OVES lNOlCATCR(S) 1 ALLOUABLE FORYARD CCNSTRAINT CbbhGE(N-DI 0.000000 ALL~S~~T~A~~~.D.~S~R~TCH~~G~~N-O~ 0. US@00 CONSTRAINT CHANGE ASKED FORlk-Cl 0.05cc00 CONSTRAINT TOLERANCEIS) 1000.000000 - _---- Ftd-Ti Nz%k-fT I E 5 0.000000 0.000000 0.000000 STEP-SIZE COEFFLCIENTS c~.oocGoo f4AXlMlJM STEP-SIZE COEFFIClENTt 51 0.000000 ._,- ____. -__- ._-____ -181.370571 45.195420 PRECICTED END POINT CHANGE(S) (TOTAL1 -774.0404RZ 3.619660 PRECICTED EN0 POINT CHANCE(S) (DIJSG) 592.669907 41.575759 PREDICTED END POINT CHANGE(S) IVICI ~- .- __-.

CHANGE IN UEI Cl-IT VIC 58.522521 CHANGE IN GAUMA R VIC 1.170450 FINAL STEP SIZE COEFFICIENT 0. oocooo 409.999996 OUSO’ FOR NEXT TRAJECTORY NEXT TRAJECTORY WILL BE A TRIAL THRUST UFT -rL!!E- wEutiLALuBlnEC~bzLauJx-E WAC- cl MACH NUMBER WE ICHT FLGW 4LPHA PITCH ANGLE -2.000 ___._ 155_8.5z?~.5.00 ..COO.~~~ 6.~.-17_0.-----~CO.NQ- C-S~.L..~~~~ e-e-2%-32?-- ____ 2?9,764 ----_-_ 421.757 0.538 25 .ooo 4.614 65.785 ,a,l_fi I 1104.493 I aE.L,tP1 44.945 2000.0~0 5949.163 131p.1362 2319.637 2697.753 1.915 25.COO 3.262 48.207 -- .- -. _.___. --- ITERATION 1 TRIAL 1 NETA= 0 LETA= 73 --. .____-- MAJORtTY VOTE TEST l ** ENTERING KCALC *+* -- -~ NAXIHUN PERMISSABLE CHANGE Ih ALTITUOE -1503.091888 CONTROLLING WITH UE IGHT STEP-SIZE COEFFI_C_IENT __~ ia.

FINAL STEP-SIZE COEFFICIENT 1.912763 l ** LEAVING KCALC l ** OUSO REWIRED FOR OBETA- 0.99495189E ‘03 G NWERATOR OF FIRST LAGRANGE PULTIPLLER 0.50510020E 03 W OUSO TOE SMALL FOR OPTlMIZATICh-CUSO RESET OUSC REWIRED FOR OBETA 0.99495189E 03 ..___-_-_--_ .-- --... _.

-. - NUMERATOR OF FIRST L4GRANGE CULTIPLIER 0.1044b995E 04 LAGRANGE MULTIPLIERS 0.20513551E 03 -0.22191217E 00 ACCURACY CHECK FOR DUSP 0.20396513E 04 ALTITUOE UE IGHT CURRENT END POINT VALUES -1_R1_SJ :. 6013 1.8 1104,493439 ----i%fi&WS VALIO STEP END POINT VALbEi lal98.317139 1066.386520 CHANGF IN ENC POINTS -45.715820 38.106918 CONSTRAINT TRAVEL INOICATOR(St -_ ~___ -----.

CONSTRAINT CROSS-OVER INDICATOR(S) 1 ALLGilABLE FORYARO CONSTRAINT ChAhGE(N-U) 4. a50000 ALLOWABLE aACKWARC CONSTRAINT CHAhGElh-Cl 1 .ccoooo ------t-ONmb-t%T -G-AT% ~ti--l=EKtfX-l o.o5coor) CONSTRAINT TCLERANCE(SI 1000.000000 NON-L1 NEAR LT IES -0.747942 -0.156841 -- -._- ~ FTFP-SIZE COEFFICIENTS 0.401no 1.912763 MAXIMUM STEP-SIZE COEFFICIENTtSI 1o.ccoooo -552.795692 PRECICTED END POINT CHANGE(S) (TOTAL) e7.597941 -- - - -~ ..___ -- PREOICTED END POINT CHANGEIS) (OUSC) -1726.432968 8.073352 PRECICTED EN0 POINT CHANGE(S) (VICI 1133.637283 79.524590 CHANGE IN YEIGHT VIC 111.939734 -- CHANGE IN GAMMA R VIC 2.238795 FtNAL STEP SILE COEFFICIENT 1.912763 OUSC FOR NEXT TRAJECTORY 2039.651367 . ..- ._.. - .-- -- -- NEXT TRAJECTCtiv.uILL RE A -V.ATfb-$-rcF -.- ---..--.-___ ----___ --- - -___--.-- .____-- _.._.- -__.

Tl !!E ‘dELGHL.-PLLLNDE GAUUA- YELQC LT.Y .RTEllSlET JR&G----- 0 YACH NUM8ER CE IGHT FLOW 4LPHA PITCH ANGLE 161.1 .p4c 5oc..~co~ -.~-..--~o~.wo. ___ 6’TO.OQO -..-- -----..62.234- _.. -.- 4921.T_zp_-..--~6674.~63.1.-- ---.- .- --49’..Cc ___._ --.- -. - -___. 421.757 0.538 L5.000 10.664 72.903 0.500 1599.440 me73 62.E71 629.861 4942.210 553.R34 420.293 461.069 0.566 25.000 a.Doa 70.879 ~.-2_5_0--~-1580,690 .-.___ -L&%892 -y’-fK&--A;.~;~ 4974.67@ 323.963 _34wdk-- .--_-__-_- 526.418 0.609 25.000 . .

2.000 1561.940 1679.385 62.8C7 727,192 5 no2.m 7 35-746 344.409 598.354 0.655 25.CtO 0.398 63.205 2./50 1543.150 2179-021 -- -_ -61.634 -- 777.954 - _.-.- .-. 5045.171 - -_____ - 3 3je,13a !! ?.a,se4__------ 674.735 0.702 25.COO -3.341 58.293 3.500 1524.440 2704.175 59.632 827.638 5082.532 -817.511 673.505 751.867 0.748 25 .OOO -7.249 52.384 4.625 1496.315 3520.164 54.843 893.304 5139.474 -1776.120 1449.194 -- _.-. -..-_ -. --- --.---_ --..------ - 854.745 0.810 25.OOC -13.853 40.990 5.750 1468.190 4324.671 47.413 932.396 5194.183 -2830.61A 2772.875 908.956 0.848 25 .OOO -20.761 26.652 5.972 1462.t34 4416.731 46.567 946.030 5204.382 829.263 77~.940 - 931.452 0.861 25.CCO 5.935 52-502 6.194 1457.075 4630.588 46.877 961.043 5214.681 1277.296 1032.216 956.783 0.075 25.COC 8.900 55.777 6.417 1451.523 4788.485 47.177 976.188 5225.147 1212.6R9 lCO~.I5@ -- 982.488 0.889 25.OOC a.229 55.4C16- 6.750 1443.150 5030.155 47.560 999.446 5241.149 1115.758 961.77@ 1022.336 0.911 25.00C 7.276 54.836 7.250 1430.690 5407.516 47.989 1035.413 5265.801 :67.?47 922.927 - 1084.772 0.945 25.CCO 5.944 53.934 8.000 1411.940 6001.e56 48.325 1091.357 5304.118 732.86P 9c1.957 1183.553 0.998 25.000 4.129 52.453 9.125 1383.815 6555.196 48.199 1178.465 5364.1‘42 359.691 946.9C2 1340.342 1.082 25.000 1.789 49.978 10.250 1355.690 7974.733 47.350 1267.674 5426.412 -34.541 ll169.472 1502.966 1.168 25.000 -0.153 47.197 11.375 1327.565 5C48.618 45.892 1357.352 5489.906 -421.016 1251.19P --- 1666.614 1.255 25.000 -1.684 44.208 DRAG ALTlTUDE __T!Y_E.-..-- -..!!E IW-.---. .-.M!!MA -R---. VEL?CITV.R -_ ._ JHII’JST ____- CI.FT_-.---- -. ___ ---.-- __-..

ALPHA PITCH ANGLE cl MACH XUMBER dE ItiHT Fcau 42.836 1490.702 5585.130 -851.944 l5b3.3Ek 13.062 1295,371 lC730.372 - 1906.894 I. 387 25.000 -2.978 39.857 39.553 1624.168 5677.865 -699.156 1771.316 - 14.750- 1243.150 12458.354 - -.--- - _--- 2142.732 1.521 25.COO -2.177 37.376 2012.363 16.437 1201.002 14236.133 38.003 1762.885 5768.072 lPOl.489 2383.9AO 1.662 25 .ooo 2.801 40.803 17.187 1182.252 1501o.L46 38.913 1 tlL7.-732 5808.590 2812.498 2775.167 -- ------ - -. .--__ 2466.939 1.720 25.COO 7.601 46.514 17.937 1163.502 15954.615 40.828 lRb0.691 5850.402 7779.515 2778.078 2511.187 1.766 25 .OC'O 7.379 48.207 2425.586 2bL6.!l,n5 18.437 LLcil.o>Z 165:s 791 41.060 1893.591 5878.99R .--- 2548.264 1. RO2 25 .OOO .6,3~i--.- 48.207 .

2(,02.734 2472.697 19.187 1132.252 i7547.8se 43.093 1947.788 5922.707 5.114 48.207 2610.777 1.860 25.000 1707.145 24PB.R75 19.860 1115.440 18462.139 43.944 2000.000 5962.500 2669.995 1.917 25.000 4.263 40.207 z ITERATION 1 TRIAL 2 NETA= 46 LETA= 73 MAJORITY VOTE TEST -I -__.. _--. - .--- .____ _--. _- --..--__ TOO NANY FUNCTIONS YITH ADVERSE TRAVEL 0.17390345E 03 OUSC REWIRED FOR OBETA -_NUnERCl_~~F_F_IR_ST-.LI1GR~GECULIlL TER Pd34cQq39E. 03 __--. ._ --.- ___ 0.11635616E 03 -0.39952482&01 LAGRANGE MlJLT IPL I ERS ACCURACY CHECK FOR OUSO 0.50991284E 03 HEIGHT ALTITUDE CURRENT END POINT VALUES 18462.139404 1115.440140 PREVIOUS VALID STEP EN0 POIhT VALUES ---- -. -..____. 10h6.3~.6520m ta19s. 317139 -___-.-...---- ---.. -_-- _..

CHANGE IN EN0 POINTS 263.022266 49.053619 CONSTRAINT TRAVEL INOICATCRISI -1 CONSTRAINT CROSS-OVER INOICATGRISI ALLOWABLE FORWARD CCNSTRAINT ChbhCElN-01 0.000000 ALLOHAELE BACKYARD CONSTRAINT CHAhGElN-Cl 0. ocoooa CONSTRAINT CHANGE ASKED FOR(N-C) -.. O.WOOO. _- .__ CONSTRAINT TOLERANCE{ Sl loco .cLoooo NON-LINEARITIES o.oocoon -0.156841 STEP-SIZE COEFFICIENTS C.@OCCOO 1.912763 -.-_-.

MAXIYUH STEP-SIZE COEFFIClEhT(S1 0.000009 w PRECICTED EN0 POINT CHANGE(S) (TGTAL) -154.957928 43.896452 -721.776566 PREOTCTEO EN0 POINT CHANGE(S) (DLS.0) 4.134153 -~-__. . _-- =:- --..-.-- __.-- - PRECICTED END POINT CHANGEIS) (VICI 566.818642 p-jq.-762i95 WE1 GHT VIC 55.969867 CHANGE IN CHANGE IN GAUUA R VIC 1.119397 --- --- FINAL STEP SIZE COEFFICIENT o.5ouooo LIUSO FOR NEXT TRAJECTORY 5C9.912R42 NEXT TRAJECTCRY HILL BE A VALIC STEP - _--._- GAYMA R THRUST DRAG TIME IL_[!TUOE VELOCITY R LIFT - ---- . .--~ .____--------___ -. .'cC IGHT ..__ ----cl----- HACH NUM8ER UElGHT FLCU ALPHA PLTCH ANGLE -c!.ooc 1555.970 5co.cco 51.119 500. COO 4921.729 424.471 336.05lJ -- 421.757 0.538 25.000 6.710 57.829 0.500 1543.470 770.210 61.290 533.312 4942.011 378.055 331.821 -- ~ -- -- 456.171 0.559 ZS.hti-0 5.4c7 66.696 1.250 1524.720 1203.t42 51.264 504.577 4974.209 280.574 343.Gll 537.785 0.515 25.olJc 3.470 64.742 2.000 1505.97c 1670.359 50.882 737.245 5008.420 142.118 363.235 ____-- - -- ..-- _ 515.17a 0.664 25.COO 1.540 52.422 3.125 1477.045 2430.178 59.581 819.194 5053.112 -151.037 451.064 740.791 0.739 25.COO -1.449 58.132 4.250 1449.720 3254.055 57.252 R99.559 5121.03@ -63C.281 653.142 --- 073.922 0.815 25.000 -4.808 52.444 5.315 1421.595 4124.881 53.592 577.405 5180.710 -1347.066 1091.385 1004.851 0.888 25.000 -8.937 44.55- 5.708 1413.262 4387.550 52.182 958.503 5198.4C7 -1531.839 1307.993 1040.442 0.9ca . G -10.455 .

6.042 1404.528 4551. EC5 51.555 1022.555 521b.OCB 934.577 903.552 1082.712 0.931 25.CCO 5.755 57.310 5.264 1399.373 4831.407 51.710 1038.944 5227.991 ea7.C3fl 9C2.@5(' -_ - 1111.425 0.946 25.500 . 1 5mi30 5.597 1391.c35 5105.a21 51.897 1063.654 5245.149 R16.657 905.229 1155.270 0.970 25.coo 4.713 55.510 7.097 1378.539 5533.209 52.084 L101.318 5273.969 709.754 919.939 ~2-zrz3 7 r----- 1.006 25.ooQ . 0 55.954 7.047 1359.789 52tl2.458 52.169 115R.852 5316.895 547.583 962.273 1325.350 1.051 25.COO 2.752 54.921 0.972 1331.564 7269.720 51.895 1246.938 5383.562 299.734 1063.712 -__ 1485.182 1.145 25.000 1.345 53.241 10.550 1209.477 9001.749 50.712 1381.293 5407.179 -72.ObC 128O.t?92 1728.448 1.277 25.000 -0.278 50.435 12.347 1247.289 1caba.727 40.705 1515.505 5592.759 -392.434 1545.421 -__-.---.-___ - -- 1954.861 1.412 25.000 -1.332 47.453 14.035 1205.1c2 12843.152 46.455 1552.085 5597.825 -511.069 1781.726 2189.885 . . - m 4.9 _ _. ------ .._--- __-_-. -___- -.-----__ TI?fE T.*nw* R VFI IlClTY R THPlfCT 1 IFT .WlCHTF RB--- a MACH NUMRER WEIGHT FLCh ALPHA PITCH ANGLE --.- _ .-. 15.??2..e.e 1162,914~_1491P.199._~_~~..44..447.__1791.8Z4___1BILO.B99 --113,222---__-l.yzq.129...__--.~..

2409.748 1.694 25.COO -0.313 44.129 17.410 1120.727 17=-265 44.322 1933.tl14 5902.731 1521.8Q3 2311.3Of1 2612.022 1.844 25.COC 3.384 48.207 4.4L9.94 ?.OOO .O-OQ _ 5951.93R --1&d 1 B- UO.O~.52?_-.. lAzl6_tl4L. ---12S7,2-3_C 23@9.tibP ______ ___..- .^_. _ 2691.904 1.915 25.000 3.213 48.207 ---- -- ITERATION 1 TRIAL 3 NETA= 36 LETA- 73 --.---- ------ MAJORITY VOTE TEST 0 - THETA TABLE 11 LINES -0.000000 67.820964 0.25cc00 67.281850 0.500000 66.696447 c.e750(‘0 65.76C060 1.250000 b4.?417C3 1.625000 63.633300 2.OO~I'300 62.422361 2.5625CO 60.429440 3.125o_qp 58.131>@ 3.@7500 55.520065 4.25Oi)l'O -.-- _- _- 52.444272 5.708333 41.726336 4.812500 49.851157 5.e75cco 49.142176 5.375000 44.655149 6.041667 57.309796 5.541667 43.236309 6.152778 57.170612 6.263889 57.0303t5 6.430555 56.819995 5.5972?? 56.609626 b.047222 56.203332 7.C97222 55.954317 7.472222 55.443790 7.A47222 54.920916 0.409722 54.1coa31 0.972222 53.240920 9.815972 51.e75994 10.659722 So.434552 ll.E03471 48.9_4.3_4_18 12.347272 47.423438 13.190_971 46.046348 14.034721 44.895663 14.878471 44.153043 15.722221 44.1291Pl 18.217864 4a.tot633 182.17EC42 40.206633 RESTART TABLE - 0 1 509.912842 0.05c000 -o.coooco 17.344538 1066.386520 1l3198.317139 1000.000000 0.050000 TRAJECTCRY SUMMARY *+*+*+*****+++*********************************** _-__.-..--_--__.- ._._ -- TIME YEIGHT ALTITUDE GAMYA R VELOCITY R THRUST LTFT DRAG 0 MACH NUM9ER UEIGHT FLOY ALPHA PITCH ANGLE 17.345 1066.387 18198.317 48.207 2000.000 5951.141 0.000 2154.912 2693.640 1.915 25.000 0.000 48.207 18.218 1100.523 13216.741 44.994 2000.000 5951.938 1297.23C 2309.Abb 2691.904 1.915 25.000 3.213 48.207 -AZ.------- --- Operating Information Program setup. -The program is divided into nine logical blocks or links, where each link performs a particular job. These links are required to ensure that the program will fit into a 32K core computer. Additional information is given in the section on program and data overlay.

Link 0: Controls the flow of the program from link to link and remains in The link consists of the following subroutines: core at all times.

LOAD LOOKlD ATMOS UNITZ BLOCK DATA Link 1: All operations related to the forward trajectory are performed.

The subroutines included in the link are: AKSTP EXEC PLAC ANPARP FPROG STEP1 ANPRTL INITCO STPl BOOM LOOK3D TLD CONTR Link 2: The numerical partial derivative check is made using the subroutine: PRTIAL Link 3: The control variable history for the nominal trajectory generation is made using the subroutine: GUIDE Link 4: The operations related to the automatic convergence logic are per- formed. The subroutines included in the link are: VALID KCALC MATRXB MATOUT UCALC VARIC MATRIX Link 5: All operations and controls related to the backward trajectory are performed. The link consists of the following subroutines: DVALB STEP2 LAMBDA STP2 Link 6: This link performs the initialization that is required once per data case. The subroutines included are: INITAL TITLES Link 7: The final output of trajectory data that is required only once per data case is performed. The subroutine called is CARDS Link 8: The last link controls the plotting of data and is called once per run, The subroutines used are PLOTZ SKALZ A set of auxiliary subroutines are required for plotting similar to the Boeing numerical plotting system (NPS).

Organization of links: The organization of the computer links and the approxi- mate core storage required are shown in figure 17. This chart shows that all links of origin one are called from the basic link 0. Links 2 and 3 (origin two) are called only by the forward trajectory, link 1.

- Stop Subroutines ---- Labeled Common Link 0 System Subroutines Link I (8500) Link Origin Two Link 4 I 1 I I W-3 Link 2 Link 3 @W (I 200)

Figure 17. PROGRAMORGANIZATION

Physically in the computer, links 1, 4, 5, 6, 7, and 8 occupy the space beginning with the same storage location (ORIGIN ONE) but at different times as Similarly links 2 and 3 occupy the stor- governed by the main program (LOAD).

age beginning with the same location (ORIGIN TWO).

Program deck: The recommended program deck setup is shown in figure 18.

The control cards required for machine operation precede link 0. The data deck for a given problem is placed following link 8.

To speed up the overlaying of links, two additional tape (or disk) units are used to store the program, The forward integration routines (links 1, 2, and 3) are stored on SYSLB3. The convergence logic routines (link 4) are stored on SYSC El. All other links are stored on SYSUT2.

Data setup. -The setup of the data is performed as detailed in “Input Data Preparation. ” The physical setup of the data deck is shown in figure 19. Note that each set of stage-dependent data must be preceded by an STG card. A STG card must be present for a stage even if no stage-dependent tables are input.

The last card in the data deck is an END card.

General Machine Operation Tape or disk requirements. -STOP may use a maximum of 10 I/O units during execution, excluding the system units reserved for input, output, punch, and normal link loading (i. e. , SYNINl, SYSOUl, SYSPPl, SYSUTB) 0 Listed in figure 20 is a brief description of all I/O units used-their prograr symbol, system name, function, mode, and buffer size. (Buffer sizes are set in subroutine UNITZ.)

Figure 21 shows the I/O units used during execution of a typical data case.

This data case is set to complete 1 iteration (NARBY = 2, NITC = l), punch the control table used for the last iteration (NC(8) = l), and plot the nominal and last iteration (NC (ll) = -1).

End of run indication. - The program will normally exit successfully by reading an end of file on the input tape (SYSINl). The monitor will print out the comment “END OF FILE READING. ” The program has the capability to run multiple data cases, and always terminates one case by trying to read the next.

[(FOR PLOT LINK) A I/ f $ INCLUDE LINK 7 INCLUDE FOVERF 11 P .

ORIGIN ONE, SYSCKI., REW LINK 3

Figure 18. PROGRAM SETUP

WEIGHTING TABLE RESTART TABLE NOMINAL GUIDANCE CONSTRAINT DEP. PARAM. / STAGE DEP. PARAMETERS INITIAL CONDITION DATA / NC CONTROLS CONTROL CARD TITLE CARD A

Figure 19. DATA DECK SETUP

Program Buffer Size System Symbol Name Function Mode (decimal) Remarks KINP SYSINl Input BCD Printed output 120 KOUT SYSOUl BCD System Units KPUN SYSPPl Punch output 22 BCD SYSUTB Link storage BIN Forward inte- SYSLB3 BIN gration link Link storage Matrix link SYSC Kl BIN KDAT Input data 256 SYSUT6 BIN storage KLAM SYSUTB Impulse response BIN 256 storage KPAR Partial storage 256 SYSUT4 BIN Scratch units required by STOP KPLT Plot data storage 20 SYSUTl BIN *KSCR Control table 20 SYSUT9 BIN storage KTAN SYSUT5 Control table BIN 256 overlay Scratch tape for 256 SYSUT7 BIN the NPS gener- Used only when ated code plotting Plot output *sysc K2 BCD 22 i * These units must be assigned to TAPE, as they may be saved

Figure 20. INPUT/OUTPUT (l/O) USAGE

Special-machine operating information. - The two I/O units, KSCR and SYSCK2, must be assigned as physical tape units and comments inserted asking that these tapes be saved.

Link I/O Unit I/O Unit Routines Read Written Purpose Loaded From Number CORE Initialize data SYSUTB 6 KINP KOUT KDAT KTAN 1 Forward integration of SYSLB3 KDAT KOUT nominal KPAR KSCR KPLT 4 Majority vote test SYSC Kl KOUT SYSUT2 5 KPAR KLAM Backward integration KOUT KDAT 4 1. Steepest ascent logic SYSC Kl KLAM KTAN 2. Build new control table KOUT KPAR 1 KOUT Forward integration of try 1 SYSLB3 KDAT KTAN 4 Analyze try 1, make another SYSC Kl KLAM KTAN try or run valid step KPAR KOUT KPAR SYSLBS 1 KDAT KOUT Forward integration of itera- KTAN KPAR tion 1 KSCR KPLT SYSUTB 7 KSCR KOUT Print and punch final output KDAT KPUN 8 KPLT KOUT Prepare plot tape SYSUT2 SYSUT7 SYSCKB Exits trying to read next SYSUTB 6 KINP data case Note: KTAN will not be used if the control table used does not exceed 1000 points.

Figure 21. EXAMPLEOF I10 USAGE

The following control cards are required for each computer run, and must precede the program deck.

SJOB 99F-61C AS2413 WATSON BCC A2 5-0934 s* ACCTTIME 5 SB PRINT ESTIMATE (2000) _ _.--.-- . . --_ s+ S+ READY A4.A5rB3rB5rCZrC3rC4 s* SATTACH A4 SYSUT7 SAS SATTACH .-. _-_ --._--. gsrrs2.- ..__.. .- ..- ..__ _ ‘- SAS SATTACH 83 SYSUT9 SAS SATTACH 85 SAS SYSLB3 SATTACH c2 sAs . . .-.. .-spsuT5-.- .-.._. _-_. -...- ..-.. -._. . ..^ SATTACH c3 SYSCKl SAS SATTACH C4 SYSUT6 SAS SALT NPS .-IBmB. .-.-.. -.-- - -.. -f- .._ .- - .- S-EXECUTE .- -.- SIBJOB GO,LOCIC~MAP~FILES~FIOCS,EXTLIB SPOOL -UNIT03-r-UNIT04-r-UNITO9-~-UN~TlO-~-UNITll-~BLK=256 SETC BUFCT=4 SSROUP -UNIT03-r-UNIT04-,-UNITO9-~-UNITlO-~-UNITll-~OPNCT=4 SETC BUFCT=4 -- SPOOL- ..--- ~;UNITOB-i’UNlTT4-,6LK=~~BU~=2 ---.- - -.-’ - - - ---- - -.

IGROUP -UNITO8-r-UNIT14-rOPNCT=2rBUFCT=2 NOTE: (1) When NC(8) # 0, a comment card must be included with the control cards, asking that a save tape be mounted on B3.

When plotting (NC( 11) # 0) , a comment card must be included with (2) the control cards stating that unit Cl will contain the plot output.

Programming Information Basic program flow. -The following flow charts are intended to give the - user a general picture of the basic program organization and the subroutines The flow charts are broken down into called by each major area of the program.

the logical areas of forward integration, reverse integration, and control logic.

-. - RUN A NOMINAL BASIC (FIRST VALID STEP) OVLKALL TERMINATE If FLOW (EXEC) NOM NAL ONLY IS DESIRED PUNCH CARDS IF REQUIRED (CARDS) EXAMINE ABORT RlJN NOMINAL IF ERKOR b IN NOMINAL (MATRIX) PUNCH CARDS IF REQUIRED (CARDS) RUN A TRIAL (FORWARD INTEGRATION) (EXEC) RlJN A VALID STEP (FORWARD INTEGRATION) (tXEC) EXAMINE VALID (CARDS) VALID STEP BASIC FORWARD INTEGRATION FLOW ENTER FORWARD INTEGRATION LEAVE FORWARD INTEGRATION CALL FPROG

I

I STPl 1.

I

J BASIC REVERSE INTEGRATION FLOW ENTER REVERSE INTEGRATION

l

CALL l I DVAL2 CALL LAMBDA I F I , CALL STP2 A LEAVE REVERSE 2 INTEGRATION I tr - CALL STEP2 z

I

z oi

b

INTEGRATE n SF A STAGE Z WV) UJ w - I CALL STP2 I BASIC CONTROL LOGIC FLOW ENTER CONTROL LOGIC

I I

REVERSE INTEGRATION VALID STEP bi ;;::, 1 CALL , MATRIX 1 TRIAL

r

CALL VARIC LEAVE CONTROL LOGIC CALL - KCALC v I CALL UCALC .

CALL MATOUT Subroutine descriptions and flow diagrams.-Each of the subroutines in STOP is described giving the purpose and the numerical method if applicable.

The approximate Other subroutines called by the one being described are listed.

core storage used is indicated. The arguments of subroutines using a calling sequence are defined in detail. Flow charts showing the organization and detail of each routine are given.

AKSTP- AK Store Purpose AKSTP defines and calculates basic and frequently used variables in the These variables include some AK’s, VAR’s, partials, and trigonome- program.

tric functions of angle of attack, latitude, flight path angle, and heading.

Subroutines Called ATMOS BOOM LOOK 1D TLD Storage Used 952 cells AKSTP Calculate XKR CK3 SK3 NO NO Colculote Colculote XKR XKR CK5 CK5 NO Calculate XKR CK6 SK6 AKSTP Calculate AK (5) t NO Calculate VAR (4) PQPH r- VAR (5) PQPV ’ !r Calculate AK (I\) AK (21) AK11 SAKll ,CAKl’l Colculato AKSTP -Analytical Partial Derivatives of Placard Equations of Motion ANPARP Purpose The analytical partial derivatives of the placard equations of motion (part of the F array) are performed in ANPARP. These derivatives of F (placard) are with respect to all the state and control variables. See appendix B for an algebraic description of these partials.

Method Derivatives of the placard equations with respect to the state variables are

1)

given as PFX(1, J) where I identifies the placard equations, and J identifies the derivative state variable; i. e. , PFX(K25, K4) = aF(K25)/dX(K4).

Derivatives of the placard equation with respect to the control variables are

2)

given as PFU(I, J) where I identifies the placard equation of the F array, and J the derivative control variable; i. e. , PFU(K25,Il) = aF(K25)/aU(Il).

All partials not calculated are set to zero.

3)

Words not in COMMON computed in ANPRTL and required by ANPARP are

4)

transmitted through the calling sequence: CALL ANPARP (PLPH, PLPV, PLPA, PLPS, PDPH, PDPV, PDPA, PDPS, PTPH, PTPV, PTPA, PTPT) where the elements in the call are defined in appendix B.

Subroutines Called BOOM Storage Used 1,371 cells Evaluate Pdialr A

n

ANPARP I Evaluate Partials I

‘r

Of Dynamic Pressure Placard F(t)

I t--

Evaluate Partials Of Qa d Placard F(M) t- Evaluate Partialr Of Resultant c Physiological Acceleration I- 1 Placard F(i) I Evaluate Partialr ANPARP Evalwte Partials Of Sonic Boom * Overpressure PlaFard prnerical Perturbation I 1 I Evaluate Sonic Boom Partials With Respect To Velocity

f-l

YES E i Evoluoto Partialr - CIfcI;h Numkr I W-0 Avai lablo w For + Growth Evaluate Partials + w of Gamma Placard W, W Avai lablo t * For Growth Avoilablo w * For Growth - ANPRTL -Analytical Partial Derivatives of the Equations of Motion Purpose The analytical partial derivatives of the equations of motion (excluding pla- card equations) are calculated in ANPRTL. These derivatives of the F array are with respect to all the state and control variables. Appendix B gives an ANPRTL is called from STPl only algebraic description of these partials.

during valid trajectories.

Method Derivatives of the equations of motion with respect to the state variables

1)

are given as PFX(I, J) where I identifies the equation of motion and J iden- tifies the derivative state variable; i. e. , PFX(K2,K4) = aF(K2)/aX(K4).

Derivatives of the equations of motion with respect to the control variables

2)

are given as PFU(I, J) where I identifies the equation of motion and J iden- tifies the derivative control variable; i. e. , PFU(K2,Il) = aF(K2)/aU(Il).

All partials not calculated are set to zero.

3)

Subroutines Called ANPARP Storage Used 1,649 cells ANPRTL

t

Compute

PLPH

QS

PQSPH PLPV PLPA PQ SPV t . 1 _ Evaluate Partials r Of Weight - Evaluate Partials Of Altitude Evaluate _ Partials Of # Relative Flight Path Angle ANPRTL YES b NO

I

/NO

/NO ANPRTL / * I I Evaluate Partials -1 Of Gravity I-L Loss I E ANPRTL Evaluate Partials : Of Thrust Vectoring Loss ANPRTL Placards > > Return ~1 Return NO T Increment Numerical Partial Calculation Indicator By I F( Return ATMOS -Atmosphere Calculations Purpose ATMOS calculates atmospheric pressure, temperature, density, and speed of sound. Any atmosphere may be used; however, the calling sequence must be the same as presently used.

Method The ATMOS subroutine normally used in the program is an analytical repre- sentation of the 1962 ARDC standard atmosphere.

The call to ATMOS is made with input and output variables transmitted through the calling sequence, i. e. , CALL ATMOS (H, VA, P, D, T) where H is the altitude, ft VA is the speed of sound, fps P is the pressure, psf D is the density, slugs/ft3 T is the temperature, deg Rankine Limitation This ATMOS subroutine is a modified version of the complete 1962 ARDC standard atmosphere, Because of a limiting altitude of 157,000 feet and other simplifications, the modified ATMOS is considerably faster than the complete atmosphere with little or no loss of accuracy within the bounds of common application. Other ATMOS subroutines are presented in reference 9.

Storage Used 250 cells ATMOS Convert Altitude To Meters

-I””

I’TREF=216.65 (OK)) .eater Than Or Equal y’A And > n 7n KM Alti tide, Al2 tae Greater Than Or E ual <20 KM. And ~~::I%-% ::d No Less Than 48 KM Less Than 32 KM.

YES

Y

T RE DT/un= ,001 (-K/M) Return HREF = 2OCW. (M) lS.l285+(PSF) IP RFF = , _ PREF = 114.346 (PSF) ’ “L’

I

I

I L Calculate Density = 0 o TM Pressure = 0.

CS (Velocity) Temp. = 489.676 P (Pressure) CS = 1084.796 t- / Return

I

Program Constants and Output Titles BLOCK - Purpose This subroutine is provided for input of numerical constants and alpha- numeric printout titles, which are required by the program and are not case These variables are set once for each computer run.

dependent.

Method All constants and titles are input by data statements.

BLOCK

1111-1-11-----

-I

I-

t

I

I Mode Indicators I I

I Accomplished

I At Loah I Time Only Set Tape Unit I Assignments !

I Set Program !

Constants !

I

Initialize

t

Arrays (DESXUK, XNAME,YNAME

I

IGUIDN ,NAM) I

I

I

------II-

a---

-Sonic Boom Overpressure and Placard Determination BOOM Purpose BOOM calculates the sonic boom overpressure and determines the equation of motion for the overpressure violation in the path of the aircraft.

Method Sonic boom overpressure on the ground is determined for the conditions along the flight path. The sonic boom penalty function is calculated by the method discussed in the analytical development section on enroute constraints. The terminal value of the penalty function (auxiliary state variable) must be con- strained to zero. Location of the shock signature on the ground is determined as a function of the longitude, latitude, and aircraft flight conditions.

Assumptions and Limitations The overpressure violation is calculated at an odd number of points so that both maximum lateral cutoff points and the midpoint are included.

The number of points is selected by setting NC(75) = number of points. If NC(75) is even, the number of points calculated is NC(75) + 1. If too few points are selected, a region where the placard should be violated may slip through the mesh without being detected.

Subroutines Called LOOK3 D PLAC Storage Used 298 cells BOOM NO + Return ..:

,”

P I

BOOM Compute Sonic Boom Placard * “olation v i’ I /

I

Return

I

Card Output Routine CARDS- Purpose CARDS is called at the end of each data case and is assigned the following tasks : 1) Recover control table from last valid trajectory from KSCR unit and print and punch control table and restart table; 2) Recover trajectory summary from KDAT unit and print trajectory summary table.

Method The control variable table is recovered by backspacing the KSCR unit to the start of the last valid trajectory and reading the control history into core. The restart table is built from convergence information in COMMON. The control table is then punched in octal format (in order to maintain sufficient accuracy) and the restart table is punched in decimal format (for ease m modification if required).

The trajectory summary is obtained from the KDAT unit which was written by LAMBDA at the end of each valid iteration.

Storage Used 681 cells CARDS Backspace KSCR Unit Read Control Variable Table From KSCR Unit + Calculate Last Point Of Control Variable Table c Print Control Variable Table CARDS

Bo

I 1

Read Trajectory Summaries From KDAT unit Set Indicator L(6) = 10 I Return

f-l

CONTR - Open Loop Control Variable Calculation Purpose CONTR performs P table lookup for the control vector (u) in the TlMEU array at each new time point as required by the integration package.

Method The control variables obtained are specified by the setting of NC’s 7’7 through 80 for 8, cp, 7, and A respectively. An NC set to 0 indicates the control vari- able is not required, set to 1 indicates the control variable is used. Values of the variables as functions of time are obtained by linear interpolation from the control variable table which is either input nr generated by the program. For the nominal trajectory, if NC(l3) = 0, the control variables are calculated as described above; but if NC(l3) # 0, the control variables 8 or r) are calculated in the GUIDE subroutine and override the value looked up by CONTR.

CONTR is called prior to GUIDE.

Assumptions and Limitations The control variables are obtained by a linear table look-up. The table look-up package does not extrapolate so the user must input data that will not permit a table overrun.

Subroutines Called GUIDE Storage Used 197 cells Set TIMEU Arra Index CONTR For lurrent Lookup * Set Control Variables Lookup Times Direct1 From InTlMEU Array TlMEUkray On Either Side Of Current Time Interpolate For Values Of Con tro I Variables save TIMEU Array Index For Next Lookup Calculate UI I , A Selected, Cl 1 Control Variable s I1 I Calculate Angle Is Bank YES UI2 A Selected < Cl 2 Control Variable SI2 NO I CONTR , YES z; = U(U) * NO I c Return + Roturn YES& Return

r-l

- Adjoint Variable Derivative Evaluation DVAL2 Purpose The derivatives of the adjoint variables and the integrated payoff sensitivities are evaluated for the reverse integration.

Method The differential equations for the adjoint variables are evaluated as discussed The. partial derivative PFX or PFU in the section on the steepest-ascent method.

and the control variables u were stored as functions of time during the forward trajectory. Because this subroutine is important for the basic understanding of the computer evaluations of the adjoint variables, weighting matrices, and I matrix, the computer language for important variables will be given.

The array of all derivatives evaluated during the backward integration, the FB array, is constructed with elements for A given by where F C PFX matrix A E XB array The derivatives for the elements of the I matrix are dI = XXDOT array dt = h’ G W-l G’ A (from equation 109) where A E XXLAM matrix (given by equation 110) G 9 PFU matrix W-1 E VINP array The derivatives of the elements for the automatic weighting are evaluated by a: - = s f (from equations 98 and 134) I I dt where s: (t) - G’ XQa (equation 130) = PFU (I, J)* XXLAM (I, 1) The elements of the UULAM matrix are the weighted impulse responses of the constraints and performance to the control variables as given by W-1 G’ A.

The call to DVAL is made through the calling statement.

CALL DVAL (UULAM) Note: When a single-dimension array is used to store the elements of a matrix, the elements of the matrix are stored by rows.

Subroutines Called LOOKlD Storage Used 1229 cells DVALZ NO Call Compute Elements LOOK10 of FB

P

Compute UUIAM’S Load Values I Into PFX Compute PFU XXDOT’S U And Store In FB NO

Y

----l- Compute Weighting Matrix Elements Of FB Calculate PFX, PFU, U Return By Inter olation

A

In TPU I; rray i - Forward Integration Flow EXEC Purpose EXEC is called at the start of every forward trajectory.

The basic purpose is to control the program flow in performing the forward trajectory.

To accomplish this, the following tasks are accomplished.

Read nonstage-dependent data from KDAT; 1) Initialize state variable array; 2) Read stage-dependent data from KDAT; 3) Position KDAT past last written trajectory summary at start of last stage 4) only; Initialize the state variable derivative array; 5) Integrate to end of stage or XSTP for last stage. Control data overlay is 6) accomplished as required by reading blocks of data from KTAN.

Method Integration of the forward trajectory utilizes the variable step Runge-Kutta option of STEPl. The final step to XSTP is iterated using the fixed step option until X(NST) is within EPSLN of XSTP.

Subroutines Called AKSTP FPROG STEP1 CONTR INITCO STPl Storage Used 665 cells EXEC (Plot Indicator) The Last Iteration NO x EXEC Read Non-Stage Dependent Table Data From KDAT Unit + Set NETA =0 L(2) = I NSTAGE =l LINES = 100 L(8) = 0 T = VAR (8) I YES set Liner = 100 I Set INST=5000 INPO = 5000

+

Read Stage Dependent Tables (If Any) For Appropriate Stage From KDAT Unit + EXEC f I 1 Summaries On EXEC Return Set Integration Mode TO 4

r-l

L I I YES

AC

ation lntegr Of First Stage Integration Of A Stage Value: Of Time And Derivatives Of EXEC Return Load Into Available YES _ . Core Storage (TIMEU Array) Control Table Dota On Data Needed KTAN Unit L I ‘r I Set INST = 5ooO INPO = 5000 T=TSTP(NSTAGE: Return

LJ

EXEC 10 Std p6a dn Greats;“ I %la n Or Equal Tar Sto&ing Co/rcKiion Lo8t Stage YE To Stopping Condition set I 1 INST = 5000 INPO = !%OO + L(8) = 0 +=-l Set Integration Indicator For lnitialirati~n be Plot c < telmct, Return

I I

EXEC F

Q

‘c

Backspace KPAR Unit Decrement The Number Of Stored Integration Points By 1 Reset Values Of State Variables And Derivatives Of State Variables I Equations of Motion (Calculation of F Array) FPROG - Purpose FPROG calculates the entire F array, which consists of the equations defining See the derivatives of the state variables and placards with respect to time.

appendix B for an algebraic definition of the F array.

Method 1) All state variable derivatives with respect to time are calculated directly in FPROG. The placard derivatives, however, are calculated in the PLAC sub- routine called by FPROG.

The value of the placard required in the evaluation of the placard derivatives 2) is obtained by a table look-up.

The array is initially zeroed so that only the calculations performed in a 3) given case may have nonzero values.

Note: The sonic boom placard, F(K35) is not determined in FPROG but by the subroutine AKSTP (which calls BOOM which, in turn, calls PLAC).

Subroutines Called PLAC Storage Used 677 cells - FPROG Be Evalu ted YES Set Index To Evaluate Al I Desired State Variables

L-

‘I

I

I

Is Th- FPROG

I 1

, Is Lati tde \ YES FPROG Evaluate Longitude Angle Equation set Dummy Time Evoluote _ Path Range Along.

Earth’s Surfuce Equation Evaluate L Ideal Relative w Delta Velocity Equotion G FPROG F c, / Evaluate ?

1 Evaluate 1 Evoluote i I Wing Sweep Plocord W, M) Avoiloble

I I

-7 :;ow th /---B

I Evaluate I

Evaluate

I 1

_ Dynamic Pressure Plocord F(M) FPROG J Eva bate Total Temperature - Placard F(t) Evaluate Normal Load _- Factor Placard - W-f, M Evaluate Altitude Placard F(M) I I

I 1

Evaluate Mach Number Placard W-Q FPROG

I Evaluate I

Gamma b / Placard I= 0-W) , 1 Available 1

*

-_I FcSowrh

Available For Growth NO NO Return Return

I I

- Closed Loop Guidance Control Variables GUIDE Purpose GUIDE calculates control variables for generating a nominal flight path.

Method The NC(l3) controls the manner in which nominal trajectory is generated.

The GUIDE subroutine is called only for the nominal trajectory, if NC(l3) # 0.

If the first stage is a tilt maneuver, GUIDE is bypassed. The number of phases used to generate the nominal is given by NC(l3). Corresponding to each of the phases, a guidance mode and phase stopping condition must be specified. The provision is made to allow the phase stopping condition to be approached from below or above. The guidance modes are discussed under “Nominal Trajectory Generation. ” Ass umpt ions and Limitations A maximum of 12 phases may be used to generate the nominal flight path.

Five tables have been reserved for guidance modes. The equations solved for the control variables are linearized as described under “Nominal Trajectory Generation. ” Subroutines Called AKSTP LOOKlD Storage Used 657 cells GUIDE kder Indicator Return NO YES I^_I a&h rrmrurn I Transfer To Selected Mode Calculote 4. Return GUldE

b

Set

ALF = Ull I Calculate + UI1 Gamma Is Function Of Time Mode

I

I

Calculate @j-- GDOT It--+@ ’ I I Gamma Dot 5 Is A function Of Time Mode GUIDE Calculate XMV GF Calculate QS i Calculate ALF IV Maximum r\l ,owoble Value

ix

p

-- R

GUIDE < DHDV Equal >E NO GUIDE Altitude Is A Function Of Mach Number Mode Calculate DHDV I-+ VDOT Equal Zerc Mode fi:~Ca~c$a~T Theta) Required = 0 VOOT Required VDOT Equal 1 Zero Mode (Calculate Throttling) F GUIDE Crdc:~~lote XM VF

0-r

Al lowablo Value Return

I

- Initial Data INITAL Purpose INITAL is called once at the beginning of a data case and performs the following tasks in the order listed below.

1) Read the title card; 2) Initialization of important program indicators and arrays; 3) Read the control card, NC array, and initial condition data; 4) Set the forward integration error limits; 5) Read the stage-dependent and constraint-dependent parameters ; 6) Print the initial conditions and stage-dependent parameters; 7) Construct heading block for printout from subroutine MATOUT; 8) Print optimization and stopping condition parameters; 9) Read and print free initial condition parameters if necessary; 10) Read and print phasing guidance parameters if necessary; 11) Print control parameters; 12) Read and print nonstage-dependent tables (including control variable and restart tables), Store part of control variable table on KTAN unit if necessary; 13) Store non&age-dependent tables (except control variable and restart tables) on KDAT unit; 14) Read and print stage-dependent tables and store on KDAT unit.

Subroutines Called TITLES Storage Used 2692 cells INITAL IO Print Data

-0

Variables And Arrays Card And NC LINET 2 NPP Read Initial Candi tion Data + Set Forward Integration Error Limits I State Variables Parameters

-_

I

INITAL Read Constraint Dependent Parameters (If There Are Constraints) Calculate GR OMEGA t Initialize VAR(28) 511 VAR(29) SI2 Cl1 Cl2 + Initialize CONST Array Calculate CONST 20) CONS? 21) I Print Note If b ntinn In Error\ Parametera And Constraint Dependent Parameters (If Any) IN ITAL onstruct coding Block Alloy For Printout From PI irlt Optimirotion Parameter And Stopping Parameter Print Error Note And Reset Partial Check Indicator Read And Print Phasing Guidance Porometers Print Control Parameters , 6, Read An rint Control Variable Table (Store In Core And nit If ,Necessary) C INITAL t Read And Print Restort Table If Core Is Being Restated I I I NO

J

NO

I

Rewind KPAR Uni And KSCR Unit L Return

T

Initialize State Variables INITCO - Purpose This subroutine is called at.the start of each stage and is assigned the following tasks: Initialize the X array at the start of stage 1; 1) Update the free elements in the X array at the start of stage 1; 2) Initialize weight at the start of every stage; Update weight at the start of each stage if weight is a free initial condition.

4) Method The X0 array is used to initialize X at the start of stage 1 and the WO array is used to initialize the weight at the start of each stage. The input WO and X0 arrays are updated in VALID using the values of DELX computed by VARIC.

Storage Used 117 cells IN ITCO Values Of X’S YES

Y

Reset Appropriate INITCO XKI = XKI + WO (NSTAGE) Zero The COF Array Calculate Calculate IXTRA 1, IXTRA 1, IXTRA 2 IXTRA 2 and IXTRA 3 and IXTRA 3 I Return

r-

KCALC -Step Size Logic Subroutine Purpose This subroutine is called following every try if MATRIX determines a non- negative majority vote. KCALC is assigned the following tasks: Compute maximum allowable constraint motion; 1) Set L(6) for next trajectory; 2) Compute step size coefficient (COEFK) for next trajectory.

3) METHOD KCALC computes the step size coefficient based on performance and con- straint linearities using the parabolic curve fitting method described in appendix C. Maximum step size coefficients are also determined based on the maximum allowable travel of each constraint. The best step size is then determined for the next trajectory.

Variables not in COMMON are transmitted between MATRIX and KCALC through the calling sequence CALL KCALC (DPSI, DPSIP, PSIBWD, PSIFWD, PSIK, PSIKTR, PSINL,

=)

where: ZZ change in constraint end points DPSI = predicted change in constraint end points DPSIP = nondimensional allowable change of a constraint in direction PSIBWD of desired constraint value = nondimensional allowable change of a constraint in direction PSIFWD away from desired constraint value = constraint step size coefficients PSIK = maximum step size coefficients of constraints due to maximum PSIKTR permissible travel = constraint nonlinearities PSINL = current end point values Storage Used 832 cells I KCALC Determine Step Size Bored On Performonce Linearity

I

t I Determine Step Size Bored On Constraint Linearities + Determine Permirrible Constraint ChClllp Store Step Sizes Based On Lineaitios In COKLIN Detormino Maximum Permirriblo Constraint Changes I Determine Maximum Stop Size Using Permissible Constraint CtKlnges A Determine Whether Or Not Each Constraint Is In Its Temporary Tolerance Bond I A KCALC Store Greatest Step Size Based On linearity In TRIALK Then Zero Out That Repeat KCALC If Control I ing With Constraint Within Tolerance Band Which Reduces Step Size Repeat Do Not Path Re eat K?AL‘ Check To Ensure Path Permissible Travel Of Constraint In Tolerance Band Is Not I> Error Coefficient Is Greatest TRIALK Obtained If Final Step Size Foils To Lie Between .5And 2. Use Appropriate Limit Return

i-I

LAMBDA - Reverse Integration Flow Controller Purpose LAMBDA controls the flow of the program during the backward integration.

Method The adjoint variables are initialized prior to performing the backward inte- gration. The automatic weighting matrix is calculated using the method described under “Automatic Weighting Matrices. ” Subroutines Called DVAL2 STP2 STEP2 Storage Used 518 cells LAMBDA On KDAT Unit Zero SENS In XB Array Zero XXLAM In XB Array Set Appropriato XB Elements To 1.0 Perform Initial Loading Of Portials lnto Allowable Core Storage ,.

’

LAMBbA Sst YES ll- VAR(B)-T -‘, c YES jaGEd) -T

I

I -TSTP (NSTAGE) Set IN= I L (5) = 1 B B LAMBDA

Q

Coil

STP2 Declemellt NSTAGE NO convet t X B EII,,I ln,licotol Decrement G9

‘a

Retw n

Fi

LOAD -Main Program Purpose LOAD controls the program flow from link to link and causes links to be overlayed.

Method each performing a specific function.

The program is broken into nine links, Link No. Purpose of Link General program flow Performs forward integration; called once per trial Numerical partial check link Closed-loop guidance link Performs the steepest-ascent calculations; called once per trial and twice per valid trajectory 5 Controls program flow for reverse integration; called once per valid trajectory Initilization, called once per data case 7 Performs final output of trajectory data; called once per data case Controls the plotting routines; called once per computer run Subroutines Called j ATMOS CARDS EXEC INITAL MATRIX LAMBDA PLOTZ Storage Used 206 cells

I - -

~ -- LOAD I NO Trajectory /

fe

A LOAD A Reset The Error Indicator To Zero Call PLOTZ

x

Return

I I

I

I

LOOKlD - One-Dimensional Table Look3 Purpose LOOKlD provides a table lookup technique for tables that have one inde- pendent variable and NDV dependent variables.

Method Given a value for the independent variable, the subroutine will furnish the dependent variable values and their partial derivatives with respect to the inde- Extrapolation is performed if the table limits are exceeded.

pendent variable.

The call is made to LOOKlD through a calling sequence.

Call LOOKlD (NDV, KK, XQ, Yl, Y2, Y3, Y4, Y5, Sl, S2, S3, S4, S5) where NDV = number of dependent variables KK = location of the first entry in the tables, NDS (table no.)

XQ = value of the independent variable Yl, Y2, Y3, Y4, Y5 = values of the dependent variables Sl, S2, S3, S4, S5 = values of the slopes of the dependent variables with respect to the independent variable Assumptions and Limitations The lookup method is based on linear interpolation.

1) The maximum number of dependent variables (NDV) is five.

2) Storage Used 229 cells LOOK I D f- Perform Lookup In Z Array Interpolate Or Extrapolate If Necessary Return

?l

LOOK3D -3D Table Lookup Purpose LOOK3D provides a method for linearly interpolating in a table to obtain a dependent variable as a function of up to three independent variables, i. e. , W = F(X,Y, Z).

Method Given the three values of the independent variables, the routine calculates the dependent variable and the slopes (derivatives) of the dependent variable with respect to the three independent variables. Linear extrapolation is performed in cases where the table limits are exceeded.

The call to LOOK3D is made through a calling sequence.

LOOK3D (NTAB, ZIND, XIND, AIND, YY, SLZ SLY, SLX) where NTAB = table number ZIND = value of the independent variable out of the plane XIND = value of the independent variable between curves on a given plane AIND = value of the independent variable that is the abscissa of the curves = value of the dependent variable YY SLZ, SLY, SLX = slopes of the dependent variable corresponding to ZIND, XIND, AND AIND, respectively Assumptions and Limitations The dependent variable is obtained from the input data by linear interpolation and the slopes by a forward difference technique. The lookup package can be used to obtain the dependent variable as a function of one or two independent variables also.

Storage Used 617 cells LOOK3D ind Curves On Either Side f Curve Desire *

F

t Find Points On Either Side Of NO

J

Find Points On Desired Point On A LOOK3D CompJto Slopes And Value Compute NO Slopes And - Value L Compute 510 es And P Vo U8 Find Curves On Either Side Of Curve Desired l

L-

Find Points On Either Side Of Desired Point Of Desired Point On Upper Curve

r-l

LOOK3D E’iher Side Of Desired Compute NO Slopes And < Value Compute Slopes And Value Return + -Matrix Output MATOUT Purpose MATOUT prints convergence information that appears following each trial The printout includes information such as current end-point values, or valid step.

constraint travel indicators, constraint changes, constraint tolerances, predicted constraint changes, step size coefficients, and other data.

Information not in COMMON is transmitted between MATRIX and MATOUT through the calling sequence CALL MATOUT (DBETA, DPSI, DPSIP, IMAJ, INDTVL, PSIBWD, PSIFWD, PSIK, PSIKTR, PSINL, XLAMDX, Xx) where: = the change in constraints modified by the end point changes due DBETA to the variation of initial conditions = actual change in constraints measured from a previous valid DPSI step DPSIP = predicted change in the constraint end points = value of the majority vote IMAJ INDTVL = constraint travel indicator = nondimensional allowable change of a constraint in the direction PSIBWD away from the desired constraint value ZZ nondimensional allowable change of a constraint in the direction PSIFWD of the desired constraint value = PSIK constraint step size coefficients = constraint step size coefficients based on maximum permissible PSIKTR constraint travel PSINL = constraint nonlinearities XLAMDX = predicted change in the constraints due to variable initial conditions.

xx = array of constraint end points storage Used 1202 cells MATOUT Print Heading Block For Control Logic Information NO rctnt I Ice Value .

ious Volid Print Current End ormance Point Values And Previous Valid J Step End Point Values Constraints i \ Per nrmance) And Print DPHIP DPHI . DPHID PHIN PHIK DPH13 Print DPH I And/Or DPSI I Print INDTVL CPSI INDSIC TTOL PSIFWD PSIBWD 1 PHINL Print PSIKTR And Or DPHIP d D SIP DBETA DPHID XLAMDX DPH13

---J

MATOUT Print Variable Initial Condition Name And tncremental Change To Be Mode To It I

I

IN0 Print Step-Size Coefficient

r-l

Print ~e,~tB~~ectory i Trial Set Indicator Print Next Trajectory Will Be A Valid Set Indicator I Return I

I

MATRIX - Basic Steepest Ascent Subroutine Purpose MATRIX is called at the end of every valid trajectory, reverse integration, and trial trajectory. After a valid trajectory, MATRIX calls VALID to examine the results (see VALID). After a reverse integration, MATRIX performs the following tasks : Recovers the I matrix from the XB array

1)

Calls MATRXB to invert 1~~ matrix 2) Computes minimum allowable DUSQ 3) Calls VARIC to compute new DELX array 4) Updates iteration counter (ITC) 5) Writes arrays required for trials on KPAR unit for temporary storage 6) Computes Lagrange multipliers

7)

Calls UCALC to update control history

8)

Calls MATOUT to print convergence information

9)

Following each trial, MATRIX performs the following tests: Reads arrays from KPAR unit written after the previous reverse integration

1)

Examines the try to see if it passes the majority vote test. If the test fails,

2)

a new trial is attempted Calls KCALC to determine the step size coefficient for the next rajectory

3)

Performs steps 7, 8, and 9 above.

4)

Method The steepest ascent method and the convergence logic used in MATRIX are explained in detail in the analytical development section and appendix C Subroutines Called KCALC MATRXB VALID MATOUT UCALC VARIC Storage Used 3454 cells IX rint Iteration

IP I

And Trial N

I NETA Ard I

Compute Change In Optimization Function Compute And Print .

XI Store Terminal Valuer Of + Constraints In XX Array Compute BVEC I Compute NO B t Compute Change In Constraint Terminal Values D

nl

A Return MATRIX Compute Ma jori ty Vote Print Too Many NO Functions With 4 1 Adverse Travel

YES

Call KCALC Adverse rrint

I

A * -erse PHI Great Acceptable Change In PHI

x

B v ,,i;: f i i MATRIX I Check Performance With No Constraints t, t, -Force Vahd -Force Vahd Step If All Step If All I I + + Constraints Counted Constraints Counted Improved But DUSQ Improved But DUSQ d d I I MATRIX Test For Bounce - If Bounce Force Valid Step Update Trial Number .

YES Compute and Scale DELX I XLAMDX

r-l

....

MATRIX Compute BINV Invert Matrix B Reset Appropriate I Calculotc And PI int Minimum I PHI PSI Transpose Store BVEC, SFISSD, B, Xl, DPSIP on KPAR Unit MATRIX YES NO I YES )

1 NO

Toleran+,Bands 1 YES YES Return

tl

Return n MATR.IX Ini tiolize Perturbation Size Set and Update Indicators Print Gradient of PHI too Negative NO Vf Decrement Iteration Counter Set NETA = 0 Compute DPSIP XIAMDX DBETA Compute DPHIP Due To DBETA G

d

G MATRIX

Q

OBETA Transpose DUSQ = Tl -= XNUM = 0 DUSQ too Small To Meet DBETA Scale DBETA To Match DUSQ -T- f H MATRIX compute * XMU(1) NO YES Set Indicator Compute DPHIP Value

I

YES Set DUSQ Bounce For Limits Optimization -

I

MAT.R IX Compute XMU’S Compute BBVEC , Compute DUSQ Check On KPAP Unit MATRIX YES Return G YES Compute DPHIP Due To DUSQ L Return

El

- Matrix Inversion Subroutine MATRX2 Purpose MATRX2 computes the solution of a set of simultaneous linear equations, the inverse of a matrix, or the value of a determinant.

Method The solution of a simultaneous linear equations follows from B = AX or solving for X X = A-l B and the determinant of A is Thus, the inverse is required for the solution, required for the inverse of A. The matrix operations are all performed in double-precision. Matrix inversion uses the Gaussian elimination method.

The call to MATRX2 is made through a calling sequence CALL MATRX2 (N, LN, LM, A, B, E, D, MM) = the order of the matrix A where N LN = number of rows in the dimension statement for matrices A and B, i.e., A is an LN x LM matrix.

LM = number of columns in matrix B = array designation for matrix A A array designation for matrix B B = E = an array used as temporary storage of a column matrix D = value of the determinant of A if D = 0 the determinant is not computed D = 1 the determinant is computed D = a scale factor, the determinant equals the value of the determinant times the scale factor MM = error indicator 1 if solution is successful 2 if overflow occurred 3 if matrix is singular The answers or the X matrix replace the original A matrix.

Assumption and Limitations The size of the problem that can be solved is limited only by machine storage space.

Storage Used 596 cells MATRX 2 Double Precision Motrix Inversion Initialize Overflow lndi cotor Reduce Matrix Bock Substitution To Find Solution Index

,&

l-r

Selectton of Element A of Maximum Absolute Value in Each ,Reduced Matrix A Interchange Rows Ploci ng Row With Maximum Return / Placard Equation Evaluation PLAC - Purpose Subroutine PLAC is called from FPROG for each selected earoute constraint and performs the following tasks.

I) Compares the constrained variable with the allowable limits from the input placard table Computes the placard equation of motion value and updates the appropriate 2) element of the F array Stores the placard limit value and the slopes of the limit value with respect 3) to its independent variables (obtained from the input placard table) in the S array.

Method The call to PLAC is made through the following calling sequence: CALL PLAC (II, ITAB, KTAB, Al, A2, A3, A4) = F array index corresponding to the selected placard.

where II II is computed by FPROG.

ITAB = first table number containing the placard values for the selected placard. ITAB is computed by FPROG.

KTAB = number of independent variables defining the input placard table. KTAB may be 1, 2, or 3, depending on the placard formulation.

Al = Value of the independent variable for a one dimensional table or the value of the “plane If for a three dimensional table.

A2 = Value of the “curve” for a 3-D table. A2 = 0. for a 1-D table.

A3 = Value of the abscissa (“point”) for a 3-D table. A3 = 0. for a 1-D table.

As an example, the calling sequence for the QLIMIT = F(t) placard is CALL PLAC (II, ITAB, 1, T, 0.) 0.3 AK(4)).

For the NLIMIT = F(H, M), the call is CALL PLAC (II, ITAB, 2, 0. , X(K2), AK(5), U13).

Note that if the 3-D option is used for a 2-D table, the value of the ‘plane” is set to zero.

PLAC will set F(II) = 0. if the limit value is not exceeded or to (LIMIT - A4)2 if the limit is exceeded.

PLAC also computes elements in the S array as follows: S (I, 2) = limit value exceeded (max. or min. ).

S(I,3) = slope of limit value with respect to Al (from the input placard table).

S(I,4) = slope of limit value with respect to A2.

S(I,5) = slope of limit value with respect to A3.

I is the internal index for the selected placard (I = 1, NPF) and is computed by PLAC.

Subroutines Called LOOKlD LOOK3D Storage Used 216 cells

x I

Return --.

!-,( LOOKBD) Set C’S YES Compute Derivative Of Placard I Return ii SET S’S -Plotting Subroutine PLOTZ Purpose PLOTZ controls the plotting operations for the Orthomat and SC 4020 plotters. Operations for scales, headings, titles, etc., are included.

Limitations The maximum number of plots is limited to seven

1)

Six curves maximum per plot 2) Multiple-data cases cannot be plotted 3) 4) X- and Y-axes annotation occurs every 2 cm Titles consist of the following: 5) LINE 1 - DEP. VAR. NAME VS. IND. VAR. NAME LINE2 - The first 72 characters taken from title card on the data deck.

Subroutines Called SKALZ Boeing numerical plotting system (NPS) subroutines called CAMERA NOS LIB TIT LEB ORTHPP AXLILI BORDER SLLILI NWPAGE SCPP PAUSE GDIJLI NOSLIL PSLILI FORM Storage Used 1436 cells PLOTZ Print The Name Of Plotter Used.

I

r

Draw Grid Draw Alignment On SC4020 Axes On Orthomat Bed NO. Of Curves NO. Of Points PLOTZ No. Of Points Per Curve

I

NO + Read Data From Read Data From KPLT For .Plot KPLT For .Plot Number 1 Number 1 4 4 . .

NO NO Read Data Read Dota From KPLT For From KPLT For The Last Plot The Last Plot Y Y ’ Read Data From ’ Read Data From KPLT For Plots KPLT For Plots Between First Between First And Last. And Last.

80 80 -0 -0 + Rewind KPLT Print The Plot No., No. Curves/ Plot And Plot Variable Names I PLOT2 Determine Max and.Min of Data for this Variable

r-l

YES lndep.

Variable 1 De-!jl: No ’ Determine Appropriate Appropriate Scale For Scale For Indep. Variable Dep. Variable

I

PLOTZ Initialize Specifications For SC4020 Annotation I w Call AXLILI Draw ORTHO Axes Initialize Specs. For Orthomat Annotation I Initialize General Specs.

For Annotation + Call NOSLIL Y Axis Annotation * Set Up Specs.

For Plotting Straight Lines Between Data Pts.

PLOTZ

Q

Initialize No. Of Data Pts For This Curve.

X ad Y Data Skips.

I--

Call SLLI LI

(3

YES PLOTZ NO PLOTZ Return Return -Numerical Partial Derivative Check PRTIAL Purpose PRTIAL calculates the partial derivatives of all of the equations of motion for state and placard variables with respect to each of the state variables and These are used for checking the analytical partial derivatives control variables.

used in the steepest-ascent method ( see appendix B ).

Method The numerical partial derivatives are evaluated using a linear approximation by perturbing the state or control variable about the present value. The pertur- bations are made both positive and negative from the mean value. For example, if z = f(y1, y2, . . . . . yn) , where yi is either a state or control variable, the numerical partial derivative is given by . . . . . yi+h, . . . . . yn) - f(y1, Y2, *---9 Yimh3 -‘*** Yn) f(Y1, y2, aZ 2h where h is the perturbation increment.

Subroutines Called AKSTP FPROG Storage Used 743 cells PRTIAL Set Stat0 Variable Increments To Bo UIod W.~;t~rration I Save Current Value Of State Variable t Increment Value Of State Variable , Of Derivative For Incremented PRTIAL Decrement Value Of State Vorioble Compute The Numerical Portiols Of Derivatives With The Perturbed State Vorioble And Reset Values Of Derivotivw + Reset The Stote Variable Volue PRTIAL

Q

.

Print Time, Analytical Portiols, And Numerical Partials For State Voriobles Save Appropriate Control Vorioble Of Derivotivcs Increment Control Variable Value

I--

PRTIAL NO t Colculote Ul2 Cl2 I- NO NO + + Set u14 = U(l4) I I I- +- Call AKSTP

P

I

I

I

The Following 4 Entries Are Made In Consecutive Made In Consecutive Order I At A Time Order I At A Time I PRTIAL Sow Derivotive For Incremented Control Vorioble Reset Control Vorioble Value t--@ r Entry 3 I Computo Numerical Partiolr Of Derivatives With The Perturbed Control Voriable Entry 4 4i2t&F+@ PRTIAL Print Analytical Partials And Numerical Partials For Control Variables

L-J

t Reset Integration Indicator Return SKALZ - Plot Scales Purpose SKALZ determines the scales for plotting.

Method Scales are based upon the number of divisions allowed and the variable range covered. The call to SKALZ is made through the calling sequence.

CALLSKALZ (DIV, DMAX, DMIN, SDIV, SMAX, SMIN, STEPP, DIST) = maximum number of subdivisions that may be used where DIV = maximum value of independent or dependent variable DMAX = minimum value of independent or dependent variable DMIN = number of divisions used for x or y axis ( < DIV ) SDIV SMAX = maximum value of independent or dependent variable (1 DMAX) SMIN = minimum value of independent or dependent variable ( 5 DMIN) step size between successive annotation on an axis STEPP = ( SMAX-SMIN/SDIV ) = distance to alter the origin of axis given in one- DIST division increments, ( DIV -SDIV )/2.

Storage Used 412 cells SKALZ Calculate Range Of Data- Calculate Range Per Division (DE LT) Find Power Of 10 (ISF) So DELvlO **ISF Lies Between 1 And 9 SKALZ Calculate NDIV, Number Of Divisions Required To Fit Data With Calculated Step Size DIV = DIV-1.

I I

T

Adjust Origin If Necessary Return -Integration Of Forward Trajectory STEP1 Purpose STEP1 integrates the equations of motion defining the state variable.

Method The integration package is a variable-step Runge-Kutta. The Runge-Kutta technique is basically a fourth-order method that is mechanized to vary the in- tegration interval based on two half-step intervals as compared with the full step.

This system evaluates the NV derivatives 11 times per integration interval.

only the basic NV equations of motion are integrated using the 11-point method.

The other equations of motion are integrated using only four evaluations of the A complete writeup of the RKVS method is derivatives to save computer time.

given in reference 8.

Assumptions and Limitations The minimum step size is specified by the user.

Subroutines Called AKSTP CONTR FPROG STPl Storage Used 1179 cells STEP2 - Adjoint Equation Integration Purpose STEP2 integrates the adjoint equations of motion backwards along the trajec- tory.

Method Integration is performed using a modification of the Runge-Kutta variable- step technique given in reference 8.

The call to STEP2 is made through the calling statement CALL STEP2 (UULAM) where UULAM is a matrix of the weighted impulse response.

Subroutines Called DVALB STP2 Storage Used 2323 cells STPl -Store and Print Forward Trajectorv Data Purpose STPl: 1) Calculates many of the AK variables, including metric conversions and transforms relative state variables to inertial coordinate system; Performs the circular satellite option calculations; 2) Stores partials and plotting data on tape; 3) 4) Tests for maximum stopping time; 5) Controls phases for nominal-guidance modes; 6) Controls printout of forward trajectory.

Subroutines Called ANPRTL Storage Used 867 cells NO I

b

Colculato Printout Variablor

u

I) YES

I

STPI NO Check For A Guidance Mode Switching Condition f Increment Phase Counter For A Chonge Of Modes NO YES

A

PI int PI int Variables Variables For Current For Current Time Time t NO STPI Increment Storage Counter Numerical YES NO \ Partiol ChecL A Nominal > Option Selected YES (%RTL ) ( fib,,,) ii:*-.

< Loop Option > ‘-- - --.. -- NO I NO L Return Plots On KPLT -Reverse Trajectory Store and Print STP2 Purpose STP2 controls the storage and printing of the influence coefficients during the reverse integration.

Method The option to print the adjoint variables is made by the user (NC (7) # 0).

The print interval is selected and printing is done only at the completion of an integration interval. Storage of the influence coefficients for use in determining the new control variable history is performed. The partial derivatives stored during the forward trajectory are called from storage and placed in core for use by DVALB.

The UULAM matrix which is computed in DVAL2 is transmitted to STP2 through the calling sequence CALL STP2 (UULAM) where UULAM = W-lG’ h Assumptions and Limitations The maximum number of points that may be stored is 1359O/Nu~l, where NUPl is the number of control variables plus 1.

Storage Used 478 cells STP2 YES t NO Print Adjoints

cl

Print Error Message i-i Store Return Influence

cl

Coeffieients On KLAM Unit .

STP2 YES W Return NO,/d;ttt$~~$ c Return

To Be Loaded Into Core l-

YES Load New Set Of Partials Into Al lowable Core From KPAR Unit Return TITLES - Output Heading Subroutine Purpose TITLES is called from INITIAL once per data case and performs the following tasks: Constructs the packed array of headings (DESXUP) from the DESXTJK array; 1) Converts the input NC(l21) through NC(l60) values to internal state variable 2) indicies ; Constructs the plot headings for the PLTNAM array; 3) Calculates the number of equations of motion (NEON), placards (NPF), 4) state variables (NX), control variables (NU), and auxiliary printout vari- ables (NK) selected by the user.

Method All information required by TITLES is available in the input NC array and the DESXTJK array compiled in subroutine BLOCK.

Storage Used 363 cells TITLES YES

t

I I

J

Convert NC , .

Array to Internal Equation Store Titles For Sequence Each Plot in , L PLTN AM Array I I I + I I Cons true t Headings f For Auxiliary Printout Variables Return

I

Thrust, Lift, and Drag Calculations TL-ET ~.

Purpose TLD calculates thrust, fuel flow, and aerodynamic coefficients as required by the options. It also calculates some VAR’s and partial derivatives as required to perform the above calculations.

Method TLD uses the NC’s input to determine the options to be used for calculating Body axes data is converted to thrust, fuel flow, and aerodynamic coefficients.

wind axes for use in the program. Thrust, fuel flow, and aerodynamic coefficients data are obtained by table lookups.

Assumptions and Limitations For thrust options 3 to 9, thrust and fuel flow data are tabulated for one engine only. If more than one engine is used, the thrust will be multiplied by the number of engines as input.

Subroutines Called LOOKlD LOOK3D Storage Used 1,419 cells TLD Zero VAR Array Transfer To - Thrust Option Selected By Input NO _ Calculate TSTAGE T = F (Alpha, •+ H, VI p ,i i TLD YES Calculate H VAR (32) VAR(31)

1 I-

Calculate VAR( I) a TLD IN0 4 :T ;) F (Alpha, I + Calculate VAR (1) Calculate YES - VAR (33) VAR (32) VAR (31) I :T i) F (Alpha, I

--I

Calculate VAR (1) TLD Calculate VAR (33) VAR (32) VAR (30 Transfer to . Fuel Option Selected by Input , Calculate VAR (3) Colculote VAR 36) VAR 35) I TLD F (ETA

I

I

Calculate YES VAR (38) VAR (37) J Calculate YES - VAR 36) I VAR 35) SFC = F (Alpha -+ V) Colculote YES VAR (38) VAR (37) TLD 2ooo Calculate I

I ’

TLD Calculate YES _ I Calculate VAR (3) e --SFFHj; F (Alpha M + b Calculate VAR (31 Calculate y i 1’ c TLD I I I Calculate VAR (3) I Call

A

TLD Scale Data Due to Number of Engines

r-l

I I I Vacuum

I I

I

I

c

t CD, CL = F W, M, A) Calculate YES L Call LOOK3D .

talculate Calculate TLD CAO, CNA, + CNA3= F (M) I t Calculate CN CL CD I t YES f%i” -i?6f PCLPM ,4 VAR (20) PCDPM TLD for CDO, PCD, CLI, CLA, CL0 Colculate EC0 Calculate VAR 20) VAR (17) YES 4 VAR 21) M&5’ I - PCLPM VAR (18) VAR 14) VAR (19) VAR 15) I I Calculate YES - Control Variable History Update UCALC Purpose UCALC is assigned the task of computing the control variable history for the next trajectory.

Method The control variable history used to generate the last valid trajectory and the impulse response history computed during the last backward integration is recovered from the KLAM unit. The Lagrange multipliers are used to compute a new control history by + [UULAM] [ks] = UOLD UNEW YS and UNEW is stored on the KPAR unit. Next the first part of uNEW that will fit in core is stored in TIMEU and the remainder is written on the KTAN unit.

STPl reads KTAN to accomplish the control history data overlay in the TIMEU array.

, computed in MATRIX, are transmitted to UCALC The Lagrange multipliers by CALL UCALC (XMU) where XMU = the array of Lagrange multipliers Storage Used 532 cells UCALC YES YES Print Error Nate And

+ NO

Set Error Indicator Read Time, Influence Cosffieientr, And Control Variables From KLAM Unit Return

I

Compute Change (fa~\.l)~ Con tro I

r-l

f t Store Time And New Value, Of Control Variables In TIMEU Array If New Control Table Too Large For Available Core Storage, Store Excess On KTAN Unit Return Valid Step Monitor VALID - Purpose VALID is called from MATRIX at the end of every valid step and is assigned the following tasks: Decides whether to accept or reject the valid step.

1) Update free initial condition state variables.

2) Tighten temporary tolerance bands.

3) Updates INDSIC array.

4) Updates CPSI array.

5) Method A detailed description of the logic used by VALID is given in Appendix C.

If VALID decides to accept the valid step, switches are triggered to allow a reverse integration. If the valid step is rejected, the KSCR and KPLT units The new valid trajectory are backspaced to the start of the previous iteration.

then will replace the rejected one.

Elements not in COMMON are transmitted between MATRIX and VALID through the calling sequence CALL VALID (DBETA, DPSI, DPSIP, INDTVL, XMU, XLAMDX, Xx) where = the change in constraints modified by the end point change DBETA due to the variation of initial conditions DPSI = actual change in constraints measured from a previous valid step = predicted change in the constraint end points DPSIP INDTVL = constraint travel indicator XMU = array of Lagrange multipliers +&] = predicted change in the constraints due to variable initial XLAMDX conditions = array of constraint end points Storage Used 735 cells VALID S0ve Values For Previous And Current Valid Iterations I Print Return < ’ Toa Many Trials - Abort Ah Initialize DBETA, DPSIP, INDTVL Initialize DPSI Establish lnitiol NO

I

VALID Print Too Many Functions With Adverm Travel Save Previous YES Reset Initial Condi ti ens NO VALID lighten Temporary Tolerances .

I Compute Number Of Iterations Since Constraint Has Been Inside Temporary Tolerance Band.

Compute Non-Dimensianol Constraint Change (Asked For) I Print

DUSQ Mod1 fi l d

YES Set DUSQ = DP2TMY I b NO r , Set Reverse lndi cator Set TOMAX L(13) = 0 I Return + -Variable Initial Conditions VARIC Purpose VARIC computes the direction and magnitude to perturb the free initial con- ditions to improve performance and also controls the maximum payload option.

This subroutine is called after every reverse integration and results in an update of the DELX array.

Method The method used in VARIC is discussed in detail in the analytical develop- ment section.

Data not in COMMON and computed in MATRIX are transmitted to VARIC by CALL VARIC (SFISSD) where SFISSD = cl 3 ’ I$ Storage Used 284 cells NO NO l NO c DELX = i/2/Payload \/ VARK Check Far Violation Of Free Initial Condition Limits Set To Approprf ate Limit When Violation Occurs

l-l

I Save Current Value Of DELX In DELXS

<~@=llP

a NO Return 2l Program Flexibility. -One of the main features of STOP is its flexibility.

The program has been modularized on both the macro and micro levels. On the macro level, for instance, the steepest-ascent block of the program is linked to the equations of motion by only one subroutine, ANPRTL. Thus, the entire equa- tions of motion block may be altered or even replaced without requiring any changes in the steepest-ascent block. Every attempt has been made to assign one task to each subroutine, and to make each subroutine as independent as pos- sible from the others. Finally, the flow has been organized so that there is one INITAL; another called once per subroutine that is called once per data case, stage, INITCO; one called once per forward trajectory, EXEC ; etc.

On the micro level, the equations of motion and partial derivatives are modularized. This feature allows the selection of any subset of the 40 available equations for a given problem. Also, equations may be added or deleted easily.

A few examples of the program flexibility are given below.

Example 1: Additional printout variable example: Suppose it is desired to print, Beta = -1 h’ h 1 , w ic is not currently available in the AK array. AK(39), not in current use, is available for growth.

In subroutine STPl add the FORTRAN statement Step 1.

AK(39) = SQRT (ABS (AK(5)**2-1 a)) Step 2. In subroutine BLOCK put the words (J/J/BE) and (TAJJJJ) in DESXUK (159) and DESXUK(lGO), respectively.

Step 3. Set NC(119) = 1 in the input data deck. This selects AK(39) to be printed in the forward trajectory output.

No other changes are required.

Example 2 : Additional placard example : Suppose it is desired to constrain the weight flow, which is not currently an available placard.

X(K37) is an avail- able state variable for placard additions.

Step 1. In subroutine FPROG, after the two FORTRAN statements: 37 CONTINUE ITAB = ITAB + 2 .

add a call to subroutine PLAC. If the W limit is a function of time, the call would be, CALL PLAC(I1, ITAB, 1, T, O., O., AK(7)) Step 2. In subroutine ANPARP, after the statement: 37 CONTINUE add the partials of F(K37) with respect to each state and control v,ariabJe. Nc$e that the call to PLAC generates F(K37) = for the portion of the trajectory where the (W - Wlimit) fuel flow exceeds the limit value(s). Some typical partials are: fi aF(K37)=2(+ -W limit) 3 H )H J& aF(K37) = 2 6 _ W limit) 3 V 1V

aw aa

2EE?k2& - --

- Wlimit) iJ 01 30 .

aF(K37) = 2 & - Wlimit) z arl or in FORTRAN TERM = 2 * (AK(7) - S(I,2)) PFX(K37,K2) = TERM * VAR(35) PFX(K37,K4) = TERM * VAR(36) PFU(K37,Il) = TERM * VAR(37) * VAR(41) PFU(K37,13) = TERM * VAR(38) Step 3. In subroutine BLOCK the heading for Wpenalty function, which consists of the words (JJW/D@) and (TJPFJJ), is placed in DESXUK( 75) and DESXUK( 76)) respectively.

Step 4. Set NC(157) = -1, 1, or 2 in the input data deck depending upon the type of limit desired, i. e., minimum, maximum, or corridor respectively.

Add a placard table to the data deck for the fuel flow limit.

Step 5.

It is a relatively straightforward Example 3 : Equation replacement example : procedure to replace the complete equation of motion block. This capability is demonstrated in the example below.

Consider a cat standing on a free-wheeling disk at point A as shown in figure 22. The cat wishes to get to point B on the disk in minimum time. The cat is constrained to the disk for the entire mission. The properties of the cat are Cat At

Figure 22. CAT ON A DISC PROBLEM

m = 0.5 slug F =lIb and the properties of the disk are r = loft r. =7ft = 70 slug ft2 I The pertinent equations of motion are: E =U 8 = V/r DEG = V2/r + F/m Cos(Y-0) fJ = -W/r + F/m Sin( Y- 8) + =o z = F/I r Sin(Y -0) DEG & =-V/r. DEG-ti i =l.

P’F = (r-r limit) Now generalizing to the STOP formulation, the state variables are: X(K.l) =r (fi) X(K2) = e Pa X(K3) = u (fps) X(K4) = V (fps) = 4 (de@ X W) X(K6) = 0 (deg/s ec) X(K7) = 6 (de@ X(K8) = t (SW X(K21) = PF (ft2/sec2) The control variable is: =Y (de@ UP) The initial conditions are: X,(Kl) = 10 X,(K2) =0 X, (K3) = 0 X. (K4) = 0 X, (K5) = 0 x, (K6) = 0 X, (K7) = 0 X, (K21) = 0 The problem statement is: R = X(K7) = 0 stopping quantity 4 = t (minimize) performance +I= r = 10 RPF=O +2= The partials required are: oE/au =l.

= -V/r2 DEG 0 tJ/ar 30 /3V = l/r DEG = -V2/r2 3 IIJ/ar = 2V/r C/a V = F/m Sin (Y- 8) RAD a he aThe = - F/m Sin (Y- 8) RAD at/Or = +W/r2 0+/a U = -V/r ;,+T/aV = -U/r a+/30 = - F/m Cos (Y-8) RAD ,ti/,Y = F/m Cos (Y-0) RAD a& /a 0 = 1.

ad;/ar = + F/I Sin (Y-8) DEG

ad/a0 = - F/I r Cos (Y-8)

ad/u = + F/I r Cos (Y-0)

ai/ar = - ai)/ar

ai/aV = - ai/aV

ag/oo = -1.

aGF/ar = 2 (r - rlimit)

The subroutines requiring changes are FPROG, ANPRTL, and ANPARP.

--....

In addition, a dummy subroutine AKSTP may be used and part of STPl may be deleted to save computer time.

As a matter of interest, the solution to the above problem is shown in figure 23. The problem solution required approximately 4 hours of engineering time and 10 minutes of computer time.

Program and data overlay. -Due to the overall size of STOP, it is not possible to fit all the subroutines and data into the IBM 7094 available core as is discussed in “Operating Information. ” It was necessary to minimize core requirements while maintaining the pro- gram’s capability, flexibility, and execution speed in that order of importance.

To this end, program overlay, data overlay, andequivalenceare used.

Program overlay: The STOP program can be broken into five major func- tional groups of subroutines: data initialization, forward integration, backward integration, convergence control, and plotting control. Each of these groups could operate separately if a main program could control their execution and inte rlinking data.

Using the IBSYS overlay feature, the program was divided into the above mentioned groups or links and three additional links. Two of these additional links were added to minimize the storage required for the forward integration subroutines, and the last link controls the final output for a data case. These nine links are described below.

Controls the program flow or execution from link to link and remains LINK 0 - in core all the time.

LINK 1 - Forward integration link Links 1, 2, & 3 LINK 2 - Numerical partial check link comprise the Nominal guidance link LINK 3 - forward integration Convergence control link LINK 4 - B(o)

Figure 23. PATH OF CAT ON DISC

--- : LINK. 5 - Backward integration link - Initialization link LINK 6 - Final output link LINK 7 Plotting link LINK 8 - Under the IBSYS overlay feature, unless otherwise specified, all links will In an effort to speed up the load- be loaded onto the I/O unit assigned as SYSUTB.

ing of links, the two most frequently used links (forward integration and conver- gence control) are loaded on separate I/O units by themselves, and these I/O units are rewound after each loading.

See “Operating Information” for a breakdown of links by subroutine. Figure 17 illustrates the program overlay and gives approximate core requirements for each link.

Data overlay: Since all of the input data for large cases will not fit in core, some of the data is written on scratch tapes and overlayed during execution.

Control table overlay: The input control variable table is normally stored in the TIMEU array, dimensioned to hold 1000 data points. If the input table exceeds 1000 points, that part of the table which won’t fit in TIMEU is written in binary format on the I/O unit referenced as KTAN. After each forward integra- tion step, the EXEC subroutine tests the current time against the highest time currently stored in the TIMEU array to see if more data should be read from KI’AN into TIMEU.

A similar procedure is used in subroutine UCALC, where the updated control table is calculated. In this subroutine, the I/O unit referenced as KPAR is used for temporary storage as the updated control table is built. Again, a maximum of 1000 .words of the control table are stored in TIMEU, and the remainder written on KTAN.

If the total number of data points in the updated control table exceeds 13,590, the data case is aborted and a comment printed.

Nonstage and stage-dependent data overlay: Except for the control and re- start tables, all other tables, both nonstage and stage dependent, are initially read into the Z array, dimensioned 7000.

The nonstage dependent data is read first and written in binary format on the I/O unit referenced as KDAT. The stage-dependent data is then read, one stage at a time and written on KDAT. If the sum of the non&age-dependent and stage- dependent data for any stage exeeds 7000 points, a message is printed and the data case aborted.

During the forward integration, the Z array contains the non&age-dependent and stage-dependent data for the current stage. At the beginning of each new stage, the data for that stage is read from KDAT into Z, replacing or overlaying the data for the preceding stage. After the last stage has been integrated, for each valid forward trajectory, the end points of that trajectory are written on KDAT.

These end points at the end of a data case are read from KDAT and printed as the trajectory summary.

At the start of each forward trajectory, trial or valid, KDAT is rewound.

The nonstage-dependent tables are read and the stage-dependent tables for the current stage are read back into the Z array. Reading the non&age-dependent tables into the Z array at the start of each iteration is necessary to allow other links to use the Z array for storage while they are executing.

-Equivalence. - The equivalence feature is used liberally throughout the pro- gram to make maximum use of available core in each link. Another purpose of the equivalence is to allow use of parts of the NC array (or its equivalent) for subscripting. The two major arrays to which other variables are equivalenced are NC and Z .

Z array equivalencing: Figure 24a illustrates the manner in which the Z array is used. The first column, represented by the shaded area, presents the usage of the Z array during forward integration. Following the weighting table, which is always in the Z array and never equivalenced, the nonstagedependent data is stored starting at Z (201). Then the stage-dependent data for any one stage is stored. Both the non&age-dependent data and the stage-dependent data are destroyed by equivalence during the backward integration and then reloaded into the Z array from KDAT in subroutine EXEC.

The Z array, from Z(201) through Z(2872), is used during the backward integration (link 5 ) for temporary storage of the state variable values, deriva- tives, integration errors, and other values needed by RKVS (Runge-Kutta var- iable step) integration package.

The Z array from Z (2900) through Z(7000) and the TJWMEU array, (l-1000), are used during the backward integration for storage of the partial derivatives read from KPAR.

NC array equivalencing: The NC array is equivalenced in subroutines throughout the program to make use of subscripted subscripts; i.e. use the NC as a subscript. Figure 24b shows the manner of NC equivalencing. From the figure it can be seen that the control variable indicators, NC (77) to NC ( 80)) are equivalenced with the control variable indices 11, 12, 13, and 14. The state var- iable and placard indicators, NC (121) through NC (160)) are equivalenced with state and placard variable indices Kl through K40.

The purpose of this section is to provide a quick refer- Trouble Shooting.- ence to aid in correcting typical problems encountered during the operation of __- WEIGHTING TABLE

r

WC9 FB ARRAY (534) XS ARRAY NC(1) ERI .

ARRAY .

NERI ARRAY NC ;7*, 1870) AMAl ARRAY 2204) OH151 ARRAY 2538) XHIST ARRAY NC (al) .

.

TPU ARRAY .

I I NC (1201

NC(l21) KI NC(I2P) K2 .

.

. .

.

.

. .

. .

l NC; 60) K40 Figure 24a. Figure 24b.

NC ARRAY EQUIVALENCING

ZARRAY EQUIVALENCING

STOP. The most common problems, as with any large, flexible computer pro- gram, are input errors. Another common problem is slow convergence and/or excessive run time. Also, STOP may abort during a run, either because of a problem it has detected or an error detected by the computing system.

It is always possible to restart STOP from the last valid iteration, whether termination was a successful exit, a STOP abort, a system-abort, or a system failure. Problems that execute more than one hour have a reasonable probability of experiencing a system failure, which makes restart capability a necessary requirement.

An input checklist and a trouble-shooting guide are provided below to assist the user in obtaining satisfactory results from STOP.

Input checklist: The following checklist will aid in data deck preparation.

It includes the most common input errors encountered by users not familiar with STOP. The checklist assumes that all keypunch errors have been removed and that all data are placed correctly in their fields. In particular, check all integer data to make sure it is right adjusted. The checklist is in the same order as the input instruction guide.

Card set 1 (title card) No problem Card set 2 (controls ) Check MSTAGE. Remember that if a tilt maneuver is used, it counts as a stage.

Card set 3 ( NC controls) 1. Check last card for blank field in columns 61-65. If last card is full, add a blank card.

2. Check NC (4 ) . It should equal the number of enroute and terminal constraints plus payoff and stopping variables.

3. Check NC (16) and NC (17). NC (16) must be 1 if the input control table was punched by STOP. NC (17) must be 1 if the restart table is included.

4. Only state variables selected for integra- Check NC (18) through NC (32).

tion may be used as the stopping condition. Three basic requirements must be examined to be sure that the stopping condition is attainable. First, XSTP must be reached before the stopping time (TSTP) of the final stage.

Second, the stopping condition must be compatible with the constraints selected (e.g. do not select altitude if the constraint on flight path angle is zero degrees ). Finally, select a stopping condition that will always be monotonic in the last stage of flight and has a large terminal time derivative.

Note that the stopping condition is only examined during the final stage.

Performance and constraint functions must be state variables selected for integration. Be sure that the sign of the performance variable index is plus for maximizing and minus for minimizing.

Constraint-dependent parameters and constraint tables must be input in the same order as they are in NC (20) to NC (32).

5. Check NC (121) -NC (160). Do not select state variables that are not needed, such as path range or the losses. Do not use latitude and heading for non- rotating-earth problems. Instead, formulate the problem in terms of equa- torial flight. Also, do not select placards unless they are -or are in danger Note that run time is almost directly proportional to of -being violated.

the number of selected state variables.

Card set 4 (initial condition data) Values on the order of 100 1. Check DUSQ if a restart table is not included.

are reasonable for airplane-type problems. For other vehicles, use an initial value of about one tenth of the controllable trajectory time.

Check FINNER and BINNER. Values of 18 are typical. The run time may 2.

be decreased by using 15’s, but convergence may be slower. For rocket problems, values of 18 seem to work well, while for SST problems 15’s are better. In general, difficult problems require settings that will result in at least 500 integration steps.

3. Check X0 ( )‘s to verify they are in the same order as the state variables selected. Also, if more than seven state variables are selected, there must be two or three cards of X0 ( )‘s 0 The X0’s are a packed array of initial state variable values. Note, for example, that X0 (7) does not necessarily equal X (K 7 ) at the starting time.

Card set 5 (stage dependent parameters) Check to be sure eight cards are present if MSTAGE is less than eight, 1.

otherwise 16 cards.

Check TSTP. The last stage requires a stopping time that will not be 2.

reached under usual conditions. It should be set to allow an abort if the stopping condition is not met. For VT0 systems the first stage contains the tilt maneuver; TSTP for this stage represents the termination of the tilt maneuver.

Check DTMIN. Too large a value will slow convergence. SST problems 3 ‘ behave badly if DTMIN is much larger than 0.1 second.

4. The first run should be a nominal only so that weight, Check KT and KM.

thrust, and other data may be checked. KM and KT may be used to adjust errors in total propellant consumption and total impulse.

Card set 6 (constraint dependent parameters ) 1. Check to be sure that one card is present for each enroute and terminal con- straint and that they are in the same order as the constraints selected by the NC’s.

Check the DPSI values. Very small values may result in slow convergence.

2.

Values of about 0.1 on angles and 100 on altitude are reasonable.

Card set 7 (free initial conditions ) 1. Check to be sure there are NC (5) cards.

If a tilt maneuver is selected, by, will be applied at end of tilt (start of

2.

stage 2). A large maximum value is required.

3. To restart with variable initial conditions: a. If weight is a variable initial condition, all positive-stage weights must be updated to be compatible with the weights from the valid step from which restart is desired. (card set 5, card subset 1). Zero and negative values are not updated.

If a variable initial condition is selected on gamma and a tilt maneuver b.

is selected, table 2 must be updated with the value of gamma at the beginning of stage 2 from the valid step from which restart is desired.

All other state variables chosen as variable initial conditions are up- C.

dated in card set 4, card 3. These values are also updated from the valid step from which restart is desired.

Card set 8 {nominal guidance data) 1. Check to be sure NC (13) cards are present.

2. Check the sign of NSW (see instructions).

Card set 9 (tables ) 1. Check table 0 if constructed by the user. This table will not be extrapolated and, therefore, a sufficiently large time must be included as the last point.

Multiple time points may be used at staging time. This table must always be present regardless of whether or not a nominal guidance is used.

2. Check the numbers of the placard tables (see notes at the end of the input instructions).

3. All tables after table 0 may be extrapolated (linearly) in any direction.

SYMPTOM DIAGNOSIS REMEDY Table card count in error a) Check card count in each table System abort during a) 1.

b) Deck or data setup error See appropriate sections of input instructions data input b) Make sure card set 4 contains at least four Error in count of state c) c) cards if more than seven state variables were variables, terminal con- straints, free initial con- selected (NC(121) - NC(160)). Make sure card ditions, or guidance modes set 5 contains 16 cards if MSTAGE is greater than seven.

Make sure card set 6 contains one card for each selected constraint (termi- nal and enroute) (NC(4) - 2 cards required).

Card set 7 must contain NC(5) cards and card set 8 must have NC(13) cards.

Check data deck d) STG or END card missing d) or misplaced NC input error Check card set 3. Make sure last card has e) e) a blank field in columns 61-65.

f) Restart error Restart table required if NC(17) = 1.

f) See output description.

Trajectory fails a) Stopping conditions not 2. a) Correct table 0 or use closed-loop options.

to meet stopping attainable on last stage See suggestions for NC(18) in checklist.

condition or b) TSTP of last stage too b) Increase TSTP of last stage. Trajectory prints note, small will always abort if this time is reached.

MAXIMUM STOPPING TIME EXCEEDED 3. I STOP abort with Input data tables exceed Reduce the amount of data or divide the note, MAXIMUM allowable storage for a stage-dependent tables into more stages ALLOWABLE stage (tables 1 through 30 TABLE STORAGE plus tables 31 through 34 EXCEEDED for a stage) SYMPTOM DIAGNOSIS REMEDY Backward integration too Lower BINNER by 3 4. STOP abort with 4 a) note, NUMBER tight OF POINTS Storage too frequent Increase ISTOR2 by 1 (Recommended b) b) STORED IN BACK- solution) WARD INTEGRA- TION EXCEEDS MAXIMUM ALLOWABLE One or more Constraints not attainable Change problem statement 5.

a) a) constraints not Solution nominal dependent b) Start from different nominal b) met after 20 Partial derivatives in Check partial derivatives (NC(12) option) Cl _ iterations error Input weighting matrix Increase weighting of applicable control

d) 4

variable problem Placard too nonlinear Nonlinearities are often caused by spikes

e) e)

in the placard table. Try redesigning placard so that the placard has a more severe violation in the spike region.

Placard corners or dis- Avoid corners or discontinuities in

f) f)

continuities placard tables.

Excessive run Unnecessary state vari- During the iterative procedure, integrate 6.

a) a)

time (Note: Run ables being integrated only the minimum number of equations time is very prob- (e. g. if the problem can be phrased in lem dependent. terms of flight around the equator, and Simple problems bank control is not required, latitude and may run several heading equations are not required). Run iterations per time is nearly proportional to the number minute, while of state variables.

very difficult Unnecessary constraints Constrain only the state variables required.

b) W

problems may selected, require 1 hour per iteration) SYMPTOM IUC ME DY DIAGNOSIS Variable initial conditions Decrease initial increment per iteration Excessive end- 7.

point oscillations allow too large a variation (card set 7) for the free initial conditions.

per iteration Slow convergence a) Weighting matrix problems a) Try decreasing the input weighting matrix 8.

elements on all control variables during the initial part of the trajectory. Also, decrease weighting cn control variables that tend to be nonlinear (such as sweep) b) Integration too loose Decrease DTMIN and increase BINNER W and FINNER. The integration should never reach minimum step size.

Error in partials Check partials (NC(12) option) c) c) d) Poor stopping condition Select a stopping variable that is mono- d) tonic and has a large time derivative at XSTP. Problems that stop on time usually behave better.

Constraint tolerances Increase DPSI’s (card set 6) e) e) too tight EPSLN too large EPSLN (card set 4) must be small enough f ) f) so that end-point motion within the EPSLN band is small compared with the end-point motion predicted by the steepest ascent logic.

STOP abort with All constraints have moved 9. See remedy for symptom 8.

note, CONVER- less than their respective GENCE FAILED final tolerance bands and the change in performance is less than or equal to lo6 times the value of perform- ance on the previous valid step.

DIAGNOSIS REMEDY SYMPTOM STOP termination All constraints are inside their Congratulations 10.

with note, THE respective final tolerance bands, LAST TRAJE C- the change in performance is less the TORY APPEARS than or equal to 10m6times TO BE AN value of performance on the EXTREMAL previous valid step, and the denominator of the first Lagrange multiplier is less than or equal to one ten thousandth of Ic + .

STOP abort with The denominator of the first See remedy for symptom 8.

11.

note, GRADIENT Lagrange multiplier plus 10 percent of 144 is negative.

OF PHI TOO NEGATIVE TO CONTINUE STOP abort with Check problem statement and table 1.

12. za;;p in inversion of I ti ti note, OVERFLOW .

OCCURRED IN MATRIX INVER- SION - ABORT Singular I# + matrix.

STOP abort with Check problem statement and table 1.

13.

note, SINGULAR MATRIX - ABORT REMEDY DIAGNOSIS SYMPTOM An excessive number of rejected See remedy for symptom 8.

STOP abort with 14.

note, NUMBER OF valid steps in combination with an excessive number of trials have TRIES PLUS RE- been run in an attempt to run an JECTED VAUD STEPS FOR CUR- accepted valid step.

RENT ITERATION EXCEEDS 10 - ABORT Plotting Information The program will plot any state, auxiliary printout, or control variable versus time or any other state, auxiliary printout, or control variable. Each Plotting is performed on either the curve represents a particular iteration.

Orthomat drafting machine or the Stromberg-Carlson (SC) 4020 microfilm record- er. Orthomat plots can be made on one of two pregridded forms: ll- by l7-inch mm vellum, or 8-l/2- by ll-inch mm vellum. SC 4020 plots may be made on 9- by g-inch vellum (hard copy) or 16-mm microfilm.

Limitations. - The following are plotting limitations: 7 plots maximum; 1) 6 curves per plot maximum 2) 3) Multiple-data cases cannot be plotted; X and Y axes annotation occurs every 2 cm. ; 4) Plot titling consists of the following 5) LINE 1 - centered, 0.2 inch below the X axis DEP. VAR. NAME VS INDEP. VAR. NAME IJNE 2 -centered, 0.2 inch below LINE 1 The first 72 characters taken from the title card on the data deck.

Input controls required for plotting. -The following values should be exam- ined in the input when plotting is desired: ORTHOMAT NC(l) = 0, output will be on ll- by 17-inch pregridded mm vellum = 1, output will be on S-1/2- by ll-inch pregridded mm vellum SC 4020 = 0, output will be 9- by g-inch vellum (hard copy) = 1, output will be 16 mm microfilm = 2, output will be 9- by g-inch vellum (hard copy) = 3, output will be both hard copy and microfilm ORTHOMAT NC (3) ---__ = 0, multiple curves will be defined with plot symbols (see multiple curve identification) = 1, multiple curves will be drawn with new colors (see multiple curve identification).

SC 4020 Not Used NC(ll) ORTHOMAT & SC 4020 = 0, No plots wanted = +N, SC 4020 plots wanted = -N, Orthomat plots wanted Note: The nominal, or first trajectory, the last trajectory, and every INI valid trajectories will be plotted, until NITC is reached or .

six trajectories plotted. Example: NITC = 10, NC(ll) = -2.

The Orthomat will plot the nominal and iterations 2, 4, 6, 8, and 10.

Each plot to be made is defined by a pair of integers that are input in the NC array, starting with NC (6 1).

Figure 11 gives the plot index for all variables.

NC(61) Dependent variable for plot 1 NC(62) Independent variable for plot 1 NC(63) Dependent variable for plot 2 NC(64) Independent variable for plot 2 .

NCi73) Dependent variable for plot 7 NC(74) Independent variable for plot 7 Note: A comment card must be included with the control cards stating that unit Cl will contain plot output.

Multiple -curve identification. -When plotting on the Orthomat, multiple curves can be identified by either plot symbols or colored ink.

(Note: Colors other than black or red do not reproduce very well).

SC 4020 multiple curves are always identified with plot symbols: Ink Color Curve No. Plot Symbol 1 None Black 2 0 Blue Green 4 0” Purple 5 h Red 6 El Black REFERENCES NASA Request for Proposal L-5347 -A Research Study for Supersonic 1.

Transport Trajectory Optimization. NASA, Langley Research Center.

Denham, W. F. : A Steepest-Ascent Method for Solving 2. Bryson, A. E . ; and Optimum Programming Problems. J. Appl. Mech., June 1962.

Problem, Optimization Formula- 3. Hague, D. S. : Three-Degree-of-Freedom tion. Wright Field FDL-TDR-64-1 Part I, Volume 3, October 1964.

Optimization Techniques. Academic Press, New York, 4. Leitmann, G. : N.Y., 1962.

5. Bryson, A. E. ; Denham, W. F. ; and Dreyfus, S. E. : Optimal Programming, Problem with Inequality Constraints I: Necessary Conditions for Extremal Solutions. AIAA Journal, November 1963.

Denham, W. F.; and Bryson, A. E. : Optimal Programming Problems 6.

with Inequality Constraints II: Solution by Steepest-Ascent. AIAA Journal, January 1964.

7. Phase II-A Comprehensive Report - Sonic Boom and Noise. Boeing Report D6-8686-7, November 1964.

8. Anderson, S. L.: Runge-Kutta Methods of Numerical Integration for Flight Simulation Programs. Boeing Report D2-139527-1, September 1966.

9. Finke, G. F.: Standard Atmospheric Properties (1962). Boeing Applied Mathematics Report AS 1772.

APPENDIX A CONTROL VARIABLE CHOICE FOR POINT MASS EQUATIONS OF MOTION The point mass equations of motion used for optimization of trajectories using angle of attack as a control variable are unstable for some types of vehicles with low thrust/weight, high lift/weight airbreathing engines. For rocket- type vehicles or high thrust/weight airbreathing vehicles, stability problems are generally not encountered. Problems arise when flying open loop with an .angle -of -attack history defined. Analysis of the linearized point mass equations of motion with angle of attack as the control variable is given herein showing The use of pitch angle as the that a lightly damped oscillatory motion results.

control variable is anlyzed; it is shown that the equations of motion are stable and that no oscillation occurs.

The linearized point mass equations of motion, assuming constant weight and flat earth, are: .

v- niV, - L,- Lo+WV Lh > cT;+VLQ

( >

- mv 5 = Q Y- - mV 0 0 v -*I ‘;-$y =o where y = perturbed flight path angle V = perturbed velocity l h = perturbed altitude v = v/v0 h = h/h, Ty= V, aT/8v Dv= VoaD/av T’h = h,aT/bh Dh = h,bD/dh Da = eD/aol Lg = h,aL/ah L,=V aL/aV Lor = &/ao!

All of the state variables listed are perturbation values and ( 10 are evaluated at the initial condition. At the initial time the conditions are assumed to be .

T Do

To= O-

- gY, m Toao+Lo w fo= mV --=

(A21

m v. 0 .

V x0 = $ y.

Taking the Laplace transform of equations Al and placing in matrix form gives: - Dv

Tv

A - Ti-i - Di s -

mV I

m v. vO Lv + w - Lo Li; -- S

W)

mV mV vO -- 0 S hO The characteristic equation of the system of equations is obtained by equating the determinant of the three by three matrix to zero giving: The characteristic equation evaluated for a typical SST in a cruise condition, .

assuming y = 0, v. = 0, ho = 50,000 feet, and V. = 3322 fps To = Do = 96,300 lb W = 346,000 lb Lo = 344,000 lb = 0.0 144 radians CY L = 2.39 x lo7 lb/rad a al-J Lo aP -= -- = - 18.98 lb/ft ah P a-h aIJ-2 Lo - = 207 lb - sec/ft av- v, BD Do 6P -= --=- 5.32 lb/ft ah P ah 6D 2Do -= - = 57.9 lb - sec/ft av VO

aT

-=-

3.18 lb/ft Bh E = 5 3 lb - sec/ft au - gives S3 + 4.89 x 1O-3 S2 + 1.96 x 1O-3 s + 4.79 x 1o-6 = 0 which has the roots s = - 2.44 x 10-3 s= - 1.22 x 10’3 f j 4.43 x 10B2 These indicate that all roots are stable, but the damping ratio of the oscillatory pair, 0.028, indicates that they represent near neutral stability. The period is 141.6 seconds. Figure Al shows the results of a digital computer run for similar conditions using the complete nonlinear point mass equation angle-of-attack constant and perturbing the initial value of the flight path angle. The motion is slightly divergent with a period of about 148 seconds. The nonlinearity of the point mass simulation could account for the slight differences as compared to the linear approximation.

The dominant term in the natural frequency for the oscillating mode is sL/6h, or more specifically, the 6p/6 h (atmospheric density gradient). This term is large for winged vehicles. The important terms in damping are: (1) BT/av, strictly an airbreather term which is destabilizing; and (2) 6D/av, a stabilizing but the desire is to have vehicles with low drag, which makes this term term - For accelerating vehicles, a high thrust/weight helps keep relatively small.

damping large, as can be seen by writing / 1.1 A low thrust/weight reduces the damping of the system.

- I II LENVELOPE FOR a = CONSTANT I 3000 4mo 5000 6GGo 7ooo 9000 ENVELOPE FOR a = CONSTANT -1 .o - * .

-1.5- 8 = CONSTANT ----- a = CONSTANT ‘I = THROTTLING = CONSTANT yo= 10 ENVELOPE FOR Q = CONSTANT I I I I I I I I I 1000 2000 3000 4000 5000 6000 7000 8000 9ooo TIME-SECONDS

Figure Al SUPERSONICCRUISEWITH 8 AND a CONTROL

Introducing 8 as the control variable, equations A3 become (upon letting

a,=& y)

T; - D; _ Th - DE D@

g + Dor

‘- mV To m V.

In vo

0 m v.

-L;+W-L To + L To + Lo LI; o!

-- : O s-t 8 (A5) mV mV mV m Vo 0 0 0 vO -- 0 S ho which has the characteristic equation r -I l- T; - D; To + Lo!

+ mV 0 m Vo LY+W-Lo g D Vo Li; ---+ -+- vo m vo ho m Vo m Vo ( L; + W - Lo v, I-i; +--- +

=o

(W

In vo

hom v,

)I

When evaluated for the same flight condition as for the cy controller, it is s3 + 0.45 s2 + 4.0 x 10-3 s + 4.79 x 10-6 = 0 (A7) where the roots are approximately s = 1.2 x 10-3 s= - 0.45 s=- 7.65 x 1O-3 Thus each root is real and negative. Figure Al shows the effect of the 13 type 8 held constant and the flight path angle perturbed). The motion control (i.e. , is overdamped and no oscillation occurs. The damping due to the flight path angle (To + La /m Vo) dominates the dynamics and eliminates the unrealistic oscillation from the cruise flight path.

The use of the 8 controller essentially introduces a fourth degree of freedom with the assumptions of neutral stability, no short-period dynamics, and a position feedback autopilot.

APPENDIX B PROGRAM EQUATIONS, VARIABLES, AND CON- STANTS DEFINED The purpose of this appendix is to define the names placed in COMMON which describe basic and frequently used equations, variables, and constants.

Equations of Motion The equations of motion defining the state variables, auxiliary state variables, and placards are defined as F( ). The argument of the functions identify the state variables as given in figure 1.

F(K1) = dW/dt = -WDOT F(K2) = dH/dt = V sin YR F(K3) = dYR/dt = [(T Sin Cv + L) COS Cp Go/W - G co6 YR 1” + O2 R cos fi (cos @ cos YR + sinp sin YR sin (I,)] /VR DEG cos YR/R + 2 w cos /k?cos eR + ‘R ) F(K4) = d VR/dt = (T cos Cy - D) Go/W - G sin YR + f&3 R cos @ (cos &? sin YR - sin /3 cos YR sin @R) F(K5) = d fi/dt = VR cos YR sin J/JR DEG/R F(K6) = d $y+/dt = - [(T sin Q + L) sin Cp G,/(W VR cos YR) sin ,8 cos tiR/(R cos 6) + 2 W (sin p cos y R + vR (I! 2 - sin YR cos B sin +R/cos YR) + R sin B cos ti cos $R/ COS YR) DEG lvR F(K7) = dX/dt = VR cos YR cos J/JR DEG/(R cos 8) F(K8) = dTD/dt = 1 F(K9) = d RNG/dt = VR cos YR Ro/(6076.103R) F(KlO)= d AHI/dt = Q VR .

F(Kll)= dA V/dt = TG,/W + ti2 R cos @ (cos @ sin YR - sin @ cos Y, sin J/JR) F(K12)= d GL/dt = G sin YR F(K13)= d DL/dt = DG,/W F(K14)= d T VL/dt = T(l - cos 01) Go/W 305 I F(K15) = dER/dt = dH/dt + VR/dt/Go F(K16) F(K17) Not defined; available for growth F(K18) WW F(K20)

F(K21) = d 8 */dt = 8 - elim(t)

c

)

F(K22) = dcp*/dt = cq - ‘Plim(t)j2 F(K23) = drl*/dt ‘((7.1 - 711imtH, M,)2 F(K24) = dh*/dt = (A - hlim(H, M))’ F(K25) = da*/dt = (a! - olim(H, M))’ F(K26) = Not defined, available for growth F(K27) = dHD*/dt = (dH/dt - dH/dtlim (H))’ F(K28) =

dQt*/dt = (Q - Qlim (t))2

F(K29) =

dQM*/dt = (Q - Qlim (NI))2

F(K30) = dQar*/dt = (Qa, - Q”lim(M))2 F(K31) = dTEMP */dt = (TEMP - TEMPlim(t))2 F(K32) = dN*/dt = (N - Nlim (H, M))2 F(K33) = dRPA*/dt = (RPA - RPAprn W) 2 F(K34) = dH*/dt = (H - ‘Hlim (M))

F(K35) = dAP */dt = AP

- Aplirn (h,@)2

( F(K36) = dM*/dt = M - M c lim(H$2 F(K37) Not defined; available for growth F(K38) = dY */dt =( Y - Y lim(H, M)j2 Not defined; available for growth F(K39) F(K40) Partial Derivatives The equations defining the partial derivatives of the equations of motion with respect to the state variables are given as PFX(A ,B) where the first argument, A, identifies the equation of motion and the second, B, identifies the derivative state variables. The partial derivative with respect to the control variables are given as PFU(A, B), where again the first argument identifies the equation of motion and the second identifies the derivative control variable.

PFX(K1, K2) = -avir/a.H PWK1, K3) = - SW/~ CY . a&a yR PFX(K1, K4) = - a’W/BVR PFU(K1, 11) = -aw/m . mhe PFU(K1, 13) = - aW/ZMj PFX(K2, K3) = VR cos y R. R&D PFX(K2, K4) = sinYR PFX(K3, Kl) = - (T sin CY+ L) COSQ Go DEG/(W2VR) PFX(K3, K2) = (BT/aH sin CL!+ aL/6H) cos Q Go/W {[ cos /3(cos p CO6 yR - ~G/BH cos yR + w2 DEG /VR - VR cos yR/R + sin @ sin yR sin 6,)

>

PFX(K3, K3) = [(T cos Q + (aT/d~! . sine! + BL/~cY) DEG). Bo!/ByR.

+ u2 R cos p (-cos p sin yR cos Q Go /W + G sin yR + sin @ cos yR sin GR) /V, - VR sin yR/R - (T sin 01 + L) cos cp Go/W - G cos yR PFX(K3, K4) = i [ + sin p sin yR sin OR)] /Vi + cc? R cos fi (cos /3 cos yR + (8T/a VR sin CY+ a L/a VR) cos cp Go/ (W VR) DEG + cos yR/R PFX(K3, K5) = k? R [cos B (-sin @ cos yR -t cos /3 sin yR sin +,) - sin /3 (cos 4 cos yR i- sin @ sin yR sin GR) /V, -2 w sin j9 cos GR y R CL? R cos fl sin @ sin PFX(K3, K6) = ‘OS $R& -2 W cos /3 sin #R = (T cosof + ( ~T/~cY sin (Y + a L/ a o!) DEG) cos cp . 0 cx/, cp.

PFU(K3, 11) Go/W vR) PFU(K3, 12) = - (T sin OL + L) sin 56 Go/(W VR) = aT/8 r) sin cx cos Q Go . DEG/(W VR) PFU(K3,13) = i3L/a A. cos Q Go DEG/(W VR) PFU(K3, 14) PFX(K4, Kl) = - (T cos a - D) Go/ W2 PFX(K4, K2) = (aT/bH cos Q - lbD/aH) Go/W - aG/b H . sin y R + O2 cos #3 (cos fi sin YR - sin fl cos YR sin Q,) 8(w/ay PFX(K4, K3) = (-T sin Q! RAD + 8T/60! cos cy - 6D/8q R. Go/W + -G cos yR + L? R cos /5 (cos B cos yR [ + sin j5 sin y R sin $,$I l RAD PFX(K4, K4) = (d!T/6VR cos Q - aD/aVR) Go/W PFX(K4, K5) = 2 R dos @ (-sin fl sin yR - cos j3 cos yR sin 0,) [ - sin fl (cos /3 sin Y - sin fl cos y R sin $R) RAD R PFX(K4, K6) = -u2 R cos 6 sin fl cos YR cos t/~, . RAD PFU(K4, 11) = (aTha! cos Cy - BD/8 o! - T sinCxRAD)B(w/a 8.Go/W = eT/t9Q cos CY Go/W PFU(K4, 13) PFU(K4, 14) = -aD/aA. Go/W PFX(K5, K2) = -VR cos YR sin QR . DEG/R2 PFX(K5, K3) = -VR sin YR sin tiR/R PFX(K5, K4) = cos YR sin GR/R PFX(K5, K6) = VR cos YR cos GR/R PFX(K6, Kl) = (T sin 01 + L) sin 56 Go DEG/(W2 VR cos YR) PFX(K6, K2) = - (8Th H sin cf + a L/BH) sin rp Go/ (W VR cos yR) C +v cos yR cos GR sin p/(R2cos 8, 2R -w sin 4 cos fl cos GR/ (VR cos yR)]. DEG .

PFX(K6, K3) = [ - (aT/6a . sin 01 + 8 L/act) DEG- T cos o!

- (T sin 01 + L) sin yR/cos yR sin tp Go/

8de YR

) p R COS y ) + (2 sin eR - w R sin B cos sin yR/VR) m (wVR yR + VR sin /3 sin y R W cos @ / cos ~0s $R/(R cm 8) PFX(K6, K4) = T sin a! + VR (-BT/aVR . sin cy - BL/BV~)] .

((1 o 2 $R /(Vi cos yR) sin q Go/W + R sin /3 cos /3 cos > - cos y R DEG cos GR sin6 / (R cos 8) > EFX(K6, K5) = - VR cos yR cos eR/(R cos2fj) - 2 W (COS j3 + sin yR sin j3 .

sin +!J~/COS YR) - a2 R (cos2B - sin2 8) cos eR/(vR cos yR) sin h sin p/(R cos p) + 2 w sin y, l P FX(K6, K6) = VR cos yR cos. /? cos #JR/cos Y’ + u2 R sin pcos Bsin GR/ CO6 YR) tvR = PFU(K6, 11) - (T cos a RAD + 6T/aa sin or + aL/m) aa/a 8 .

[ sin CPGo/W VR COB %) ] DEG PFU(K6, 12) = - (T sin Q + L) co8 Cp G,/(W VR cos %) =- bT/Bq sin a sin 4. R

PFU(K6, 13) Go DEG/(W VR COS YR)

=- PFU(K6, 14) 6L/Bh sin Cp Go DEG/(W VR cos YR) = PFX(K7, K2) - VR co6 YR CoS #R DEG/(R2 cos 8) =- PFX(K7, K3) vR sin y R co6 GR DEG/(R co6 p) PFX(K7, K4) = cos Y R cos #R DEG/(R cos 8) cos GR sin J3 /(R cos 2 PFX(K7, KS) = VR cos YR 8) =- sin tjR/(R cos 6) PFX(K7, K6) VR cos YR = - vR cos y, R,/(6076.103 R2) PFX(K9, K2) = - VR Sin yR R,/(6076.103 R) PFX(K9, K3) = cos YR Ro/(6076. 103 R) PFX(K9, K4) = VR t9QhH PFX(K1O, K2) = 3Q PFX(K1O, K4) = - T Go/W2 PFX(K1l, Kl) = 6T/6H Go/W + u2 cos fi (cos fi sin YR PFX(K1l, K2) - Sin p cos YR sin $-$ = 8T/6a* s&r, Go/W + CL’ 2R co6 B (CO6 B COS YR PFX(K11, K3) + sin fl sin yR sin $,) RAD = 8T/8VR Go/W PFX(Kl1, K4) PFX(K11, K5) = k-2 R cos @ (- sin 6 sin yR - cos f3 co8 yR sin h) C - sin B (cos fi sin & - sin ficos yR sin-+&] RAD = - a2 PFX(K11, K6) R cos fl sin/3 cos YR cos (IrR RAD = aT/aa* Ba/B&G,/W PFU(K1l, 11) = 8T/t9?? Go/W PFU(K1l, 13) PFX(K12, K2) = 8G/8 H sin YR PFX(K12, K3) = G cos YR RAD a PFX(K13, Kl) = - D Go/W2 = aD/aH Go/W PFX(K13, K2) PFX(K13, K3) = bD/Ba @a/byR Go/W PFX(K13, K4) = aD/6VR Go/W PFU(K13, 11) = sD/iYQI iYa!/B8 Go/W = 8D/8h Go/W PFU(K13, 14) PFX(K14, Kl) = - T (1 - cos a) Go/W2 PFX(K14, K2) = BT/6H (1 - cos a) Go/W PFX(K14, K3) = (BT/Bcr (1 - cos a) + T sin a. RAD) 8~ /s % Go/W PFX(K14, K4) = BT/BVR (1 - cos a) Go/W @T/act (1 - cos a) + T sin Cy RAD) a cy/a 8 Go/W PFU(K14, 11) = PFU(K14, 13) = aT/ a t) (1 - cos a) Go/W PFX(K15, Kl) = VRGoPFX(K4, Kl) PFX(K15, K2) = VRG,PFX(K4, K2) = PFX(K2, K3) + (VRGoPFX(K4, K3) PFX(K15, K3) PFX(K2, K4) + VRG,PFX(K4, K4) + dVR/dt/Go PFX(K15, K4) = = VRGoPFX(K4, K5) PFX(K15, K5) = VRG,PFX(K4, KG) PFX(K15, K6) PFU(K15, 11) = VRG,PFU(K4, 11) = VRGoPFU(K4, 12) PFU(K15, 12) PFU(Kl5, 13) = vRGoPFu(K4, 13) PFU(K15, 14) = VRG,PFU(K4, 14) PFU(K21, 11) = 2(e -8 limb PFU(K22, 12) = 2( ~ -16 lim) = - 2(77 /s H + sM/‘H .arlIim/‘M) PFX(K23, K2) - qlim)(ar)lim.

PFX(K23, K4) = - 2(71 - 7jlim.) (aqlim /aM .aM/aV) PFU(K23, 13) = 2(rl - Vlim) PFX(K24, K2) = - 2 (h - h & (aAlim /dH + 8h lim/BM 6M/BH) Iim./aM .6M/aV PFX(K24, K4) = - 2 (A - ‘Iim)” PFU(K24, 14)

= VA - Ali,!

PFX(K25, K2) = - 2 (o! - ~ Iim.) (acw Iim/aH + eM/aH. 6cYIim /‘M) = 2(a! - (Ylim) PFX(K25, K3) 8ddYR PFX(K25, K4) = - 2 (” - align) a8!Iim/dM . aM/6V

PFU(K25, 11) = 2(a - slim) 80he

PFX(K27, K2) = - 2 F(K2) - HDIim aHDIim/6H [ - HDlim PFX(K2, K3) PFX(K27, K3) = 2 [FW PFX(K27, K4) = 2 F(K2) - HDIim PFX(K2, K4) [ PFX(K28, K2)

= 2 (Q - Q lim) eQ/aH

PFX(K28, K4)

= 2 (Q - Qlim.) aQ/aVR

P FX(K29, K2)

= 2 (Q - Qlim) (aQ/aH - 8WaH * BQlim/@@

PFX(K29, K4) = 2 (Q - Q lim) (eQ/eVR - aMfaVR * 8QIim/6M)

P FX(K30, K2) ) (a 6 &/6H - 8M/8H . SQa!Iim/BM)

= 2 (Q(u - Q~ lirn

PFX(K30, K3) = 2 (&a! - Qclr lim) (Q . ‘~/‘YR

PFX(K30, K4) = 2 (Qcu - &a liml (a 6 Q/eVR - BM/BVR . aQa!Iim/BM) PFU(K30, 11) ) Q8d8@ = 2 t&a - QQ lim = 2 (TEMP - TEMP Iim + 0.2 M2) 8TEM/BH PFX(K31, K2) )[U t9M/BH] + 0.4M . TEM .

= 2 (TEMP - TEMP Iim) (0.4 M . TEM . 6M/8VR) PFX(K31, K4) PFX(K32, Kl) = - 2 (N - Nlim) N/W = 2 (N PFX(K32, K2) - N Iim)[(BL/BH cos 01 + 6D/6H . sin Q)/W - (BNlim/BH + tZ’M/BH * 8N lim/aW] PFX(K32, K3) = 2 (N - sin o! + D cos a) RAD Nlim)[t- L + aL/BcY cos 0 + @D/&Y sin (y.(act/6yR)/W PFX(K32, K4) = 2 (N - Nlim) (8L/av~ co6 8! + aD/aVR sin a)/W [ - aM/BVR l aNlim/aM] PFU(K32, 11) = 2 (N -NIim) (-L sin@ +Dcosg)RAD + 8L/60! cos ‘a + ~D/BQ! sin ct

@a/W/w

PFU(K32, 14) = 2 (N - Nlim) (0L/aA cos 8! + tYD/BA sin a)/W = - 2 (RPA) - RPA PFX(K33, Kl) lim) RPA/W PFX(K33, K2) = B(RPA - RPA lim) G2/(W2 RPA) [(L + T sin CV) a (8L/8H + 8T/aH l sin a) + (D - T cos a) (a~/aH - aT/aH . cos Q ) = 2 (RPA T RPA lim ) Go2/(W2 RPA)[(L + T sin a) (aL/t?a’ PFX(K33, K3) + aT/‘a sin ct + T co6 CY . RAD) + (D - T cos a) (aD/a Cl!

aT/dcu cos (Y + T sin CY .

RAD)] ~cY/~YR = Go2/(W2 RPA) [(L +-T sin 0) (aL/aVR PFX(K33, K4) 2 (RPA - RPqim) + sin a) + (D - T CO6 a) (aD/aVR - aT/ VR .

8 T/6VR * ,

cos Q) J

= 2 (RPA - RPA lim PFU(K33, 11) ) Go2/(W2 RPA)[ (L + T sin CY) (aL/aa!

+ aT/aol sin Q + T cos CY . RAD) + (D - T cos a) (aD/aa ~T/~cY cos o! + T sin CY * RAD)]do/ae = ) Go2/(W2 RPA) [(L + T sin Cu) 2 (RPA - RPA Iim PFU(K33, 13) gI’/ar) . sin CY - (D - T cos CY) aT/ar) . COS 8!

= ) Go2/(W2 RPA) [(L + T sin cw) . aL/aA PFU(K33, 14) 2 (RPA lim + T cos CY) aD/aa] 0’ - ZZ 2 (H - Hlim) (1 - aM/8H . ‘Hlim/6M) PFX(K34, K2) = 2 (H - H Iim) (BM/‘VR . ‘H lim/‘M) PFX(K34, K4) = F(K35) F(K35)H] /lo0 PFX(K35, K2) [ H + 100 - ZZ F(K35) (ecu/a yR)/O. 1 F(K35) PFX(K35, K3) I cu’O.1 - = F(K3;VR] /lo PFX(K35, K4) WWVR + 1o - C = F(K35)B] /O. 1 PFX(K35, K5) WWB + o. 1 - C = PFX(K35, K6) F(K35) + + o. 1 - F(K35$ /O. 1 E z PFX(K35, K7) F(K35)X] /O . 1

FtKWA +o, 1-

[

= PFU(K35, 11) (a ah 8)/O. 1 F(K35&] WWa + o 1 - [ .

= F(K35)* ] /O . 1 PFU(K35, 14) WWA + o, 1 - [ = PFX(K36, K2) 2 (M - M lim) (8M/a H - aM Iim/8H) = PFX(K36, K4) 2 (M - MIim) BM/aVR = PFX(K38, K2) -2 (YR - + aYRlim/ah aM/BH) YRIim) (a %Iim’aH = PFX(K38, K3) 2 (‘R - YRlim) = PFX(K38, K4) -2 (YR -7

Rlim) 3 Y~li~/ aM aM/b VR

f / Auxiliary Printout Variables ), are functions that are desired for The auxiliary printout variables, AK( printout in addition to the state and control variables. The argument identifies State variables and printout variables in the variable as given in figure 13.

metric units are also defined as AK’s.

= T = (T) Tref XKT - PAe/144) Cos 6c AW) = L = CLQSref AK@) = D = CDQSref A K(3) =Q = 0.5PVR2 A J-W) =M=V /V AW4 R A = GCR = cos -’ sin p0 sin B + cos B,, cos B cos (x-a] 60 DEG AK(G) [: = WDOT = Tref XKM/IS,,f + WDOTinert

AK(T)

= VI= (R 0 cos /3 t- VR cos yR cos QR )2 + sin2 y AK(g) [. R 2 l/2 + (VR cm YR cos QR ) = yI = sin -’ [ VR sin YR/VI ] DEG

AK(g)

AK(lO) = 4 = sin -’ [VR cos YR sin tiR/(VI sin Y$] DEG AK(ll) =“=e-YR AK(12) = TEMP = TEM (1 + 0.2M2) AK(13) =VA ’ AK(14) = P From ATMOS Routine AK(15) = p AK(16) = TEM M = 4.448221615 XKl AK(17) = W AK(18) = HM = 0.3048 x (K2) l/2 AK(19) = RPA = (L + T sin S)2 + (D - T cos a)2 Go/W [ AK(20) = N = (L cos Q + D sin Or)/W AK(21) =QCt=Qa AK(22) = G =c, /R2 AK(23) =R=H+Ro

l/2 ~ (M2 _ I)1/8,H3/4

AK(24) =AP =A PJ (pop) From tables 33 and 34 AK(25) = CL From tables 23,33, and 34 AK(26) = CD AK(27) = L/D = L/D AK(28) = SFC = 3600/ISref M = 0.3048 X (K4) AK(29) = VR = 4.448221615 AK( 1) AK(30) = T M = 4.448221615 AK(2) AK(31) = L M = DM = 4.448221615 AK(3) AK(32) = QM = 47. 8802589 AK(4) AJw3) = WDCT”= 0.45387006 AK(7) Aww = 0.3048 AK(8) AK(35) =VI M AK(36) =&a M = 0.8356681645 AK(21) = 47.8802589 AK(24)

AK(37) =A P M

M = 367.322574 AK(28) AK(38) = SFC Not defined; available for growth AK( 40) I VAR Array Many variables, functions, and terms are placed in a COMMON array called VAR( ) because of convenience or frequent usage. This array allows for the addition of variables to the program without changing the COMMON block.

VAR(1) = Tref i VAR(2) =A From Tables 31 and 32 e ) VAR( 3) =‘JSref VAR( 4) = 6Va/BH VAR(5) = 6P/8H From numerical differentiation of ATMOS Routine = BP/6H VAR(6) VAR(7) = 6 TEM/6 H VAR( 8) = Initial trajectory time (TSTART) From Table 30 VAR(9) = APJ VAR(lO) Not defined; available for growth VAR(ll) VAR(12) VAR(13) = CD f = (H - Href) 6 CHf/B H Where Href and BCDf/bH are from Table 23 VAR(14) = ac,hH VAR(15) = ND/t3VR VAR(16) = KDhat VAR(17) = ec,hh VAR(18) = 6CL/6H VAR(19) = 6CL/dVR VAR(20) = 8CL/aQ VAR(21) = aC,/sh VAR(22) = (cos fi cos yR + sin 6 sin yR sin I,!),) VAR(23) = (COS B sin yR - sin 6 cos yR sin I,!J~) 2 l/2 VAR(24) = VcsR = (ROW cos fi - VcsL cos $I)2 + (VcsL sin 4) C Not defined; available for growth VAR(25) Last-stage final weight VAR(26) = Wf VAR(27) = WI Last-stage initial weight

VAR(28) =A?’ tilt

VAR(29) = y0 VAR( 30) = MLSP (Maximum last-stage propellant, input) VAR( 31) = BTref/6 H VAR(32) = c?Tref/UR VAR(33) = aTref/aO1 /a r) (= 0 for current thrust options) VAR(34) = a Tref VAR(35) = aWDOT/a H VAR( 36) = a WDOT/d VR = a WDOT/dQ VAR( 37) = aWDOT/ar) VAR(38) VAR(39) =pG=8+(xsin$R+ycos$R)/R VAR(40) = AG = x + (x cos QR - y sin IjR)/(Ro cosfi) VAR(41) = PAPT = a&J 8 = 1 VAR(42) = PAPG = ea/ByR = -1 VAR(43) = PAPF = ea/egb = 0 Note: VAR(44) through VAR(50) not defined at present CONST Array A number of constants or functions of input constants are in the COMMON Again the array allows adding constants to ) .

array identified as CONST ( COMMON without modifying the common block.

= sin PO CONST(1) CONST(2) = cos PO = cos 6 CONST( 3) Cl = cos 6 CONST(4) c2 = cos 6 CONST ( 5 ) c3 = cos 6 CONST( 6) c4 CONST(7) =cos bc5 = cos 6 CONST( 8) c6 = cos 6 CONST( 9) c7 CONST(lO) = cos bc8 CONST(11) = cos 6cg CONST(12) = cos bclo CONST(13) = cos 6c11 CONST(14) = cos bc12 CONST(15) = cos bc13 CONST(16) =cos bc14 CONST(17) = 6076.103 CONST(18) = ~ = 1.9 = (Sea level pressure) CONST(19) = PO CONST( 20) = RD= (Desired final radius) [ I’ /RD]1’2 CONST(21) = VcsI = = (Sea level speed of sound) CONST(22) = VA0 CONST(23) = L = (Boom PF reference length) Note: CONST(24) through CONST(30) not defined at present Miscellaneous COMMON Words in COMMON or in calling sequences that are not defined by an array are listed below. These are primarily defined for convenience for the programmer.

= Sti/8 H = VAR(35) PWPH = i&8 VR = VAR(36) PWPV = 8ti/8a! = VAR(37) PWPA = ati/ q = VAEt(38) PWPT = 8T/8H = (q . XKT l BTre{8H - Ae8P/8H/144) COs hc PTPH = 8T/8 VR = ?j l XKT . 8Tref/8VR l COs 6c PTPV = 8~/8a! = 7 l XKT . 8Tre48” * COs 6-C PTPA = 8T/8 ?j = XKT l Tref l ~0s 6, PTPT PQSPH = 8Q/8H l Sref = PQPH l Sref PQSPV = 8Q/8VR l Sref = PQPV l Sref PLPH = 8 L/8H = 8CL/8 H l QS + Cz’ PQSPH l QS + CL’PQSPV PLPV = 8L/dVR = 8CL/8VR = 8L/aa = 8CL/8a l QS PLPA PLPS = SL/aA = 8CL/aAy QS PDPH = 8D/8H = 8CD/8H l QS + CD*PQSPH PDPV = aD/8VR = 8CD/8VR l QS + CtiPQSPH = 8D/8a = 6CD/862 l QS PDPA 8D/8h= aC,/81\* QS PDPS = PMPH = (8 M/8 H) vR = - VR/VA2 PMPV = (8M/8VR)H = l/VA = (8Q/8H) VR = 0.5 VR2 l 8P/8H PQPH

PQPV = (8&/8V), = P VR

= R. = Earth Radius = 20,902,992 ft RO = Degrees/Radian DEG RAD = l/DEG = cos 8 CIl SIl = sin 8 = cos cg CI2 CK6 = cos l/JR SI2 =sin(b CK3 = cos y R SK3 = siny R CK5 = cos j3 SK5 =&l/3 = sin Ic, SK6 R CAKll = cos Q =sinQ SAKll =w XKl ZZ III1 U12 =cp UI3 =r) =A UI4 =Q AKll GO = Go= Mass-To-Weight Conversion, 32.174 GR = 1 = Gravitional Constant, 14.081718 x 1015 OMEGA = o= Earth Rotation Rate, 0.72921152 x 10m4 PGPH = 8G/8H = -2 G/R Internal Program Indicators The following indicators (L array) are used for controlling internal operations in the program.

Number of integration points stored during integration of last valid tra- L(1) jectory (Set in subroutine LAMBDA). Used in the recovery of the control variable history.

Subroutine C&iTR lookup memory index during forward integration and L(2) subroutine DVAL2 lookup memory index during backward integration.

Maximum travel indicator L(3) = 0 No constraints traveled more than maximum allowable.

= 1 One or more of the constraints traveled more than maximum allowable.

Minimum DUSQ indicator

L(4)

= 0 DUSQ > minimum = 1 DUSQ = minimum Subroutine DVAL2 lookup control

L(5)

= 0 lookup if T (time ) not equal to previous time.

= 1 lookup Current iteration phase

L(6)

= -1 valid step, nominal only (set when NARBY = 0) = 0 trial = 1 valid step = 2 reverse integration =lO endof run Not used

L(7)

Forward trajectory last integration step indicator

L(8)

= 0 Intermediate integration step = 1 Last integration step Plot indicator

J-J(9)

= 0 Don’t plot current iteration = 1 Plot current iteration Gamma tilt variable initial condition indicator

U10)

= 0 Gamma tilt not selected as variable initial condition = 1 Gamma tilt selected as variable initial condition Current number of records of partial derivatives remaining on KPAR J411) (partial) unit during reverse integration.

Half-step indicator J-412) = 0 Not a half step = 1 Half step - do not print Iteration phase indicator

U13)

= 0 Next forward trajectory will be a valid step = 1 Next forward trajectory will be a trial Secondary use Valid step rejected indicator = 0 Valid step not rejected = 1 Valid step rejected Numerical partial check indicator (incremented by 1 every stored integration step and reset to 0 whenever numerical partials are computed) = NC (12 ) Compute numerical partials # NC (12 ) Don’t compute numerical partials - --

APPENDIX C PROGRAM CONTROL LOGIC The following tests made after different phases of an iteration are in no way connected with the steepest-ascent logic. They are present only to ensure the convergence to a solution in a reasonable length of time. All of the tests are strictly empirical and all have at some time or another played an important part in the convergence of particular cases. Experience is the only answer that can be given to the question as to the need of some of the tests. The control logic as presented in its current form is by no means able to solve all problems with which it is confronted. The control logic has, however, successfully caused the convergence of a great number of cases ranging from high-speed reentry vehicles to supersonic transports, thus showing its flexibility and capability.

Accept or reject valid step tests. -Majority vote test: The constraints out- side their respective tolerance bands on the previous valid step are examined with respect to their travel on this valid step. The constraints with adverse the constraints with favorable travel are assign- travel are assigned a value of -1; ed a value of +l. The algebraic sum of these assigned values is then taken. If the value of DUSQ is large enough to allow optimization, the performance function is treated in a similar way and added to the sum. The direction of change in per- formance is not meaningful when the program is not optimizing and therefore is not included in the majority vote. The change in performance is just the result of the constraint motion in this case. If the majority vote (algebraic sum) is negative, the valid step just calculated will be rejected and another valid step will later be attempted. If the majority vote is either zero or positive, the valid step will normally be accepted and another iteration will be run if desired. It should be noted that the majority vote is neither computed nor tested after the first valid step (nominal or first trajectory of a restart) of a run.

Excessive adverse performance test: An allowable adverse change in per- formance is computed as the larger of ten percent of the value of performance on the previous valid step, and five percent of the greatest value of performance attained on any previous trial or valid step (including the current valid step).

If the change in performance is in the wrong direction and it is greater than the allowable adverse change, a step-size coefficient based on adverse performance is computed and the valid step is rejected. The nominal or first trajectory of a restart cannot be rejected because of an adverse performance change.

It should be noted that, in the above two tests, the valid step will not be rejected if the decision to attempt the valid step was made because too many trials were made or DUSQ is at its minimum value. A valid step is run after a specified maximum number of trials or when DUSQ is set to the minimum allowable value to obtain a new set of influence coefficients, and is run under the assump- tion that many of the constraints and possibly performance will travel in an unfavor- able direction.

Abort test after a valid step. -Exit on too many valid step rejections: Under certain conditions the program control logic cannot make a successful valid step and will continue to reject consecutive valid steps. If the number of trials taken before the attempted valid step plus the number of rejected valid steps ever exceeds 10, the case being run is aborted with an appropriate message.

Abort tests following reverse integration. -Matrix inversion: If the inverse ._~~ of I## matrix cannot be accomplished because of overflow during matrix inver- sion or the matrix is singular, an appropriate message is printed and the run is aborted.

Check for optimum: If all of the constraints are within their respective input tolerance bands, the change in performance is less than or equal to lo-6 times the value of performance on the previous valid step, and the denominator of the first Lagrange multiplier is less than or equal to one ten thousandth of IO@, the case is terminated and a message stating that the last valid step is an optimum is printed.

Detect convergence failure: If all of the constraint changes are less than their respective input tolerances and the change in performance is less than or equal to 10m6 times the value of performance on the previous valid step, the case is terminated and a message stating convergence failure is printed. This indi- cates that the constraints are not within the prescribed band and that improve- ment in performance and constraints cannot be made satisfactorily.

Gradient of Q is too negative ‘test: If the denominator of the first Lagrange multiplier plus ten percent of IQ<p is less than zero, the case is terminated and a message stating that the gradient of Cpis too negative to continue is printed.

Tests following a trial causing another trial to be made. -Majority vote test: This test is similar to the majority vote test previously mentioned, except in two areas. The constraints being examined are the ones outside of their res- pective temporary tolerance bands on the current valid step and if all of the constraints being counted improved, an indicator is set. If the majority vote is negative, another trial will almost always be made (only other test examined after a negative majority vote is the bounce test, which is described later). A zero or positive majority vote will cause the step-size coefficients to be calculated.

Excessive adverse performance test: This test is identical to the adverse performance test previously described. If there is excessive adverse perform- ance on the current trial, another trial is almost always run. (Only other test examined after excessive adverse performance is detected is the bounce test. ) Tests after a trial that will cause a valid step to be run. -End point search: If the step-size coefficient equals two and the control logic is in an end point search mode (i.e., dU2< d - ’ I -I d - ), the current trial will be accepted and P&b P a valid step will be run. An appropriate message will also be printed.

Penalty function test: The penalty functions are examined (if any). If one or more have been violated after having no value for at least one iteration, a valid step is run - provided that the step-size coefficient is greater than or equal to one. An appropriate message is printed.

Performance test (no constraints): If at least one trial has already been run and the case has no constraints, the performance change is examined. If the performance change is favorable and the step-size coefficient is equal to one half, the current trial is accepted and a valid step is run. An appropriate message is printed.

All constraints favorable test: If all of the constraints counted in the majority vote test improved, the step-size coefficient was not reduced due to a constraint change that was larger than its maximum permissible travel, and the step size (DUSQ) is reduced over the previous step size on the last valid step. then the current trial is accepted and a valid step is run. An appropriate message is printed.

Bounce test: If the trend after a set of trials following a reverse integration is either a monotonically increasing or decreasing step size, and on the current trial, the trend is reversed, the last consistent trial is accepted and a valid step is run. An appropriate message is printed.

Too many trials: A valid step will be run if the following conditions are met: six or more trials have been run since the last accepted-valid step, the changes in constraints being asked for are not too great (d $4 d (1,max), and the control logic is not in an end point search (d U2 > dD’ I$#-’ dg ). An appropriate mes- sage is printed.

Too many trials: This test has no other conditions on it except that if eight trials have been run, the current trial is accepted and a valid step is run.

The above two tests (too many trials) are not satisfied under normal running conditions. The program normally runs one or two trials before making a valid step.

Factors affecting the selection of step size coefficient, - The step-size coefficient calculation will greatly affect the decision of the control logic in decid- ing to run another trial or valid step. This calculation is thedominant factor in deciding between a trial and a valid step.

A step-size coefficientbased on linearity is computed for each constraint and per- formance. The maximum step-size coefficients are then calculated for each constraint using its respective maximum permissible constraint change in the parabolic approxi- mation. The largest step-size coefficient based on constraint and performance linearities is then selected. The control logic is then said to be controlling with the function corresponding to the largest step-size coefficient. The largest step-size coefficient is compared with the maximum step-size coefficients of the constraints outside their respective temporary tolerance bands. If the smallest value of the maximum step-size coefficients being examined is less than the largest step-size coefficient, then it is chosen as the new step-size coefficient.

If at this point an indicator shows that the step-size coefficient calculation has not been repeated (explained later) and the step-size coefficient is less than 0.5, the control logic will return from this calculation with the step-size coefficient set This will cause the running of another trial. If the above conditions equal to 0.5.

are not met, tests are made to determine if the step-size calculation should be Conditions under which the calculations can be repeated are: the repeated.

logic is controlling with a constraint that is within its temporary tolerance band, and the maximum step-size coefficient corresponding to this constraint is less than the current step-size coefficient.

The maximum permissible travel of the controlling constraint is checked to ensure that it is not greater than its error. If the travel is less than the error, the logic returns from the calculation without recalculating the step-size coef- ficient. If the travel, however, is greater than the error, the step-size calcula- tion is repeated. Repeating the calculation means that the next largest step-size coefficient is chosen (based on linearities) and the logic is then said to be con- trolling with this function (another constraint or performance). The logic then follows the same path as described above.

It is possible to repeat the calculation many times before choosing a final step-size coefficient.

The final step-size coefficient is then compared with the limits 0.5 and 2. If the step-size coefficient lies outside the limits, it is set to the closest limit.

A final step-size coefficient of 0.5 or 2 will probably cause the running of another trial, whereas a step-size coefficient within the limits will cause a valid step to be attempted.

Minimum allowable value of DUSQ. -The value of DUSQ required for d F is calculated after a valid step as d$” x I## -’ x d 5. 0.25 percent of this value is the minimum allowable value of DUSQ. After a trial, reverse integration, or rejected valid step, the value of DUSQ is compared with the minimum allowable.

If the current value of DUSQ is less than the minimum allowable, DUSQ is reset to the minimum allowable, recalculations are made with the new value of DUSQ, and a valid step is run.

The valid step under this mode of operation cannot be rejected. It should be noted that no trials are run while the program is in the minimum DUSQ mode.

A need for a minimum allowable DUSQ resulted from a number of difficult cases, which reduced the value of DUSQ to such a point that any changes in con- straint or performance values were due strictly to computer round-off error and noise. Control logic predictions under these conditions were not reliable and therefore caused an even further reduction in the value of DUSQ. It is obvious that once this mode of operation is entered it is extremely difficult for the program to break out of it. The minimum DUSQ concept successfully bypasses this prob- lem, yet causes another. It is possible for the program to cause the divergence of a case when it is in the minimum DUSQ mode because the control logic is being forced to use a value of DUSQ greater than it would otherwise choose. Experience has shown that difficult cases will run in the minimum DUSQ mode for a number of iterations and break out of it and run in the normal program mode. Examining these cases in detail has shown that a very slight amount of divergence has occur- red during the iterations run in the minimum DUSQ mode. This occurred at a point when the program found it extremely difficult to shape the flight path properly.

Once this was overcome the program stepped out of this mode and successfully optimized.

Minimum DUSQ for optimization. -The value of the minimum DUSQ for opti- mization is computed as (2.1-FXTRAl)*Tl where FXTRAl is the nondimensional amount of the constraint error to be eliminated in this iteration (increased by five percent every iteration) and Tl is the amount of DUSQ required to meet the cur- rent A B .

After a trial, reverse integration, or rejected valid step, the current value of DUSQ is compared with the minimum DUSQ for optimization. If the current value is less than the minimum value, the value of DUSQ is reset to this minimum value and recalculations are made with the new value of DUSQ. The control logic continues with no other change. At the start of a case the minimum DUSQ required for optimization is more than twice as much as that needed to meet the amount of constraint error being eliminated. This essentially forces the program to opti- mize from the beginning of a run. Without this logic, too many cases used the entire DUSQ to meet the constraints and thus never optimized. Once the con- straints were met, the program would then turn to the problem of optimization, which resulted in the convergence to a local optimum. The current logic will cause the program to optimize first and then meet the constraints.

Scaling of dB. -If the denominator of the first Lagrange multiplier minus 10s8 x IQVis less than or equal to zero (1st test) and the numerator of the first Lagrange multiplier is less than zero (2nd test), DBETA (d-), Tl, DELX, B XLAMDX, and DPSIP are scaled by the factor (DUSQ/T1)l 2. If the first test above is satisfied, the program is said to be in the end point search mode. If the second test above is not satisfied, scaling does not occur and DUSQ is set equal to Tl. If the first test above is not satisfied and the numerator of the first La- grange multiplier is negative, scaling does occur and the program is said to be asking for the constraints too quickly. The program will also operate in the end point search mode when the user selects the boundary value mode.

Nondimensional (1, change.

-The nondimensional amount of constraint error to be eliminated for the new iteration is computed after every valid step. The amount of constraint error asked for is normally increased by five percent each iteration. If, however, the program is not in its optimization mode and a constraint is inside its tolerance band, its respective value of nondimensional (II change is not increased every iteration and is set to 10 percent. This causes the program to pay more attention to constraints outside their respective toler- ance bands when in an end point search mode, minimum DTJSQ mode, etc. The nondimensional ti change for a constraint is reset once it is outside its tolerance band. It should be noted that, only after a minimum of 20 iterations, will a full 100 percent of the constraint error be asked for.

Temporary tolerance bands. -Initial values of the temporary bands are equal to 10 times the final tolerance bands (input by the user). The constraints are examined after a valid step to check if they are inside their respective tem- porary tolerance bands. If all of the constraints are inside their respective temporary tolerance bands, all of the tolerances are reduced to one half their value. The only exception to this is that a temporary tolerance band cannot be less than its final input tolerance band.

Nondimensional permissible constraint change. - Two nondimensional per- missible constraint changes are computed after a valid step for each constraint.

One is for motion toward the desired constraint value and one is away from the desired constraint value. The values are a function of the nondimensional con- - straint error asked for (CPSI). As successive iterations are run, the non- dimensional permissible constraint changes are decreased. The two formulas are: PSIFWD (forward motion) = 5.0 - 3.0 x CPSI -, i PSIBWD (backward motion) = 1. 05 - CPSI Tolerance band indicators. -The number of consecutive iterations that a constraint has been outside its respective tolerance band is accumulated and stored in the computer word INDSIC. This value is reset to zero whenever a constraint is inside its tolerance band.

CR-793 - 21