\ . A Weak Hamiltonian Finite Element Method 

for Optimal Control Problems 

Dewey H. Hodges’ and Robert R. Bless2 

Georgia Institute of Technology, Atlanta, Georgia 

Abstract 

A temporal finite element method based on a mised form of the Hamiltonian n-eal; 

principle is developed for dynamics and optimal control problems. The mixed form of 

Hamilton’s weak principle contains both displacements and momenta as primary variables 

that are expanded in terms of nodal values and simple polynomial shape functions. Unlike other forms of Hamilton’s principle, however, time derivatives of the nionient a and 

displacements do not appear therein; instead, only the virtual momenta and virtual displaceinents are differentiated with respect to time. Based on the duality that is obser\-cd 

to esist between the mixed form of Hamilton’s weak principle and variational principles 

governing classical optiinal control problems, a temporal finite element forinulation of the 

latter can be developed in a rather straightforward manner. Several n-ell-l;nowl problelxs 

in dynamics and optimal control are illustrated. The esample dynamics problem inr-olws a 

time-marching problem. As optimal control esaniples, elementary trajectory optimizat ion 

problems are treated. 

Introduction 

This paper examines a finite element approach to addressing optimal control problenis. 

Hamilton’s principle has traditionally been used in anaZyticaZ mechanics as a method of 

