1

Optimal Control Computation via Evolution
Partial Differential Equation with Arbitrary
Definite Conditions
Sheng ZHANG, En-Mi YONG, and Wei-Qi QIAN
(2017.11)
Abstract: The compact Variation Evolving Method (VEM) that originates from the continuous-time dynamics stability theory seeks the
optimal solutions with variation evolution principle. It is further developed to be more flexible in solving the Optimal Control Problems
(OCPs), by relaxing the definite conditions from a feasible solution to an arbitrary one for the derived Evolution Partial Differential
Equation (EPDE). To guarantee the validity, an unconstrained Lyapunov functional that has the same minimum as the original OCP is
constructed, and it ensures the evolution towards the optimal solution from infeasible solutions. With the semi-discrete method, the
EPDE is transformed to the finite-dimensional Initial-value Problem (IVP), and then solved with common Ordinary Differential
Equation (ODE) numerical integration methods. Illustrative examples are presented to show the effectiveness of the proposed method.

Key words: Optimal control, dynamics stability, variation evolution, evolution partial differential equation, initial-value problem

I. INTRODUCTION
Optimal control theory aims to determine the inputs to a dynamic system that optimize a specified performance index while
satisfying constraints on the motion of the system. It is closely related to engineering and has been widely studied [1]. Because of
the complexity, Optimal Control Problems (OCPs) are usually solved with numerical methods. Various numerical methods are
developed and generally they are divided into two classes, namely, the direct methods and the indirect methods [2]. The direct
methods discretize the control or/and state variables to obtain the Nonlinear Programming (NLP) problem, for example, the
widely-used direct shooting method [3] and the classic collocation method [4]. These methods are easy to apply, whereas the
results obtained are usually suboptimal [5], and the optimal may be infinitely approached. The indirect methods transform the OCP
to a Boundary-value Problem (BVP) through the optimality conditions. Typical methods of this type include the well-known
indirect shooting method [2] and the novel symplectic method [6]. Although be more precise, the indirect methods often suffer
from the significant numerical difficulty due to the ill-conditioning of the Hamiltonian dynamics, that is, the stability of costates
dynamics is adverse to that of the states dynamics [7]. The recent development, representatively the Pseudo-spectral (PS) method
[8], blends the two types of methods, as it unifies the NLP and the BVP in a dualization view [9]. Such methods inherit the
advantages of both types and blur their difference.
Theories in the control field often enlighten strategies for the optimal control computation, for example, the non-linear variable
transformation to reduce the variables [10]. Recently, a Variation Evolving Method (VEM), which is enlightened by the states
evolution within the stable continuous-time dynamic system, is proposed for the optimal control computation [11][12][13]. The
The authors are with the Computational Aerodynamics Institution, China Aerodynamics Research and Development Center, Mianyang, 621000, China. (e-mail:
zszhangshengzs@hotmail.com).

2
VEM also synthesizes the direct and indirect methods, but from a new standpoint. The Partial Differential Equation (PDE), which
describes the evolution of variables towards the optimal solution, is derived and the optimality conditions will be satisfied under
this frame. In Ref. [11], the Augmented Evolution PDE (AEPDE) is developed with the employment of the costates. In Ref. [12],
an effective form of the VEM is proposed to obtain the AEPDE for the classic time-optimal control problem with control constraint.
In Ref. [13], a compact VEM that uses only the original variables is developed to reduce the complexity of the computation. The
costate-free optimality conditions are derived and the corresponding Evolution PDE (EPDE) is established. However, the definite
conditions for the EPDE are required to be feasible solutions that satisfy the state equations. In this paper, the compact VEM is
further developed, and a modified EPDE that uses arbitrary definite conditions but still seeks the optimal solution is developed,
which brings extra flexibility for the computation.
Throughout the paper, our work is built upon the assumption that the solution for the optimization problem exists. We do not
describe the existing conditions for the purpose of brevity. Relevant researches such as the Filippov-Cesari theorem are
documented in [14]. In the following, first the principle of the VEM is reviewed. Then the compact VEM is further developed to
establish the modified EPDE that accommodates arbitrary definite conditions. Later illustrative examples are solved to verify the
effectiveness of the proposed method.
II. PRINCIPLE OF VEM
The VEM is a newly developed method for the optimal solutions. It originates from the continuous-time dynamics stability
theory in the control field.
Lemma 1 [15] (with small adaptation): For a continuous-time autonomous dynamic system like
x = f ( x )
where x ∈ \ n is the state, x =

(1)

dx
is its time derivative, and f : \ n → \ n is a vector function. Let x̂ , contained within the
dt

domain D , be an equilibrium point that satisfies f ( xˆ ) = 0 and D . If there exists a continuously differentiable function
V : D → \ such that

i) V ( xˆ ) = c and V ( x ) > c in D / { xˆ } .
ii) V ( x ) ≤ 0 in D and V ( x ) < 0 in D / { xˆ } .
where c is a constant. Then x = xˆ is an asymptotically stable point.

Lemma 1 aims to the dynamic system with finite-dimensional states, and it may be directly generalized to the
infinite-dimensional case as

Lemma 2: For an infinite-dimensional dynamic system described by

δ y ( x)
= f ( y, x)
δt

(2)

∂y ( x, t )
= f ( y, x)
∂t

