Variation Evolving for Optimal Control Computation, a Compact Way

1 Variation Evolv ing fo r Optimal Control Computation, a Compact W ay  Sheng ZHANG , Jiang-Tao H UANG , Kai-Feng HE , and Fei LIAO 1 Abstract A co mpact version of the variation evolving m ethod (VEM) is developed in the primal variable space for o ptimal control computation. Follo wing the idea that o riginates from the Lyapunov continuou s -time d ynamics stabil ity theo ry in t he control field, t he opti mal solu tion is analog ized to the stable equilibrium p oint of a d ynamic s ystem a nd ob tained asymptotically through the variation motion. With t he i ntroduction of a virtual di mension, namely th e variation time, the evolution partial differential equatio n (EPDE) , which seek s the opti mal sol ution with a theoretical g uarantee, is d eveloped for the op timal control pr oblem (OCP) with free ter minal states, and the equivalent optimality conditio ns with no employment of co states are established in the primal space. These conditions show that the optimal feedback control law is generally not an alyticall y available because the optimal control is relate d to the future state s. Since the derived EPDE is suitable to be computed with the semi-discrete method in the field of PDE numerical calcu latio n, the optimal solution may be obtained by solving the resulting fi nite-di mensional initial-value pro blem (IVP ). Keywords Optimal co ntrol • Lyapuno v d ynamics stabili ty • E volution p artial differential equation • I nitial-value problem • Costate-free opti mality condition AMS Cla ssification 49 M05 • 49M29 • 35B 35 1 Introduction Optimal con trol is a discipline that studies the opti mization on d ynamic systems. I t is clo sely related to the e ngineering and has been widely studied [1]. Because of the ir co mplexity, op timal control prob lems (OCPs) are usually solv ed numerical ly . T he prevailing n umerical method s may be divided into t wo types, namely, the d irect method s and the indirect methods [2]. The d irect methods discr etize the control or/and state variable s to obtain a nonli near progra mming (N LP) pr oblem [ 3-5] . T hese methods are easy to apply, whereas the results obtained are usually suboptimal [6]. T he indirect methods transfor m the OCP to a boundary-value problem (B VP) through the o ptimality con ditions [2 , 7] . T hey may produce more precise results, b ut often suffe r from significant numerical difficult y due to the ill -co nditioning of t he Ha miltonian dynamics; that is, the stability of t he costate dynamics is adverse to t hat of the state d ynamics [8], and this makes t he co mputation di fficult without a good initial guess [9]. ____________ __________ 1 The authors are with the Computational Aerodynamics Insti tution, China Aerod ynamics Research and Develop ment Center, Mianyang, 6210 00, China. (e- mail: zszhangshengzs@ hotmail.com). 2 Practical effective ways to ad dress the difficulty may be e mploying the evolutionary algorithms and metaheuristics [9], or combining the homotop y approach [ 10 , 11 ]. Researchers, from a philosophical perspective o f unity, try to uncover the connection between the direct and the indirect methods. Po sitive studies include the direct colloca tion method [6], the Runge-Kutta discretization m ethod [ 12 ], and t he pseudo-spectral (P S) method [ 13 - 15 ], which is popular in recent decades . A s such met hods blend the two t ypes of methods in a dualization view, the y inherit advantage s of both types and blur their differe nce [ 16 ] . Theories in the dynamics and control field often enlighten strategies for the numerical optimization. The notion of differential flatness, arising f ro m the non-linear control theor y, is advanced for the optimal control co mputation since its proposal in Ref. [ 17 ]. With the n o n-linear coordinate transformation to reduce th e variables , the OC Ps are simplified and the s olutions may be generated rapidly [ 18 - 20 ]. In the field of NLP, which is the static counterp art to the OCP, the dynamic m ethods based on the continuous-time ordinary differential eq uation (ODE) ar e prosperously de veloped [ 21 , 22 ]. With these met hods, NLP p roblems are tra nsformed to the initial-value proble ms ( IVP s) to be solved . T hrough the control parameterization technique , t he OCP may be discretized and then solved by the d ynamic method for NLP [ 23 ] . Interestingl y, add ressing the OCPs directl y i n a co ntinuous d ynamical manner will result in eq uations with the partial dif ferential equation (PDE) form, which determines the evolution dynamics of v ariab les . In Ref. [ 24 ] , f rom t he view o f “d ifferential distance” , which might be m ore accurate ly called “variational dista nce” , the evolution P DE (EPDE) for OCP s w ith M ay er performance index is pr opos ed through the “method of gradient” . In Ref. [ 25 ], a variation evolving method (VEM), inspired by the Lyapunov continuo us-time dynamics stabilit y theo ry, i s developed for the optimal control computation, and the EPDE, upo n a reco nstructed unconstrained functio nal, is d erived from the vie wpoint of variation motion in typical O CPs. Usi ng the well-kno wn semi -discrete method f or PDE numerical calculation [ 26 ], the finite-dimensional IVP s obtained may be solved with the mature ODE integration methods . The VEM also synthesizes the direct and indirect m ethods, but fro m a new sta ndpoint. B ecause the extre mum of t he OCP is theoretically guaranteed the eq uilibrium point of the d educed dynamic system, the op timal solution will be gradua lly appro ached. However, in the work of Re f. [ 25 ], b esides the states a nd the controls, the costates ar e also introduced, which increases the complexity o f the co mputation. In this p aper, a compact version of the VEM that uses onl y the pr imal variab les is devel oped and the corresp onding EPDE is formulated. Bo th the VEM and th e “m ethod of g radient” in [ 24 ] take advantage of the variation motion for the solutio n o f the O CPs. However, herei n the more gen eral Bolza performance index is considered, w hich avoids the introduction of extra variab le w hen transforming to the Mayer form. Therefore, the EP DE, derived from a di fferent approach in the following, is more common. In particular, the costate-free optimality conditions are established with the primal variables , an d they characterize the fixed point of the PDE d ynamic system. More importantly, under the frame of t he infinite-dimensio nal Lyapunov principle , the theoretical con vergence to the o ptimal solutio n is proved . Throughout the paper, our work is built upon the assumption that the solution for the optimizatio n pro blem exists. We d o not 3 describe the existing condition s for the purpose of brevity. Relevant studies may be found in [ 27 , 28 ] and references therein. In the following, f irst preli minaries that state the inf inite -dimensional Lyapunov theor y and the pr inciple of the VEM ar e presen ted. Then the co mpact VEM to solve t he OCP s with free ter minal states is d eveloped. During this co urse, the costate-free optimality conditions are derived , an d are proved equivalent to the class ic conditions. After that, illustrative examples are solved to verify th e effectiveness of the method. Besides , in -depth comments are presented before concluding the p aper. 2 Preliminaries 2. 1 Infinite-Dimensional Lyapunov Theory The VEM is a newly develop ed method for op timal solutio ns. It is enlightened fro m t he inverse consideration o f the Lyapunov continuous-time dynamics stab ility theory in the control field [ 29 ]. As the theoretical foundation of this method, the in finite-dimensio nal Lyapunov principle is introduced as follows. Definition 2. 1 Consider an infinite-dimensional d ynamic syste m described by ( , ) ( , , ) xt x t    y f y p , (1) d ( ) ( , ) d t t  p Γ yp , (2) where t is t he t ime, x  is t he i ndependent variable , ( ) ( ) ( ) n x x x   y is the function vector o f x , m   p is a state vector, : ( ) ( ) n xx    f and : ( ) m x  Γ are vector f unctionals, () x is a certain function set an d is a certain nu mber set . If ˆ ˆ ( ( ), ) ( ( ) ) x x   yp satisfies   ˆ ˆ ( ), , xx  f y p 0 and   ˆ ˆ ( ), x  Γ y p 0 , the n ˆ ˆ ( ( ), ) x yp is called an equilibrium sol ution. Note that De finition 2. 1 allo ws f to be a functional, a n exte nsion of t he function type used in Ref. [ 25 ]. Mathematically, an equilibrium solution o f the dynamic system (1) and (2) is the fixed point of their right functional s. Definition 2. 2 The equilibrium solu tion ˆ ˆ ( ( ), ) x yp is an as ymptotically stab le equilibri um solut ion in ( ( ) ) x  if for an y initial conditions 0 ( , ) ( ) ( ) t x t x x    yy and 0 () t t    pp , there is ˆ ˆ l im ( ( , ), ( )) ( ( ), ) 0 t x t t x    y p y p , where 1/ 2 1/ 2 22 11 ( ( ), ) : m ax sup , nm ii x ii x y p                     yp is the defined supr emum norm. 4 With the same proo f in Ref. [ 25 ], the following lemma may b e established . Lemma 2. 1 F or the infinite -dimensional dynamic system (1) and (2) , if there exists a co ntinuously differentia ble function al : ( ) V x  such that i)   ˆ ˆ ( ), V x c  yp and   ( ), V x c  yp in ˆ ˆ ( ( ) ) / { ( ( ), )} xx  y p , ii)   ( ), 0 V x  yp in ( ( ) ) x  and   ( ), 0 V x  yp in ˆ ˆ ( ( ) ) / { ( ( ), )} xx  y p , where c is a constan t, then ˆ ˆ ( ( ), ) ( ( ), ) xx  y p y p is an asymptotically stable equ ilibrium solution in ( ( ) ) x  . Lemma 2. 1 is a generalizatio n of the Lyapuno v stability theore m, without requiri ng V vanishing at the eq uilibrium. W hen employing Lemma 2. 1 in th e following, note that th e variables are functions of t (not x ) and t he time indicates th e variation time  , a virtual dimension i ntroduced to describ e the variation evolution pr ocess. 2. 2 Principle of VEM The VEM analogizes the optimal solution of an OCP to the asymptotically stable equilibriu m p oint of an in finite -dimensional dynamic system, and seeks such dynamics to minimize a specific performance index that acts as the Lyapunov functional. To implement the idea, the variation time  is introduced to describe the process that a variable evolves to th e optimal solution under the d ynamics gover ned by the variation d ynamic evolutio n equations, which may be presented in the for m of the EP DE and the evolution differential equat ion (EDE) [ 25 ]. Figure 1 illustrates the variation evolutio n of the v ariable () xt in the VEM to solve the OCP. Through the variation motion, the initial guess of () xt will evolve to the opti mal solution, and the optimality conditions will be gradually achieved. Fig. 1 The illustr ation of the variable evolution along the variation time  in the VEM [ 25 ] For example, consider the following calculus- of -variations problem, which may be regarded as OCPs w ith integrator dynamics. 5 Problem 1 For the follo wing functional depending on the va riable vector () n t  y   0 ( ), ( ) , d f t t J F t t t t   yy , (3 ) where t  is th e time. The elements of y belong to 2 0 [ , ] f C t t , w hich denotes th e set of variables w ith continuo us second-order derivatives . The fun ction : nn F    and its first-order and second-order partial derivatives are continuous with respect to y , its time deriva tive d d t  y y and t . 0 t and f t are th e fixed initial and term inal times respecti vely, and the boundary conditions are prescribed as 00 () t  yy and () ff t  yy . Find the optimal solution ˆ y that minimizes J , i.e., ˆ arg min( ) J  y . (4 ) Following the idea of d ynamic evolution to reduce some performance index , we anticipate that any initial guess of () t y , whose elements belong to 2 0 [ , ] f C t t , w ill evolve to the m inimu m solution through the variation moti on . Like the decrease of the Lyapunov function for a stable d ynamic s ystem, if J in Eq . (3) dec reases with re spect to the variatio n time  , i.e., 0 J    , we may eventually obtain the optimal solution. Differentiating Eq. (3) w ith respect to  (even if  does n ot ex plicitl y exist) produces 0 0 0 TT T TT d d ( ) d d f f f t t t t tt J F F t F F F F t t                               yy y y y y yy y y y , (5) where the column vecto rs F F    y y and F F    y y are the shorthand notation s of partial derivatives . T he superscript “ T ” denotes the transpo se operator. “ f t ” and “ 0 t ” m ean “ evaluated at f t ” and “ evaluated at 0 t ”, respectively. From an initial g uess of () t y that satisfies the boundar y conditions at 0 t and f t , then by enforcing 0 J    , we may set t hat 0 d ( ) , ] , [ d f F F t t t t          yy y K , (6) where K is an nn  dimensional positive-definite gain matrix. T he va riation dynamic evolution equation (6) may be considered from the view of the PDE formulation, by r ep lacing the variation operation “  ” and the differential operator “ d ” with the partial differential operator “  ” as ( , ) () t FF t           yy y K . (7) 6 In addition, the i nitial condition is 0 ( , ) ( ) t t     yy , (8) where 2 0 ( ) [ , ] f t C t t  y is an arbitrary sol ution that satisfies the bo undary conditions. Theorem 2.1 Solving the IVP defined by Eq s . (7) and (8) , when    , () t y will s atisfy the optimality conditions of Problem 1. Proof For the infinite-di mensional dynamics governed by Eq . (7) , Eq . (3) may ac t as the Lyapunov functio nal in 2 0 [ , ] f C t t . According to Lemma 2. 1, the minimum solutio n of Problem 1 is an as ymptotically stable equilibrium so lution. Thus, w hen    , () t y will meet the opti mality condition , namel y, the Euler-Lagrange equa tion [ 30 , 31 ] d () d FF t  yy 0 . (9) QED. □ For the EPDE (7), its right f unction depends only on the time t and is suitable to be solved with the semi-discrete method in the field of PDE n umerical calculation . With the discretization along the n or mal t ime dimension, Eq s . (7) and ( 8) are tran sformed to an IVP with finite states; then t he mature ODE numerical inte gration methods may be used to get the opti mal solution. Note that the resulting IVP is defined with respect to the variation time  , not the nor mal time t . In the previous work [ 25 ] , a demonstrative example is solved t o verif y the results. 3 The Compact V EM 3.1 Pro blem Definition In this paper, we consider the following clas s of OCP that is defined as Problem 2 Consider the perfor mance index of Bolza form   0 ( ( ), ) ( ), ( ), d f t ff t J t t L t t t t    x x u , ( 10 ) subject to the d ynamic equation ( , ) t  x f x, u , ( 11 ) where t  is the time. n  x is the state vector , with its ele ments belonging to 2 0 [ , ] f C t t . m  u is the control vec tor, with its ele ments belongi ng to 1 0 [ , ] f C t t . T he functio n : nm L    and its first -order partial der ivatives are continuous with respect to x , u and t . The function : m   and its first-order and second-order partial derivatives are continuous with 7 respect to x and t . The vector function : n m n    f and its first-order partial derivatives are continu ous and L ip schitz in x , u and t . The initial time 0 t is fixed and the ter minal time f t is free. The initial bou ndary conditions are p rescribed as 00 () t  xx , ( 12 ) and the terminal states ar e free. Find the o ptimal solution ˆ ˆ ( , ) xu that minimizes J , i.e., ˆ ˆ ( , ) arg m in( ) J  xu . ( 13 ) 3.2 Deriv ation of EPDE In Ref. [ 25 ], the problem is ad dressed by constructing an unconstrained functional proble m that has the same extrem um. This operation may b e practical bu t it introduces the costate vector, wh ich has the same dime nsion a s t he sta te vector . T o avoid the resulting complexity, here we will address P roblem 2 directly in the primal space, as we treated Problem 1 in the p receding . Consider the p roblem within the feasible solution domain o , in which any sol ution satis fies Eqs. ( 11 ) and ( 12 ). First, w e transform the ge neral Bolza perfor mance index ( 10 ) to the eq uivalent Lagrange t ype, i.e.,         0 00 0 T 00 T 00 ( ( ), ) ( ), ( ), d ( ( ), ) ( , ) d ( ), ( ), d ( ( ), ) ( , ) ( , , ) d f ff f t ff t tt t tt t t t J t t L t t t t t t t t L t t t t t t t L t t                    x x x x u x f x, u x u x f x, u x u , ( 14 ) where t  and  x are the p artial derivatives. No te that the term 00 ( ( ), ) tt  x is co nstant and may be neglected in the performance index. Similarly, we differentiate Eq. ( 14 ) with respect to the variation time  to obtain 0 T T T T T T T ( ) ( ) ( ) d f f t f t t u t t t J L L L t                                x x xx x x x x u xu f f f f , ( 15 ) where t  x and  xx are th e s eco nd-order partial d erivatives in th e form of (column) vector and matrix respectively, and x f and u f are the Jacob i an matrixes. Different from the calc ulus- of -variations pro blems,   x and   u are related because the profiles of x are determined b y u . For the solutions in o , they need to satisf y the following variatio n equation          xu x x u ff , ( 16 ) wi th the i nitial condition 0 t    x 0 . Note th at x f and u f are tim e-depende nt m atrix es linearized at the f eas ible solution () t x and () t u . Eq. ( 16 ) is linear with respect to the variables   x and   u , and has a zero initial value. T hus, it satisfies the su perposition principle, and its sol ution may b e explicitly expr essed. 8 Lemma 3.1 [ 32 ] For the linear time-va rying system ( ) ( ) t t   x A x B u , ( 17 ) with the initial value 0 () t  x0 , where n  x is the state vector, m  u is the control vector, and () t A and () t B are the righ t dimensional coefficient ma trix es , the solution is 0 ( ) ( , ) ( ) d t t t t s s s   x H u , ( 18 ) where ( , ) ts H is the nm  dimensiona l impulse respon se function t hat satisfies ( , ) ( ) ( , ) nm t s s t s ts st        Φ B H 0 , ( 19 ) and ( , ) ts Φ is the nn  dimensiona l state transition matrix fo r the system from time point s to time poin t t , which satisfies ( , ) ( ) ( , ) ts t t s t    Φ A Φ , ( 20 ) ( , ) ( , ) ( ) ts t s s s    Φ Φ A , ( 21 ) ( , ) nn tt   Φ 1 , ( 22 ) where nn  1 is the nn  dimensiona l identity matrix. According to Lemma 3.1, Eq . ( 16 ) has the solution 0 ( , ) ( ) d t o t t s s s        xu H , ( 23 ) where ( , ) o ts H is the impulse response function corresponding to the specific () t x f and () t u f . S ub stitute Eq . ( 23 ) into Eq . ( 15 ); there is     00 0 0 0 T T T T T T T T T T T T T T ( ) ( , ) ( ) d ( ) d ( ) ( , ) ( ) d d ( ) d f f ff f tt f t t o u t tt t t t f t t o u t t t t t J L L t s s s L t t L L t s s s t L t                                                              x x xx x x x x u x x xx x x x x u uu f f f H f uu f f f H f  . ( 24 ) By exchanging the ord er in the double integral, we may derive the following transfor mation as       00 0 T T T T T T T T ( ) ( ) ( ) ( ) ( ) ( ) ( , ) ( ) d d ( ) ( ) ( ) ( ) ( ) ( ) ( , ) d ( ) d f ff tt to tt tt to ts t t t t t L t t s s s t t t t t t L t t s t s s                       x xx x x x x xx x x x u f f H u f f H . ( 25 ) 9 To make the result clear er , the variable symbol s are changed as t   and st  , without altering the final result of Eq. ( 25 ), to produce         0 0 T T T T T T T T ( ) ( ) ( ) ( ) ( ) ( ) ( , ) d ( ) d ( ) ( ) ( ) ( ) ( ) ( ) ( , ) d ( ) d ff ff tt to ts tt to tt t t t t t L t t s t s s L t t t                            x xx x x x x xx x x x u f f H u f f H . ( 26 ) Thus, Eq . ( 24 ) may be refor mulated as       0 T T T T T T T ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( , ) d ( ) ( ) d ff f tt f t t o u t tt t J L L t L t t                              x x xx x x x x u u f f f H f . ( 27 ) Now to achieve 0 J    , we may derive the variation d ynamic evolution equa tion for the control u as       T T T T ( , ) ( ) ( ) ( ) ( ) ( ) ( ) d f t u o t t L t L                       u x x x xx x x u K f H f f , ( 28 ) and for the terminal ti me f t as   T f f f tt t t kL        x f , ( 29 ) where K is an mm  dimensional posi tive-definite gai n matrix and f t k is a scalar positive gain parameter . Re mark 3. 1 With the pro perty of the state tra nsition matrix, it may be obtained that     T T T T TT ( , ) ( ) ( ) ( ) ( ) ( ) d d = ( , ) ( ) d d = ( , ) ( ) ( ) ( ) f f t ot t t o t o f f u t t t t t t t                     x xx x x x xx H f f H Hf . ( 30 ) Th u s, Eqs. ( 27 ) and ( 28 ) may be equivalently refor mulated as follows.     0 T T T T ( ) ( ) ( , ) d ( ) ( , ) ( ) d ff f tt f t o f o f t tt t J L L t t t t L t t                   x x x u u f H H , ( 31 )   TT + ( , ) ( ) ( , ) ( ) d f t o f f o t L t t t t L           u x x u K H H . ( 32 ) Similarly, using the partial differential operator “  ” and the differential operator “ d ” to reformulate the variation dynamic evolution equatio ns ( 23 ) , ( 28 ) and ( 29 ), we may obtain the f ollowing EPDE and EDE as       0 T T T T ( , ) ( , ) d ( , ) ( , ) ( , ) ( ) ( ) ( ) ( ) ( ) ( ) d f t o t t u o t t s t s s t t L t L                                        u x x x xx x x u H x u K f H f f , ( 33 ) 10   T d d f f f tt t t kL       x f . ( 34 ) Put into this perspective ; the initial condit ions are 0 ( , ) ( ) ( , ) ( ) t t t t                  xx uu , ( 35 ) 0 f f tt    , ( 36 ) where () t x and () t u are arbitra ry feasible solution s that satisfy Eqs. ( 11 ) and ( 12 ), and f t is the initial guess of f t . 3.3 Establishment of Costa te-F ree Opti m ality C onditions The equilibriu m solution of the infinite -dimensio nal dynamic system given b y Eqs. ( 33 ) and ( 34 ) will satisfy     T T T T ( , ) ( ) ( ) ( ) ( ) ( ) ( ) d f t u o t t L t L                    u x x x xx x x f H f f 0 , ( 37 ) T ( ) ( ) ( ) ( ) 0 f t f f f L t t t t     x f . ( 38 ) Actually, Eqs. ( 37 ) and ( 38 ) are the first-order optimality condition s for Prob lem 2 without the employment o f costates. We will show that they are equivalent to the classic one s with costates [ 33 ], i.e., T HL      x x x λ λ f λ 0 , ( 39 ) T HL    u u u f λ 0 , ( 40 ) a nd the tran sversality conditio ns ( ) ( ) ff tt   x λ , ( 41 ) ( ) 0 f ft Ht   , ( 42 ) where T HL  λ f is the Ha miltonian and λ is the costate vector. Theorem 3.1 Fo r Problem 2, the optimal ity conditions given by E qs. ( 37 ) and ( 38 ) are equivalent to the optima lity con ditions given by Eqs. ( 39 ) - ( 42 ) . Proof W ith Eq. ( 19 ), Eq. ( 37 ) may be reformulated as     T T T T T ( , ) ( ) ( ) ( ) ( ) ( ) ( ) d f t u u o t t L t L                    u x x x x x x x ff Φ f f 0 , ( 43 ) where ( , ) o t  Φ is the state transitio n matrix related to ( , ) o t  H . Define a quantity () t γ as   T T T ( ) ( ) ( , ) ( ) ( ) ( ) ( ) ( ) ( ) d f t ot t t t t L                   x x x xx x x γ Φ ff . ( 44 ) 11 Then Eq. ( 43 ) is si mplified as T u L  u f γ 0 . ( 45 ) Obviously, fro m Eq. ( 44 ) , when f tt  , there is ( ) ( ) ff tt   x γ . (46) Differentiate () t γ with respect to t . In the process, we will use the Leibniz rule [ 34 ]       ( ) ( ) ( ) ( ) d d d ( , ) d ( ), ( ) ( ), ( ) ( , ) d d d d a t a t t b t b t h t h a t t a t h b t t b t h t t t t         , ( 47 ) and the proper ty of ( , ) o t  Φ (see Lemma 3.1) ( , ) ( , ) ( ) o o t tt t      x Φ Φ f , ( 48 ) ( , ) o n n tt   Φ 1 . ( 49 ) Then we have         T T T T T T T T T T T T d ( ) ( , ) ( ) ( ) ( ) ( ) ( ) ( ) d d ( ) ( , ) ( ) ( ) ( ) ( ) ( ) ( ) d () f f t t t o t t t ot t t L t L t L t t L Lt                                                   x xx x x xx x x x x x xx x x x x x x x xx x x xx γ f f f f Φ ff f Φ ff f γ .( 50 ) From Eqs. ( 50 ) and ( 46 ) , it is found that () t γ conforms to the same dynamics and boundary conditions as the costates () t λ . Thus, we can conclude th at ( ) ( ) tt  γ λ . Then Eq. ( 45 ) and Eq. ( 40 ) are identical . This m eans that Eq. ( 37 ) is equivalent to Eqs. ( 39 ) , ( 40 ) and ( 41 ) . On the o th er ha nd, w ith E q. ( 41 ) , it is e asy to s ho w t hat Eq s. ( 38 ) a nd ( 42 ) a re t he s a me. □ According to the above ana lysis, the costate-free Hamiltonian may be presented as     T T T T ( ) ( , ) ( ) ( ) ( ) ( ) ( ) ( ) d f t ot t H L t t L                    x x x xx x x f Φ ff . ( 51 ) Then the EPDE for the co ntrol variables is simply refor mulated as ( , ) t H      u u K . ( 52 ) An interesting trial on the costate -free optimality co nditions is to solve the optimal co ntrol from these theoretical expressions. B y investigating Eq. ( 37 ), it is found th at the optimal control is related to the future state s. This means generally the optimal feedback control la w i n the analytic form can not be solved. However, for the linear quadr atic r eg ulator (LQR) prob lem with the state equation given b y Eq. ( 17 ) and the per formance index given by 12 0 TT 1 ( ) d 2 t Jt    x Qx u Ru , ( 53 ) where Q and R are positive-de finite matrixes, the op timal contro l law may be derived as 1 T T ( ) ( ) ( , ) ( ) d t t t t        u R B Φ Qx . ( 54 ) subject to 1 T T ( ) ( ) ( ) ( , ) ( ) d t t t t t         x A x B R B Φ Qx , ( 55 ) and the closed-for m analytic solution may be obtained upon the follo wing identity T ( , ) ( ) d ( ) ( ) t t t t       Φ Qx P x , ( 56 ) where the matrix () t P is the solutio n of the Riccati di fferential eq uation T 1 T , ( ) nn          P A P PA P BR B P Q P 0 . ( 57 ) This may be verified by eli minating () t x in Eq. ( 56 ) and then imple menting the differen tiation . When the state equation ( 17 ) is time-invariant, then the clo sed- form analytic solution may be obtained upo n the following identit y T () ( ) d ( ) t t et       A Qx Px , ( 58 ) where the matrix P is the solutio n of the Riccati alg ebr ai c equation T 1 T nn       A P PA PB R B P Q 0 . ( 59 ) After proving the equivalence of op timality conditions, now the variables’ evolving direction using the VEM is easy to determine. Theorem 3.2 Solving the IVP defined by Eqs . ( 33 )-( 36 ) , when    , ( , , ) f t xu will satisfy the optimality cond itions of Pro blem 2. Proof Upon the i nfinite-dime nsional dynamics governed by Eqs. ( 33 ) and ( 34 ) , w ith the in itial conditions of Eq s . ( 35 ) and ( 36 ), we have 0 J    and 0 J    when J reaches its minimum within t he feasible solution d omain o . That is, the functional ( 10 ) is the Lyapunov functional in o . According to Lemma 2. 1, the minimum solution o f Eq. ( 10 ) is an asymptoticall y stable equilibriu m solution. Thus, the solution o f the IVP will as ymptotically meet the o ptimality conditions, namely , Eqs. ( 37 ) and ( 38 ). □ 13 Note that althou gh only the first-ord er op timality conditio ns ar e explicitl y guaranteed, the solution of the IVP, defined by E qs. ( 33 )-( 36 ), will not halt at a maximum or a sad dle point of Pr oblem 2 (unle ss they are the initial guesses), since t hose solutions are not asymptotically stable eq uilibrium solution s. Presume that we alrea dy ha ve a feasible i nitial solution of () t x and () t u ; then theoretically Eqs. ( 33 ) and ( 34 ) may be used to obt ain the optimal solution of Prob lem 2 . Recall Fig. 1, Eqs. ( 33 ) and ( 34 ) realize the anticipated variable evolution a long the variation time  . T he initial co nditions of ( , ) t  x and ( , ) t  u at 0   belong to the feasible solutio n do main and their value at    is op timal for the OCP. The right part of the EPDE ( 33 ) is also only a vector function of time t . T hus we may app ly the semi-discrete method to discretize it along the normal time dimension an d further use ODE integration methods to get the numerical solutio n. Th e OCP d efined in Problem 2 has a free ter minal time f t . When f t is fixed , the evolut ion equatio ns for the variables x and u remain the same, while the equation re garding the terminal time f t is no longer necessar y. 3.4 Co mputation of Impulse Response F unction In em plo ying the EPDE ( 33 ) t o get the solution, one important issue is to determine the impulse response function ( , ) o ts H , w hich is related to the current solution of () t x and () t u . We have ( , ) o n m t s when t s   H0 , ( 60 ) and for ts  , there is d ( , ) ( ) ( , ), ( , ) ( ) d o o o t s t t s t s s wh en t s t    xu H f H H f . ( 61 ) When applying the semi-discrete method for the initial-value problem of the PDE system, the numerical computation of ( , ) o ts H , along the d iscretization of the ti me coordinate , i.e. , ( 1 , 2, ..., ) i t i N  with N being the number of d iscretizatio n time points, is performed as follo ws. Given a fixed time point j t , for an y ij tt  , w e may use an implicit integration format, which brings higher preci si on and a larger stability region t han those o f the explicit format, to get     1 1 1 1 21 2 1 2 2 1 1 1 0.5 ( ) ( ) ( , ) ( ) ( , ) ( , ) ( , ) ( , ) ( , ) 0.5 ( ) ( ) ( , ) ( ) ( , ) ... ... .. ( , ) ( , ) j j j o j j j o j j o j j o j j o j j o j j j j j o j j j o j j o N j o N j t t t t t t t t t t t t t t t t t t t t t t t t t t t t                                            xx xx f H f H HH HH f H f H HH   1 1 1 . 0.5 ( ) ( ) ( , ) ( ) ( , ) N N N o N j N o N j t t t t t t t t             xx f H f H . ( 62 ) Then we may deri ve that 14     1 1 1 1 1 1 1 1 2 3 4 0.5( ) ( ) ( , ) ... o j j j n n j j j j j nn o j j nn t t t tt                   x H f1 0 Mtx C C C C H 0 . ( 63 ) where 1 2 1 ( , ) ( , ) ... ( , ) o o j j o j j j o N j tt tt tt            H H H Mtx H , 1 ... ... nn n n n n j n n n n             1 11 C 11 , 3 ... nn n n n n j nn n n n n            1 11 C 1 11 , 1 21 2 1 0.5( ) 0.5( ) ... 0.5( ) j j n n j j n n j N N n n tt tt tt                1 1 C 1 , and 1 2 4 () () ... () j j j N t t t          x x x f f C f . Note that 1 1 1 1 1 2 3 4 () j j j j      C C C C is a triangular m atrix and it s inverse al ways exists f or a sufficie ntly small discretiza tion granularity. Therefore, the values of the impulse response function at different time point s, with respect to j t , are     1 1 1 1 1 1 1 1 1 2 3 4 ( , ) ... ... ( , ) 0.5( ) ( ) ( ) ( , ) ... ... ( , ) nn oj nn o j j j j j n j o j j j j j j nn o N j nn tt tt t t t t tt tt                                             xu 0 H 1 H f 1 f H 0 C C C C H 0 . ( 64 ) 4 Examples First, a linear e xample taken fro m X i e [ 35 ] is con sidered. Example 4. 1 Consider the following dynamic sys tem u  x Ax b , where 1 2 x x     x , 01 00     A , and 0 1     b . Find the solution that mini mizes the perfor mance index   0 T T 2 11 ( ) ( ) d 22 f t ff t J t t Ru t     x Fx x Qx , 15 with t he i nitial boundary conditio ns 0 1 () 1 t     x , where the in itial time 0 0 t  and the ter minal ti me 3 f t  are fixed. The w eighted matrixes are 10 02     F , 21 14     Q and 1 2 R  . In solving this exa mple using the VEM , the EPDE der ived wa s       0 T () T T ( ) T ( ) d ( ) ( ) ( ) d f t ts t t t t u e s s u K Ru e u                               A A b x b Fx b Q FA A F x Fb , with t he one-dimensional gain K set as 2 2 10 K   . T he initial conditio ns o f t he EP DE , with the cont rol input ( ) 0 u t  , were 1 () 1 t t       x . Using the se mi-discrete me thod, the time horizon 0 [ , ] f tt was discretized uniformly with 61 po ints. Th us , a dynamic system with 1 83 states wa s obtained and the OCP wa s transfor med to a finite-dimensional IVP . The ODE integrator “ ode45 ” in MAT LAB, with a default relative error to lerance of 1 × 10 -3 and a default absolute erro r tolerance of 1 × 10 -6 , wa s employed to solve th e IVP . Even for this simple example, the an alytic solution is not easy to obtain. For com pariso n, w e com p uted the optimal solution with GPOP S- II [ 36 ], a Radau PS method -based OCP solver . Fig ure s 2, 3 and 4 show the e volution proce ss of 1 () xt , 2 () xt and () ut solutions to their opti mal val ues, respectively. A t  = 300 s, the n umerical solutions are indistinguishable from the optim al, which shows the effectiveness of the VEM. Figure 5 plots the profile of the p erfor mance index value a gainst t he v ar iation ti me. This i ndex decli nes rapi dly at first and is close to t he minimum when  = 50 s. Then it approach es the minimum monotonica lly and is almost fixed at the op timal value. 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 3.5 4 t x 1 ( t ) The op ti m al s ol ut i o n N um eri cal s o l u ti o ns wi th the V EM  = 1.0s  = 1.8s  = 0s  = 0.4s  = 300s  = 20.7s Fig. 2 The evo lution of numerical solut ions of 1 x to the optimal so lution 16 0 0 .5 1 1.5 2 2 .5 3 -0 .5 0 0.5 1 t x 2 ( t ) The o p ti m al so l uti on N um er i c al s o l ut i o ns wi th th e V EM  = 0s  = 1.0s  = 0.4s  = 20.7s  = 1.8s  = 300s Fig. 3 The evo lution of numerical solut ions of 2 x to the optimal so lution -0 .5 0 0 .5 1 1.5 2 2.5 3 -6 -5 -4 -3 -2 -1 0 1 t u ( t ) The o p ti m al so l u ti on N um er i c al s o l uti ons wi th th e V E M  = 0.4s  = 0s  = 1.8s  = 300 s  = 20.7s  = 70.9s Fig. 4 The evo lution of numerical solut ions of u to the optimal so lution 0 50 100 150 200 250 300 0 10 20 30 40 50  (s) J M i ni m u m o f J Functi onal va l ue o f th e n um eri cal s o l uti on 0 2 4 6 0 20 40  (s) J Fig. 5 The approa ch to the minim um performance index 17 Now we consider a nonlinear example with free terminal time f t , the 2-dimensiona l homing missile p roble m ad apted f rom Hull [ 37 ]. Example 4. 2 Consider the pr oblem of a constant spee d missile intercepting a co nstant speed tar get moving in a s traight line. T he dynamic equatio ns are ( , ) u  x f x , where M x y        x and cos( ) cos( ) sin( ) sin( ) T T M M T T M M M VV VV u V            f . Her e, x and y are th e abscissa and ordinate o f th e relative po sition , respectively, M  is the azimuth of the missile, T  =30 deg is the azi muth of the tar get, M V =1000 m/s is t he constan t speed of the missile , T V =500 m/s is t he constant speed of the target , and u is the missile nor mal acceleration. To intercept the target and penalize too lar ge normal acceleration , the performance index to be minimized i s defined as 0 T2 11 ( ) ( ) d 22 f t ff t J t t R u t   x Fx , where the weighted matrixes are 2 2 1 10 0 0 0 2 10 0 0 0 0          F and 4 5 10 R   . T he initial b oundary conditions ar e 0 0 10000 m 5000m 0 deg M t x y                        and the terminal time f t is free. Be fore th e computation, the s tate and the control v ariables were scaled to improve the n umerical efficie ncy. In the specific form of the EPDE ( 33 ) and the EDE ( 34 ), the gain para meters K and f t k were set to be 6 1 .5 10   and 4 1 10   , respectivel y . The initial conditions for the evol ution e quations, i.e., 0 ( , ) ( , ) () f t ut t           x , were obtained b y numerical integratio n at the time horizon of [0, 2 5]s with t he co ntrol inp ut ( ) 0 u t  . We also discretized the time horizon 0 [ , ] f tt uniformly, with 51 points . Th erefore, an IVP w ith 205 states (i ncluding the ter minal ti me) wa s obtained. W e still emplo y ed “ ode 45 ” in M ATLAB for the numerical inte gration. In t he integrator setting, the d efault r elative er ror tolerance and the abso lute erro r tolerance were 1 × 10 -3 and 1 × 10 -6 , respec tively. For comparison , the opti mal solution w a s again co mputed with GPOPS - II. 18 Fig ure 6 gives the states curve in the xy relative position coordinate p lane , showing that the numerical res ults tend to the optimal sol ution o ver ti me. For the optimal so lution, the missile will intercep t the target w ith a fairly small position error. The control solutions are plotted in Fig. 7, and th e asymptotical approac h of th e numerical results to the optimal s o lution is demonstrated. In Fig. 8, the terminal time profile against the variation time  is plotted. The results of f t os cillate sharply at first and then gradually ap proach the o ptimal interception ti me, and the value only chan ges slightly after  = 40 s. At  = 3 00 s, we compute that f t = 23.52 s from the VEM , the sa me as the result from GPOPS - II . -5 00 0 0 5000 10000 0 2000 4000 6000 8000 10000 12000 x ( m ) y ( m ) The o p ti m al s ol u ti o n N um eri cal s o l uti ons wi th th e V E M  = 3.1s  = 2.1s  = 1.2s  = 0.6s  = 0s  = 300s Optimal t ermina l re lativ e coord intates Initial relativ e c oor dintates Fig. 6 The ev olution of numerical solut ions in the xy coordinate plane to the optimal solutio n 0 5 10 15 20 25 -1 0 0 10 20 30 40 50 60 70 t u ( t ) The op ti m al s ol ut i o n N um eri cal so l u ti o ns wi th the V EM  = 2.1s  = 1.2s  = 0.6s  = 0.3s  = 0s  = 300s Fig. 7 The evo lution of numerical solut ions of u to the optimal so lution 19 0 50 100 150 200 250 300 23 .4 23 .6 23 .8 24 24 .2 24 .4 24 .6 24 .8 25  (s) t f (s ) M i ni m u m decl i ne ti m e E vol vi ng p rofi l e o f t f -2 0 2 4 6 23 .5 24 24 .5 25  (s) t f (s ) Fig. 8 The evo lution profile of f t to the optimal intercept ion time 5 Further Comments In contrast to th e augmented EPDE (AEPDE) that introduced th e costates i n Ref. [ 25 ], this paper develops a compact version of the EPDE that uses the primal v ariab les only . Discussion between the tw o formulations is worthwhile. For convenience, the A EPDE is again presented. ( , ) H t t t H H t t t H                                                            x yy x u λ x f y K x λ f 0 , ( 65 ) wher e       x y λ u , T HL  λ f is th e Hamiltonian, and K is a positive-de finite g ain matrix. By comparison, the advantages of the EPDE ( 33 ) over the AE PDE ( 65 ) may be listed as follows. i) The costates are not included and the evaluation on the second-ord er derivatives of H is avoided, which will relieve the computation burden. ii) The solution of the EPDE is m ore capable of reaching a minimum, while the sol ution of the AEPDE may halt at a saddle. iii) Generall y, t he EPDE req uires the in tegration, and the differentiation, as displayed in the AEPDE , may be avoided. This is ad vantageous to red uce the numerical erro r in seekin g solutions. Regarding the d isadvantages , i) c urrently, the EPDE may address OCPs w ith f ree terminal states only . W he n the terminal states are co nstrained, which is co mmon in t he OCP for mulation, it is not d irectly applicab le. One may penalize the terminal boundary conditions in the per formance index. However , this is not a satisfactor y solution. ii) Mo reover, the EPDE requires the initial solution to b e feasible, which is also inflexible, e ven if a feasible solutio n for the problem focused herein may be generated by the high-precision numerical integration. 20 The mechanism t hat the classic iterative methods use to solv e the OCPs ma y be generally descr ibed as   1 ( ( ), ( )) ( ( ), ( )) ii t t t t   x u M x u , ( 66 ) where M represents so me mapping operator and i denotes the number of iteration . A v ital conce pt introduced in the VEM is the variation time  . Interpreted from t he view of the n umerical computation dimension, it may b e associated with the iterations in the iterative m ethod s. Namely, the numerical iterations are th e manifestation of the discrete dynamics along the virtual time dimension  . C onsider in this way; the VEM proposed in this paper may be regarded as a continuous version of the iterative gradient method [ 33 ], while the AEPDE that in cludes the seco nd-order Hessian derivative is a continuous realization o f the Newto n t ype iterati ve mechanism, in particular , with co nvergence theoretica lly guaranteed. Regarding the computation, through the semi-discrete m ethod, th e IVP s on the infinite-dimensional dynamics may b e solved with the mature ODE integration methods. Particularly d uring t he discretization, techniques such as the Legendre-Gauss (LG), Legendre-Gauss -Radau ( LGR) and Legendre -Gauss-Lobatto (LG L) discretizatio n, w hich have stronger appr oximation capacity, may be emplo yed to improve the accurac y and efficiency. In add ition, si nce the i ntegration can be achieved with a simp le analog circuit, t he EP DE may pro vide a promising way of opti mal control in engineering. 6 Conclusions A co mpact version of the variation evolvin g method (VEM) for the optimal control computation is developed in the primal space . This method introduces no e xtra variab les in transfor ming the optimal control pr oblems (OCPs) to the i nitial-value pr obl em s (IVPs). T he costate-free o ptimality conditions are established and t hey ar e proved to be equivalent to the c lassic optim alit y conditions. One d irect result found by analyzing the optimality condition s is that the analytic op timal feedback co ntrol law is generally not algebraica lly solvable b ecause the control is r elated to the future state s. Regarding the computatio n, with the mature ordinary differential eq uation (ODE) integration methods, the solution of the resulting IVP may reach the minimum of the OCP effectively . In particular, under the fram e of IVP s upon continu o us-time dynamics, the daunting task of searching for a reasonable step size and t he annoying oscillatio n pheno menon around the minimum, as occ urs in the d iscrete iterative numerical method, are eliminated. Ho wever, currently t he method can solve the OCP with free ter minal states onl y, and this restricts its application. Further studies will be carried out to address this issue. R EFERENCES 1. P e s c h , H . J . , P l a i l , M. : Th e m a x i m u m p r i n c i p l e o f o p t i m a l c o n t r o l : A h i s t o r y o f i n g e n i o u s i d ea s a n d m i s s e d o p p o r t u n i t i e s . Co n t r o l C y b e r n . 38 ( 4 ) , 973 - 995 ( 2 0 0 9 ) 2. B e t t s , J . T. : S u r v e y o f n u m e r i c a l m e t h od s f o r t r a j e c t o r y o p t i m i z a t i o n . J . G u i d . C o n t r o l D y n . 21 ( 2 ) , 1 9 3 - 206 ( 1 9 9 8 ) 21 3. L i n , Q . , L o x t o n , R . , T e o , K . L . : T h e c o n t r o l p a r a m e t e r i z a t i o n m e t h o d f o r n o n l i n e a r o p t i m a l c o n t r o l : A s u r v e y . J . I n d . M a n a g . O p t i m . 1 0 ( 1 ) , 275 - 309 ( 2 0 1 4 ) 4. H a r g r a v e s , C . , P a r i s , W. : D i r ec t t r a j e c t o r y o p t i m i z a t i o n u s i n g n o n l i n ea r p r o g r a m m i n g a n d c o l l o c a t i o n . J . G u i d . C o n t r o l D y n . 10 ( 4 ), 338 - 3 4 2 ( 1987 ) 5. M a r t i n , C . S . , C o n w a y , B . A. : A n e w n u m e r i c a l o p t i m i z a t i o n m e t h o d b a s e d o n T a y l o r s e r i e s . I n P r o c e e d i n g s o f 2 1 s t A A S / A I A A S p a c e f l i g h t M e c h a n i c s M e e t i n g , N e w O r l e a n s , 2 0 1 1 , A AS P a p e r 1 1 - 1 5 6 ( 2 0 1 1 ) 6. S t r y k , O . V . , B u l i r s c h , R. : D i r e c t a n d i n d i r e c t m e t h o d s f o r t r a j e c t o r y o p t i m i z a t i o n . A n n . O p e r . R e s . 37 ( 1 ), 357 - 373 ( 1 9 9 2 ) 7. P e n g , H . J . , Ga o , Q . , Wu , Z . G . , Z h o n g , W . X . : S y m p l e c t i c a p p r o a c h e s f o r s o l v i n g t w o -p o i n t b o u n d a r y - v a l u e p r o b l e m s . J . G u i d . C o n t r o l D yn . 35 ( 2 ), 653 - 658 ( 2 0 1 2 ) 8. Rao , A . V. : A s u r v e y of n u m e r i c a l m e t h o d s f o r o p t i m a l c o n t r o l . In P r o c e e d i n g s o f AA S / A I A A A s t r o d yn a m i c s Sp e c i a l i s t C o n f e r e n c e , P i t t s b u r g h , P A , 2 0 0 9 , A A S P a p e r 0 9 - 334 ( 2 0 0 9 ) 9. C o n wa y , B . A . : A s u r v e y o f m e t h o d s a v a i l a b l e f o r t h e n u m e r i c a l o p t i m i z a t i o n o f c o n t i n u o u s d yn a m i c s y s t e m s . J . O p t i m . T h e o r y A p p l . 1 5 2 , 271 - 3 0 6 ( 2 0 1 2 ) 10. B e r t r a n d , R . , E p e n o y , R . : Ne w s m o o t h i n g t e c h n i q u e s f o r s o l v i n g ba n g - b a n g o p t i m a l c o n t r o l p r o b l e m s - n u m e r i c a l r e s u l t s a n d s t a t i s t i c a l i n t e r p r e t a t i o n . O p t i m . C o n t r o l A p p l . M e t h . 2 3 ( 4 ) , 1 7 1 - 1 9 7 ( 2 0 0 2 ) 11. P a n , B . F . , Lu , P. , P a n , X. , Ma , Y . Y . : D o u b l e - h o m o t o p y m e t h o d fo r s o l v i n g o p t i m a l co n t r o l p r o b l e m s . J. Gu i d . Co n t r o l D y n . 3 9 ( 8 ) , 1 7 0 6 -1 7 2 0 ( 2 0 1 6 ) 12. H a g e r , W. W . : R u n g e - Ku t t a m e t h o d s i n o p t i m a l c o n t r o l a n d t h e t r a n s f o r m e d a d j o i n t s y s t e m . N u m e r . M a t h . 8 7 , 2 4 7 - 282 ( 2 0 0 0 ) 13. G a r g , D . , Pa t t e r s o n , M . A . , H a g e r , W. W . , R a o , A . V . , e t a l . : A Un i f i e d fr a m e w o r k fo r t h e n u m e r i c a l so l u t i o n of o p t i m a l c o n t r o l pr o b l e m s us i n g p s e u d o s p e c t r a l m e t h o d s . A u t o m a t i c a , 46 ( 11 ) , 1 8 4 3 - 1 8 5 1 ( 2010 ) 14. G o n g , Q . , R o s s , I . M. , Ka n g , W. , F a h r o o , F . : C o n n e c t i o n s be t w e e n t h e c o v e c t o r ma p p i n g th e o r e m a n d co n v e r g e n c e o f p s e u d o s p e c t r a l m e t h o d s f o r o p t i m a l c o n t r o l . C o m p u t . O p t i m . A ppl . 4 1 ( 3 ) , 307 - 3 3 5 ( 2 0 0 8 ) 15. H a g e r , W. W . , Li u , J. , M o h a p a t r a , S. , R a o , A. V . , Wa n g , X . S . : Co n v e r g e n c e r a t e f o r a G a u s s c o l l o c a t i o n me t h o d a p p l i e d to c o n s t r a i n e d op t i m a l c o n t r o l . S I A M J . C o n t r o l O p t i m . 5 6 , 1 3 8 6 - 1 4 1 1 ( 2 0 1 8 ) 16. R o s s , I. M ., Fa h r o o , F. : A p e r s p e c t i v e o n m e t h o d s f o r t r a j e c t o r y o p t i m i z a t i o n . I n P r o c e e d i n g o f A I A A / A A S A s t r o d y n a m i c s C o n f e r e n c e , M o n t e r e y , C A , 2 0 0 2 , A I A A P a p e r N o . 2 0 0 2 - 4727 ( 2 0 0 2 ) 17. F l i e s s , M. , Lé v i n e , J. , M a r t i n , P . , R o u c h o n , P. : F l a t n e s s a n d d e f e c t o f n o n l i n e a r s y s t e m s : In t r o d u c t o r y t h e o r y a n d ex a m p l e s . In t . J . C o n t r . 6 1 ( 6 ) , 1327 - 1 3 6 1 ( 1 9 9 5 ) 18. R o s s , I. M . , F a h r o o , F. : P s eu d o s p e c t r a l m e t h o d s f o r o p t i m a l m o t i o n p l a n n i n g o f d i f f e r e n t i a l l y f l a t s y s t e m s . I E E E T r a n s . A u t o m . C o n t r o l . 49 ( 8 ), 1410 - 1 4 1 3 ( 2 0 0 4 ) 19. M i l a m , M . B. : R e a l - t i m e o p t i m a l t r a j e c t o r y g e n e r a t i o n f o r c o n s t r a i n e d d y n a m i c a l s y s t e m s . P h . D . t h e s i s , C a l i f o r n i a In s t i t u t e o f T e c h n o l o g y, P a s a d e n a , C A ( 2003 ) 20. Z h a n g , S . , Z h a o , Q . , H u a n g , H . B . , T a n g , G . J .: Ra p i d p a t h p l a n n i n g f o r z e r o p r o p e l l a n t m a n e u v e r s b a s ed o n f l a t n e s s . J . A e r o s p a c e E n g . 3 0 ( 5 ) , 0 4 0 1 7 0 5 4 ( 2017 ) 21. S c h r o p p , J. : A d y n a m i c a l s y s t e m s a p p r o a c h t o c o n s t r a i n e d m i n i m i z a t i o n . N u m e r . M a t h . 21 , 5 3 7 - 551 ( 2000 ) 22. M e r t i k o p o u l o s , P . , S t a u d i g l , M. : O n t h e c o n v e r g e n c e o f g r a d i e n t - l i k e f l o w s w i t h n o i s y g r a d i e n t i n p u t . S I A M J . O p t i m . 28 , 163 - 1 9 7 ( 2 0 1 8 ) 23. S n y m a n , J . A . , F r a n g o s , C . , Y a v i n , Y. : P e n a l t y fu n c t i o n s o l u t i o n s t o op t i m a l c o n t r o l p r o b l e m s w i t h g e n e r a l co n s t r a i n t s v i a a d y n a m i c o p t i m i s a t i o n m e t h o d . C o m p u t . M a t h . A ppl . 23 ( 23 ), 47 - 55 ( 1992 ) 22 24. K e l l e y , H . J . : M e t h o d s o f g r a d i e n t s . In : L e i t m a n n , G. : ( e d . ) : O p t i m i z a t i o n T e c h n i q u e s , p p . 2 1 6 - 2 3 2 . Ac a d e m i c P r e s s , L o n d o n ( 1 9 6 2 ) 25. Z h a n g , S . , Y o n g , E . M . , Q i a n , W . Q . , He , K. F . : A v a r i a t i o n ev o l v i n g m e t h o d f o r o p t i m a l c o n t r o l c o m p u t a t i o n . J . Op t i m . T h e o r y Ap p l . 1 8 3 , 246 - 270 ( 2 0 1 9 ) 26. Z h a n g , H . X . , S h e n , M . Y. : C o m p u t a t i o n a l F l u i d D y n a m i c s — Fu n d a m e n t a l s a n d A p p l i c a t i o n s o f F i n i t e D i f f e r e n c e M e t h o d s . N a t i o n a l D e f e n s e In d u s t r y P r e s s , B e i j i n g ( 2003 ) 27. H a r t l , R . F . , S e t h i , S . P . , V i c k s o n , R . G. : A s u r v e y o f t h e m a x i m u m p r i n c i p l e s f o r o p t i m a l c o n t r o l p r o b l e m s w i t h s t a t e c o n s t r a i n t . S IA M R e v . 37 ( 2 ), 181 - 218 ( 1 9 9 5 ) 28. L o u , H . W . : A n a l y s i s o f t h e o p t i m a l r e l a x e d c o n t r o l t o a n o p t i m a l c o n t r o l p r o b l e m . A p p l . M a t h . O p t i m . 5 9 ( 1 ) , 75 - 97 ( 2 0 0 9 ) 29. K h a l i l , H . K . : N o n l i n e a r S y s t e m s , 3 rd e d n . P r e n t i c e H a l l , U p p e r S a d d l e R i v e r ( 2002 ) 30. Wu , D . G. : V a r i a t i o n M e t h o d . H i g h e r E d u c a t i o n P r e s s , B e i j i n g ( 1987 ) 31. C a s s e l , K . W. : V a r i a t i o n a l M e t h o d s w i t h A p p l i c a t i o n s i n S c i e n c e a n d E n g i n ee r i n g . C a m b r i d g e U n i v e r s i t y P r e s s , C a m b r i d g e ( 2013 ) 32. Z h e n , D . Z. : Li n e a r S y s t e m T h e o r y . T s i n g h u a U n i v e r s i t y P r e s s , B e i j i n g ( 2002 ) 33. B r y s o n , A . E . , Ho , Y . C. : A p p l i e d O p t i m a l C o n t r o l : O p t i m i z a t i o n , E s t i m a t i o n , a n d C o n t r o l . H e m i s p h e r e , W a s h i n g t o n , D C ( 1975 ) 34. W a n g , H . , L u o , J . S . , Z h u , J . M .: A d v a n c e d M a t h e m a t i c s . N a t i o n a l U n i v e r s i t y o f D e f e n s e T e c h n o l o g y P r e s s , C h a n g s h a ( 200 0 ) 35. X i e , X . S. : O p t i m a l C o n t r o l Th e o r y a n d A p p l i c a t i o n . Ts i n g h u a U n i v e r s i t y P r e s s , B e i j i n g ( 1986 ) 36. P a t t e r s o n , M . A . , R a o , A . V. : GP O P S - I I : AM A T L A B s o f t w a r e f o r s o l v i n g m u l t i p l e - p h a s e o p t i m a l c o n t r o l p r o b l e m s u s i n g hp - a d a p t i v e Ga u s s i a n q u a d r a t u r e c o l l o c a t i o n m e t h o d s a nd s p a r s e n o n l i n e a r p r o g r a m m i n g . A C M Tr a n s . M a t h . S o f t wa r e 41 ( 1 ), 1 - 37 ( 2014 ) 37. H u l l , D . G. : O p t i m a l C o n t r o l T h e o r y f o r Ap p l i c a t i o n s . S p r i n g e r , N e w Y o r k ( 2003 )