1

Variation Evolving for Optimal Control
Computation, A Compact Way
Sheng ZHANG and Kai-Feng HE
(2017.08)
Abstract: A compact version of the Variation Evolving Method (VEM) is developed for the optimal control computation. It follows the
idea that originates from the continuous-time dynamics stability theory in the control field. The optimal solution is analogized to the
equilibrium point of a dynamic system and is anticipated to be obtained in an asymptotically evolving way. With the introduction of a
virtual dimension——the variation time, the Evolution Partial Differential Equation (EPDE), which describes the variation motion
towards the optimal solution, is deduced from the Optimal Control Problem (OCP), and the equivalent optimality conditions with no
employment of costates are established. In particular, it is found that theoretically the analytic feedback optimal control law does not
exist for general OCPs because the optimal control is related to the future state. Since the derived EPDE is suitable to be solved with the
semi-discrete method in the field of PDE numerical calculation, the resulting Initial-value Problems (IVPs) may be solved with mature
Ordinary Differential Equation (ODE) numerical integration methods.

Key words: Optimal control, dynamics stability, variation evolution, initial-value problem, costate-free optimality condition.

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 [2] and the classic collocation method [3]. These methods are easy to apply, whereas the
results obtained are usually suboptimal [4], 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 [5]. 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 [6]. The recent development, representatively the Pseudo-spectral (PS) method
[7], blends the two types of methods, as it unifies the NLP and the BVP in a dualization view [8]. 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 [9]. Recently, a new Variation Evolving Method (VEM), which is enlightened by the states
The authors are with the Computational Aerodynamics Institution, China Aerodynamics Research and Development Center, Mianyang, 621000, China. (e-mail:
zszhangshengzs@hotmail.com).

2
evolution within the stable continuous-time dynamic system, is proposed for the optimal control computation [10]. The 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 extremal solution, is derived from the viewpoint of variation motion in typical OCPs. Using
the well-known semi-discrete method in the field of PDE numerical calculation [11], the PDEs are transformed to the
finite-dimensional Initial-value Problems (IVPs) to be solved, with the mature Ordinary Differential Equation (ODE) integration
methods. Because the extremum is guaranteed be the equilibrium point of the deduced dynamic system, the optimal solution will
be gradually approached. However, in the work of Ref. [10], besides the states and the controls, the costates are also introduced,
which increases the complexity of the computation. In this paper, a compact version of the VEM that uses only the original
variables is developed. The corresponding variation dynamic evolution equations, which may be re-presented in the PDE
formulation, are derived.
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 [12]. In the following, first preliminaries that state the inspiration of the VEM are presented. Then the foundational
VEM bred under this idea is demonstrated for the unconstrained calculus-of-variations problem. Next the compact VEM to solve
the OCPs is established. During this course, the costate-free optimality conditions are derived, and it is proved equivalent to the
traditional conditions. Later illustrative examples are solved to verify the effectiveness of the method. Besides, comparison
between the evolution equation and that in Ref. [10] is presented at the end.
II. PRELIMINARIES
The VEM is a newly developed method for the optimal solutions. For a better understanding, its motivations are recalled. 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. Suppose that x̂ is a
dt

asymptotically stable equilibrium point of system (1) that satisfies f ( xˆ ) = 0 , then from any initial condition x (t ) t = 0 = x0 within
the stability domain D that contains x̂ , the state x will tend to x̂ over time t [13]. According to the Lyapunov theory, there is a
continuously differentiable function V : D → \ such that
i) V ( xˆ ) = 0 and V ( x ) > 0 in D / { xˆ } .
ii) V ( x ) ≤ 0
For example, maybe f ( x ) satisfies ( x − xˆ )T f ( x ) ≤ 0 , and then a feasible Lyapunov function can be constructed as
V=

1
( x − xˆ )T ( x − xˆ )
2

(2)