(3)

or presented equivalently in the PDE form as

where “ δ ” denotes the variation operator and “ ∂ ” denotes the partial differential operator. x ∈ \ is the independent
variable, y ( x) ∈ \ n ( x) is the function of x , and f : \ n ( x) × \ → \ n ( x) is a vector function. Let yˆ ( x) , contained within a

3

certain function set D ( x) , is an equilibrium function that satisfies f ( yˆ ( x), x) = 0 . If there exists a continuously differentiable
functional V : D ( x) → \ such that
i) V ( yˆ ( x) ) = c and V ( y ( x) ) > c in D ( x) / { yˆ ( x)} .
ii) V ( y ( x) ) ≤ 0 in D ( x) and V ( y ( x) ) < 0 in D ( x) / { yˆ ( x)} .
where c is a constant. Then y ( x) = yˆ ( x) is an asymptotically stable solution.
In the system dynamics theory, from the stable dynamics, we may construct a monotonously decreasing function (or functional)
V , which will achieve its minimum when the equilibrium is reached. Inspired by it, now we consider its inverse problem, that is,

from a performance index function to derive the dynamics that minimize this performance index, and optimization problems are
just the right platform for practice. Under this thought, the optimal solution is analogized to the equilibrium of a dynamic system
and is anticipated to be obtained in an asymptotically evolving way. Accordingly, a virtual dimension, the variation time τ , is
introduced to implement the idea that a variable x (t ) evolves to the optimal solution to minimize the performance index within the
dynamics governed by the variation dynamic evolution equations. Fig. 1 illustrates the variation evolution process of the VEM in
solving the OCP. Through the variation motion, the initial guess of variable will evolve to the optimal solution.

Fig. 1. The illustration of the variable evolution along the variation time

τ

in the VEM.

For example, consider the calculus-of-variations problems which may be regarded as OCPs with integrator dynamics.
J = ∫ F ( y (t ), y (t ), t ) dt
tf

t0

(4)

where t ∈ \ is the time. The elements of the variable vector y (t ) ∈ \ n belong to C 2 [t0 , t f ] . t0 and t f are the fixed initial
and terminal time, and the boundary conditions are prescribed as y (t0 ) = y0 and y(t f ) = y f . Within the following variation
dynamic evolution equation

δy
d
⎛
⎞
= − K ⎜ Fy − ( Fy ) ⎟ ，t ∈ (t0 , t f )
t
δτ
d
⎝
⎠
where the column vectors Fy =

(5)

∂F
∂F
and Fy =
are the shorthand notations of partial derivatives, and K is a n × n
∂y
∂y

dimensional positive-definite matrix. Starting from a feasible solution y (t ) that satisfies the boundary conditions, Eq. (5) drives
the performance index J to decrease and the variables y approaches the optimal solution asymptotically. This is because that by
differentiating Eq. (4) with respect to τ (even τ does not explicitly exist) and substituting Eq. (5) in, we have

4
t
δJ
=∫
t
δτ

0

f

δ y ⎞
⎛ Tδy
dt
+ Fy T
⎜ Fy
δτ
δτ ⎟⎠
⎝

= Fy T

t
δy
δy
− Fy T
+∫
t
δτ t
δτ t

0

f

0

f

T
⎛⎡
d
⎤ δ y⎞
⎜ ⎢ Fy − ( Fy ) ⎥
⎟ dt
⎜⎣
dt
⎦ δτ ⎟⎠
⎝

(6)

tf ⎛ ⎡
d
d
⎤
⎡
⎤⎞
= − ∫ ⎜ ⎢ Fy − ( Fy ) ⎥ K ⎢ Fy − ( Fy ) ⎥ ⎟ dt
t0 ⎜ ⎣
dt
dt
⎦
⎣
⎦ ⎟⎠
⎝
≤0
T

where the superscript “ T ” denotes the transpose operator, and when

δJ
= 0, y (t ) will satisfy the optimal conditions, namely, the
δτ

Euler-Lagrange equation [16][17]

Fy −

d
( Fy ) = 0
dt

(7)

Generally, the variation dynamic evolution equations for typical OCPs may be reformulated as the EPDE and the Evolution
Differential Equation (EDE), by replacing the variation operator “ δ ” with the partial differential operator “ ∂ ” and differential
operator “ d ”. For instance, Eq. (5) may be re-presented as
∂y
∂
⎛
⎞
= − K ⎜ Fy − ( Fy ) ⎟
∂τ
∂
t
⎝
⎠

(8)

Since the right function of the EPDE (like Eq. (8)) only depends on the time t , it is suitable to be solved with the well-known
semi-discrete method in the field of PDE numerical calculation [18]. With the discretization along the normal time dimension, the
EPDE is transformed to be IVPs with finite states. Note that the resulting IVP is defined with respect to the variation time τ , not
the normal time t . Then, we may use mature Ordinary Differential Equation (ODE) integration methods to get the numerical
solution.
III. MODIFIED EVOLUTION PARTIAL DIFFERENTIAL EQUATION
As Ref. [13], we still consider the following class of OCP that is defined as

Problem 1：Consider performance index of Bolza form
J = ϕ ( x (t f ), t f ) + ∫ L ( x (t ), u(t ), t ) dt
tf

t0

(9)

subject to the dynamic equation
x = f ( x,u, t )