obtaining the equations of motion for dynaiiiical systems. Bailey [ 11 follon-ed by sm-eid 

others (see [‘2 - 41 for esample) obtained direct solutions to dynamics problems I.-. ’WlS 2 

form of Hamilton’s principle known as the la.sv of varying action thus opening the door for 

its use in computational mechanics. , 

More recently, it has been shown that expression of Hamilton‘s law as a n-eak form 

(commdhly referred to as Hamilton’s weak principle or HWP) provides a powerful alternative to nunierical solution of ordinary differential equations in the time domain i3. Gi. 

The accuracy of the time-marching procedure derived in [5, 61 is competitix-e xitli st andard ordinary differential equation solvers. However, in order to derive an uiicondition,?!!~- 

stable algoritlini in [5]$reduced/selective integration had to be used. Further col~li~u~~- 

tional advantages may be obtcned in so-called mixed formulations of HUT in 11-hich the 

Professor, School of Aerospace Engineering. Associate Fellou-, -41-A-4. 

Graduate Research Assistant, School of Aerospace Engineering. Student. \lei:1b~r. 

AI-4 A. 

(NASA-CR-185336) A HEAK HAHILTONIAN FINITE 

ELEfiENT HETHOU FOR OP’IIllAL CONTROL PROBLEMS 

(Georgia Inst. of Tech.) 30 p CSCL 12a 

N89 -248 73 

Unclas 

63/64 0217251 
generalized coordinates and momenta appear as independent unknowns [ 71. Therein, an 

unconditionally stablk algorithm emerges for the linear oscillator with exact integration. 

HUT’ also has shown to be an ideal tool for obtaining periodic solutions for autonomous 

systems, as well as finding the corresponding transition matrix for perturbations about the 

periodic solution [SI. These are complex two-point boundary value problems; its utility for 

these problems and its superior performance in mised form strongly suggest that it could 

be used in optimal control problems. 

In this paper n7e show that HWP in mixed form provides a useful parallel to optimal 

coiitrol problems. Although finite elements have been appliecl to optimal control problems 

[9], the present formulation is believed to offer advantages over existing formulations. For 

example, it allows simple choices for shape functions in the finite element formulation [7] 

which simplifies element quadrature. 

TVe begin by summarizing HIVP in both displacement and mised forms. The mised 

form is then applied to initial value problems in dynamics. Then, making use of the 11-ell1;nown analogy between dynamics and optimal control, we develop our Hamiltonian weal; 

form for optimal control problems. Finally, we apply it to two relatively simple optimal 

control problems to illustrate its power. 

Hamilton’s Weak Principle for Dynamics 

In [l - 81, one can see the potential of obtaining a direct solution in the time domain 

17ery much analogous to obtaining the solution of a beam deflection problem n-it11 the 

beam axial coordinate brolien into several segments or finite elements. In the present 

case, however, it is the time interval which is broken into segments; thus, the term “finite 

elements in time” has been adopted by several investigators. 

Only recently has a mised formulation of HTW been investigated as a computational 

tool for finite elements in time [7]. In this section we will formulate the mixed form of 

HWP and illustrate its application to dynamics problems. 

G e ne r a1 D eve 1 o p me nt 

To this aim, let us consider an arbitrary holonomic inechanical s3;stem. Tlie configuration is completely defined by a set of generalized coordinates q. Further, let us denow 

with L(q, q, t) the Lagrangean of the system, Q the set of nonconservative generalized 

forces applied to the system, and p = aL/aq the set of generalized momenta. Then tlie 

following variational equation, 1;nown as HWP [5], describes the real motion of tlie spem 

between the two known tiines to and tf: 

This particular variational equation is said to be in displacement form because it only involves the variation of q. Although this formulation has been shown to be of practical use in 
dynamics [5, GI, an even more useful formulation may be derived if independent variations 

in both displacements and momenta are allowed, resulting in a mised formulation. 

To derive the mixed formulation, we begin by defining the Hamiltonian as 

Talcing the variation of Eq. (2) and substituting for 6L in Eq. (1) results in 

Integrating the first term of the integrand by parts yields 

This is called a mixed formulation because it contains independent variations of q and p. 

It is also a weak form in the sense that all boundary conditions are of tlie natural or weal; 

type. 

There are two main advantages of the mixed forinulation over tlie displacement formulation. The first advantage is that tlie mised formulation generally provides a more 

accurate solution for a given level of computational effort than does the displacement formulation. The second advantage is that a simpler choice of shape functions is allo~ed. 

Note in Eq. (4) that time derivatives of Sq and Sp are present. Hon-ever, no time derivatives of q and p exist. Therefore, it is possible to implement linear shape functions for 6q 

and Sp and constant shape functions for q and I:, within each eleiiieiit. 

Let us break tlie time interval from to to tf into AT equally spaced elements. The 

nodal values of these elements are ti for i = 1,. . . ,N+ 1 n-here to = tl and tf = tx+l. We 

define a nondimensional elemental time T as 

The linear shape functions for the virtual coordinates and momenta are 

For tlie generalized coordinates and momenta 

3 
and 

When these shape functions are substituted into Eq. (4), one can either generate an implicit time-marching procedure for nonlinear problems or apply standard finite element 

assembly procedures to solve periodic or two-point boundary value problems [SI. When 

this forinulation is applied to the linear oscillator, a time-marching algori thin emerges that 

is unconditionally stable. Higher-order (so-called p-version) elements could be del-eloped, 

and they would certainly be attractive for linear problems or for nonlinear problems n-ith 

nonlinearities of low order. For nonlinear problems in general. use of the crude shape 

functions allowable with the mixed inethod would seem to be more efficient than use of 

higher-order shape functions in a p-version. The reason for this is that, xith the esception 

of tlie term involving Q, all element quadrature can lse done by inspection regardless of 

the order of the nonlinearities. Detailed comparison of these methods is beyond the scope 

of tlie present paper but is being undertaken by the first author at this time. 

Example 1: Noiiliiiear Initial-Value Problems 

Applying the shape functions of Eqs. (6) - (8) to Ey. (4) for an initial value problem, 

we obtain a recursive set of nonlinear algebraic eqdations of the form 

where n is four times the number of degrees of freedom of tlie system. Eq. (9) can be 

solved by a Newton-hplison method yielding an implicit time-marching procedure. The 

key advantage of using finite elements and a weal; variational approach over nuiiierical 

integration is that the solution (for linear problems) is stable for dl time steps. In other 

words, no matter how large a time step is used, a finite approximation of the solut.ion \-,-ill 

be obtained. This unconditional stability is obtained without ad hoc procedures sucli as 

selective or reduced integration which are necessary in displacement formulations. 

We also point out that 2qi = Gi +@;+I and 2Fi =$i +$i+l. Thus, it is possible to cut 

tlie number of equations and unknowns in half. This can be very useful for a multi-degree 

of freedom problem. 

4 
i 

Consider a simple pendulum composed of a lumped mass rn and a weightless bar of 

length C (see Fig. 1). The single generalized coordinate q is the angular displacement of the 

bar from the vertical. Denoting the kinetic energy of the system with K and the potential 

energy with V, then we may define the following: 

1 2 -3 

v = 17229e(~ - cosq) 

I{ = -me q- 2 

L = I< - v 

There are no nonconservative forces Q applied to this system. 

defined in Eqs. (6) - (8) into Ecl. (4) we obtain 

Substituting t = t; + rAt, along with Eq. (lo), and substituting the shape functions 

Carrying out the integration, we obtain the four independent equations 

(12) 

There are six unknowns; however, for an initial-value problem, vie will specify 6, and 5, 

and solve for the remaining unlmowns. Thus, Eq. (12) is of the form of Eq. (9). 

For this simple pendulum example, we will nondiliiellsiollalize the variables as follon-s. 

If we define w2 = g/l, then a dimensionless time step At may he defined as At = --At. 

Also, instead of solving directly for p, we will solve for the dimensionless p/772C22. 

= GO” and = 0.0. Thc equations will be sol\-ed 

for A? = 0.4, 0.8, and 1.G. Graphs of the solutions are shown in Fig. 2 and Fig. 3 and 

We will start our pendulum at 

5 
1 

At = 0.4 gives acceptable results for both displacement and angular velocity. Also, note 

that even tlie large 1.6 time step yields a finite approximation of the exact solution. 

Weak Principle for Optimal Control 

A definite analogy exists between the mixed formulation of HWP in dynamics a.nd 

tlie first variation of tlie performance index in optimal control theory. Specifically, there 

is a similarity between tlie generalized coordinates and generalized momenta in dynainics 

and tlie states and costates in optimal control theory. 

Since the mixed foriiiulation has proven to be so valuable in dynamics, we will formulate a weal; variational formulation of the performance index. When deriving this 

formulation, we will lceep two things in mind. First, the resulting formulation must satisfy 

the Euler-Lagrange equations and boundary conditions that have already been established 

in optimal control theory [lo]. Second, whatever terms are necessary 1\7ill be added to 

the forinulation to transfer all strong boundary conditions to natural or weal; boundary 

conditions. 

Gene 1- a1 D eve lop iiie 11 t 

We start with a performance index talcen from Eq. (2.8.4) of [lo]. The first variation 

of the performance index will be talcen in a standard manner, except that states, costates. 

and controls will have arbitrary variations. Rather than setting the first variation equal 

to zero, however, it will be set equal to an expression which contains the terms that are 

necessary to transform all boundary conditions to the natural or “weak” type. The final 

weak form is then obtained by integration of this equation by parts in such a way that no 

time derivatives of states or costates appear. 

Consider a system defined by a set of n states J: and a set of 172 controls u. Furthermore. 

let the system be governed by a set of state equations of the form i = f(z, u, t), We ma?- 

denote elements of the performance index, J, with an integrand L(I, u, t) and a discrete 

function of the states and time at tlie final time d[x(:(tj),t,]. In addition, any terminal 

constraints placed on tlie states may be placed in the set of q functions $[~:(tj)~tf) and 

adjoined to the performance index by a set of q discrete Lagrange multipliers v. Finally, we 

may adjoin the state equations to tlie performance index with a set of Lagrange niultiplier 

functions X(t) which will be referred to as costates. This yields a performance index of the 

form 

(13; 

Taking the first variation of J and setting it equal to an expression chosen so that all 

boundary conditions are of tlie weak type, one obtains 
I :: = 6tf (A%) It, - sgio - SAT (3: - x) 

where 

The right hand side of Eq. (14) contains ternis necessary to form all of the proper boundary 

conditions as natural ones. The quantities 3: and X are discrete values of the states and 

costates at the initial (with subscript 0) and final times (with subscript f), Depending on 

the problem, these values will either be specified or left as unlmowns; they need not coincide 

with the values of the states or costates taken on within the elements at the beginning and 

at the end of the time interval. This is yet another indication of the “weakness” of tlie 

formulation. 

From Eq. (14), we niay directly write don711 a weak formulation. However, as stated 

earlier, one of the requirements of the formulation is that it satisfies the Euler-Lagrange 

equations and boundary conditions. To sl10~7 this, we will integrate the SiT term in 

Eq. (14) by parts, expand the variation of L, substitute Eq. (15), and group ternis yielding 

lot’ {SAT(? - f) - su [(E) + (Z) ”I 

-SzT (gp+i]}dt 

+ 6x7 (Af - if) - sxr (Ao - io) 

+ 6x7 (Gf - “j) - SX, T (20 - zo) = 0 

wliei-e xO,XO,xf, and Ai represent tlie values of those functions at the initial and final 

times, respectively. The coefficients of SAT, 6zT, and SuT in the integrand are the t.lree 

7 
correct Euler-Lagrange equations, Eqs. 2.8.15 - 2.8.17 from [lo]. There are also six trailing 

terms in Eq. (16) from which the boundary conditions can be determined. The first fsur 

of these equations and the sixth equation correspond to the correct boundary conditions 

in [lo]. Namely, the requirement for the coefficient of 6vT to vanish yields Eq. (2.8.21). 

The requirement for tlie coefficient of 6tf to vanish is equivalent to Eq. (2.8.20). The 

requirement for the coefficient of 65; to vanish sliows that the final value of x approaclies 

as given in Eq. (15), which corresponds to Eq. (2.8.19). If 10 is chosen as zero, the 

requirement for the coefficient of 62: to vanish requires the initial value of x to approacli 

zero as the accuracy of the approximation sclienie is increased (such as from adding more 

finite elements); on the other hand, tlie requirement for the coefficient of 6Xr to vanish 

requires the initial value of x to approach io, in accordance with Eq. (2.8.18). Finally, the 

requirement for the coefficient of SA: to vanish demands that the final value of z approach 

the discrete value if; this has no counterpart in [lo] since the elements of i?f are usually 

unknown. 

Having satisfied our requirement that none of the fundamental equations are altered, 

we may now derive our weal; formulation from Eq. (14). Recall from the mised formulation 

for dynamics that we do not want time derivatives of q or p to appear in the formulation. 

Similarly, we do not wish for the tinie derivatives of 5 and X to appear in the present 

weal; forniulation for optimal control. Therefore, we will integrate the k term by parts in 

Eq. (14). The resulting equation is 

This is the governing equation for tlie weal; Hamiltonian method for optimal cont,rol problems of the form specified. It will serve as tlie basis for the finite element discretization 

described below. It should be noted that normally one 11-ill encounter various types of inequality constraints on the states and controls in problems that deal with optimal control. 

This aspect has not been treated yet and will be tlie subject of future research. 

As in dynamics, we may choose linear shape functions for 63: and SA. We may choose 

piecewise constant shape functions for 2 and A. Thus, we will be working with shape 

functions similar to those of Eqs. (6) - (8). In addition, note that tlie time derivatives of 

u and Su do not appear in the formulation. This, we let 
u = u; 

su = sui 

Plugging in tlie shape functions described for z,X, and u, substituting T for t, and 

carrying out the simple integration from 0 to 1, we obtain 

This is the general algebraic forin of our Hamiltonian weak form for optimal control 

problems of the form specified. Eq. (19) is a system of algebraic equations. In fact. for 

AT elements, there are 2n(N + 1) + mhi + q + 1 equations and 2n(N + 2) + mN + q + 1 

unknowns. Therefore, 2n of the 4n endpoint values for tlie states and costates (20: Xo, ?j. 

and if) must be specified. In general, io (the initial conditions) is known in accordance 

with physical constraints. Also, if can be specified in terms of other unknowns with the 

use of Eq. (15). NOW we have the same nuniber of equations as unknowns. 

,. 

Normally Eq. (19) can be solved by writing ai explicit Jacobian and using a Ne\\-tonRaphson solution procedure. For the example problems which follow, the iteration procedure will converge quickly for a small number of elements with a trivial initial guess. 

Then, the answers obtained for a small nuinber of eleineiits can be used to generate initial 

guesses for a higher number of elements. Thus, a large number of elements can be solved 

with a very efficient run-time on the computer. 

Unfortunately, for some highly nonlinear problems, trivial initial guesses ma:- not be 

adequate. This problem may be overcome by noting that the 6zz and SS,;~ equations in 

Eq. (19) are n(N +1) equations that happen to be linear in the n(N+l) unlcnown costates. 

Now, we need n(N + 1) fewer initial guesses and we can solve for tlie costates in ternis 

of the other unlmoxviis. This is particularly useful since generating initial guesses for tlie 

costates can be difficult. This has proven to be a very useful way of obtaining ansn-ers for 

highly nonlinear equations. 
Example 2: A Fixed-Final-Time Problem 

As the first optimal control problem, we will examine the transfer of a particle to a 

rectilinear path (see Fig. 4). This is an example talien froin [lo], article 2.4. We will let 

z1 and 22 denote the position of the particle at a given time and 53 and 24 denote the 

particle's velocity at a given time. The thrust angle u is the control and the particle has 

niass nz and a constant acceleration a. 

The state equations are defined as 

0010 0 

0001 

a cos u 

0000 a sin u 

(20) 

Thc final time T is fixed and we would like to maximize the final horizontal component 

of velocity-. Thus, L = 0 and 4 = x3(T). There are also two terminal constraints on t.he 

states. These are tliat the particle arrive with a fixed final height (h) and that the final 

vertical component of vcloci ty be zero. We do not care what the filial horizontal component 

of position is. These constraints can be stated analytically as 

$ = { x22; h} (21) 

The initial conditions are s(0) = 20 = 10 0 0 OJT. Finally, we will eliniiimte the .. un1;nown Xf's by writing it in terms of other unknowns. In accordance with Eq. (15) 

(22) 

With L = 0, and tf fised, along with all zero initial conditions on the states, then t.he 

general formulation of Eq. (19) takes the following form. 

- SUI [ At ("','Xi]} - au 1 -6V?'4lt, 

(23) 

10 
1 

This produces a system of nonlinear algebraic equations whose size depends on the 

number of elenient s N. 

These equations axe solved by writing doilrn an explicit Jacobian and using a NewtonRa,plison algorithm. Trivial initial guesses (that are never changed regardless of input 

parameters) are used for N = 2. These results are then used to obtain tlie initial guesses 

for a.rbitrary AT by simple interpolation. In all results obtained to date for this problem, 

no a.dditiona1 steps are necessary to obtain results as accurate as desired. 

Representative numerical results for all four states versus dimensionless time t/T are 

presented in Figs. 5 - 8. For this example, we have taken h = 100.0, T = 20.0, and 

41t/aT2 = O.SS97. The results for 2, 4, and 8 elements are plotted against the esact 

solution available in [lo]. It can easily 5e seen that AT = 8 gives acceptable results for all 

the states. 

In Fig. 9, the control angle IL versus dimensionless time t/T is presented. Once again, 

the results are seen to be escellent for AT = 8. 

Three of the four costatcs are constants for all time and this method j-ields tn-o of 

these exactly. Tlie third costate is very close to tlie exact answer. The fourth costate 

corresponding to the vertical component of velocity A4 is sho\1711 in Fig. 10. Tlie results 

compare nicely with the exact results. 

-4 plot of the relative error of the performance index J = 53(T) and the endpoint 

multiplier u1 versus the numnber of eknieiits is shown in Fig. 11. It is seen to be nearly a 

straight line on a log-log scale. Tlie slope of tlie line is about -2 which indicates that the 

error varies inversely with the square of N. 

Notice in Fig. 11 that there is a bend in the endpoint multiplier curve. It is a common 

characteristic of niised formulations for an error curve not to be nionotonically decrcasing. 

Example 3: A Free-Final-Time Probleim 

The second optimal control problem is similar to Example 2 escept that now the 

final time is free and we would like to obtain a given horizontal coniponcnt of velocity 

(U) in the minimum time (see problem 9, article 2.7 of [lo]). Our algorithm from tlie 

preceding esaiiiple is readily modified to fit this problem by noting the folloxing changes. 

The performance index is now tlie final time T; so 

additional endpoint constraint on the states; namely 

have 

4 = 0 and L = 1. Also, there is an 

that 53 = U. With these changes: we 

0 

= { .) 

V3 

(24) 

and 

11 
(25) 

Along with these changes to our equations, we also pick up the additional Stf ecluation. 

The new system of equations is solved in tlie same manner described previously. Again, 

trivial initial guesses are satisfactory for N = 2, and these answer are used to obtain initial 

guesses for arbitrary N. 

Representative nuinerical results for all four states versus dimensionless time t/T are 

presented in Figs. 12-15 for a case with ah/U2 = 0.75. The results for 2, 4, aiid 8 elements 

are plotted against the exact solution available in [lo]. It can easily be seen that AT = S 

gives acceptable results for all the states. 

The control angle u versus dimensionless time t/T is presented in Fig. 16. Once again, 

tlie results are seen to be excellent €or AT = 8. 

-4s with the fixed time problem, tliree of tlie four costates are constants. The costate 

results are all as accurate or better then tlie costate depicted in Fig. 10. 

A plot of the relative error of the perforinance indes J = T versus the number of 

elements AT is shown in Fig. 17. As before, it is seen to be nearly a straight line on n loglog scale. The slope of the line is about -2 which indicates that tlie error varies inversely 

with the square of N. 

It should be noted that the computer time on a Cyber 990 is only about 2 seconds 

for AT = 2, N = 4, and N = 8 and 3 seconds for N = 16. Thus, tlie run tiiir,. is relatively 

insensitive to AT. 

Conclusions and Future Work 

In this report, a mixed form of Hamilton's Weak Principle has been >tated foi- dynamics problems. Finite elements in time were applied to this forinulatioii and a sii?iple 

initial-value problem has been used to demonstrate the principles involi-cc1 -4 ten:poral 

finite element method based on a mised form of the Hamiltonian weak priiit.iple \[-a< then 

developed for optimal control problems from the dynamics principles. It 1i:1-. been d( .: 11011- 

strated that the mixed form allows for a simple choice of shape functions an(\ is esseii:iall)- 

self-starting. Two simple optimal control problems have been examined and tlie resulrs are 

seen to compare excellently with tlie exact solutions for even a very few elements. Overall. 

the method provides very accurate results for the problems investigated to date n-itli only 

a few elements aiid for minimal coiiiputational effort. 

Future research will be done in applying this method to practical prolilerns such as 

development of on-board trajectory optimization algorithms for launch vehicles [ 111 ~ Such 

applications will require the relialility, efficiency, and self-starting characteristics illustrated in the present approach. Hon-ever, more work needs to lie done to solve the equations as efficiently as possible. One area with great potential for malting the method more 

12 
efficient is to take advantage of the sparsity of the Jacobian. The method will also have to 

be estended to allow for inequality constraints on the states and controls. 

Acknowledgements 

This work was supported in part by NASA Grant NAG-1-939 of which Dr. Daniel D. 

Moerder is the technical monitor. Technical discussions with Profs. A. J. Calise and 14. 

Borri are also gratefully acknowledged. 

References 

1. Bailey, C. D., “Application of Hamilton’s Law of Varying Action,” AIAA Journal, 

2. Sinilcins, T. E., ‘[Unconstrained Variational Statements for Initial and Boundary-Value 

Problems,” AIAA Journal, Vol. 16, 1978, pp. 559 - 563. 

3. Hitzl, D. L., and Levinson, D. .4., “Application of Hamilton’s Law of Varj-inS -Action 

to the Restricted Three-Body Problem,” Celestial Mechanics, 1701. 22, 19S0, pp. ‘253 

4. Baruch, M., and Riff, R., “Hamilton’s Principle, Hamilton’s Law, G” Correct Forniu5. Borri, M. et al, “Dynamic Response of h4eclianical Systems by a Weak Hamiltonian 

Formulation,” Computers and Structures, Vol. 20, No. 1 - 3, 19S5, pp. 495 - ;>OS. 

6. Peters, David A., and Izadpanali, Amir, “hp-Version Finite Elements for the SpaceTime Domain,” Computational Mechanics, Vol. 3, 19238, pp. 7.3 - 88. 

7. Borri, M., et ak., “Primal and Mixed Forms of Hamilton’s Principle for Constrained 

and Flesible Dynamical Systems: Numerical Studies,” ARO/-4FOSR TVorlisho1) 0x1 

Nonlinear Dynamics, Virginia Poly technic Institute and Stat e Uni versi t J-, B laclcsburg . 

Virginia, June 1 - 3, 198s. 

8. Borri, M., “Helicopter Rotor Dynamics by Finite Element Time -4pprosimat ion.-‘ 

Computers and Mathematics with Applications, Vol. 12-4, No. 1, 19SG, pp. 1-19 - 1GO. 

9. Pstten, William Neff, “Near Optimal Feedback Control for Nonlinear -4erod~-namic 

Systems with an Application to the High-Angle-of-Attack TVing Rock Problem.” -\I-\-\ 

Paper 8S-4052-CP, 19S8. 

10. Byson, Arthur E. Jr., and Ho, Yu-Chi, Applied Optimal Control, Blaisdell Publishin2 

11. Hardtla, J. W., Piehler, M. J., and Bradt, J. E., “Guidance Requirements for FLIture Launch Vehicles,” Proceedings of the AIAA Guidance, Navigation, and Co~troI 

Conference, Paper No. 87-2463, Aug. 19S7. 

Vol. 13, 1975, PI). 433 - 451. 

- 266. 

latioils,” AIAA Journal, Vol. 20, 1982, pp. GS7 - 692. 

I. 

Company, Waltham, Massachusetts, 1969, Chapter 2. 

13 
5 - - . . . - ... 

Fig. 1: Nomenclature for Example 1. This is a simple pendulum composed of a lump mass 

ni arid a weightless bar of length L 
60 

40 

20 

til 

a I 

v 

h 

w 

E: 

W 

F 

V 

c; 

", 

m 

zo 

L--.( 

.- 

R 

L 

3 

co 

E: 

2 -20 

4 

-40 

-60 I ' -1 

0 1 2 3 4 5 6 7 

D i mens ion less Ti me 

Fig. 2: Angular displacement q versus dimensionless time for t.hree values of the time step 

A?. These solutions are being compared to the exact elliptic integral solution. Note that 

even a large time step yields a stable solution. 

- .*- . .. 
l'Ol 

0.8 - 

0.6 - 

0.4 - c 

t: 

3 

u 

-- 

-0.4 

Exact 

A 1.6 

I 0.8 

0 0.4 

t' 

2 3 4 5 6 7 0 1 

Dimensionless Time 

Fig. 3: DiiiieIisionless momentum p/nd2w versus dimensionless time for thee values of the 

the step A?. These solutions are being compared to the exact elliptic integral solution. 

Note that even a large time step yields a stable solution. 
a 

f 

Fig. 4: Nomenclat.ure for Example 2. Transfer of a particle to a rectilinear path with a 

fmed final time. The problem is to maximize the final horizontal component of velocity 

- scbjzct. to specified terminal conditions on the states. 

-, ~ - 

I 

I -: 

i 

I 
1 .z 

1 .o 

0.5 

0.0 

A N=2 

R N=4 

N=8 

EXACT 

0.2 0.4 0.6 0.8 1 .o 0.0 

Dimensionless Time 

Fig. 5: Dimensionless liorizontal position zl/h versus t/T. Note that, the final horizontal 

component of posit.ion is not specified. 
1 .c 

0.E 

0.6 

0.4 

0.2 

0.0 

A 

81 

6 

0.0 0.2 0.4 0.6 0.8 1 .o 

Dimensionless Time 

Fig. 6: Dimensionless vertical position z2/h versus t/T. The final height is constrained to 

be h at the final time. 

..-. . ._ . 
2.5 

2.0 

1.5 

1 .o 

0.5 

0.0 

A 

A 

I I I 1 1 

A 

N=2 

N=4 

N=8 

0.2 0.4 0.6 0.8 1 .o 0.0 

Dimensionless Time 

Fig. 7: Dimensionless horizontal velocity 23 T/h versus t/T. Note that the periormance 

index J = zS(T). 

- 

I 

i 
1 .E 

1.5 

1.2 

0.9 

0.6 

0.3 

0.0 

A N=2 

N=4 

* N=S 

EXACT 

0.2 0.4 0.6 0.8 1 .o 0.0 

Dimensionless Time 

Fig. 8: Dimensionless vertical velocity z4T/h. versus t./T. The final vertical component of 

- ve1odt.y is constrained to be zero at the final time. 

.,- . . _.._ 
. a 

90 

60 

30 

(A 

a 

-60 

-90 

0.0 0.2 0.4 0.6 

N=8 

EXACT 

Dirnensionless Time 

0.8 1 .o 

Fig. 9: Control angle 'u versus t/T- The control is chosen so as to maximize the performance 

index J = z3(T). 
411 

-r 

-3 - 

-4 - 

2 3i 1 

1 I 

cs 

v 

z h 

P c -1 c.. 

-2 1 

a N=4 

N=S 

- Exact 

a 

-5 ! I 1 I 1 1 

0.0 0.2 0.4 0.6 0.8 1 .O 

D i me n s i o 11 I e s s Ti 111 e 

Fig. 10: Vertical velocity costate X;i versus t/T. The results for this c0stat.e are the least 

~ 2cccrat.e of all the costate.. 
\ 

-Performance Index 

---- End Point Multiplier 

\ 

\ 

\ 

\ \ \\ \ \ 

\ 

\ 

\ \ 

I 

100 

, \ , , , ,,'1 

Number of Elentents 

Fig. 11: Relat've error of the performance index 23 and the endpoint. multip 

A'. The error varies inversely with the square of N. 

1 

ier u1 versus 
0.0 

0.0 

A N=2 

N=4 

e N=8 P 

- exact 

A 

i 

A 

0.2 0.4 0.6 0.8 1 .o 

Dimensionless Time 

Fig. 12: Dimensionless horizontal position z1 /h. versus t/T. Note that. t.he final horizontal 

component. of position is not specified. 

. .,- --_. . _. 
1 .o 

0.23 

c 0.6 - 

\ 

c 

0 .- 

.c1 .- 

v) 

0 

CL 

c: 

V 

L 

I 

.- 

u 

0.4 $ 

0.2 

0.0 

A N=2 

N=4 

N=S 

- exact 

0.0 0.2 0.4 0.6 0.8 1.0 

Dimensionless Time 

Fig. 13: Dimeasionless vertical position x?/h versus i/T. The final height is constrained 

to be h at the final time. 

~ -. -. 
A N=2 

N=4 

N=8 

- exact 

A 

0.0 0.2 0.4 0.6 0.8 1 .o 

Ditnensionless Time 

Fig. 14: Dimensionless horizontal velocity x3 /lJ versus t/T. The final horizontal coniponent. of ve1ocit.y is constrained to be U. The problem is to reach this velocity in the 

Dinurnurn time T. 

.>-. .. --_ 
0.8 

0.G 

0.2 

0.0 

N=2 

- exact 

I I I I h 

0.0 0.2 0.4 0.6 0.8 1 .o 

Dimensionless Time 

Fig. 15: Dinlensionless vertical velocity 24 /rJ versus f./T. The final vertical component of 

velocity is constrained to be zero at the final time. 

I 

I 
90 

60 

30 

0 

-30 

-60 

-90 

0.0 

A N=2 

N=4 

€D N=8 

- exact 

I I I 1 i 

0.2 0.4 0.6 0.8 1.0 

Dimensionless Titne 

Fig. 16: Control angle u versus i/T. The accuracy of the initial control u(tO) deternines 

t.he accuracy of the costates. 
.1 

.o 1 

-00 1 

-000 1 

1 10 100 

Number of Elements 

Fig. 17: Reiat.ive error of the performance index T versus N. The error vzries inverselJwith the square of N.