where the superscript “ T ” denotes the transpose operator. The dynamics given by f ( x ) determines that V ≤ 0 and x will
converge to the equilibrium x̂ . Fig. 1 sketches the trajectory of some state in the stable dynamic system and the corresponding
Lyapunov function value. No matter what the initial condition x0 is, as long as it falls into the stability domain D , the state x
will approaches the equilibrium x̂ gradually.

3

Fig. 1. The sketch for the state trajectory and the Lyapunov function value profile.

In the system dynamics theory, from the stable dynamics of state x , we may construct a monotonously decreasing function

V ( x ) , which will achieve its minimum when x reaches x̂ . 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. Consider the parameter optimization
problem with performance index
J = h(θ )

(3)

where θ is the optimization parameter vector and h : \ n → \ is a scalar function. To find the optimal value θ̂ that minimizes J ,
we make the analogy to the Lyapunov function and differentiate J , i.e., function h here, with respect to a virtual time τ , which is
used to describe the derived dynamics.
dJ dh
∂h dθ
=
= ( )T
dτ dτ
∂θ dτ
To guarantee that J decreases with respect to τ , i.e.,

(4)

δJ
≤ 0 , we may set
δτ
dθ
∂h
= −Kh
dτ
∂θ

(5)

where K h is a positive-definite matrix. Under this dynamics, h will decrease until it reaches a extremum, and θ will approaches
θ̂ , the equilibrium point of system (5), which satisfies

∂h
∂θ

= 0 . This equilibrium condition is exactly the first-order optimality
θ =θˆ

condition for the optimization problem (3).
A bolder idea further arises hereafter. If the optimization parameter θ can approach its optimal under the dynamics given by Eq.
(5), we can imagine a variable x (t ) might also evolve to the optimal solution to minimize some performance index within certain
dynamics. Fig. 2 illustrates the idea of the VEM in solving the OCP, and we formally introduce the virtual variation time τ , a new
dimension orthogonal to the normal time t , to describe the variation evolution process. Through the variation motion, the initial
guess of variable will evolve to the optimal solution. We will detail the implementation of the idea in the following.

4

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

τ

in the VEM.

III. THE FOUNDATIONAL VARIATION EVOLVING METHOD

Calculus-of-variations problems may be regarded as OCPs with integrator dynamics. To start with, the foundational VEM,
which was first demonstrated in Ref. [10], is again presented for the unconstrained calculus-of-variations problem defined as
Problem 1: For the following functional depending on variable vector y (t ) ∈ \ n
J = ∫ F ( y (t ), y (t ), t ) dt
tf

(6)

t0

where t ∈ \ is the time. The elements of y belong to C 2 [t0 , t f ] , which denotes the set of variables with continuous second-order
derivatives (indicated by the superscript). The function F : \ n × \ n × \ → \ and its first-order and second-order partial
dy
derivatives are continuous with respect to y , its time derivative y =
and t . t0 and t f are the fixed initial and terminal time,
dt
and the boundary conditions are prescribed as y (t0 ) = y0 and y(t f ) = y f . Find the extremal solution ŷ that minimizes J , i.e.

yˆ = arg min( J )

(7)

Follow the idea of dynamics evolution to reduce some performance index, we anticipate that any initial guess of y (t ) , whose
elements belong to C 2 [t0 , t f ] , will evolve to the minimum along the variation dimension. Like the decrease of a Lyapunov

function, if J in Eq. (6) decreases with respect to the variation time τ , i.e.,

δJ
≤ 0 , we may finally obtain the optimal solution.
δτ

Differentiating (6) with respect to τ (even τ does not explicitly exist) produces
t
δJ
=
δτ ∫t

0

f

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

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

0

f

where the column vectors Fy =

0

f

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

(8)

∂F
∂F
and Fy =
are the shorthand notations of partial derivatives. From an initial guess y (t )
∂y
∂y

that satisfies the boundary conditions at t0 and t f , then by enforcing

δJ
≤ 0 , we may set that
δτ

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

(9)

5
where K is a n × n dimensional positive-definite matrix. The variation dynamic evolution equation (9) describes the variation
motion of y (t ) starting from y (t ) , and it is proved that the motion is directed to the extremum [10]. It drives the performance
index J to decrease until