(10)

where t ∈ \ is the time. x ∈ \ n is the state vector and its elements belong to C 2 [t0 , t f ] . u ∈ \ m is the control vector and its
elements belong to C1[t0 , t f ] . The function L : \ n × \ m × \ → \ and its first-order partial derivatives are continuous with
respect to x , u and t . The function ϕ : \ m × \ → \ and its first-order and second-order partial derivatives are continuous with
respect to x and t . The vector function f : \ n × \ m × \ → \ n and its first-order partial derivatives are continuous and Lipschitz
in x , u and t . The initial time t0 is fixed and the terminal time t f is free. The initial boundary conditions are prescribed as

x (t0 ) = x0
and the terminal states are free. Find the optimal solution ( xˆ , uˆ ) that minimizes J , i.e.

(11)

5
( xˆ , uˆ ) = arg min( J )

(12)

A. Former results
In Ref. [13], by differentiating Eq. (9) and considering the dynamic constraint of Eq. (10) within the feasible solution domain, we
derived the following variation dynamic evolution equations as
t
δx
δu
= ∫ Φo (t , s ) f u ( s )
(s) d s
t
δτ
δτ

(13)

0

δu
= − Kpu
δτ
δtf
δτ

(14)

(

= − kt f L + ϕt + ϕ x T f

)

(15)
tf

where

pu = Lu + fu Tϕ x + fu T

(∫

tf

t

(

)

Φo T (σ , t ) Lx (σ ) + ϕtx (σ ) + ϕ xx T (σ ) f (σ ) + f x (σ )T ϕ x (σ ) d σ

)

(16)

K is the m × m dimensional positive-definite matrix and kt f is a positive constant. Φo (t , s ) is the n × n dimensional state

transition matrix from time point s to time point t , which satisfies
∂
Φo (t , s ) = f x (t )Φo (t , s)
∂t

(17)

Note that the necessary conditions, that Eqs. (13)-(15) seek the optimal solution successfully, require that the initial value of

⎡ x (t ) ⎤
⎡ x (t ) ⎤
variables ⎢
=⎢
⎥
⎥ is feasible, that is, satisfying Eqs. (10) and (11). From such initial variables, Eqs. (13)-(15) will
⎣ u(t ) ⎦ τ =0 ⎣ u (t ) ⎦
evolve the solution to the optimal that satisfies the costate-free optimality conditions, i.e.
pu = 0

(18)

L(t f ) + φt (t f ) + ϕ x T (t f ) f (t f ) = 0

(19)

Use the partial differential operator “ ∂ ” and the differential operator “ d ” to reformulate the variation dynamic evolution
equations (13)-(15), we may get the following EPDE and ED as
∂u
⎡ t
⎤
∂ ⎡ x ⎤ ⎢ ∫t Φo (t , s ) f u ( s ) ( s ) d s ⎥
=
0
∂
τ
⎥
∂τ ⎢⎣ u ⎥⎦ ⎢
− Kpu
⎣⎢
⎦⎥

dt f
dτ

(

= − kt f L + ϕ t + ϕ x T f

with the definite conditions including the initial guess of t f , i.e., t f

τ =0

)

(20)

(21)
tf

⎡ x (t ,τ ) ⎤
⎡ x (t ) ⎤
= tf , and ⎢
=⎢
⎥
⎥ . In this formulation, the
⎣ u(t ,τ ) ⎦ τ =0 ⎣ u (t ) ⎦

initial conditions of x (t ,τ ) and u(t ,τ ) at τ = 0 belong to the feasible solution domain and they will achieve the optimal solution
of the OCP at τ = +∞ .

6
B. Modification of variation dynamic evolution equations
To increase the flexibility of computation, we hope to eliminate the drawback in Ref. [13] that a feasible initial solution of
variables is required. Now the problem of how to drive an infeasible solution to be feasible arises, that is, from a solution that
violates Eqs. (10) and (11) to the one that meets them. Intuitively, we may drive an error variable e(t ) (and a parameter trivially) to
be zero in terms of the following equations

δ e(t )
= − ke(t )
δτ

(22)

where k is a positive constant. Through this idea, we redefine the dynamics constraint as
x − f ( x , u, t ) = e f (t )

(23)

where e f is the dynamics error variable. Correspondingly, the equation describing the variation evolution of Eq. (23) is

δ x
δx
δu δef
= fx
+ fu
+
δτ
δτ
δτ δτ

(24)

To eliminate the dynamics error, we artificially introduce the relation during the solution evolution as

δef
δτ

= −K f e f

(25)

where K f is a positive-definite matrix. Also, the evolution equation for x (t0 ) may be established as

δ x (t0 )
= − K x ex
δτ
0

(26)

0

where K x0 is a positive-definite matrix and e x0 = x (t0 ) − x0 is the initial state error. By the linear system theory [19], we may get
the solution for Eq. (24), which is the adaptation of Eq. (13), as

δef
⎛
δ x (t0 ) t
δ x (t )
δu
= Φo (t , t0 )
+ ∫ Φo (t , s ) ⎜ f u ( s )
( s) +
t
δτ
δτ
δτ
δτ
⎝
0

= −Φo (t , t0 ) K x0 e x0 + ∫

t

t0

⎞
⎟d s
⎠

(27)

