Home

as a PDF

Contents

1. 0 1 0 2 0 3 0 4 0 5 0 6 0 7 time Figure B9 Constrained double integrator Hamiltonian 41 0 8 0 9 H uu EigenValues vs time 2 T T T T T T T s gt 5 m 3 T 0 l l l l l l l l l 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 time Figure B10 Constrained double integrator eigenvalues of Huu Random numbers will work for initial estimates The results for this case are not included herein How does the user know whether the trajectory touches or rides the constraint Always run the unconstrained problem first to see whether the constraint limits are violated First order constraints always have a constrained arc whereas second order constraints frequently have touch point solutions and constrained arcs Finally if a touch point solution is assumed and the actual solution rides the constraint then somewhere there 1s a constraint violation Control Constrained Problem This example is taken from section 3 8 of reference 7 The problem is to minimize loce cd ire J T dt sry 5 v where T 10 z and u are scalars and the initial condition is 0 19 945596 The state equation is amp h t u with A ft 1 t at Two control inequality constraints are imposed to enforce u 1 g u 1 lt 0 ga u 1 lt 0 42 The following mac file produces the needed plant mex4 file This is the fixed time control const
2. Terms for dynamic equations alt h hf rho rhoO exp alt betar faxial ca 0 b rho v vf 2 sref fnorm cna 0 b rho v vf 2 sref alpha 180 Api drag fnorm sin alpha tfaxial cos alpha xlift fnorm cos alpha faxial sin alpha p pO exp alt betap thri tvac engi engi ae p thr2 tvac eng2 eng2 ae p DYNAMICS EQUATIONS PHASE 1 rmassdoti tvac engi rmassO grav spimpi hdoti v sin gam vf hf vdoti thri cos alpha drag rmass rmassO vf earmu sin gam vf retalt 2 46 gamdoti thri sin alpha txlift rmass v rmassO vf t COrf v reralt earmu v vf retalt 2 cos gam rmassdoti i00 rmassdoti eta 2 hdoti 100 hdoti eta 2 vdoti 100 vdoti eta 2 gamdoti i100 gamdoti eta 2 DYNAMICS EQUATIONS PHASE 2 rmassdot2 tvac eng2 rmassO grav spimp2 hdot2 v sin gam vf hf vdot2 thr2 cos alpha drag rmass rmassO vf earmu sin gam vf retalt 2 gamdot2 thr2 sin alpha xlift rmass v rmassO vf vftv retalt earmu v vf retalt 2 cos gam rmassdot2 100 rmassdot2 eta 2 hdot2 100 hdot2 eta 2 vdot2 100 vdot2 eta 2 gamdot2 100 gamdot2 eta 2 DYNAMICS EQUATIONS LIST delist rmassdoti hdoti vdoti gamdoti1 0 rmassdot2 hdot2 vdot2 gamdot2 0 Several features are important in this mac file 1 Notice that rmass not mass is used for the state because the FORTRAN files treat mass as an integer 2 Note that the states in psilist are
3. few of the commercially available programs for solving optimal control problems are mentioned below The first two programs solve general MPBVP s whereas the last two are particularly designed to optimize flight vehicle trajectories The Chebyshev Trajectory Optimization Program CTOP is useful in several practical applications ref 2 This program solves problems directly and parameterizes the functions using Chebyshev polynomials Penalty functions enforce the equations of motion and path constraints The Nonlinear Programming for Direct Optimization of Trajectories NPDOT package uses piecewise polynomials and collocation to satisfy the differential equations Results presented in reference 1 show that the NPDOT runs more quickly than the CTOP does Both programs are generic optimization programs that are not limited to aerospace problems The Program to Optimize Simulated Trajectories POST targets and optimizes point mass trajectories for a powered or unpowered vehicle that operates near a rotating oblate planet ref 18 The POST allows the solution of a wide range of flight problems that include aircraft performance orbital maneuvers and injection into orbit The user can select the optimization variable the dependent variables and the independent variables from a list of more than 400 program variables The POST also operates on several computer systems Another useful program is Optimal Trajectories by Implicit Simulation OTIS The
4. 0 so that 7 is a constant Its value is equal to the final time as yet unknown In this case to prevent the time scale from becoming negative set n t t t 7 Now dx dx o 7 dyn dt Therefore all the differential equations are multiplied by z2 Similarly the VTOTS can also solve nonautonomous problems In this instance the time t becomes a state with the additional boundary condition that this new state has an initial condition of 0 the corresponding differential equation is f 1 Multiphase problems can be handled by a straightforward extension of this technique Examples of time scaling are given in the Additional Examples section The File Vtotsinfo m In addition to the files namcom nml and plant mex4 created in the MACSYMA environ ment the user must supply a MATLAB file called vtotsinfo m Because the VTOTS uses an iterative method to solve MPBVP s an initial guess is required The file vtotsinfo m stores this initial guess with several other optional variables Some variables are common to both algorithms and some are specific to either the finite element or the shooting algorithm All the variables are discussed below Variables common to the finite element and shooting algorithms The following variables may be defined in vtotsinfo m if not defined they are not used prob name comment about the current problem placed in single quotes timestate integer between 1 and number of states in pro
5. 1 The next variable is the column vector t This variable defines the location on the discretized time line where the estimates are needed Recall that estimates are needed at the endpoints of the phase and at the midpoints of the elements fig 1 Next the initial guess of the solution is put in a column vector with the reshape command Finally estimates for the discrete multipliers and time are added Because psilist has a length of six and tsilist has a length of one seven estimates of 1 for the discrete multipliers are given Also because the problem has been scaled to run from 0 to 1 the estimate for the final time is 1 Finally 25 the variable scale scales the output quantities See the section entitled Variables Common to the Finite Element and Shooting Algorithms Because only the time line is scaled the first number is 20 0 the actual final time and the next 8 values states and costates are 1 0 because these were not scaled Another example of the use of scale is given in the section entitled A Two Stage Rocket Problem Running the MACSYMA commands in figure Al creates the plant mex4 file After the plant files plant mex4 vtotsinfo m and namcom nml are created VTOTS is ready to run After the user enters MAT LAB typing vtots at the MATLAB prompt returns the following Variational Trajectory Optimization Tool SET VTOTS v 2 0 PROBLEM SPECIFIC INFORMATION BHO FIX B amp H Fixed time prob Number of stat
6. 9 J e x 5 3 0 1 Les t t S 2 I 0 l l l 4 l l l 0 5 10 15 20 0 5 10 15 20 time time Figure A4 Costate histories 2 ul vs time zo L l fi 0 5 10 15 20 time Figure A5 Control history 31 Hamiltonian vs time 173 709 T T T T T T T T 86 855 Hamiltonian 0 l l l l l l l l l 0 2 4 6 8 10 12 14 16 18 20 time Figure A6 Hamiltonian history integral cost plus adjoined dynamics measure of global convergence of algorithm 32 genValues H uu Ei H_uu EigenValues vs time 90 l l l l l 0 2 4 6 8 10 12 14 16 18 time Figure A7 Eigenvalues of second partial of Hamiltonian with respect to control second order sufficient condition 33 20 Appendix B Additional Example Files This appendix presents several example problems for use with the VTOTS The Unconstrained Double Integrator As a first example consider the simple double integrator defined by two states x and v with differential equations x v v u and boundary conditions v 0 1 v 1 2 1 x 0 0 x 1 0 The problem is to find the condition u t that minimizes 1 1 J 5 u dt 2 Jo All the information needed to set up the appropriate mac file to produce the plant mex4 file is shown above The mac file is This is the fixed time double integrator problem stlist x v ctlist ul glist L 11 glist L l1l ellist 0 5 u 2 phi 0 thi 0 psilist x 1 1 x0 v 1 1
7. Appl amp Methods vol 10 no 2 Apr June 1989 pp 147 171 Menon P K A and Lehman L L A Parallel Quasi Linearization Algorithm for Air Vehicle Trajectory Optimization J Guid Control amp Dyn vol 9 no 1 Jan Feb 1986 pp 119 121 Roberts S M and Shipman J S Multipoint Solution of Two Point Boundary Value Problems J Optim Theory Appl vol 7 no 4 Apr 1971 pp 301 318 Bless Robert R Time Domain Finite Elements in Optimal Control With Application to Launch Vehicle Guidance NASA CR 4376 1991 Roberts Sanford M and Shipman Jerome S Two Point Boundary Value Problems Shooting Methods American Elsevier Publ Co 1972 Subrahmanyam M B A Computational Method for the Solution of Time Optimal Control Problems by Newton s Method nt J Control vol 44 no 5 Nov 1986 pp 1233 1243 Brauer G L Cornick D E Habeger A R Petersen F M and Stevenson R Program To Optimize Simulated Trajectories POST Volume Formulation Manual NASA CR 132689 1975 Macsyma Reference Manual Version 13 Doc No 5MI0500030 013 Symbolics Inc Nov 1988 MATLAB for Sun Workstations User s Guide Math Works Inc Jan 31 1990 Kane Thomas R and Levinson David A Dynamics Theory and Applications McGraw Hill c 1985 53 REPORT DOCUMENTATION PAGE E Public reporting burden for this collection of infomation is estimated to average 1 hour per response including
8. by expressing terminal values of x and A their values at the end of phases as functions of initial values their values at the beginning of phases This conversion is achieved by integrating the solution of the ordinary differential equations eqs 7 and 8 from the initial values to the terminal conditions For simplicity consider the case with no internal events so that the boundary conditions of the problem are v Xo Xy 0 19 96 Ob T T Ao Oxo u xy 0 20 96 Ow T T Ag Ox v ax 0 21 where xo x 0 Ao A 0 xy x ts Ap A tp The variables Xy and af are evaluated as 1 xir f x d 7 etg 22 0 l H Ap ZA0 72 3 0M de 7 etg 23 where H is the Hamiltonian that is defined in equation 6 and is evaluated along x t A t and u t and is a root of PE x 0 24 which is obtained by numerical solution of Hy 0 in terms of x and A at each instant The result is that appears as u x A in the calculations The partial derivatives and Wy are Up Hoa Hur 25 9 uy Hay Hux 26 where Huu is assumed to have full rank The variable 7 in equations 22 and 23 is a parameter that scales the dummy independent variable teg 0 lt lt 1 27 In the implementation of the VTOTS shooting algorithm 7 is appended to x as an additional state variable with 7 0 28 and is solved with boundary conditions appropriate to the free or fixe
9. dynamics in equation 8 will have positive real parts Direct numerical differentiation of W V o would require perturbation of Vo an action that could excite modes corresponding to unstable eigenvalues This problem is avoided through the use of equations 32 and 33 Although the shooting method is slower than the finite element method the shooting method solution is as numerically accurate as the integrator used to propagate the state and costate equations Concluding Remarks This report provides a technical overview and a brief description of the algorithms that comprise a new software package for solving optimal control problems Although many excellent programs are available for this purpose the Variational Trajectory Optimization Tool Set VTOTS offers some new features 1 The VTOTS provides two algorithms based on indirect methods most available programs are based on direct methods 2 The VTOTS provides a finite element algorithm for approximate solutions and a shooting algorithm for numerically accurate solutions 3 An optimal control problem from any discipline may be solved when properly formatted however this flexibility requires a VTOTS user to supply application specific code The appendixes contain a complete user s manual that includes a detailed example and helpful hints Additional examples even those using very few elements in the finite element algorithm show that a good approximation to a solution is pos
10. getu m jacobi m makepsi m psiend m psist m rhs m salvo m ushape m varstr m calculates the Hamiltonian H and the eigenvalues of 0 H 0u for plotting purposes solves for the optimal control using a Newton iteration calculates an analytical Jacobian matrix needed by rhs m calculates the error vector W used to solve for the initial values with a Newton iteration calculates the OW OX matrix that is part of the Newton step to find the initial values calculates the QW O X matrix that is part of the Newton step to find the initial values calculates the right side of the differential equations integrated by the ode45 m integrator the driver m file for the shooting code all integrations and error calcula tions are done in this file conditions a control guess for getu m saves variables so that fewer globals are needed 52 References 1 10 11 12 13 14 15 16 17 18 19 20 21 Hargraves C R and Paris S W Direct Trajectory Optimization Using Nonlinear Programming and Collocation J Guid Control amp Dyn vol 10 July Aug 1987 pp 338 342 Hargraves Charles Johnson Forrester Paris Stephen and Rettie lan Numerical Computation of Optimal Atmospheric Trajectories J Guid amp Control vol 4 no 4 July Aug 1981 pp 406 414 Vlassenbroeck Jacques A Chebyshev Polynomial Method for Optimal Control With State Constraints Automatica
11. improve the initial estimate or try another method The shooting algorithm runs only a Newton method In general a finite element solution should be obtained before the shooting algorithm 1s attempted Output After a successful finite element run is executed the user is prompted to save a variable called yout This variable is the same length as the user supplied yin and contains the converged values of the solution vector To save this variable use the command save yout dat yout ascii 20 PROBLEM MAC all f files MACSYMA all o files COMPILER VTOTSINFO M namcom nml plant f plant mex4 FMEX USER NAMCOM NML SUPPLIED VALUES MATLAB VIOTS ALGORITHM Figure A2 Flowchart of problem setup 21 After completion of a finite element solution the user is always prompted to run another problem with a different number of elements The number of elements is usually increased to obtain better accuracy but the number of elements may be decreased The user must input the number of elements as a vector of a length that corresponds to the number of phases When the number of elements is increased or decreased code convergence is not guaranteed After completion of all finite element or shooting runs the program stores a matrix of values in yall This matrix is used for plotting and it can be saved in the same way as yout except the user is not prompted to do so The
12. order necessary conditions for optimality are derived by requiring J to be stationary The conditions are ref 7 KO go x wl 7 2 Ou __ OH _ _ 1 A 3x Hx 8 aH Ho 9 QuO 7 T where each equation holds for t lt t lt tj4 and i 1 2 N The boundary conditions are 4 0 10 Od AUD ep L j22 3 N 1 11 8x G 1 t i Y Od NOME 21 2 N 12 ti dx t i 2 N 12 and the transversality conditions are Od A ot HY t1 0 13 oP HHD t HO ru 0 i 2 3 N 14 9 dt HW tva 0 15 The optimal control problem defined above is a nonlinear MPBVP The solution to the MPBVP yields a stationary point of J1 or a candidate optimal solution The problem can now be extended to include control and state inequality constraints on the system Control constraints see a standard optimal control text such as ref 7 for details are defined as a function of the states and the control where the control appears explicitly but the states may not of the form g x u lt 0 16 To solve this problem the constraint g 1s adjoined to the cost function with a Lagrange multiplier function p t This augmentation is equivalent to redefining the Hamiltonian of the system H as H L A f p g 17 The necessary conditions in equations 7 through 15 remain unchanged when the new definition of H is used However the multiplier p requires additional necessary con
13. scaled 3 No boundary conditions on eta occur in psilist because eta has a different unknown constant value in each phase 4 When multiplied by the differential equations the variable eta 1s squared only to ensure a positive value Therefore the returned value of eta is the square root of the length of the phase 9 The differential equation for eta 1s 0 because eta is a constant 6 Because each phase has been scaled the final time of the first phase is 1 and the final time of the second phase is 2 as indicated in tsilist The name list file which defines the values of the vehicle parameters and physical constants is namcom tvac 2594963 0d 00 spimpi 430 55d 00 spimp2 430 55d 00 earmu 3 98601d14 re 6 378145d6 grav 9 81d 00 rmass0 890149 09d 00 47 ho 0 0d 00 vO 20 0d 00 gamO 0 157d 01 hf 148011 1d 00 vf 7854 0d 00 gamf 0 0d 00 dropma 29920 0d 00 rhoO 1 35924d 00 sref 5 518d 01 betar 0 1405594d 03 ae 3 823d 00 pO 97136 2d 00 betap 0 141864 03 engi 5 0d 00 eng2 1 0d 00 ca 0 35d 00 cna 0 045d 00 end The initial estimate from trial and error is loaded into the vtotsinfo m file with the load command Also optional variables that may be defined in the vtotsinfo m file are timestate and scale These variables are for plotting only The vtotsinfo m file is prob_name NLS jocv 32 32 load yout3232 dat yin yout3232 m0 8901
14. solution to the optimization problem A user s manual some examples with results and a brief description of the individual subroutines are included in this report Introduction Background The optimal control problem featured in this report is described as follows Consider a dynamical system defined by a finite dimensional set of ordinary differential equations and assume a finite dimensional vector of time varying control variables The optimization problem is to choose the control variables to satisfy the given boundary conditions while a given performance index or cost functional is minimized or maximized Methods available for the numerical solution of optimal control problems generally fall into two distinct categories direct and indirect Direct methods which involve parameter optimization directly minimize the performance index by varying the values of the parameters Indirect methods on the other hand minimize the performance index indirectly by satisfying first order necessary conditions for optimality that are established from the calculus of variations The direct approach to the solution of optimal control problems requires parameterization of the control and state time histories and results in a nonlinear programming problem to solve The choice of parameterization schemes is not unique and success of the direct methods has been achieved with Hermite polynomials ref 1 Chebyshev polynomials refs 2 and 3 single te
15. the time line Remember that this is a three phase problem with coincident nodes defined at 0 2 and 0 7 sec These times are Just estimates as to when the constrained arc starts and ends After the tablel routine is used estimates for the multipliers m are included as the last column of yin Note that because the constraint is assumed to be inactive in the first and third phases the multiplier is necessarily 0 The reshape command is used to produce a column vector Finally estimates are made for the discrete multipliers v and the final times of each phase The state costate control Hamiltonian and 07H u eigenvalue histories are shown in figures B6 through B10 respectively for the case jbev 4 4 4 Notice that in figure B6 the state xl x does not violate the given constraint of 0 1 Also v which is denoted with x2 in figure B6 and u in figure B8 both remain at 0 during the constrained phase Figure B7 shows that both costates have discontinuities at the start of the second phase due to the tangency conditions specified in psilist these discontinuities are part of the necessary conditions listed in the section entitled Generalized Optimal Control Problem Finally in figure B9 the Hamiltonian is not constant 1n the first and third phases This lack of consistency indicates that the exact solution has not been found as expected The Hamiltonian becomes constant as more elements are used Shooting cannot be used on
16. the number of states nz the number of controls nu the number of control constraints np the number of state constraints nq the number of state boundary conditions length of psi m c the number of time boundary conditions length of tsilist tbc and the number of phases nph Also the variable jbev defines the number of elements per phase The formula for determining the number of initial guesses for single phase problems is 2nz nu 2np nq Jbev 2 mbc tbc 1 For example a single phase problem with three states two controls one control constraint zero state constraints four state boundary conditions one time boundary condition and five elements would require an initial guess file of length 8 3 2 42 5 2 4 1 1 16 For multiple phase problems there is an obvious extension to this formula Unknowns occur at the midpoint of the elements in each phase and at the endpoints of each phase Two coincident nodes appear at the juncture of phases Although these nodes occur at the same instant the values of the variables states costates and controls may be different In fact this discontinuity in one or more variables often requires introduction of the additional phases The assembly of the initial guess vector 1s similar to the single phase process Sets of unknowns for the first phase are assembled as described above for the single phase problem Next before the values of v are added sets of unknowns are added for the second
17. the time for reviewing instructions searching existing data sources gathering and maintaining the data needed and completing and reviewing the collection of information Send comments regarding this burden estimate or any other aspect of this collection of information including suggestions for reducing this burden to Washington Headquarters Services Directorate for Information Operations and Reports 1215 Jefferson Davis Highway Suite 1204 Arlington VA 22202 4302 and to the Office of Management and Budget Paperwork Reduction Project 0704 0 188 Washington DC 20503 1 AGENCY USE ONLY Leave blank 2 REPORT DATE 3 REPORT TYPE AND DATES COVERED July 1993 Technical Memorandum 4 TITLE AND SUBTITLE 5 FUNDING NUMBERS Variational Trajectory Optimization Tool Set Technical Description and User s Manual WU 946 01 00 82 AUTHOR S Robert R Bless Eric M Queen Michael D Cavanaugh Todd A Wetzel and Daniel D Moerder PERFORMING ORGANIZATION NAME S AND ADDRESS ES 8 PERFORMING ORGANIZATION NASA Langley Research Center REPORT NUMBER Hampton VA 23681 0001 L 17166 SPONSORING MONITORING AGENCY NAME S AND ADDRESS ES 10 SPONSORING MONITORING National Aeronautics and Space Administration AGENCY REPORT NUMBER Washington DC 20546 0001 NASA TM 4442 11 SUPPLEMENTARY NOTES Bless Lockheed Engineering amp Sciences Co Hampton VA Queen and Moerder Langley Research Center Hampton VA Cavanaugh George Washington Univ
18. this problem because a constraint is imposed The assumption that this problem is composed of three arcs is true only for certain values of For a larger value of l for example 0 2 the trajectory only touches the constraint limit In that case the optimal trajectory would consist of only two phases with no tangency conditions The MACSYMA setup file to solve this problem for 0 2 would be This is a fixed final time second order state constraint problem with a touch point solution Section 3 11 Bryson and Ho 39 stlist x v ctlist u qlist 11 ellist 0 5 u 2 0 5 u 2 phi 0 thi 0 psilist x 1 1 x0 v 1 1 vO x 2 1 ellim x 1 2 x 2 1 v 1 2 v 2 1 x 2 2 xf v 2 2 vf tsilist thyme 2 1 delist Llv ul v u namcom x0 v0 xf vf ellim 0 1 xl vs time 1 x2 vs time 0 08 gt 0 06 7 x S 0 04 0 02 0 i 1 0 0 5 1 0 0 5 1 time time Figure B6 Constrained double integrator state histories 40 lambdal vs time 10 lambda vs time 20r 5F E 3 Ea 0r 5 2 OF E E E E 20 _ 5r 40 10 0 0 5 1 0 0 5 1 time time Figure B7 Constrained double integrator costate histories 40 Hamiltonian 0 4 0 35 0 3 0 25 0 2 0 15 0 1 0 05 0 ul vs time 0 0 5 1 time Figure B8 Constrained double integrator control history Hamiltonian vs time
19. trajectory and only one state constraint 1s active at a time Furthermore for problems with state constraints the control is assumed to be continuous across Junction points of constrained and unconstrained arcs Assuming continuity of the control is tantamount to saying that the Hamiltonian of the problem is regular that is a unique optimal control exists for a given state and costate time history The user should be aware of these assumptions and carefully study solutions obtained from the VTOTS package especially for constrained problems In general the user should be aware that with the finite element algorithm or any discretization algorithm the output is only a candidate solution to an extremal For problems with control constraints the user is not required to specify the switching structure of the constraint in other words the user does not need to know or specify in the problem setup when the constraints will be active or inactive However for problems with state constraints the user must know the order in which the constrained and unconstrained arcs occur Further if the program has active control and state constraints a switching structure must be 3 specified only for the state constraints Details and examples of handling constrained optimal control problems are presented in subsequent sections Finally neither the finite element algorithm nor the shooting algorithm handles optimal control problems with singular arcs Pu
20. unknown Because f varies monotonically from 0 to ty r varies monotonically from 0 to 1 Also note that dx dx dy dt Thus the differential equations for any fixed final time problem can be scaled from 0 to 1 by multiplying each equation by the desired known final time 15 files problem mac allell mac allf mac allg mac allq mac allphi mac allthi mac allpsi mac alltsi mac plant mac USER supplied VTOTS supplied MACSYMA Environment commands batch problem mac gentranin plant mac plant f gentranin allf mac allf f gentranin allell mac allell f gentranin gentranin gentranin gentranin gentranin allphi mac allphi f allthi mac allthi f allg mac allg f allg mac allq f allpsi mac allpsi f gentranin alltsi mac alltsi f quit files plantg f plant f allf f allell f allphi f allthi f allg f allq f allpsi f alltsi f VTOTS supplied MACSYMA generated commands 77 c all f 8 fmex plant f plantg f all o FORTRAN MATLAB file plant mex4 used by VTOTS for problem specific information Figure Al Commands for creating plant mex4 16 This method can be used even if the final time is not known a priori For a free final time problem define an extra state for example 7 to be solved by the system The differential equation for 7 is T
21. vO x 1 2 xf v 1 2 vf tsilist Lthyme 1 1 delist v ull namcom x0 vO xf vf The name list file namcom nml is namcom XO 0 0d 00 VO 1 0d 00 XF 0 0d 00 VF 1 0d 00 end 34 Finally vtotsinfo m is set up by the user as prob name unconstrained double integrator jbcv 2 yin rand 26 1 yin 26 1 0 Note that the last estimate in yin is the final time which is known to be 1 0 To run the shooting algorithm directly change vtotsinfo m to prob name unconstrained double integrator yin 1 1 1 1 ynu 1 1 1 1 utrial 2 nnode 0 time 0 1 Results for the finite element case with jbev 8 are shown in figures B1 through B5 The state histories for x and v are denoted by x1 and x2 in figure B1 the corresponding costate histories are in figure B2 and the control history is in figure B3 The nonsmoothness of the curves results from the use of a relatively coarse grid with eight elements The fact that the Hamiltonian in figure B4 is constant indicates that an accurate solution has been found Finally the eigenvalues of 92 H Ou are shown in figure B5 This graph is important because its value is always positive for all time and it provides a second order sufficient condition that a local minimum has been found State Constrained Double Integrator The problem described in the previous section 1s solved again this time with a state constraint of the form S x v x l
22. with first and second order time derivatives S x v x v S xvu v u Because the control u first appears in the second time derivative of S this parameter is a second order state constraint In order for VTOTS to handle this problem the user must decide on a switching structure for the constraint From the results of the unconstrained problem one might assume the solution 1s composed of an unconstrained arc followed by a constrained arc followed by an unconstrained arc For certain values of this solution is correct if so the MACSYMA setup file would be This is a fixed final time second order state constraint problem Section 3 11 Bryson and Ho stlist x v ctlist u glist C2 E11 qlist LL ul 1 35 ellist 0 5 u 2 0 5 u 2 0 5 u 2 phi 0 thi 0 psilist x 1 1 x0 v 1 1 vO x 2 1 ellim x 1 2 x 2 1 v 2 1 v 1 2 v 2 1 x 2 2 x 3 1 v 2 2 v 3 1 x 3 2 xf v 3 2 vf tsilist thyme 3 1 delist v ul v u v u namcom x0 v0 xf vf ellim 0 25 x vs time 1 x2 vs time m 0 5 0 15 x of 0 1 0 05 0 5 J 0 l j 0 5 1 0 0 5 l time time Figure B1 Unconstrained double integrator state histories 2 178 x10 15 lambdal vs time 4 lambda2 vs time S e 2 1 089 So Xl E E E 0 0 0 0 5 1 0 0 5 I time time Figure B2 Unconstrained double integrator costate histories 36 Hamiltoni
23. 11 Control constrained problem state histories 44 0 871 lambdal Vs time 0 lambda2 vs time 20 436 E t 5 E 0 2 5 i 10 0 5 10 time cing Figure B12 Control constrained problem costate histories 1 ul vs time 0 5 J s 0 0 5 4 1 0 5 10 time Figure B13 Control constrained problem control history A time state eta has been introduced so that each phase has a duration of 1 See section entitled Time Scaling The effect of the time state is multiplication of each differential equation and ellist by eta The resulting setup file is macsyma script problem mac Simplified NLS model nondimensionalized with aerodynamics analytical two phase stage problem no fairing drop mass drop at staging staging determined at optimal time PARAMETERS in namcom nml namcom tvac spimpi spimp2 earmu re grav rmassO hO vO gamO hf vf gamf dropma rhoO sref betar ae pO betap engi eng2 ca cna 45 STATE LIST stlist rmass h v gam eta CONTROL LIST ctlist alpha INTEGRAL COST LIST ellist 0 0 TERMINAL COST phi rmass 2 2 thi 0 PHASE BOUNDARY CONSTRAINTS LIST psilist rmass 1 1 1 h 1 1 hO hf v 1 1 vO vf gam 1 1 gamO rmass 2 1 rmass 1 2 dropma rmass0 h 1 2 h 2 1 v 1 2 v 2 1 gam 1 2 gam 2 1 h 2 2 1 v 2 2 1 gam 2 2 gamf 1 tsilist thyme 1 1 thyme 2 2
24. 49 09 hf 148011 1 vf 7854 0 scale 100 0 m0 hf vf 1 1 1 m0 hf m0 vf m0 1 100 0 m0 hf vf 1 1 1 m0 hf m0 vf m0 1 timestate 5 By defining timestate 5 the z axes of the output plots are scaled to the true lengths of each phase The variable scale redimensionalizes the states and costates Scaling the states automatically scales the costates Plots of the states and costates except for the piecewise constant time state and the corresponding time costate are shown in figures B14 and B15 for a finite element run of 32 elements in each phase The control history is shown in figure B16 48 10 xl 8000 6000 2 4000 2000 x105 x1 vs time 15 x105 x2 vs time os 4 1 N gt 0 5 0 500 0 500 time time x3 vs time 2 x4 vs time L 4 1 5r 7 L J 10 J 2 0 5 2 0 500 0 500 time time Figure B14 State histories for two stage rocket problem 49 lambdal lambda3 lambdal vs time lambda vs 0 0 time epp e 5 1 5 0 500 time time 20 lambda3 vs time 5000 lambda4 vs time 30 J 0 40 3 Z 5000 4 50 J 5 60 1 10000 2 70 15000 0 500 0 500 time time Figure B15 Costate histories for two stage rocket problem 04 ul vs time 0 2 H 4 of 0 2 F 2 0 4 0 500 time Figure B16 Control history for two stage rocket problem 50 Ap
25. NASA Technical Memorandum 4442 Variational Trajectory Optimization Tool Set Technical Description and User s Manual Robert R Bless Eric M Queen Michael D Cavanaugh Todd A Wetzel and Daniel D Moerder J ULY 1993 NASA Technical Memorandum 4442 Variational Trajectory Optimization Tool Set Technical Description and User s Manual Robert R Bless Lockheed Engineering amp Sciences Company Hampton Virginia Eric M Queen Langley Research Center Hampton Virginia Michael D Cavanaugh George Washington University Hampton Virginia Todd A Wetzel lowa State University Ames lowa Daniel D Moerder Langley Research Center Hampton Virginia The use of trademarks or names of manufacturers in this report is for accurate reporting and does not constitute an official endorsement either expressed or implied of such products or manufacturers by the National Aeronautics and Space Administration Contents List of Figures Abstract Introduction Background What Is the VTOTS V TOTS Software Capabilities Purpose and Overview of Report Symbols Technical Description of Methods Generalized Optimal Control Problem Finite Element Method Shooting Method Concluding Remarks Appendix A User s Manual Using MACSYMA for Problem Setup The setup file os Variable names to avoid Example setup file problem mac Creating the MATLAB plant module plant med Time Scaling The File Vtotsinfo m Variabl
26. OTIS is a three degree of freedom point mass optimization program that includes a six degree of freedom and multiple vehicle simulation ref 1 The user can simulate aircraft missiles reentry vehicles and hypervelocity vehicles The methods used were chosen to improve speed convergence and applicability of the OTIS over existing performance programs Both the POST and the OTIS are reliable and accurate programs but they specifically target aerospace applications What Is the VTOTS The VTOTS package 1s a set of optimal control algorithms each based on a common problem specific user setup and interface The two methods for solving optimal control problems are a 2 finite element and a shooting method Each method uses a symbolic mathematics package to organize the system equations and to calculate system Jacobians The VTOTS package also uses the finite element algorithm to obtain initial estimates for the more accurate shooting code Combining the finite element results with a shooting initial condition provides a fast solution technique for nonlinear optimal control problems The VTOTS package was designed to minimize the user programming needed to solve optimal control problems and still provide a quick accurate solution procedure Three software packages that are used by the VTOTS are described in the next section VTOTS Software The VTOTS optimal control algorithms use three computer languages 1 MACSYMA MACSYMA is a sym
27. able is optional In order to use the finite element method estimates must be provided for all unknowns Consider a single phase problem A set of unknowns occurs at the midpoint of each element denoted by z and also at the beginning and end of each phase z Each set of unknowns consists of in the following order the states x1 Xn the costates 34 An the controls uj Um the multipliers for each control constraint 44 ppp the slack variables for each control constraint k1 knp and the multipliers for the state constraints 71 nq There may not be any constraints therefore no multipliers or slack variables are re quired After these estimates have been assembled several more estimates are added to the end These estimates correspond to the discrete Lagrange multipliers v that adjoin the boundary conditions held in 4 and to the discrete multipliers v that adjoin the boundary conditions in tsi Finally an estimate for the final time is made after the multiplier estimates For brevity the set of unknowns for a problem with three states two controls one control constraint one state constraint four state boundary conditions and one time boundary condition 1s z X1 X2 X3 AL Ag AZ Uj u2 44 ki m and the format of the initial estimates for jbcv 5 is Z1 Zi Z2 Za Z4 Z5 22 Vi v2 V3 V4 vy tu A general formula can be defined for the size of the initial estimate file Name
28. an ul vs time T ul t 0 0 5 1 time Figure B3 Unconstrained double integrator control history Hamiltonian vs time 4 T T T T T T T 0 l l l l l l l l 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 time Figure B4 Unconstrained double integrator Hamiltonian 37 H uu EigenValues vs time 2 T T T T T T T s gt 5 m 3 T 0 l l l l l l l 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 time Figure B5 Unconstrained double integrator eigenvalues of Huu The problem has now been constructed as a three phase problem The same differential equations hold for each phase The variable glist that holds the state constraint information consists of three parts one for each phase For the first phase qlist is empty an indication that no constraints are active In the second phase the state constraint is assumed to be active therefore S is put in qlist Finally the third phase is unconstrained so qlist is empty again The user must specify the boundary conditions for this problem in psilist Recall that the tangency conditions S 0 and 0 discussed in the section entitled Generalized Optimal Control Problem must be satisfied at the start of the second phase These conditions are listed as the third and fifth conditions in psilist The fourth and sixth conditions specify that the states x and v are continuous that is the values at the end of the first phase equa
29. and subsequent phases At the juncture of phases the sets of unknowns may have identical values When all phases have been assembled 18 one v for each boundary condition in v and tsi and estimates for the final times of each phase are added to the end of the initial estimate vector The general formula nph 2na nu 2np nq p Jbcv 2 zl mbe tbc nph 1 1 may be used to calculate the length of the initial estimate file Shooting variables The VTOTS provides a shooting algorithm that may be run directly or automatically without user interface after a successful finite element run The setup outlined in this section describes how to run the shooting algorithm directly The VTOTS initializes the shooting startup automatically when the finite element shooting method is operating so that no additional setup beyond the finite element initialization is required As with the finite element method starting estimates must be provided for all shooting method unknowns which are the state and costate values at the beginning of each phase and at any user specified interior phase points nodes In addition this method requires Lagrange multipliers and a control estimate A summary of these estimates and the variables that specify the number and frequency of nodes is shown below and must be included in the file vtotsinfo m yin initial estimates for each phase and node this column vector must contain the states and costates of the fi
30. arc These boundary conditions are placed in because of these conditions the user must define the switching structure of the constrained arc Thus the user must decide when the trajectory enters and leaves the constraint boundary because each independent arc of the trajectory is a new phase with corresponding boundary conditions Without loss of generality all constrained problems can be set up as minimizing problems with the constraints defined as less than or equal to zero The VTOTS also requires this constraint format Finite Element Method The finite element method converts the continuous time necessary conditions into nonlinear algebraic equations The process for generating the algebraic equations is outlined below Full details of the method are described in reference 15 For simplicity the finite element method 1s outlined for a one phase problem that 1s one with no internal events To begin the derivation of the finite element equations the first variation of an augmented performance index 1s taken the resulting expression 1s integrated by parts so that no time derivatives of the states x or costates A appear Instead one time derivative of the variational states 6x and variational costates 6A appears This appearance identifies the simple choice of approximating functions Next shape functions or approximating functions for the states costates and controls are chosen With the expression that is developed for the first v
31. ard c3 MACSYMA command line variable storage di MACSYMA display line variable storage ez MACSYMA internal variable sequence emq variable name string reserved by VTOTS The user must avoid the variables sin cos log and exp because these strings are treated as the corresponding function names Also the user must avoid using the tangent function in the setup file because MACSYMA does not successfully convert this function to FORTRAN The user is responsible for ensuring that each variable name used in the MACSYMA problem setup does not begin with the letters 1 j k 1 m or n because these letters are reserved for integers in FORTRAN Do not use thyme as a variable except in thi and tsilist Also any MACSYMA keyword that 1s used as a variable name leads to unpredictable results The user should always check the MACSYMA output for error messages Example setup file problem mac This example will help the reader understand how the MACSYMA setup file is defined A complete optimal control problem example is presented in the section entitled A Detailed Example Consider this linear quadratic optimal control problem find u t to minimize the scalar performance index J where ie W t dt subject to x t x t u t with boundary conditions 13 The following file problem mac loads this problem into MACSYMA stlist x defines the state variable names ctlist u defines the control variable names glist LCI
32. ariation the simplest possible shape functions are chosen for the states costates and controls namely plecewise constant functions To begin the discretization scheme associated with this finite element method a time line is broken into a series of equal segments known as elements The length of each element is At ti to M where M is the number of elements The endpoints of each element are referred to as nodes We will denote the values of the states costates and controls at the element midpoints as barred symbols Similarly values at the nodes will be symbols with carets Figure 1 is an example of a time line that is broken into five elements only the state variables are labeled Nodal values at the beginning and end of a phase and at all midpoint values are treated as unknowns for the states costates and controls The remaining unknowns are the discrete multipliers v and the event times t See appendix A x X5 X3 X4 Xs Figure 1 Discretized time line The state differential equations that are discretized become X1 XI Afi x f x u gt 0 x 4f x34 AE 1 2 M T xm Sify 3 where f denotes the value of f at midpoint i The costate differential equations become H A x i BE StH X A U H 1 EY Athy Athy ae Ox 0 24Ai Hx Aig 5 Hx 21 2 M 1 Au S Ay A where H denotes the value of H at midpoint i The algebraic optimality condition be
33. ble Integrator for tips on how to get switching structure for state constrained problems When solving problem with control constraints do not choose zero as an initial guess for the multiplier and slack variable this choice causes a singular matrix Remember that all constrained problems must be minimization problems Any maxi mization problem can be transformed into minimization problem by multiplying the performance index by 1 In general avoid an estimate of zero for unknowns Remember to list all known continuity conditions on states for problems with multiple phases V TOTS cannot directly handle boundary conditions that contain states and time If this situation occurs introduce another state that corresponds to the time as shown in the section entitled Control Constrained Problem 23 Detailed Example Consider the transfer of a particle to a rectilinear path as described in section 2 4 of reference 7 The particle has constant acceleration a The problem is defined in terms of four states z coordinate y coordinate velocity in z direction velocity in y direction G gc M and one control B angle of acceleration vector measured positive from x axis he differential equations are given by ru you a cos 8 v asin 0 The goal is to maximize the velocity in the z direction after 20 sec All states are initially zero and the final velocity in the y direction is zero The final
34. blem which corresponds to position of time state timestate defined only for plotting purposes defining timestate automatically scales the x axis of the plots to the correct values of the independent variable scale matrix of scaling factors n by nph where n is the number of states plus costates and nph is the number of phases each row of the matrix scales the states and costates in the corresponding phase the z 7 element of scale multiplies the 7th state in the jth phase for a problem that has been nondimensionalized scale will dimensionalize the problem as an example see the section entitled A Two Stage Rocket Problem in appendix B Finite element variables The following variables are defined by the user in vtotsinfo m if the finite element algorithm is run jbev vector of number of elements in each phase vector length is equal to number of phases jbev determines the mesh density of the solution in each phase jbcv gt 1 this variable is required yin vector of initial estimates for all unknowns size and order of the initial guess are defined below this variable is required 17 converge variable that defines the convergence criterion default value is 1 x 107 sometimes useful to raise this convergence value if the code approaches a solution but does not reach it raising convergence value allows the user to look at the answer before full convergence 1s reached to see if the solution is being approached or not this vari
35. bolic mathematics package that computes analyt ical derivatives of mathematical expressions A V TOTS preprocessor was written in MACSYMA a language developed by Symbolics Inc to organize and calculate expressions needed by the VTOTS algorithms The preprocessor then translates these mathematical expressions into FORTRA N 2 FORTRAN The result of the VTOTS preprocessor is a series of FORTRAN sub routines that are written to disk Each subroutine is generated by the MACSYMA algorithm 3 MATLAB MATLAB is a computer language that specializes in matrix manipu lation and vector analysis The VTOTS program and associated algo rithms are written in MATLAB a language developed by Mathworks Inc The FORTRAN subroutines supplied by MACSYMA are compiled into a single problem specific module using MATLAB compiler The plant module 1s then accessed by MATLAB function files Capabilities The VTOTS package provides solutions to a variety of optimal control problems with both the finite element and shooting algorithms Both algorithms can solve nonlinear optimal control problems with multiple state or state rate discontinuities Also the boundary conditions can be any nonlinear function of the states The finite element algorithm but not the shooting algorithm solves problems with control and or state constraints The number of control constraints 1s arbitrary however it 1s assumed that the same number of constraints acts over the entire
36. cobian feocbvp m the main driver routine for the finite element code determines the prob lem parameters and prompts the user for a solution method fepsi m defines the error vector and Jacobian for the boundary conditions held in psilist fesolv m the driver subroutine to fill the error vector and Jacobian also solves the linearized system if appropriate geth m calculates the Hamiltonian H and the eigenvalues of 92 H Ou for plotting purposes inphas m defines the error vector and Jacobian for the elements on the interior of a phase jacob m finds an error vector for use by fsolve m morenode m uses the MATLAB linear interpolation routine to generate new initial estimates for feocbvp m 51 nodal m solve m stphas m timcond m unod m extracts nodal values of states costates and controls assembles these values with the appropriate time vector in a matrix called yall for plotting by plotter m the user may save yall and call plotter m directly if desired called by feocbv p m when the Newton method 1s chosen by the user determines the step size logic and convergence criteria defines the error vector and Jacobian for the equations at the beginning of a phase defines the error vector and Jacobian that corresponds to the boundary conditions held in tsilist and the boundary conditions on the Hamiltonian uses a Newton method to determine the nodal values of the control called by nodal m Shooting Method geths m
37. comes Hy x A u p n 0 Hy Xj Ay j p 7 0 1 2 wr eg M The remaining equations involve the state and costate boundary conditions and the transver sality conditions The same number of equations as unknowns appears in this formulation 8 Additional algebraic conditions are associated with control constraints The finite element algorithm handles the control inequality constraints g x u lt 0 by introducing a positive function k such that g k 0 The function k is referred to as a slack variable and becomes an unknown Note that when on the constraint g 0 therefore k 0 Additional unknowns associated with state constraints are listed in appendix A A finite element method yields an approximate solution to the optimal control problem From numerical experience the accuracy of the solution or closeness to the exact answer increases quadratically with the number of elements ref 15 however for a numerically accurate answer a shooting method is available Shooting Method The VTOTS includes shooting algorithm for solving the necessary conditions in equa tions 7 through 15 The solution technique converts the MPBVP for the Hamiltonian system eq 6 subject to equation 7 and boundary conditions eqs 10 through 12 into an alge braic root finding problem in the values taken on by x A and at the initial and terminal points of the trajectory and at internal events The procedure 1s accomplished
38. d time problem The costate A that corresponds to 7 is appended to A and is evaluated at t t with the appropriate modification of equation 23 The x 7 A and A variables and their propagation expressions eqs 22 23 and 28 are concatenated to form the system 1 Vp Votr A F V dC 29 vic Ix T XU ar FT If 0 HF 27H which satisfies the equation Y Vo Vy 20 30 where V is a concatenation of equations 19 through 21 reexpressed in components of Vo and V y Equation 30 is solved by expressing V y as V Vo with equation 29 and employing a Newton Raphson iteration to obtain Vo The jth iteration is dy NO Vo Vo Ga v vo G 0 1 31 The value Vo o the initial guess for the iteration is usually provided by boundary values from a converged finite element run For problems addressed to date with the VTOTS these values have proved to be sufficiently close to the shooting solution so that no line search was necessary in equation 31 The Jacobian matrix dY dVo in equation 31 is dv v aw dV 2 dVo 8Vy OV dVo 92 where A i f gt fi dF dV l d dV t 1 Iv qv 39 The use of equations 32 and 33 to obtain dW dVo rather than the use of direct numerical differentiation with respect to Vo is motivated by concern for numerical stability in integrating 10 V When the plant states x contain dissipative effects some eigenvalues of the adjoint
39. ditions For a minimizing problem the conditions are as follows a multiplier of zero when the constraint is not active g lt 0 and a nonnegative multiplier when the constraint is active g 0 Consider problems with state inequality constraints of the form S x lt 0 One of several methods available to solve problems with state constraints 1s to take total time derivatives of the constraint until the control appears explicitly this method requires substitution of the differential equations for the state rates If q time derivatives are required for the control to appear explicitly then the constraint 1s referred to as a qth order state inequality constraint Now the gth time derivative of the constraint plays the same role as the control constraint g x u above After a Lagrange multiplier function m t is introduced the Hamiltonian is ras H L 2 f n qu 18 where the following statements apply 1 The multiplier y 0 when the constraint is not violated S lt 0 2 The value d S dt 0 when the constraint is active S 0 3 The multiplier 7 gt 0 when the constraint is active if minimizing cost In addition to taking time derivatives of the constraint tangency conditions must be met at the point of entry onto a constrained arc These conditions are that S and the first q 1 time derivatives of S are zero at the beginning of a constrained arc Also the states must be continuous at the beginning and end of each
40. e interval to to t y discontinuities in the states as well as in the state equations may occur at interior points i e times between to and ty These interior initial and final points are referred to as events and the intervals between events are referred to as phases The time of event 7 is denoted as and the states and controls in phase z are denoted as x and ul The optimal control problem of interest in this report is to choose a control history that minimizes a performance index J and satisfies the state equations x O x wl and boundary conditions Elements of a performance index may be denoted with an integrand LO wl which is continuous and differentiable within each phase and a discrete function 5 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 N ft pug Tab x t5 x 2 x 2 x 3 s uva states stv Ml LO x a9 di 1 i21 subject to the state equation constraints ut t lt t lt tari 12 N 2 with boundary conditions specified as v x ty 63 xO 65 xO 68 aw itat sta O 3 With the introduction of Lagrange multiplier functions A t referred to as costates and discrete Lagrange multipliers v an augmented performance index J may be defined as Ji i p gt na L 4 AG ro 0 dt 4 i21 For convenience and H are defined as d u 5 and H LO AO O i21 2 N 6 The first
41. e name is also helpful For example K gt gt save bhofix5 dat yout ascii double saves the output in the file bhofix5 dat named for the Bryson and Ho fixed time problem with jbcv 5 The format is double precision ascii Usually a good procedure is to start with a small value for jbev and build up until the desired resolution is reached After saving the output type return lt cr gt at the K prompt to continue Next the code gives the option of changing jbev 27 Change JBCV y or n 2 y INPUT JBCV 10 Enter 1 to run a continuation method Enter 2 to run fsolve otherwise cr to run a Newton method A carriage return here restarts the program with an initial estimate automatically generated by linear interpolation of the preceding solution the initial error is 0 13357 Step size 1 the error is 0 0148538 the iteration number is 1 U step size the error is 2 86589e 05 the iteration number is 2 U step size the error is 2 56208e 10 the iteration number is 3 CONVERGED Total run time is 22 0878 Now is your chance to save yout K gt gt The initial error is smaller than it was in the first run and the problem converges in fewer iterations The run time 1s about the same because a larger system of equations 1s solved at every iteration After saving this new solution the user is again prompted to change Jbcv If the reply is n this time then the s
42. ell f allf mac creates allf f allg mac creates allg f allq mac creates allq f allphi mac creates allphi f allthi mac creates allthi f allpsi mac creates allpsi f alltsi mac creates alltsi f plant mac creates plant f and namcom nml The FORTRAN file allell f evaluates the cost integrand L for all phases Similarly allf f evaluates the right side of the differential equations for all phases allg f and allq f evaluate the constraints allphi f and allthi f evaluate the discrete cost terms allpsi f and alltsi f evaluate the boundary conditions plant f is the master routine that coordinates calls to the other FORTRAN files and namcom nml contains the variables in namcom 14 One additional file plantg f is required to construct the plant module This file is supplied and does not require changes by the user The file plantg f is a gateway file to pass information between the FORTRAN routines and MATLAB The commands for running the MACSYMA preprocessor are batch problem mac gentranin plant mac plant f gentranin allf mac allf f gentranin allell mac allell f gentranin allphi mac allphi f gentranin allthi mac allthi f gentranin allg mac allg f gentranin allq mac allq f gentranin allpsi mac allpsi f gentranin alltsi mac alltsi f quit This sequence of commands can also be placed in a file batchfile mac for example and batched at the system level comma
43. equations that form a multipoint boundary value problem MPBVP which can be solved with simpler root finding techniques Analytical solutions to a MPBVP are generally unobtainable except for the simplest problems hence numerical methods are usually employed The two main methods for solving a nonlinear MPBVP are shooting and quasi linearization methods Shooting methods refs 10 through 12 are frequently used and can be described as follows The differential equations and the known initial conditions are satisfied at each stage of the process but the final conditions are not satisfied A nominal solution is generated by guessing the missing initial conditions and integrating the differential equations forward to reduce the error in the final conditions at each iteration Quasi linearization methods refs 7 and 13 are used to choose 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 and costate equations are 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 methods include the method of adjoints ref 14 and finite element methods ref 15 The system of equations in these methods is typically solved by Newton Raphson ref 16 or continuation algorithms ref 17
44. ersity Hampton VA Wetzel Towa State University Ames TA 12a DISTRIBUTION AVAILABILITY STATEMENT 12b DISTRIBUTION CODE Unclassified U nlimited Subject Category 18 13 ABSTRACT Maximum 200 words This report briefly describes the algorithms that comprise the Variational Trajectory Optimization Tool Set VTOTS package The VTOTS is a software package for solving nonlinear constrained optimal control problems from a wide range of engineering and scientific disciplines The WTOTS package was specifically designed to minimize the amount of user programming in fact for problems that may be expressed in terms of analytical functions the user needs only to define the problem in terms of symbolic variables This version of the VTOTS does not support tabular data thus problems must be expressed in terms of analytical functions The VTOTS package consists of two methods for solving nonlinear optimal control problems a time domain finite element algorithm and a multiple shooting algorithm These two algorithms under the VTOTS package may be run independently or jointly The finite element algorithm generates approximate solutions whereas the shooting algorithm provides a more accurate solution to the optimization problem user s manual some examples with results and a brief description of the individual subroutines are included in this report 14 SUBJECT TERMS 15 NUMBER OF PAGES Optimal control algorithm Finite elements Shooting m
45. es nx 4 Number of controls nu 1 Number of control constraints np 0 Max No of active state constraints nq 0 Number of state b c s mbc 6 Number of time b c s the 1 Number of phases nph 1 Number of elements in phase number 1 jbcv 1 5 1 Finite Elements READY 2 Shooting 00707 3 F E Shooting READY 4 Exit Program READY INPUT SELECTION gt gt The comment line is displayed and is followed by a brief summary of the important parameters for this problem Next a list of options appears with a READY message indicating the options that can be selected Because the different methods require different initial guess forms not all methods are available at once In this example the most powerful option for solving unconstrained optimal control problems is demonstrated Starting from the initial guess defined in vtotsinfo m the finite element algorithm is successfully run Initial estimates can then be generated for the shooting algorithm to produce an essentially exact answer to the problem The finite element shooting option is method 3 and is selected by typing 3 cr at the prompt This selection leads to the following message Enter 1 to run a continuation method Enter 2 to run fsolve otherwise l
46. es common to the finite element and shooting algorithms Finite element variables Shooting variables Overview of Problem Setup Solution Method Options Output Program Diagnostics Helpful Hints Detailed Example Appendix B Additional Example Files The Unconstrained Double Integrator State Constrained Double Integrator Control Constrained Problem A Two Stage Rocket Problem Appendix C Programmer File Reference List VTOTS Driver Subroutines Finite Element Method Shooting Method References ili List of Figures 1 Al A2 A3 A4 A5 A6 AT B1 B2 B3 B4 B5 B6 BT B8 B9 B10 Bll B12 B13 B14 B15 B16 Discretized time line Commands for creating plant mex4 Flowchart of problem setup State histories Costate histories Control history Hamiltonian history integral cost plus adjoined dynamics measure of global convergence of algorithm Eigenvalues of the second partial of Hamiltonian with respect to control second order sufficient condition nconstrained double integrator state histories nconstrained double integrator costate histories nconstrained double integrator Hamiltonian U U Unconstrained double integrator control history U U nconstrained double integrator eigenvalues of Huu Constrained double integrator state histories Constrained double integrator costate histories Constrained double integrator control history Constrained double integra
47. ethods 56 16 PRICE CODE AQ4 17 SECURITY CLASSIFICATION 18 SECURITY CLASSIFICATION 19 SECURITY CLASSIFICATION 20 LIMITATION OF REPORT OF THIS PAGE OF ABSTRACT OF ABSTRACT Unclassified Unclassified NSN 7540 01 280 5500 Standard Form 298 Rev 2 89 Prescribed by ANSI Std 739 18 298 102
48. hooting algorithm is started and the convergence is displayed on the screen SHOOTING iteration 1 abs error 2 73e 02 iteration 2 abs error 6 89e 04 iteration 3 abs error 4 18e 07 iteration 4 abs error 5 26e 13 total time 165 74 seconds The user is next given the option to plot any of the following states costates controls Hamiltonian and eigenvalues of the second partial of the Hamiltonian with respect to the controls do you wish to plot the results y or n y do you wish to see the state histories plotted y or n y 28 After each plot of four histories the following message appears type print lt cr gt to get a hardcopy of graphs type return lt cr gt to see rest of states K gt gt Typing print cr would print the graphics window of four plots to the MATLAB printer In this case typing return lt cr gt goes to the option for displaying the costate histories because only four states are available Next follow the other plot options The plots produced by VTOTS for this example are shown in figures A3 through AT do you wish to see the costate histories plotted y or n y do you wish to see the control histories plotted y or n y do you wish to see the Hamiltonian plotted y or n y do you wish to see Eigenvalues of H uu plotted y or n y These and all remaining plots in this report reflect the plots produced automatically by the VTOTS The states labeled x1 x2 xn in figure A3 correspond
49. ith the variable names but without the values To run the problem an initial guess must be supplied in vtotsinfo m prob_name BHO FIX B amp H Fixed time prob jbcv 5 tab 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 0 1 1 0 100 0 100 0 10 0 0 0 1 0 1 0 1 0 1 0 0 1 t 0 1 3 5 7 9 1 0 yin table1 tab t yin reshape yin 63 1 yin yin ones 7 1 1 0 scale 20 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 The first line of vtotsinfo m gives a comment that is displayed when the problem is run This comment is an optional declaration by the user The second line defines the number of elements in the first finite element run This variable must be defined if the finite element algorithm is used The user has the option of increasing or decreasing Jbcv if the first run is successful The next four lines demonstrate the ability of the MATLAB tablel function to create an initial guess from linear interpolation The matrix tab consists of 2 rows and 10 columns The first column is an independent variable that starts at 0 and ends at 1 The next nine columns are in order estimates at the states four columns estimates at the costates four columns and estimates at the control one column In this example only crude estimates for the beginning and ending values of the variables are made and tablel draws straight lines between them For example in the third column of tab the estimated value of the second state y is 0 at t 0 and 100 at t
50. l the values at the start of the second phase Continuity conditions are also listed at the junction node of the second and third phases No tangency conditions are required at the end of a constrained arc Finally the final time of the third phase is specified as 1 but no information is known about when the first and second phases end These times which are estimated by the user in yin are determined by the V TO TS The name list file for this problem is namcom XO 0 0d 00 VO 1 0d 00 XF 0 0d 00 VF 1 0d 00 38 ELLIM 0 1d 00 end and vtotsinfo m is as follows prob name state constrained double integrator jbev 1 1 1 load yall8 dat yin ya118 1 3 t 0 1 2 2 45 7 7 855 1 0 yin tablei yin t eta 0 0 0 1 1 1 0 0 0 yin yin ones 9 3 eta yin reshape yin 54 1 yin yin ones 11 1 yin 66 0 2 yin 67 0 7 yin 68 1 0 The answer from the unconstrained problem saved in the variable yall8 dat has been used to generate initial estimates of the states for the constrained problem Usually the costate and control histories change drastically as compared with the unconstrained case and are not useful for estimates The matrix yall is loaded into vtotsinfo m and then the variable yin is defined as the matrix containing all rows and the first three columns of yall These columns are the time and the two states Next a new variable t is defined to locate the points of unknowns along
51. n n multiplier vector NT time scales A costate vector T multiplier vector U vector of Lagrange multipliers discrete portion of augmented performance index a discrete portion of performance index ap vector of boundary condition expressions Abbreviations CTOP Chebyshev Trajectory Optimization Program I identity matrix MPBVP multipoint boundary value problem NPDOT Nonlinear Programming for Direct Optimization of Trajectories OTIS Optimal Trajectories by Implicit Simulation POST Program to Optimize Simulated Trajectories VTOTS Variational Trajectory Optimization Tool Set Technical Description of Methods In this section a nonlinear constrained optimal control problem is defined Then a brief description of a finite element method and a shooting method is presented to solve the optimal control problem Further details of these methods are given in the cited references Generalized Optimal Control Problem An optimal control problem is defined below First the notation is defined and the first order necessary conditions for unconstrained problems are derived Then the inclusion of constraints on the system is considered and the additional conditions for optimality are defined Consider a system that is defined from initial time to to final time ty 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 th
52. nd prompt by typing the following batch command macsyma lt batchfile mac gt std out amp where std out will contain MACSYMA run time information and error messages These files must then be compiled with the system FORTRAN compiler On the Sun systems the commands are as follows 77 c all f amp Then the plant f and plantg f files must be compiled with a MATLAB compiler and linked to the other object code with the command fmex plant f plantg f all o The result is the plant mex4 file which can be moved to a convenient working directory and accessed by MATLAB routines in much the same way that functions are called Figure 1 shows a summary of the commands for creating the plant module plant mex4 Any plant module that is acceptable for use with the shooting algorithm will also work for the finite element algorithm however the converse may not be true For example a plant module that includes constraints will work with the finite element algorithm but not with the shooting algorithm Time Scaling The finite element algorithm in the VTOTS does not require special scaling of the time parameter However in order to run the shooting algorithm in the VTOTS the user must scale the time of each phase to a length of one This procedure requires the conversion of free time problems to fixed time problems The variable 7 is defined such that t ty where t is the independent variable and tp is the final time possibly
53. nitial estimate file is not the proper length The word READY appears next to each option if the initial estimate is the proper size Choosing an option without a READY results in errors When the option for a finite element algorithm is chosen the user must decide between three different solution methods to solve the algebraic equations The user is prompted to choose between a continuation method MATLAB s fsolve algorithm ref 20 and a Newton method The continuation method is a simple type of homotopy described in reference 21 This option is the most robust of the three methods that is it allows for the least accurate initial estimate and still finds a solution but it is also the slowest After the continuation method is completed the Newton method is automatically called to obtain the solution In certain cases the integrator for the continuation method is interrupted and gives an error message like Singularity likely at t 0 456 The Newton method is called at this time and may converge on the solution in such a case the message can be ignored The fsolve algorithm in MATLAB is a Newton method with a line searching algorithm The fsolve algorithm is generally not as robust as the continuation method but it does run faster The Newton method is the fastest of the three solution methods but it requires the best initial estimate Generally the Newton method should be attempted first If the program does not converge then either
54. no control constraints specified qlist CC no state constraints specified phi 0 no discrete cost on states defined thi 0 no discrete cost on times defined ellist x 24u 2 quadratic cost function psilist Notice the indices for boundary conditions x 1 1 x0O 1 1 1st phase initial time x 1 2 xf 1 2 1st phase final time 1 The same index scheme is used in phi tsilist thyme 1 1 1 final time of the first phase delist LL x u 11 differential equations namcom x0 xf these variables are found in the FORTRAN namelist The MACSYMA comments delimited by and on the right do not need to appear Creating the MATLAB plant module plant mex4 In this section a listing of file names and UNIX commands is given to show how to use the MACSYMA preprocessor and how the MACSYMA produced files are compiled into a single problem specific module Several versions of MACSYMA and FORTRAN are available and these vary from one machine to another The existing versions of MACSYMA version 417 100 FORTRAN version 1 4 and MATLAB version 3 51 described in this report are specific to Sun SPARCstation IPC and IPX workstations After the problem specific information has been set up in a file such as problem mac the MACSYMA preprocessor can be run The MACSYMA preprocessor consists of the following nine MACSYMA script files that create FORTRAN files allell mac creates all
55. pendix C Programmer File Reference List Below is a list and brief description of all the MATLAB m files in the VTOTS First a handful of files used by both algorithms is listed Then the m files specifically used by the finite element code are listed followed by the shooting code m files VTOTS Driver Subroutines alert m issues error or warning messages as prompted by vtots m fems m produces the initial estimate for shooting after a finite element run plotter m produces plots of the output after a successful finite element or shooting run called by vtots m but may be called by the user directly printfull m modified print m file that prints the plots produced by plotter m in landscape mode type printfull instead of print at the MATLAB prompt vtots m the main driver routine for VTOTS reads vtotsinfo m checks for a proper initial guess calls the finite element and shooting algorithms and calls the plotter m routine Finite Element Method enphas m defines the error vector and Jacobian at the end of a phase errorvec m finds an error vector for use by fsolve m febc m defines the error vector and Jacobian for the costate boundary conditions at the beginning and end of each phase fecontin m provides MATLAB s ode45 m integrator with the differential equa tions needed to solve the system of equations with a simple continuation method fejac m calls stphas m inphas m and enphas m to fill in most of the error vector or Ja
56. raint problem Section 3 8 Bryson and Ho glist u ulimu utulim1 glist L l1l stlist x t ctlist Lu ellist 0 5 u 2 phi 0 5 x 1 2 2 thi 0 psilist x 1 1 x0 t 1 1 tsilist thyme 1 10 delist u 14t 3 t 2 17 111 namcom xO ulimu uliml and the corresponding namcom nml file is namcom XO 19 945596d 00 ULIMU 1 0d 00 ULIML 1 0d 00 end Because the state equation is an explicit function of time nonautonomous an extra state is introduced This extra state t imitates an independent variable because it runs from 0 to 10 The MATLAB tablel function generates initial estimates when information is known about some variables See the section entitled A Detailed Example One vtotsinfo m file that worked is listed below prob_name BHO FIX B amp H control constraint prob jJbcv 5 tab 0 00 19 9 0 0 1 1 1 1 1 1 1 10 0 0 0 10 0 1 1 1 1 1 1 11 t 0 1 3 5 7 9 10 yin tablel tab t yin reshape yin 63 1 yin yin ones 3 1 10 0 Results for the states costates and control are shown in figures B11 through B13 Notice that the control history does not violate the specified constraints A Two Stage Rocket Problem For one last example using a time state to set up a shooting problem consider the following model of a two stage rocket The states chosen are mass m altitude h velocity V and flight path angle y the control is the angle of attack a so the d
57. ree files in this procedure are provided by the user problem mac namcom nml and vtotsinfo m The MACSYMA file in this case problem mac can have any name however the other two files must be named namcom nml and vtotsinfo m The first step is to process the MACSYMA file as described in the section entitled Using MACSYMA for Problem Setup to produce the files plant f and namcom nml and eight other FORTRAN files At this point the name list file namcom nml has a list of parameters with no values The user must edit this file and input the parameter values Next the all f files should be compiled to form object files The object files the plant f file and the plantg f file are combined into a file called plant mex4 through the use of the MATLAB fmex utility At this point the command vtots in MATLAB causes VTOTS to access the files plant mex4 vtotsinfo m and namcom nml The user is prompted for several options which are discussed in the next section Solution Method Options After the plant mex4 namcom nml and vtotsinfo m files are created the user is ready to start MATLAB and run the VI OTS by typing in vtots at the MATLAB prompt A menu appears that lists four solution method options The user can choose the finite element algorithm the shooting algorithm or the finite element algorithm followed by the shooting algorithm The fourth option is to exit the program a useful choice if the name list file is not set properly or if the i
58. rm Walsh series ref 4 and splines ref 5 After the parameterization scheme is chosen a parameter optimization algorithm is used to determine the free parameters These algorithms are In common use today and include variable metric techniques or quasi Newton methods ref 6 and variations on gradient methods Gradient methods refs 7 and 8 were developed to surmount the initial guess difficulty associated with other methods such as Newton algorithms The gradient methods are characterized by iterative algorithms for improving estimates of the state and control time histories First order gradient methods rapidly improve the state and control histories during early iterations when sufficiently far from the optimal solution however these methods exhibit only linear convergence close to the solution Second order gradient methods provide quadratic convergence but are more sensitive to initial guesses Conjugate gradient methods exploit the approximately quadratic variation of the objective function near the solution to accelerate convergence Reference 9 contains a thorough description of the gradient method and other algorithmic methods in optimal control Because the direct method is presented as a nonlinear programming problem the solution is much more difficult to obtain especially from a software standpoint Conversely when the indirect method satisfies the first order necessary conditions the problem is converted into a system of
59. rpose and Overview of Report This report describes the finite element and shooting algorithms and explains how to solve optimal control problems with the VTOTS The section Technical Description of Methods defines an optimal control problem and provides a technical description of the finite element and shooting algorithms A brief discussion of each algorithm and the VITOTS package is then presented Concluding Remarks summarizes the unique features of the VTOTS Appendix A is auser s manual for solving optimal control problems with the VTOTS and includes an example and some helpful hints Appendix B contains several additional example files and output for problems that are solvable with the VIOTS Finally appendix C briefly describes the VTOTS MATLAB files Symbols F vector of right sides for state and costate equations f right side of differential equations g state and control constraints H Hamiltonian J scalar performance index Jy scalar augmented performance index k slack variable L integral portion of performance index M number of elements m number of controls N number of phases n number of states q order of state inequality constraint S state inequality constraints tr final time t time at ith event u control vector V vector containing states and costates x state vector x state vector at event points x state vector at midpoints x state time derlvative vector x state variation GA costate variatio
60. rst phase followed by the states and costates of the first node etc length of yin 2nz nnode nph this variable is required utrial control estimate for the system at the initial time this variable is required nnode column vector that contains number of nodes in each phase the first element in the vector specifies the number of nodes in the first phase etc a 0 is needed if the phase does not contain nodes this variable is required ynu column vector containing the Lagrange multipliers length of ynu mbc this variable is required time matrix in which each column holds node times for each phase including a 0 to start the phase and a 1 to end it shorter columns fewer nodes in a particular phase must be padded with 0 s to make the matrix square for a single phase problem the vector time must be a column vector this variable is required err specifies the integrator error default is 1 x 1079 this variable is optional For example consider a two state two phase problem with two nodes in the first phase at times 0 2 and 0 6 and one node in the second phase at time 0 5 Three boundary conditions exist nnode 2 1 time 0 2 6 1 0 5 10 err 1e 6 utrial 5 yin 12345678 9 10 11 12 13 14 15 16 17 18 19 20 ynu 1 2 3 Notice that the trailing 0 1n the time variable makes the matrix rectangular 19 Overview of Problem Setup The problem setup procedure is illustrated in figure A2 Th
61. save command may be evoked after completion of the plotting The matrix yall contains the following columns of data the time the states the costates the controls the Hamiltonian and the eigenvalues of the second partial derivative of the Hamiltonian with respect to the controls Because the Hamiltonian is constant for each phase at the exact solution the value of the Hamiltonian should be on the order of 1 x 107 for the shooting code which uses an integrator with an error tolerance of 1 x 107 The finite element algorithm is not as accurate unless the number of elements jbcv is large The eigenvalues are important because they serve as a second order necessary condition for a minimum or maximum The eigenvalues should be positive everywhere for a minimization problem and negative everywhere for a maximization problem Although the multipliers for the constraints are not available in yall these values are available in the vector yout The plotter routine may be called directly by the user if yall is saved To call the plotter enter plotter nx nu yall at the MATLAB prompt Each of these arguments should be in the workspace after a successful run by either the finite element or the shooting algorithm Type help plotter for more information Program Diagnostics The following list shows some potential errors that can occur 1 Common MACSYMA mistakes are a Use of an equal sign instead of a colon b Not ending a line with a
62. semicolon the result of which is usually a MACSYMA error message stating that some variable is not an Infix operator c Use of wrong number of brackets when defining MACSYMA variables The variables delist glist and qlist are lists of lists that require an extra set of brackets Incorrect number of brackets usually results in the message part fell off end d Failure to compile an indication of a mistake in the MACSYMA setup file 2 Segmentation violation during a call to plant mex4 is caused by a mistake in the MACSYMA setup file 3 No READY light by any of the solution options except 4 Exit Program indicates that the initial estimate 1s not the correct length Choosing the desired option should point to the discrepancy 4 Failure to provide values for the name list can produce strange results These values are held in the file namcom nml 5 A warning that a matrix is singular or badly scaled given during a Newton method means that the Jacobian matrix 1s singular and cannot be inverted by MATLAB In this case either the initial estimate leads to a singular matrix the problem is poorly defined or the 22 problem is singular at the solution Fixing this problem requires remodeling the problem or changing the initial estimate file no converge in unod mindicates that one of the control values during an interpolation routine was not found This condition is generally caused by a bad solution vector altho
63. sible This approximation may be used to start the shooting algorithm Finally a brief description of the VTOTS MATLAB files is included NASA Langley Research Center Hampton VA 23681 0001 April 19 1993 11 Appendix A User s Manual This appendix describes how to set up run and solve optimal control problems with the VTOTS In particular the development of three files that are needed to run VTOTS is described These files are a plant module plant mex4 a name list file namcom nml and a MATLAB file vtotsinfo m The first stage in using the VTOTS system is to set up the optimal control problem in a MACSYMA readable form this step is the creation of a file that defines specific MACSYMA variable names equation lists cost expressions and lists of parameters that define the problem The MACSYMA setup file and commands for producing the MATLAB FORTRAN interface are described in the next section The section entitled Time Scaling discusses how and when to scale the independent variable of the problem The Vtotsinfo m section describes the user supplied MATLAB file that is read in by the VI OTS That section includes a discussion of the initial guess vector that is required for the finite element and shooting algorithms An overview of the steps required to set up the VT OTS files is provided in the Overview of Problem Setup section The solution methods available to the user are described in the section entitled Solution Me
64. silist will be zeroed in the solution of the problem tsilist list of boundary conditions on time for each phase may be empty glist list of state and control inequality expressions qlist list defining the switching structure and the qth time derivative of a state constraint 12 namcom list of scalar FORTRAN variables for placement in the parameter name list useful for parameters that vary across a family of problems for example an initial condition could be put in namcom and then changed without having to rerun MACSYMA namarray list of lists of variables appearing in the name list that need to be dimensioned in FORTRAN variables are expressed in namarray with the correct dimension for example namarray a 3 b 4 7 dimensions a at 3 b at 4 f at 7 namarray is optional The variables phi and psilist have a common convention for defining event conditions In these variables a state name followed by two indices is used The first index denotes the phase number and the second denotes the initial or final time of the phase 1 for initial and 2 for final The variables delist glist and qlist are lists of sublists They contain one sublist for each phase Refer to the section entitled Additional Example Files for further clarification Variable names to avoid The following variables cause errors that may not be detectable by the MACSYMA preprocessor In the following denotes a number and denotes a wild c
65. t cr gt to run a Newton method 26 Hitting a carriage return at this prompt begins execution of the Newton method A sample of the execution follows the initial error is 8 14201 Step size 1 the error is 3 66657 the iteration number is 1 step size 1 the error is 0 508319 the iteration number is 2 step size 1 the error is 0 220737 the iteration number is 3 step size 1 the error is 0 0100028 the iteration number is 4 step size 1 the error is 7 62849e 06 the iteration number is 5 step size 1 the error is 2 0227e 12 the iteration number is 6 CONVERGED Total run time is 20 5638 Now is your chance to save yout K gt gt The given initial estimate converged with the Newton method in six iterations and required 20 5638 sec of run time The step size listed in the left column refers to the Newton method line search step size If the error is not reduced with the current step size then the step size is reduced The step size is continually reduced until the error improves decreases In this case all iterations improved the error Convergence is obtained when the error is less than 1 x 107 unless another value is set by the user in vtotsinfo m with the variable converge See the section entitled Finite Element Variables The output after a problem converges should be saved because many successful initial estimate files are generated by slight changes to output files from similar problems A descriptiv
66. thod Options The Output section describes the output that is available to the user when a VI OTS run is successfully completed Some program diagnostics and helpful hints are provided Finally a detailed example of the use of the VTOTS to solve an optimal control problem 1s presented Using MACSYMA for Problem Setup The first step in solving an optimal control problem with the VTOTS is to set up the problem in MACSYMA readable form This process separates the dynamics boundary conditions and performance index of an optimal control problem and assigns these expressions to MACSYMA specific variables A general problem statement for an optimal control problem was given previously The setup file The following list of MACSYMA variable names must be loaded into the MACSYMA memory stack These variables must be loaded into the file problem mac Standard MACSYMA syntax must be followed when these expressions are set up See the MACSYMA Reference Manual ref 19 for details stlist list of state variable names ctlist list of control variable names phi scalar cost expression that is a function of states at events only thi scalar cost expression that is a function of time at events only ellist list of integral cost terms corresponds to L in the performance index delist list of differential equations corresponds to f in the problem statement psilist list of boundary conditions that are a function of states at events only each term in p
67. to the states in stlist defined by the user in the setup file The costates labeled lambdal lambda2 lambdan in figure A4 correspond to the costates in the same order as the states The controls labeled ul u2 un in figure A5 correspond to the controls listed in ctlist defined by the user in the setup file Also the second partial derivative of the Hamiltonian H with respect to u is denoted by H uu in figure AT Finally the VI OTS makes no provisions for units on the plots because the units will change from problem to problem After the last plot the VTOTS is finished the user is then returned to the MATLAB prompt Note that in figure A3 all boundary conditions specified in psilist are satisfied Also the z axis on the graphs runs from 0 to 20 because the scale variable is used in vtotsinfo m Without that variable the z axis would run from 0 to 1 29 xl x3 150 100 50 15 10 xl vs time T x2 vs time T 10 time X3 vs time T 10 time 100 80 4 60 J N gt 40 4 20 F 4 l 0 l l l 15 20 0 5 10 15 20 time 10 x4 vs time m gt 15 20 0 5 10 15 20 time Figure A3 State histories 30 6 924 x10 24 ambdal vs time 0 746 lambda2 vs time E E 3 462 0 373 E E 0 l l l 0 l l l 0 5 10 15 20 0 5 10 15 20 time time 2 lambda3 vs time 4 lambda4 vs time
68. tor Hamiltonian Constrained double integrator eigenvalues of Huu Control constrained problem state histories Control constrained problem costate histories Control constrained problem control history State histories for two stage rocket problem Costate histories for two stage rocket problem Control history for two stage rocket problem iv 16 21 30 31 31 32 33 36 36 37 37 38 40 40 41 41 42 44 45 45 49 50 50 Abstract This report briefly describes the algorithms that comprise the Variational Trajectory Optimization Tool Set VTOTS package The VTOTS ts a software package for solving nonlinear constrained optimal control problems from a wide range of engineering and sci entific disciplines The VTOTS package was specifically designed to minimize the amount of user programming in fact for problems that may be expressed in terms of analytical functions the user needs only to define the problem in terms of symbolic variables This verston of the VT OTS does not support tabular data thus problems must be expressed in terms of analytical functions The VTOTS package consists of two methods for solving nonlinear optimal control prob lems a time domain finite element algorithm and a multiple shoot ing algorithm These two algorithms under the VTOTS package may be run independently or jointly The finite element algorithm generates approaimate solutions whereas the shooting algorithm pro vides a more accurate
69. ugh convergence was obtained Commonly a state or control that 1s an angle assumes a value in the wrong quadrant no converge during a shooting run generally indicates that the initial estimate provided by the user 1s too far from the solution A warning during compilation that a do loop is not executed in alltsi f may be ignored This warning occurs whenever tsilist is empty Helpful Hints In this section helpful hints are suggested for obtaining a solution to an optimal control problem It is assumed that the plant mex4 file is bug free and the name list file is complete 1 10 11 12 A finite element solution is almost always easier to obtain than a shooting solution therefore start with finite elements When using finite elements start with a small number of elements and increase in general the initial estimate does not need to be as accurate for a small number of elements as for larger numbers of elements When increasing the number of elements it is not necessary to increase the elements in each phase Avoid the use of numbers in the MACSYMA setup file Define these constants as variables in the name list Make sure that namcom nml is filled in properly in double precision A name list that is not filled in could lead to a singularity in the Jacobian Avoid the use of variables starting with 1 J k 1 m or n See the example in the section entitled State Constrained Dou
70. vol 24 July 1988 pp 499 506 Zhu Jian Min and Lu Yong Zai Hierarchical Strategy for Non Linear Optimal Control Systems Via STWS Approach Int J Control vol 47 no 6 June 1988 pp 1837 1848 Prenter P M Splines and Variational Methods John Wiley amp Sons Inc c 1975 Kelley C T and Sachs E W Quasi Newton Methods and Unconstrained Optimal Control Problems STAM J Control amp Optim vol 25 no 6 Nov 1987 pp 1503 1516 Bryson Arthur E Jr and Ho Yu Chi Applied Optimal Control Revised printing Hemisphere Publ Corp c 1975 Kuo Chung Feng and Kuo Chen Yuan Improved Gradient Type Algorithms for Zero Terminal Gradient Optimal Control Problems J Dyn Syst Meas amp Control vol 109 Dec 1987 pp 355 362 Gruver W and Sachs E Algorithmic Methods in Optimal Control Pitman Publ 1981 Oberle H J Numerical Treatment of Minimax Optimal Control Problems With Application to the Reentry Flight Path Problem J Astronaut Sci vol 36 nos 1 2 Jan June 1988 pp 159 178 Pesch Hans Josef Real Time Computation of Feedback Controls for Constrained Optimal Control Problems Part 1 Neighbouring Extremals Opt Control Appl amp Methods vol 10 no 2 Apr June 1989 pp 129 145 Pesch Hans Josef Real Time Computation of Feedback Controls for Constrained Optimal Control Problems Part 2 A Correction Method Based on Multiple Shooting Opt Control
71. y coordinate is 100 There is no integral cost and no constraints are imposed In order to demonstrate both the finite element algorithm and the shooting algorithm the problem is scaled so that the phase has a duration of one as required by the shooting algorithm The differential equations are multiplied by the final time to achieve the scaling See the section entitled Time Scaling Several constants are used in this problem the magnitude of the acceleration a the final time and the specified initial and final conditions on the states These constants are assigned values in the file namcom nml and can be changed between VTOTS runs without repeating the MACSYMA process For this problem the MACSYMA input file is as follows This is the fixed time trajectory optimization problem Section 2 4 Bryson and Ho stlist x y u vl ctlist betal glist L 11 glist L l1l ellist 0 phi u 1 2 thi 0 psilist x 1 1 x0 y 1 1 yo u 1 1 u0 v 1 1 vO y 1 2 yf v 1 2 vf 24 tsilist thyme 1 1 delist timtu tim v a tim cos beta a tim sin beta namcom x0 y0 u0 v0 yf vf a tim The name list file namcom nml is namcom XO 0 0d 00 YO 0 0d 00 UO 0 0d 00 VO 0 0d 00 YF 100 0d 00 VF 0 0d 00 A 1 12397d 00 TIM 20 0d 00 end The name list starts with a dollar sign in the second column and no data are entered in the first column The MACSYMA scripts produce this file w
72. ynamic equations are Tyac G1 sp 43 h V sin y Tfcose D usiny V m r2 a Saat ae V H E mV ro rV EY where T Tac Aep dyac 18 the thrust in a vacuum Ae is the nozzle exit area p is the pressure Isp is the specific impulse g is the acceleration due to gravity at sea level u is used here as the Earth s gravitational constant and r is the distance from the center of the Earth where Re h is the radius of the Earth The drag D and the lift L are composed of axial and normal components _ 1 2 Fa qSCa Fy qSC Na D Fy sina Fa cosa L Fy cosa Fa sin o where Fa and C are the axial force and coefficient Fy and Cy are the normal force and coefficient p is the density S is the reference area and q is the dynamic pressure The performance index is J mle and the final time ty is free The initial conditions specified are m 0 mo 890 149 09 kg h 0 0 m V 0 20 0 m sec and 0 1 57 rad The drop mass of the booster is 29920 kg The final velocity and altitude are V ty Vy 7854 m sec and h ty hy 148011 1 m Other constant values are listed in the name list file From a numerical standpoint all variables should be of the same order of magnitude Therefore the equations have been nondimensionalized by defining h hfhy V V Vr 0 xl VS time 10 x2 vs time 5L 8r jJ 6 L 4 x 10H 7 8 4L J 15 F 4 2L 20 0 0 5 10 0 5 10 time time Figure B

Download Pdf Manuals

Related Search

Related Contents

  BAM 50-51  取扱説明書 - Panasonic    Mode d`emploi  Grundig GMN 3540 hair dryer  Télécharger manuel et vues éclatées - Solano  Keor T 80-120 kVA  Genius G-Shot HD575T  ASUS A55BM-E T8564 User's Manual  

Copyright © All rights reserved. Failed to retrieve file