δJ
δJ
= 0, this determines the optimal conditions, namely, the Euler-Lagrange
= 0 , and when
δτ
δτ

equation [14][15]
Fy −

d
( Fy ) = 0
dt

(10)

The variation dynamic evolution equation (9) may be considered from the view of PDE formulation, by replacing the variation
operation “ δ ” and the differential operator “ d ” with the partial differential operator “ ∂ ” as
∂y
∂
⎛
⎞
= − K ⎜ Fy − ( Fy ) ⎟
∂τ
∂t
⎝
⎠

(11)

For this PDE, its right function only depends on the time t . Thus it is suitable to be solved with the semi-discrete method in the
field of PDE numerical calculation. With the discretization along the normal time dimension, Eq. (11) 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 . In the previous
work [10], a demonstrative example is solved to verify the result.

IV. THE COMPACT VARIATION EVOLVING METHOD
A. Problem definition
In this paper, we consider the following class of OCP that is defined as

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

t0

(12)

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

(13)

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

(14)

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

(15)

B. Derivation of variation dynamic evolution equations
In Ref. [10], the problem is addressed by constructing an equivalent unconstrained functional problem that has the same
extremum. This operation is practical but it introduces the costate vector, which has the same dimension as the state vector. To
avoid the resulting complexity, here we will address Problem 2 directly. Consider the problem within the feasible solution domain

6

D o , in which any solution satisfies Eqs. (13) and (14). First we transform the Bolza performance index to the equivalent Lagrange
type, i.e.

J =∫

tf

t0

(ϕ + ϕ
t

x

T

)

f ( x,u, t ) + L( x , u, t ) dt

(16)

where ϕt and ϕ x are the partial derivatives notated as before. Similarly we differentiate Eq. (16) with respect to the variation time

τ to obtain
δtf
t
δJ
= (ϕt + ϕ x T f + L)
+∫
t
t δτ
δτ

f

0

f

δu⎞
⎛
T
T
T
T δx
dt
+ (ϕ x T fu + Lu T )
⎜ (ϕtx + f ϕ xx + ϕ x f x + Lx )
δτ
δτ ⎟⎠
⎝

(17)

where ϕtx and ϕ xx are second-order partial derivatives in the form of (column) vector and matrix, and f x and f u are the Jacobi
matrixes. Different from the calculus-of-variations problems,

δx
δu
and
are related because the profiles of x are determined by
δτ
δτ

u . For the solutions in D o , they need to satisfies the following variation equation

δ x
δx
δu
= fx
+ fu
δτ
δτ
δτ
with the initial condition

(18)

δx
= 0 . Note that f x and f u are time-dependent matrixes linearized at the feasible solution x (t ) and
δτ t
0

u(t ) . Eq. (18) is linear with respect to the variables

δx
δu
and
, and has a zero initial value. Thus it satisfies the superposition
δτ
δτ

principle, and its solution may be explicitly expressed.

Lemma 1 [16]: For the linear system

x = A(t ) x + B (t )u

(19)

with the initial value x (t0 ) = 0 , where x ∈ \ n is the state vector, u ∈ \ m is the control vector, A(t ) and B (t ) are the right
dimensional coefficient matrixes, the solution of system (19) is
t

x (t ) = ∫ H (t , s ) u( s ) d s
t0

(20)

where is the n × m dimensional H (t , s) is the impulse response function that satisfies

⎧Φ (t , s) B ( s )
H (t , s) = ⎨
⎩0

t≥s
s <t

(21)

and Φ (t , s) is the n × n dimensional state transition matrix for the system from time point s to time point t .

According to Lemma 1, it may be derived that Eq. (18) has the solution
t
δx
δu
= ∫ H o (t , s )
(s) d s
t
δτ
δτ

(22)

0

where H o (t , s) is the impulse response function corresponding to the specific f x (t ) and f u (t ) . Substitute Eq. (22) into Eq. (17),
there is

7