t
δu
Φo (t , s ) f u ( s )
( s ) d s − ∫ Φo (t , s ) K f e f d s
t
δτ
0

It is worthwhile to reinvestigate how the performance index given by Eq. (9) varies when starting from an infeasible solution. See

δJ δ
=
δτ δτ

( ∫ (ϕ + ϕ
tf

t0

t

x

T

)

x + L( x, u, t) dt

)

δt f

δef
tf ⎛
δx
δu⎞
+ ∫ ⎜ (ϕtxT + x Tϕxx + ϕxT f x + LxT ) + ϕxT
+ (ϕxT fu + LuT ) ⎟ dt
t0
t f δτ
δτ
δτ
δτ ⎠
⎝
δt f tf ⎛ T δ u
δef
δx ⎞
= (ϕt + ϕxT x + L)
+ ∫ ⎜ pu
+ pf T
+ px0 T 0 ⎟ dt
t0
t f δτ
δτ
δτ
δτ ⎠
⎝
= (ϕt + ϕxT x + L)

(28)

where
pu = Lu + fu Tϕ x + fu T
tf

(∫

tf

t

(

)

Φo T (σ , t ) Lx (σ ) + ϕtx (σ ) + ϕ xx T (σ ) x + f x (σ )T ϕ x (σ ) d σ

(

)

p f = ϕ x (t ) + ∫ Φo T (σ , t ) Lx (σ ) + ϕtx (σ ) + ϕ xx T (σ ) x (σ ) + f x (σ )T ϕ x (σ ) d σ
t

(

p x0 = Φo T (t , t0 ) Lx + ϕtx + ϕ xx T x + f x Tϕ x

)

By investigating Eq. (28), it is found that a better way to reduce the value of J is to adapt Eqs. (14) and (15) as

)

(29)
(30)
(31)

7

δu
= − Kpu
δτ
δtf
δτ

(

= −kt f L + ϕt + ϕ x T x

(32)

)

(33)
tf

However, we still cannot expect J will monotonously decrease, because the term

and the sign of

tf

∫t

0

⎛ T δef
δ x0
+ p x0 T
⎜ pf
δ
τ
δτ
⎝

⎞
⎟ dt may be positive,
⎠

δJ
is uncertain.
δτ

With these treatment, we anticipate that the modified variation dynamic evolution equations (27), (32) and (33) will evolve an
arbitrary initial guess of solutions to the optimal, by achieving the feasibility and optimality simultaneously. However, although the
modification is practically intuitive, lacking the convergence guarantee by Lemma 2 (with Eq. (9) as the Lyapunov functional), one
may argue that
i) Is it true that the optimal solution of Problem 1 is the equilibrium solution of variation dynamic evolution equations (27), (32)
and (33)?
ii) Is it ensured that under the dynamics governed by the evolution equations (27), (32) and (33), the variables will approach the
equilibrium solution from arbitrary initial value, instead of converging to the limit cycle as the Van der Pol oscillator [15]?
For the first question, it is easy to verify that the solution that satisfies Eqs. (10), (11), (18) and (19) is the equilibrium solution of
variation dynamic evolution equations (27), (32) and (33). Now we will answer the second question with rigorous mathematic
argument as follows.
C. Mathematic validation
Before we carry out the mathematic analysis, certain assumption is presented.

Assumption 1: The variables p f (t ) and p x0 (t ) , defined in Eqs. (30) and (31) respectively, are bounded within the time horizon
[t0 , t f ] as

p f (t )

2

≤ d1

(34)

p x0 (t )

2

≤ d2

(35)

where ⋅ 2 denotes the 2-norm of vector.

Lemma 3: For the unconstrained functional
V = e x0 T e x0 + ∫

tf

t0

1
e f T e f dt + c1 J + c2 (e f T e f )
tf
2

(36)

where e x0 = x (t0 ) − x0 , e f = x − f ( x , u, t ) , J is defined in Eq. (9), there exists certain positive constants c1 and c2
⎛
⎞
min(eig( K x0 ))
min(eig( K f ))
⎟
,
c1 < min ⎜
2
⎜ d 2 max(eig( K x ))(t f − t0 ) d1 max(eig( K f ))(t f − t0 ) ⎟
0
⎝
⎠
c2 >

kt f
2c1 min(eig( K f ))

(37)

(38)

such that Eq.(36) is a Lyapunov functional for the variation dynamic evolution equations (27), (32) and (33). Here eig(⋅) is the

8
function of eigenvalue.
Proof: First we show that the minimum solution of Problem 1, denoted by ( xˆ , uˆ ) , is also the minimum solution of the
unconstrained functional (36). If the variables x and u are located within the feasible domain, in which the variables meet Eqs.
(10) and (11), then we have
V = c1 J

(39)

Obviously for this case the minimum of Problem 1 is the minimum of the unconstrained functional (36). When the variables lies in
the infeasible domain, we consider the neighborhood around the minimum solution ( xˆ , uˆ ) . Since ( xˆ , uˆ ) satisfies Eqs. (10), (11),
(18) and (19), we have the first order variation of functional (36) at ( xˆ , uˆ ) as

δV = δ ex

0

2

+∫

tf

δef

t0

tf

tf

t0

t0

dt + c1 ∫ p f Tδ e f dt + c1 ∫ p x0 Tδ e x0 dt

