- Email: [email protected]
A Finite Element Method for the Solution of Optimal Control Problems
A Finite E l e m e n t M e t h o d for t h e Solution of Optimal Control P r o b l e m s
Robert R. Bless Lockheed Engineering and Science Company Hampton, Virginia, 23666
Dewey H. Hodges School of Aerospace Engineering Georgia Institute of Technology Atlanta, Georgia 30332-0150
I. I N T R O D U C T I O N This chapter deals with a new type of numerical method based on finite elements in time for solving optimal control problems. The optimal control problem of interest in this chapter may be described as follows. Consider a system that is completely defined by a finite number of states, i.e., quantitles that describe the current status of the system. The status of the system is determined by a set of first-order ordinary differential equations referred to as state equations. The states are influenced by a finite number of control variables. The optimization problem is to choose the control variables to satisfy the given state equations, boundary conditions, and any constraints imposed on the states and/or controls, while minimizing or maximizing a given performance index, or cost functional. Use of the calculus of variations results in a multi-point boundary-value problem. Due to nonlinearities occurring in the state equations and the complex nature of constraints, not many analytical solutions to these types of optimal control problems have been found. Thus, the analyst must resort to numerical methods in order to solve the problem. CONTROL AND DYNAMIC SYSTEMS, VOL. 72 Copyright © 1995 by Academic Press, Inc. All rights of reproduction in any form reserved.
183
184
ROBERT R. BLESS AND DEWEY H. HODGES
A. METHODS OF SOLUTION Methods available for the solution of optimal control problems generally fall into two distinct categories: direct and indirect. Among the direct methods are those that transcribe the infinite-dimensional continuous problem to a finite-dimensional nonlinear programming (NLP) problem by some parameterization of the control variables and, possibly, the state variables. A discussion of NLP problems and the solution of these problems is given in [1]. Continuing advances in NLP algorithms and related software have made these the methods of choice in many applications [2, 3]. Indirect techniques, on the other hand, seek to minimize the performance index indirectly by satisfying the first-order necessary conditions for optimality as established from the calculus of variations (see, for example, [4, 5, 6]). Below is given a brief description of some of the indirect and direct methods now being used to solve optimal control problems. Since finite element methods may be derived either directly or indirectly, a separate discussion of them is given. 1. Direct Methods The direct approach to the solution of optimal control problems requires parameterization of the control and/or state time histories. The choice of parameterization schemes is not unique, and success of the direct methods has been achieved using schemes such as cubic polynomials [7] and Chebychev polynomials [8, 9]. Once the parameterization scheme is chosen, a parameter optimization algorithm is then used to improve the initial guess of the free parameters. These algorithms are in common use today and include quasi-Newton methods [10], sequential-quadratic programming (SQP) methods [1], and gradient methods [4]. Gradient methods were developed to surmount the "initial guess" difficulty associated with other methods such as NewtonRaphson. They are characterized by iterative algorithms for improving estimates of the control parameters in order to come closer to satisfying the stationarity conditions. First-order gradient methods usually show rapid improvements when sufficiently far from the optimal solution. However, the rate of convergence drastically decreases in the neighborhood of the solution. Second-order gradient methods have excellent convergence characteristics near the optimal solution, similar to a Newton-Raphson method. Conjugate gradient methods are very powerful because they combine the first-order and second-order gradient methods. A thorough description of the gradient method and many other algorithmic methods in optimal con-
A FINITE ELEMENT METHOD
185
trol may be found in [11]. Below are mentioned three of the widely used codes for the direct solution of optimal control problems. The first one is a FORTRAN code called NPSOL [12]. NPSOL is used to solve NLP problems and also serves as the solver for another code, OTIS, for solving optimal control problems. OTIS, or Optimal Trajectories by Implicit Simulation [13], is a three-degree-offreedom (point-mass) simulation program for multiple vehicles. The user can simulate a wide variety of vehicles such as aircraft, missiles, re-entry vehicles, and hypervelocity vehicles. A third code is POST, or Program to Optimize Simulated Trajectories [14]. POST provides the capability to target and optimize point-mass trajectories for a powered or unpowered vehicle operating near a rotating oblate planet. POST offers the solution to a wide range of flight problems including aircraft performance, orbital maneuvers, and injection into orbit. 2. Indirect Methods Indirect methods are based on finding the solution of a boundary-value problem which results from the first-order necessary conditions of optimal control. For many practical optimization problems the boundary-value problems are quite difficult; but, here again, modern numerical algorithms and associated software have enlarged the class of solvable problems significantly [15, 16]. Two methods for solving nonlinear multi-point boundary-value problems are shooting methods and quasilinearization methods. Shooting methods [17, 18] are frequently used and can be described as follows: The initial conditions and the differential equations are satisfied at each stage of the process while the final conditions are sacrificed somewhat. A nominal solution is generated by guessing the missing initial conditions and forward integrating the differential equations. The intent is to reduce the error in the final conditions at each iteration. Quasilinearization techniques [4, 19] involve choosing nominal functions for the states and costates that satisfy as many of the boundary conditions as possible. The control is then found by using the optimality conditions. The system equations and costate equations are then linearized about the nominal values, and a succession of nonhomogeneous, linear two-point boundary value problems are solved to modify the solution until the desired accuracy is obtained. Other indirect techniques include Newton methods [2.0], steepest descent methods [21], finite difference [22], and the method of adjoints and continuation methods
[lS]o
The virtue of the indirect methods is the high precision they offer and their rapid convergence in the immediate neighborhood of the optimal so-
186
ROBERTR. BLESSANDDEWEYH. HODGES
lution. These are important features, for example, when one is conducting sensitivity studies [23]. In addition, optimal trajectories are often composed of sequences of arcs along which various state or control constraints are alternately active and inactive, referred to as the switching structure [24]. Indirect methods require that the analyst find the correct switching structure, generally by trial-and-error, in order to obtain a converged solution. Even very short arcs must be anticipated. With NLP-based approaches, switching structures do not have to be known; however, the analyst can have difficulties in correctly identifying some optimal switching structures, particularly those with very short arcs that are missed by the parameterization scheme chosen with the direct methods. 3. Finite Element Methods Finite element methods originated in the field of structural mechanics in the 1950s, and the number and types of applications have since been greatly expanded. They are formulated by considering a set of spatial subdomains, called finite elements, which make up the domain of interest. The behavior of the set of finite elements, with certain inter-element continuity conditions enforced, is intended to approximate the behavior of the domain. The behavior of each element is analyzed based on a restricted set of modes, called shape functions. Usually the shape functions are low-order polynomials and are chosen so that the spatial degrees of freedom are either displacements or rotations at specified points, called nodes, or are a combination of nodal degrees of freedom and coefficients of polynomials that vanish at the nodes. We consider finite element methods separately from the direct and indirect methods, because finite element methods can be developed either from energy principles (directly) or from application of Galerkin's method to the governing partial differential equations (indirectly). In some instances the resulting spatially discretized equations are exactly the same, regardless of which path is taken in the derivation. Finite element methods have also been developed for the time domain. The earliest works seem to be those of Fried [25] and Argyris and Scharpf [26], although some of the later work, including the work presented in this chapter, has its roots in the mid-1970s with the development of direct time-domain methods for dynamics based on Hamilton's principle. Hamilton's principle had traditionally been used in analytical mechanics only as a method of obtaining the equations of motion for dynamical systems. Bailey [27] followed by several others [28, 29, 30] obtained direct solutions to dynamics problems using a form of Hamilton's principle known as the law of varying action, thus opening the door for its use in computational mechanics. Recasting of Hamilton's law as a weak form, referred to as Hamilton's
A FINITE ELEMENT METHOD
187
weak principle or HWP [31], was shown to provide computational characteristics that are superior to those of the form used by Bailey. Bailey's form was later shown to occasionally exhibit convergence problems [32]. The accuracy of the time-marching procedures derived in [31, 32] is competitive with standard ordinary differential equation solvers and provides a powerful alternative to numerical solution of ordinary differential equations in the time domain. In what appears to be one of the first papers using finite elements to solve optimal control problems [33] the authors applied a modified RitzTrefftz direct method to state-regulator control problems. Finite element methods based on Ritz and Galerkin methods [34, 35] as well as the method of collocation (see, for example, [36]) have also been used to solve optimal control problems. Other finite element methods based on the method of Ritz can be found in [37, 38], and one based on collocation is found in [7]. A direct method for solution of inverse response problems is presented in [39]. The method is set up to "march" in time using a direct representation of the dynamics via Hamilton's law of varying action instead of using the standard state-space ODE representation of the system dynamics. B. PRESENT APPROACH The methods described above, although very accurate and useful, suffer from some computational challenges. For example, the approximating functions must satisfy all strong boundary conditions, which means that one set of functions may not suffice for all types of problems or that certain equations and unknowns need to be removed depending on the particular problem. This turns the development of a general-purpose algorithm into an even more complex task. Also element quadrature in the previous methods must be done by numerical means, which greatly increases the computational effort. This chapter describes in detail a computationally efficient, generalpurpose time-domain finite element approach, taken from [40, 41, 42], for solving optimal control problems based on a weak formulation of the firstorder necessary conditions. The above problems are avoided by reformulating the variational problem so that all boundary conditions appear as natural boundary conditions. This way, the simplest possible approximating functions are allowed for all cases so that all element quadrature may be done by inspection, and the Jacobian matrix involved in the solution procedure is very sparse. The equations can be set up for solution prior to specifying a particular problem, making the algorithm ideal for generalpurpose computer software. The background of the present method of treating the time domain
188
ROBERT R. BLESS AND DEWEY H. HODGES
lies in extensions of HWP [31]. HWP was 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 [43]. These are complex two-point boundary value problems, and the utility of HWP for these problems strongly suggests that it could be used in optimal control problems. Finally, an important step was the development of a mixed form of HWP [44], in which the generalized coordinates and momenta appear as independent unknowns. The mixed form provides additional computational advantages over the original displacement form [31], in which only generalized coordinates appear. These advantages include the possibility of using simpler shape functions and its ability to yield an unconditionally stable algorithm for the linear oscillator with exact element quadrature. Although the early developments with Hamilton's law and HWP are direct methods, the mixed form of HWP is really just a weak form of the governing equations, including boundary conditions. Therefore, a finite element method based on it is an indirect method. The finite element method to be described herein is also an indirect method. Of the indirect methods listed above, shooting methods yield numerically exact solutions that are as accurate as the integrator used for the differential equations. Shooting methods are, unfortunately, very sensitive to initial guesses; hence, it may be difficult to obtain a converged solution without good initial guesses for the costates. In addition, shooting methods are not well-suited for real-time applications due to the computationally intensive operations associated with repetitive numerical integrations needed to calculate the J acobian. Therefore, the purpose of this chapter is to derive and evaluate another type of indirect method which will 1. rapidly obtain an approximate solution for possible real-time implementation; 2. identify the optimal switching structure for constrained problems; and 3. produce initial guesses for a shooting method. This chapter describes a finite element method based on a weak formulation of the first-order necessary conditions of optimality to achieve these three goals. The remainder of this section will mathematically define the optimal control problem. Then, Section II develops the weak principle for optimal control theory for unconstrained problems. Two examples are given to illustrate the use of the weak principle. In one example, the states and controls are continuous. The second example involves a problem with discontinuities in the states and system equations. Further extensions of the weak
A FINITE ELEMENT METHOD
189
principle are given in Sections III and IV with the inclusion of control and state inequality constraints. Several example problems are given in these two sections to highlight the behavior of the finite element method. Some features of the weak formulation and conclusions are discussed in Section V. C. O P T I M A L C O N T R O L PROBLEM D E F I N I T I O N Consider a system defined over a time interval from the initial time to to the final time t I by a set of n states, x, and a set of m controls, u. The states of the system are governed by a set of first-order differential equations referred to as state equations. During the time interval from to to t f, there may be discontinuities in the states, or discontinuities in the state equations at interior points (times between to and t l ) , or control and/or state constraints to satisfy. These interior points, along with the initial and final points, will be referred to as events, and the time interval between events will be referred to as phases. The time of event i will be denoted with ti and times just before and after ti will be denoted with t~and t +, respectively. All event times will be considered unknown except that of the zeroth event, coinciding with the initial time to. Elements of a performance index may be denoted with an integrand L(i)[x(t), u(t)] over each phase i and a discrete function ¢ of the states and/or times at any of the events. A general class of such problems with N phases involves choosing u(t) to minimize J
--
¢[x(t+o),x(tl),X(t+l),...,X(tN);tl,t2,...,tN]
+
~
L(i)[x(t),u(t)] dt
iX'~
(1)
¢--1
subject to the state equation constraints
Jc - f(i)Ix(t), u(t)];
(t+_l < t < t~-; i - 1 , . . . , N)
(2)
with boundary conditions on the states and times specified as
),...
,tN] - o
(3)
We introduce Lagrange multiplier functions )~(t), referred to as costates, and discrete Lagrange multipliers v. Then, for convenience, (I) and H are defined as
¢
-
¢+~¢
H(i)
_
L(i) + ATf(i)
(i -- 1 , . . . , N )
(4) (5)
190
ROBERT R. BLESS A N D D E W E Y H. H O D G E S
The augmented cost function obtained by adjoining Eqs. (2) and (3) to Eq. (1) using the costates A(t) and multipliers u, respectively, is
Ja = ~ + ~~1=
-1 H(~) _
)~T~
) dt
(6)
The derivation of the first-order necessary conditions for optimality begins by expanding the increment of the augmented cost function of Eq. (6) (i.e., the difference between the cost function evaluated at a solution neighboring to the optimal and the cost function evaluated at the optimal solution) into a Taylor series about the optimal solution. The first-order term of the Taylor series is referred to as the first variation. Provided that the variations to the optimal solution are small enough so that the sign of the increment of the cost function is determined by the first variation, then it is well-known [45] that the optimal solution furnishes a stationary point to Ja. The first variation is ;
+ 6uT
Ou
dt + du T¢
(7)
N-1
+
Z i=1
+
-~i + L(i) (t:~ ) -
(t +) dti + -~N + L(N) (tN) dtg
[ ]
Zi =~1 Ox(t;)°~ dx(t;)+ Zi--0 Ox(t+)
where 6x(t) signifies the perturbation of x when t is considered fixed and dx(t +) and dx(tT~) signify the perturbation of x due to perturbations of ti, as well as due to changes in x at t + and t~-, respectively. These are related as
dx(t +) dx(tT ) -
6x(t +) + 2(t + )dti 6x(t 7) + ~(t 7 )dti
(8) (9)
Note that since x can be discontinuous at ti, so can the perturbation of x. After integrating the 6~TA term by parts and using Eqs. (8) - (9) to eliminate 6x(t +) and 6x(t:~) in terms of dx(t +), dx(t:~ ), and dti, we obtain
i=1
,-
A FINITE ELEMENT METHOD
5uT
OU
191
(10)
dt + duT¢
N-1
+
E i=1
-~i + H(i)(t[) -
(t +) dti + -~N + H(N)(tN) N-1
0e
Ox(t+~------~+
dtN
~ (t+)] ex(t+ )
The first-order necessary conditions of optimal control require that 5Ja - 0
(11)
for all independently chosen perturbations dtl, dt2, ..., dtN, du, dx(t+), dX(tl) , dx(t+), ..., dX(tN) , 5x(t), 5,~(t), and 5u(t). Hence, the discrete boundary condition terms are c0(I) H(i+I) Ot---~+ H(O (t;) (t +)
-
0
O0 Ot---N +
H(N)(tN)
-
0
(13)
~[x(t0+),x(t 1), . . . , x(tN); tl, . . . , tN]
--
0
(14)
-
0
(i- 1,2,...,N)
(15)
+ ~(t+~ ) -
0
(i - 0, 1 , . . . , N -
1)(16)
0O o~(t;)
OO
o~(t+~ )
-
;~(t~-)
(i- 1,2,...,N-
1)(12)
where Eq. (12), for instance, is obtained by setting all perturbations equal zero except for dti, which in turn forces the coefficient of dti to be zero. The Euler-Lagrange equations that come from Eq. (10) for each phase i are :~_f(0 (OH(i))
~+
=
0
(17)
-
o
(18)
=
0
(19)
T
Ox
OH(O Ou
Equations (12) - (19) define a multi-point boundary-value problem. Since we have only considered necessary conditions, not sufficient, the solution of these equations may yield an optimal value of the cost function. In Section II, a weak formulation of these equations is derived to obtain an approximate, candidate optimal solution.
192
ROBERT R. BLESS AND DEWEY H. HODGES
II. UNCONSTRAINED PROBLEMS In this section, a variational formulation for solving unconstrained optimal control problems is developed. It will be shown that the formulation has no strong boundary conditions (i.e., boundary conditions requiring the virtual, or variational, quantities to be zero) but only natural, or weak, boundary conditions (i.e., those determined by setting the coefficient of a virtual quantity to zero). This allows the test functions (which approximate the virtual quantities) and shape functions (which approximate the dependent variables) to be chosen from a less restrictive class of functions; indeed, the same test and shape functions may be chosen for every optimal control problem. Thus, a general set of algebraic equations will be developed which do not have to be altered in terms of numbers of equations with any changes in the boundary conditions. Hence, this formulation is also referred to as a weak formulation. A. GENERAL DEVELOPMENT The derivation of the weak formulation begins by using the conditions in Eqs. (12) - (16) to simplify Eq. (10), yielding
+ 5uT
Ou
dt
(20)
To put Eq. (20) into the weak form for the finite element discretization, we will integrate by parts so that no time derivatives of x or ~ appear in the weak formulation. As will be shown, this allows the simplest possible shape functions to be chosen and for all element quadrature to be done in closed form. Integrating by parts yields
5uT
OU
dt + E [5xT(t:~ )£(t:~ ) -5£T(t~)x(t~)] (21) i--1
+
N-1
Z [ ~ (t~+)~(t~+) - ~ (t~+)~(t~+)] i=0
A FINITE ELEMENT METHOD
i to
--
At(1 ) t (1)
I
t (2)
tl
II --
I
t (3)
193
I
t (5)
t (4)
I
t2 -- t (6)
Figure 1" Discretized time line with nodes labeled Equation (21), when taken with Eqs. (12) - ( 1 6 ) , is called the weak formulation of the first-order necessary conditions for the optimal control problem defined in Eqs. ( 1 ) - (3). Equations ( 1 2 ) - (16), and Eq. (21) will be used for the finite element discretization scheme described next. B. F I N I T E E L E M E N T DISCRETIZATION
i-
Let the time interval from t +_ 1 to t~- be broken into Mi elements where 1, 2 , . . . , N. For convenience, define i
m-ZMj
for
i-
1,2,...,N
(22)
j=l
and define M0 = 0. This yields a subdivision of the original time interval from to to t f into M N subintervals. The boundaries of these subintervals are called nodes and are denoted by t (i) for i - 1 , 2 , . . . , M N + 1. Note that to = t (1), ti = t (M{+I), and tl - tN : t (MN+I). A nondimensional elemental time T is defined as T
=
t - t (i) t(i+1) - t(~)
=
t - t (i) At(0
(23)
SO that 0 <_ T < 1 within each element. Figure 1 gives a pictorial representation of a two-phase time line with M1 = 2 and M2 - 3. Since no derivatives of x, A, or u appear in Eq. (21), then the simplest possible shape functions in each element, namely piecewise constants, may be used. The benefit of choosing piecewise constant shape functions is that element integration may be done by inspection. The shape functions are
X-
{
~+ if T -- 0; ~i if 0 < T < 1; X~-+I i f T - - 1
(24)
194
ROBERT R. BLESS AND DEWEY H. HODGES
^
~1+
X5
I t (1)
t (2)
1
t (3)
I ~
I t (1)
t (4)
II x~
t(2)
t(3)
I x4
t(4)
I
I
t (6)
t(5)
I
A5
t (6)
t(5)
^
~1+
I t(I)
_
x6
~5
Ul
t (2)
t (3)
t (4)
t (5)
_
u6
I
t (6)
Figure 2: Location of unknowns for x, A, and u
{ {
Aand u-
A+ A~ ~-+1
if "/- - 0; if 0 < T < 1; if T - - 1
(25)
~.+ fi~
if ~- - 0; if0
(26)
Iti---F1
Since first-order time derivatives of 5x and 5A appear in Eq. (21), linear shape functions are chosen for 5x and 5A in each element. Since no time derivatives of 5u appear, then piecewise constant shape functions may be chosen in each element. Again, by choosing the simplest functions allowed, element integration is done by inspection. The functions chosen were 5x
-
5x+(1 - T)
+ (~XhlT
(27)
5A
--
5 A + ( 1 -- T) + (~i-+l T
(28)
and, with 5D defined as the Dirac delta function ~U -- ~U~D(T)-t-~Ui
nL- (~U~+I(~D( T -- 1)
(29)
A FINITE ELEMENT METHOD
195
(~Xl+
(~X2
(~X3 (~X3+
(~X4
(~X5
(~X6
t(l)
t(21
t(a)
t(4)
t(5)
t(6)
I
I
t(1)
I !
II
II
I
I
I
I
I
I
t(2)
t(3)
t(4)
t(5)
t(6)
t(2)
t(3)
t(4)
t(5)
t(6)
+
t(1)
Figure 3: Location of independent variational quantities for 5x, 5A, and 5u
The superscripted "-" and "+" signs in Eqs. (24) - (29) signify values just before and after a particular nodal time t(0. For all nodes except event nodes, the values for x and A, as well as for 5x and 5A, are equal on either side of the node. In other words, only at an event are the states, costates, and their variations allowed to "jump" in value. Thus any discontinuities in the states or costates must be known in advance and modeled with an appropriate event and phase. The examples throughout the chapter will clarify the specifics of defining multi-phase problems. The Dirac delta functions appearing in the discretization of 5u have the effect that the coefficient of the ~u term in the integrand of Eq. (21) is forced to zero pointwise wherever the Delta function has a nonzero value. More specifically, the Dirac delta functions enforce OH/Ou = 0 at the left and right limits of each node. We note that the control could be discontinuous at an internal node even though the states and costates must be continuous. For an example case with two elements in phase one and three elements in phase two, the quantities defined in Eqs. (24) - (26) that will appear in the algebraic weak form are shown on a timeline in Figure 2. The quantities defined in Eqs. (27) - (29) are shown in Figure 3. After changing the variable of integration from t to T with Eq. (23) and
196
ROBERT R. BLESS AND DEWEY H. HODGES
inserting the shape and test functions from Eqs. (24) - (29) into Eq. (21), the element quadratures can be carried out analytically yielding an algebraic equation. Noting that j indexes the phase number and i indexes the element number within a phase, the equation is 0
-+
+
N<
..[_T ( ,, ..~_ + ~A'M'j_ lq--1 -- X'M'j_ 1-1--1 :E'M'j-1-'b1
E j=l
+~ (~X " ~ j _
[~+
~+1
~j
_t..T ~£t~Mj_lnt-1
--1
+1
"~
-
E
2
"M'j- 1"[-1
M'j-1 +1
Atj -
{~)~T (--Xi-1
Atj )
i=Mj_I+2
-3I- ~XT [~i_l +
-1-t-'1
-
A~j_ 1 +1
Mj
+
2
Mj Z
i=Mj_I+I
+
&~j+l
+
(~XMj+1
e<
[(v:):] Atj 2
--~MJ 2
0/:/ -+- (SZ~Mj + 1
(30)
2- -g~- ~-1
"-~
Cv, + ^-
x
XMj +1
-gg
~,
)
- X:~,+I
+1]>
In Eq. (30), note that the superscripted "-" and "+" signs are dropped except at the event nodes where discontinuities may take place in both the state and costate quantities and the variational state and variational costate quantities. A l s o , / ~ - H(2, fi, A), f - f(~, ~), a n d / 2 / _ H(~, fi, A). Additionally, we have taken the elements within each phase j to be of constant width Atj. Finally, in the above equation, 5~+ and 5fi- terms should appear at every node and not just the event nodes; however, the 5fi terms at internal nodes are decoupled from the other ones and may be used to solve for internal nodal values of the control after the other unknowns have been solved.
A FINITE ELEMENT METHOD
197
When the shape functions for x, A, and u as defined in Eqs. (24) - (26) are substituted into the boundary conditions of Eqs. (12) - (16), then the following equations are obtained.
+'
Ot--~+ H~M i + I - /~-~i +1 O0 +[-ti)tN MN+I ^
_
^ - -
Ox(t:~ ) 0(1)
-- )~ i+I--M"
+
^~_T
=
0 (i = 1 , 2 , . . . , N -
=
o
(32)
=
0
(33)
0 (i - 1 , 2 , . . . , N )
(34)
(i = 0, 1 , . . . , N -
1)(35)
--
=
0
1)(31)
Now, in Eq. (30), the coefficient of each of the arbitrary virtual quantities, 5x, 5)~, and 5u, must be set equal to zero in order to satisfy the equation. When the coefficients are set to zero and these equations are combined with Eqs. (31) - (35), the result is a sparse system of nonlinear equations whose size depends on the number of elements. Solution of these algebraic equations yields an approximate candidate optimal solution. A summary of the number of equations and unknown is given next. C. SUMMARY OF EQUATIONS AND UNKNOWNS Consider a problem with N phases, n states (and, hence, n costates), m controls, and q boundary conditions in ¢. For any given phase i, the unknowns appearing in Eq. (30) and Eqs. (31) - (35) are the barred quantities for each element for x, ~, and u, along with the hatted (nodal) quantities at the beginning and end of each phase for x, ~, and u (see Figure 2). Thus, there are a total of (2n + m)(Mi + 2) unknowns for the ith phase. Summing these over N phases results in (2n + m ) ( M y + 2N) total barred and hatted quantities for x, A, and u. In addition to these unknowns, there are q unknown multipliers v and the N unknown times tl, t2, ..., tg. Therefore, the total number of variables to solve for is (2n + m ) ( M g + 2N) + q + Y. For any given phase i, the equations appearing in Eq. (30) and Eqs. (31) - (35) are the coefficients of each independent 5x, 5)~, and 5u when set equal to zero. Since the 5x and 5~ quantities appear at each node (see Figure 3), there are 2n(Mi + 1) equations for the ith phase. Since the 5u quantities appear at the end nodes of each phase, and also at the midpoints of each element (see Figure 3), there are m(Mi + 2) equations. The total number
198
ROBERT R. BLESS AND DEWEY H. HODGES
of equations then from phase i is m + (2n + m)(M~ + 1). Summing these over N phases results in m N + (2n + m)(My + N) total equations from the coefficients of 5x, hA, and 5u. Additionally, Eqs. (31) - (32) provide N equations, Eq. (33) provides q equations, and Eqs. (34)- (35) provide 2nN equations. Thus the total number of equations is (2n + m)(MN + 2N) + q + N, which agrees with the number of unknowns. The algebraic system of equations is typically solved by expressing the Jacobian explicitly and using a Newton-Raphson solution procedure. For the example problems which follow, the iteration procedure will converge quickly for a small number of elements with a relatively simple initial guess. Then, the answers obtained for a small number of elements can be used to generate initial guesses for a higher number of elements. When the sparsity of the Jacobian is exploited, it is possible to efficiently generate a reasonably accurate solution. Although the nonevent nodal values for xi and Ai do not appear in the algebraic equations, their values can be easily recovered after the solution is found. This is most easily seen by looking at the following ordinary differential equation multiplied by a test (or weighting) function 5A and integrated over some time interval where the integral makes sense.
ftl 6A[f (x, t) - 5] dt
(36)
0
After an integration by parts, using the linear shape function for 6A defined in Eq. (28), using the piecewise constant shape function for x defined in Eq. (24), and substituting T for t as given in Eq. (23), then the following equation is obtained from Eq. (36). ~/~1
--:~1 ~- "-~'Y nt- Xl
-{- (~/~2
:~1 -~- - ~ - f -- :E2
-- 0
(37)
With arbitrary ~A1 and ~A2, the coefficients must vanish, forming two equations of the form AtAt-
-
0
-
0 o
Now, by subtracting the second equation from the first, it is seen that _ = ~1 + ~2 Xl
2
(38)
or, in words, that the interior value (the bar value) is simply the average of the surrounding nodal (or hatted) quantities. Once a solution to the
A FINITE ELEMENT METHOD
199
algebraic equations is found, all the midpoint values and the end nodal values are known, and thus all other nodal values can be recovered by repeatedly using Eq. (38) for both x and A. In view of the redundancy, only the nodal values are plotted for the states and costates in this chapter. As for the control, interior nodal values are available once the nodal values for the states and costates are found by using the optimality condition (OH/Ou = 0) to solve for u at nodes. Recall that these equations appeared in the weak formulation but were decoupled from the other equations. We note that the midpoint value for the control is not the average value of the surrounding nodal values; thus, both midpoint and nodal values are plotted in the figures in this chapter. D. EXAMPLE: A SINGLE-PHASE P R O B L E M As the first example problem, the transfer of a particle to a rectilinear path will be examined. This example is taken from section 2.4 of [4]. The purposes of this example are to (1) explicitly state the boundary conditions and algebraic equations of the weak formulation for a one-phase problem, and (2) present some results which show the accuracy of the method for this problem. Let x(1) and x(2) denote the horizontal and vertical position coordinates, respectively, of the particle at a given time and x(3) and x(4) denote the particle's velocity components at a given time. (A subscripted number in parentheses refers to the state index to avoid confusion with the element index.) The thrust angle u is the control and the particle has mass m and a constant acceleration a. The state equations are defined as X(3) --
X(4)
a cos u a sin u
= f
(39)
The problem is to maximize the final horizontal component of velocity in a fixed amount of time t l = T. Thus, L = 0
and
¢ = x(3)(tl)
(40)
As for boundary conditions on the states, the initial conditions are zero and there are two terminal constraints. These are that the particle arrive with a fixed final height h and that the final vertical component of velocity be zero. The final horizontal component of position is free.
200
ROBERT R. BLESS AND DEWEY H. HODGES
For single-phase problems, Eqs. (30) - (35) are greatly simplified. Note that since there are no internal events, the "+" and "-" superscripts may be dropped. Looking at the boundary conditions first, Eq. (a3) becomes x(1)(to)
x(2)(to)
x(3) (t0) x(4) (to)
=0
x(2)(tl)-h
(41)
x(4) (tl) tl - T
where the conditions in ¢ are adjoined to the cost with multipliers Vl, /22, ..., vT. Equation (31) is no longer applicable and EQ. (32) becomes ^
/27 + H~'a+ 1 - 0
(42)
The costate boundary conditions in Eqs. (34) and (35) become 0
/21 1 nt-
-- 0
/22 /23 /24
^ ~+1
and
/25 1
-
-0
(43)
/26
Finally, Eq. (30) is simplified to
0
=
e~l~
q-
/22"= ~ / T
+
exT
i--1
-~1 + ~1
(
~-fl
At-
Xi-1----2-
At
I~-1 + ~ - - V ~
--:~/-1
-~
i-1
-
(
at ( 0 ~ )
)
+ &xT~+I --~M~ ~-fM~+ a?Mx+l
+ &5~+1 XMI-- 5- -0~ M1--
)
'~MI+I]
(44)
A FINITE ELEMENT METHOD
+
~uTI q-i
-~
201
Mlq-I
where At = (tl - to)/M1 = tl/M1 = T/M1 for all elements. These equations are solved by expressing the Jacobian explicitly and using a Newton-Raphson algorithm. For M1 - 2, suitable initial guesses for the nonlinear iterative procedure can be found by simply choosing element values that are not too different from the boundary conditions. The results from solving the M1 = 2 equations are then used to obtain the initial guesses for arbitrary M1 by linear interpolation. Representative numerical results for all four states versus dimensionless time t / T are presented in Figures 4 - 7. For this example, h = 100, T = 20, and 4 h / a T 2 - 0.8897018. This last number is chosen to yield a value of 75 ° for the initial control angle of the exact solution available in [4]. The results for 2 and 8 elements are plotted against the exact solution. It can easily be seen that M1 -- 8 gives acceptable results for all the states. Even the very crude 2-element mesh yields a decent approximation to the answer and led to convergent initial guesses for the 8-element run. In Figure 8, the thrust angle u versus dimensionless time t / T is presented. Once again, the results are seen to be excellent for M1 - 8. Note that instead of just nine data points (circles) for the nine nodal values which were plotted for the states, there are 17 data points available for the control since both midpoint and nodal values are obtained from repeated use of the optimality condition. This is of great value since it is the control variable which is generally of the most interest. Three of the four costates are constants for all time and this method yields two of these exactly. The third costate is very close to the exact answer. The fourth costate, A(4), corresponding to the vertical component of velocity is shown in Figure 9. The results compare nicely with the exact results. A plot of the relative error of the performance index J - x(3)(T) and the initial control angle u(0) - Ul versus the number of elements is shown in Figure 10. For 8 elements or more, the lines are nearly straight on a loglog scale. The slope of the line is about - 2 which indicates that the error varies inversely with the square of M1, similar to a-posteriori error bounds as formulated in usual finite element applications [46]. Notice in Figure 10 that there is a bend in the initial thrust angle curve. It is not unusual for mixed formulations to have an error curve that is not monotonically decreasing. It should be noted that developments of mathematical proofs for convergence and expressions for error bounds are not state-of-the-art for mixed methods.
202
ROBERT R. BLESS AND DEWEY H. HODGES
1.2-
1.0c-
FE(2)
•
FE(8)
.
/
Exac
C i
•
O n
{D
0.8-
O 0_
-
n
C O N O
7"
0.6_ O.4-
0.2-
0.0
q
0.0
0.2
0.4
0.6
0.8
Dimensionless Time Figure 4" Normalized horizontal position versus
t/T
1.0
A FINITE ELEMENT METHOD
203
1.00.90.8-
v ,,
FE(8)
-
-
Exact
0.7 -
P,
:,~
FE(2) /
J
~
0.6-
13_ 0.5o4>
0.3 0.2 0.1
O0 1 0.0
'
I
0.2
'
I
0.4
'
I
'
0.6
I
0.8
Dimensionless Time Figure 5: Normalized vertical position versus
t/T
'
I
1.0
204
ROBERT R. BLESS AND DEWEY H. HODGES
3.0-
2.5-, O
• •
~'~ 2.0-
~
FE(2) FE(8)
//~ ~
•
1.5-
N
1.0-
-r"
0.5-
0.0
l"
I
0.0
'
I
0.2
'
I
0.4
'
I
0.6
'
I
0.8
Dimensionless Time Figure 6: Normalizedhorizontal velocity versus t/T
'
I
1.0
Vertical Velocity/(h/T) 0"q
0
cb
0
--,,1
0
0
~..,o
0
CI)
CL .,¢ c--F C~
o o I.-,L"F
.<
C'D
CI)
o
b0
0
o ±
0
~
~
0
I
J
~
0
I
i
~
0
I
~
b0 I
~
J
o
~
I
t
ro I
~
;
.~
I
206
ROBERT R. BLESS AND DEWEY H. HODGES
80-
60-
40-
v
20-
m
O
m
-20 -
-40 -
•
FE(8) Exact
~
-60 ip
-80
'
0.0
I
0.2
'
I
0.4
'
I
0.6
'
I
0.8
Dimensionless Time Figure 8: Thrust angle versus t/T
'
I
1.0
A FINITE ELEMENT METHOD
4
~
207
!
3
2 0
1
"~
o
-2
-3
9
-4
'
0.0
I
0.2
'
I
0.4
'
I
'
0.6
J
0.8
Dimensionless Time Figure 9: Vertical velocity costate versus
t/T
'
I
1.0
208
R O B E R T R. B L E S S A N D D E W E Y H. H O D G E S
100 -
1 0 "1 - _
0
I,..
LLI .
(1) >
.
Cost
1 0 -2 - -
m
n
rc 1 0 .3 - -
Initial Thrust Angle
10 .4
I
1
I
I
I
I
I
I I ]
1
I
I
I
I
I I I I
10
Number of Elements Figure 10: Relative error for cost and initial thrust angle versus number of elements
100
A FINITE ELEMENT METHOD
209
E. EXAMPLE" MULTI-PHASE ASCENT TO ORBIT A two-stage launch vehicle will now be examined with the three purposes of (1) providing an example of a multi-phase problem definition, (2) examining the finite element accuracy for highly nonlinear problems, and (3) more closely examining how the Newton method behaves when using linear interpolation to increase the number of elements. Consider the following model of a rocket flying over a spherical, nonrotating Earth with the states chosen to be mass m, altitude h, speed V, and flight-path angle 7, and the control to be the angle of attack a. The dynamical equations are approximated to be
Tvac gI~p
=
h l)"
-
(45)
ysin7 Tcosc~-D m Tsina+L mV
=
(46) #sin7 r2 +
(V
r
(47) it ) r2V cos 7
(48)
where T - Tva c - A e p , T,a~ is the thrust in a vacuum, Ae is the exit area of the nozzles, p is the local atmospheric pressure, Isp is the specific impulse, g is acceleration due to gravity at sea level, # is the Earth's gravitational constant, and r is the distance from the center of the Earth to the vehicle, given by the radius of the Earth Re + the altitude h. Defining p as the density, S as the reference area and Ca and CN as the axial and normal force coefficients, respectively, then the dynamic pressure, q, the drag, D, and the lift, L, are defined as
q-
1 -~pV 2
(49)
D
-
aS(Ca 2t- CNOL 2)
(50)
L
-
q S ( C N - C~)a
(51)
The numerical values for the physical constants are # - 3.986 x 1014 m3/s 2, Re = 6378000 m, and g - 9.81 m/s 2. The propulsive model used was
Tvac
=
1.2975 z 107 N before staging
(52)
Tvac Ae
=
2.5950 × 106 N after staging
(53)
=
19.115 m 2 before staging
(54)
210
ROBERT R. BLESS AND DEWEY H. HODGES
Ae
=
3.823 m 2 after staging
I~p
=
430.55 s
(55) (56)
and the atmospheric and aerodynamic model was kg/m 3
(57)
p
=
1.225exp(-h/8600)
p
=
101325exp(-h/7600) Pa
(58)
S
=
55.18 m 2
(59)
Ca
--
0.3
(60)
CN
--
3.1
(61)
We note that constant values for the aerodynamic coefficients are unrealistic; however, this assumption does simplify the numerical operations needed to solve the problem. Also, for numerical simplicity, all aerodynamic and atmospheric terms were neglected after staging, since they are extremely small by this point in the trajectory. The performance index is defined as J = ¢ = - m ( t f ) , and the final time tf is free. Defining m0 = 890150 k.g as the initial mass, hf = 148160 m as the final altitude, and VI = 7854 m/s as the final speed, then all the state and time boundary conditions required for the finite element formulation can be stated as re(to)
-
mo
h(to) - 60
V(to)-25 ")'(to)- 1.5 re(t+) - r e ( t 1 ) +
¢ -
29920
h(t+) - h ( t l )
- 0
v(t+l)- V(tl) ~(t+l ) - ~ ( t l ) h(t~) - h~ v(t~) - vI ~(t~)
(62)
tl -- 195 Note that (1) t l is the specified time of staging, (2) there is a discontinuity in mass of 29920 kg at staging, and (3) continuity of the other three states is explicitly stated. By adjoining the conditions in ~ with multipliers ui for i = 1, 2 , . . . , 12, then Eqs. (31)and (32) become
~1~ + / t ~ + 1 -/~+~,+1
= 0
(63)
--
(64)
^
H-~2 + 1
0
A FINITE ELEMENT METHOD
211
The costate b o u n d a r y conditions in Eqs. (34) and (35) become Pl
+ +
u2
- 0
/23
and
--/]5 ~Y6
^
A~+I
-
-
P4 and
^ q-
nl-
-
0
(65)
= 0
(66)
--/]8
u5
AMI+I
-
--/27
/26
u7
-1 --
0
and
^ A~M2+ 1 --
/28
v9 Ul0 /211
Since the internal boundary condition was based on a time of staging instead of a particular state value (such as mass), there are no jumps in the costates. An example with multiple stages and jumps in the costates will be given in Section IV. For numerical conditioning, the states were normalized in a simple manner by defining ~ - m/mo, h - h / h f , and V - V / VI. Normalizing the states automatically scales the costates. The numerical results which follow were generated with the code VTOTS which is described in [16]. Obtaining guesses for a small number of elements is generally easier than for a large number of elements. Additionally, once a solution has been found for a few elements, more dements may be added readily in most cases. For example, Table I was produced after a 2-element-per-phase solution was found (a very crude approximation). By linearly interpolating the solution, new initial guesses were computed for a 4-element-per-phase solution. These new guesses yielded an average initial error of the algebraic equations of 2.60015 x l0 -2, as displayed in the second column. The last column of Table I displays the total number of iterations, using a Newton method, required to obtain a converged solution. The "4:4" solution was then used to obtain initial guesses for the "8:8" case, and similarly, the table was completed. It is seen that the initial error monotonically drops, as does the number of iterations required for convergence. Doubling the number of elements has proven to be a quick and easy way of obtaining reasonably accurate answers. Table I. Build up of finite element values Elements Per Phase 4:4 8:8 16:16 32:32
Average Initial Error 2.60015 x 10 -2 9.93419 x 10 -3 2.35328 x10 -3 4.81494 x 10 -4
Number of Iterations 5 4 3 ' 3
212
ROBERT R. BLESS AND DEWEY H. HODGES
Some of the results are shown in Figures 11 - 18. In all of the figures, finite element results and shooting results are shown. The finite element results are taken with 2 elements per phase and 8 elements per phase and are denoted with "FE(2:2)" and "FE(8:8)", respectively. Figure 11 shows the normalized altitude profile, and Figure 12 shows the normalized velocity. It is seen that the 2-element case only gives a rough approximation to the answer, and that the 8-element case lies essentially on top of the shooting curves. Figure 13 shows the flight-path angle, V, versus time. The vehicle is seen to roll over quickly to take advantage of the gravity field. Notice that the finite element results show some error for this case. The error is due to an insufficient number of elements being used to model the rapidly varying 7 profile at the initial time. A finer mesh increases the accuracy in this area. Figures 14 - 17 give the four costate histories for the vehicle. The finite element results behave similarly to the state histories; namely, that the 2-element results are a rough approximation, and that the 8-element results are quite accurate except in areas where the variables change more quickly than the finite element mesh is able to model. Examples of this are seen in Figures 15 and 17 for the altitude and flight-path angle costate, respectively. The angle of attack history is shown in Figure 18. The control experiences a very rapid change from a negative angle of attack to a positive one at the initial time. The 2-element case yields ridiculous results initially with an initial angle of attack of -81.9°! However, since the rest of the 2element control history is reasonable, it is still possible to use that answer to obtain guesses for a finer mesh. The 8-element case is seen to be much more accurate with an initial control angle o f - 2 3 . 5 ° , however, it is still not very accurate in the first stage where the control varies rapidly. Figure 19 shows the relative error of the midpoint values of control for varying meshes. In all cases, we have used 16 elements in the second phase, and varied the number of elements in the first phase. It is seen that the control error typically decreases with more elements, although there are singular places where this is not true; and that the largest error is at the initial time. These examples conclusively demonstrate the accuracy and efficiency of the present methodology for unconstrained problems. The method provides virtual self-starting capability and a rapid convergence of states, costates, and controls as the mesh is refined. These are both important for potential real-time applications. III. T R E A T M E N T
OF CONTROL
CONSTRAINTS
This section is concerned with the treatment of optimal control problems
A FINITE ELEMENT METHOD
213
-
1.0~
.~
v,
~-
0.8
.~~ o.6 N ¢
0.4
:2) FE(8:8)
•
t~ng
Z 0.2
0.0
"
._...~ I
I
100
'
I
200
'
I
300
Time (s)
Figure 11" Normalized altitude versus time
'
I
400
'
I
500
214
ROBERT R. BLESS AND DEWEY H. HODGES
1.00.90.8•~,
a
m
0
0.70.6-
"~
0.5-
~
0.4-
Z
0.3-
N
-
•
FE(2:2)
•
FE(8:8) Shooting
0.20.10.0 T 0
'
I
100
'
U
200
~
I
300
Time (s) Figure 12" Normalizedvelocityversus time
'
I
400
'
I
50O
A FINITE ELEMENT METHOD
215
90-
80-
70-
ED (D
60-
_.m
50-
v
ED c"
<
J.4=,.J :
Q. !
.
tED
.
m
40FE(2:2)
30-
FE(8:8) Shooting
20-
m
LL 10m
-10 0
'
I
100
'
I
200
'
I
300
Time (s) Figure 13: Flight-path angle versus time
I
400
I
50O
216
R O B E R T R. B L E S S A N D D E W E Y H. H O D G E S
o
-0.1 -~ -0.2 -0.3
(D a3
r.D
O
O
-0.4 -0.5
{D
09
-0.6 FE(2:2) FE(8:8) Shooting
-0.7-0.8-0.9-1.0
'
0
I
100
I
I
2O0
'
I
I
3OO
4O0
Time (s) Figure 14" Mass costate versus time
I
v
I
5O0
A FINITE ELEMENT METHOD
217
0.0-0.02 -0.04 -0.06 -0.08
0
-o.1
"0 :3
-0.12-
<:::
-o.14-
•
FE(2:2)
•
FE(8:8) Shooting
-0.16-0.181
-0.2 '
0
I
100
'
I
200
I
I
300
Time (s) Figure 15: Altitude costate versus time
i
i¸
I
400
5OO
218
ROBERT R. BLESS AND DEWEY H. HODGES
-0.1 -0.15 -0.2
(D
{D O
(O m
m
-0.25 -0.3 -0.35
•
O (D
>
FE(2:2) FE(8:8)
•
m
Shooting
-0.4 -0.45 -0.5 -0.55 -0.6
I
0
'
I
100
'
I
200
'
I
300
Time (s) Figure 16: Velocity costate versus time
'
I
400
I
5OO
A FINITE ELEMENT METHOD
219
.10 -~ 1.0-
0.5-
FE(2:2) FE(8:8) Shooting
0.0
0
(1)
-o.5
~C:D -1.0
t'3 !
-1.5-
c" C:)') •--2.0IJ_ -2.5-
-3.0 0
'
I
100
'
I
200
'
I
300
Time (s) Figure 17: Flight-path-angle costate versus time
'
I
400
'
I
500
220
ROBERT R. BLESS AND DEWEY H. HODGES
20-
O
m
"0
v
0
-20 -
< 0
FE(2:2) FE(8:8) Shooting
-40 -
m
c-
<
-60 -
-80 -
i
0
I
IO0
i
I
20O
i
I
3OO
Time (s) F i g u r e 18" A n g l e of a t t a c k v e r s u s t i m e
I
400
I
5O0
A FINITE ELEMENT METHOD
221
100 8 Elements I
1
0
"1
t
-_
"",,,
%
~ \
",,
%%
\
_
16 E l e m e n t s
"',,,
64 Elements
',
\ ,,
_
~L_
\ ,,
0
UJ (D >
32 E l e m e n t s
""
\ \, o° i \
1 0 "2 - _
. , m
_
\
,s
\ %~ \
\
\
\
\
%%% ",~,% i - - - - . . .
°°
s "
,~
I
\
I
m
(D rr
oo
oo
\
\
,"
l 1 0 -3 - -
\\
I I
10"
' 0
I 50
'
I 1 O0
\
'
I 150
Time (s) Figure 19: Relative error in angle of attack versus time
'
\\\
I 200
222
ROBERT R. BLESS AND DEWEY H. HODGES
that are subject to inequality constraints which are explicit functions of the controls. A. GENERAL D E V E L O P M E N T Consider the same system as defined in Eqs. (1) - (3). For simplicity, it is assumed that the problem has only one phase. Suppose now that the p constraints are written as
gi(x, u) ~_ 0
for
i = 1, 2 , . . . , p
(67)
One way of handling inequality constraints is to use a "slack" variable [47]. The idea is that if gi _< 0 then gi plus some nonnegative number, the ith slack variable k 2, is equal to zero. Thus, we form p x 1 column matrices for K and 6K
K6K
Lk~ k~ =
[2~1 (~kl
...
k]j T
2]{2(~2
...
(68)
2kp(~kpJT
(69)
and a p x 1 column matrix g(x, u), the elements of which are gi(x, u). Now, from Eq. (67) g(x, + K : 0 (70) Eq. (70) will be adjoined to the augmented cost function Ja in Eq. (6) by using p Lagrange multiplier functions #(t). After dropping the superscripted phase number on H, we obtain
J~ - • +
fl 1 [ H -
IT9 + #T (g + K)] dt
(71)
By redefining the Hamiltonian in Eq. (5) to be
(72)
H =_ L + )~Tf + #T(g + K)
then, after the usual integration by parts and making use of Eqs. (8) - (9), we obtain the first variation of Ja as
61T (f -- x) + ~xT
--
+
5u T
~
-~x
+spT(g+K)+6KT#
+ i
dt
(73)
A FINITE ELEMENT METHOD
+
+
223
~ 1 + H ( t l ) d r 1 + dp T ¢ [ 0(I)
OX(tl)
)~T(tl)]
dx(tl) + Ox(to) + AT(t°) dx(to)
Setting 5Ja = 0 in Eq. (73) results in the same necessary conditions (for a single-phase problem) as are listed in Eqs. (12) - (19). There are also the additional conditions that gi+k 2-0 #iki-O
~ J
for i - 1
'
2,
"" " P
(74)
These conditions, when satisfied, ensure that the constraint will not be violated at element midpoints and nodes and that either the constraint or the multiplier will be zero. For example, if the constraint gi equals zero, then ki must be zero, and so #i need not be zero. On the other hand, if gi is less than zero (inactive), then ki must be nonzero, which in turn forces pi to zero. An additional requirement for a minimizing problem is that the multipliers p must be nonnegative for all time. This is enforced computationally, as will be shown in the example, by squaring the value of p that appears in the costate equations and optimality conditions. There may be numerical problems associated with the introduction of slack variables. As pointed out in [1], there can be poor conditioning or even singularities in the Jacobian occurring during the solution procedure if both #i and ki were zero at the same time. The advantage of using a slack variable, however, is that control-constrained problems may be treated as single-phase problems. Since the control is allowed to come on and off the constraint as required, there is no need to define a switching structure (i.e., the sequence of constrained and unconstrained arcs) for the problem. To derive the weak formulation, Eq. (73) is simplified by using the boundary conditions in Eqs. (12) - (16). Then an integration by parts to eliminate any time derivatives of x or A from appearing in the formulation yields
0
--
+ +
5)~Tf + 5~Tx + 5xT
5u T
~
~X
+spT(g+K)+SKT# - e,
r(tl)z(tl) +
-- (~2T/~
dt
(75)
ar(to)z(to) - e z r ( t o ) a ( t o )
Equations (12)- (16) and Eq. (75) define the weak formulation for optimal control problems with control constraints. Equation (75) can be generalized in a straightforward manner to multiple phases, similar to Eq. (21).
224
ROBERT R. BLESS AND DEWEY H. HODGES
The discretization of the weak form is similar to what is defined in Section II, B.
B. FINITE ELEMENT DISCRETIZATION To the discretization defined in Section II, B, we define shape functions for 5K and 5# of the form (with 5D defined as the Dirac delta function) 5K
-
2k+Sk+~D(T)+ 2ki6ki + 2ki--+l(~ki--+l(~D(T- 1)
~lZ
--
6ft+6D(T)
nt- ~fti + (~i-+l(~D(T
--
1)
(76)
(77)
where the subscript i now refers to the element number. The Dirac delta functions appearing in the discretization of 5K and 5# have the effect that the associated parts of the integrand of Eq. (73) are forced to zero pointwise wherever the Delta functions have nonzero values. The shape functions for K and # are of the form K-
]c + ~2
if T = O; if 0 < ~- < 1;
(78)
"2--
ki+l if W-- 1 and
if # =
/~ fti-+l
-o;
if 0 < T < I;
(79)
ifT--1
The superscripted "-" and "+" signs in Eqs. (76) - (79) signify values just before and after nodal times. By substituting the shape and test functions defined in Eqs. (76) - (79) and the shape and test functions defined in Eqs. (24) - (29), and carrying out the element quadrature over 7 from 0 to 1, the general algebraic form of the weak principle is obtained. The equation is similar to Eq. (44) with the addition of the appropriate 6# and 6K equations. These equations will be written out explicitly for the example problem. Note that if the time t does not appear explicitly in the problem formulation then all integration is exact and can be done by inspection. If t does appear explicitly, then t may be approximated by its midpoint value over each element and the integration may still be done by inspection. This is equivalent to redefining t to be a new state. The latter case occurs in the example problem presented next. In addition to the unknowns and equations listed in Section II, C, there are an additional 2p(M1 + 2) equations corresponding to the coefficients of
A FINITE ELEMENT METHOD
225
~p and 5k and an additional 2p(M1 + 2) unknowns to solve for the k, ]~, #, and/2 quantities. Thus, there are still the same number of equations as unknowns. Again, once the nodal values for the states and costates are found, then one may use the optimality condition (OH/Ou = 0), the constraint equations (g+K = 0), and the condition that either k or p be zero (k# = 0) to solve for u, K, and p at a node. This procedure is used in the following example problem. C. EXAMPLE This example is taken from [4], section 3.8. The problem is to minimize J-~
1
x(T) 2 +~
1 fT
]0
u2dt
(80)
where T = 8, x and u are scalars, and the initial condition is x(0) - ( 5 + 8/~). The state equation is
2. = h(t)u
h(t) - v/2sin
with
-~-
The following two control inequality constraints are imposed. gl
=
u-1
g2
--
-(u+
50
(81)
1) _< 0
(82)
The exact solution is found to be x(T) - - 1 and for the control
-x(T)h(t) 1
u(t)-
-x(T)h(t) -1
-x(T)h(t)
for0_
(83)
The algebraic equations which come from the weak principle can be verified to be Xl -+- 5 -t- 8/71"
--
0
(84)
t l -- 8
--
0
(85)
for the state and time boundary conditions (in ¢), ~1 q-/Jl
--
0
(86)
226
ROBERT R. BLESS AND DEWEY H. HODGES ~1 -- X1
--
0
(87)
Xi -- Xi+I
--
0
'~M1 -- '~MI-t"I
--
0
(89)
'~MI-}-I -- X M I + I
--
0
(90)
for i - 1, 2 , . . . ,
M1
-
1
(88)
for the costate equations (obtained from the coefficients of 5x set equal to zero) and boundary conditions, ~ti "-~ ~ih ( t i ) + fitli -- fit2i ~-- 0
k1~#1~ - 0
k+#2~ = 0 fii - 1 ~2 _ 0 -~i-l+k!: 0
for i - l , 2 , . . . , M ~
(91)
for the coefficients of ~fi, 5k, and ~#, (note that the equations formed from setting the coefficients of 55, 5k, and ~ at the first and last node of the phase equal to zero are similar to these) and :~1 X_i + l
_ -- Xi
At
~--h (t-l)Ul
-
Xl
(92)
A t [h (t-i)ui + h ( t i + l ) Ui+I]
--
0
(93)
-
1,2,...,M1-1
--
0
for i At --:~M, -- --~-h (t-M1)UM, ~- :EMI-t-1
(94)
for the 5A coefficients. Note that t-i is an average time value for the ith element and (if At - Ati - tl/M1 for all i) can be expressed as 2i- 1 t-i = ~ A t 2
for i - 1, 2 , . . . , M1
(95)
Recall from [4] that one of the additional necessary conditions for problems with control constraints is that the multipliers be nonnegative for a minimizing problem. Therefore, in practice, the multipliers p appearing in the first of Eqs. (91) are squared to ensure that they remain nonnegative. Further recall that if the constraint is not violated, then p - 0. This condition is satisfied by the second and third equation in Eq. (91) which implies that either k or p is zero for each element. It is readily apparent from Eqs. (86) - (90) that all the costate variables are equal to ~?M1+1 -- x(T). Therefore, these equations were eliminated and all costates that occurred in the remaining equations were replaced with i Ml+l- Also, the state and time boundary conditions in ¢ were solved. The remaining 6M1 + 1 algebraic equations in Eqs. (91) - (94) were solved
A FINITEELEMENTMETHOD
227
using a Newton-Raphson method and a FORTRAN code written on a SUN 3/260. The sparse, linearized equations were solved using subroutine MA28 from the Harwell subroutine library [48]. This subroutine takes advantage of sparsity which leads to great computational savings. Table II shows the convergence rate of XMl+l = x(T), the elapsed computer time for the first 5 iterations, and the percentage of zeroes in the Jacobian (i.e., the sparsity) versus the number of elements. The x(T) column shows that the 32-element case has almost converged on the exact solution. (Recall that the exact solution is x(T) = - 1 . ) Note further that the approximate x(T) is not an upper or lower bound of the exact value, which is common in mixed formulations. The third column of Table II gives the elapsed computer time for five iterations. It is easily seen that there is a modest increase in computer time with an increase in the number of elements. Note that in some cases acceptable convergence criteria were satisfied in fewer than five iterations. This is because the answers obtained from a small number of elements (say 2 or 4) may be interpolated to generate initial guesses for a higher number of elements. Finally, the extremely sparse structure of the Jacobian is demonstrated in the last column. This strongly encourages the use of a smart sparse matrix solver such as MA28. This subroutine leads to quicker solutions and tremendous savings in memory allocation since only the nonzero elements of the J acobian need be stored. T a b l e II. x(T), e l a p s e d c o m p u t e r t i m e a n d p e r c e n t s p a r s i t y of J a c o b i a n v e r s u s t h e n u m b e r of e l e m e n t s M1
M1
x(T)
1 2 4 8 16 32 Exact
-7.5465 -0.4439 -0.8385 -1.0685 -1.0141 -1.0034 - 1
Time (sec) 0.42 0.44 0.66 0.76 1.03 1.52
Sparsity (%) 65.3 80.5 89.6 94.6 97.3 98.6
Results for the state and control are shown in Figures 20 and 21 for 4 and 8 elements and the exact solution. Note from Figure 21 that the 4-element case does not define the constraint boundaries very accurately; however, it is accurate enough to generate guesses for the 8-element case. Thus, in a problem with many constrained and unconstrained arcs, a small number of elements could still be used to generate guesses for a higher number of elements.
228
ROBERT R. BLESS AND DEWEY H. HODGES
One of the constraint multipliers, #1, is shown in Figure 22. Note t h a t #1 = 0 when u < 1, so only half the multiplier history is shown. It can be seen t h a t both the 4- and 8-element cases yield very poor approximations for the multiplier, indicating the need for more elements. However, from the control history plot, the 8-element answer seemed very accurate. This highlights the need to check all results and verify tbnt the multiplier histories seem reasonable. Note also that, due to the s y m m e t r y in this problem, the nonzero segment of #2 is identical to #1.
IV. TREATMENT OF STATE CONSTRAINTS A more difficult constrained problem is one in which the constraints are functions of the state variables only. As will be seen, state constraints produce internal boundary conditions on the states leading to discontinuities in the costates. Choosing the number and type of internal b o u n d a r y conditions can be a nontrivial, trial-and-error task for the analyst. A. G E N E R A L C O M M E N T S Consider now the case where there is a state constraint of the form
S(x) < 0
(96)
imposed on the problem. W i t h the imposed constraint, the optimal trajectory m a y (1) never violate the constraint, (2) touch the state constraint at one or more instances in time (touch points), (3) remain on the curve S = 0 for one or more finite amounts of time (a b o u n d a r y arc), or (4) have some combination of boundary arcs and touch points. The sequence of constrained arcs, unconstrained arcs, and touch points will be referred to as the switching structure. To define the order of the constraint in Eq. (96), successive total time derivatives of S are taken and f(x, u) is substituted for ~ until an expression explicitly dependent on the control u is obtained. If p total-time derivatives are required, then S is called a pth-order state-variable inequality constraint 1. One problem facing the analyst trying to solve a state-constrained problem is finding the switching structure of the constrained arcs. Some general guidelines are given in [49] and [50]. Assuming t h a t the hodograph of the 1One could allow S - S(x, u) so that p = 0. This constraint could then be treated in the framework of this section as opposed to that of Section III. This may provide some advantages in certain situations.
A FINITE ELEMENT METHOD
O
229
~
-1-
-2-
C) t
•
FE(4)
•
FE(8)
-3-
m
"1"(D
-4-
,4..~
09
-5-
-6-
-7-
-8 0
1
2
3
4
Time Figure 20: State history versus time
5
6
7
8
Control History s
6
b
i,--, •
crq
b
6
b
6
~
6
b
o
b
o
b
o
o
~
b
('b
bO
0
Do
0 i--,. ¢-P
o
Or3
=
Or3
CD
c.o
--I
3
j~
oi
--~
o~ -4
Go
-
•
m
x o
mm
m oo
m j~
A FINITE ELEMENT METHOD
231
0.6-
0.5-
FE(4) FE(8) FE(16) Exact
~ " 0.4 0
•
ll/ /
\n \
•
. , m
"r-
~- 0.3 a
~
0.2
0.1
0.0 0
1
2
Time Figure 22" Half of the multiplier history versus time
3
4
232
ROBERT R. BLESS AND DEWEY H. HODGES
problem is strictly convex [51] (also referred to as having a regular Hamiltonian), then it has been shown that the control must be continuous across a junction point (i.e., points where constrained and unconstrained arcs are joined). In this case, for a scalar state constraint, it is proved in [49] that only touch-point solutions are possible for odd-order constraints of order greater than two. Most problems with practical applications, however, involve first- or second-order constraints. In [50], it is shown that first-order constraints will only exhibit boundary-arc solutions, except under three special conditions. One condition is the trivial case, where the constraint is identically satisfied by the unconstrained solution. The second condition arises when there is a corner in the hodograph. The third case is if the pth derivative of the constraint is a global extremum when evaluated at the optimal state and control values. This unusual case is demonstrated in the third example of this section. These three conditions are quite uncommon so the analyst can reasonably assume, at least initially, that there will be a boundary-arc solution for first-order constraints. For a second-order constraint, it is not uncommon that there will be some touch-point cases and some boundary-arc cases depending on the constraint limit. This behavior is demonstrated in the first example. 1. Touch-Point Case For the special case where the trajectory only touches the curve S = 0 at discrete instances of time, then the problem defined in Eqs. (1) - (6) is of the proper form. For example, if there were only one touch point occurring in the solution at some unknown time, then the conditions that S(tl) - 0 and the continuity of all states at time tl, (i.e., x ( t ~ ) - x(t+l) - 0), would be added to the vector ¢ of Eq. (3). Note also that the problem would now have one additional phase due to the extra event at t l. 2. Boundary-Arc Case Consider the situation where the solution has one unconstrained arc between to and t l, followed by a constrained arc between t l and t2, and then another unconstrained arc between t2 and tf (or t3). The augmented cost function defined in Eq. (6) is modified in the second phase to include the state constraint. As is done in [4] and [52], successive time derivatives of the state constraint are taken until the control appears explicitly, say in the pth derivative. Then, to enforce the condition
dPS dtP
=0
(97)
A FINITE ELEMENT METHOD
233
for time between t l and t2, the term
dPS ~7 dtP
(98)
is added to the integrand of Eq. (1) and the augmented cost function of Eq. (6) takes the form
N~I: ( Ja - (~-~-i~ -, H(i)
__ , ~ T ~ +
dPS'~
? ~ . ~ ) dt
(99)
Here, r/is a new multiplier function of time which is nonnegative for tl+ _< t _< t 2 and is identically zero for t + _< t < t I and t + _< t < t 3. By redefining the Hamiltonian in Eq. (5) to be
dPS H - L + )~Tf + rl dtP
(100)
then we see an analogy between the development of the weak formulation for state constraints and that for control constraints as is given in Section III, A. The test and shape functions for 57/and r/will be comparable to those given for 5# and # in Eqs. (77) and (79), respectively. Continuity conditions for the states at t l and t2, namely X(tl ) - x(tt )
-
o
x ( t 2 ) - x(t +)
-
0
(101)
and tangency conditions on the state constraint S
T[x(t+)] = 0 where
T - IS
~d--t
"'"
~,,-1~]~ dtP -1
(102)
(103)
must be satisfied. These internal boundary conditions are enforced by adding these conditions, Eqs. (101) - (102) to the constraint vector ¢ of Eq. (3). B. EXAMPLE: A SECOND-ORDER PROBLEM Consider the double integrator problem in Section 3.11 of [4]. Let x and v define the position and velocity of a particle in rectilinear motion. The governing state equations are
234
ROBERTR. BLESSAND DEWEYH. HODGES
~)
-
v
(104)
-
u
(105)
where the control u is the force, normalized by the mass, imparted to the particle. The problem is to minimize
1/01
u2dt
(106)
However, complications occur when a state inequality constraint is added. Let S - x - L < 0 where L is a constant. The first and second total time derivative of S yield
and S-~) -u Thus, this is a second-order constraint. For certain choices of L, the solution exhibits a touch-point behavior. For other, lower L, the solution rides a boundary arc. First, the unconstrained problem is solved. The state and time boundary conditions are defined to be
x(to)
v(to)- 1 -
X(tl)
(107)
-0
v(tl) + 1 tl-1 Exact and numerical results both show that the state constraint S < 0 is not active for L >_ .25. (The estimated unconstrained value of the limit L is usually a function of the number of elements; this problem is an exception.) For values of L that yield a touch-point solution, the boundary conditions are defined to be
x(t0)
v(to)- 1 x(t+)-L x(t+)-x(tl)
x(t )
v(t2) + I t2- I
-
o
(lo8)
A FINITE ELEMENT METHOD
235
For the constrained boundary-arc case, the boundary conditions are defined to be
x(to)
v(to)- 1 x(t+l ) - L
v(t+l )
x(t+)-x(tl) (t+l ) x(t+ ) - x(t v(t+ ) x(t )
-0
(109)
) )
v(t3) + 1
t3-1 where the tangency conditions are enforced at the unknown entry time, t l, and continuity conditions have been specified at t l and the unknown exit time, t2. In both the unconstrained and constrained cases, the costate boundary conditions are given by Eqs. (15) - (16). Note that the analyst must decide what the switching structure is based on the value of L chosen. Figure 23 shows the displacement of the particle, and Figure 24 shows the control history. In the figures, the exact unconstrained solution is graphed versus a 4-element solution. Also, a touch-point case (with L = 0.2) is graphed versus a finite element run with 4 elements in each phase, denoted with FE(4:4). Finally, a value of L - 0.1 was chosen to demonstrate the exact solution versus the finite element solution (4 elements per phase) for a boundary-arc case. It is seen for this simple example that there is excellent agreement between all results for both the displacement and the control. Of greater interest in this problem is the accuracy with which the finite element method determines the switching structures (i.e., whether the solution is unconstrained, or has a touch-point or boundary-arc solution). For this particular example, the finite element approach identifies the exact unconstrained value of L = 0.25 for the displacement. However, the estimated value of L for which there is a touch-point solution or a boundary-arc solution is a function of the number of elements. The numerical values displayed in Table III below were generated by assuming a touch-point solution. Table III shows the smallest value that L can take and still yield satisfactory results. For instance, when running just 2 elements per phase (the first row of the table), L could be set as low as 0.125 and still yield an approximate answer that satisfies all the discretized necessary conditions. If, however, L - 0.124, then the discretized answer violates the constraint limit before and after (due to the symmetry of the problem) the touch point, which always remains at t = 0.5. As the number of elements is increased, the
236
ROBERT R. BLESS AND DEWEY H. HODGES
0.25
-
0.2-
c" (~
E (:D
0.15-
~J
m
c~
a
l i r a
0.1-
0.05 -
0.0
'
0.0
I
0.2
'
I
'
0.4
I
0.6
Time Figure 23: Displacement versus time
'
I
0.8
~
T
1.0
C~
o
~.~°
o
c-h
o
o0
C~
0~
~.~o
CD
3
0
•
u
b0-
0
0
0
0
0 0
q
•
..s
6"
v
...O --" "" ""
o 0
C ~
I
v
~
O~
o 0
C ~
"1"I rT] m ×
....
Ko Or) ~
0
r" II
m m m m ~ m ~ m ×
. . . . . . . .
". Ko
0
".
0
rII
rII
I i
I
I
,
r" II
0
I
o
..,.. " - " "
.6,
•- - , O
U"
"
i
"~
sS
n
,S
"""~,,
Control
oo
nji o oo
o°
N
•
,
o o"
"",u.
oZ
r"
.n,
"'1, oo
,',II
.............-- "" g ""
° So°
,%
'. n
o I
I
qli
!
~i,
I
,I,
I
- -~
..,..,.- .-- .,,., -,-
-,,..,.....
,
238
ROBERT R. BLESS AND DEWEY H. HODGES
1 allowable values for L increases towards the exact value of ~.
T a b l e I I I . A l l o w a b l e v a l u e s of c o n s t r a i n t limit for t o u c h - p o i n t s o l u t i o n v e r s u s n u m b e r of e l e m e n t s Elements Minimum Value of L for Touch-Point Solution Per Phase 2:2 0.125 (0.124") 4:4 0.146 (0.145") 8:8 0.157 (0.156") 16:16 0.162 (0.161") 32:32 0.165 (0.164") Exact 0.16666666 * denotes the value of L where the constraint is violated Table IV was generated under the assumption that a boundary-arc case existed. This table displays the maximum value of L versus the number of elements where a boundary-arc solution exists. For example, with 2 elements per phase, we could set L = 0.1875. At this value of L, the entry and exit times are almost equal, approximately 0.5 for each. At L = 0.1876, the code produces a solution such that the exit time (t2) is less than the entry time (tl), a nonsense solution. T a b l e I V . A l l o w a b l e v a l u e s of c o n s t r a i n t l i m i t for b o u n d a r y - a r c s o l u t i o n v e r s u s n u m b e r of e l e m e n t s Maximum Value of L for Elements Boundary-arc Solution Per Phase 2:2:2 0.1875 (0.1876") 4:4:4 0.1718 (0.1719") 8:8:8 0.1679 (0.1680") 16:16:16 0.1669 (0.1670") 32:32:32 0.1667 (0.1668") Exact 0.16666666 * denotes the value of L where t2 < t l D
By comparing Tables III and IV, it appears that fewer elements per phase are required to approximate the limiting value between a touchpoint solution and a boundary-arc solution (L = 1/6) when assuming the boundary-arc solution. It is not known if this is true in general.
A FINITE E L E M E N T M E T H O D
239
C. EXAMPLE" B R A C H I S T O C H R O N E P R O B L E M We next consider the brachistochrone problem described in Section 3.11 of [4]. Let x and y define the horizontal and vertical (positive downward) position of a particle, respectively. The governing state equations are
=
2 x / ~ cos u
(110)
-
2 V ~ sin u
(111)
where g is the acceleration due to gravity (= 9.81) and the control u is the angle between the tangent to the path and the horizontal. The problem is to minimize the time it takes for the particle to move from the origin to any point on the line x - t~. A first-order state inequality constraint of the form S = y - x tan 0 - h < 0 is also imposed where 0 and h are constants. The first total time derivative of S yields _ 2v/~ysin(u - O) cosO = 0
(112)
or u = 0 along the constraint boundary. Since this is a first-order constraint, it is reasonable to assume that the constrained solution contains a b o u n d a r y arc. The boundary conditions on the states for this problem may be placed in ~ and the costate boundary conditions are then found from Eq. (15) (16). For the unconstrained case, we have
¢-
{
y(t0)-y0 X(tl)
--
}
--0
(113)
and for the constrained case, we have
x(to) y(to)-
yo
y(t +) - x(t +) tan x(t+)-x(tl) y(t+l ) - y ( t l ) • - x(t ) y(t+ ) - y(t )
0- h -0
(114)
x(t3) where the tangency condition is enforced at the unknown entry time. Notice t h a t the initial value for y has been left as a variable. In the original statement of the problem y0 - 0 for which the Jacobian becomes
240
ROBERTR. BLESSAND DEWEYH. HODGES
singular. Also, the exact solution yields an initial value of the y-costate of - c ~ so that the initial control angle is 90 °. As will be shown, the large and rapid changes in the y-costate value decrease the accuracy of the finite element solution. Figures 25 and 26 show the results for the unconstrained and constrained cases with Y0 - 0.001. The curves are the shooting (essentially exact) results and the discrete points are the finite element values. The constrained case was run with h = 1, t~ = 10, and t a n 0 = 1/2. For the unconstrained case, the finite elements were run with 8 elements. The points are a little off. The reason is clearly seen in Figure 26 (the control history) where it is seen that the final time has been poorly approximated. Also, the initial control angle is not very accurate either. Both of these problems occur because the small initial y0 leads to an extremely rapidly varying costate at the initial time. This is shown in Figure 27 where the unconstrained ycostate history is shown for y0 = 0.001 and Y0 = 0.1. A finer discretization is required to accurately model this costate well. Although the general shape of the state trajectory is found, it is the entry, exit, and final times that are poorly approximated, as is shown in the next two figures. We note that it may be possible to reformulate the problem so that the costate does not tend to - 0 o . With the shooting results taken to be exact, Figures 28 and 29 give relative error of the approximated answers. The relative errors for the entry, exit, and final times, as well as the initial control angle, versus the total number of elements (sum of the elements in each phase) are shown. The slopes of these curves indicate how the relative error varies in proportion to the number of elements. Figure 28 shows the relative errors for Y0 = 0.001. Given that all four curves are nearly parallel, and that the slope continually improves, we can calculate the slopes of the last three segments (the best ones) to be about -1.1, -1.4, and -1.6, respectively. This indicates that the relative error varies proportionally with the inverse of the total number of elements raised to the 1.1, 1.4, or 1.6 powers, and that the error is heading toward second-order accuracy. This accuracy is less than what is typically seen for the finite element method, and is caused by the large gradient in the y-costate value at the initial time. This large gradient in the costate changes as the initial value for y is made larger, as was seen in Figure 27. For y0 = 0.1, the relative errors are plotted in Figure 29 and are seen to be an order of magnitude better. The slopes of the curves are more constant after the initial "bend" in the curve, and the slopes of the last three segments are -1.9, - 2 , and - 2 indicating that the relative error varies approximately inversely with the square of the number of elements. This "second-order" nature of the present finite element method appears to be typical for most problems.
A FINITE ELEMENT METHOD
241
X 0 0 i
1 I
2 I
3 I
4 I
5 I
6 I
7 I
8 I
9 I
10 .I
0
1
,~,,, I
2
• ""o,,,,,,
•
, %,%%,
l
>4
~
1
6
~
7
~
Figure 25: Trajectory versus time
Unconstrained; Shooting Unconstrained; FE(8) h/l=O. 1; Shooting h/l=O. 1; FE(4:4:4)
242
ROBERT R. BLESS A N D D E W E Y H. H O D G E S
90
80
• i
70-
o',,
60-
'\
'\
h/l=0.1;
FE(4:4:4)
• i
50\
m
O L_
O
FE(8)
•
o',
"~ E O
Unconstrained;
• \
-O
Shooting
h/l=O. 1; S h o o t i n g
\
o\
v
Unconstrained;
-
\ i
i
i i
40-
i
i i
_
\
• \ i
30-
\ t
i
\
• ~--e-,
20-
•
O"',
10-
i 0.0
i
i
i
I 0.5
i
,
,
,
I
'
'
1.0
Time (s) Figure 26" Control history versus time
'
'
I 1.5
'
" "
'
I 2.0
A FINITE ELEMENT METHOD
243
-2
-3 y(O) = 0.001; Shooting 0 0 >,,
•
-4
y(O) = 0.001; FE(8) y(O) = O. 1; Shooting
•
-5
y(O) = 0.1; FE(8)
-6
-7
-8
I
0.0
~
w
I
,
I
0.5
~
1
I
~
I
~
~
1.0
Time (s) Figure 27: Costate history versus time
~
'
I
1.5
'
I
2.0
244
ROBERT R. BLESS AND DEWEY H. HODGES
10 0 -
t.__
1 0 "1
0 t._
:: ::::::::i
t._
LU
i
>
,.m
•~
"....
~
m
n-'
1 0 -2
•
.........
• •
"~
"... "\
,, °°°iooo ~,"
"".,,.~.............
Error in Entry T i m e
.................. Error in Exit Time
'.,,~ .......
".,
Error in Final T i m e
10 .3
I
I
I
w
I
I
I
I
i
i
I
10 ~
Number of Elements Figure 28" Error for Yo - 0 . 0 0 1 versus time
i
i
I
I
I
I
10 ~
A FINITE ELEMENT METHOD
245
1 0 o-
,%. '%,,~
] 0 "1 ~ '"'"'"'.. ~.~,
t.__
O
"" """ %'" '~'%'",, ""'......
•~.
•~.
I=.
ILl >
%, %',%
"'....
,~
....
"%.. ~.
%. "%
•~.
1 0 "~
•
'%
".... "~.
, . m
,%
m
~ •~
"'..... "'.... •~ .... "~ "..
rc
%%,
,%. "'... \. "'-.. \
\
"... "'.. "'"'".. \.\. "'".....
u'~''~''~
~. "'% \ "..
1 0 -3 - -
"N.\"'.......
E r r o r in Initial C o n t r o l
",,
"\'N "'"'"'.
•~ "%.
E r r o r in E n t r y T i m e
~
.................. E r r o r in Exit T i m e E r r o r in Final T i m e 10'
,
l
J
,
,
, i I
,
,
I
Number of Elements 29:
Error
f o r yo = 0 . 1 v e r s u s
I
I I , I 10 2
101
Figure
I
"... "\.\. %...... \ "% "\.\ "'.o.., "\
time
246
ROBERT R. BLESS AND DEWEY H. HODGES
D. EXAMPLE: MULTI-PHASE ASCENT TO ORBIT We now turn to a more nonlinear problem to evaluate the usefulness of the finite element method in the presence of state constraints. The same launch vehicle model as used in Section II, E is used, along with a dynamic pressure constraint of the form S ( h , V ) - q - qt~m <_ 0
(115)
This is a first-order state constraint as is seen by taking the first total time derivative and obtaining Op hV 2 -(7=
Oh
2
~- p V ~ /
(116)
By substituting in values for p and h, and simplifying the resulting expression, we obtain the following control-dependent algebraic constraint that must be enforced during the time interval that q = qlim T cos a - D ( a )
p sin 3'
r2
V 2 sin 3' =0 2(8600)
(117)
This constraint will be enforced as one boundary arc occurring before the staging time. Thus, the first phase of the unconstrained problem will be broken into three phases consisting of an unconstrained, constrained, and unconstrained arc. This is followed by the final phase from staging to the final time, for a total of four phases. The boundary condition vector must also be modified to include the tangency conditions and extra internal boundary conditions. The vector is
A FINITE ELEMENT METHOD
247
re(to) - ,no h(to) - 60
V(to)-25 7 ( t 0 ) - 1.5
q(tl+) - q . ~
m(t~+) - m ( t ~ )
) -h(tl) ) - V(tl) ) -7(tl) m(t+~ ) - m(t~) h(t+~ ) - h(t~) V(t+~ ) - V(t~) ~(t~+) - 7(t~) h(t+l V(t+~ ~(t+l
~b
=0
(118)
m ( t +) - r e ( t 3 ) + 29920 h(t +) - h(t3 )
v(t+~ ) - v(t~) .y(t+~ ) - ~(t~) h ( t 4 ) - hy V ( t 4 ) - Vf ~'(t4) t 3 - 195 Some of the results are shown in Figures 30 - 37. The maximum value of the dynamic pressure q is approximately 55.3 kPa. Shooting results and finite element results were obtained for cases where the maximum allowable value for the dynamic pressure is 45 kPa, 35 kPa, and 25 kPa. The dynamic pressure histories are shown in Figure 30. However, for clarity in the figures, results will only be presented for the 25 kPa case. In all of the figures, the unconstrained shooting results are shown, for comparison purposes. The finite element results are for 4 elements in each of the first three phases and 8 elements in the last phase (the one after staging) and are denoted with "FE(4:4:4:8)". Figure 31 shows the normalized altitude profile versus the normalized velocity. As expected, the vehicle climbs higher, more quickly, as the constraint becomes stricter. Note also that even with the coarse discretization presented, the finite element results lie essentially on top of the exact curves. Figure 32 shows the flight-path angle 7 versus time. The vehicle is seen to roll over more slowly with the strict q limit in order to climb out of the atmosphere more quickly. Recall from Figure 13 that the finite element results show some error for the unconstrained results. However, with the addition of finite elements and the less rapidly varying flight-path angle history, the answers are much more accurate for the constrained case.
248
ROBERT R. BLESS AND DEWEY H. HODGES
Figures 3 3 - 36 give four costate histories for the vehicle. Due to the internal boundary condition on q, a jump occurs in the altitude costate and the velocity costate. Note also that the finite elements have solved for this jump and the initial costate values accurately. These values are, of course, critical to the shooting code. The angle of attack history is shown in Figure 37. The control experiences large variations due to the changing constraint limits. Once again the finite element results are quite accurate, although less so at the initial time where the angle of attack changes rapidly. As mentioned at the beginning of this section, an interesting touchpoint solution exists for this problem. The control is positive along the constrained arc for certain values of the dynamic pressure limit, and is negative along the constrained arc for other values. At a particular value of qtim, the control will be zero, and a touch-point solution occurs. A first-order touch-point solution is possible because, from Eq. (117)
0b
0v
0--~ = 0 a = 0
at
a - 0
(119)
Figure 38 shows the length of the constrained arc and final mass loss as a function of dynamic pressure limit. It is interesting to note that the length of the constrained arc is extremely small (less than 0.5 seconds) until after the limit is set below the touch-point limit, which occurs just under the 35 kPa arc, at which point the length of the constrained arc grows rapidly as the limit is decreased. Also of interest is the fact that the loss of final mass becomes increasingly steeper with a stricter dynamic pressure limit. The combination of rapidly obtained finite element solutions and accurate shooting solutions allowed for the discovery of this anomalous first-order state-constraint behavior.
C~
=
CD C~
t~
(1) oo 0° o
O~
m v
3
--I
0 0
0 0
0
~ o
o
o
0
0
L
I
0
I
I
I
I~ 0
I I I I t
0 ~'.
o
~.
oo
o 0
I
II 0"I ~-D
II 07 ~-D
o1 ~-D
II
x .o
0 0 =. ~
0 0 ~
=.
OOO0007"I
x .Q
x .0
0 0 =. ~
"'"
66
"~ ..
~
O1 ~. "D
II
x c~
•
I
.~ 0
""--~)
..~.. ~ ~ -
--..:........ " " ..--. .-...,,,,...
Z.-.--._-.ZT-.
................................
I
..... •.
...................T..." .... .....
, ' , 3 3 3 3
I I I i
I
0.~ 0
I
Dynamic Pressure (kPa) I
01 0
. ..,...,.,.... ~ "
-'.--.~..~
I
I
0"~ 0
250
ROBERT R. BLESS AND DEWEY H. HODGES
1.0~
A
A
_
..~ 0.8
_~~o,6 N 0.4
Z
0.0
:4:4:8) max q = 25kPa; Shooting oting
/,,,,"
0.2
•
0.0
'
I
0.2
'
I
0.4
'
I
0.6
'
I
0.8
NormalizedVelocity
Figure 31" Normalized altitude versus normalized velocity
~
I
1.0
A FINITE ELEMENT METHOD
251
80
70
¢D ,<: _
5o
i "
40
",,,,, \
•
max q = 25kPa; FE(4:4:4:8)
-
max q = 25kPa; Shooting
.&, r--
.~
2o-
IJ _
-
100-
-10
'
0
I
100
'
I
'
200
I
300
Time (s) Figure
32" F l i g h t - p a t h
angle versus
time
'
I
400
'
I
500
252
ROBERT R. BLESS AND DEWEY H. HODGES
o.o-0.1
"'"" ..... .
-0.2 -0.3 -O.4
8
0
-o.5
~
-0.6 -0.7 -0.8 -0.9 -1.0
I 0
100
200
300
Time (s) Figure 33" Mass costate versus time
400
500
A FINITE ELEMENT METHOD
253
0.0 zooOO° ° ° ° °
I
-0.05
oo
o
(D
-0.1
0
0
(D "0
m
m
< m
-0.15
max q = 25kPa; FE(4:4:4:8)
-0.2
max q = 25kPa; Shooting q unconstrained; Shooting
-0.25
-0.3
I
0
'
I
100
'
I
200
'
I
300
Time (s) Figure 34: Altitude costate versus time
'
I
I
400
5OO
254
ROBERT R. BLESS AND DEWEY H. HODGES
00-] -0.1
-0.2 .b..,0
-0.3
0
0
>.,
i
m
-0.4
m
0
0 (D
>
-0.5
max q = 25kPa; FE(4:4:4:8) max q = 25kPa; Shooting q unconstrained; Shooting
-0.6
-0.7 I
-0.8
'
0
I
100
'
I
200
'
I
300
Time (s) Figure 35: Velocity costate versus time
'
I
i
400
5OO
A FINITE ELEMENT METHOD
255
"10 ~
0.2
max q - 25kPa; FE(4:4:4:8)
0.04
max q = 25kPa; Shooting q unconstrained; Shooting
i i
\
-0.2 -
i
0O
0
0
i i
\
-0.4-
\
i
i
(D
133
/
/
#
-0.6 /
<
t-
/
-0.8-
II
/
\, \
u
',,,,
-1.0-
c-"
m
m
O3 m
U..
/
",,
-1.2-1.4-
~,,
ooo1°
,,/
•
•
-1.6'
0
I
100
'
I
2OO
'
I
300
Time (s) Figure 36" Flight-path-angle costate versus time
'
U
400
'
I
500
256
ROBERT R. BLESS AND DEWEY H. HODGES
3025-
max q - 25kPa; FE(4:4:4:8) m a x q - 2 5 k P a ; Shooting
20
(D "O
t ! t
i
,
15-
;i
q unconstrained; \
",\ i
i
i
\ i
v
rD
Shooting
10-
i
i
i II i i
< O (D
%%%
O
_
m~=..=
¢-
<
-5-10-15-20
' 0
I 100
I 200
'
I 300
Time (s) Figure 37" Angle of attack versus time
'
I 400
'
I 500
0
0
r~ r~
0
r/l
.°
V
m .
3
m .
r-
"u
0
m l
3
01
O
0
O
I
t
O Q
i
I
~
0
t 0 0 Q
I
'
i
I
T 01 0 Q
(/)
E
..,..
"1"1 --,.,,.
0
0
,.,..
¢.. (-)
:3:3
i
Mass Loss (kg)
::)>
p-..q
0 :::I U)
0
:3"
i
I
0 Q O
i
Length of Constrained Arc (s) I
,
0
01
I
258
ROBERT R. BLESS AND DEWEY H. HODGES
V. SUMMARY This chapter was concerned with a method for the solution of optimal control problems. The method, called the weak principle for optimal control, is based on time-domain finite elements. The goals of this work, as stated in Section I, B, were to develop an efficient and accurate algorithm that would (1) rapidly obtain an approximate solution for possible real-time implementation; (2) identify the optimal switching structure for constrained problems; and (3) produce initial guesses for a shooting method. With respect to these three goals, it has been observed that: 1. It is helpful to obtain a converged answer ments and then add elements using guesses converged solution. As was demonstrated II, E (see Table I), this proved to be a fast ing accurate answers which could then be
for a small number of eleobtained from a previously in the example in Section and reliable way of obtainused in a shooting code.
2. After the dynamics of the system have been properly modeled with some minimum number of elements, then the finite element method possesses second-order accuracy. It was shown in Section IV, C that there is a degradation in accuracy if the variables change more rapidly than the finite element discretization is able to model. 3. As was shown in Tables III and IV of Section IV, B, this formulation can be useful in determining the switching structure for an optimal trajectory. 4. The speed of solution is directly related to the speed of solving a sparse linear system. The sparsity of the Jacobian was highlighted in Table II of Section III, C. With the use of an intelligent sparse solver, this finite element method can be very efficient in terms of run-time and memory requirements. There is now a patent pending on a real-time guidance algorithm based on the present methodology
[53].
Some of the characteristics and features of the weak principle are summarized below. 1. The necessary conditions of optimality are discretized, and all strong boundary conditions are transformed into weak ones. The weak principle finds candidate extremal solutions, i.e., ones that satisfy all the discretized necessary conditions. 2. Because all strong boundary conditions are cast in the form of natural or weak boundary conditions, then the same test and shape functions may be chosen for every optimal control problem.
A FINITE ELEMENT METHOD
259
3. The choice of shape functions allowed by the weak principle permits all integration to be done by inspection, regardless of the degree of nonlinearity in the problem. Thus, no errors are introduced as the result of numerical quadrature. 4. One advantage of the integration being done by inspection is that algebraic equations may be derived prior to specifying the problem to be solved. Because the form of the algebraic equations is known, it was possible to create a general algorithm, VTOTS, for the solution of optimal control problems. 5. The derivations in this chapter lay the groundwork for higher-order test and shape functions to be used. Although this introduces the need for element quadrature, greater accuracy is obtained with fewer elements [54].
VI. A C K N O W L E D G E M E N T S Technical discussions with Drs. Marco Borri, Anthony J. Calise, J. Eric Corban, Daniel D. Moerder, David A. Peters, J. V. R. Prasad, and Hans Seywald are gratefully acknowledged. The authors also thank Michael S. Warner for editorial comments and a careful reading of the manuscript. This work was partially supported by grants from NASA Langley Research Center, the technical monitor of which was Daniel D. Moerder.
References [1] P. E. Gill, W. Murray, and M. H. Wright, Practical Optimization, Academic Press, New York (1981). [2] J. Betts and W. Huffman, "The Application of Sparse Nonlinear Programming to Trajectory Optimization," Proceedings of the AIAA Guidance, Navigation, and Control Conference, Portland, Oregon, AIAA Paper 90-3448, pp. 1198-1218 (1990). [3] J. E. Rader and D. J. Hull, "Computation of Optimal Aircraft Trajectories Using Parameter Optimization Methods," Journal of Aircraft 12(11), pp. 864-866 (1975). [4] A. E. Bryson, Jr. and Y.-C. Ho, Applied Optimal Control, Hemisphere Publishing Corporation, New York (1975).
260
ROBERT R. BLESS AND DEWEY H. HODGES
[5] E. B. Lee and L. Markus, Foundations of Optimal Control Theory, John Wiley & Sons, Inc., New York (1967).
[6]
L.D. Berkovitz, Optimal Control Theory, Springer-Verlag, New York (1974).
[7]
C. R. Hargraves and S. W. Paris, "Direct Trajectory Optimization Using Nonlinear Programming and Collocation," Journal of Guidance 10(4), pp. 338-342 (1987).
[8] C. Hargraves, F. Johnson, S. Paris, and I. Rettie, "Numerical Com-
putation of Optimal Atmospheric Trajectories," Journal of Guidance and Control 4(4), pp. 406-414 (1981).
[9] J. Vlassenbroeck, "A Chebychev Polynomial Method for Optimal Control with State Constraints," Automatica 24(4), pp. 499-506 (1988). [10] C. T. Kelley and E. W. Sachs, "Quasi-Newton Methods and Unconstrained Optimal Control Problems," SIAM Journal on Control and Optimization 25(6), pp. 1503-1516 (1987). [11] W. A. Gruver and E. Sachs, Algorithmic Methods in Optimal Control, Pitman Publishing, Boston, Massachusetts (1981).
[12] P. E. Gill, W. Murray, M. A. Saunders, and M. H. Wright, "User's
Guide for NPSOL: A FORTRAN Package for Nonlinear Programming," Technical Report SOL 86-2, Department of Operations Research, Stanford University (1986).
[13] S. W. Paris and C. R. Hargraves, "Optimal Trajectories by Implicit Simulation: Volume I-Formulation Manual," AFWAL-TR-883057 (1988).
[14] G. L. Brauer, et al., "Final Report: Program To Optimize Simulated Trajectories (POST): Volume I-Formulation Manual," NASA CR132689 (1975). [15] H. J. Oberle and W. Grimm, "BNDSCO: A Program for the Numerical Solution of Optimal Control Problems," English Translation of DLRMitt. 85-05, Oberpfaffenhofen, Germany. [16] R. R. Bless, et al., "Variational Trajectory Optimization Tool Set; Technical Description and User's Manual," NASA TM 4442 (1993). [17] S. M. Roberts and J. S. Shipman, "Multipoint Solution of Two-Point Boundary-Value Problems," Journal of Optimization Theory and Applications 7(4), pp. 301-318 (1971).
A FINITEELEMENTMETHOD
261
[18] 5. M. Roberts and J. S. Shipman, Two-Point Boundary Value Problems: Shooting Methods, American Elsevier Publishing Company, Inc., New York, New York (1972). [19] P. K. A. Menon and L. L. Lehman, "A Parallel Quasi-Linearization Algorithm for Air Vehicle Trajectory Optimization," Journal of Guidance 9(1), pp. 119-121 (1986). [20] M. B. Subrahmanyam, "A Computational Method for the Solution of Time-Optimal Control Problems by Newton's Method," International Journal of Control 44(5), pp. 1233-1243 (1986). [21] D. E. Kirk, Optimal Control Theory: An Introduction, Prentice-Hall Inc., Englewood Cliffs, New Jersey (1970). [22] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, Cambridge (1986). [23] H. G. Visser, H. J. Kelley, and E. M. Cliff, "Energy Management of Three-Dimensional Minimum-Time Intercept," Journal of Guidance, Control and Dynamics 10(6), pp. 574-580 (1987). [24] H. Seywald, "Optimal Control Problems with Switching Points," NASA c a 4393 (1991). [25] I. Fried, "Finite Element Analysis of Time Dependent Phenomena," AIAA Journal 7, pp. 1170-1173 (1969). [26] J. H. Argyris and D. W. Scharpf, "Finite Elements in Time and Space," Aeronautical Journal of the Royal Society 73, pp. 1041-1044 (1969). [27] C. D. Bailey, "Application of Hamilton's Law of Varying Action," AIAA Journal 13(9), pp. 1154-1157 (1975). [28] T. E. Simkins, "Unconstrained Variational Statements for Initial and Boundary Value Problems," AIAA Journal 16, pp. 559-563 (1978). [29] D. L. Hitzl and D. A. Levinson, "Application of Hamilton's Law of Varying Action to the Restricted Three-Body Problem," Celestial Mechanics 22, pp. 255-266 (1980). [30] M. Baruch and R. Rift, "Hamilton's Principle, Hamilton's Law, 6n Correct Formulations," AIAA Journal 20, pp. 687-692 (1982). [31] M. Borri, et al., "Dynamic Response of Mechanical Systems by a Weak Hamiltonian Formulation," Computers and Structures, 20(1-3), pp. 495-508 (1985).
262
ROBERTR. BLESSANDDEWEYH. HODGES
[32] D. A. Peters and A. Izadpanah, "hp-Version Finite Elements for the Space-Time Domain," Computational Mechanics 3, pp. 73-88 (1988). [33] W. E. Bosarge. Jr. and O. G. Johnson, "Error Bounds of High Order Accuracy for the State Regulator Problem Via Piecewise Polynomial Approximations," SIAM Journal on Control 9(1), pp. 15-28 (1971). [34] P. M. Prenter, Splines and Variational Methods, John Wiley & Sons, Inc., New York, New York (1975). [35] L. J. Segerlind, "Weighted Residual Solutions in the Time-Domain," International Journal for Numerical Methods in Engineering 28, pp. 679-685 (1989). [36] W. N. Patten, "Near Optimal Feedback Control for Nonlinear Aerodynamic Systems with an Application to the High-Angle-of-Attack Wing Rock Problem," AIAA Paper 88-4052-CP (1988). [37] F. G. Mathis and G. W. Reddien, "Ritz-Trefftz Approximations in Optimal Control," SIAM Journal on Control and Optimization 17(2), pp. 307-310 (1979). [38] B. Asselmeyer, "A Ritz-Type Optimization Method for Optimal Control Problems and Its Application to Hierarchical Final-Value Controllers," Control and Dynamic Systems: Advances in Theory and Application 23, pp. 295-318 (1986). [39] H. Oz and A. Raffle, "Inverse Response Problem (Control) of Dynamic Systems Via Hamilton's Law," Computer Methods in Applied Mechanics and Engineering 62, pp. 17-26 (1987). [40] D. H. Hodges and R. R. Bless, "Weak Hamiltonian Finite Element Method for Optimal Control Problems," Journal of Guidance, Control, and Dynamics 14(1), pp. 148-156 (1991). [41] D. H. Hodges, R. R. Bless, A. J. Calise, and M. Leung, "Finite Element Method for Optimal Guidance of an Advanced Launch Vehicle," Journal of Guidance, Control, and Dynamics 15(3), pp. 664-671 (1992). [42] R. R. Bless and D. H. Hodges, "Finite Element Solution of Optimal Control Problems with State-Control Inequality Constraints," Journal of Guidance, Control, and Dynamics 15(4), pp. 1029-1032 (1992). [43] M. Borri, "Helicopter Rotor Dynamics by Finite Element Time Approximation," Comp. and Mathematics with Applications 12A(1), pp. 149-160 (1986).
A FINITE ELEMENT METHOD
263
[44] M. Borri, F. Mello, and S. N. Atluri, "Primal and Mixed Forms of Hamilton's Principle for Constrained Rigid Body Systems: Numerical Studies," Computational Mechanics 7, pp. 205-220 (1991). [45] I. M. Gelfand and S. V. Fomin, Calculus of Variations, Prentice-Hall, Inc., Englewood Cliffs, New Jersey (1963).
[46]
E. B. Becker, G. F. Carey, and J. T. Oden, Finite Elements: An Introduction, Volume 1, Prentice-Hall, Inc., Englewood Cliffs, New Jersey (1981).
[47]
F. A. Valentine, "The Problem of Lagrange with Differential Inequalities as Added Side Conditions," Contributions to the Calculus of Variations, Chicago, IL, Chicago University Press, pp. 407-448 (1937).
[48]
I. S. Duff, Harwell Subroutine Library, Chapter M, Computer Science and Systems Division, Harwell Laboratory, Oxfordshire, England (1988).
[49]
D. H. Jacobson, M. M. Lele, and J. L. Speyer, "New Necessary Conditions of Optimality for Control Problems with State-Variable Inequality Constraints," Journal of Mathematical Analysis and Applications 35, pp. 255-284 (1971).
[50]
H. Seywald, "On the Existence of Touch Points for First-order State Inequality Constraints," Proceedings of the AIAA Guidance, Navigation, and Control Conference, Monterey, California, AIAA Paper 93-3743CP, pp. 372-376 (1993).
[51]
E. M. Cliff, H. Seywald, and R. R. Bless, "Hodograph Analysis in Aircraft Trajectory Optimization," Proceedings of the AIAA Guidance, Navigation, and Control Conference, Monterey, California, AIAA Paper 93-3742-CP, pp. 363-371 (1993).
[52]
A. E. Bryson, Jr., W. F. Denham, and S. E. Dreyfus, "Optimal Programming Problems with Inequality Constraints, I: Necessary Conditions for Extremal Solutions," AIAA Journal 1(11), pp. 2544-2550 (1963).
[53]
G. T. Haskins, M. G. Johnson, and D. H. Hodges, "Real Time Missile Guidance System," United States Patent Pending (1993).
[54] D. H. Hodges and L. J. Hou, "Shape Functions for Mixed p-Version Finite Elements in the Time Domain," J. Sound and Vibration 145(2), pp. 169-178 (1991).
Robert R. Bless Lockheed Engineering and Science Company Hampton, Virginia, 23666
Dewey H. Hodges School of Aerospace Engineering Georgia Institute of Technology Atlanta, Georgia 30332-0150
I. I N T R O D U C T I O N This chapter deals with a new type of numerical method based on finite elements in time for solving optimal control problems. The optimal control problem of interest in this chapter may be described as follows. Consider a system that is completely defined by a finite number of states, i.e., quantitles that describe the current status of the system. The status of the system is determined by a set of first-order ordinary differential equations referred to as state equations. The states are influenced by a finite number of control variables. The optimization problem is to choose the control variables to satisfy the given state equations, boundary conditions, and any constraints imposed on the states and/or controls, while minimizing or maximizing a given performance index, or cost functional. Use of the calculus of variations results in a multi-point boundary-value problem. Due to nonlinearities occurring in the state equations and the complex nature of constraints, not many analytical solutions to these types of optimal control problems have been found. Thus, the analyst must resort to numerical methods in order to solve the problem. CONTROL AND DYNAMIC SYSTEMS, VOL. 72 Copyright © 1995 by Academic Press, Inc. All rights of reproduction in any form reserved.
183
184
ROBERT R. BLESS AND DEWEY H. HODGES
A. METHODS OF SOLUTION Methods available for the solution of optimal control problems generally fall into two distinct categories: direct and indirect. Among the direct methods are those that transcribe the infinite-dimensional continuous problem to a finite-dimensional nonlinear programming (NLP) problem by some parameterization of the control variables and, possibly, the state variables. A discussion of NLP problems and the solution of these problems is given in [1]. Continuing advances in NLP algorithms and related software have made these the methods of choice in many applications [2, 3]. Indirect techniques, on the other hand, seek to minimize the performance index indirectly by satisfying the first-order necessary conditions for optimality as established from the calculus of variations (see, for example, [4, 5, 6]). Below is given a brief description of some of the indirect and direct methods now being used to solve optimal control problems. Since finite element methods may be derived either directly or indirectly, a separate discussion of them is given. 1. Direct Methods The direct approach to the solution of optimal control problems requires parameterization of the control and/or state time histories. The choice of parameterization schemes is not unique, and success of the direct methods has been achieved using schemes such as cubic polynomials [7] and Chebychev polynomials [8, 9]. Once the parameterization scheme is chosen, a parameter optimization algorithm is then used to improve the initial guess of the free parameters. These algorithms are in common use today and include quasi-Newton methods [10], sequential-quadratic programming (SQP) methods [1], and gradient methods [4]. Gradient methods were developed to surmount the "initial guess" difficulty associated with other methods such as NewtonRaphson. They are characterized by iterative algorithms for improving estimates of the control parameters in order to come closer to satisfying the stationarity conditions. First-order gradient methods usually show rapid improvements when sufficiently far from the optimal solution. However, the rate of convergence drastically decreases in the neighborhood of the solution. Second-order gradient methods have excellent convergence characteristics near the optimal solution, similar to a Newton-Raphson method. Conjugate gradient methods are very powerful because they combine the first-order and second-order gradient methods. A thorough description of the gradient method and many other algorithmic methods in optimal con-
A FINITE ELEMENT METHOD
185
trol may be found in [11]. Below are mentioned three of the widely used codes for the direct solution of optimal control problems. The first one is a FORTRAN code called NPSOL [12]. NPSOL is used to solve NLP problems and also serves as the solver for another code, OTIS, for solving optimal control problems. OTIS, or Optimal Trajectories by Implicit Simulation [13], is a three-degree-offreedom (point-mass) simulation program for multiple vehicles. The user can simulate a wide variety of vehicles such as aircraft, missiles, re-entry vehicles, and hypervelocity vehicles. A third code is POST, or Program to Optimize Simulated Trajectories [14]. POST provides the capability to target and optimize point-mass trajectories for a powered or unpowered vehicle operating near a rotating oblate planet. POST offers the solution to a wide range of flight problems including aircraft performance, orbital maneuvers, and injection into orbit. 2. Indirect Methods Indirect methods are based on finding the solution of a boundary-value problem which results from the first-order necessary conditions of optimal control. For many practical optimization problems the boundary-value problems are quite difficult; but, here again, modern numerical algorithms and associated software have enlarged the class of solvable problems significantly [15, 16]. Two methods for solving nonlinear multi-point boundary-value problems are shooting methods and quasilinearization methods. Shooting methods [17, 18] are frequently used and can be described as follows: The initial conditions and the differential equations are satisfied at each stage of the process while the final conditions are sacrificed somewhat. A nominal solution is generated by guessing the missing initial conditions and forward integrating the differential equations. The intent is to reduce the error in the final conditions at each iteration. Quasilinearization techniques [4, 19] involve choosing nominal functions for the states and costates that satisfy as many of the boundary conditions as possible. The control is then found by using the optimality conditions. The system equations and costate equations are then linearized about the nominal values, and a succession of nonhomogeneous, linear two-point boundary value problems are solved to modify the solution until the desired accuracy is obtained. Other indirect techniques include Newton methods [2.0], steepest descent methods [21], finite difference [22], and the method of adjoints and continuation methods
[lS]o
The virtue of the indirect methods is the high precision they offer and their rapid convergence in the immediate neighborhood of the optimal so-
186
ROBERTR. BLESSANDDEWEYH. HODGES
lution. These are important features, for example, when one is conducting sensitivity studies [23]. In addition, optimal trajectories are often composed of sequences of arcs along which various state or control constraints are alternately active and inactive, referred to as the switching structure [24]. Indirect methods require that the analyst find the correct switching structure, generally by trial-and-error, in order to obtain a converged solution. Even very short arcs must be anticipated. With NLP-based approaches, switching structures do not have to be known; however, the analyst can have difficulties in correctly identifying some optimal switching structures, particularly those with very short arcs that are missed by the parameterization scheme chosen with the direct methods. 3. Finite Element Methods Finite element methods originated in the field of structural mechanics in the 1950s, and the number and types of applications have since been greatly expanded. They are formulated by considering a set of spatial subdomains, called finite elements, which make up the domain of interest. The behavior of the set of finite elements, with certain inter-element continuity conditions enforced, is intended to approximate the behavior of the domain. The behavior of each element is analyzed based on a restricted set of modes, called shape functions. Usually the shape functions are low-order polynomials and are chosen so that the spatial degrees of freedom are either displacements or rotations at specified points, called nodes, or are a combination of nodal degrees of freedom and coefficients of polynomials that vanish at the nodes. We consider finite element methods separately from the direct and indirect methods, because finite element methods can be developed either from energy principles (directly) or from application of Galerkin's method to the governing partial differential equations (indirectly). In some instances the resulting spatially discretized equations are exactly the same, regardless of which path is taken in the derivation. Finite element methods have also been developed for the time domain. The earliest works seem to be those of Fried [25] and Argyris and Scharpf [26], although some of the later work, including the work presented in this chapter, has its roots in the mid-1970s with the development of direct time-domain methods for dynamics based on Hamilton's principle. Hamilton's principle had traditionally been used in analytical mechanics only as a method of obtaining the equations of motion for dynamical systems. Bailey [27] followed by several others [28, 29, 30] obtained direct solutions to dynamics problems using a form of Hamilton's principle known as the law of varying action, thus opening the door for its use in computational mechanics. Recasting of Hamilton's law as a weak form, referred to as Hamilton's
A FINITE ELEMENT METHOD
187
weak principle or HWP [31], was shown to provide computational characteristics that are superior to those of the form used by Bailey. Bailey's form was later shown to occasionally exhibit convergence problems [32]. The accuracy of the time-marching procedures derived in [31, 32] is competitive with standard ordinary differential equation solvers and provides a powerful alternative to numerical solution of ordinary differential equations in the time domain. In what appears to be one of the first papers using finite elements to solve optimal control problems [33] the authors applied a modified RitzTrefftz direct method to state-regulator control problems. Finite element methods based on Ritz and Galerkin methods [34, 35] as well as the method of collocation (see, for example, [36]) have also been used to solve optimal control problems. Other finite element methods based on the method of Ritz can be found in [37, 38], and one based on collocation is found in [7]. A direct method for solution of inverse response problems is presented in [39]. The method is set up to "march" in time using a direct representation of the dynamics via Hamilton's law of varying action instead of using the standard state-space ODE representation of the system dynamics. B. PRESENT APPROACH The methods described above, although very accurate and useful, suffer from some computational challenges. For example, the approximating functions must satisfy all strong boundary conditions, which means that one set of functions may not suffice for all types of problems or that certain equations and unknowns need to be removed depending on the particular problem. This turns the development of a general-purpose algorithm into an even more complex task. Also element quadrature in the previous methods must be done by numerical means, which greatly increases the computational effort. This chapter describes in detail a computationally efficient, generalpurpose time-domain finite element approach, taken from [40, 41, 42], for solving optimal control problems based on a weak formulation of the firstorder necessary conditions. The above problems are avoided by reformulating the variational problem so that all boundary conditions appear as natural boundary conditions. This way, the simplest possible approximating functions are allowed for all cases so that all element quadrature may be done by inspection, and the Jacobian matrix involved in the solution procedure is very sparse. The equations can be set up for solution prior to specifying a particular problem, making the algorithm ideal for generalpurpose computer software. The background of the present method of treating the time domain
188
ROBERT R. BLESS AND DEWEY H. HODGES
lies in extensions of HWP [31]. HWP was 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 [43]. These are complex two-point boundary value problems, and the utility of HWP for these problems strongly suggests that it could be used in optimal control problems. Finally, an important step was the development of a mixed form of HWP [44], in which the generalized coordinates and momenta appear as independent unknowns. The mixed form provides additional computational advantages over the original displacement form [31], in which only generalized coordinates appear. These advantages include the possibility of using simpler shape functions and its ability to yield an unconditionally stable algorithm for the linear oscillator with exact element quadrature. Although the early developments with Hamilton's law and HWP are direct methods, the mixed form of HWP is really just a weak form of the governing equations, including boundary conditions. Therefore, a finite element method based on it is an indirect method. The finite element method to be described herein is also an indirect method. Of the indirect methods listed above, shooting methods yield numerically exact solutions that are as accurate as the integrator used for the differential equations. Shooting methods are, unfortunately, very sensitive to initial guesses; hence, it may be difficult to obtain a converged solution without good initial guesses for the costates. In addition, shooting methods are not well-suited for real-time applications due to the computationally intensive operations associated with repetitive numerical integrations needed to calculate the J acobian. Therefore, the purpose of this chapter is to derive and evaluate another type of indirect method which will 1. rapidly obtain an approximate solution for possible real-time implementation; 2. identify the optimal switching structure for constrained problems; and 3. produce initial guesses for a shooting method. This chapter describes a finite element method based on a weak formulation of the first-order necessary conditions of optimality to achieve these three goals. The remainder of this section will mathematically define the optimal control problem. Then, Section II develops the weak principle for optimal control theory for unconstrained problems. Two examples are given to illustrate the use of the weak principle. In one example, the states and controls are continuous. The second example involves a problem with discontinuities in the states and system equations. Further extensions of the weak
A FINITE ELEMENT METHOD
189
principle are given in Sections III and IV with the inclusion of control and state inequality constraints. Several example problems are given in these two sections to highlight the behavior of the finite element method. Some features of the weak formulation and conclusions are discussed in Section V. C. O P T I M A L C O N T R O L PROBLEM D E F I N I T I O N Consider a system defined over a time interval from the initial time to to the final time t I by a set of n states, x, and a set of m controls, u. The states of the system are governed by a set of first-order differential equations referred to as state equations. During the time interval from to to t f, there may be discontinuities in the states, or discontinuities in the state equations at interior points (times between to and t l ) , or control and/or state constraints to satisfy. These interior points, along with the initial and final points, will be referred to as events, and the time interval between events will be referred to as phases. The time of event i will be denoted with ti and times just before and after ti will be denoted with t~and t +, respectively. All event times will be considered unknown except that of the zeroth event, coinciding with the initial time to. Elements of a performance index may be denoted with an integrand L(i)[x(t), u(t)] over each phase i and a discrete function ¢ of the states and/or times at any of the events. A general class of such problems with N phases involves choosing u(t) to minimize J
--
¢[x(t+o),x(tl),X(t+l),...,X(tN);tl,t2,...,tN]
+
~
L(i)[x(t),u(t)] dt
iX'~
(1)
¢--1
subject to the state equation constraints
Jc - f(i)Ix(t), u(t)];
(t+_l < t < t~-; i - 1 , . . . , N)
(2)
with boundary conditions on the states and times specified as
),...
,tN] - o
(3)
We introduce Lagrange multiplier functions )~(t), referred to as costates, and discrete Lagrange multipliers v. Then, for convenience, (I) and H are defined as
¢
-
¢+~¢
H(i)
_
L(i) + ATf(i)
(i -- 1 , . . . , N )
(4) (5)
190
ROBERT R. BLESS A N D D E W E Y H. H O D G E S
The augmented cost function obtained by adjoining Eqs. (2) and (3) to Eq. (1) using the costates A(t) and multipliers u, respectively, is
Ja = ~ + ~~1=
-1 H(~) _
)~T~
) dt
(6)
The derivation of the first-order necessary conditions for optimality begins by expanding the increment of the augmented cost function of Eq. (6) (i.e., the difference between the cost function evaluated at a solution neighboring to the optimal and the cost function evaluated at the optimal solution) into a Taylor series about the optimal solution. The first-order term of the Taylor series is referred to as the first variation. Provided that the variations to the optimal solution are small enough so that the sign of the increment of the cost function is determined by the first variation, then it is well-known [45] that the optimal solution furnishes a stationary point to Ja. The first variation is ;
+ 6uT
Ou
dt + du T¢
(7)
N-1
+
Z i=1
+
-~i + L(i) (t:~ ) -
(t +) dti + -~N + L(N) (tN) dtg
[ ]
Zi =~1 Ox(t;)°~ dx(t;)+ Zi--0 Ox(t+)
where 6x(t) signifies the perturbation of x when t is considered fixed and dx(t +) and dx(tT~) signify the perturbation of x due to perturbations of ti, as well as due to changes in x at t + and t~-, respectively. These are related as
dx(t +) dx(tT ) -
6x(t +) + 2(t + )dti 6x(t 7) + ~(t 7 )dti
(8) (9)
Note that since x can be discontinuous at ti, so can the perturbation of x. After integrating the 6~TA term by parts and using Eqs. (8) - (9) to eliminate 6x(t +) and 6x(t:~) in terms of dx(t +), dx(t:~ ), and dti, we obtain
i=1
,-
A FINITE ELEMENT METHOD
5uT
OU
191
(10)
dt + duT¢
N-1
+
E i=1
-~i + H(i)(t[) -
(t +) dti + -~N + H(N)(tN) N-1
0e
Ox(t+~------~+
dtN
~ (t+)] ex(t+ )
The first-order necessary conditions of optimal control require that 5Ja - 0
(11)
for all independently chosen perturbations dtl, dt2, ..., dtN, du, dx(t+), dX(tl) , dx(t+), ..., dX(tN) , 5x(t), 5,~(t), and 5u(t). Hence, the discrete boundary condition terms are c0(I) H(i+I) Ot---~+ H(O (t;) (t +)
-
0
O0 Ot---N +
H(N)(tN)
-
0
(13)
~[x(t0+),x(t 1), . . . , x(tN); tl, . . . , tN]
--
0
(14)
-
0
(i- 1,2,...,N)
(15)
+ ~(t+~ ) -
0
(i - 0, 1 , . . . , N -
1)(16)
0O o~(t;)
OO
o~(t+~ )
-
;~(t~-)
(i- 1,2,...,N-
1)(12)
where Eq. (12), for instance, is obtained by setting all perturbations equal zero except for dti, which in turn forces the coefficient of dti to be zero. The Euler-Lagrange equations that come from Eq. (10) for each phase i are :~_f(0 (OH(i))
~+
=
0
(17)
-
o
(18)
=
0
(19)
T
Ox
OH(O Ou
Equations (12) - (19) define a multi-point boundary-value problem. Since we have only considered necessary conditions, not sufficient, the solution of these equations may yield an optimal value of the cost function. In Section II, a weak formulation of these equations is derived to obtain an approximate, candidate optimal solution.
192
ROBERT R. BLESS AND DEWEY H. HODGES
II. UNCONSTRAINED PROBLEMS In this section, a variational formulation for solving unconstrained optimal control problems is developed. It will be shown that the formulation has no strong boundary conditions (i.e., boundary conditions requiring the virtual, or variational, quantities to be zero) but only natural, or weak, boundary conditions (i.e., those determined by setting the coefficient of a virtual quantity to zero). This allows the test functions (which approximate the virtual quantities) and shape functions (which approximate the dependent variables) to be chosen from a less restrictive class of functions; indeed, the same test and shape functions may be chosen for every optimal control problem. Thus, a general set of algebraic equations will be developed which do not have to be altered in terms of numbers of equations with any changes in the boundary conditions. Hence, this formulation is also referred to as a weak formulation. A. GENERAL DEVELOPMENT The derivation of the weak formulation begins by using the conditions in Eqs. (12) - (16) to simplify Eq. (10), yielding
+ 5uT
Ou
dt
(20)
To put Eq. (20) into the weak form for the finite element discretization, we will integrate by parts so that no time derivatives of x or ~ appear in the weak formulation. As will be shown, this allows the simplest possible shape functions to be chosen and for all element quadrature to be done in closed form. Integrating by parts yields
5uT
OU
dt + E [5xT(t:~ )£(t:~ ) -5£T(t~)x(t~)] (21) i--1
+
N-1
Z [ ~ (t~+)~(t~+) - ~ (t~+)~(t~+)] i=0
A FINITE ELEMENT METHOD
i to
--
At(1 ) t (1)
I
t (2)
tl
II --
I
t (3)
193
I
t (5)
t (4)
I
t2 -- t (6)
Figure 1" Discretized time line with nodes labeled Equation (21), when taken with Eqs. (12) - ( 1 6 ) , is called the weak formulation of the first-order necessary conditions for the optimal control problem defined in Eqs. ( 1 ) - (3). Equations ( 1 2 ) - (16), and Eq. (21) will be used for the finite element discretization scheme described next. B. F I N I T E E L E M E N T DISCRETIZATION
i-
Let the time interval from t +_ 1 to t~- be broken into Mi elements where 1, 2 , . . . , N. For convenience, define i
m-ZMj
for
i-
1,2,...,N
(22)
j=l
and define M0 = 0. This yields a subdivision of the original time interval from to to t f into M N subintervals. The boundaries of these subintervals are called nodes and are denoted by t (i) for i - 1 , 2 , . . . , M N + 1. Note that to = t (1), ti = t (M{+I), and tl - tN : t (MN+I). A nondimensional elemental time T is defined as T
=
t - t (i) t(i+1) - t(~)
=
t - t (i) At(0
(23)
SO that 0 <_ T < 1 within each element. Figure 1 gives a pictorial representation of a two-phase time line with M1 = 2 and M2 - 3. Since no derivatives of x, A, or u appear in Eq. (21), then the simplest possible shape functions in each element, namely piecewise constants, may be used. The benefit of choosing piecewise constant shape functions is that element integration may be done by inspection. The shape functions are
X-
{
~+ if T -- 0; ~i if 0 < T < 1; X~-+I i f T - - 1
(24)
194
ROBERT R. BLESS AND DEWEY H. HODGES
^
~1+
X5
I t (1)
t (2)
1
t (3)
I ~
I t (1)
t (4)
II x~
t(2)
t(3)
I x4
t(4)
I
I
t (6)
t(5)
I
A5
t (6)
t(5)
^
~1+
I t(I)
_
x6
~5
Ul
t (2)
t (3)
t (4)
t (5)
_
u6
I
t (6)
Figure 2: Location of unknowns for x, A, and u
{ {
Aand u-
A+ A~ ~-+1
if "/- - 0; if 0 < T < 1; if T - - 1
(25)
~.+ fi~
if ~- - 0; if0
(26)
Iti---F1
Since first-order time derivatives of 5x and 5A appear in Eq. (21), linear shape functions are chosen for 5x and 5A in each element. Since no time derivatives of 5u appear, then piecewise constant shape functions may be chosen in each element. Again, by choosing the simplest functions allowed, element integration is done by inspection. The functions chosen were 5x
-
5x+(1 - T)
+ (~XhlT
(27)
5A
--
5 A + ( 1 -- T) + (~i-+l T
(28)
and, with 5D defined as the Dirac delta function ~U -- ~U~D(T)-t-~Ui
nL- (~U~+I(~D( T -- 1)
(29)
A FINITE ELEMENT METHOD
195
(~Xl+
(~X2
(~X3 (~X3+
(~X4
(~X5
(~X6
t(l)
t(21
t(a)
t(4)
t(5)
t(6)
I
I
t(1)
I !
II
II
I
I
I
I
I
I
t(2)
t(3)
t(4)
t(5)
t(6)
t(2)
t(3)
t(4)
t(5)
t(6)
+
t(1)
Figure 3: Location of independent variational quantities for 5x, 5A, and 5u
The superscripted "-" and "+" signs in Eqs. (24) - (29) signify values just before and after a particular nodal time t(0. For all nodes except event nodes, the values for x and A, as well as for 5x and 5A, are equal on either side of the node. In other words, only at an event are the states, costates, and their variations allowed to "jump" in value. Thus any discontinuities in the states or costates must be known in advance and modeled with an appropriate event and phase. The examples throughout the chapter will clarify the specifics of defining multi-phase problems. The Dirac delta functions appearing in the discretization of 5u have the effect that the coefficient of the ~u term in the integrand of Eq. (21) is forced to zero pointwise wherever the Delta function has a nonzero value. More specifically, the Dirac delta functions enforce OH/Ou = 0 at the left and right limits of each node. We note that the control could be discontinuous at an internal node even though the states and costates must be continuous. For an example case with two elements in phase one and three elements in phase two, the quantities defined in Eqs. (24) - (26) that will appear in the algebraic weak form are shown on a timeline in Figure 2. The quantities defined in Eqs. (27) - (29) are shown in Figure 3. After changing the variable of integration from t to T with Eq. (23) and
196
ROBERT R. BLESS AND DEWEY H. HODGES
inserting the shape and test functions from Eqs. (24) - (29) into Eq. (21), the element quadratures can be carried out analytically yielding an algebraic equation. Noting that j indexes the phase number and i indexes the element number within a phase, the equation is 0
-+
+
N<
..[_T ( ,, ..~_ + ~A'M'j_ lq--1 -- X'M'j_ 1-1--1 :E'M'j-1-'b1
E j=l
+~ (~X " ~ j _
[~+
~+1
~j
_t..T ~£t~Mj_lnt-1
--1
+1
"~
-
E
2
"M'j- 1"[-1
M'j-1 +1
Atj -
{~)~T (--Xi-1
Atj )
i=Mj_I+2
-3I- ~XT [~i_l +
-1-t-'1
-
A~j_ 1 +1
Mj
+
2
Mj Z
i=Mj_I+I
+
&~j+l
+
(~XMj+1
e<
[(v:):] Atj 2
--~MJ 2
0/:/ -+- (SZ~Mj + 1
(30)
2- -g~- ~-1
"-~
Cv, + ^-
x
XMj +1
-gg
~,
)
- X:~,+I
+1]>
In Eq. (30), note that the superscripted "-" and "+" signs are dropped except at the event nodes where discontinuities may take place in both the state and costate quantities and the variational state and variational costate quantities. A l s o , / ~ - H(2, fi, A), f - f(~, ~), a n d / 2 / _ H(~, fi, A). Additionally, we have taken the elements within each phase j to be of constant width Atj. Finally, in the above equation, 5~+ and 5fi- terms should appear at every node and not just the event nodes; however, the 5fi terms at internal nodes are decoupled from the other ones and may be used to solve for internal nodal values of the control after the other unknowns have been solved.
A FINITE ELEMENT METHOD
197
When the shape functions for x, A, and u as defined in Eqs. (24) - (26) are substituted into the boundary conditions of Eqs. (12) - (16), then the following equations are obtained.
+'
Ot--~+ H~M i + I - /~-~i +1 O0 +[-ti)tN MN+I ^
_
^ - -
Ox(t:~ ) 0(1)
-- )~ i+I--M"
+
^~_T
=
0 (i = 1 , 2 , . . . , N -
=
o
(32)
=
0
(33)
0 (i - 1 , 2 , . . . , N )
(34)
(i = 0, 1 , . . . , N -
1)(35)
--
=
0
1)(31)
Now, in Eq. (30), the coefficient of each of the arbitrary virtual quantities, 5x, 5)~, and 5u, must be set equal to zero in order to satisfy the equation. When the coefficients are set to zero and these equations are combined with Eqs. (31) - (35), the result is a sparse system of nonlinear equations whose size depends on the number of elements. Solution of these algebraic equations yields an approximate candidate optimal solution. A summary of the number of equations and unknown is given next. C. SUMMARY OF EQUATIONS AND UNKNOWNS Consider a problem with N phases, n states (and, hence, n costates), m controls, and q boundary conditions in ¢. For any given phase i, the unknowns appearing in Eq. (30) and Eqs. (31) - (35) are the barred quantities for each element for x, ~, and u, along with the hatted (nodal) quantities at the beginning and end of each phase for x, ~, and u (see Figure 2). Thus, there are a total of (2n + m)(Mi + 2) unknowns for the ith phase. Summing these over N phases results in (2n + m ) ( M y + 2N) total barred and hatted quantities for x, A, and u. In addition to these unknowns, there are q unknown multipliers v and the N unknown times tl, t2, ..., tg. Therefore, the total number of variables to solve for is (2n + m ) ( M g + 2N) + q + Y. For any given phase i, the equations appearing in Eq. (30) and Eqs. (31) - (35) are the coefficients of each independent 5x, 5)~, and 5u when set equal to zero. Since the 5x and 5~ quantities appear at each node (see Figure 3), there are 2n(Mi + 1) equations for the ith phase. Since the 5u quantities appear at the end nodes of each phase, and also at the midpoints of each element (see Figure 3), there are m(Mi + 2) equations. The total number
198
ROBERT R. BLESS AND DEWEY H. HODGES
of equations then from phase i is m + (2n + m)(M~ + 1). Summing these over N phases results in m N + (2n + m)(My + N) total equations from the coefficients of 5x, hA, and 5u. Additionally, Eqs. (31) - (32) provide N equations, Eq. (33) provides q equations, and Eqs. (34)- (35) provide 2nN equations. Thus the total number of equations is (2n + m)(MN + 2N) + q + N, which agrees with the number of unknowns. The algebraic system of equations is typically solved by expressing the Jacobian explicitly and using a Newton-Raphson solution procedure. For the example problems which follow, the iteration procedure will converge quickly for a small number of elements with a relatively simple initial guess. Then, the answers obtained for a small number of elements can be used to generate initial guesses for a higher number of elements. When the sparsity of the Jacobian is exploited, it is possible to efficiently generate a reasonably accurate solution. Although the nonevent nodal values for xi and Ai do not appear in the algebraic equations, their values can be easily recovered after the solution is found. This is most easily seen by looking at the following ordinary differential equation multiplied by a test (or weighting) function 5A and integrated over some time interval where the integral makes sense.
ftl 6A[f (x, t) - 5] dt
(36)
0
After an integration by parts, using the linear shape function for 6A defined in Eq. (28), using the piecewise constant shape function for x defined in Eq. (24), and substituting T for t as given in Eq. (23), then the following equation is obtained from Eq. (36). ~/~1
--:~1 ~- "-~'Y nt- Xl
-{- (~/~2
:~1 -~- - ~ - f -- :E2
-- 0
(37)
With arbitrary ~A1 and ~A2, the coefficients must vanish, forming two equations of the form AtAt-
-
0
-
0 o
Now, by subtracting the second equation from the first, it is seen that _ = ~1 + ~2 Xl
2
(38)
or, in words, that the interior value (the bar value) is simply the average of the surrounding nodal (or hatted) quantities. Once a solution to the
A FINITE ELEMENT METHOD
199
algebraic equations is found, all the midpoint values and the end nodal values are known, and thus all other nodal values can be recovered by repeatedly using Eq. (38) for both x and A. In view of the redundancy, only the nodal values are plotted for the states and costates in this chapter. As for the control, interior nodal values are available once the nodal values for the states and costates are found by using the optimality condition (OH/Ou = 0) to solve for u at nodes. Recall that these equations appeared in the weak formulation but were decoupled from the other equations. We note that the midpoint value for the control is not the average value of the surrounding nodal values; thus, both midpoint and nodal values are plotted in the figures in this chapter. D. EXAMPLE: A SINGLE-PHASE P R O B L E M As the first example problem, the transfer of a particle to a rectilinear path will be examined. This example is taken from section 2.4 of [4]. The purposes of this example are to (1) explicitly state the boundary conditions and algebraic equations of the weak formulation for a one-phase problem, and (2) present some results which show the accuracy of the method for this problem. Let x(1) and x(2) denote the horizontal and vertical position coordinates, respectively, of the particle at a given time and x(3) and x(4) denote the particle's velocity components at a given time. (A subscripted number in parentheses refers to the state index to avoid confusion with the element index.) The thrust angle u is the control and the particle has mass m and a constant acceleration a. The state equations are defined as X(3) --
X(4)
a cos u a sin u
= f
(39)
The problem is to maximize the final horizontal component of velocity in a fixed amount of time t l = T. Thus, L = 0
and
¢ = x(3)(tl)
(40)
As for boundary conditions on the states, the initial conditions are zero and there are two terminal constraints. These are that the particle arrive with a fixed final height h and that the final vertical component of velocity be zero. The final horizontal component of position is free.
200
ROBERT R. BLESS AND DEWEY H. HODGES
For single-phase problems, Eqs. (30) - (35) are greatly simplified. Note that since there are no internal events, the "+" and "-" superscripts may be dropped. Looking at the boundary conditions first, Eq. (a3) becomes x(1)(to)
x(2)(to)
x(3) (t0) x(4) (to)
=0
x(2)(tl)-h
(41)
x(4) (tl) tl - T
where the conditions in ¢ are adjoined to the cost with multipliers Vl, /22, ..., vT. Equation (31) is no longer applicable and EQ. (32) becomes ^
/27 + H~'a+ 1 - 0
(42)
The costate boundary conditions in Eqs. (34) and (35) become 0
/21 1 nt-
-- 0
/22 /23 /24
^ ~+1
and
/25 1
-
-0
(43)
/26
Finally, Eq. (30) is simplified to
0
=
e~l~
q-
/22"= ~ / T
+
exT
i--1
-~1 + ~1
(
~-fl
At-
Xi-1----2-
At
I~-1 + ~ - - V ~
--:~/-1
-~
i-1
-
(
at ( 0 ~ )
)
+ &xT~+I --~M~ ~-fM~+ a?Mx+l
+ &5~+1 XMI-- 5- -0~ M1--
)
'~MI+I]
(44)
A FINITE ELEMENT METHOD
+
~uTI q-i
-~
201
Mlq-I
where At = (tl - to)/M1 = tl/M1 = T/M1 for all elements. These equations are solved by expressing the Jacobian explicitly and using a Newton-Raphson algorithm. For M1 - 2, suitable initial guesses for the nonlinear iterative procedure can be found by simply choosing element values that are not too different from the boundary conditions. The results from solving the M1 = 2 equations are then used to obtain the initial guesses for arbitrary M1 by linear interpolation. Representative numerical results for all four states versus dimensionless time t / T are presented in Figures 4 - 7. For this example, h = 100, T = 20, and 4 h / a T 2 - 0.8897018. This last number is chosen to yield a value of 75 ° for the initial control angle of the exact solution available in [4]. The results for 2 and 8 elements are plotted against the exact solution. It can easily be seen that M1 -- 8 gives acceptable results for all the states. Even the very crude 2-element mesh yields a decent approximation to the answer and led to convergent initial guesses for the 8-element run. In Figure 8, the thrust angle u versus dimensionless time t / T is presented. Once again, the results are seen to be excellent for M1 - 8. Note that instead of just nine data points (circles) for the nine nodal values which were plotted for the states, there are 17 data points available for the control since both midpoint and nodal values are obtained from repeated use of the optimality condition. This is of great value since it is the control variable which is generally of the most interest. Three of the four costates are constants for all time and this method yields two of these exactly. The third costate is very close to the exact answer. The fourth costate, A(4), corresponding to the vertical component of velocity is shown in Figure 9. The results compare nicely with the exact results. A plot of the relative error of the performance index J - x(3)(T) and the initial control angle u(0) - Ul versus the number of elements is shown in Figure 10. For 8 elements or more, the lines are nearly straight on a loglog scale. The slope of the line is about - 2 which indicates that the error varies inversely with the square of M1, similar to a-posteriori error bounds as formulated in usual finite element applications [46]. Notice in Figure 10 that there is a bend in the initial thrust angle curve. It is not unusual for mixed formulations to have an error curve that is not monotonically decreasing. It should be noted that developments of mathematical proofs for convergence and expressions for error bounds are not state-of-the-art for mixed methods.
202
ROBERT R. BLESS AND DEWEY H. HODGES
1.2-
1.0c-
FE(2)
•
FE(8)
.
/
Exac
C i
•
O n
{D
0.8-
O 0_
-
n
C O N O
7"
0.6_ O.4-
0.2-
0.0
q
0.0
0.2
0.4
0.6
0.8
Dimensionless Time Figure 4" Normalized horizontal position versus
t/T
1.0
A FINITE ELEMENT METHOD
203
1.00.90.8-
v ,,
FE(8)
-
-
Exact
0.7 -
P,
:,~
FE(2) /
J
~
0.6-
13_ 0.5o4>
0.3 0.2 0.1
O0 1 0.0
'
I
0.2
'
I
0.4
'
I
'
0.6
I
0.8
Dimensionless Time Figure 5: Normalized vertical position versus
t/T
'
I
1.0
204
ROBERT R. BLESS AND DEWEY H. HODGES
3.0-
2.5-, O
• •
~'~ 2.0-
~
FE(2) FE(8)
//~ ~
•
1.5-
N
1.0-
-r"
0.5-
0.0
l"
I
0.0
'
I
0.2
'
I
0.4
'
I
0.6
'
I
0.8
Dimensionless Time Figure 6: Normalizedhorizontal velocity versus t/T
'
I
1.0
Vertical Velocity/(h/T) 0"q
0
cb
0
--,,1
0
0
~..,o
0
CI)
CL .,¢ c--F C~
o o I.-,L"F
.<
C'D
CI)
o
b0
0
o ±
0
~
~
0
I
J
~
0
I
i
~
0
I
~
b0 I
~
J
o
~
I
t
ro I
~
;
.~
I
206
ROBERT R. BLESS AND DEWEY H. HODGES
80-
60-
40-
v
20-
m
O
m
-20 -
-40 -
•
FE(8) Exact
~
-60 ip
-80
'
0.0
I
0.2
'
I
0.4
'
I
0.6
'
I
0.8
Dimensionless Time Figure 8: Thrust angle versus t/T
'
I
1.0
A FINITE ELEMENT METHOD
4
~
207
!
3
2 0
1
"~
o
-2
-3
9
-4
'
0.0
I
0.2
'
I
0.4
'
I
'
0.6
J
0.8
Dimensionless Time Figure 9: Vertical velocity costate versus
t/T
'
I
1.0
208
R O B E R T R. B L E S S A N D D E W E Y H. H O D G E S
100 -
1 0 "1 - _
0
I,..
LLI .
(1) >
.
Cost
1 0 -2 - -
m
n
rc 1 0 .3 - -
Initial Thrust Angle
10 .4
I
1
I
I
I
I
I
I I ]
1
I
I
I
I
I I I I
10
Number of Elements Figure 10: Relative error for cost and initial thrust angle versus number of elements
100
A FINITE ELEMENT METHOD
209
E. EXAMPLE" MULTI-PHASE ASCENT TO ORBIT A two-stage launch vehicle will now be examined with the three purposes of (1) providing an example of a multi-phase problem definition, (2) examining the finite element accuracy for highly nonlinear problems, and (3) more closely examining how the Newton method behaves when using linear interpolation to increase the number of elements. Consider the following model of a rocket flying over a spherical, nonrotating Earth with the states chosen to be mass m, altitude h, speed V, and flight-path angle 7, and the control to be the angle of attack a. The dynamical equations are approximated to be
Tvac gI~p
=
h l)"
-
(45)
ysin7 Tcosc~-D m Tsina+L mV
=
(46) #sin7 r2 +
(V
r
(47) it ) r2V cos 7
(48)
where T - Tva c - A e p , T,a~ is the thrust in a vacuum, Ae is the exit area of the nozzles, p is the local atmospheric pressure, Isp is the specific impulse, g is acceleration due to gravity at sea level, # is the Earth's gravitational constant, and r is the distance from the center of the Earth to the vehicle, given by the radius of the Earth Re + the altitude h. Defining p as the density, S as the reference area and Ca and CN as the axial and normal force coefficients, respectively, then the dynamic pressure, q, the drag, D, and the lift, L, are defined as
q-
1 -~pV 2
(49)
D
-
aS(Ca 2t- CNOL 2)
(50)
L
-
q S ( C N - C~)a
(51)
The numerical values for the physical constants are # - 3.986 x 1014 m3/s 2, Re = 6378000 m, and g - 9.81 m/s 2. The propulsive model used was
Tvac
=
1.2975 z 107 N before staging
(52)
Tvac Ae
=
2.5950 × 106 N after staging
(53)
=
19.115 m 2 before staging
(54)
210
ROBERT R. BLESS AND DEWEY H. HODGES
Ae
=
3.823 m 2 after staging
I~p
=
430.55 s
(55) (56)
and the atmospheric and aerodynamic model was kg/m 3
(57)
p
=
1.225exp(-h/8600)
p
=
101325exp(-h/7600) Pa
(58)
S
=
55.18 m 2
(59)
Ca
--
0.3
(60)
CN
--
3.1
(61)
We note that constant values for the aerodynamic coefficients are unrealistic; however, this assumption does simplify the numerical operations needed to solve the problem. Also, for numerical simplicity, all aerodynamic and atmospheric terms were neglected after staging, since they are extremely small by this point in the trajectory. The performance index is defined as J = ¢ = - m ( t f ) , and the final time tf is free. Defining m0 = 890150 k.g as the initial mass, hf = 148160 m as the final altitude, and VI = 7854 m/s as the final speed, then all the state and time boundary conditions required for the finite element formulation can be stated as re(to)
-
mo
h(to) - 60
V(to)-25 ")'(to)- 1.5 re(t+) - r e ( t 1 ) +
¢ -
29920
h(t+) - h ( t l )
- 0
v(t+l)- V(tl) ~(t+l ) - ~ ( t l ) h(t~) - h~ v(t~) - vI ~(t~)
(62)
tl -- 195 Note that (1) t l is the specified time of staging, (2) there is a discontinuity in mass of 29920 kg at staging, and (3) continuity of the other three states is explicitly stated. By adjoining the conditions in ~ with multipliers ui for i = 1, 2 , . . . , 12, then Eqs. (31)and (32) become
~1~ + / t ~ + 1 -/~+~,+1
= 0
(63)
--
(64)
^
H-~2 + 1
0
A FINITE ELEMENT METHOD
211
The costate b o u n d a r y conditions in Eqs. (34) and (35) become Pl
+ +
u2
- 0
/23
and
--/]5 ~Y6
^
A~+I
-
-
P4 and
^ q-
nl-
-
0
(65)
= 0
(66)
--/]8
u5
AMI+I
-
--/27
/26
u7
-1 --
0
and
^ A~M2+ 1 --
/28
v9 Ul0 /211
Since the internal boundary condition was based on a time of staging instead of a particular state value (such as mass), there are no jumps in the costates. An example with multiple stages and jumps in the costates will be given in Section IV. For numerical conditioning, the states were normalized in a simple manner by defining ~ - m/mo, h - h / h f , and V - V / VI. Normalizing the states automatically scales the costates. The numerical results which follow were generated with the code VTOTS which is described in [16]. Obtaining guesses for a small number of elements is generally easier than for a large number of elements. Additionally, once a solution has been found for a few elements, more dements may be added readily in most cases. For example, Table I was produced after a 2-element-per-phase solution was found (a very crude approximation). By linearly interpolating the solution, new initial guesses were computed for a 4-element-per-phase solution. These new guesses yielded an average initial error of the algebraic equations of 2.60015 x l0 -2, as displayed in the second column. The last column of Table I displays the total number of iterations, using a Newton method, required to obtain a converged solution. The "4:4" solution was then used to obtain initial guesses for the "8:8" case, and similarly, the table was completed. It is seen that the initial error monotonically drops, as does the number of iterations required for convergence. Doubling the number of elements has proven to be a quick and easy way of obtaining reasonably accurate answers. Table I. Build up of finite element values Elements Per Phase 4:4 8:8 16:16 32:32
Average Initial Error 2.60015 x 10 -2 9.93419 x 10 -3 2.35328 x10 -3 4.81494 x 10 -4
Number of Iterations 5 4 3 ' 3
212
ROBERT R. BLESS AND DEWEY H. HODGES
Some of the results are shown in Figures 11 - 18. In all of the figures, finite element results and shooting results are shown. The finite element results are taken with 2 elements per phase and 8 elements per phase and are denoted with "FE(2:2)" and "FE(8:8)", respectively. Figure 11 shows the normalized altitude profile, and Figure 12 shows the normalized velocity. It is seen that the 2-element case only gives a rough approximation to the answer, and that the 8-element case lies essentially on top of the shooting curves. Figure 13 shows the flight-path angle, V, versus time. The vehicle is seen to roll over quickly to take advantage of the gravity field. Notice that the finite element results show some error for this case. The error is due to an insufficient number of elements being used to model the rapidly varying 7 profile at the initial time. A finer mesh increases the accuracy in this area. Figures 14 - 17 give the four costate histories for the vehicle. The finite element results behave similarly to the state histories; namely, that the 2-element results are a rough approximation, and that the 8-element results are quite accurate except in areas where the variables change more quickly than the finite element mesh is able to model. Examples of this are seen in Figures 15 and 17 for the altitude and flight-path angle costate, respectively. The angle of attack history is shown in Figure 18. The control experiences a very rapid change from a negative angle of attack to a positive one at the initial time. The 2-element case yields ridiculous results initially with an initial angle of attack of -81.9°! However, since the rest of the 2element control history is reasonable, it is still possible to use that answer to obtain guesses for a finer mesh. The 8-element case is seen to be much more accurate with an initial control angle o f - 2 3 . 5 ° , however, it is still not very accurate in the first stage where the control varies rapidly. Figure 19 shows the relative error of the midpoint values of control for varying meshes. In all cases, we have used 16 elements in the second phase, and varied the number of elements in the first phase. It is seen that the control error typically decreases with more elements, although there are singular places where this is not true; and that the largest error is at the initial time. These examples conclusively demonstrate the accuracy and efficiency of the present methodology for unconstrained problems. The method provides virtual self-starting capability and a rapid convergence of states, costates, and controls as the mesh is refined. These are both important for potential real-time applications. III. T R E A T M E N T
OF CONTROL
CONSTRAINTS
This section is concerned with the treatment of optimal control problems
A FINITE ELEMENT METHOD
213
-
1.0~
.~
v,
~-
0.8
.~~ o.6 N ¢
0.4
:2) FE(8:8)
•
t~ng
Z 0.2
0.0
"
._...~ I
I
100
'
I
200
'
I
300
Time (s)
Figure 11" Normalized altitude versus time
'
I
400
'
I
500
214
ROBERT R. BLESS AND DEWEY H. HODGES
1.00.90.8•~,
a
m
0
0.70.6-
"~
0.5-
~
0.4-
Z
0.3-
N
-
•
FE(2:2)
•
FE(8:8) Shooting
0.20.10.0 T 0
'
I
100
'
U
200
~
I
300
Time (s) Figure 12" Normalizedvelocityversus time
'
I
400
'
I
50O
A FINITE ELEMENT METHOD
215
90-
80-
70-
ED (D
60-
_.m
50-
v
ED c"
<
J.4=,.J :
Q. !
.
tED
.
m
40FE(2:2)
30-
FE(8:8) Shooting
20-
m
LL 10m
-10 0
'
I
100
'
I
200
'
I
300
Time (s) Figure 13: Flight-path angle versus time
I
400
I
50O
216
R O B E R T R. B L E S S A N D D E W E Y H. H O D G E S
o
-0.1 -~ -0.2 -0.3
(D a3
r.D
O
O
-0.4 -0.5
{D
09
-0.6 FE(2:2) FE(8:8) Shooting
-0.7-0.8-0.9-1.0
'
0
I
100
I
I
2O0
'
I
I
3OO
4O0
Time (s) Figure 14" Mass costate versus time
I
v
I
5O0
A FINITE ELEMENT METHOD
217
0.0-0.02 -0.04 -0.06 -0.08
0
-o.1
"0 :3
-0.12-
<:::
-o.14-
•
FE(2:2)
•
FE(8:8) Shooting
-0.16-0.181
-0.2 '
0
I
100
'
I
200
I
I
300
Time (s) Figure 15: Altitude costate versus time
i
i¸
I
400
5OO
218
ROBERT R. BLESS AND DEWEY H. HODGES
-0.1 -0.15 -0.2
(D
{D O
(O m
m
-0.25 -0.3 -0.35
•
O (D
>
FE(2:2) FE(8:8)
•
m
Shooting
-0.4 -0.45 -0.5 -0.55 -0.6
I
0
'
I
100
'
I
200
'
I
300
Time (s) Figure 16: Velocity costate versus time
'
I
400
I
5OO
A FINITE ELEMENT METHOD
219
.10 -~ 1.0-
0.5-
FE(2:2) FE(8:8) Shooting
0.0
0
(1)
-o.5
~C:D -1.0
t'3 !
-1.5-
c" C:)') •--2.0IJ_ -2.5-
-3.0 0
'
I
100
'
I
200
'
I
300
Time (s) Figure 17: Flight-path-angle costate versus time
'
I
400
'
I
500
220
ROBERT R. BLESS AND DEWEY H. HODGES
20-
O
m
"0
v
0
-20 -
< 0
FE(2:2) FE(8:8) Shooting
-40 -
m
c-
<
-60 -
-80 -
i
0
I
IO0
i
I
20O
i
I
3OO
Time (s) F i g u r e 18" A n g l e of a t t a c k v e r s u s t i m e
I
400
I
5O0
A FINITE ELEMENT METHOD
221
100 8 Elements I
1
0
"1
t
-_
"",,,
%
~ \
",,
%%
\
_
16 E l e m e n t s
"',,,
64 Elements
',
\ ,,
_
~L_
\ ,,
0
UJ (D >
32 E l e m e n t s
""
\ \, o° i \
1 0 "2 - _
. , m
_
\
,s
\ %~ \
\
\
\
\
%%% ",~,% i - - - - . . .
°°
s "
,~
I
\
I
m
(D rr
oo
oo
\
\
,"
l 1 0 -3 - -
\\
I I
10"
' 0
I 50
'
I 1 O0
\
'
I 150
Time (s) Figure 19: Relative error in angle of attack versus time
'
\\\
I 200
222
ROBERT R. BLESS AND DEWEY H. HODGES
that are subject to inequality constraints which are explicit functions of the controls. A. GENERAL D E V E L O P M E N T Consider the same system as defined in Eqs. (1) - (3). For simplicity, it is assumed that the problem has only one phase. Suppose now that the p constraints are written as
gi(x, u) ~_ 0
for
i = 1, 2 , . . . , p
(67)
One way of handling inequality constraints is to use a "slack" variable [47]. The idea is that if gi _< 0 then gi plus some nonnegative number, the ith slack variable k 2, is equal to zero. Thus, we form p x 1 column matrices for K and 6K
K6K
Lk~ k~ =
[2~1 (~kl
...
k]j T
2]{2(~2
...
(68)
2kp(~kpJT
(69)
and a p x 1 column matrix g(x, u), the elements of which are gi(x, u). Now, from Eq. (67) g(x, + K : 0 (70) Eq. (70) will be adjoined to the augmented cost function Ja in Eq. (6) by using p Lagrange multiplier functions #(t). After dropping the superscripted phase number on H, we obtain
J~ - • +
fl 1 [ H -
IT9 + #T (g + K)] dt
(71)
By redefining the Hamiltonian in Eq. (5) to be
(72)
H =_ L + )~Tf + #T(g + K)
then, after the usual integration by parts and making use of Eqs. (8) - (9), we obtain the first variation of Ja as
61T (f -- x) + ~xT
--
+
5u T
~
-~x
+spT(g+K)+6KT#
+ i
dt
(73)
A FINITE ELEMENT METHOD
+
+
223
~ 1 + H ( t l ) d r 1 + dp T ¢ [ 0(I)
OX(tl)
)~T(tl)]
dx(tl) + Ox(to) + AT(t°) dx(to)
Setting 5Ja = 0 in Eq. (73) results in the same necessary conditions (for a single-phase problem) as are listed in Eqs. (12) - (19). There are also the additional conditions that gi+k 2-0 #iki-O
~ J
for i - 1
'
2,
"" " P
(74)
These conditions, when satisfied, ensure that the constraint will not be violated at element midpoints and nodes and that either the constraint or the multiplier will be zero. For example, if the constraint gi equals zero, then ki must be zero, and so #i need not be zero. On the other hand, if gi is less than zero (inactive), then ki must be nonzero, which in turn forces pi to zero. An additional requirement for a minimizing problem is that the multipliers p must be nonnegative for all time. This is enforced computationally, as will be shown in the example, by squaring the value of p that appears in the costate equations and optimality conditions. There may be numerical problems associated with the introduction of slack variables. As pointed out in [1], there can be poor conditioning or even singularities in the Jacobian occurring during the solution procedure if both #i and ki were zero at the same time. The advantage of using a slack variable, however, is that control-constrained problems may be treated as single-phase problems. Since the control is allowed to come on and off the constraint as required, there is no need to define a switching structure (i.e., the sequence of constrained and unconstrained arcs) for the problem. To derive the weak formulation, Eq. (73) is simplified by using the boundary conditions in Eqs. (12) - (16). Then an integration by parts to eliminate any time derivatives of x or A from appearing in the formulation yields
0
--
+ +
5)~Tf + 5~Tx + 5xT
5u T
~
~X
+spT(g+K)+SKT# - e,
r(tl)z(tl) +
-- (~2T/~
dt
(75)
ar(to)z(to) - e z r ( t o ) a ( t o )
Equations (12)- (16) and Eq. (75) define the weak formulation for optimal control problems with control constraints. Equation (75) can be generalized in a straightforward manner to multiple phases, similar to Eq. (21).
224
ROBERT R. BLESS AND DEWEY H. HODGES
The discretization of the weak form is similar to what is defined in Section II, B.
B. FINITE ELEMENT DISCRETIZATION To the discretization defined in Section II, B, we define shape functions for 5K and 5# of the form (with 5D defined as the Dirac delta function) 5K
-
2k+Sk+~D(T)+ 2ki6ki + 2ki--+l(~ki--+l(~D(T- 1)
~lZ
--
6ft+6D(T)
nt- ~fti + (~i-+l(~D(T
--
1)
(76)
(77)
where the subscript i now refers to the element number. The Dirac delta functions appearing in the discretization of 5K and 5# have the effect that the associated parts of the integrand of Eq. (73) are forced to zero pointwise wherever the Delta functions have nonzero values. The shape functions for K and # are of the form K-
]c + ~2
if T = O; if 0 < ~- < 1;
(78)
"2--
ki+l if W-- 1 and
if # =
/~ fti-+l
-o;
if 0 < T < I;
(79)
ifT--1
The superscripted "-" and "+" signs in Eqs. (76) - (79) signify values just before and after nodal times. By substituting the shape and test functions defined in Eqs. (76) - (79) and the shape and test functions defined in Eqs. (24) - (29), and carrying out the element quadrature over 7 from 0 to 1, the general algebraic form of the weak principle is obtained. The equation is similar to Eq. (44) with the addition of the appropriate 6# and 6K equations. These equations will be written out explicitly for the example problem. Note that if the time t does not appear explicitly in the problem formulation then all integration is exact and can be done by inspection. If t does appear explicitly, then t may be approximated by its midpoint value over each element and the integration may still be done by inspection. This is equivalent to redefining t to be a new state. The latter case occurs in the example problem presented next. In addition to the unknowns and equations listed in Section II, C, there are an additional 2p(M1 + 2) equations corresponding to the coefficients of
A FINITE ELEMENT METHOD
225
~p and 5k and an additional 2p(M1 + 2) unknowns to solve for the k, ]~, #, and/2 quantities. Thus, there are still the same number of equations as unknowns. Again, once the nodal values for the states and costates are found, then one may use the optimality condition (OH/Ou = 0), the constraint equations (g+K = 0), and the condition that either k or p be zero (k# = 0) to solve for u, K, and p at a node. This procedure is used in the following example problem. C. EXAMPLE This example is taken from [4], section 3.8. The problem is to minimize J-~
1
x(T) 2 +~
1 fT
]0
u2dt
(80)
where T = 8, x and u are scalars, and the initial condition is x(0) - ( 5 + 8/~). The state equation is
2. = h(t)u
h(t) - v/2sin
with
-~-
The following two control inequality constraints are imposed. gl
=
u-1
g2
--
-(u+
50
(81)
1) _< 0
(82)
The exact solution is found to be x(T) - - 1 and for the control
-x(T)h(t) 1
u(t)-
-x(T)h(t) -1
-x(T)h(t)
for0_
(83)
The algebraic equations which come from the weak principle can be verified to be Xl -+- 5 -t- 8/71"
--
0
(84)
t l -- 8
--
0
(85)
for the state and time boundary conditions (in ¢), ~1 q-/Jl
--
0
(86)
226
ROBERT R. BLESS AND DEWEY H. HODGES ~1 -- X1
--
0
(87)
Xi -- Xi+I
--
0
'~M1 -- '~MI-t"I
--
0
(89)
'~MI-}-I -- X M I + I
--
0
(90)
for i - 1, 2 , . . . ,
M1
-
1
(88)
for the costate equations (obtained from the coefficients of 5x set equal to zero) and boundary conditions, ~ti "-~ ~ih ( t i ) + fitli -- fit2i ~-- 0
k1~#1~ - 0
k+#2~ = 0 fii - 1 ~2 _ 0 -~i-l+k!: 0
for i - l , 2 , . . . , M ~
(91)
for the coefficients of ~fi, 5k, and ~#, (note that the equations formed from setting the coefficients of 55, 5k, and ~ at the first and last node of the phase equal to zero are similar to these) and :~1 X_i + l
_ -- Xi
At
~--h (t-l)Ul
-
Xl
(92)
A t [h (t-i)ui + h ( t i + l ) Ui+I]
--
0
(93)
-
1,2,...,M1-1
--
0
for i At --:~M, -- --~-h (t-M1)UM, ~- :EMI-t-1
(94)
for the 5A coefficients. Note that t-i is an average time value for the ith element and (if At - Ati - tl/M1 for all i) can be expressed as 2i- 1 t-i = ~ A t 2
for i - 1, 2 , . . . , M1
(95)
Recall from [4] that one of the additional necessary conditions for problems with control constraints is that the multipliers be nonnegative for a minimizing problem. Therefore, in practice, the multipliers p appearing in the first of Eqs. (91) are squared to ensure that they remain nonnegative. Further recall that if the constraint is not violated, then p - 0. This condition is satisfied by the second and third equation in Eq. (91) which implies that either k or p is zero for each element. It is readily apparent from Eqs. (86) - (90) that all the costate variables are equal to ~?M1+1 -- x(T). Therefore, these equations were eliminated and all costates that occurred in the remaining equations were replaced with i Ml+l- Also, the state and time boundary conditions in ¢ were solved. The remaining 6M1 + 1 algebraic equations in Eqs. (91) - (94) were solved
A FINITEELEMENTMETHOD
227
using a Newton-Raphson method and a FORTRAN code written on a SUN 3/260. The sparse, linearized equations were solved using subroutine MA28 from the Harwell subroutine library [48]. This subroutine takes advantage of sparsity which leads to great computational savings. Table II shows the convergence rate of XMl+l = x(T), the elapsed computer time for the first 5 iterations, and the percentage of zeroes in the Jacobian (i.e., the sparsity) versus the number of elements. The x(T) column shows that the 32-element case has almost converged on the exact solution. (Recall that the exact solution is x(T) = - 1 . ) Note further that the approximate x(T) is not an upper or lower bound of the exact value, which is common in mixed formulations. The third column of Table II gives the elapsed computer time for five iterations. It is easily seen that there is a modest increase in computer time with an increase in the number of elements. Note that in some cases acceptable convergence criteria were satisfied in fewer than five iterations. This is because the answers obtained from a small number of elements (say 2 or 4) may be interpolated to generate initial guesses for a higher number of elements. Finally, the extremely sparse structure of the Jacobian is demonstrated in the last column. This strongly encourages the use of a smart sparse matrix solver such as MA28. This subroutine leads to quicker solutions and tremendous savings in memory allocation since only the nonzero elements of the J acobian need be stored. T a b l e II. x(T), e l a p s e d c o m p u t e r t i m e a n d p e r c e n t s p a r s i t y of J a c o b i a n v e r s u s t h e n u m b e r of e l e m e n t s M1
M1
x(T)
1 2 4 8 16 32 Exact
-7.5465 -0.4439 -0.8385 -1.0685 -1.0141 -1.0034 - 1
Time (sec) 0.42 0.44 0.66 0.76 1.03 1.52
Sparsity (%) 65.3 80.5 89.6 94.6 97.3 98.6
Results for the state and control are shown in Figures 20 and 21 for 4 and 8 elements and the exact solution. Note from Figure 21 that the 4-element case does not define the constraint boundaries very accurately; however, it is accurate enough to generate guesses for the 8-element case. Thus, in a problem with many constrained and unconstrained arcs, a small number of elements could still be used to generate guesses for a higher number of elements.
228
ROBERT R. BLESS AND DEWEY H. HODGES
One of the constraint multipliers, #1, is shown in Figure 22. Note t h a t #1 = 0 when u < 1, so only half the multiplier history is shown. It can be seen t h a t both the 4- and 8-element cases yield very poor approximations for the multiplier, indicating the need for more elements. However, from the control history plot, the 8-element answer seemed very accurate. This highlights the need to check all results and verify tbnt the multiplier histories seem reasonable. Note also that, due to the s y m m e t r y in this problem, the nonzero segment of #2 is identical to #1.
IV. TREATMENT OF STATE CONSTRAINTS A more difficult constrained problem is one in which the constraints are functions of the state variables only. As will be seen, state constraints produce internal boundary conditions on the states leading to discontinuities in the costates. Choosing the number and type of internal b o u n d a r y conditions can be a nontrivial, trial-and-error task for the analyst. A. G E N E R A L C O M M E N T S Consider now the case where there is a state constraint of the form
S(x) < 0
(96)
imposed on the problem. W i t h the imposed constraint, the optimal trajectory m a y (1) never violate the constraint, (2) touch the state constraint at one or more instances in time (touch points), (3) remain on the curve S = 0 for one or more finite amounts of time (a b o u n d a r y arc), or (4) have some combination of boundary arcs and touch points. The sequence of constrained arcs, unconstrained arcs, and touch points will be referred to as the switching structure. To define the order of the constraint in Eq. (96), successive total time derivatives of S are taken and f(x, u) is substituted for ~ until an expression explicitly dependent on the control u is obtained. If p total-time derivatives are required, then S is called a pth-order state-variable inequality constraint 1. One problem facing the analyst trying to solve a state-constrained problem is finding the switching structure of the constrained arcs. Some general guidelines are given in [49] and [50]. Assuming t h a t the hodograph of the 1One could allow S - S(x, u) so that p = 0. This constraint could then be treated in the framework of this section as opposed to that of Section III. This may provide some advantages in certain situations.
A FINITE ELEMENT METHOD
O
229
~
-1-
-2-
C) t
•
FE(4)
•
FE(8)
-3-
m
"1"(D
-4-
,4..~
09
-5-
-6-
-7-
-8 0
1
2
3
4
Time Figure 20: State history versus time
5
6
7
8
Control History s
6
b
i,--, •
crq
b
6
b
6
~
6
b
o
b
o
b
o
o
~
b
('b
bO
0
Do
0 i--,. ¢-P
o
Or3
=
Or3
CD
c.o
--I
3
j~
oi
--~
o~ -4
Go
-
•
m
x o
mm
m oo
m j~
A FINITE ELEMENT METHOD
231
0.6-
0.5-
FE(4) FE(8) FE(16) Exact
~ " 0.4 0
•
ll/ /
\n \
•
. , m
"r-
~- 0.3 a
~
0.2
0.1
0.0 0
1
2
Time Figure 22" Half of the multiplier history versus time
3
4
232
ROBERT R. BLESS AND DEWEY H. HODGES
problem is strictly convex [51] (also referred to as having a regular Hamiltonian), then it has been shown that the control must be continuous across a junction point (i.e., points where constrained and unconstrained arcs are joined). In this case, for a scalar state constraint, it is proved in [49] that only touch-point solutions are possible for odd-order constraints of order greater than two. Most problems with practical applications, however, involve first- or second-order constraints. In [50], it is shown that first-order constraints will only exhibit boundary-arc solutions, except under three special conditions. One condition is the trivial case, where the constraint is identically satisfied by the unconstrained solution. The second condition arises when there is a corner in the hodograph. The third case is if the pth derivative of the constraint is a global extremum when evaluated at the optimal state and control values. This unusual case is demonstrated in the third example of this section. These three conditions are quite uncommon so the analyst can reasonably assume, at least initially, that there will be a boundary-arc solution for first-order constraints. For a second-order constraint, it is not uncommon that there will be some touch-point cases and some boundary-arc cases depending on the constraint limit. This behavior is demonstrated in the first example. 1. Touch-Point Case For the special case where the trajectory only touches the curve S = 0 at discrete instances of time, then the problem defined in Eqs. (1) - (6) is of the proper form. For example, if there were only one touch point occurring in the solution at some unknown time, then the conditions that S(tl) - 0 and the continuity of all states at time tl, (i.e., x ( t ~ ) - x(t+l) - 0), would be added to the vector ¢ of Eq. (3). Note also that the problem would now have one additional phase due to the extra event at t l. 2. Boundary-Arc Case Consider the situation where the solution has one unconstrained arc between to and t l, followed by a constrained arc between t l and t2, and then another unconstrained arc between t2 and tf (or t3). The augmented cost function defined in Eq. (6) is modified in the second phase to include the state constraint. As is done in [4] and [52], successive time derivatives of the state constraint are taken until the control appears explicitly, say in the pth derivative. Then, to enforce the condition
dPS dtP
=0
(97)
A FINITE ELEMENT METHOD
233
for time between t l and t2, the term
dPS ~7 dtP
(98)
is added to the integrand of Eq. (1) and the augmented cost function of Eq. (6) takes the form
N~I: ( Ja - (~-~-i~ -, H(i)
__ , ~ T ~ +
dPS'~
? ~ . ~ ) dt
(99)
Here, r/is a new multiplier function of time which is nonnegative for tl+ _< t _< t 2 and is identically zero for t + _< t < t I and t + _< t < t 3. By redefining the Hamiltonian in Eq. (5) to be
dPS H - L + )~Tf + rl dtP
(100)
then we see an analogy between the development of the weak formulation for state constraints and that for control constraints as is given in Section III, A. The test and shape functions for 57/and r/will be comparable to those given for 5# and # in Eqs. (77) and (79), respectively. Continuity conditions for the states at t l and t2, namely X(tl ) - x(tt )
-
o
x ( t 2 ) - x(t +)
-
0
(101)
and tangency conditions on the state constraint S
T[x(t+)] = 0 where
T - IS
~d--t
"'"
~,,-1~]~ dtP -1
(102)
(103)
must be satisfied. These internal boundary conditions are enforced by adding these conditions, Eqs. (101) - (102) to the constraint vector ¢ of Eq. (3). B. EXAMPLE: A SECOND-ORDER PROBLEM Consider the double integrator problem in Section 3.11 of [4]. Let x and v define the position and velocity of a particle in rectilinear motion. The governing state equations are
234
ROBERTR. BLESSAND DEWEYH. HODGES
~)
-
v
(104)
-
u
(105)
where the control u is the force, normalized by the mass, imparted to the particle. The problem is to minimize
1/01
u2dt
(106)
However, complications occur when a state inequality constraint is added. Let S - x - L < 0 where L is a constant. The first and second total time derivative of S yield
and S-~) -u Thus, this is a second-order constraint. For certain choices of L, the solution exhibits a touch-point behavior. For other, lower L, the solution rides a boundary arc. First, the unconstrained problem is solved. The state and time boundary conditions are defined to be
x(to)
v(to)- 1 -
X(tl)
(107)
-0
v(tl) + 1 tl-1 Exact and numerical results both show that the state constraint S < 0 is not active for L >_ .25. (The estimated unconstrained value of the limit L is usually a function of the number of elements; this problem is an exception.) For values of L that yield a touch-point solution, the boundary conditions are defined to be
x(t0)
v(to)- 1 x(t+)-L x(t+)-x(tl)
x(t )
v(t2) + I t2- I
-
o
(lo8)
A FINITE ELEMENT METHOD
235
For the constrained boundary-arc case, the boundary conditions are defined to be
x(to)
v(to)- 1 x(t+l ) - L
v(t+l )
x(t+)-x(tl) (t+l ) x(t+ ) - x(t v(t+ ) x(t )
-0
(109)
) )
v(t3) + 1
t3-1 where the tangency conditions are enforced at the unknown entry time, t l, and continuity conditions have been specified at t l and the unknown exit time, t2. In both the unconstrained and constrained cases, the costate boundary conditions are given by Eqs. (15) - (16). Note that the analyst must decide what the switching structure is based on the value of L chosen. Figure 23 shows the displacement of the particle, and Figure 24 shows the control history. In the figures, the exact unconstrained solution is graphed versus a 4-element solution. Also, a touch-point case (with L = 0.2) is graphed versus a finite element run with 4 elements in each phase, denoted with FE(4:4). Finally, a value of L - 0.1 was chosen to demonstrate the exact solution versus the finite element solution (4 elements per phase) for a boundary-arc case. It is seen for this simple example that there is excellent agreement between all results for both the displacement and the control. Of greater interest in this problem is the accuracy with which the finite element method determines the switching structures (i.e., whether the solution is unconstrained, or has a touch-point or boundary-arc solution). For this particular example, the finite element approach identifies the exact unconstrained value of L = 0.25 for the displacement. However, the estimated value of L for which there is a touch-point solution or a boundary-arc solution is a function of the number of elements. The numerical values displayed in Table III below were generated by assuming a touch-point solution. Table III shows the smallest value that L can take and still yield satisfactory results. For instance, when running just 2 elements per phase (the first row of the table), L could be set as low as 0.125 and still yield an approximate answer that satisfies all the discretized necessary conditions. If, however, L - 0.124, then the discretized answer violates the constraint limit before and after (due to the symmetry of the problem) the touch point, which always remains at t = 0.5. As the number of elements is increased, the
236
ROBERT R. BLESS AND DEWEY H. HODGES
0.25
-
0.2-
c" (~
E (:D
0.15-
~J
m
c~
a
l i r a
0.1-
0.05 -
0.0
'
0.0
I
0.2
'
I
'
0.4
I
0.6
Time Figure 23: Displacement versus time
'
I
0.8
~
T
1.0
C~
o
~.~°
o
c-h
o
o0
C~
0~
~.~o
CD
3
0
•
u
b0-
0
0
0
0
0 0
q
•
..s
6"
v
...O --" "" ""
o 0
C ~
I
v
~
O~
o 0
C ~
"1"I rT] m ×
....
Ko Or) ~
0
r" II
m m m m ~ m ~ m ×
. . . . . . . .
". Ko
0
".
0
rII
rII
I i
I
I
,
r" II
0
I
o
..,.. " - " "
.6,
•- - , O
U"
"
i
"~
sS
n
,S
"""~,,
Control
oo
nji o oo
o°
N
•
,
o o"
"",u.
oZ
r"
.n,
"'1, oo
,',II
.............-- "" g ""
° So°
,%
'. n
o I
I
qli
!
~i,
I
,I,
I
- -~
..,..,.- .-- .,,., -,-
-,,..,.....
,
238
ROBERT R. BLESS AND DEWEY H. HODGES
1 allowable values for L increases towards the exact value of ~.
T a b l e I I I . A l l o w a b l e v a l u e s of c o n s t r a i n t limit for t o u c h - p o i n t s o l u t i o n v e r s u s n u m b e r of e l e m e n t s Elements Minimum Value of L for Touch-Point Solution Per Phase 2:2 0.125 (0.124") 4:4 0.146 (0.145") 8:8 0.157 (0.156") 16:16 0.162 (0.161") 32:32 0.165 (0.164") Exact 0.16666666 * denotes the value of L where the constraint is violated Table IV was generated under the assumption that a boundary-arc case existed. This table displays the maximum value of L versus the number of elements where a boundary-arc solution exists. For example, with 2 elements per phase, we could set L = 0.1875. At this value of L, the entry and exit times are almost equal, approximately 0.5 for each. At L = 0.1876, the code produces a solution such that the exit time (t2) is less than the entry time (tl), a nonsense solution. T a b l e I V . A l l o w a b l e v a l u e s of c o n s t r a i n t l i m i t for b o u n d a r y - a r c s o l u t i o n v e r s u s n u m b e r of e l e m e n t s Maximum Value of L for Elements Boundary-arc Solution Per Phase 2:2:2 0.1875 (0.1876") 4:4:4 0.1718 (0.1719") 8:8:8 0.1679 (0.1680") 16:16:16 0.1669 (0.1670") 32:32:32 0.1667 (0.1668") Exact 0.16666666 * denotes the value of L where t2 < t l D
By comparing Tables III and IV, it appears that fewer elements per phase are required to approximate the limiting value between a touchpoint solution and a boundary-arc solution (L = 1/6) when assuming the boundary-arc solution. It is not known if this is true in general.
A FINITE E L E M E N T M E T H O D
239
C. EXAMPLE" B R A C H I S T O C H R O N E P R O B L E M We next consider the brachistochrone problem described in Section 3.11 of [4]. Let x and y define the horizontal and vertical (positive downward) position of a particle, respectively. The governing state equations are
=
2 x / ~ cos u
(110)
-
2 V ~ sin u
(111)
where g is the acceleration due to gravity (= 9.81) and the control u is the angle between the tangent to the path and the horizontal. The problem is to minimize the time it takes for the particle to move from the origin to any point on the line x - t~. A first-order state inequality constraint of the form S = y - x tan 0 - h < 0 is also imposed where 0 and h are constants. The first total time derivative of S yields _ 2v/~ysin(u - O) cosO = 0
(112)
or u = 0 along the constraint boundary. Since this is a first-order constraint, it is reasonable to assume that the constrained solution contains a b o u n d a r y arc. The boundary conditions on the states for this problem may be placed in ~ and the costate boundary conditions are then found from Eq. (15) (16). For the unconstrained case, we have
¢-
{
y(t0)-y0 X(tl)
--
}
--0
(113)
and for the constrained case, we have
x(to) y(to)-
yo
y(t +) - x(t +) tan x(t+)-x(tl) y(t+l ) - y ( t l ) • - x(t ) y(t+ ) - y(t )
0- h -0
(114)
x(t3) where the tangency condition is enforced at the unknown entry time. Notice t h a t the initial value for y has been left as a variable. In the original statement of the problem y0 - 0 for which the Jacobian becomes
240
ROBERTR. BLESSAND DEWEYH. HODGES
singular. Also, the exact solution yields an initial value of the y-costate of - c ~ so that the initial control angle is 90 °. As will be shown, the large and rapid changes in the y-costate value decrease the accuracy of the finite element solution. Figures 25 and 26 show the results for the unconstrained and constrained cases with Y0 - 0.001. The curves are the shooting (essentially exact) results and the discrete points are the finite element values. The constrained case was run with h = 1, t~ = 10, and t a n 0 = 1/2. For the unconstrained case, the finite elements were run with 8 elements. The points are a little off. The reason is clearly seen in Figure 26 (the control history) where it is seen that the final time has been poorly approximated. Also, the initial control angle is not very accurate either. Both of these problems occur because the small initial y0 leads to an extremely rapidly varying costate at the initial time. This is shown in Figure 27 where the unconstrained ycostate history is shown for y0 = 0.001 and Y0 = 0.1. A finer discretization is required to accurately model this costate well. Although the general shape of the state trajectory is found, it is the entry, exit, and final times that are poorly approximated, as is shown in the next two figures. We note that it may be possible to reformulate the problem so that the costate does not tend to - 0 o . With the shooting results taken to be exact, Figures 28 and 29 give relative error of the approximated answers. The relative errors for the entry, exit, and final times, as well as the initial control angle, versus the total number of elements (sum of the elements in each phase) are shown. The slopes of these curves indicate how the relative error varies in proportion to the number of elements. Figure 28 shows the relative errors for Y0 = 0.001. Given that all four curves are nearly parallel, and that the slope continually improves, we can calculate the slopes of the last three segments (the best ones) to be about -1.1, -1.4, and -1.6, respectively. This indicates that the relative error varies proportionally with the inverse of the total number of elements raised to the 1.1, 1.4, or 1.6 powers, and that the error is heading toward second-order accuracy. This accuracy is less than what is typically seen for the finite element method, and is caused by the large gradient in the y-costate value at the initial time. This large gradient in the costate changes as the initial value for y is made larger, as was seen in Figure 27. For y0 = 0.1, the relative errors are plotted in Figure 29 and are seen to be an order of magnitude better. The slopes of the curves are more constant after the initial "bend" in the curve, and the slopes of the last three segments are -1.9, - 2 , and - 2 indicating that the relative error varies approximately inversely with the square of the number of elements. This "second-order" nature of the present finite element method appears to be typical for most problems.
A FINITE ELEMENT METHOD
241
X 0 0 i
1 I
2 I
3 I
4 I
5 I
6 I
7 I
8 I
9 I
10 .I
0
1
,~,,, I
2
• ""o,,,,,,
•
, %,%%,
l
>4
~
1
6
~
7
~
Figure 25: Trajectory versus time
Unconstrained; Shooting Unconstrained; FE(8) h/l=O. 1; Shooting h/l=O. 1; FE(4:4:4)
242
ROBERT R. BLESS A N D D E W E Y H. H O D G E S
90
80
• i
70-
o',,
60-
'\
'\
h/l=0.1;
FE(4:4:4)
• i
50\
m
O L_
O
FE(8)
•
o',
"~ E O
Unconstrained;
• \
-O
Shooting
h/l=O. 1; S h o o t i n g
\
o\
v
Unconstrained;
-
\ i
i
i i
40-
i
i i
_
\
• \ i
30-
\ t
i
\
• ~--e-,
20-
•
O"',
10-
i 0.0
i
i
i
I 0.5
i
,
,
,
I
'
'
1.0
Time (s) Figure 26" Control history versus time
'
'
I 1.5
'
" "
'
I 2.0
A FINITE ELEMENT METHOD
243
-2
-3 y(O) = 0.001; Shooting 0 0 >,,
•
-4
y(O) = 0.001; FE(8) y(O) = O. 1; Shooting
•
-5
y(O) = 0.1; FE(8)
-6
-7
-8
I
0.0
~
w
I
,
I
0.5
~
1
I
~
I
~
~
1.0
Time (s) Figure 27: Costate history versus time
~
'
I
1.5
'
I
2.0
244
ROBERT R. BLESS AND DEWEY H. HODGES
10 0 -
t.__
1 0 "1
0 t._
:: ::::::::i
t._
LU
i
>
,.m
•~
"....
~
m
n-'
1 0 -2
•
.........
• •
"~
"... "\
,, °°°iooo ~,"
"".,,.~.............
Error in Entry T i m e
.................. Error in Exit Time
'.,,~ .......
".,
Error in Final T i m e
10 .3
I
I
I
w
I
I
I
I
i
i
I
10 ~
Number of Elements Figure 28" Error for Yo - 0 . 0 0 1 versus time
i
i
I
I
I
I
10 ~
A FINITE ELEMENT METHOD
245
1 0 o-
,%. '%,,~
] 0 "1 ~ '"'"'"'.. ~.~,
t.__
O
"" """ %'" '~'%'",, ""'......
•~.
•~.
I=.
ILl >
%, %',%
"'....
,~
....
"%.. ~.
%. "%
•~.
1 0 "~
•
'%
".... "~.
, . m
,%
m
~ •~
"'..... "'.... •~ .... "~ "..
rc
%%,
,%. "'... \. "'-.. \
\
"... "'.. "'"'".. \.\. "'".....
u'~''~''~
~. "'% \ "..
1 0 -3 - -
"N.\"'.......
E r r o r in Initial C o n t r o l
",,
"\'N "'"'"'.
•~ "%.
E r r o r in E n t r y T i m e
~
.................. E r r o r in Exit T i m e E r r o r in Final T i m e 10'
,
l
J
,
,
, i I
,
,
I
Number of Elements 29:
Error
f o r yo = 0 . 1 v e r s u s
I
I I , I 10 2
101
Figure
I
"... "\.\. %...... \ "% "\.\ "'.o.., "\
time
246
ROBERT R. BLESS AND DEWEY H. HODGES
D. EXAMPLE: MULTI-PHASE ASCENT TO ORBIT We now turn to a more nonlinear problem to evaluate the usefulness of the finite element method in the presence of state constraints. The same launch vehicle model as used in Section II, E is used, along with a dynamic pressure constraint of the form S ( h , V ) - q - qt~m <_ 0
(115)
This is a first-order state constraint as is seen by taking the first total time derivative and obtaining Op hV 2 -(7=
Oh
2
~- p V ~ /
(116)
By substituting in values for p and h, and simplifying the resulting expression, we obtain the following control-dependent algebraic constraint that must be enforced during the time interval that q = qlim T cos a - D ( a )
p sin 3'
r2
V 2 sin 3' =0 2(8600)
(117)
This constraint will be enforced as one boundary arc occurring before the staging time. Thus, the first phase of the unconstrained problem will be broken into three phases consisting of an unconstrained, constrained, and unconstrained arc. This is followed by the final phase from staging to the final time, for a total of four phases. The boundary condition vector must also be modified to include the tangency conditions and extra internal boundary conditions. The vector is
A FINITE ELEMENT METHOD
247
re(to) - ,no h(to) - 60
V(to)-25 7 ( t 0 ) - 1.5
q(tl+) - q . ~
m(t~+) - m ( t ~ )
) -h(tl) ) - V(tl) ) -7(tl) m(t+~ ) - m(t~) h(t+~ ) - h(t~) V(t+~ ) - V(t~) ~(t~+) - 7(t~) h(t+l V(t+~ ~(t+l
~b
=0
(118)
m ( t +) - r e ( t 3 ) + 29920 h(t +) - h(t3 )
v(t+~ ) - v(t~) .y(t+~ ) - ~(t~) h ( t 4 ) - hy V ( t 4 ) - Vf ~'(t4) t 3 - 195 Some of the results are shown in Figures 30 - 37. The maximum value of the dynamic pressure q is approximately 55.3 kPa. Shooting results and finite element results were obtained for cases where the maximum allowable value for the dynamic pressure is 45 kPa, 35 kPa, and 25 kPa. The dynamic pressure histories are shown in Figure 30. However, for clarity in the figures, results will only be presented for the 25 kPa case. In all of the figures, the unconstrained shooting results are shown, for comparison purposes. The finite element results are for 4 elements in each of the first three phases and 8 elements in the last phase (the one after staging) and are denoted with "FE(4:4:4:8)". Figure 31 shows the normalized altitude profile versus the normalized velocity. As expected, the vehicle climbs higher, more quickly, as the constraint becomes stricter. Note also that even with the coarse discretization presented, the finite element results lie essentially on top of the exact curves. Figure 32 shows the flight-path angle 7 versus time. The vehicle is seen to roll over more slowly with the strict q limit in order to climb out of the atmosphere more quickly. Recall from Figure 13 that the finite element results show some error for the unconstrained results. However, with the addition of finite elements and the less rapidly varying flight-path angle history, the answers are much more accurate for the constrained case.
248
ROBERT R. BLESS AND DEWEY H. HODGES
Figures 3 3 - 36 give four costate histories for the vehicle. Due to the internal boundary condition on q, a jump occurs in the altitude costate and the velocity costate. Note also that the finite elements have solved for this jump and the initial costate values accurately. These values are, of course, critical to the shooting code. The angle of attack history is shown in Figure 37. The control experiences large variations due to the changing constraint limits. Once again the finite element results are quite accurate, although less so at the initial time where the angle of attack changes rapidly. As mentioned at the beginning of this section, an interesting touchpoint solution exists for this problem. The control is positive along the constrained arc for certain values of the dynamic pressure limit, and is negative along the constrained arc for other values. At a particular value of qtim, the control will be zero, and a touch-point solution occurs. A first-order touch-point solution is possible because, from Eq. (117)
0b
0v
0--~ = 0 a = 0
at
a - 0
(119)
Figure 38 shows the length of the constrained arc and final mass loss as a function of dynamic pressure limit. It is interesting to note that the length of the constrained arc is extremely small (less than 0.5 seconds) until after the limit is set below the touch-point limit, which occurs just under the 35 kPa arc, at which point the length of the constrained arc grows rapidly as the limit is decreased. Also of interest is the fact that the loss of final mass becomes increasingly steeper with a stricter dynamic pressure limit. The combination of rapidly obtained finite element solutions and accurate shooting solutions allowed for the discovery of this anomalous first-order state-constraint behavior.
C~
=
CD C~
t~
(1) oo 0° o
O~
m v
3
--I
0 0
0 0
0
~ o
o
o
0
0
L
I
0
I
I
I
I~ 0
I I I I t
0 ~'.
o
~.
oo
o 0
I
II 0"I ~-D
II 07 ~-D
o1 ~-D
II
x .o
0 0 =. ~
0 0 ~
=.
OOO0007"I
x .Q
x .0
0 0 =. ~
"'"
66
"~ ..
~
O1 ~. "D
II
x c~
•
I
.~ 0
""--~)
..~.. ~ ~ -
--..:........ " " ..--. .-...,,,,...
Z.-.--._-.ZT-.
................................
I
..... •.
...................T..." .... .....
, ' , 3 3 3 3
I I I i
I
0.~ 0
I
Dynamic Pressure (kPa) I
01 0
. ..,...,.,.... ~ "
-'.--.~..~
I
I
0"~ 0
250
ROBERT R. BLESS AND DEWEY H. HODGES
1.0~
A
A
_
..~ 0.8
_~~o,6 N 0.4
Z
0.0
:4:4:8) max q = 25kPa; Shooting oting
/,,,,"
0.2
•
0.0
'
I
0.2
'
I
0.4
'
I
0.6
'
I
0.8
NormalizedVelocity
Figure 31" Normalized altitude versus normalized velocity
~
I
1.0
A FINITE ELEMENT METHOD
251
80
70
¢D ,<: _
5o
i "
40
",,,,, \
•
max q = 25kPa; FE(4:4:4:8)
-
max q = 25kPa; Shooting
.&, r--
.~
2o-
IJ _
-
100-
-10
'
0
I
100
'
I
'
200
I
300
Time (s) Figure
32" F l i g h t - p a t h
angle versus
time
'
I
400
'
I
500
252
ROBERT R. BLESS AND DEWEY H. HODGES
o.o-0.1
"'"" ..... .
-0.2 -0.3 -O.4
8
0
-o.5
~
-0.6 -0.7 -0.8 -0.9 -1.0
I 0
100
200
300
Time (s) Figure 33" Mass costate versus time
400
500
A FINITE ELEMENT METHOD
253
0.0 zooOO° ° ° ° °
I
-0.05
oo
o
(D
-0.1
0
0
(D "0
m
m
< m
-0.15
max q = 25kPa; FE(4:4:4:8)
-0.2
max q = 25kPa; Shooting q unconstrained; Shooting
-0.25
-0.3
I
0
'
I
100
'
I
200
'
I
300
Time (s) Figure 34: Altitude costate versus time
'
I
I
400
5OO
254
ROBERT R. BLESS AND DEWEY H. HODGES
00-] -0.1
-0.2 .b..,0
-0.3
0
0
>.,
i
m
-0.4
m
0
0 (D
>
-0.5
max q = 25kPa; FE(4:4:4:8) max q = 25kPa; Shooting q unconstrained; Shooting
-0.6
-0.7 I
-0.8
'
0
I
100
'
I
200
'
I
300
Time (s) Figure 35: Velocity costate versus time
'
I
i
400
5OO
A FINITE ELEMENT METHOD
255
"10 ~
0.2
max q - 25kPa; FE(4:4:4:8)
0.04
max q = 25kPa; Shooting q unconstrained; Shooting
i i
\
-0.2 -
i
0O
0
0
i i
\
-0.4-
\
i
i
(D
133
/
/
#
-0.6 /
<
t-
/
-0.8-
II
/
\, \
u
',,,,
-1.0-
c-"
m
m
O3 m
U..
/
",,
-1.2-1.4-
~,,
ooo1°
,,/
•
•
-1.6'
0
I
100
'
I
2OO
'
I
300
Time (s) Figure 36" Flight-path-angle costate versus time
'
U
400
'
I
500
256
ROBERT R. BLESS AND DEWEY H. HODGES
3025-
max q - 25kPa; FE(4:4:4:8) m a x q - 2 5 k P a ; Shooting
20
(D "O
t ! t
i
,
15-
;i
q unconstrained; \
",\ i
i
i
\ i
v
rD
Shooting
10-
i
i
i II i i
< O (D
%%%
O
_
m~=..=
¢-
<
-5-10-15-20
' 0
I 100
I 200
'
I 300
Time (s) Figure 37" Angle of attack versus time
'
I 400
'
I 500
0
0
r~ r~
0
r/l
.°
V
m .
3
m .
r-
"u
0
m l
3
01
O
0
O
I
t
O Q
i
I
~
0
t 0 0 Q
I
'
i
I
T 01 0 Q
(/)
E
..,..
"1"1 --,.,,.
0
0
,.,..
¢.. (-)
:3:3
i
Mass Loss (kg)
::)>
p-..q
0 :::I U)
0
:3"
i
I
0 Q O
i
Length of Constrained Arc (s) I
,
0
01
I
258
ROBERT R. BLESS AND DEWEY H. HODGES
V. SUMMARY This chapter was concerned with a method for the solution of optimal control problems. The method, called the weak principle for optimal control, is based on time-domain finite elements. The goals of this work, as stated in Section I, B, were to develop an efficient and accurate algorithm that would (1) rapidly obtain an approximate solution for possible real-time implementation; (2) identify the optimal switching structure for constrained problems; and (3) produce initial guesses for a shooting method. With respect to these three goals, it has been observed that: 1. It is helpful to obtain a converged answer ments and then add elements using guesses converged solution. As was demonstrated II, E (see Table I), this proved to be a fast ing accurate answers which could then be
for a small number of eleobtained from a previously in the example in Section and reliable way of obtainused in a shooting code.
2. After the dynamics of the system have been properly modeled with some minimum number of elements, then the finite element method possesses second-order accuracy. It was shown in Section IV, C that there is a degradation in accuracy if the variables change more rapidly than the finite element discretization is able to model. 3. As was shown in Tables III and IV of Section IV, B, this formulation can be useful in determining the switching structure for an optimal trajectory. 4. The speed of solution is directly related to the speed of solving a sparse linear system. The sparsity of the Jacobian was highlighted in Table II of Section III, C. With the use of an intelligent sparse solver, this finite element method can be very efficient in terms of run-time and memory requirements. There is now a patent pending on a real-time guidance algorithm based on the present methodology
[53].
Some of the characteristics and features of the weak principle are summarized below. 1. The necessary conditions of optimality are discretized, and all strong boundary conditions are transformed into weak ones. The weak principle finds candidate extremal solutions, i.e., ones that satisfy all the discretized necessary conditions. 2. Because all strong boundary conditions are cast in the form of natural or weak boundary conditions, then the same test and shape functions may be chosen for every optimal control problem.
A FINITE ELEMENT METHOD
259
3. The choice of shape functions allowed by the weak principle permits all integration to be done by inspection, regardless of the degree of nonlinearity in the problem. Thus, no errors are introduced as the result of numerical quadrature. 4. One advantage of the integration being done by inspection is that algebraic equations may be derived prior to specifying the problem to be solved. Because the form of the algebraic equations is known, it was possible to create a general algorithm, VTOTS, for the solution of optimal control problems. 5. The derivations in this chapter lay the groundwork for higher-order test and shape functions to be used. Although this introduces the need for element quadrature, greater accuracy is obtained with fewer elements [54].
VI. A C K N O W L E D G E M E N T S Technical discussions with Drs. Marco Borri, Anthony J. Calise, J. Eric Corban, Daniel D. Moerder, David A. Peters, J. V. R. Prasad, and Hans Seywald are gratefully acknowledged. The authors also thank Michael S. Warner for editorial comments and a careful reading of the manuscript. This work was partially supported by grants from NASA Langley Research Center, the technical monitor of which was Daniel D. Moerder.
References [1] P. E. Gill, W. Murray, and M. H. Wright, Practical Optimization, Academic Press, New York (1981). [2] J. Betts and W. Huffman, "The Application of Sparse Nonlinear Programming to Trajectory Optimization," Proceedings of the AIAA Guidance, Navigation, and Control Conference, Portland, Oregon, AIAA Paper 90-3448, pp. 1198-1218 (1990). [3] J. E. Rader and D. J. Hull, "Computation of Optimal Aircraft Trajectories Using Parameter Optimization Methods," Journal of Aircraft 12(11), pp. 864-866 (1975). [4] A. E. Bryson, Jr. and Y.-C. Ho, Applied Optimal Control, Hemisphere Publishing Corporation, New York (1975).
260
ROBERT R. BLESS AND DEWEY H. HODGES
[5] E. B. Lee and L. Markus, Foundations of Optimal Control Theory, John Wiley & Sons, Inc., New York (1967).
[6]
L.D. Berkovitz, Optimal Control Theory, Springer-Verlag, New York (1974).
[7]
C. R. Hargraves and S. W. Paris, "Direct Trajectory Optimization Using Nonlinear Programming and Collocation," Journal of Guidance 10(4), pp. 338-342 (1987).
[8] C. Hargraves, F. Johnson, S. Paris, and I. Rettie, "Numerical Com-
putation of Optimal Atmospheric Trajectories," Journal of Guidance and Control 4(4), pp. 406-414 (1981).
[9] J. Vlassenbroeck, "A Chebychev Polynomial Method for Optimal Control with State Constraints," Automatica 24(4), pp. 499-506 (1988). [10] C. T. Kelley and E. W. Sachs, "Quasi-Newton Methods and Unconstrained Optimal Control Problems," SIAM Journal on Control and Optimization 25(6), pp. 1503-1516 (1987). [11] W. A. Gruver and E. Sachs, Algorithmic Methods in Optimal Control, Pitman Publishing, Boston, Massachusetts (1981).
[12] P. E. Gill, W. Murray, M. A. Saunders, and M. H. Wright, "User's
Guide for NPSOL: A FORTRAN Package for Nonlinear Programming," Technical Report SOL 86-2, Department of Operations Research, Stanford University (1986).
[13] S. W. Paris and C. R. Hargraves, "Optimal Trajectories by Implicit Simulation: Volume I-Formulation Manual," AFWAL-TR-883057 (1988).
[14] G. L. Brauer, et al., "Final Report: Program To Optimize Simulated Trajectories (POST): Volume I-Formulation Manual," NASA CR132689 (1975). [15] H. J. Oberle and W. Grimm, "BNDSCO: A Program for the Numerical Solution of Optimal Control Problems," English Translation of DLRMitt. 85-05, Oberpfaffenhofen, Germany. [16] R. R. Bless, et al., "Variational Trajectory Optimization Tool Set; Technical Description and User's Manual," NASA TM 4442 (1993). [17] S. M. Roberts and J. S. Shipman, "Multipoint Solution of Two-Point Boundary-Value Problems," Journal of Optimization Theory and Applications 7(4), pp. 301-318 (1971).
A FINITEELEMENTMETHOD
261
[18] 5. M. Roberts and J. S. Shipman, Two-Point Boundary Value Problems: Shooting Methods, American Elsevier Publishing Company, Inc., New York, New York (1972). [19] P. K. A. Menon and L. L. Lehman, "A Parallel Quasi-Linearization Algorithm for Air Vehicle Trajectory Optimization," Journal of Guidance 9(1), pp. 119-121 (1986). [20] M. B. Subrahmanyam, "A Computational Method for the Solution of Time-Optimal Control Problems by Newton's Method," International Journal of Control 44(5), pp. 1233-1243 (1986). [21] D. E. Kirk, Optimal Control Theory: An Introduction, Prentice-Hall Inc., Englewood Cliffs, New Jersey (1970). [22] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, Cambridge (1986). [23] H. G. Visser, H. J. Kelley, and E. M. Cliff, "Energy Management of Three-Dimensional Minimum-Time Intercept," Journal of Guidance, Control and Dynamics 10(6), pp. 574-580 (1987). [24] H. Seywald, "Optimal Control Problems with Switching Points," NASA c a 4393 (1991). [25] I. Fried, "Finite Element Analysis of Time Dependent Phenomena," AIAA Journal 7, pp. 1170-1173 (1969). [26] J. H. Argyris and D. W. Scharpf, "Finite Elements in Time and Space," Aeronautical Journal of the Royal Society 73, pp. 1041-1044 (1969). [27] C. D. Bailey, "Application of Hamilton's Law of Varying Action," AIAA Journal 13(9), pp. 1154-1157 (1975). [28] T. E. Simkins, "Unconstrained Variational Statements for Initial and Boundary Value Problems," AIAA Journal 16, pp. 559-563 (1978). [29] D. L. Hitzl and D. A. Levinson, "Application of Hamilton's Law of Varying Action to the Restricted Three-Body Problem," Celestial Mechanics 22, pp. 255-266 (1980). [30] M. Baruch and R. Rift, "Hamilton's Principle, Hamilton's Law, 6n Correct Formulations," AIAA Journal 20, pp. 687-692 (1982). [31] M. Borri, et al., "Dynamic Response of Mechanical Systems by a Weak Hamiltonian Formulation," Computers and Structures, 20(1-3), pp. 495-508 (1985).
262
ROBERTR. BLESSANDDEWEYH. HODGES
[32] D. A. Peters and A. Izadpanah, "hp-Version Finite Elements for the Space-Time Domain," Computational Mechanics 3, pp. 73-88 (1988). [33] W. E. Bosarge. Jr. and O. G. Johnson, "Error Bounds of High Order Accuracy for the State Regulator Problem Via Piecewise Polynomial Approximations," SIAM Journal on Control 9(1), pp. 15-28 (1971). [34] P. M. Prenter, Splines and Variational Methods, John Wiley & Sons, Inc., New York, New York (1975). [35] L. J. Segerlind, "Weighted Residual Solutions in the Time-Domain," International Journal for Numerical Methods in Engineering 28, pp. 679-685 (1989). [36] W. N. Patten, "Near Optimal Feedback Control for Nonlinear Aerodynamic Systems with an Application to the High-Angle-of-Attack Wing Rock Problem," AIAA Paper 88-4052-CP (1988). [37] F. G. Mathis and G. W. Reddien, "Ritz-Trefftz Approximations in Optimal Control," SIAM Journal on Control and Optimization 17(2), pp. 307-310 (1979). [38] B. Asselmeyer, "A Ritz-Type Optimization Method for Optimal Control Problems and Its Application to Hierarchical Final-Value Controllers," Control and Dynamic Systems: Advances in Theory and Application 23, pp. 295-318 (1986). [39] H. Oz and A. Raffle, "Inverse Response Problem (Control) of Dynamic Systems Via Hamilton's Law," Computer Methods in Applied Mechanics and Engineering 62, pp. 17-26 (1987). [40] D. H. Hodges and R. R. Bless, "Weak Hamiltonian Finite Element Method for Optimal Control Problems," Journal of Guidance, Control, and Dynamics 14(1), pp. 148-156 (1991). [41] D. H. Hodges, R. R. Bless, A. J. Calise, and M. Leung, "Finite Element Method for Optimal Guidance of an Advanced Launch Vehicle," Journal of Guidance, Control, and Dynamics 15(3), pp. 664-671 (1992). [42] R. R. Bless and D. H. Hodges, "Finite Element Solution of Optimal Control Problems with State-Control Inequality Constraints," Journal of Guidance, Control, and Dynamics 15(4), pp. 1029-1032 (1992). [43] M. Borri, "Helicopter Rotor Dynamics by Finite Element Time Approximation," Comp. and Mathematics with Applications 12A(1), pp. 149-160 (1986).
A FINITE ELEMENT METHOD
263
[44] M. Borri, F. Mello, and S. N. Atluri, "Primal and Mixed Forms of Hamilton's Principle for Constrained Rigid Body Systems: Numerical Studies," Computational Mechanics 7, pp. 205-220 (1991). [45] I. M. Gelfand and S. V. Fomin, Calculus of Variations, Prentice-Hall, Inc., Englewood Cliffs, New Jersey (1963).
[46]
E. B. Becker, G. F. Carey, and J. T. Oden, Finite Elements: An Introduction, Volume 1, Prentice-Hall, Inc., Englewood Cliffs, New Jersey (1981).
[47]
F. A. Valentine, "The Problem of Lagrange with Differential Inequalities as Added Side Conditions," Contributions to the Calculus of Variations, Chicago, IL, Chicago University Press, pp. 407-448 (1937).
[48]
I. S. Duff, Harwell Subroutine Library, Chapter M, Computer Science and Systems Division, Harwell Laboratory, Oxfordshire, England (1988).
[49]
D. H. Jacobson, M. M. Lele, and J. L. Speyer, "New Necessary Conditions of Optimality for Control Problems with State-Variable Inequality Constraints," Journal of Mathematical Analysis and Applications 35, pp. 255-284 (1971).
[50]
H. Seywald, "On the Existence of Touch Points for First-order State Inequality Constraints," Proceedings of the AIAA Guidance, Navigation, and Control Conference, Monterey, California, AIAA Paper 93-3743CP, pp. 372-376 (1993).
[51]
E. M. Cliff, H. Seywald, and R. R. Bless, "Hodograph Analysis in Aircraft Trajectory Optimization," Proceedings of the AIAA Guidance, Navigation, and Control Conference, Monterey, California, AIAA Paper 93-3742-CP, pp. 363-371 (1993).
[52]
A. E. Bryson, Jr., W. F. Denham, and S. E. Dreyfus, "Optimal Programming Problems with Inequality Constraints, I: Necessary Conditions for Extremal Solutions," AIAA Journal 1(11), pp. 2544-2550 (1963).
[53]
G. T. Haskins, M. G. Johnson, and D. H. Hodges, "Real Time Missile Guidance System," United States Patent Pending (1993).
[54] D. H. Hodges and L. J. Hou, "Shape Functions for Mixed p-Version Finite Elements in the Time Domain," J. Sound and Vibration 145(2), pp. 169-178 (1991).