δtf
t ⎧
δJ
δu
δ u⎫
⎛ t
⎞
( s ) d s ⎟ + (φ x T f u + Lu T ) ⎬ dt
= (ϕt + ϕ x T f + L)
+ ∫ ⎨ ϕtx T + f Tϕ xx + ϕ x T f x + Lx T ⎜ ∫ H o (t , s )
t
t
t
δτ
δτ
δτ
δτ ⎭
⎝
⎠
⎩
δtf
t
t ⎧
t ⎛
δu ⎫
δu⎞
( s ) ⎬ d s d t + ∫ ⎜ (ϕ x T f u + Lu T )
dt
= (ϕt + ϕ x T f + L)
+ ∫ ∫ ⎨ ϕtx T + f Tϕ xx + ϕ x T f x + Lx T H o (t , s )
t
t ⎩
t
t δτ
δτ
δτ ⎠⎟
⎭
⎝
f

0

f

(

)

(

f

0

f

0

0

)

(23)

f

0

By exchanging the order in the double integral, we may derive the following transformation as

δu ⎫
(s)⎬ d s d t
δτ
⎭
δu
T
T
T
T
( s) d s
ϕtx (t ) + f (t )ϕ xx (t ) + ϕ x (t ) f x (t ) + Lx (t ) H o (t , s) d t
δτ

∫t ∫t ⎨⎩(ϕtx
=∫

tf

t

0

0

tf

t0

⎧

(∫

tf

s

(

T

)

(t ) + f T (t )ϕ xx (t ) + ϕ x T (t ) f x (t ) + Lx T (t ) H o (t , s)

(24)

)

)

To make the result clear, the change of variable symbol t → σ and s → t , without changing the final result of Eq. (24), produce
tf
tf
δu
T
T
T
T
∫t0 ∫s ϕtx (t ) + f (t )ϕ xx (t ) + ϕ x (t ) f x (t ) + Lx (t ) H o (t , s) d t δτ (s) d s
(25)
tf
tf
δu
= ∫ ∫ ϕtx T (σ ) + f T (σ )ϕ xx (σ ) + ϕ x T (σ ) f x (σ ) + Lx T (σ ) H o (σ , t ) d σ
(t ) d t

(
(

t0

(

)

)

(

t

) δτ

)

Thus, Eq. (23) may be reformulated as
δtf
δJ
= (ϕt + ϕ x T f + L)

δτ

+∫

tf

t0

Now to achieve

{( ∫

tf

tf

t

(ϕ

tx

T

δτ

)

)

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

T

T

T

δJ
≤ 0 , we may derive the variation dynamic evolution equation for the control u as
δτ
t
δu
= − K Lu + f u Tϕ x + ∫ H o T (σ , t ) Lx (σ ) + ϕtx (σ ) + ϕ xx T (σ ) f (σ ) + f x (σ )T ϕ x (σ ) d σ
t
δτ

{

(

(

f

}

δu
(t ) d t
fu + Lu )
δτ
T

)

)}

(26)

(27)

and for the terminal time as

δtf
δτ

(

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

)

(28)
tf

where K is the m × m dimensional positive-definite matrix and kt f is a positive constant. Presume that we already have a feasible
initial solution x (t ) and u (t ) , then Eqs. (22), (27) and (28) drive it to the one that satisfies
Lu + fu Tϕ x +

(∫

tf

t

(

)

)

H o T (σ , t ) Lx (σ ) + ϕtx (σ ) + ϕ xx T (σ ) f (σ ) + f x (σ )T ϕ x (σ ) d σ = 0

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

(29)
(30)

C. Equivalence to the classic optimality conditions
Actually, Eqs. (29) and (30) are the first-order optimality conditions for Problem 2 without the employment of costates. We will
show that they are equivalent to the traditional ones with costates [17]

λ + H x = λ + Lx + f x T λ = 0
H u = Lu + f u T λ = 0

(31)
(32)

and transversality conditions

λ(t f ) = ϕ x (t f )

(33)

H (t f ) + ϕt f = 0

(34)

8
where H = L + λT f is the Hamiltonian and λ is the costate vector.