2

(40)

According to Assumptions 1, and with the Holder’s inequality, there is
− d1 (t f − t0 ) ∫

tf

t0

δef

2

− d 2 (t f − t0 ) 2 δ e x0

dt ≤ − ∫

tf

pf

t0

2

≤ −∫

tf

px0

t0

2

2

dt ∫

tf

t0

dt ∫

tf

δ ex

t0

tf

2

dt ≤ ∫ p f Tδ e f dt

2

dt ≤ ∫ px0 Tδ e x0 dt

δef

(41)

t0

tf

0

(42)

t0

Then we have

(

)

δ V ≥ 1 − c1d 2 (t f − t0 )2 δ e x

2

0

(

+ 1 − c1d1 (t f − t0 )

) ∫t

tf

δef

0

2

(43)

dt

According to Eq. (37), we have δ V > 0 . Especially, since c1 may be arbitrarily small, the infeasible domain where the minimum
maintains may be arbitrarily large. In summary, the solution ( xˆ , uˆ ) determines a minimum for the functional (36).
Now we consider the derivative of V with respect to the variation time τ . Differentiating Eq. (36) produces
T

ex
δ x (t0 ) t
δV
=
+∫
t
δτ
δτ
ex Tex
0

e f T e f δτ

0

0

δef

efT

f

0

dt + c1

δef
δtf
δJ
+ c2 (e f T
) + ( e f Te f )
δτ
δτ t
t δτ

(44)

f

f

Substitute Eq. (28) in and use Eqs. (27), (32) and (33), we have
T

ex
t
δV
=−
K x ex − ∫
t
δτ
ex Tex
0

0

0

− c1 ∫

tf

t0

0

efT

f

K f e f d t − c2 ( e f T K f e f )

T

ef ef

0

0

p f T K f e f dt − c1 ∫

tf

t0

(

− kt f c1 ϕ t + ϕ x T x + L

tf

(

p x0 T K x0 e x0 dt − kt f c1 e f T e f L + ϕt + ϕ x T x

)

)

2
tf

− c1 ∫

tf

t0

pu T Kpu dt

(45)

tf

With the Young’s inequality, there is

(

− kt f e f T e f L + ϕt + ϕ x T x

)

kt f c

≤

2

tf

( L + ϕ + ϕ x )
t

x

T

2

+
tf

kt f
2c

( ef

2

(46)

)

2 t
f

According to Assumptions 1, and with the Holder’s inequality, there are
tf

− ∫ p f T K f e f dt ≤ max(eig( K f ))
t0

(∫

tf

t0

tf

pf

2

− ∫ p x0 T K x0 e x0 dt ≤ max(eig( K x0 ))(t f − t0 )
t0

Substituting the inequalities (46)-(48) into Eq. (45) gives

dt

(∫

)∫

tf

t0

tf

t0

p x0

ef

2

2

dt ≤ d1 max(eig( K f ))(t f − t0 ) ∫

tf

t0

)

dt e x 0

2

ef

2

(47)

dt

≤ d 2 max(eig( K x0 ))(t f − t0 ) 2 e x0

2

(48)

9
δV
≤ − min(eig( K x )) − c1d 2 max(eig( K x ))(t f − t0 ) 2 e x
δτ

(

0

kt ⎞
⎛
− ⎜ c2 min(eig( K f )) − f ⎟ ( e f
⎜
2c1 ⎠⎟
⎝

)

0

2

)

tf

2 t
f

0

2

− c1 ∫ pu Kpu dt −
T

t0

With the value of c1 and c2 set by Eqs. (37) and (38), we have that
and (33), and

(

− min(eig( K f )) − c1d1 max(eig( K f ))(t f − t0 )
kt f c1
2

(ϕ

+ ϕ x x + L
T

t

)

) ∫t

tf

ef

0

2

dt

(49)

2
tf

δV
≤ 0 hold under the dynamics governed by Eqs. (27), (32)
δτ

δV
= 0 when Eqs. (10), (11), (18) and (19) are satisfied.
δτ

■

Theorem 1: Solving the IVP with respect to τ , defined by the variation dynamic evolution equations (27), (32) and (33) with any
initial solution, when τ → +∞ , ( x, u) will satisfy the feasibility conditions and the optimality conditions of Problem 1.
Proof: The proof is a direct application of Lemmas 2 and 3. From Lemmas 3, the functional (36) is ensured a Lyapunov
functional for the dynamics system (27), (32) and (33) around the equilibrium that meets (10), (11), (18) and (19). According to
Lemma 2, the equilibrium solution is an asymptotically stable solution and ( x, u) will satisfy the feasibility conditions (10), (11)
and the optimality conditions (18), (19) of Problem 1 when τ → +∞ .

■

