Solve the General Constrained Optimal Control Problem with Common Integration Method

1  Sheng Z HANG and Jin-Mei GAO (2017. 12 ) Abstract : Com putation of general state- and/or contro l-constrained Op timal Con trol Problems (OCPs) is difficu lt for v arious constraints, especially the intractable path constraint. For su ch problems, the theoretical convergence of numerical algorithms is us ually not guaranteed, and the r ight solution may not be suc cessfully obtained. W ith the recently proposed V ariation Evo lving Method (VEM), the evolution equations, which guarantee the convergence towards the optima l so lution in th eory even fo r the gene ral constr ained O CPs, ar e derived. In particular, the costate-free optimality conditions are established. Besides the analytic expressions of the costates and the Lagrange multipliers adjoining the ter m inal constraint, the integral equation that determines the Karush-Kuhn-Tucker (K KT) multiplier variable is also derived . Upon the work in this paper, the general constrained OCPs may be tr ansfor med to the Initial-value Problems (IVPs) to be solved, with common Ordinary Differential Eq uation (ODE) numerical integration methods. Key words : Optim al control, path constraint, dynamics stability , variation evolution , th e second evolution equation, costate-free optimality condition, Karush-Kuhn-Tucker (KKT) condition , initial-value problem. I. I NTRODUC TION Optimal control th eory ai ms to deter mine the inputs to a d ynamic syste m that optimize a specified performance index while satisfying constraint s on the mo tion of the system. It is clo sely related to en gineering and has b een widely studied [1] . Because of the co mplexity, Optimal Control Problems (OCP s) are usually sol ved with numerical methods. Various numerical meth od s are developed and generally t hey are divided into two classes, namel y, t he direct methods and the i ndirect methods [2] . The direct methods d iscretize the control or /and state variab les to o btain the Nonlinear Programming (N LP ) problem, for example, the widely-used direct shooting met hod [3] and the classic collocation method [4]. These met hods are easy to app ly, whereas the results obtained are usually suboptimal [5], and th e optimal may be infinitely approached. The in direct m ethods transform the OCP to a B oundary-value Pr oblem (BVP) through the o ptimality condition s. Typical met hods of this type incl ude the wel l - known indirect shooting meth od [2 ] and the novel symplectic met hod [6 ]. Although be more pr ecise, the indirect methods often su ffer from the significant numerical difficulty due to the ill -conditioning of the Ha miltonian dyna mics, that is, the stabilit y o f costates dynamics is ad verse to that of the state s d ynamics [7] . The recent develop ment, repr esentatively the Pseudo -spectral (PS) method [8], blends the t w o types o f methods, as it unifies the NLP and the BVP in a dualizati on vie w [9 ]. Such methods inherit the advantages of both types a nd blur their difference. Even if th ere are m any numerical m ethods available, the computation of general state- and/or control-constrained OCPs is s till a tough topic tod ay. T he complex i nequality p ath constrai nts result in the da unting optimality co nditions in t he classic optimal Th e au thors are with the Computa tional Ae rodynamics I nstitution, China Aerodynamics R esearch and Devel opment Center, Mianyang , 621000, China. (e -mail: zszhangshengz s1@outlook.com). Solve the General Constrained Optim al Control Problem with Common Integration Method 2 control theory [10]. The indirect multiple shooting method is once the mainstream approach for th e optimal solutions of constrained OCPs. Ho wever, the utilization of the complex optimality conditions with costates is very user-unfrie ndly. Moreo ver, determination of s uch op timality conditions is still not co mpletely solved [11 ] . As the dev elopment o f the direct methods, p eople prefer to use them because they are much more understandable. Also, the greatly improved efficiency has ma de th e direct m ethods popular nowadays. Ho wever, the theoretic convergence of such algorithms is not well guaranteed for constrained OCPs, even i f varieties of pro mising work has been finished [12] [13]. Theories in the control fie ld often enlighten strate gies for the optimal contro l computation, for example, the no n -linear va riable transformation to reduce the variables [14 ]. Recently, a Variation Evolving Metho d (VEM), which is inspired by th e continuous-time dynamics stability t heory, is p roposed for the op timal contro l co mputation [1 5]-[21] . T he VEM also synthesize s the direct and i ndirect method s, but from a new sta ndpoint. T he Evolution Partia l Differential Eq uation ( EPDE), which d escribes the evol ution of variables towards the o ptimal solution, is derived from t he vie wpoint of variation motion, and the op timality conditions will be gradually met u nder this frame. I n Refs . [15] and [16], besides the states and the controls, the costates are also employed in developing the EPDE, and this increases the c omplexity o f the computation. In p articular, the time-opti mal control problems with control constraint are addressed through that thread. However, it is n ot widel y applicable to the general constrained OCPs. In Ref. [17], a compact version of t he VEM that u ses o nly the o riginal variables is proposed . T he costate-free optimality conditions are established and the cor responding EP DE is derived f or th e OCPs w it h free term inal state s. In Refs. [18] an d [19 ], the compact VEM is furthered develop ed to address the O CPs with termi nal Equalit y Co nstraints (ECs) and Inequal ity Constraints (IECs). Normally, u nder the fr ame o f t he compact VEM, the definite co nditions f o r th e EPDE are req uired to be feasible solutions. In Ref s . [2 0] and [21], the Modified EP DE (MEP DE) that is valid even in the infeasible so lution do main is proposed to facilitate the computation o f the OCPs. In this pap er, we further devel op the VEM to solve the general co nstrained OCPs. Throughout the pap er, our work is built upon the assu mption that the solution for the optimization problem exists. We do not describe the e xisting conditions for the p urpose o f brevity. Relevant re searches s uch as the Filipp ov -Cesari t heorem are documented in Ref. [11]. In the following, first the principle of the VEM and the attr ibutes of IEC s in optimization prob lems ar e reviewed. Then th e VEM for t he g eneral constrained OCPs is developed. During this course, the costate-free o ptimality conditions are established , which uncover the anal ytic relatio n of t he co states and multipliers i n the classic treat ment to the state a nd control variables. L ater illustrative e xamples are solved to verify the effectiveness o f the m e thod. B esides, co mparison between the derived EPDE and that in Ref. [15] is pr esented at the end. II. P RELIMINARIES A. Principle of VEM The VEM is a newly d eveloped method for the o ptimal sol utions. It is enlightened b y the states evolution within the stab le continuous-time d ynamic syste m in the control field. Lemma 1 [22 ] (with small adap tation): For a co ntinuous-time autonomous d ynamic system like ()  x f x (1) where n  x is the state, d d t  x x is its time der ivative, and : nn  f is a vector function . Let ˆ x , contained within the domain , be an equilibrium point that s atisfies ˆ ()  f x 0 . If there exists a co ntinuously differentiable function : V  such that i) ˆ () Vc  x and () Vc  x in ˆ / { } x . 3 ii) ( ) 0 V  x in and ( ) 0 V  x in ˆ / { } x . where c is a constant. T hen ˆ  xx is an asymptoticall y stable point in . Lemma 1 aims to the d ynamic system with finite-dimensional states, and it m ay be directly generalized to the infinite-dimensio nal case as Lemma 2 : Fo r an infinite -dimensional d ynamic system desc ribed by () ( , ) x x t    y fy (2) or presented eq uivalently in the P artial Differential Equation (PDE) form as ( , ) ( , ) xt x t    y fy (3) where “  ” denotes the va riation operator and “  ” denotes the partial different ial operator . t  is the tim e. x  is the independent variable, ( ) ( ) n x x  y is the function vector of x , and : ( ) ( ) nn xx  f is a vector function. Let ˆ () x y , contained within a c ertain f unction se t () x , is an equil ibrium function that satisfi es ˆ ( ( ), ) xx  f y 0 . If there exists a continuously different iable function al : ( ) V x  such that i)   ˆ () V x c  y and   () V x c  y in ˆ ( ) / { ( )} x x y . ii)   ( ) 0 V x  y in () x and   ( ) 0 V x  y in ˆ ( ) / { ( )} x x y . where c is a constant. T hen ˆ ( ) ( ) xx  yy is an asymptotically stable so lution in () x . In th e system d ynamics theory, from the stable dynamics, we may construct a monotonously decreasing function (or fu nctional) V , which will achie ve its minimu m when the equilibrium is reached . Inspired by it, no w we consider its inverse proble m, that is, from a perfor mance i ndex funct ion to derive the dynamics that minimize this p erformance index, and opti mization problems are just the right platform for practice. T he optimal solution is analogized to the stable equilibriu m of a d ynamic s ystem a nd is anticipated to be obtained in an asymptotical ly evolvi ng way. A cco rdingly , a virtual dime nsion , the variation ti me  , is introduced to implement the i dea that a variable () t x evolves to th e optimal solution to m i nimize the perfo rmance index within t he dynamics governed by the variation d ynamic evolution equations (in the form o f Eq. (2)) . Fig . 1 illustrates the variation evolution process of the VEM in solving the O CP. T hrough the variation motion, the initial guess of variable s will evolve to the opti mal solution. Fig. 1 The illustratio n of the var iable ev olving along the variation time  in the VEM . 4 The VEM bred under this idea is de monstrated for the u nconstrained calc ulus- of -variations probl em s first [ 15] [17]. The variation dynamic evolutio n equations, derived under the frame of the VEM, ma y be reformulated as the EPDE and the Evolution Differential Equation (EDE ), by r eplacing the v ar iation operation “  ” w ith the partial differential operator “  ” and the differential operato r “ d ” . Under the dynamics gover ned by the EPDE , the variables w il l achieve the opti mality co nditions gradually. For exa mple, consider the c alculus- of -variations pr oblems defined as   0 ( ) , ( ) , d f t t J F t t t t   yy (4) where the elements of the variable vect or ( ) ( ) n t t  y belong to 2 0 [ , ] f C t t . 0 t and f t are the fixed i nitial and ter minal ti me, and the boundary co nditions ar e prescribed as 00 () t  yy and () ff t  yy . T he variation d ynamic evolution eq uation ob tained with the VEM is d () d FF t         yy y K (5) where the column vectors F F    y y and F F    y y are the shorthand notations of partial derivatives, and K is a nn  dimensional positive -definite gain matri x. Co rrespondin gly, the refor mulated EPDE is () FF t          yy y K (6) The equilibriu m solution of th e EPDE (6) will satisfy the optimality co ndition, i.e., the Euler -Lagrange equatio n [23][24] d () d FF t  yy 0 ( 7) Since the right function of the E PDE o nly dep ends on the ti me t , it is suitable to be solved via the well-known semi-disc rete method in the field of P DE numerical ca lculation [25] . W ith the discr etization alon g the nor mal time dimension, the EPDE is transformed to the fi nite-dimensional Initial -value P roblem ( IVP) to be solved, with c ommon Ordinar y Differential Equation (ODE) integration methods. Note that the resulting IVP is defined with respect to the variat ion time  , not the nor mal time t . B. Active IEC and inactive IEC Optimization pr oblems with IECs ar e more intractable. To search the right evol ution equations, Ref. [19] investigated th e attributes of ECs and I ECs in opti mization prob lems and uncover ed the intrin s ic relatio ns to th e multipliers. C onsider the following generalized opti mization problem with the performance inde x as   ( ), gg J J t  yp ( 8) subject to   ( ), , E t t t  g y p 0 (9)   ( ), , I t t t  C y p 0 (10) where ( ) ( ) y n t t  y is the opti mization variable vector and p n  p is the opti mization p arameter vector. Eq. (9) represents the ECs acti ng in the time set E and Eq. ( 10 ) refers to the IECs acting in the time set I . Find the optimal solution   ˆ ˆ ( ), t yp that minimizes g J , i.e.   ˆ ˆ ( ) , ar g m in( ) g tJ  yp ( 11 ) 5 In this ge neral formulation o f optimization prob lems, t he ECs (9) are ca tegorized according to their influence to the optimal performance index, denoted by ˆ g J , and the IECs ( 10 ) are classified according to their acti veness at the o ptimal solution, na mely Definition 1 : Co nsider a specific time po int E E t  and reformulate E q. (9) as   ( ), , EE tt  g y p a (12) where a is a right d imensional vector. For the i th co mponent, if there is 0 ˆ d 0 d i g i a J a   , then E i t g is cat egorized as a Positive-effect EC; if 0 ˆ d 0 d i g i a J a   , then E i t g is categorized as a Negative-effect EC. Here in Definitio n 1, “ 0 i a  ” mean s “evaluated at 0 i a  ” and “ E t ” mean s “evaluated at E t ” . In the follo wing, there are similar implications on “ 0 t ” and “ f t ” . Definition 2 : Co nsider a specific time po int I I t  , for an I EC   ( ), , 0 iI C t t  yp (13) it is said to be an active IEC if   ˆ ˆ ( ), , 0 iI C t t  yp (14) and said to be an inactive IEC if   ˆ ˆ ( ), , 0 iI C t t  yp (15) Note that an inacti ve I EC may be activated for some () t y and p during the optimization process, but we will not call it an active IEC in the pap er. From Definition 2, it is readily to find that strengt hening an IEC ( 13 ) to be an EC as   ( ), , 0 iI C t t  yp (16) t he op timal solution w ill not be changed if t his IEC is a n active IEC. Als o , removing a n inactive IEC from t he optimization pro blem, t he opti mal solutio n will not b e changed either. In ord er to effectively distinguish the IECs within an o ptimization p roblem , we discovered their relatio ns to the strengthened ECs and p rovided a feasible way. See Theorem 1 [1 9]: The IEC ( 13 ) is an ac tive IEC if and only i f the strengthened EC ( 16 ) is a Positive-effect EC, and the IEC ( 13 ) is an inactive IEC if and o nly if the EC ( 16 ) is a Negati ve -e ffect EC. Theorem 2 [19]: Consider a specific ti me point E E t  and use the Lagrange multiplier π to adjo in Eq. (9) with the p erformance index (8). T hen E i t g is a Positive-effect EC if and only if 0 i   . Also, E i t g is a Negative-effect EC i f and o nly if 0 i   . Therefore , we may deter mine the type of an IEC from the multiplier infor mation o f it s s trengthened EC, without t he n eed of substituting optimized solutions into the IEC for verification. In p ractice, w e may first strengthen all I ECs to get the corresponding Lagrange multiplier s, and then use T heorem 2 to determine their types. 6 III. VEM FOR THE G ENERAL CONSTRAINED OCP S A. Problem definition In this pap er , w e consider the general constrai ned OCPs that are defined as Problem 1 ： Co nsider performance index o f Bolza form   0 ( ( ) , ) ( ), ( ), d f t ff t J t t L t t t t    x x u ( 17 ) subject to the d ynamic equation ( , ) t  x f x, u ( 18 ) where t  is the ti me. n  x are the states and each ele ment is piece wise differentiable. m  u are the control inputs and each ele ment is p iecewise di ffe rentiable. T he functio n : nm L    and its first -order par tial derivatives are continuous with respect to x , u and t . The function : n   an d its first-order par tial derivatives are co ntinuous with respect to x and t . T he vector fu nction : n m n    f and its first-order par tial der ivatives are continu ous an d Lipschitz in x , u and t . T he initial time 0 t is fixed a nd the ter minal time f t is free. The initial and ter minal boundary co nditions ar e respective ly prescribed as 00 () t  xx ( 19 )   ( ), ff tt  g x 0 ( 20 ) where : nq  g is a q dimensional vector function with continuous first-order partial derivatives . T he path constraints are describ ed by   ( ), ( ), t t t  C x u 0 (21) where : n m r    C is a r dimensional vector function with continuou s first-order partial derivatives in its augment s. Find the opti mal solution ˆ ˆ ( , ) xu that minimizes J , i.e. ˆ ˆ ( , ) arg min( ) J  xu ( 22 ) The definition of Problem 1 rep resents a large class of OCPs in engineering. Besides t he general B olza form performance index , for the proble ms with no terminal constraints or p ath constr aints, the y are just the degraded cases of P roblem 1 ; for the situatio ns wh ere the terminal time is fixed, they actually bec ome simpler because now t he stud y req uires no deter mination of the ter minal time. Thus, the method develo ped for Pro blem 1 may be widely applied and relevant results are of general mea ning. B. Derivation of varia tion dynamic evolution equation s W e first consider th e problem w it hin the f easible solution do m ain o , in which any solution satisfies Eq s. ( 18 )-( 21 ) . According to the Lyapunov principle, d ifferentiat ing Eq. ( 17 ) with respect to the variation t ime  gives 0 0 T T T T T T T () ( ) ( ) d () ( ) d f ff ff f f f f f t f f f f t tt t t ff t t t t t t t J L L L t tt L L L t                                                   x x u x x x u x xu f x xu f ( 23 ) 7 where : f f tt t   and : f f t   xx . t t      and      x x are the par tial d erivatives , in the form of scalar and (column) vector , respectively. L x and L u are the partial d erivatives. “ T ” is the transpose op erator. For the solutions in o ,   x and   u are related because of Eq. ( 18 ), and they need to sati sfy the followin g variation equatio n as          xu x x u ff ( 24 ) with t he i nitial condit ion 0 t    x 0 . Note that the Jac obi matrixes x f and u f , linearized at the feasible solution () t x and () t u , are time-dependent. Eq . ( 24 ) is a linear time -varying equat ion and has a zero initial valu e. T hus according to the li near s ystem theory [26] , its solution may be explicitly expressed as 0 ( , ) ( ) d t o t t s s s        xu H ( 25 ) where ( , ) o ts H is the nm  dimensional i mpulse response function corr esponding to the specific () t x f and () t u f , name ly ( , ) ( ) ( , ) o o t s s t s ts st       u Φ f H 0 ( 26 ) and ( , ) o ts Φ is the nn  dimensional state t ransition matrix from time point s to time point t , which satisfi es ( , ) ( ) ( , ) oo t s t t s t    x Φ f Φ ( 27 ) In particular, 0 () ( , ) ( ) d f t f of t t t s s s       x u H ( 28 ) Use Eq s . ( 25 ) and ( 28 ), and follo w the similar d erivation as Ref. [1 7]; we may obtain 0 TT ( ) d f ff f t f t t t t J Lt              xu u fp ( 29 ) where   TT ( ) ( , ) ( , ) ( ) d f f t o f o t t L t t t L         u u x x p H H ( 30 ) Now the question of how to find fea sible equations f or   u and f t   arises, which n o t only guara ntee 0 J    but also satisfy the variation equation of the ter minal EC s ( 20 ) as   () f f f f ff t t tt          xx x g g g f g 0 ( 31 ) and the variation motion allowed b y the path constraint s ( 21 ) as TT ( ) ( ) , 1 , 2, ... , 0 p i i i i C C C r t i               xu xu ( 32 ) where : f f tt t  gg and : f f t  xx gg . T he time set p i is defined for the i th path constraint as 0 { | [ ) ] , } ( , 0 , i p i f Ct t t t t    x, u (33) 8 Before answeri ng this questio n, we introduce the Feasibilit y-preserving Evolution Op timization Pr oblem (FPEOP) that is d efined as FPEOP : 3 1 2 11 min 22 .. , 1 , 2, . , 0 .. t t t p i i p J J J st C i r t           g 0 ( 34) where 0 TT 1 ( ) d f ff f t f tt t t t J L t            xu u fp (35) 0 2 2 T1 ( 11 ( ) d 22 ) f f t f t t t t J k t         uu K ( 36 ) with   u being the optimization variable and f t   being the optimization parameter . K is a mm  dimensional positive-definite matrix and f t k is a positive constant. T he time set pp i is a subset of p i defined as 0 ( , ) 0, 0 is a n a cti v e I { | [ , ]} E C, pp i f i i C t t tt C t        x, u (37) Now the question o f determining right   u and f t   will be ans w ered by the following theore m. Theorem 3 : The follo wing variation d ynamic evolution eq uations guarantee tha t the solution stays i n the feasible domain and the change of perfor mance index 0 J    pc    u u Kp ( 38 )   TT () f f f f f f f t t t t t kL          xx f π g f g (39) where T T T T T ( ) ( , ) ( ) ( , ) ( ) ( )d f f t pc o f o t t t t t s t s s s      u u x u x pp Η g π C μ Η C μ ( 40 ) K is the mm  dimensional p ositive-definite gain matrix, f t k is a positive gain co nstant , and u p is defined in Eq. ( 30 ). Assume that t he d ynamic s ystem sati sfies the controllability requirement (See Ref. [27]), then the par ameter vector q  π is calculated by 1   π M r (41) The qq  dimensional matrix M and the q dimension al vector r are   0 T TT ( , ) ( , ) ( )( ) d f f f f f f f ff t o o t t t t ff t t t k t t t      x x x x M g H H g g f g g f g K (42) 9 0 21 ( ) ( )d f f t t t t t    x g r h μ h (43) where   0 T 1 ( , ) d ( ) ( ) f f f f f f f f t o f t t t t t t t t k L        x u x x h g H Kp g f g f (44)   0 T T T 2 ( ) ( , ) ( , ) d ( , ) t o f o o f t t t s t s s t t   xu h H KH C H KC (45) The variable vector   T 12 ( ) ( ) ( ) .. . ( ) () r r t t t t t      μ is determined by 00 TT ( 1 , 2, ..., ) ( ) 0 ( 1 , 2, ..., ) ( ) ( ) ( , ) ( )d ( , ) ( ) d ( , ) ( ) d ( ) 0 ff i t t t W L R A i i i i i t pp i tt pp i when t i r t ti wh r C t t t n d e tt                        u KC μ d μ d μ d μ u ( 46 ) with   0 TT T T 1 T T 1 22 ( , ) ( ) ( , ) ( , ) d ( ) ( ) ( , ) ( ) f f f f t W ii i o o f o f t CC t t t s t s s t t t                       x x x x d H KH g M g h KH g M g h xu (47)   0 TT T T T ( , ) ( ) ( , ) ( , ) d ( ) ( ) ( , ) ( ) L ii i o o o t CC t t t s s s t t                       xu d H K Η C H KC xx ( 48)   0 TT T T T T ( , ) ( ) ( , ) ( , ) d ( ) ( ) ( , ) ( ) t R ii i o o o t CC t t t s s s t t                      xx d H K Η CK Η C xu ( 49)   00 T T T T 1 T T T T 1 11 ( ) ( ) ( , ) d ( ) ( , ) ( , ) d ( ) ( ) ( , ) ff tt A i i i i i o o o f o f tt C C C C d t t s s t s t s s t t               u x u x H Kp H KH g M h Kp KH g M h x x u u (50) Moreover, under the evolution eq uations ( 38 ) and ( 39 ) , 0 J    occurs onl y when () pc t  u p0 (51)   TT ( ) 0 f f f f f tt t L       xx f π g f g ( 52 ) Proof : We will derive Eqs. ( 38 ) and ( 39 ) though the optimization theory. Reformulate Eq. ( 29 ) as a constrained opti mization problem (use 1 t J to denote the per formance i ndex as defined in Eq. ( 35 )) subj ect to constraints ( 31 ) and ( 32 ) . Note that now   u is the optimization variable and f t   is the optimization p arameter . Ho wever, since the mini mum of this opti mization problem may be negative in finity , to penalize too large optimization variab le (par ameter), we i ntroduce anot her perfor mance index 2 t J as defined in Eq. ( 36 ) to formulate a Multi -objective Optimization Pr oblem (MOP) as 12 min( , ) 0 .. , 1 , 2, ..., p i tt i JJ st C i r t          g 0 (53) 10 We use the weig hting method to solve the P areto optimal solution o f this MOP , and the resulti ng performance index is 3 1 2 t t t J aJ b J  ( 54 ) where 0 a  , 0 b  and 1 ab  . W hen 1 , 0 ab  , w e get a solution that minimizes 1 t J . When 0, 1 ab  , w e g et a solutio n that minimizes 2 t J . Otherwise, we get a compromising solution . For this MOP, o bviously in the case of 0, 1 ab  , the Pareto optimal solutio n is t hat 0 f t    and    u 0 , and now the value o f perfor mance indexes are 1 0 t J  and 2 0 t J  . For any other cases, the compromising solution guarantees t hat 1 0 t J  . Set 1 2 a  , 1 2 b  , and because the inactive IECs ma y be removed and the active IECs may be strengthened without changin g the optimal solution, t hen we have the FP EOP defined in Eq. ( 34 ). Introduce the Lagra nge multiplier parameter q  π and the KKT m ultiplier variable r  μ to adjoin the constraints, we m a y get the unconstrained op timization problem from the FPEOP as 0 0 0 0 TT 4 1 2 TT 12 1 1 1 1 d 2 2 2 2 1 1 1 1 ( , ) ( ) d ( ) ( , ) ( ) d d 2 2 2 2 f ff f f f f t t t t t t t t ff t t o f t o t t t t J J J t tt J J t t t t t s s s t                                  x x x u gC π μ u u u π g H g f g μ C H C (55) with ( ) 0 ( 1 , 2, ..., ) p i i p when t t i r     (56) By exchanging the ord er in the d ouble integral and the s ymbols t s , there is   00 0 TT dd ( , ) ( ) d ( , ) ( ) d f ff t tt t t oo tt t t s s s s t t ts t           xx μ CH μ CH uu (57) Thus   00 T T T 4 1 2 1 1 1 1 ( , ) ( ) d ( ) ( , ) d ( ) d 2 2 2 2 f f f f f f f t t t ff t t t o f t o t t t t tt J J J t t t t s t s t t                    x x u x uu π g H g f g μ C μ CH ( 58 ) Use the fist-order optimality conditions, i.e. , 4 t J         u 0 and 4 0 t f J t          , we may get Eq s. ( 38 ) and ( 39 ) . Substitute Eqs . ( 38 ) and ( 39 ) into Eq. ( 31 ); w e have       0 T T T T T TT ( , ) ( , ) ( , ) ( ) ( ) d ( ) d () ff ff f f f f f f f f tt o f o f o tt t t t t t t t t t s t s s s t t kL             x u x x u x x x g H K p H g π HC μ C μ g f g f π g f g 0 (59) Again with the tec hnique of exchanging t he order in the double integral a nd the symbols t s , there is     0 00 T T T T ( ) ( )d ( , ) ( , ) d ( ) ( , ) ( ) d ( , )d ff f tt o o o f o tt f tt tt s s s t t s t s t s s t t tt t     xx H H C μ H H C μ KK ( 60) With further deductio n, we have   0 0 T T T 21 ( , ) ( , ) d ( )( ) ( ) ( )d f f f f f f f f f f f t o f o f t t t t t t t t t t t t k t t t           x x x x x g H KH g g f g g f g π gh μ h ( 61) 11 where 1 h and 2 () t h are defined in Eqs. ( 44 ) and ( 45 ), respectively. Thus, with the definition of M in Eq. ( 42 ) and r in Eq. ( 43 ), Eq. ( 61 ) that determines π is simplified as  M π r (62) Regarding this li near equation, assuming that the control satisfies the controllab ility r equirement [27], then the solution i s guaranteed. T hus the parameter π may be calculated b y Eq. ( 41 ). For ( 1 , 2, ..., ) pp i t i r   , there is 0 i C    , i.e. 0 TT ( , ) ( ) d () 0 () ii t o t CC t s s s          u xu u H ( 63 ) Substituting Eq. ( 38 ) into Eq. ( 63 ) gives     0 T T T T T T T T T T T T ( ) ( , ) ( , ) ( , ) ( ) ( )d ( ) d ( ) ( , ) ( , ) ( ) ( )d ( ) 0 f f f f tt i o o f o ts t i o f o t C t s t s s s s C t t s t s s s t                   u x x u u x x u H K p H g π HC μ C μ x K p H g π HC μ C μ u (64) Substitute Eq. ( 41 ) in and use       0 0 0 0 TT T T T T ( , ) ( , ) ( ) ( )d d ( , ) ( , ) d ( ) ( )d ( , ) ( , ) d ( ) ( )d f f tt oo ts t t t o o o o t t t t t s s s t s s s t s s s                    x xx H K H C μ H KH C μ H KH C μ (65)     0 0 0 0 T T 1 T T 1 22 ( , ) ( , ) ( ) ( )d d ( , ) ( , ) d ( ) ( )d ff f f f f t t t t o o f o o f t t t t t s t s s t s t s s             x x x x H KH g M g h μ H KH g M g h μ (66) Then we have 00 TT ( ) ( ) ( , ) ( )d ( , ) ( ) d ( , ) ( )d ( ) 0 ff t t t W L R A i i i i i t t t C t t t t d t                    u KC μ d μ d μ d μ u (67) where the 1 m  dimensional matri xes ( , ) W i t  d , ( , ) L i t  d , ( , ) R i t  d and the scalar quantit y () A i dt are defined in Eqs. ( 47 )-( 50 ) . Furthermore, Eq . ( 29 ) may be r eformulated as     0 00 0 T TT T T ( ) d ( , ) ( ) d d ( , ) d ( ) f f f f f f f f f f f f t f pc tt t t tt o tt t f o f t t t t J Lt t s s s t t t t t                                   x x u xu xx u f π g f g p uu μ C H C u π g H g f g (68) Since under Eqs. ( 38 ) and ( 39 ), the last two terms in the right par t of Eq. ( 68 ) vanish. Then   0 2 T T T ( ) ( ) d f f f f f f f t pc pc t t t t t J k L t            x x u u f π g f g p Kp (69) T his mea ns 0 J    a nd 0 J    o cc urs on ly whe n Eq s. ( 51 ) a nd ( 52 ) ho ld . ■ Remark 1 : For the opti mal sol ution, t here is ( 1 , 2, ... ) , pp p ii i r   . T he o ptimal value of () t μ (correspo nding to the right pp i ) and the optimal value of π satisfy the gain-i ndependent equatio ns as 12 0 TT ( 1 , 2, ..., ) ( ) 0 ( 1 , 2, ..., ) ˆ ˆ ˆ ˆ ( ) ( ) ( , ) ( ) d ( , ) ( )d ( ) ( ) 0 f i tt L R A i i i i i tt pp i pp i t i r t t i r C t t t d t when whe t n                   π u C μ d μ d μ d π u (70) and 0 1 1 2 2 2 ˆ ( ) ( )d f f s t s t s s t t t                     x M g 0 h μ r r π M (71) where   0 TT T T T ˆ ( , ) ( ) ( , ) ( , ) d ( ) ( ) ( , ) ( ) L ii i o o o t CC t t t s s s t t                       xu dH Η C H C xx ( 72)   0 TT T T T T ˆ ( , ) ( ) ( , ) ( , ) d ( ) ( ) ( , ) ( ) t R ii i o o o t CC t t t s s s t t                      xx dH Η C Η C xu ( 73) 0 TT ˆ ( ) ( ) ( , ) d ( ) t A ii io t CC d t t s s     uu H p p xu (74) 0 T T T T T ˆ ( ) ( ) ( , ) ( , ) d ( ) ( , ) f t ii i o o f o f t CC t t s t s s t t        π x d H H H g xu (75) The qq  dimensional matrixes 1 s M , 2 s M and the q dimension al vectors 1 s r , 2 s r are   0 TT 1 ( , ) ( , ) d f ff t s o f o f t t t t t t   xx M g H H g ( 76 ) T 2 ( )( ) f f f f f s t t t    xx M g f g g f g ( 77 ) 0 1 ( , ) d f f t s o f t t t t   xu r g H p (78) T 2 ( ) ( ) f f f f f s t t t L      xx r g f g f (79)   0 T T T 2 ˆ ( ) ( , ) ( , ) d ( , ) t o f o o f t t t s t s s t t   xu h H H C H C (80) Proof : Regarding the ar gument that ( 1 , 2, ... ) , pp p ii i r   for the op timal solutio n o f Pr oblem 1 , thi s i s bec ause a ny ti me point t in p i also belongs to pp i ultimately, or this activated inequalit y path constraint at time t will become inactive. For the op timal values o f () t μ and π , since K may b e arb itrary right -dimensional positive -definite matrix, we set K = 1 in Eq. ( 64 ) , wh er e 1 is the mm  dimensional identit y matrix. Then we get E q. ( 70 ) after appropriate simplification. From Eq . ( 61 ) , again bec ause K may be ar bitrary ri ght-dimensio nal positive-definite matri x and f t k may b e arbitrary positive co nstant, w e consider three case s, t hat are, i) K = 1 , 1 f t k  , ii) K = 2 1 , 1 f t k  , and iii) K = 1 , 2 f t k  . B y comparing the three cases o f substituting t he specific values into Eq. ( 61 ), we may obtain Eq. ( 71 ), which is irrelevant to the specific v alue of K and f t k . ■ 13 C. Equivalen ce to the classic optimality con ditions Actually, Eqs. ( 51 ) and ( 52 ) are the f irst-order optimality conditions f or Problem 1 without the employment of co states. We will show that they are eq uivalent to the classic o nes with costa tes [28] . By the d irect adjo ining method [11] , we may co nstruct the augmented functional as     0 T T T ( ( ), ) ( ), ( ) d f t f f f t J t t t t L t         x π gx λ fx μ C (81) and ( ) 0 ( 1 , 2, ..., ) ( ) 0 ( 1 , 2, .. ., ) p i i i i p i t t i r t C wh e tr n wh n i e         (82) where n  λ is the co state variable , q  π is Lagrange multiplier parameter, and () r t  μ is the KKT multiplier variable . Then the first-ord er variation may be derived as       0 T TT T T T T ( ) ( ) ( ) ( ) d f f f f f f t t f f f t t t J H t t t H H H t                       xx λ xu π g λ g π x x λ λ x u C μ (83) where TT HL    λ f μ C is the augmented Hamiltonian. T hrough 0 J   , we have T T HL       x x x x λ λ f λ C μ 0 ( 84 ) T T HL     u u u u f λ C μ 0 ( 85 ) and the transversalit y conditions TT ( ) 0 ff f tt t L      λ f π g ( 86 ) T () ff f t     xx λ g π 0 (87) and the KKT conditio n ( 82 ). Theorem 4 : For Pr oblem 1, the optim alit y conditions given by Eq s. ( 51 ) and ( 52 ) are equivalent to the optimality co nditions given by Eqs . ( 82 ) , ( 84 ) -( 86 ). Proof : Define a variable () t γ as T T T T T T ( ) ( , ) ( , ) ( , ) ( ) d ( , ) ( ) ( )d ff ff tt o f o f o o tt t t t t t t L t              x x x x γ Φ Φ g π Φ Φ C μ ( 88 ) Then Eq. ( 51 ) is simplified as TT u L    uu f γ C μ 0 ( 89 ) Obviously, when f tt  , there is T () ff f t   xx γ g π (90) Under Eqs. ( 38 ) and ( 39 ) , Eq. ( 68 ) may be simplified as     0 T TT ( ) d f f f f f f t f pc tt t t t J Lt              x x u u f π g f g p (91) which holds in the feasible solution domain o . Further combined with Eq. ( 90 ) and ignore  , we have 14     0 T TT ( ) d f ff f t pc t t f f t t J L t t t           u π g γ f p u (92) Eq. ( 92 ) obviously hold s at th e optim al solution. C o mpare Eq . ( 92 ) w it h Eq. ( 83 ) and consider the v ar iation of th e term inal tim e. T o achieve the extre mal condition, Eqs. ( 52 ) and ( 86 ) should be same, i.e.     T T T T () f f f f ff f t t t t tt L t L         γ f π g λ f π g (93) Since Eq. ( 93 ) generally hold for arbitrary g and f , we can conclude that  π π (94 ) ( ) ( ) ff tt  γ λ ( 95 ) Therefor E q. ( 90 ) is same to Eq. ( 87 ) . Consider the variation on the control w i thin the f easible solution do main in Eqs . ( 92 ) an d ( 83 ). There should be the conclusion that Eq s . ( 89 ) and ( 85 ) are identical, namely T T T T u L L      u u u u u f f γ λ C μ C μ (96) Because ( 1 , 2, ... , ) pp p ii i r   for the optimal soluti on, Eq. ( 96 ) implies ( ) ( ) tt  μ μ (97) ( ) ( ) tt  λ γ (98) Differentiate () t γ , as defined in Eq. ( 88 ) , with respect to t . In the proce ss, we will use the Leibniz rule [29]       ( ) ( ) ( ) ( ) 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         ( 99 ) and the property of ( , ) o t  Φ [26] ( , ) ( , ) ( ) o o t tt t      x Φ Φ f (100) ( , ) o tt  Φ 1 ( 101 ) where 1 is the nn  dimensional identit y matrix. Then we have   T T T T T T T T T T T T T T T T T T T d ( ) ( , ) ( , ) ( , ) ( ) d ( ) ( ) ( , ) ( ) ( )d d ( ) ( ) ( , ) ( , ) ( , ) ( ) d ( , ) ( ) ( )d ff ff ff ff tt o f o f o o tt tt o f o f o o tt t t t t t L t L t t L t t t t t L t                                      x x x x x x x x x x x x x x x x x γ f Φ f Φ g π f Φ C μ f Φ C μ C μ f Φ Φ g π Φ Φ C μ TT ( ) ( ) ( ) Lt   x x x f γ C μ (102) With Eq. ( 97 ) , t his means () t γ confor ms to the same dynamics as the costates () t λ , and Eqs. ( 102 ) and ( 84 ) are exactly the same . T hus the t heo re m i s pr ove d. ■ From T heorems 3 a nd 4 , w e have got the explicit analytic relations o f the costates λ , the L a grange multiplier para meter s π and the KKT multiplier variables μ for the classic treatment i n Eq. ( 81 ) to the original (state a nd control) variables , which formerly can only be obtained numerica lly by solving the BVP . To present the se results more clearly, they are summarized in T able 1. 15 Table 1 The class ic optimality conditions and the costate-fre e optimality conditions for Pr oblem 1 Feasibility condit i ons i) ( , ) t  x f x, u ; ii) 00 () t  xx ; iii)   ( ), ff tt  g x 0 ; i v)   ( ), ( ), t t t  C x u 0 . Classic optimal ity conditions i) TT L     x x x λ f λ C μ 0 ; ii) T T L   u uu λ C μ f0 ; iii) T T () 0 f f f t t t L      π g λ f ; iv) T () f f f t     xx λ π g0 ; v) ( ) 0 ( 1 , 2, ..., ) ( ) 0 ( 1 , 2, ..., ) p i i i i p i t t i r t C whe tr n wh n i e         . Costate-free optimality conditions i) pc  u p0 in which π is calculated by Eq . ( 41 ) and () t μ is determined by Eq. ( 46 ) ; ii) TT ( ) 0 ff ff f tt t L       xx f π g f g where π is calculated by Eq . ( 41 ). Analytic rel ations i)     T T T T ( ) ( , ) ( , ) ( ) ( ) ( ) d f f f t o f o t t t t t L            x x x x λ Φ g π Φ C μ ; ii) 1   π M r where M is given in Eq. ( 42 ) and r is given in Eq. ( 43 ) ; iii) 00 TT ( ) ( ) ( , ) ( ) d ( , ) ( ) d ( , ) ( )d ( ) 0 ( ) 0 ff t t t W L R i i i i i t t t A i C t t t t t for t d                      u KC μ d μ d μ d μ u wit h ( , ) W i t  d , ( , ) L i t  d , ( , ) R i t  d and () A i dt defined in Eqs. ( 47 )-( 50 ) . After the pr oof of Theorem 4, n ow the v ar iables evolving direction using t he VE M is easy to determ ine and t he optimal solutio n of Problem 1 will be sought with theor etical guarantee. Theorem 5 : Solving the IVP with respect to  , defined b y th e variatio n d ynamic e volution equations ( 25 ) , ( 38 ) and ( 39 ) fro m a feasible initial solutio n, when     , ( , ) xu will satisfy the opti mality conditions of Prob lem 1. Proof : By Lemma 2 and with Eq. ( 17 ) as the Lyapunov functional, we may claim that the minimum solution of Prob lem 1 is an asymptotically stable solution within the feasibility do main o for the i nfinite-dimensional d ynamics governed by Eqs. ( 25 ) , ( 38 ) and ( 39 ). Fro m a feasible initial solution, any evolution und er these d ynamics maintains the feasibility o f the variables, and they also guarantee 0 J    . The functional J will decrease u ntil 0 J    , which occurs w hen     due to the asymptotical approach. When J   =0, this deter mines the op timal co nditions, na mely, Eqs. ( 51 ) and ( 52 ). ■ D. Formulation o f EPDE Use t he partial di fferential o perator “  ” and the di fferential operator “ d ” to refor mulate t he variation d ynamic e volution equations, we ma y get the EPDE and the EDE as 0 ( , ) ( , ) d ( , ) ( , ) t o t pc s t s s t t                      u u H x u Kp ( 103)   TT d () d f f f f f f f t t t t t kL         xx f π g f g ( 104 ) Put into this perspective, t he d efinite conditions are the initial guess of f t , i.e., 0 f f tt    , and 16 0 ( , ) ( ) ( , ) ( ) t t t t                  xx uu (105) where () t x and () t u are the feasible initial solutions. Eqs. ( 103 ) and ( 104 ) realize the a nticipated variable evolving along the variation time  as d epicted in Fi g. 1. T he initial conditions of ( , ) t  x and ( , ) t  u at 0   belong to the f easible sol ution domain and their values at    represent the optimal solution of t he OCP. T he right par t of the EP DE (103) is also only a v ec tor function of time t . Thus w e m a y appl y the semi-discrete method to discretize it alo ng the nor mal time dimension an d further use ODE i ntegration methods to get the n umerical solutio n. Meanwhile, t he Lagrange multiplier parameters π and the KKT multiplier variables () t μ need to be solved during the evol ution process. To determine th e right functions of Eqs. (103) and ( 104 ) at any variat ion time  , gen er ally w e compute () t μ through Eq. ( 46 ) first, and then u se Eq . ( 41 ) to get π . For the integral equation in Eq. ( 46 ), it may be solved numerically at t he d iscretization time points th at b elong to ( 1 , 2, ..., ) pp i i r  . Note that whether a specific time point b elongs to ( 1 , 2, ..., ) pp i i r  may b e determined in light of T heorems 1 and 2 . IV. D ISCUSS ION A. Various path constraints The path constraints ( 21 ) defined in P roblem 1 take the mixed state-control for m. Act ually they can also include other par allel forms. For example, i f the path constraint is in the form of pure -state type , i.e.   ( ), tt  C x 0 (106) t hen  u C0 and all ter m s in the evolution equations releva nt to u C vanish . When the path co nstraint takes the pure -control form as   ( ), tt  C u 0 (107) then  x C0 and the ter ms relevant to x C disappear in the evolution equations. In par ticular, the integral equation in Eq. ( 46 ) that determines () t μ is now simplified as 0 TT ( ) ( ) ( , ) ( )d ( ) 0 f t WA i ii t C t t d t          u KC μ d μ u (108) which is a typical Fred holm integral equatio n. Besides the inequality ty p e path constraint, we may also investigate the equality type. Consider an extreme ca se that all th e path constraints in ( 21 ) take the equality form, suc h as   ( ), ( ), t t t  C x u 0 (109) For such case, the treatme nt in deriving the e volution eq uations i s still same, w hile now the determination of the mu l tiplier v ariable μ becomes simpler, j ust through solvin g the following inte gral equation for the whole time horizon 0 [ , ] f tt as 00 T (t) ( , ) ( )d ( , ) ( ) d ( , ) ( )d ff t t t W L R A t t t t t t                  uu CC μ D μ K μ D μ D d 0 (1 10) where   0 T T 1 T T 1 22 ( , ) ( ) ( , ) ( , ) d ( ) ( ) ( , ) ( ) f f f f t W o o f o f t t t t s t s s t t t         x x x u x x D C H KH g M g h C KH g M g h (1 11) 17   0 T T T ( , ) ( ) ( , ) ( , ) d ( ) ( ) ( , ) ( ) L o o o t t t t s s s t t         x x x u D C H K Η C C H KC (1 12)   0 T T T T ( , ) ( ) ( , ) ( , ) d ( ) ( ) ( , ) ( ) t R o o o t t t t s s s t t        x x u x D C H K Η C C K Η C (113)   00 T T 1 T T 1 11 ( ) ( , ) d ( , ) ( , ) d ( , ) ff tt A o o o f o f tt t t s s t s t s s t t       x u x x u u u x d C H Kp C H KH g M h C Kp C KH g M h (1 14) When t here are both inequality and equality path constraints in th e O CP, th en Eqs. ( 46 ) and ( 110 ) need to be s ynthesized in seeking the multiplier variables. B. Simpler OCP formulation s Since being derived for the general constrai ned OCPs, the equations obtained in t his pap er are generally applicable . Now we consider simpler OCP formulations . If the terminal time f t is fixed in the OCP , then Eqs. ( 103 ) and ( 104 ) may be directly applied b y setting 0 f t k  . If there is no path co nstraint ( 21 ) in the OCP, then we do not n ee d to solve () t μ and it m a y b e set as () t  μ 0 in Eqs. (103) and ( 104 ) to get the optimal solutio n. When there is no ter minal con straint ( 20 ) , then we may let  π 0 in Eqs. (103) and ( 104 ), and the KKT multiplier variables () t μ is now deter mined by 0 TT ( 1 , 2, ..., ) ( ) 0 ( 1 , 2, ..., ) ( ) ( ) ( , ) ( ) d ( , ) ( )d ( ) 0 f i tt L R A i p i p i i t pp i t i t i r t t i r C t t t d t when when                  u KC μ d μ d μ u ( 115) w it h ( , ) L i t  d defined in Eq. ( 48 ) , ( , ) R i t  d defined in Eq. ( 49 ) , and () A i dt being 0 TT ( ) ( ) ( , ) d ( ) t A ii io t CC d t t s s     uu H Kp Kp xu (1 16) C. Generation o f a feasible initial solution The d erivation in Sec. III starts from the p remise that the i nitial so lutions () t x and () t u are f ea s ib le, na mely, satisfying Eqs. ( 18 ) -( 21 ). For the general constrained OC Ps, usually it is n ot an ea sy ta sk to find a f easible solution. We may eit her fu rther develop the VEM to be valid in the infeasible solution d omain as Refs. [20] and [21] did, or w e find an a lternative to determine a f ea sible in itial solution. Here we choose the second way. Upon the work in Ref. [21], which solves the OCPs with arbitrary initializatio n, we can generate a feasible i nitial solution b y solving the following Feasible Solutio n Searching Opti mization Problem (FSSOP ) as FS S OP :     0 1 00 ( ), ( ), d ( , ) () () .. , min f r t fs i i t i ff J w C t t t t t t t st t           xu x f x , u xx g x 0 ( 117) 18 where i w is the w eig ht co efficient for th e i th path constraint i C in ( 21 ). T he dynamics constraint a nd the bo undary conditions ar e same to those in Pr oblem 1, while the terminal time f t may be different. Through solving thi s FSSOP with t he method in Ref. [2 1], we may obtai n a solution that is feasible for P roblem 1 . D. Numerical soft barrie r Theoretically, the e volution equations will precisely see k th e o ptimal solution. During the variable evolution process, once the inequality path constraint at a specific ti me point is activated, the corresponding variation constraint (in E q. ( 32 ) ) w ill be triggered immediately to maintain the feasibility of sol utions. However, since we reso rt to the numerical method for the s olution, concretely by u sing the ODE integrati on methods to solve the tra nsformed finite-dime nsional IVPs, the numerical erro r is unavoidable, and this may lead to the v io lation of the path constraints. Refer to the strategy to eliminate the violatio ns on th e ter minal IECs [19], the numerical soft barrier technique is again emplo yed to remove the possible numerical error on the path constraints, by ada pting the FPEOP ( 34 ) as Adapted FPEOP : 3 1 2 11 min 22 .. , 1 , 2, ..., 0 t t t i Ci pp i J J J st C k C i r t             g 0 (1 18) where C k is a positive co nstant and now the ti me set pp i is defined as 0 ( , ) 0 { | [ , , 0 is a n a cti v ]} e I EC , i i C i pp i f C C t k t t t t C         x, u (11 9) Through solving t he adapted FPE OP, the evolution equatio ns derived are still similar except () A i dt in Eq. ( 46 ) is modified as     00 T T T T 1 T T T T 1 11 ( ) ( ) ( , ) d ( ) ( , ) ( , ) d ( ) ( ) ( , ) ff tt A i i i i i o o o f o f C i tt C C C C d t t s s t s t s s t t k C                u x u x H Kp H KH g M h K p KH g M h x x u u (120) In this way, the possible violations on the path constraints du e to the numerical erro r will be eli minated gradually. V. I LLUSTRATIV E EXAMPL ES First a linear exa mple taken from Z h a n g [30 ] is solved. Example 1 : Consider t he following dynamic system u  x A x b where 1 2 x x     x , 01 00     A , and 0 1     b . Find the solution that mini mizes the perfor mance index f Jt  wi th t he control constraint 11 u    and the boundary co nditions 0 1 () 1 t     x , 0 () 0 f t     x 19 where the initial ti me 0 0 t  is fixed. In solving this exa mple using the VEM , the pat h constrai nt is re formulated as 2 10 u  T he n the EPDE derived is   0 () T () T ( ) d 2 ( ) f t ts t tt e s s u K e u t                         A A u b x b π where π is solved by Eq. ( 41 ) and the scalar KKT multiplier variable () t  is deter mined by Eq . ( 46 ) . The o ne-dimensional gain matrix K wa s 0. 2 K  and the scalar f t k was 0.1 f t k  . T he barrier parameter C k in Eq . (120 ) was set to be 0.1. The definite conditions of the EP DE, i.e., the feasible initial guess of the states () t x and the control () u t , were obtained by solvin g the following FS S OP as 0 2 0 0 1 d 2 10 ( ) , ( ) min .. 10 0, 8 f t fs t f f J u t u tt st tt                  x A x b xx Note that f t =8s is also the initial g uess of the terminal time f t for this example. Using the se mi-discrete method, the time horizon 0 [ , ] f tt was discretized unifor mly with 41 points. T hus, a dynamic s ys tem with 1 24 states (including the ter minal time) wa s obtained and the OCP wa s transformed to a finite-di mensional I VP. The ODE integrator “ode45” in Matlab , with default relative error toler ance 1 × 10 - 3 and default absolute erro r toler ance 1× 10 -6 , wa s e m ployed to solve the IVP . For comparison, t he analytic solution by solving the B VP is also presented. 2 2 1 1 22 11 22 [ 1 ( 6 / 2) , 1 6 ˆ 0.5 ( ) 3.5 ˆ 0.5 1 ˆ ˆ 11 ˆ ˆ 6 / 3 6 / 3 ˆ ˆ ( 6 / 3) ( 6 / 3) 1 ( 6 / 3) ( 6 ] [0, 1 ( 6 / 2 ) / 3 ) 1 ) 1 6 6 ˆ ˆ 1 6 1 x t t x t t x t t t u t t x t u                                                 Fig s. 2 and 3 present the evolution process o f 1 () xt and () ut towards the analytic solutio ns, sh owing the asymptotically approach of the numerical res ults to the opti mal. A t  = 300s, they are ver y close to the analytic sol utions, and this de monstrates the effectivene ss o f the VE M . For the con trol results p lotted in Fi g. 3, it is shown that the control s witch is acc urately capt ured from the close- up . In Fig. 4 , the stat es results are again compared with t he analytic sol ution i n the state plane, illustrating t he evolution process of the states from a different an gle . T he pr ofile for the terminal t ime is give n in Fig 5 . It m onoto nously decreases from f t = 8s and is almost unc hanged after  = 100s. At  = 300 s, we co mpute tha t f t = 3.44 s, very clo se to t he analytic result. Regarding t he 20 Lagrange multipliers t hat adj oin the terminal constraint, we computed that 1.00 0.83        π at  = 300s . From t he a nalytic relation to the costates in T able 1, we have   T () T 0.83 0.83 0.83 1. 10 ( ) ( 1.00 8 , 1 6 ) f tt of f t t t t t e t                  A λ Φ π π This is close to the analytic s olution of 1 2 ˆ ˆ       . In F ig . 6 , the numerical solution of KKT multiplier variable () t  at  = 3 00 s is presented, and the s harp angle of the curve at the control s witch time point is clearl y shown. 0 1 2 3 4 5 6 7 8 0 0.5 1 1.5 2 t x 1 ( t ) The a nal yt i c s o l ut i o n N um eri c al s o l uti o ns wi th VE M  = 30 0s  = 50 .1s  = 31 .1s  = 22.9s  = 9.3s  = 0s Fig. 2 The evolution of numer i cal sol u tions of 1 x to the analy t ic solution. 0 1 2 3 4 5 6 7 8 -1 . 5 -1 -0 . 5 0 0.5 1 1.5 t u ( t ) The a nal yti c s o l uti o n N um eri c al s ol uti o ns w i th V EM 2 2.5 -1 0 1  = 9.3s  = 22.9s  = 31.1s  = 50.1s  = 0s  = 300 s  = 30 0s Fig. 3 The ev olution of numerical solutions of u to the analy ti c solution. 21 -0 . 5 0 0.5 1 1.5 2 -1 . 5 -1 -0 . 5 0 0.5 1 x 1 ( t ) x 2 ( t ) The an al yti c s o l uti o n N um eri c a l s ol uti o ns w i th V EM  = 0s  = 9.3s  = 50.1s  = 300 s  = 22 .9s  = 31 .1s Fig. 4 The ev oluti on of numeric al solutions in 12 xx state plane to the analy t ic solution. 0 50 100 150 200 250 300 3 4 5 6 7 8  (s) t f (s ) M i ni mu m co nt rol ti me E vo l vi ng pro fi l e o f t f Fig. 5 The ev olution profile of f t to the analy t ic result. 0 0 .5 1 1.5 2 2.5 3 3.5 0 0.2 0.4 0.6 0.8 1 t  ( t ) K KT mu l ti p l i er v ari abl e at  = 30 0s Fig. 6 The KKT mu ltiplier v a riable profil e for th e optimal so lution with VE M . 22 Now we consider a nonlinear example with pure-state path c onstraint , the constrained Brachistochrone pr oblem [28], which describes the motion curve of the fastes t descending with position constraint. Example 2 : Consid er the following d ynamic system ( , ) u  x f x where x y V       x , sin( ) cos( ) cos( ) Vu Vu gu       f , and 10 g  is the gravity constant. Find the solution that minimizes the per formance index f Jt  with the boundary co nditions 0 0 0 0 0 t x y V             , 2 f t x  In addition, during the descending, t he position states are constrained b y a slope as   ( ), ( ) 0.5 ( ) ( ) 0.35 0 C x t y t x t y t      In the sp ecific form of t he EP DE (103) and the EDE (104), the gain p arameters K and f t k were set to b e 0.1 and 0.05 , respectively . T he barrier parameter C k in Eq. (120) was set to be 0 .2. The d efinite conditions, i.e ., 0 ( , ) ( , ) () f t ut t           x , were ob tained from a physical motio n along a straight line that co nnects the initial positio n to the terminal positio n o f 2 1     , i.e. 22 a rc t a n ( 2) 2 2 1 5 f u x t t t yV t       We also d iscretized the time horizon 0 [ , ] f tt uniformly, with 101 points . T hus, a large IVP with 405 states (including the terminal time) wa s ob tained. W e still e mploy ed “ ode 45 ” in Matlab fo r the nu merical inte gration. In the i ntegrator setting, the default r elative error tolerance and the absolute error tolerance were 1 × 10 -3 and 1 × 10 -6 , respectively. For comparison, we computed the optimal solution with GP OPS- II [31], a Radau PS method based OCP solver . Fig . 7 gives t he state s cur ve in the xy coordinate plane, s howing that the n umerical resul ts starting fro m the straight line approach the optimal solution o ver time, and the op timal d escending curve is constrained by t he slope. T he control solutions are plotted in Fig. 8 . T he asymptotical appr oach of the numerical results is demonstrated, and the restriction effect from t he slope on the control is clearl y s hown . I n Fig. 9, the ter minal time profile aga inst the variatio n time  is plotted . T he result of f t declines rapidly at fir st and then graduall y ap proaches the minimum decline ti me, and it only changes slightly after  = 50s. At  = 300s, we compute th at f t = 0.80 01s from the VEM , very close to the r esult of 0 .7999s f r om GPOPS- II . Fig. 10 p resents the pro files of the path constrai nt   ( ), ( ) C x t y t and the KKT multiplier variable () t  at  = 300s, with respect to the x position co ordinate . T he path constraint is acti ve within the co ordinate interval 0 .5 6 1 .0 6? x  , and the correspo nding positive KKT multiplier variable is obviously shown . Different from t he results in Fig. 6 , t he ir values for the active path constraint oscillate due to th e nu merical error arising from discretizatio n. 23 0 0 . 5 1 1.5 2 2.5 -1 . 5 -1 -0 . 5 0 x ( t ) y ( t ) The op ti m al s o l uti o n N um eri c al s o l uti o ns wi th VE M  = 2.5s  = 12.2s  = 6.8s  = 4.0s  = 0s  = 300s The slope Fig. 7 The ev oluti on of numeric al solutions in t he xy coordinate plane to the optimal solution. 0 0 .2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 t u ( t ) The op ti m al s o l uti o n N um eri c al s o l uti o ns wi th V EM 0.4 0.5 0 .6 1.0 5 1.1  = 0s  = 2.5s  = 4.0s  = 6.8s  = 12.2s  = 300s  = 300s Fig. 8. The evo lut ion of numeri cal solutions of u to the optimal so lution. 0 50 100 150 200 250 300 0.7 5 0.8 0.8 5 0.9 0.9 5 1  (s) t f (s ) M i ni mu m decl i ne ti me E vo l vi ng p rofi l e o f t f Fig. 9 The ev olution profile of f t to the minim um decline time. 24 0 0.5 1 1.5 2 -0 .4 -0 .2 0 0.2 x C Pro fi l e of the p ath c o nstra i nt at  = 30 0s 0 0.5 1 1.5 2 0 0.5 1 x  K KT mu l ti p l i er v ari ab l e at  = 30 0s 0.5 1 -4 -2 0 x 1 0 -3 x C Fig. 10 The p rofile s of the path constr a int and t he KKT multiplier with VEM. VI. F URTHER C OMMEN TS The EPDE derived under the compact VEM is first presented in Ref. [17 ], and we gave an immature discussion then between it and the a ugmented EPDE deriv ed in R ef. [15], which wa s cal led the first evolution eq uation or iginally. Since w e have dee pened the study along the thread of the c ompact VEM and achieved s ystema tic results, it makes sense to give a revie w . To disting uish, t he EPDE (103) derived here is n amed the second evolution equation. For convenience, the fir st evolution equation is again presented. ( , ) H t t t H H t t t H                                                          x yy x u λ x f y K x λ f 0 (121) where       x y λ u , T HL  λ f is the Ha miltonian, and λ is the costate variable s. K is a (2 ) ( 2 ) n m n m    dimensional positive-definite matrix. Both the first a nd the second evolution equations or iginate from the cont inuous-ti me dynamics stabilit y theor y, and their solutions are guarantee d to ultimately meet t he o ptimality co nditions. The right parts of bot h equations ar e on ly v ecto r f unctions of time t . T his makes t hem suitab le to be solved with t he semi-d iscrete method in t he field of PDE nu merical calc ulation. T hen the numerical solution may be obtained with the commo n ODE integration methods. However, there exist obvious di fferences. The f irst evol ution equa tion is derived from a co nstructed un co nstrained f u nctional v i a em p loying the classic op timality conditions with costates. It may handle typical OCPs w ith terminal co nstraint [15]. Ho wever, ev en if it has solved the time-op tima l control problem with control constrain t [16], it is not applicable to the g eneral constrained OC Ps (at least for now ). The introduction o f the costates als o co mplicates the formula and intensifies the co mputation burden. In particular, its solution may halt at a saddle po int since it cannot d ifferentiate the minimu m from the saddle upon t heir first- order optimality conditions. The second evolution equatio n searches the mi nimum so lution from t he pri mary proble m, and the equivale nt costa te-free optimality co nditions are established meanwhile, which unco ver the analytic r elations bet w ee n the original variables and the augmented quantities, including the costates a nd the multiplier s . It has been shown that the second evolution eq ua tion may solve 25 general constrained O CPs and various t ypical OCPs, and it ma y be modified to be valid in the infeasible solutio n domain as Refs. [20] and [21] show. In principle, the second evolution equatio n requires th e integration, and the differentiation, as displayed in the first evolution equatio n, may be avoided . T his is advantageous to reduce the nu merical error in seeking opti mal solutions. VII. C ONCLUSION The Var iation Evolving Metho d ( VEM) is d eveloped to solve the general state- and \or control-constrained Optimal C ontrol Problems (OCPs) . In der iving t he evolution equations, the costate-free optimality conditions are established, and the anal ytic relations bet ween the original variab les and the c ostates, the KKT multiplier variables , the Lagrange m ultiplier p arameters , introduced in the classic treatment, are uncovered. T he se results are authenticated between the VEM and th e adjoining method, an d are helpfu l to deepen the understanding on t he o ptimal control theor y . In our w ork, the studies of th e VEM are carried out upon the assumption that the solu tion o f the OCP exists. Actually this may ofte n be ascertained through the p hysical analysis. Once the existence of t he solution is se cured, the VEM th eoretically guara ntees the co nvergence to the optimal solution . For the user, this method allows automatically generated initial guess, thus it may be an initial-guess free method for th e users. Moreover , the VEM mainly r equires co mmon Ord inary Dif ferential Eq uation ( ODE) nu merical integratio n to get the solution. Althou gh we did not highlight the small ti me consu mption, the solutions ar e usuall y o btained fairly fast. In its application, since co mplex numerical computations are avoided, the inte gration may be i mplemented with the simple analo g cir cuit . As an outlook, t hese merits might make the achieveme nt of more reliable and practical on -lin e optimal control in engineering p ossible. R EFERENCES [ 1 ] H . J . P e s c h a n d M. P l a i l , “ 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 , ” C o n t r o l & C y b e r n e t i c s , v o l. 38 , n o . 4 , p p . 9 7 3 - 995 , 2 0 0 9 . [ 2 ] J . T . B e t t s , “ S u r v e y o f n u m e r i c a l 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 , ” J . G u i d . C o n t r o l D y n a m . , v o l . 2 1 , n o . 2 , p p . 1 9 3 - 2 0 6 , 1 9 9 8 . [ 3 ] Q . L i n , R . Lo x t o n , a n d K. L. T e o , " 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 o u r n a l o f In d u s t r i a l a n d M a n a g e m e n t O p t i m i z a t i o n , v o l . 1 0 , n o . 1 , p p . 2 7 5 - 3 0 9 , 2 0 1 4 [ 4 ] C . H a r g r a v e s a n d W . P a r i s , “ D i r e c 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 e a 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 a m . , vo l . 1 0 , n o . 4 , p p . 3 3 8 - 3 4 2 , 1 9 8 7 . [ 5 ] O . V . S t r yk a n d R . B u l i r s c h , “ 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 . , v o l . 3 7 , n o . 1 , p p . 3 5 7 - 3 7 3 , 1 9 9 2 . [ 6 ] H . J . P e n g , Q . Ga o , Z . G. W u , a n d W . X . Zh o n g , “ S ym p l e c t i c a p p r o a c h es 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 y n a m . , v o l . 3 5 , n o . 2 , p p . 6 5 3 - 6 5 8 , 2 0 1 2 . [ 7 ] A. V. Rao , “ A s u r v e y of n u m e r i c a l m et h o d s f o r o p t i m a l c o n t r o l , ” in Pr o c . A A S / A I A A A s t r o d y nam . Sp e c i a l i s t C o n f . , Pi t t s b u r g h , PA , 2 0 0 9 , AA S P a p e r 0 9 - 334. [ 8 ] D . Ga r g , M . A. P a t t e r s o n , W . W. H a g e r , A. V . R a o , e t a l , A U n i f i e d f r a m e w o r k f o r t h e n u m e r i c a l s o l u t i o n o f 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 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 , v o l . 4 6 , n o . 1 1 , p p . 1 8 4 3 - 1 8 5 1 , 2 0 1 0 . [ 9 ] I. M. Ro s s a n d F. F a h r o o , “ A p e r s p e c t i v e on m e t h o d s f o r tr a j e c t o r y o p t i m i z a t i o n , ” in P r o c . AI A A / A A S A s t r o d y n a m . C o n f . , 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 - 4 7 2 7 . [ 1 0 ] H . J . P es c h , “ R e a l - t i m e c o m p u t a t i o n o f f e ed b a c k c o n t r o l s f o r c o n s t r a i n e d o p t i m a l c o n t r o l p r o b l e m s , p a r t 1 : n ei g h b o r i n g e x t r e m a l s , ” O p t i m a l C o n t r o l A p p l i c a t i o n s & M e t h o d s , v o l . 10 , n o . 2 , p p . 129 - 1 4 5 , 1 9 8 9 . [ 1 1 ] R . F . H a r t l , S . P . S e t h i , a n d R . G. V i c k s o n , “ 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 I A M R e v . , v o l . 3 7 , n o . 2 , p p . 1 8 1 - 2 1 8 , 1 9 9 5 . 26 [ 1 2 ] Q. G o n g , W. K a n g , a n d I. M. R o s s , “ A p s e u d o s p e c t r a l m e t h o d f o r t h e o p t i m a l c o n t r o l o f c o n s t r a i n ed f e e d b a c k l i n e a r i z a b l e 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 , v o l. 51 , n o . 7, pp. 1 1 1 5 - 1129 , 2 0 0 6 . [ 1 3 ] Q . G o n g , I. M . R o s s , W . K a n g , a n d F . Fa h r o o , “ C o n n e c t i o n s b e t w e en t h e c o v e c t o r m a p p i n g t h e o r e m a n d c o 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 p p l . , v o l . 41 , p p . 307 - 335 , 2 0 0 8 . [ 1 4 ] I. M . R o s s a n d F . F a h r o o, “ 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 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 , v o l . 49 , n o . 8 , p p . 1410 - 1 4 1 3 , 2 0 0 4 . [ 1 5 ] S . Z h a n g , E . M . Y o n g , W . Q . Q i a n , a n d K . F . He . “ A v a r i a t i o n e v 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 , ” a r X i v : 1 7 0 3 . 1 0 2 6 3 [ c s . S Y ] . [ 1 6 ] S . Z h a n g a n d W . Q . Q i a n , “ C o m p u t a t i o n o f t i m e - o p t i m a l c o n t r o l p r o b l e m w i t h v a r i a t i o n e v o l u t i o n p r i n c i p l e , a r X i v : 1 7 1 1 . 0 2 9 9 8 [ c s . SY ] . [ 1 7 ] S . Z h a n g a n d K . F . He , “ V a r i a t i o n e v o l vi n g 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 , a c o m p a c t w ay , ” a r X i v : 1 7 0 9 . 0 2 2 4 2 [ c s . S Y ] . [ 1 8 ] S . Z h a n g , B . L i a o , a n d F . L i a o , “ C o m p u t a t i o n o f 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 t er m i n a l c o n s t r a i n t v i a v a r i a t i o n e v o l u t i o n , ” a r X i v : 1 8 0 1 . 0 1 3 8 3 [ c s . S Y ] . [ 1 9 ] S . Z h a n g , Y . Q . C h e n g , a n d W . Q . Q i a n , “ O n t h e c o m p u t a t i o n o f o p t i m a l c o n t r o l p r o b l e m s wi t h t er m i n a l i n e q u a li t y c o n s t r a i n t v i a v a r i a t i o n e v o l u t i o n , ” a r X i v : 1 8 0 1 . 0 7 3 9 5 [ c s . S Y ] . [ 2 0 ] S . Z h a n g , E . M . Y o n g , a n d W . Q . Q i a n , “ 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 v i a e v o l u t i o n p a r t i a l d i f f e r e n t i a l e q u a t i o n w i t h a r b i t r a r y d e f i n i t e c o n d i t i o n s , ” a r X i v : 1 7 1 2 . 0 9 7 0 2 [ c s . S Y ] . [ 2 1 ] S. Zh a n g , K . F . H e , a n d F . Li a o , “ C o m p u t a t i o n of o p t i m a l c o n t r o l p r o b l e m s wi t h t er m i n a l c on s t r a i n t vi a m o d i f i e d e v o l u t i o n p a r t i a l di f f er e n t i a l e q u a t i o n , ” a r X i v : 1 8 0 1 . 1 0 4 8 6 [ c s . S Y ] . [ 2 2 ] H . K . Kh a l i l , N o n l i n e a r S y s t e m s . N e w J e r s e y , U S A : Pr e n t i c e H a l l , 2 0 0 2 , p p . 1 1 1 - 181. [ 2 3 ] D . G . Wu , V a r i a t i o n M e t h o d . B e i j i n g , C h i n a : H i g h e r E d u c a t i o n P r e s s , 1 9 8 7 , p p . 1 2 - 68. [ 2 4 ] K. W . C a s s e l , 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 en c e a nd E n g i n e e r i n g . N e w Y o r k : C a m b r i d g e U n i v e r s i t y P r e s s , 2 0 1 3 , p p . 2 8 - 81 [ 2 5 ] H . X. Z h a n g a n d M . Y . S h e n . 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 — F u n d a m e n t a l s a n d Ap p l i c a t i o n s o f F i n i t e Di f f e r e n c e M e t h o d s . B e i j i n g , C h i n a : 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 , 2 0 0 3 , p p . 7 6 - 78. [ 2 6 ] D. Z . Z h e n . L i n e a r S y s t e m T h e o r y . B e i j i n g , C h i n a : T s i n g h u a U n i v e r s i t y P r e s s , 2 0 0 2 , p p . 8 5 - 134. [ 2 7 ] D . G . H u l l , O p t i m a l C o n t r o l T h e o r y f o r A p p l i c a t i o n s , N e w Y o r k , U S A : S p r i n g e r , 2 0 0 3 , p p . 8 9 - 7 1 . [ 2 8 ] A. E. B r y s o n a n d Y . C . H o , 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 . Wa s h i n g t o n , D C , U S A : H e m i s p h e r e , 1975 , p p . 42 - 125. [ 2 9 ] H . W a n g , J . S . L u o , a n d J . M . Z h u . A d v a n c e d M a t h e m a t i c s . C h a n gs h a , C h i n a : 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 , 200 0 , p p . 203 - 210. [ 3 0 ] G . C . Z h a n g , 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 M e t h o d s . C h e n g d u , C h i n a : C h e n g d u S c i e n c e a n d T e c h n o l o g y U n i v e r s i t y P r e s s , 1991 , p p . 8 4 - 86 . [ 3 1 ] M . A . P a t t e r s o n a n d A . V . R a o , “ G P O P S - I I : A M 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 h p - 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 et h o d s a n d 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 T r a n s . M a t h . S o f t w a r e , v o l . 41 , n o . 1, pp. 1 - 3 7 , 2 0 1 4 .