Theorem 1: For Problem 2, the optimal conditions given by Eqs. (29) and (30) are equivalent to the optimality conditions given by
(31)-(34).
Proof: With Eq. (21), Eq. (29) may be reformulated as
Lu + fu Tϕ x + f u T

(∫

tf

t

(

)

)

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

where Φo (σ , t ) is the state transition matrix related with H o (σ , t ) . Define a quantity γ (t ) as

(

tf

(35)

)

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

(36)

Then Eq. (35) is simplified as
Lu + fu T γ = 0

(37)

γ (t f ) = φ x (t f )

(38)

Obviously, from Eq. (36), when t = t f , there is

Differentiate γ (t ) with respect to t . In the process, we will use the Leibniz rule [18]
d
dt

(∫

a (t )

b (t )

)

h(σ , t ) d σ = h ( a(t ), t )

a (t )
d
d
a(t ) − h ( b(t ), t ) b(t ) + ∫ ht (σ , t ) d σ
(t )
b
dt
dt

(39)

and the property of Φo (σ , t ) [16]
∂Φo (σ , t )
= −Φo (σ , t ) f x (t )
∂t

(40)

Φo (t , t ) = 1

(41)

where 1 is the n × n dimensional identity matrix. Then we have

(

)

(

)

tf
d
γ (t ) = ϕtx + ϕ xx T f − Lx + ϕtx + ϕ xx T f + f x Tϕ x − f x T ∫ Φo T (σ , t ) Lx (σ ) + ϕtx (σ ) + ϕ xx T (σ ) f (σ ) + f x (σ )T ϕ x (σ ) d σ
t
dt

(

tf

(

)

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

)

(42)

= − Lx − f x T γ (t )

From Eqs. (42) and (38), it is found that γ (t ) conforms to the same dynamics and boundary conditions as the costates λ(t ) . Thus,
we can conclude that γ (t ) = λ(t ) . Then Eq. (37) and Eq. (32) are identical. This means that Eq. (29) is equivalent to Eqs. (31), (32)
and (33).
On the other hand, with Eq. (33), it is readily to show that Eqs. (30) and (34) are same.

■

By investigating the optimality condition (29), it is found that the optimal control are related to the future state, and this means
the optimal feedback control law in the analytic form does not exists for general OCPs. However, optimal control aided with
numerical computation is still possible. After proving the equivalence of optimality conditions, now the variables evolving
direction using the VEM is easy to determine.

Theorem 2: Solving the IVP with respect to τ , defined by the variation dynamic evolution equations (22), (27) and (28), from a
feasible initial solution, when τ → +∞ , ( x, u) will satisfy the optimality conditions of Problem 2.

9
Proof: From a feasible initial solution, any evolution under the dynamics of Eq. (22) guarantees the feasibility of solution.
Substitute Eqs. (27) and (28) into Eq. (26), we have

τ → +∞ due to the asymptotical approach. When

δJ
δJ
≤ 0 . The functional J will decrease until
= 0 , which occurs when
δτ
δτ

δJ
=0, this determines the optimal conditions, namely, Eqs. (29) and (30).
δτ

■

Presume that we already have a feasible initial solution x (t ) and u (t ) that satisfy Eqs. (13) and (14), then theoretically the
infinite-dimensional variation dynamic evolution equations (22), (27) and (28) may be used to obtain the optimal solution that
minimizes Eq. (12). The corn during the computation is to solve H o (t , s ) which may be evaluated upon the current solution of
x (t ) and u(t ) .
D. Formulation of Evolution PDE
Similarly, we may use the partial differential operator “ ∂ ” and the differential operator “ d ” to reformulate the variation
dynamic equations to get the following Evolution PDE (EPDE) and Evolution Differential Equation (EDE) as
⎡
∂ ⎡ x⎤ ⎢
=⎢
∂τ ⎣⎢ u ⎦⎥ ⎢
− K Lu + f u Tϕ x +
⎢⎣