D. Formulation of modified EPDE
Similarly, we may use the partial differential operator “ ∂ ” and the differential operator “ d ” to reformulate the variation
dynamic evolution equations to get the modified EPDE and EDE as
t
t
⎡
∂u
⎤
⎛ ∂x
⎞
⎢ −Φo (t , t0 ) K x0 ( x (t0 ) − x0 ) + ∫t0 Φo (t , s) f u ( s ) ∂τ ( s ) d s − ∫t0 Φo (t , s ) K f ⎜ ∂t ( s ) − f ( s ) ⎟ d s ⎥
⎝
⎠
∂ ⎡ x⎤ ⎢
⎥
=
⎢
⎥
⎢
⎧
∂τ ⎣ u ⎦
∂x
⎞ ⎫⎥
⎛
⎞
T
T ⎛ tf
T
T
T
⎢ − K ⎨ Lu + fu ϕ x + fu ⎜ ∫ Φo (σ , t ) ⎜ Lx (σ ) + ϕtx (σ ) + ϕ xx (σ ) (σ ) + f x (σ ) ϕ x (σ ) ⎟ d σ ⎟ ⎬ ⎥
∂t
⎝
⎠
⎢⎣
⎝ t
⎠ ⎭ ⎦⎥
⎩

∂x ⎞
⎛
= − kt f ⎜ L + ϕt + ϕ x T
⎟
∂t ⎠ t
dτ
⎝
f

dt f

In particular, the definite conditions including t f

τ =0

(50)

(51)

⎡ x (t ,τ ) ⎤
⎡ x (t ) ⎤
= tf and ⎢
=⎢
⎥
⎥ may be infeasible solutions. Recall the
τ
(
t
,
)
u
⎣
⎦ τ =0 ⎣ u (t ) ⎦

anticipated variable evolution along the variation time τ illustrated in Fig. 1, the initial conditions of x (t ,τ ) and u(t ,τ ) at τ = 0
may be arbitrary and their value at τ = +∞ will be the optimal solution of the OCP. The right part of the EPDE (50) is also only a
vector function of time t . Thus we may apply the semi-discrete method to discretize it along the normal time dimension and
further use ODE integration methods to get the numerical solution.
IV. ILLUSTRATIVE EXAMPLES

First a linear example taken from Xie [20] is considered.
Example 1: Consider the following dynamic system
x = Ax + bu

⎡x ⎤
⎡0⎤
⎡0 1 ⎤
, b = ⎢ ⎥ . Find the solution that minimizes the performance index
where x = ⎢ 1 ⎥ , A = ⎢
⎥
⎣1 ⎦
⎣0 0⎦
⎣ x2 ⎦

10
J=

(

)

1
1 tf
x (t f )T Fx (t f ) + ∫ x T Qx + Ru 2 dt
2
2 t0

⎡1⎤
with the initial boundary conditions x (t0 ) = ⎢ ⎥ , where the initial time t0 = 0 and the terminal time t f = 3 are fixed. The
⎣1⎦
⎡1 0 ⎤
⎡2 1⎤
1
weighted matrixes are F = ⎢
⎥ , Q = ⎢1 4 ⎥ and R = 2 .
0
2
⎣
⎦
⎣
⎦
In solving this example using the VEM, the EPDE derived is
t
t
⎡ A ( t − t0 )
∂u
⎛ ∂x
⎞ ⎤
K x0 ( x (t0 ) − x0 ) + ∫ e A(t − s ) b ( s) d s − ∫ e A(t − s ) K f ⎜ ( s) − Ax ( s) − bu ( s) ⎟ d s ⎥
⎢ −e
t
t
0
0
∂
∂
τ
t
⎝
⎠ ⎥
∂ ⎡ x⎤ ⎢
=
⎥
T
t
⎧
⎫
∂τ ⎣⎢ u ⎦⎥ ⎢
x
⎧
∂
⎫
⎛
⎞
f
(σ ) + AT Fx (σ ) ⎟ d σ ⎬⎬
− K ⎨ Ru + bT Fx + bT ⎨ ∫ e A(σ −t ) ⎜ Qx (σ ) + F
⎢
⎥
t
∂t
⎝
⎠ ⎭⎭
⎢⎣
⎥⎦
⎩
⎩

(

)

⎡ 0.1 0 ⎤
The one-dimensional matrix K is K = 2 × 10−2 . The 2 × 2 dimensional matrixes K x0 and K f are both ⎢
⎥ . The definite
⎣ 0 0.1⎦

conditions of the EPDE, i.e., the initial guess of the states x (t ) and the control u (t ) , were directly set as x (t ) = 0 and u (t ) = 0 .
Using the semi-discrete method, the time horizon [t0 , t f ] was discretized uniformly, with 61 points. Thus, a dynamic system with
183 states was obtained and the OCP was transformed to an IVP. Regarding the time derivative

∂x
, it is computed with the finite
∂t

difference method at the discretization points. The ODE integrator “ode45” in Matlab, with default relative error tolerance 1×10-3
and default absolute error tolerance 1×10-6, was employed to solve the IVP. Even for this simple example, there is no analytic
solution. For comparison, we computed the optimal solution with GPOPS-II [21], a Radau PS method based OCP solver.
Figs. 2, 3 and 4 show the evolving process of x1 (t ) , x2 (t ) and u (t ) solutions to the optimal, respectively. At τ = 300s, the
numerical solutions are indistinguishable from the optimal, and this shows the effectiveness of the VEM. It is interesting to
investigate the performance index profile for the variable evolution starting from an infeasible solution. Fig. 5 plots the curve of
performance index value against the variation time. Since beginning from a zero solution, its value first increases from zero. When
it gradually turns to be feasible, it declines to approach the minimum of the OCP.

The optimal solution
Numerical solutions with VEM

1.2
1