{

(∫

tf

t

∂u
(s) d s
∂τ

⎤
⎥
⎥
H o T (σ , t ) Lx (σ ) + ϕtx (σ ) + ϕ xx T (σ ) f (σ ) + f x (σ )T ϕ x (σ ) d σ ⎥
⎦⎥
t

∫t

H o (t , s )

0

(

dt f
dτ

)

(

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

)

)}

(43)

(44)
tf

Put into this perspective, the definite conditions are the initial guess of t f , i.e., t f

τ =0

= tf , and

⎡ x (t ,τ ) ⎤
⎡ x (t ) ⎤
=⎢
⎢ u(t ,τ ) ⎥
⎥
⎣
⎦ τ = 0 ⎣ u (t ) ⎦

(45)

where x (t ) and u (t ) are the initial feasible solution. Recall Fig. 2, Eqs. (43) and (44) realize the anticipated variable evolution
along the variation time τ . The initial conditions of x (t ,τ ) and u(t ,τ ) at τ = 0 belong to the feasible solution domain and their
value at τ = +∞ may be the optimal solution of the OCP. The right part of the EPDE (43) 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.
This paper considers the OCPs defined with free terminal time, yet the results obtained are also applicable to the simpler case
with fixed terminal time. For those OCPs, the evolution equations regarding state variable x and control variable u are same,
while the equation regarding the terminal time t f is not necessary anymore.

V. ILLUSTRATIVE EXAMPLES
First a linear example taken from Xie [19] is considered.

Example 1: Consider the following dynamic system
x = Ax + bu
⎡x ⎤
⎡0 1 ⎤
⎡0⎤
where x = ⎢ 1 ⎥ , A = ⎢
, b = ⎢ ⎥ . Find the solution that minimizes the performance index
⎥
⎣0 0⎦
⎣1 ⎦
⎣ 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 A ( t − s ) ∂u
⎡
⎤
∫t0 e b ∂τ (s) d s
⎥
∂ ⎡ x⎤ ⎢
⎥
=⎢
⎢
⎥
∂τ ⎣ u ⎦ ⎢
⎧
⎞ ⎫⎥
T
T ⎛ tf
T
A (σ −t ) T
−
+
+
+
+
+
K
Ru
e
u
σ
σ
σ
b
Fx
b
(
Q
FA
A
F
)
x
(
)
Fb
(
)
d
⎨
⎜
⎟⎬
⎢⎣
⎝ ∫t
⎠ ⎭⎥⎦
⎩

(

)(

)

with the one-dimensional matrix K = 2 × 10−2 . The definite conditions of the EPDE, i.e., the initial guess of the states x (t ) and the
control u (t ) , were obtained by numerical integration with control input 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. 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 [20], a Radau PS method based OCP solver.
Figs. 3, 4 and 5 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. Fig. 6 plots the profile of
performance index value against the variation time. It declines rapidly at first and almost reaches the minimum when τ = 50s. Then
it keeps approaching the minimum monotonously.
4
3.5

The optimal solution
Numerical solutions with VEM
τ = 0s

3

x1(t)

2.5
τ = 0.4s

2

τ = 1.0s

1.5
τ = 20.7s

1
0.5
0
0

τ = 1.8s
τ = 300s

0.5

1

1.5

2

2.5

3

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

11

The optimal solution
Numerical solutions with VEM

1.2
1

τ = 0s

x2(t)

0.8

τ = 0.4s

0.6
0.4

τ = 1.0s

0.2

τ = 1.8s

0
-0.2

τ = 20.7s

-0.4

τ = 300s

0

0.5

1

1.5

2

2.5

3

t
Fig. 4 The evolution of numerical solutions of

1
0
-1

x2

to the optimal solution.

τ = 0s
τ = 0.4s
τ = 1.8s

u (t)

τ = 20.7s

-2
-3

τ = 70.9s

-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. 5 The evolution of numerical solutions of

u

to the optimal solution.

50
Minimum of J
Funcitonal value of numerical solution
40
40
20

J

J

30

20

0
0

10

0
0

50

100

150

τ (s)

200

2

4 6

τ (s)

250

300

Fig. 6 The approach to the minimum of performance index.

Now we consider a nonlinear example with free terminal time t f , the homing missile problem adapted from Hull [21].

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 ⎥
. x and y are the relative abscissa and ordinate, respectively. θ M is the
where x = ⎢⎢ y ⎥⎥ , and f = ⎢ T
⎢
⎥
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 ⎥⎦
⎢⎣ 0 deg ⎥⎦
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 (43) and the EDE (44), the parameters K and kt f were set to be 1.5 × 10−6 and 1× 10−4 , respectively. The definite
⎡ x (t ,τ ) ⎤
⎢
⎥
conditions, i.e., ⎢u (t ,τ ) ⎥
, were obtained by numerical integration at time horizon [0, 25]s with control input u (t ) = 0 . We also
⎢ t f (τ ) ⎥
⎣
⎦ τ =0
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. 7 gives the states curve in the x y relative coordinate plane, showing that the numerical results approach the optimal
solution over time. For the optimal solution, the missile will intercept the target with a fairly small position error. The control
solutions are plotted in Fig. 8, and the asymptotical approach of the numerical results are demonstrated. In Fig. 9, the terminal time
profile against the variation time τ is plotted. The result of t f oscillates sharply at first and then gradually approaches to the
optimal interception time, and it only changes slightly after τ = 40s. At τ = 300s, we compute that t f = 23.52s from the VEM,
same to the result from GPOPS-II.

13

12000
The optimal solution
Numerical solutions with VEM

10000
τ = 0s

y (m)

8000
τ = 0.6s

6000
4000

Initial relative coordintates

τ = 1.2s

τ = 300s

2000

τ = 2.1s
τ = 3.1s

0
-5000

Optimal terminal relative coordintates

0

5000

10000

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

xy

relative coordinate plane to the optimal solution.

70
The optimal solution
Numerical solutions with VEM

τ = 300s

60
50

u (t)

40

τ = 3.1s
τ = 2.1s

30

τ = 1.2s

20
10

τ = 0.6s

0
-10
0

τ = 0s

5

10

15

20

25

t
Fig. 8. The evolution of numerical solutions of

u

to the optimal solution.

25
Minimum decline time
Evolving profile of tf

24.8
24.6

25

tf (s)

24.2

f

t (s)

24.4

24

24.5
24
23.5
-2 0 2 4 6

23.8

τ (s)

23.6
23.4
0

50

100

150

τ (s)

200

250

300

Fig. 9 The evolution profile of t f to the optimal interception time.

14

VI. FURTHER COMMENTS
Through the substitution of operator, the variation dynamic evolution equations may be reformulated as PDE. Discussions
between the EPDE (43) from the compact VEM and the Augmented Evolution PDE (AEPDE) derived in [10], which was called
the ZS equation originally, is worthwhile. For convenience, the AEPDE is again presented here.
⎛
⎡⎛
∂λ ⎞ ⎤
⎡ ⎛ ∂x
⎞ ⎤⎞
⎜
⎢⎜ H x + ∂t ⎟ ⎥
⎢ ⎜ ∂t − f ⎟ ⎥ ⎟
⎠⎥
⎝
⎠ ⎥⎟
⎜
⎢⎝
∂y
∂ ⎢
= −2 K ⎜ H yy ⎢ ⎛
∂x ⎞ ⎥ − ⎢⎛ ∂λ
⎞⎥ ⎟
⎜
∂τ
⎢ ⎜ f − ⎟ ⎥ ∂t ⎢⎜ + H x ⎟ ⎥ ⎟
∂t ⎠ ⎥
⎠⎥ ⎟
⎜
⎢⎝
⎢⎝ ∂t
⎜
⎢⎣
⎥
⎢
⎥⎦ ⎟⎠
H
0
⎣
u
⎦
⎝

(46)

⎡ x⎤
where y = ⎢⎢ λ ⎥⎥ and H = L + λT f is the Hamiltonian.
⎢⎣ u ⎥⎦