τ = 300s
τ = 32.3s

x1(t)

0.8

τ = 15.3s

0.6

τ = 3.3s

0.4

τ = 1.6s

0.2
τ = 0s

0
-0.2
0

0.5

1

1.5

2

2.5

3

t
Fig. 2 The evolution of numerical solutions of x1 to the optimal solution.

11

The optimal solution
Numerical solutions with VEM

1.2
1
0.8

x2(t)

0.6

τ = 0s

0.4
τ = 1.6s

0.2

τ = 4.6s

0
τ = 15.3s

-0.2

τ = 32.3s

-0.4

τ = 300s

0

0.5

1

1.5

2

2.5

3

t
Fig. 3 The evolution of numerical solutions of

x2

to the optimal solution.

1
τ = 0s

0

u (t)

-1
-2

τ = 4.6s
τ = 15.3s

τ = 32.3s

-3
τ = 59.3s

-4

τ = 300s

The optimal solution
Numerical solutions with VEM

-5
-6
-0.5

0

0.5

1

1.5

2

2.5

3

t
Fig. 4 The evolution of numerical solutions of

u

to the optimal solution.

3.5
3
2.5

J

2
1.5
1

Minimum of J
Funcitonal value during the evolution

0.5
0
0

50

100

150

τ (s)

200

250

300

Fig. 5 The evolution profile of the performance index value

Now we consider a nonlinear example with free terminal time t f , the homing missile problem adapted from Hull [22].
Example 2: Consider the problem of a constant speed missile intercepting a constant speed target moving in a straight line. The

dynamic equations are

12
x = f ( x , u )
⎡ VT sin(θT ) − VM sin(θ M ) ⎤
⎡ x ⎤
⎢V cos(θ ) − V cos(θ ) ⎥
T
M
M ⎥
where x = ⎢⎢ y ⎥⎥ , and f = ⎢ T
. x and y are the relative abscissa and ordinate, respectively. θ M is the
⎢
⎥
u
⎢⎣θ M ⎥⎦
⎢
⎥
VM
⎣
⎦

azimuth of missile. θT =30 deg is the azimuth of the target. VM =1000 m/s is the constant speed of missile. VT =500 m/s is the
constant speed of the target. u is the missile normal acceleration. To intercept the target and penalize too large normal acceleration,
the performance index to be minimized is defined as
J=
⎡1× 10−2
⎢
where the weighted matrixes are F = ⎢ 0
⎢ 0
⎢⎣

1
1 tf
x (t f )T Fx (t f ) + ∫ Ru 2 dt
2
2 t0
0
2 × 10−2
0

0⎤
⎥
0 ⎥ and R = 5 × 10−4 . The initial boundary conditions are
0 ⎥⎥
⎦

⎡ x ⎤
⎡10000m ⎤
⎢ y ⎥
= ⎢⎢ 5000m ⎥⎥ and the terminal time t f is free.
⎢ ⎥
⎢⎣θ M ⎥⎦
⎢⎣ 0deg ⎥⎦
t0 = 0
Before the computation, the states and the control variables were scaled to improve the numerical efficiency. In the specific form
of the EPDE (50) and the EDE (51), the parameters K and kt f were set to be 2 × 10−6 and 2 × 10−4 , respectively. The 3 × 3

dimensional matrixes K x0 and K f

0⎤
⎡ 0.1 0
⎢
were both ⎢ 0 0.1 0 ⎥⎥ . The definite conditions, i.e.,
⎢⎣ 0
0 0.1⎥⎦

⎡ x (t ,τ ) ⎤
⎢
⎥
, were directly set to be
⎢ u(t ,τ ) ⎥
⎢ t f (τ ) ⎥
⎣
⎦ τ =0

x (t ) = 0 and u (t ) = 0 with tf = 25s . We also discretized the time horizon [t0 , t f ] uniformly, with 51 points. Thus, a large IVP
with 205 states (including the terminal time) is obtained. We still employed “ode45” in Matlab for the numerical integration. In the
integrator setting, the default relative error tolerance and the absolute error tolerance are 1×10-3 and 1×10-6, respectively. For
comparison, the optimal solution is again computed with GPOPS-II.
Fig. 6 gives the states curve in the x y relative coordinate plane, showing that the numerical results approach the optimal
solution and gradually meet the boundary conditions over time. Note that the solution at τ = 0s is just a point located at the origin.
For the optimal solution, the missile will intercept the target with a fairly small position error. The control solutions are plotted in
Fig. 7, and the asymptotical approach of the numerical results are demonstrated. In Fig. 8 the terminal time profile against the
variation time τ is plotted. The result of t f oscillates at first and then gradually approaches to the optimal interception time, and it
only changes slightly after τ = 100s. At τ = 300s, we compute that t f = 23.51s from the VEM, very close to the result of 23.52s
from GPOPS-II.

13

6000
5000

τ = 300s

y(m)

4000
τ = 15.8s
τ = 12.6s

3000

τ = 9.9s
τ = 8.1s

2000
1000

The optimal solution
Numerical solutions with VEM

0
-1000
-2000

0

2000

4000

6000

8000

10000

x(m)
Fig. 6 The evolution of numerical solutions in the

xy

relative coordinate plane to the optimal solution.

70
The optimal solution
Numerical solutions with VEM