Advantages of the EPDE over the AEPDE are listed as follows. i) The costates are not included and the evaluation on the
second-order derivatives of H is avoided, which will relieve the computation burden. ii) The solution of the EPDE is more
capable to reach a minimum, while the solution of AEPDE may halt at a saddle extremum. iii) Generally the EPDE requires the
integration, and the differentiation, as displayed in the AEPDE, may be avoided. This is advantageous to reduce the numerical error
in seeking solutions.
Regarding the disadvantages, i) currently the EPDE may only address OCPs with free terminal states. When the terminal states
are constrained, which is common in the OCP formulation, it is not directly applicable. One may penalize the terminal boundary
conditions in the performance index. However, this is not a satisfactory solution. ii) Moreover, it requires the initial solution to be
feasible, which is also inflexible.
About the computation, both the right parts of the PDEs are only the functions of time t . This makes them suitable to be solved
with the semi-discrete method in the field of PDE numerical calculation. Then they may be solved with the mature ODE integration
methods. During the numerical calculation, techniques such as the Legendre-Gauss (LG), Legendre-Gauss-Radau (LGR) and
Legendre-Gauss-Lobatto (LGL) discretization, which have stronger approximation capacity, may further improve the accuracy
and efficiency. The common deficiency that both equations face is the incapability to address the general state- and
control-constrained OCPs. These are interesting topics for future studies. Return to the principle, since complex numerical
computations are avoided, and the integration may be achieved with the simple analog circuit, the equations may provide a
promising way of optimal control in the engineering.

VII. CONCLUSION
A compact version of the Variation Evolving Method (VEM) for the optimal control computation is developed. It introduces no
extra variables in transforming the Optimal Control Problems (OCPs) to the Initial-value Problems (IVPs). The optimality
conditions that use no information of the costates are derived and they are proved to be equivalent to the traditional optimality
conditions. One important conclusion drawn by analyzing the optimality conditions is that generally the analytic optimal feedback
control law does not exists because the control is related to the future state. With the mature Ordinary Differential Equation (ODE)
integration methods, the solution of the resulting IVP may reach the minimum of the OCP. In particular, under the frame of IVPs
upon continuous-time dynamics, daunting task of searching reasonable step size and annoying oscillation phenomenon around the
extremum, as occurs in the discrete numerical method, are eliminated. However, currently the method can only solve the OCP with
free terminal states, and this restricts its application. Further studies will be carried out to address this issue.

15
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] 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.
[4] O. V. Stryk and R. Bulirsch, “Direct and indirect methods for trajectory optimization,” Ann. Oper. Res., vol. 37, no. 1, pp. 357-373, 1992.
[5] 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.
[6] 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.
[7] 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.
[8] 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.
[9] 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.
[10] S. Zhang, E. M. Yong, W. Q. Qian, and K. F. He, “A Variation Evolving Method for Optimal Control,” arXiv: 1703.10263 [cs.SY]
[11] 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.
[12] 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.
[13] H. K. Khalil, Nonlinear Systems. New Jersey, USA: Prentice Hall, 2002, pp. 111-181.
[14] D. G. Wu, Variation Method. Beijing, China: Higher Education Press, 1987, pp. 12-68.
[15] K. W. Cassel, Variational Methods with Applications in Science and Engineering. New York: Cambridge University Press, 2013, pp. 28-81
[16] D.Z. Zhen. Linear System Theory. Beijing, China: Tsinghua University Press, 2002, pp. 85-134.
[17] A. E. Bryson and Y. C. Ho, Applied Optimal Control: Optimization, Estimation, and Control. Washington, DC, USA: Hemisphere, 1975, pp.
42-125.
[18] H. Wang, J. S. Luo, and J. M. Zhu. Advanced Mathematics. Changsha, China: National University of Defense Technology Press, 2000, pp.
203-210.
[19] X. S. Xie, Optimal Control Theory and Application. Beijing, China: Tsinghua Univ. Press, 1986, pp. 167-215.
[20] 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.
[21] D. G. Hull, Optimal Control Theory for Applications, New York, USA: Springer, 2003, pp. 89-71.