60
50
τ = 300s

u (t)

40
30
20

τ = 12.6s
τ = 8.1s
τ = 4.9s

10
τ = 2.6s

0
τ = 0s

-10
0

5

10

15

20

25

t
Fig. 7. The evolution of numerical solutions of

u

to the optimal solution.

25
Minimum decline time
Evolving profile of tf

tf (s)

24.5

24

23.5

23

22.5
0

50

100

150

τ (s)

200

250

300

Fig. 8 The evolving profile of t f to the optimal interception time.

14

V. CONCLUSION
The Variation Evolving Method (VEM) is further developed to be more flexible in solving the Optimal Control Problems
(OCPs). In computing the Evolution Partial Differential Equation (EPDE) to seek the optimal solution, the requirement for a
feasible definite condition is relaxed to an arbitrary one, and then the transformed Initial-value Problems (IVPs) may be initialized
with any initial values of variables. During the mathematic validation on the modified evolution equations, an unconstrained
Lyapunov functional that has the same minimum as the original OCP is constructed. However, this unconstrained functional only
has theoretical meaning, and it is not practical for the computation of the OCP.
REFERENCES
[1] H. J. Pesch and M. Plail, “The maximum principle of optimal control: A history of ingenious ideas and missed opportunities,” Control &
Cybernetics, vol. 38, no. 4, pp. 973-995, 2009.
[2] J. T. Betts, “Survey of numerical methods for trajectory optimization,” J. Guid. Control Dynam., vol. 21, no. 2, pp. 193-206, 1998.
[3] Q. Lin, R. Loxton, and K. L. Teo, "The control parameterization method for nonlinear optimal control: a survey," Journal of Industrial and
Management Optimization, vol. 10, no. 1, pp. 275- 309, 2014
[4] C. Hargraves and W. Paris, “Direct trajectory optimization using nonlinear programming and collocation,” J. Guid. Control Dynam., vol. 10,
no. 4, pp. 338-342, 1987.
[5] O. V. Stryk and R. Bulirsch, “Direct and indirect methods for trajectory optimization,” Ann. Oper. Res., vol. 37, no. 1, pp. 357-373, 1992.
[6] H. J. Peng, Q. Gao, Z. G. Wu, and W. X. Zhong, “Symplectic approaches for solving two-point boundary-value problems,” J. Guid. Control
Dynam., vol. 35, no. 2, pp. 653-658, 2012.
[7] A. V. Rao, “A survey of numerical methods for optimal control,” in Proc. AAS/AIAA Astrodynam. Specialist Conf., Pittsburgh, PA, 2009, AAS
Paper 09-334.
[8] D. Garg, M. A. Patterson, W. W. Hager, A. V. Rao, et al, A Unified framework for the numerical solution of optimal control problems using
pseudospectral methods,” Automatica, vol. 46, no. 11, pp. 1843-1851, 2010.
[9] I. M. Ross and F. Fahroo, “A perspective on methods for trajectory optimization,” in Proc. AIAA/AAS Astrodynam. Conf., Monterey, CA, 2002,
AIAA Paper No. 2002-4727.
[10] I. M. Ross and F. Fahroo, “Pseudospectral methods for optimal motion planning of differentially flat systems,” IEEE Trans. Autom. Control,
vol. 49, no. 8, pp. 1410-1413, 2004.
[11] S. Zhang, E. M. Yong, W. Q. Qian, and K. F. He, “A Variation Evolving Method for Optimal Control,” arXiv: 1703.10263 [].
[12] S. Zhang, and W. Q. Qian, “Computation of time-optimal control problem with variation evolution principle,” arXiv:1711.02998 [].
[13] S. Zhang, and K. F. He, “Variation Evolving for Optimal Control Computation, A Compact Way,” arXiv:1709.02242 [].
[14] R. F. Hartl, S. P. Sethi, and R. G. Vickson, “A survey of the maximum principles for optimal control problems with state constraint,” SIAM
Rev., vol. 37, no. 2, pp. 181-218, 1995.
[15] H. K. Khalil, Nonlinear Systems. New Jersey, USA: Prentice Hall, 2002, pp. 111-181.
[16] D. G. Wu, Variation Method. Beijing, China: Higher Education Press, 1987, pp. 12-68.
[17] K. W. Cassel, Variational Methods with Applications in Science and Engineering. New York: Cambridge University Press, 2013, pp. 28-81
[18] H. X. Zhang and M. Y. Shen. Computational Fluid Dynamics — Fundamentals and Applications of Finite Difference Methods. Beijing, China:
National Defense Industry Press, 2003, pp. 76-78.
[19] D.Z. Zhen. Linear System Theory. Beijing, China: Tsinghua University Press, 2002, pp. 85-134.

15
[20] X. S. Xie, Optimal Control Theory and Application. Beijing, China: Tsinghua Univ. Press, 1986, pp. 167-215.
[21] M. A. Patterson and A. V. Rao, “GPOPS-II: AMATLAB software for solving multiple-phase optimal control problems using hp-adaptive
Gaussian quadrature collocation methods and sparse nonlinear programming,” ACM Trans. Math. Software, vol. 41, no. 1, pp. 1-37, 2014.
[22] D. G. Hull, Optimal Control Theory for Applications, New York, USA: Springer, 2003, pp. 89-71.

