# Survey of Direct Transcription for Low-Thrust Space Trajectory Optimization with Applications

Authors: F. Topputo; C. Zhang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Space trajectory design is usually addressed as an optimal control problem. Although it relies on the classic theory of optimal control, this branch possesses some peculiarities that led to the development of ad hoc techniques, which can be grouped into two categories: direct and indirect methods. This paper gives an overview of the principal techniques belonging to the direct methods. The technique known as “direct transcription and collocation” is illustrated by considering Hermite-Simpson, high-order Gauss-Lobatto, and pseudospectral methods. Practical examples are given, and several hints to improve efficiency and robustness are implemented.


---

## Body

## 1. Introduction

In January 1959, the Soviet Union launched Luna 1, a lunar probe that became the first spacecraft placed in heliocentric orbit (http://en.wikipedia.org/wiki/Luna_1). This event marked a new era for deep space exploration. In over half a century, various probes visited the major planets of the Solar System, and they provided a remarkable scientific return. Nowadays, many space agencies have established ambitious space programs, whose accomplishment stimulated the development of new technologies.Electric propulsion is one of those technologies that improve the efficiency of space transport. Electric motors have a specific impulse that is approximately ten times that of the chemical engines, and therefore they allow a considerable saving of propellant mass, which in turn makes it possible to embark heavier instruments, so increasing the scientific return of a mission. Electric propulsion also involves long working hours, extended launch windows, and flexible and precise control abilities. Electric propulsion has been used in NASA’s Deep Space 1 in 1998 [1], in JAXA’s Hayabusa in 2003 [2], and in ESA’s SMART-1 in 2003 [3]. Many future space missions foresee the use of this technology. If compared to that of the chemical engines, electric propulsion produces low levels of thrust, which cannot be modelled as producing instantaneous velocity changes [4]. This is also the case of solar sails, an emerging and prominent technology demonstrated by the recent success of JAXA’s IKAROS [5] and NASA’s NanoSail-D2 [6]. Solar sails are not propellant-constrained, and the effect of the solar radiation pressure on a large surface can be modelled as producing a low-thrust acceleration [7, 8]. The combination of low-thrust propulsion with gravity assists is among the most promising techniques for deep space explorations, which can be utilized to save fuel and shorten the flight time [9–12].Low-thrust space trajectories are studied as a specialization of the optimal control problem for continuous time systems. Unfortunately, no analytic solutions exist for this problem even under the simple two-body dynamics. Thus, numerical methods must be used. Low-thrust trajectory optimization involves determining the control law (thrust magnitude and direction) and the associated transfer orbit while minimizing a given performance index (propellant mass or time-of-flight) and satisfying boundary conditions (departure and arrival orbits), midpoint conditions (at patching points), and path constraints (thrust saturation). In this branch of optimal control, two main solution methods have been devised. Within direct methods, the state and control variables are discretized, and the optimal control problem is converted into a nonlinear programming (NLP) problem [13, 14]. This process is called direct transcription. Direct methods are generally robust and can easily accommodate path constraints, but they often require much computational effort especially for multispiral trajectories. Indirect methods rely on the calculus of variations. The necessary conditions of optimality require the solution of a two-point boundary value problem (TPBVP). This method ensures rapid convergence of good starting guesses, but most of the difficulties are related to the small convergence radius, the high sensitivity to the initial costates, and their lack of physical meaning.This paper complements previous surveys on optimal control of space trajectories [15–18] by focussing on direct methods. In particular, the direct transcription and collocation is studied. The collocation schemes considered implementing different integration methods (Hermite-Simpson, high-order Gauss-Lobatto, and pseudospectral methods). These techniques are implemented to solve practical examples, and the behavior of each method is assessed against the implementation of several issues aimed at improving the computational efficiency and the methods robustness.The remainder of the paper is organized as follows. In Section2, the optimal control problem in space flight mechanics is stated, and the NLP problem as well as the direct transcription concept is recalled. In Section 3 direct transcription and collocation is treated. Computational issues are discussed in Section 4 and practical cases are solved in Section 5. Concluding remarks are given in Section 6.


## 2. The Optimal Control Problem in Space Flight Mechanics

Given a set ofn first-order differential equations
(1)
x
˙
=
f
(
x
,
u
,
t
)
,
where x
(
t
)
∈
R
n is a vector of state variables, u
(
t
)
∈
R
m denotes vector of control variables, and t represents the independent time variable, t
∈
[
t
i
t
f
], the following performance index has to be minimized
(2)
J
=
φ
(
x
(
t
f
)
,
t
f
)
+
∫
t
i
t
f
‍
L
(
x
,
u
,
t
)
d
t
,
while satisfying q-dimensional final boundary conditions
(3)
ψ
(
x
(
t
f
)
,
u
(
t
f
)
,
t
f
)
=
0
.
The solution to this problem is derived by the calculus of variations, which leads to the derivation of the Euler-Lagrange equations. Let ν be a q-dimensional constant vector of multipliers of the final boundary constraints, and let λ be the n-dimensional variable vector of adjoint or costate multipliers of the dynamics. The augmented performance index is defined as
(4)
J
¯
=
φ
(
x
(
t
f
)
,
t
f
)
+
ν
T
ψ
(
x
(
t
f
)
,
u
(
t
f
)
,
t
f
)
+
∫
t
i
t
f
‍
[
L
(
x
,
u
,
t
)
+
λ
T
(
f
(
x
,
u
,
t
)
-
x
˙
)
]
d
t
.
The augmented performance index (4) embeds the dynamics (1) as well as the final boundary conditions (3). The problem consists in deriving the necessary conditions for a stationary point of J
¯ [19]. This is achieved by imposing that its first variation is zero, namely, δ
J
¯
=
0. In order to write the necessary conditions in a compact form, it is convenient to define the Hamiltonian
(5)
H
(
x
,
λ
,
u
,
t
)
=
L
(
x
,
u
,
t
)
+
λ
T
f
(
x
,
u
,
t
)
.
The necessary conditions for optimality, also referred to as the Euler-Lagrange equations [20, 21], are
(6)
x
˙
=
H
λ
,
λ
˙
=
-
H
x
,
0
=
H
u
,
where the subscript denotes partial derivation. The first of (6) is equivalent to (1), the second describes the dynamics of the costates, and the third is an algebraic equation for the control functions. This differential-algebraic system must be solved together with the final boundary conditions (3) and the following transversality conditions:
(7)
λ
(
t
f
)
=
[
φ
x
+
ν
T
ψ
x
]
t
=
t
f
.
With x
(
t
i
)
=
x
i given, the problem represents a two-point boundary value problem. The last of (6) is an application of the Pontryagin maximum principle [22]. A more general expression is in fact
(8)
u
=
arg
min
⁡
u
∈
U
H
(
x
,
λ
,
u
,
t
)
,
where U defines the domain of feasible controls. The maximum principle states that the control variables must be chosen to optimize the Hamiltonian at every instant of time: the solution of the optimal control problem is an extremum for H. In essence, the maximum principle is a constrained optimization problem in the function u
(
t
) at all values of t.


### 2.1. The Optimal Trajectory Design Problem

An optimal trajectory design problem is a specialization of the classical optimal control problem above. In space flight mechanics, the equations of motion have the following form [16]:
(9)
x
˙
=
{
r
˙
v
˙
}
=
{
v
g
(
r
)
+
a
c
u
^
}
.
In (9), r and v are the spacecraft position and velocity vectors, respectively, g
(
r
) is the gravitational vector field, a
c is the thrust acceleration magnitude, and u
^ is the thrust direction unit vector. The control variables are a
c and u
^. The control acceleration magnitude is upper bounded for technological reasons; that is, 0
≤
a
c
≤
a
c
max
⁡. To minimize the total velocity change, and therefore the propellant mass, the objective function (2) is such that L
=
a
c, φ
=
0. The Hamiltonian (5) is, therefore,
(10)
H
=
a
c
+
λ
r
·
v
+
λ
v
·
[
g
(
r
)
+
a
c
u
^
]
=
a
c
[
1
+
λ
v
·
u
^
]
+
λ
r
·
v
+
λ
v
·
g
(
r
)
,
where λ
r and λ
v are the costate vectors associated to the position and velocity vectors, respectively. As a
c is bounded, the last of (6) cannot be applied, and (8) has to be used instead. Since the Hamiltonian has to be minimized at any time, the following observations can be made.(1)
The thrust unit vectoru
^ has to be parallel and opposite to λ
v. That is, u
^
=
-
α
λ
v with α
>
0; this explains why λ
v is also referred to as theprimer vector.
(2)
The thrust magnitudea
c has to be chosen according to the sign of theswitching function:
(11)
S
=
1
+
λ
v
·
u
^
.
In particular, a
c
=
a
c
max
⁡ when S
<
0 and a
c
=
0 when S
>
0. This makes the function a
c
(
t
) have abang-bang structure; that is, it is piece-wise discontinuous and it is either zero or maximum. This property is of great importance to evaluate the optimal control profile a posteriori.The necessary conditions only guarantee that the optimal trajectory is an extremum for the Hamiltonian. Thus, to assess the optimality of the solutions, one is supposed to check the second-order conditions. However, due to the nature of the space trajectory problem, there is no upper bound to the propellant that can be consumed in one trajectory, so one may be confident that a solution that satisfies the necessary conditions is a local minimum and not a local maximum [16].In trajectory optimization problems, the initial state,x
(
t
i
), is generally given. If the initial costate, λ
(
t
i
), was given as well, the optimal solution could be obtained by integrating the Euler-Lagrange equations with an implicit “bang-bang” thrusting structure [23]. Thus, the low-thrust optimal trajectory design can be converted into a TPBVP, which consists in finding the unknown initial costate vector. This is the essence of indirect methods. Another philosophy consists instead in translating the continuous optimal control problem into a NLP problem and solving for a finite set of variables [15, 24–27]. This procedure is the direct transcription and the approach is said direct method.


#### 2.1.1. The General Optimal Trajectory Design Problem

It is convenient to state the optimal trajectory design problem in a more general fashion, which copes with the direct approach. The dynamics are let to incorporate a number of constant parametersp; that is,
(12)
x
˙
=
f
(
x
(
t
)
,
u
(
t
)
,
p
,
t
)
,
and initial and final conditions can be defined within some prescribed lower and upper bounds
(13)
ψ
i
,
l
≤
ψ
i
(
x
(
t
i
)
,
u
(
t
i
)
,
p
,
t
i
)
≤
ψ
i
,
u
,
ψ
f
,
l
≤
ψ
f
(
x
(
t
f
)
,
u
(
t
f
)
,
p
,
t
f
)
≤
ψ
f
,
u
.
In addition, the solution can be subject to path constraints of the form
(14)
g
l
≤
g
(
x
(
t
)
,
u
(
t
)
,
p
,
t
)
≤
g
u
,
as well as simple bounds on the state variables
(15)
x
l
≤
x
(
t
)
≤
x
u
,
and on the control variables
(16)
u
l
≤
u
(
t
)
≤
u
u
.
The basic problem is to determine the control vectors u
(
t
) to minimize the performance index
(17)
J
=
ϕ
(
x
(
t
f
)
,
t
f
)
,
which is written in the Mayer form [15].


### 2.2. The Nonlinear Programming Problem

Essentially, any numerical method for solving the trajectory optimization problem incorporates some type of Newton method to solve for a finite set of unknowns. In Section3 it is shown how an optimal control problem can be transformed into a NLP problem [16, 27]. A NLP problem is a decisional problem concerning a scalar objective function and a vector of constraints. As opposite to the optimal control problem, no dynamics is involved in a NLP problem. Suppose that the n variables x must be chosen to solve
(18)
min
⁡
x
F
(
x
)
,
subject to the m equality constraints
(19)
c
(
x
)
=
0
,
where m
≤
n. The Lagrangian of this problem is
(20)
L
(
x
,
λ
)
=
F
(
x
)
-
λ
T
c
(
x
)
,
which is a scalar function of the n variables x and the m Lagrange multipliers λ. The necessary conditions for a point (x
*
,
λ
*) to be a constrained optimum require solving the following system:
(21)
∇
x
L
(
x
,
λ
)
=
g
(
x
)
-
G
T
(
x
)
λ
=
0
,
∇
λ
L
(
x
,
λ
)
=
-
c
(
x
)
=
0
,
where g
=
∇
x
F and G are the gradient of the objective function and the Jacobian of the equality constraint vector, respectively. The system (21) can be solved via a Newton method to find the (
n
+
m
) variables (
x
*
,
λ
*
). Given a generic initial guess (
x
,
λ
), its corrections (
Δ
x
,
Δ
λ
) to construct the new solution (x
+
Δ
x
,
λ
+
Δ
λ) are given by solving the linear system
(22)
[
H
L
-
G
T
G
0
]
{
Δ
x
Δ
λ
}
=
{
-
g
-
c
}
,
also referred to as Karush-Kuhn-Tucker (KKT) system. In (22), H
L is the Hessian of (20) in x; namely,
(23)
H
L
=
∇
x
2
F
-
∑
i
=
1
m
‍
λ
i
∇
x
2
c
i
.
It is important to observe that an equivalent way to define the search direction Δ
x is to minimize the quadratic form
(24)
1
2
Δ
x
T
H
L
Δ
x
+
g
T
Δ
x
subject to the linear constraints
(25)
G
Δ
x
=
-
c
.
This is the reason why this problem is also referred to as a quadratic programming (QP) problem.The NLP problem formulated above can be generalized to the case that occur when inequality constraints are imposed; them constraints are of the form
(26)
c
(
x
)
≥
0
.
Constraints that are strictly satisfied, that is, those for which c
i
(
x
)
>
0, are called inactive; the remaining active set of constraints are on their bounds; that is, c
j
(
x
)
=
0. If the active set of constraints is known, the inactive constraints are ignored and the problem is simply solved using the method for an equality constrained problem discussed above.In summary, the general NLP problem requires finding then vectors to solve
(27)
min
⁡
x
F
(
x
)
,
subject to the m constraints
(28)
c
L
≤
c
(
x
)
≤
c
U
,
and bounds
(29)
x
L
≤
x
≤
x
U
.
In this formulation equality constraints can be imposed by setting c
j
,
L
=
c
j
,
U.


### 2.3. Direct Transcription

With a direct approach, the solution to the optimal control problem is strictly connected to the numerical integration of the differential equations. The core of this method consists in the way the dynamics are handled; the set of differential equations governing the motion of a spacecraft can be transcribed into a finite set of equality constraints. If these are respected, then the original optimal trajectory design problem is solved within the degree of accuracy of the numerical scheme used [16, 18, 27–30].In problem (12)–(17), the time domain can be uniformly discretized as
(30)
t
i
=
t
1
<
t
2
<
⋯
<
t
N
=
t
f
,
where the time labels are referred to as mesh points or nodes; t
i and t
f are the initial and final time, respectively, and h
=
(
t
N
-
t
1
)
/
(
N
-
1
) is the fixed step size of the discretization (a nonuniform time discretization is also possible). The states and the controls can be discretized over the mesh (30) by defining x
k
=
x
(
t
k
) and u
k
=
u
(
t
k
). The discretized states and the controls are now ready to be treated as a set of NLP variables. The whole variable vector of the problem is
(31)
y
=
{
x
1
,
u
1
,
…
,
x
N
,
u
N
}
T
.
The differential equations are replaced by a finite set of defects constraints derived by the numerical integration scheme. For instance, if a forward Euler scheme is used, the defects are of the form
(32)
ζ
k
≡
x
k
+
1
-
x
k
-
h
f
k
,
where f
k
=
f
(
x
k
,
u
k
,
p
,
t
k
). As a result of the transcription, the optimal control constraints (13)-(14) are replaced by the NLP constraints
(33)
c
L
≤
c
(
y
)
≤
c
U
,
where
(34)
c
(
y
)
≡
{
ζ
1
,
ζ
2
,
…
,
ζ
N
-
1
,
ψ
1
,
ψ
N
,
g
1
,
g
2
,
…
,
g
N
}
T
,
with ψ
1
=
ψ
i, ψ
N
=
ψ
f, and
(35)
c
L
≡
{
0
,
…
,
0
,
g
1
,
L
,
…
,
g
N
,
L
}
T
,
c
U
≡
{
0
,
…
,
0
,
g
1
,
U
,
…
,
g
N
,
U
}
T
.The firstn
(
N
-
1
) equality constraints in (34) require that the defect vectors ζ
k, k
=
1
,
…
,
N
-
1, are zero, thereby they satisfy the differential equations (12) within the accuracy of the numerical integration. The boundary conditions (13) are enforced directly by the equality constraints on ψ
1 and ψ
N, and the nonlinear path constraints (14) are imposed at the grid points. In a similar fashion the objective function, either in the form (2) or (17), can be written in terms of y; namely, F
=
F
(
y
). The optimal trajectory design problem is so translated into the forms (27)–(29) and it can be solved as a standard NLP problem through (18)–(23).In this derivation, the simple forward Euler scheme is used. More accurate methods are likely to be used for practical applications. This is done to keepN within reasonable values, avoiding out-of-memory problems. The integration schemes influence the robustness of the method and the solution accuracy. In Section 3, several implicit quadrature methods are shown.


## 2.1. The Optimal Trajectory Design Problem

An optimal trajectory design problem is a specialization of the classical optimal control problem above. In space flight mechanics, the equations of motion have the following form [16]:
(9)
x
˙
=
{
r
˙
v
˙
}
=
{
v
g
(
r
)
+
a
c
u
^
}
.
In (9), r and v are the spacecraft position and velocity vectors, respectively, g
(
r
) is the gravitational vector field, a
c is the thrust acceleration magnitude, and u
^ is the thrust direction unit vector. The control variables are a
c and u
^. The control acceleration magnitude is upper bounded for technological reasons; that is, 0
≤
a
c
≤
a
c
max
⁡. To minimize the total velocity change, and therefore the propellant mass, the objective function (2) is such that L
=
a
c, φ
=
0. The Hamiltonian (5) is, therefore,
(10)
H
=
a
c
+
λ
r
·
v
+
λ
v
·
[
g
(
r
)
+
a
c
u
^
]
=
a
c
[
1
+
λ
v
·
u
^
]
+
λ
r
·
v
+
λ
v
·
g
(
r
)
,
where λ
r and λ
v are the costate vectors associated to the position and velocity vectors, respectively. As a
c is bounded, the last of (6) cannot be applied, and (8) has to be used instead. Since the Hamiltonian has to be minimized at any time, the following observations can be made.(1)
The thrust unit vectoru
^ has to be parallel and opposite to λ
v. That is, u
^
=
-
α
λ
v with α
>
0; this explains why λ
v is also referred to as theprimer vector.
(2)
The thrust magnitudea
c has to be chosen according to the sign of theswitching function:
(11)
S
=
1
+
λ
v
·
u
^
.
In particular, a
c
=
a
c
max
⁡ when S
<
0 and a
c
=
0 when S
>
0. This makes the function a
c
(
t
) have abang-bang structure; that is, it is piece-wise discontinuous and it is either zero or maximum. This property is of great importance to evaluate the optimal control profile a posteriori.The necessary conditions only guarantee that the optimal trajectory is an extremum for the Hamiltonian. Thus, to assess the optimality of the solutions, one is supposed to check the second-order conditions. However, due to the nature of the space trajectory problem, there is no upper bound to the propellant that can be consumed in one trajectory, so one may be confident that a solution that satisfies the necessary conditions is a local minimum and not a local maximum [16].In trajectory optimization problems, the initial state,x
(
t
i
), is generally given. If the initial costate, λ
(
t
i
), was given as well, the optimal solution could be obtained by integrating the Euler-Lagrange equations with an implicit “bang-bang” thrusting structure [23]. Thus, the low-thrust optimal trajectory design can be converted into a TPBVP, which consists in finding the unknown initial costate vector. This is the essence of indirect methods. Another philosophy consists instead in translating the continuous optimal control problem into a NLP problem and solving for a finite set of variables [15, 24–27]. This procedure is the direct transcription and the approach is said direct method.


### 2.1.1. The General Optimal Trajectory Design Problem

It is convenient to state the optimal trajectory design problem in a more general fashion, which copes with the direct approach. The dynamics are let to incorporate a number of constant parametersp; that is,
(12)
x
˙
=
f
(
x
(
t
)
,
u
(
t
)
,
p
,
t
)
,
and initial and final conditions can be defined within some prescribed lower and upper bounds
(13)
ψ
i
,
l
≤
ψ
i
(
x
(
t
i
)
,
u
(
t
i
)
,
p
,
t
i
)
≤
ψ
i
,
u
,
ψ
f
,
l
≤
ψ
f
(
x
(
t
f
)
,
u
(
t
f
)
,
p
,
t
f
)
≤
ψ
f
,
u
.
In addition, the solution can be subject to path constraints of the form
(14)
g
l
≤
g
(
x
(
t
)
,
u
(
t
)
,
p
,
t
)
≤
g
u
,
as well as simple bounds on the state variables
(15)
x
l
≤
x
(
t
)
≤
x
u
,
and on the control variables
(16)
u
l
≤
u
(
t
)
≤
u
u
.
The basic problem is to determine the control vectors u
(
t
) to minimize the performance index
(17)
J
=
ϕ
(
x
(
t
f
)
,
t
f
)
,
which is written in the Mayer form [15].


## 2.1.1. The General Optimal Trajectory Design Problem

It is convenient to state the optimal trajectory design problem in a more general fashion, which copes with the direct approach. The dynamics are let to incorporate a number of constant parametersp; that is,
(12)
x
˙
=
f
(
x
(
t
)
,
u
(
t
)
,
p
,
t
)
,
and initial and final conditions can be defined within some prescribed lower and upper bounds
(13)
ψ
i
,
l
≤
ψ
i
(
x
(
t
i
)
,
u
(
t
i
)
,
p
,
t
i
)
≤
ψ
i
,
u
,
ψ
f
,
l
≤
ψ
f
(
x
(
t
f
)
,
u
(
t
f
)
,
p
,
t
f
)
≤
ψ
f
,
u
.
In addition, the solution can be subject to path constraints of the form
(14)
g
l
≤
g
(
x
(
t
)
,
u
(
t
)
,
p
,
t
)
≤
g
u
,
as well as simple bounds on the state variables
(15)
x
l
≤
x
(
t
)
≤
x
u
,
and on the control variables
(16)
u
l
≤
u
(
t
)
≤
u
u
.
The basic problem is to determine the control vectors u
(
t
) to minimize the performance index
(17)
J
=
ϕ
(
x
(
t
f
)
,
t
f
)
,
which is written in the Mayer form [15].


## 2.2. The Nonlinear Programming Problem

Essentially, any numerical method for solving the trajectory optimization problem incorporates some type of Newton method to solve for a finite set of unknowns. In Section3 it is shown how an optimal control problem can be transformed into a NLP problem [16, 27]. A NLP problem is a decisional problem concerning a scalar objective function and a vector of constraints. As opposite to the optimal control problem, no dynamics is involved in a NLP problem. Suppose that the n variables x must be chosen to solve
(18)
min
⁡
x
F
(
x
)
,
subject to the m equality constraints
(19)
c
(
x
)
=
0
,
where m
≤
n. The Lagrangian of this problem is
(20)
L
(
x
,
λ
)
=
F
(
x
)
-
λ
T
c
(
x
)
,
which is a scalar function of the n variables x and the m Lagrange multipliers λ. The necessary conditions for a point (x
*
,
λ
*) to be a constrained optimum require solving the following system:
(21)
∇
x
L
(
x
,
λ
)
=
g
(
x
)
-
G
T
(
x
)
λ
=
0
,
∇
λ
L
(
x
,
λ
)
=
-
c
(
x
)
=
0
,
where g
=
∇
x
F and G are the gradient of the objective function and the Jacobian of the equality constraint vector, respectively. The system (21) can be solved via a Newton method to find the (
n
+
m
) variables (
x
*
,
λ
*
). Given a generic initial guess (
x
,
λ
), its corrections (
Δ
x
,
Δ
λ
) to construct the new solution (x
+
Δ
x
,
λ
+
Δ
λ) are given by solving the linear system
(22)
[
H
L
-
G
T
G
0
]
{
Δ
x
Δ
λ
}
=
{
-
g
-
c
}
,
also referred to as Karush-Kuhn-Tucker (KKT) system. In (22), H
L is the Hessian of (20) in x; namely,
(23)
H
L
=
∇
x
2
F
-
∑
i
=
1
m
‍
λ
i
∇
x
2
c
i
.
It is important to observe that an equivalent way to define the search direction Δ
x is to minimize the quadratic form
(24)
1
2
Δ
x
T
H
L
Δ
x
+
g
T
Δ
x
subject to the linear constraints
(25)
G
Δ
x
=
-
c
.
This is the reason why this problem is also referred to as a quadratic programming (QP) problem.The NLP problem formulated above can be generalized to the case that occur when inequality constraints are imposed; them constraints are of the form
(26)
c
(
x
)
≥
0
.
Constraints that are strictly satisfied, that is, those for which c
i
(
x
)
>
0, are called inactive; the remaining active set of constraints are on their bounds; that is, c
j
(
x
)
=
0. If the active set of constraints is known, the inactive constraints are ignored and the problem is simply solved using the method for an equality constrained problem discussed above.In summary, the general NLP problem requires finding then vectors to solve
(27)
min
⁡
x
F
(
x
)
,
subject to the m constraints
(28)
c
L
≤
c
(
x
)
≤
c
U
,
and bounds
(29)
x
L
≤
x
≤
x
U
.
In this formulation equality constraints can be imposed by setting c
j
,
L
=
c
j
,
U.


## 2.3. Direct Transcription

With a direct approach, the solution to the optimal control problem is strictly connected to the numerical integration of the differential equations. The core of this method consists in the way the dynamics are handled; the set of differential equations governing the motion of a spacecraft can be transcribed into a finite set of equality constraints. If these are respected, then the original optimal trajectory design problem is solved within the degree of accuracy of the numerical scheme used [16, 18, 27–30].In problem (12)–(17), the time domain can be uniformly discretized as
(30)
t
i
=
t
1
<
t
2
<
⋯
<
t
N
=
t
f
,
where the time labels are referred to as mesh points or nodes; t
i and t
f are the initial and final time, respectively, and h
=
(
t
N
-
t
1
)
/
(
N
-
1
) is the fixed step size of the discretization (a nonuniform time discretization is also possible). The states and the controls can be discretized over the mesh (30) by defining x
k
=
x
(
t
k
) and u
k
=
u
(
t
k
). The discretized states and the controls are now ready to be treated as a set of NLP variables. The whole variable vector of the problem is
(31)
y
=
{
x
1
,
u
1
,
…
,
x
N
,
u
N
}
T
.
The differential equations are replaced by a finite set of defects constraints derived by the numerical integration scheme. For instance, if a forward Euler scheme is used, the defects are of the form
(32)
ζ
k
≡
x
k
+
1
-
x
k
-
h
f
k
,
where f
k
=
f
(
x
k
,
u
k
,
p
,
t
k
). As a result of the transcription, the optimal control constraints (13)-(14) are replaced by the NLP constraints
(33)
c
L
≤
c
(
y
)
≤
c
U
,
where
(34)
c
(
y
)
≡
{
ζ
1
,
ζ
2
,
…
,
ζ
N
-
1
,
ψ
1
,
ψ
N
,
g
1
,
g
2
,
…
,
g
N
}
T
,
with ψ
1
=
ψ
i, ψ
N
=
ψ
f, and
(35)
c
L
≡
{
0
,
…
,
0
,
g
1
,
L
,
…
,
g
N
,
L
}
T
,
c
U
≡
{
0
,
…
,
0
,
g
1
,
U
,
…
,
g
N
,
U
}
T
.The firstn
(
N
-
1
) equality constraints in (34) require that the defect vectors ζ
k, k
=
1
,
…
,
N
-
1, are zero, thereby they satisfy the differential equations (12) within the accuracy of the numerical integration. The boundary conditions (13) are enforced directly by the equality constraints on ψ
1 and ψ
N, and the nonlinear path constraints (14) are imposed at the grid points. In a similar fashion the objective function, either in the form (2) or (17), can be written in terms of y; namely, F
=
F
(
y
). The optimal trajectory design problem is so translated into the forms (27)–(29) and it can be solved as a standard NLP problem through (18)–(23).In this derivation, the simple forward Euler scheme is used. More accurate methods are likely to be used for practical applications. This is done to keepN within reasonable values, avoiding out-of-memory problems. The integration schemes influence the robustness of the method and the solution accuracy. In Section 3, several implicit quadrature methods are shown.


## 3. Direct Transcription and Collocation

Collocation is used to transcribe differential dynamic constraints into a set of algebraic constraints. The basic idea is to choose a polynomial up to a certain degree with a number of points in the time domain (collocation points), and to enforce the polynomials to satisfy the equations of motion at the collocation points. The peculiarity of each collocation method relies on the way the state and control variables are discretized and how the dynamic constraints are satisfied [17].


### 3.1. Hermite-Simpson Method

A basic form of collocation is the Hermite-Simpson method [24], illustrated in Figure 1. For each of the segments [
t
k
,
t
k
+
1
], the two end points, denoted as “nodes” (blue dots), represent the corresponding state and control NLP variables; that is, [
x
k
,
u
k
,
x
k
+
1
,
u
k
+
1
] (scalar states and controls are used from now on to ease the notation). The dynamics are used to provide time derivative values at the two nodes, so the four pieces of information [
x
k
,
x
k
+
1
,
f
(
x
k
,
u
k
)
,
f
(
x
k
+
1
,
u
k
+
1
)
] can be used to construct a third-order Hermite interpolate polynomial. This interpolate polynomial cannot satisfy the equations of motion at any time within [
t
k
,
t
k
+
1
], because it does that only at the nodes. Let [
x
c
,
u
c
] be the state and control at t
c, the middle point of [
t
k
,
t
k
+
1
]; this is called “collocation point” (red diamond). Enforcing Δ
=
x
˙
c
-
f
(
x
c
,
u
c
)
=
0 makes it possible to have a polynomial that not only satisfies the dynamics at the two nodes but also does that at the collocation point. If a large number of intervals are used, the state motion approaches the real dynamics within the whole time domain.Figure 1
Hermite-Simpson collocation method.The detailed procedure can be derived as follows [16]. Let the state function x
(
t
) be represented on each segment [
t
k
,
t
k
+
1
] through a cubic polynomial of the form
(36)
x
(
t
)
=
a
0
+
a
1
t
+
a
2
t
2
+
a
3
t
3
,
which yields
(37)
x
˙
(
t
)
=
a
1
+
2
a
2
t
+
3
a
3
t
2
,
where [
a
0
,
a
1
,
a
2
,
a
3
] are the coefficients of the polynomial. In order to simplify the argument, the time domain is transformed such that t
∈
[
0
,
h
] (h is the time interval of the segment). Let x
(
0
)
=
x
k, x
(
h
)
=
x
k
+
1, x
˙
(
0
)
=
x
˙
k, and x
˙
(
h
)
=
x
˙
k
+
1. Evaluating (36)-(37) at t
=
0 and t
=
h yields
(38)
[
x
(
0
)
x
˙
(
0
)
x
(
h
)
x
˙
(
h
)
]
=
[
1
0
0
0
0
1
0
0
1
h
h
2
h
3
0
1
2
h
3
h
2
]
[
a
0
a
1
a
2
a
3
]
,
which allows us to compute the four coefficients of interpolate polynomial
(39)
[
a
0
a
1
a
2
a
3
]
=
[
1
0
0
0
0
1
0
0
-
3
h
2
-
2
h
3
h
2
-
1
h
2
h
3
1
h
2
-
2
h
3
1
h
2
]
[
x
(
0
)
x
˙
(
0
)
x
(
h
)
x
˙
(
h
)
]
.
Substituting [
a
0
,
a
1
,
a
2
,
a
3
] into (36)-(37) allows us to compute the collocation point as
(40)
x
c
=
x
(
h
2
)
=
1
2
(
x
k
+
x
k
+
1
)
+
h
8
[
f
(
x
k
,
u
k
)
-
f
(
x
k
+
1
,
u
k
+
1
)
]
,
as well as its time derivative
(41)
x
˙
c
=
x
˙
(
h
2
)
=
-
3
2
h
(
x
k
-
x
k
+
1
)
-
1
4
[
f
(
x
k
,
u
k
)
+
f
(
x
k
+
1
,
u
k
+
1
)
]
.
The control variable at the collocation point can be computed by simple linear interpolation; that is,
(42)
u
c
=
u
k
+
u
k
+
1
2
.
The difference between the interpolated and calculated derivatives at the collocation point defines the integration defect
(43)
Δ
=
x
˙
c
-
f
(
x
c
,
u
c
)
=
-
3
2
h
(
x
k
-
x
k
+
1
)
-
1
4
[
f
(
x
k
,
u
k
)
+
f
(
x
k
+
1
,
u
k
+
1
)
]
-
f
(
x
c
,
u
c
)
=
x
k
-
x
k
+
1
+
h
6
[
f
(
x
k
,
u
k
)
+
4
f
(
x
c
,
u
c
)
+
f
(
x
k
+
1
,
u
k
+
1
)
]
.
The NLP solver will select [
x
k
,
u
k
,
x
k
+
1
,
u
k
+
1
] to drive Δ to zero and in this way the interpolating polynomial will approximate the true dynamics within the accuracy of the numerical integration. Note that the last row of (43) is actually an implicit Hermite integration. Thus, if the collocation constraints are satisfied, the system is said to be “implicitly” integrated.


### 3.2. High Order Gauss-Lobatto Method

In the Hermite-Simpson method, only two nodes are used to construct a third-order polynomial. As a general rule of thumb, a lower number of segments can be handled if higher order integration is performed. Herman and Conway [28] demonstrated that higher order Gauss-Lobatto methods are more robust and more efficient than the lower-order Hermite-Simpson scheme. Figure 2 illustrates the collocation constraints of fifth-order Gauss-Lobatto method. Similar to Hermite-Simpson method, in each segment, six pieces of information of nodes are used to construct a fifth-order Hermite polynomial to approximate the state time history. Then the resulting interpolation polynomial is used to evaluate the states at the remaining two collocation points.Figure 2
Fifth-order Gauss-Lobatto constraint formulation.In the fifth-order method, the collocation points are [28]
(44)
x
1
=
1
686
{
(
39
21
+
231
)
x
k
+
224
x
c
iiiiiiiiiiiii
+
(
-
39
21
+
231
)
x
k
+
1
iiiiiiiiiiiii
+
h
[
(
3
21
+
21
)
f
k
-
16
21
f
c
iiiiiiiiiiiiiiiiiiiii
+
(
3
21
-
21
)
f
k
+
1
]
}
,
x
2
=
1
686
{
(
-
39
21
+
231
)
x
k
+
224
x
c
iiiiiiiiiiiii
+
(
39
21
+
231
)
x
k
+
1
iiiiiiiiiiiii
+
h
[
(
-
3
21
+
21
)
f
k
+
16
21
f
c
iiiiiiiiiiiiiiiiiiiii
+
(
-
3
21
-
21
)
f
k
+
1
]
}
,
and the two collocation constraints defects to zero are
(45)
Δ
1
=
1
360
{
(
32
21
+
180
)
x
k
-
64
21
x
c
iiiiiiiiiiiii
+
(
32
21
-
180
)
x
k
+
1
iiiiiiiiiiiii
+
h
[
(
9
+
21
)
f
k
+
98
f
1
+
64
f
c
iiiiiiiiiiiiiiiiiiiii
+
(
9
-
21
)
f
k
+
1
]
}
,
Δ
2
=
1
360
{
(
-
32
21
+
180
)
x
k
+
64
21
x
c
iiiiiiiiiiiii
+
(
-
32
21
-
180
)
x
k
+
1
iiiiiiiiiiiii
+
h
[
(
9
-
21
)
f
k
+
98
f
2
+
64
f
c
iiiiiiiiiiiiiiiiiiiii
+
(
9
+
21
)
f
k
+
1
]
}
.A detailed derivation of the formulas above is required when the argument is extended to arbitrary higher orders. In [31], an alternative framework for unifying arbitrary higher-order methods is presented. In this approach, the Legendre-Gauss-Lobatto (LGL) discrete points ξ
j are used to improve both interpolating precision and quadrature performance. Figure 3 graphically illustrates the location of three, five, and seven LGL points. Table 1 lists the corresponding position of these LGL points (the time interval is [
-
1,1
]).Table 1
Locations three, five, and seven Legendre-Gauss-Lobatto points.
Number of points
Location
3
0
±
1
5
0
±
3
7
±
1
7
0
±
495
-
66
15
33
±
495
+
66
15
33
±
1Figure 3
Three, five, and seven Legendre-Gauss-Lobatto points.With reference to Figure3, the blue points are the generic nodes, which are defined by
(46)
τ
j
=
ξ
2
j
-
1
,
j
=
1
,
…
,
(
n
+
1
)
2
,
while the red diamonds are the generic collocation points, denoted as
(47)
ζ
j
=
ξ
2
j
,
j
=
1
,
…
,
(
n
-
1
)
2
,
such that there is one collocation point between every two adjacent nodes [31]. The nodes are used for constructing the interpolation polynomial, while the collocation points are used to formulate the defect constraints (it is noted that nodes and collocation points should not overlap in Gauss-Lobatto method). With three points, the Gauss-Lobatto method degenerates to the Hermite-Simpson method. The state in the ith subinterval is approximated by the nth degree Hermite interpolating polynomial
(48)
x
(
τ
)
≃
a
0
+
a
1
τ
+
a
2
τ
2
+
a
3
τ
3
+
⋯
+
a
n
τ
n
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
τ
∈
[
-
1,1
]
,
where the coefficients of Hermite interpolating polynomial a
=
[
a
0
,
a
1
,
a
2
,
…
,
a
n
]
T are determined by using the values of the states and vector field at the points τ
j; that is,
(49)
[
1
τ
1
1
τ
1
2
⋯
τ
1
n
-
1
τ
1
n
1
τ
2
1
τ
2
2
⋯
τ
2
n
-
1
τ
2
n
⋮
⋮
⋮
⋮
⋮
⋮
1
τ
(
n
+
1
)
/
2
1
τ
(
n
+
1
)
/
2
2
⋯
τ
(
n
+
1
)
/
2
n
-
1
τ
(
n
+
1
)
/
2
n
0
1
2
τ
1
1
⋯
(
n
-
1
)
τ
1
n
-
2
n
τ
1
n
-
1
⋮
⋮
⋮
⋮
⋮
⋮
0
1
2
τ
(
n
+
1
)
/
2
1
⋯
(
n
-
1
)
τ
(
n
+
1
)
/
2
n
-
2
n
τ
(
n
+
1
)
/
2
n
-
1
]
×
[
a
0
a
1
a
2
⋮
a
n
-
1
a
n
]
=
[
x
(
τ
1
)
x
(
τ
2
)
⋮
x
(
τ
(
n
+
1
)
/
2
)
h
2
f
(
τ
1
)
⋮
h
2
f
(
τ
(
n
+
1
)
/
2
)
]
.
Let (49) be written as A
a
=
b. The term A depends only on the number and location of node points, whereas b is a concatenated vector of state and vector field values. Thus, given a vector of states, the coefficients of the Hermite interpolation in the ith subinterval, a, can be found by a
=
[
A
]
-
1
b. The values of the states located at collocation points are then
(50)
x
(
ζ
j
)
=
[
1
ζ
1
ζ
1
2
⋯
ζ
1
n
1
ζ
2
ζ
2
2
⋯
ζ
2
n
⋮
⋮
⋮
⋯
⋮
1
ζ
(
n
-
1
)
/
2
ζ
(
n
-
1
)
/
2
2
⋯
ζ
(
n
-
1
)
/
2
n
]
[
A
]
-
1
b
=
Φ
b
,
j
=
1,2
,
…
,
(
n
-
1
)
2
.
The derivatives of the Hermite interpolating polynomial located at the collocation points are given by
(51)
x
˙
(
ζ
j
)
=
[
0
1
2
ζ
1
⋯
n
ζ
1
n
-
1
0
1
2
ζ
2
⋯
n
ζ
2
n
-
1
⋮
⋮
⋮
⋯
⋮
0
1
2
ζ
(
n
-
1
)
/
2
⋯
n
ζ
(
n
-
1
)
/
2
n
-
1
]
[
A
]
-
1
b
=
Φ
′
b
,
j
=
1,2
,
…
,
(
n
-
1
)
2
.
In this form the matrices Φ and Φ
′ are constants, thus a system of (
n
-
1
)
/
2 constraints per interval is obtained as follows:
(52)
Δ
i
=
x
˙
(
ζ
j
)
-
h
2
f
(
x
(
ζ
j
)
,
u
(
ζ
j
)
)
=
Φ
′
b
-
h
2
f
(
Φ
b
,
u
(
ζ
j
)
)
,
and the analytic Jacobian for the defect constraints can be derived by chain rule as
(53)
∂
Δ
i
∂
[
x
(
τ
j
)
,
u
(
τ
j
)
]
i
=
Φ
′
∂
b
∂
[
x
(
τ
j
)
,
u
(
τ
j
)
]
i
-
h
2
[
∂
f
(
x
(
ζ
j
)
,
u
(
ζ
j
)
)
∂
x
(
ζ
j
)
Φ
∂
b
∂
[
x
(
τ
j
)
,
u
(
τ
j
)
]
i
iiiiiiiiiiiiiii
+
∂
f
(
x
(
ζ
j
)
,
u
(
ζ
j
)
)
∂
u
(
ζ
j
)
∂
u
(
ζ
j
)
∂
[
x
(
τ
j
)
,
u
(
τ
j
)
]
i
]
.
If the final transfer time t
f is a variable, within each segment, the analytic Jacobian of t
f is written as
(54)
∂
Δ
i
∂
t
f
=
Φ
′
∂
b
∂
h
∂
h
∂
t
f
-
[
1
2
∂
h
∂
t
f
f
(
x
(
ζ
j
)
,
u
(
ζ
j
)
)
∂
f
(
x
(
ζ
j
)
,
u
(
ζ
j
)
)
∂
u
(
ζ
j
)
iiiiii
+
h
2
∂
f
(
x
(
ζ
j
)
,
u
(
ζ
j
)
)
∂
x
(
ζ
j
)
Φ
∂
b
∂
h
∂
h
∂
t
f
]
.
A graphical illustration of the Jacobian structure with free final time t
f is given in Figure 4. It can be seen that each segment has the same Jacobian module; thus (53)-(54) can be calculated offline and stored before solving the NLP. In this way, the computational time can be considerably reduced compared to the use of finite difference approximation.Figure 4
Jacobian structure for third-order Gauss-Lobatto method.


### 3.3. Pseudospectral Method

Recently, new direct methods have been applied to low-thrust trajectory optimization, the most popular one being the pseudospectral scheme [32–35]. The major difference between Gauss-Lobatto and pseudospectral collocation schemes is the way in which the interpolation polynomial is constructed and the defect constraints are defined. When the pseudospectral method is employed, a global Lagrange interpolation polynomial is constructed to approximate the state profile. The polynomial is then differentiated and evaluated at all the nodes to compute interpolated derivative values of the states, and equations of motion are used to provide physical time derivatives. The differences between the two sets of time derivatives form the defects. The formulation of defect constraints is sketched in Figure 5.Figure 5
Formulation of the collocation constraints in the pseudospectral method.The defect constraints of the pseudospectral method can be derived as follows. The state curvep
(
t
) is approximated by using values of x at the discrete time points t
i and the corresponding Lagrange interpolating polynomials L
i
(
t
), (
i
=
1
,
…
,
n
)
(55)
p
(
t
)
=
∑
i
=
1
n
‍
L
i
(
t
)
x
(
t
i
)
,
where L
i
(
t
)
(
i
=
1
,
…
,
n
) are defined as
(56)
L
i
(
t
)
=
∏
j
=
1
,
j
≠
i
n
‍
t
-
t
j
t
i
-
t
j
.
Then, the derivative of p
(
t
) is given:
(57)
p
′
(
t
)
=
∑
i
=
1
n
‍
L
i
′
(
t
)
x
(
t
i
)
.
The time derivative of the polynomial can be expressed in compact form through
(58)
p
′
=
D
x
,
where D is differentiation matrix. This only depends on the chosen node spacing, so it can be calculated offline and stored. Then the defect constraints can be written as
(59)
Δ
=
p
′
-
f
(
x
,
u
)
=
D
x
-
f
(
x
,
u
)
.
A variety of pseudospectral schemes emerged and the major differences are the way in which the nodes are selected. The most general form is Legendre-Gauss-Lobatto discretization [36], but Legendre-Gauss [34, 37] and Legendre-Gauss-Radau [35] are also commonly used. Figure 6 illustrates three different kinds of Gauss points. Table 2 lists the position of equidistant points with three kinds of Gaussian points. We can see that LG points are located on the open interval (
-
1,1
), in contrast, LGR points are located on half-open interval [
-
1,1
), and LGL points are located on closed interval [
-
1,1
].Table 2
Position of equidistant, LG, LGR, and LGL with five points within[−1, 1].
Case
Type
x
1
x
2
x
3
x
4
x
5
1
Equidistant
−1
−0.5
0
0.5
1
2
LG points
−0.906
−0.538
0
0.538
0.906
3
LGR points
−1
−0.720
−0.167
0.446
0.886
4
LGL points
−1
−0.654
0
0.655
1Figure 6
Legendre-Gauss (LG), Legendre-Gauss-Radau (LGR), and Legendre-Gauss-Lobatto (LGL) with five points.The main characteristic of the Gauss-Lobatto method is the combination of reasonable accuracy with highly sparse constraint Jacobians and Hessians. The pseudospectral method has a more elegant form of Jacobian and it offers spectral accuracy for smooth problems, but the constraint Jacobian is much denser. So there is a balance between accuracy and efficiency. In general, the Gauss-Lobatto method use a limited number of nodes per segment (n
<
4) with many segments, while pseudospectral methods use many nodes per segment (n
>
10) with a limited number of segments.


### 3.4. Control Parameterization Method

The control parameterization method is a further direct transcription technique for solving optimal control problems [38–40]. In this method, only the controls are discretized and represented by a linear combination of basis functions. This approximation scheme yields a suboptimal control for the original problem. The optimal control problem is converted to an approximate nonlinear optimization problem with a finite number of decision variables, which can be solved by nonlinear programming techniques. The advantage of piecewise-constant approximation scheme is due to its simplicity and convergence [41].To solve optimal control problem using the control parameterization method, we approximate the control profileu as (see Figure 7)
(60)
u
(
t
)
≈
u
p
(
t
)
=
σ
k
,
t
∈
[
τ
k
-
1
,
τ
k
)
,
k
=
1
,
…
,
p
,
where p
≥
1 is a given integer, τ
k, k
=
0
,
…
,
p are knot points, and σ
k
∈
R
r, k
=
1
,
…
,
p are vectors containing the approximate control values. The knot points satisfy
(61)
0
=
τ
0
≤
τ
1
≤
τ
2
≤
⋯
≤
τ
p
-
1
≤
τ
p
=
T
.
The approximate control u
p can be written as
(62)
u
p
(
t
)
=
∑
k
=
1
p
-
1
‍
σ
k
χ
[
τ
k
-
1
,
τ
k
)
(
t
)
+
σ
p
χ
[
τ
p
-
1
,
τ
p
)
(
t
)
,
where for a given subinterval I
⊂
[
0
,
T
], the characteristic function χ
I
:
R
→
R is defined by
(63)
χ
I
∶
=
{
1
,
if
t
∈
I
,
0
,
otherwise
.
Note that u
p is a piecewise-constant function with potential discontinuities at the points t
=
τ
k, k
=
1
,
…
,
p
-
1. These points are calledswitching times. Reference [40] surveys the key developments in the control parameterization.Figure 7
Piecewise-constant control approximation with equidistant knot points.


## 3.1. Hermite-Simpson Method

A basic form of collocation is the Hermite-Simpson method [24], illustrated in Figure 1. For each of the segments [
t
k
,
t
k
+
1
], the two end points, denoted as “nodes” (blue dots), represent the corresponding state and control NLP variables; that is, [
x
k
,
u
k
,
x
k
+
1
,
u
k
+
1
] (scalar states and controls are used from now on to ease the notation). The dynamics are used to provide time derivative values at the two nodes, so the four pieces of information [
x
k
,
x
k
+
1
,
f
(
x
k
,
u
k
)
,
f
(
x
k
+
1
,
u
k
+
1
)
] can be used to construct a third-order Hermite interpolate polynomial. This interpolate polynomial cannot satisfy the equations of motion at any time within [
t
k
,
t
k
+
1
], because it does that only at the nodes. Let [
x
c
,
u
c
] be the state and control at t
c, the middle point of [
t
k
,
t
k
+
1
]; this is called “collocation point” (red diamond). Enforcing Δ
=
x
˙
c
-
f
(
x
c
,
u
c
)
=
0 makes it possible to have a polynomial that not only satisfies the dynamics at the two nodes but also does that at the collocation point. If a large number of intervals are used, the state motion approaches the real dynamics within the whole time domain.Figure 1
Hermite-Simpson collocation method.The detailed procedure can be derived as follows [16]. Let the state function x
(
t
) be represented on each segment [
t
k
,
t
k
+
1
] through a cubic polynomial of the form
(36)
x
(
t
)
=
a
0
+
a
1
t
+
a
2
t
2
+
a
3
t
3
,
which yields
(37)
x
˙
(
t
)
=
a
1
+
2
a
2
t
+
3
a
3
t
2
,
where [
a
0
,
a
1
,
a
2
,
a
3
] are the coefficients of the polynomial. In order to simplify the argument, the time domain is transformed such that t
∈
[
0
,
h
] (h is the time interval of the segment). Let x
(
0
)
=
x
k, x
(
h
)
=
x
k
+
1, x
˙
(
0
)
=
x
˙
k, and x
˙
(
h
)
=
x
˙
k
+
1. Evaluating (36)-(37) at t
=
0 and t
=
h yields
(38)
[
x
(
0
)
x
˙
(
0
)
x
(
h
)
x
˙
(
h
)
]
=
[
1
0
0
0
0
1
0
0
1
h
h
2
h
3
0
1
2
h
3
h
2
]
[
a
0
a
1
a
2
a
3
]
,
which allows us to compute the four coefficients of interpolate polynomial
(39)
[
a
0
a
1
a
2
a
3
]
=
[
1
0
0
0
0
1
0
0
-
3
h
2
-
2
h
3
h
2
-
1
h
2
h
3
1
h
2
-
2
h
3
1
h
2
]
[
x
(
0
)
x
˙
(
0
)
x
(
h
)
x
˙
(
h
)
]
.
Substituting [
a
0
,
a
1
,
a
2
,
a
3
] into (36)-(37) allows us to compute the collocation point as
(40)
x
c
=
x
(
h
2
)
=
1
2
(
x
k
+
x
k
+
1
)
+
h
8
[
f
(
x
k
,
u
k
)
-
f
(
x
k
+
1
,
u
k
+
1
)
]
,
as well as its time derivative
(41)
x
˙
c
=
x
˙
(
h
2
)
=
-
3
2
h
(
x
k
-
x
k
+
1
)
-
1
4
[
f
(
x
k
,
u
k
)
+
f
(
x
k
+
1
,
u
k
+
1
)
]
.
The control variable at the collocation point can be computed by simple linear interpolation; that is,
(42)
u
c
=
u
k
+
u
k
+
1
2
.
The difference between the interpolated and calculated derivatives at the collocation point defines the integration defect
(43)
Δ
=
x
˙
c
-
f
(
x
c
,
u
c
)
=
-
3
2
h
(
x
k
-
x
k
+
1
)
-
1
4
[
f
(
x
k
,
u
k
)
+
f
(
x
k
+
1
,
u
k
+
1
)
]
-
f
(
x
c
,
u
c
)
=
x
k
-
x
k
+
1
+
h
6
[
f
(
x
k
,
u
k
)
+
4
f
(
x
c
,
u
c
)
+
f
(
x
k
+
1
,
u
k
+
1
)
]
.
The NLP solver will select [
x
k
,
u
k
,
x
k
+
1
,
u
k
+
1
] to drive Δ to zero and in this way the interpolating polynomial will approximate the true dynamics within the accuracy of the numerical integration. Note that the last row of (43) is actually an implicit Hermite integration. Thus, if the collocation constraints are satisfied, the system is said to be “implicitly” integrated.


## 3.2. High Order Gauss-Lobatto Method

In the Hermite-Simpson method, only two nodes are used to construct a third-order polynomial. As a general rule of thumb, a lower number of segments can be handled if higher order integration is performed. Herman and Conway [28] demonstrated that higher order Gauss-Lobatto methods are more robust and more efficient than the lower-order Hermite-Simpson scheme. Figure 2 illustrates the collocation constraints of fifth-order Gauss-Lobatto method. Similar to Hermite-Simpson method, in each segment, six pieces of information of nodes are used to construct a fifth-order Hermite polynomial to approximate the state time history. Then the resulting interpolation polynomial is used to evaluate the states at the remaining two collocation points.Figure 2
Fifth-order Gauss-Lobatto constraint formulation.In the fifth-order method, the collocation points are [28]
(44)
x
1
=
1
686
{
(
39
21
+
231
)
x
k
+
224
x
c
iiiiiiiiiiiii
+
(
-
39
21
+
231
)
x
k
+
1
iiiiiiiiiiiii
+
h
[
(
3
21
+
21
)
f
k
-
16
21
f
c
iiiiiiiiiiiiiiiiiiiii
+
(
3
21
-
21
)
f
k
+
1
]
}
,
x
2
=
1
686
{
(
-
39
21
+
231
)
x
k
+
224
x
c
iiiiiiiiiiiii
+
(
39
21
+
231
)
x
k
+
1
iiiiiiiiiiiii
+
h
[
(
-
3
21
+
21
)
f
k
+
16
21
f
c
iiiiiiiiiiiiiiiiiiiii
+
(
-
3
21
-
21
)
f
k
+
1
]
}
,
and the two collocation constraints defects to zero are
(45)
Δ
1
=
1
360
{
(
32
21
+
180
)
x
k
-
64
21
x
c
iiiiiiiiiiiii
+
(
32
21
-
180
)
x
k
+
1
iiiiiiiiiiiii
+
h
[
(
9
+
21
)
f
k
+
98
f
1
+
64
f
c
iiiiiiiiiiiiiiiiiiiii
+
(
9
-
21
)
f
k
+
1
]
}
,
Δ
2
=
1
360
{
(
-
32
21
+
180
)
x
k
+
64
21
x
c
iiiiiiiiiiiii
+
(
-
32
21
-
180
)
x
k
+
1
iiiiiiiiiiiii
+
h
[
(
9
-
21
)
f
k
+
98
f
2
+
64
f
c
iiiiiiiiiiiiiiiiiiiii
+
(
9
+
21
)
f
k
+
1
]
}
.A detailed derivation of the formulas above is required when the argument is extended to arbitrary higher orders. In [31], an alternative framework for unifying arbitrary higher-order methods is presented. In this approach, the Legendre-Gauss-Lobatto (LGL) discrete points ξ
j are used to improve both interpolating precision and quadrature performance. Figure 3 graphically illustrates the location of three, five, and seven LGL points. Table 1 lists the corresponding position of these LGL points (the time interval is [
-
1,1
]).Table 1
Locations three, five, and seven Legendre-Gauss-Lobatto points.
Number of points
Location
3
0
±
1
5
0
±
3
7
±
1
7
0
±
495
-
66
15
33
±
495
+
66
15
33
±
1Figure 3
Three, five, and seven Legendre-Gauss-Lobatto points.With reference to Figure3, the blue points are the generic nodes, which are defined by
(46)
τ
j
=
ξ
2
j
-
1
,
j
=
1
,
…
,
(
n
+
1
)
2
,
while the red diamonds are the generic collocation points, denoted as
(47)
ζ
j
=
ξ
2
j
,
j
=
1
,
…
,
(
n
-
1
)
2
,
such that there is one collocation point between every two adjacent nodes [31]. The nodes are used for constructing the interpolation polynomial, while the collocation points are used to formulate the defect constraints (it is noted that nodes and collocation points should not overlap in Gauss-Lobatto method). With three points, the Gauss-Lobatto method degenerates to the Hermite-Simpson method. The state in the ith subinterval is approximated by the nth degree Hermite interpolating polynomial
(48)
x
(
τ
)
≃
a
0
+
a
1
τ
+
a
2
τ
2
+
a
3
τ
3
+
⋯
+
a
n
τ
n
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
τ
∈
[
-
1,1
]
,
where the coefficients of Hermite interpolating polynomial a
=
[
a
0
,
a
1
,
a
2
,
…
,
a
n
]
T are determined by using the values of the states and vector field at the points τ
j; that is,
(49)
[
1
τ
1
1
τ
1
2
⋯
τ
1
n
-
1
τ
1
n
1
τ
2
1
τ
2
2
⋯
τ
2
n
-
1
τ
2
n
⋮
⋮
⋮
⋮
⋮
⋮
1
τ
(
n
+
1
)
/
2
1
τ
(
n
+
1
)
/
2
2
⋯
τ
(
n
+
1
)
/
2
n
-
1
τ
(
n
+
1
)
/
2
n
0
1
2
τ
1
1
⋯
(
n
-
1
)
τ
1
n
-
2
n
τ
1
n
-
1
⋮
⋮
⋮
⋮
⋮
⋮
0
1
2
τ
(
n
+
1
)
/
2
1
⋯
(
n
-
1
)
τ
(
n
+
1
)
/
2
n
-
2
n
τ
(
n
+
1
)
/
2
n
-
1
]
×
[
a
0
a
1
a
2
⋮
a
n
-
1
a
n
]
=
[
x
(
τ
1
)
x
(
τ
2
)
⋮
x
(
τ
(
n
+
1
)
/
2
)
h
2
f
(
τ
1
)
⋮
h
2
f
(
τ
(
n
+
1
)
/
2
)
]
.
Let (49) be written as A
a
=
b. The term A depends only on the number and location of node points, whereas b is a concatenated vector of state and vector field values. Thus, given a vector of states, the coefficients of the Hermite interpolation in the ith subinterval, a, can be found by a
=
[
A
]
-
1
b. The values of the states located at collocation points are then
(50)
x
(
ζ
j
)
=
[
1
ζ
1
ζ
1
2
⋯
ζ
1
n
1
ζ
2
ζ
2
2
⋯
ζ
2
n
⋮
⋮
⋮
⋯
⋮
1
ζ
(
n
-
1
)
/
2
ζ
(
n
-
1
)
/
2
2
⋯
ζ
(
n
-
1
)
/
2
n
]
[
A
]
-
1
b
=
Φ
b
,
j
=
1,2
,
…
,
(
n
-
1
)
2
.
The derivatives of the Hermite interpolating polynomial located at the collocation points are given by
(51)
x
˙
(
ζ
j
)
=
[
0
1
2
ζ
1
⋯
n
ζ
1
n
-
1
0
1
2
ζ
2
⋯
n
ζ
2
n
-
1
⋮
⋮
⋮
⋯
⋮
0
1
2
ζ
(
n
-
1
)
/
2
⋯
n
ζ
(
n
-
1
)
/
2
n
-
1
]
[
A
]
-
1
b
=
Φ
′
b
,
j
=
1,2
,
…
,
(
n
-
1
)
2
.
In this form the matrices Φ and Φ
′ are constants, thus a system of (
n
-
1
)
/
2 constraints per interval is obtained as follows:
(52)
Δ
i
=
x
˙
(
ζ
j
)
-
h
2
f
(
x
(
ζ
j
)
,
u
(
ζ
j
)
)
=
Φ
′
b
-
h
2
f
(
Φ
b
,
u
(
ζ
j
)
)
,
and the analytic Jacobian for the defect constraints can be derived by chain rule as
(53)
∂
Δ
i
∂
[
x
(
τ
j
)
,
u
(
τ
j
)
]
i
=
Φ
′
∂
b
∂
[
x
(
τ
j
)
,
u
(
τ
j
)
]
i
-
h
2
[
∂
f
(
x
(
ζ
j
)
,
u
(
ζ
j
)
)
∂
x
(
ζ
j
)
Φ
∂
b
∂
[
x
(
τ
j
)
,
u
(
τ
j
)
]
i
iiiiiiiiiiiiiii
+
∂
f
(
x
(
ζ
j
)
,
u
(
ζ
j
)
)
∂
u
(
ζ
j
)
∂
u
(
ζ
j
)
∂
[
x
(
τ
j
)
,
u
(
τ
j
)
]
i
]
.
If the final transfer time t
f is a variable, within each segment, the analytic Jacobian of t
f is written as
(54)
∂
Δ
i
∂
t
f
=
Φ
′
∂
b
∂
h
∂
h
∂
t
f
-
[
1
2
∂
h
∂
t
f
f
(
x
(
ζ
j
)
,
u
(
ζ
j
)
)
∂
f
(
x
(
ζ
j
)
,
u
(
ζ
j
)
)
∂
u
(
ζ
j
)
iiiiii
+
h
2
∂
f
(
x
(
ζ
j
)
,
u
(
ζ
j
)
)
∂
x
(
ζ
j
)
Φ
∂
b
∂
h
∂
h
∂
t
f
]
.
A graphical illustration of the Jacobian structure with free final time t
f is given in Figure 4. It can be seen that each segment has the same Jacobian module; thus (53)-(54) can be calculated offline and stored before solving the NLP. In this way, the computational time can be considerably reduced compared to the use of finite difference approximation.Figure 4
Jacobian structure for third-order Gauss-Lobatto method.


## 3.3. Pseudospectral Method

Recently, new direct methods have been applied to low-thrust trajectory optimization, the most popular one being the pseudospectral scheme [32–35]. The major difference between Gauss-Lobatto and pseudospectral collocation schemes is the way in which the interpolation polynomial is constructed and the defect constraints are defined. When the pseudospectral method is employed, a global Lagrange interpolation polynomial is constructed to approximate the state profile. The polynomial is then differentiated and evaluated at all the nodes to compute interpolated derivative values of the states, and equations of motion are used to provide physical time derivatives. The differences between the two sets of time derivatives form the defects. The formulation of defect constraints is sketched in Figure 5.Figure 5
Formulation of the collocation constraints in the pseudospectral method.The defect constraints of the pseudospectral method can be derived as follows. The state curvep
(
t
) is approximated by using values of x at the discrete time points t
i and the corresponding Lagrange interpolating polynomials L
i
(
t
), (
i
=
1
,
…
,
n
)
(55)
p
(
t
)
=
∑
i
=
1
n
‍
L
i
(
t
)
x
(
t
i
)
,
where L
i
(
t
)
(
i
=
1
,
…
,
n
) are defined as
(56)
L
i
(
t
)
=
∏
j
=
1
,
j
≠
i
n
‍
t
-
t
j
t
i
-
t
j
.
Then, the derivative of p
(
t
) is given:
(57)
p
′
(
t
)
=
∑
i
=
1
n
‍
L
i
′
(
t
)
x
(
t
i
)
.
The time derivative of the polynomial can be expressed in compact form through
(58)
p
′
=
D
x
,
where D is differentiation matrix. This only depends on the chosen node spacing, so it can be calculated offline and stored. Then the defect constraints can be written as
(59)
Δ
=
p
′
-
f
(
x
,
u
)
=
D
x
-
f
(
x
,
u
)
.
A variety of pseudospectral schemes emerged and the major differences are the way in which the nodes are selected. The most general form is Legendre-Gauss-Lobatto discretization [36], but Legendre-Gauss [34, 37] and Legendre-Gauss-Radau [35] are also commonly used. Figure 6 illustrates three different kinds of Gauss points. Table 2 lists the position of equidistant points with three kinds of Gaussian points. We can see that LG points are located on the open interval (
-
1,1
), in contrast, LGR points are located on half-open interval [
-
1,1
), and LGL points are located on closed interval [
-
1,1
].Table 2
Position of equidistant, LG, LGR, and LGL with five points within[−1, 1].
Case
Type
x
1
x
2
x
3
x
4
x
5
1
Equidistant
−1
−0.5
0
0.5
1
2
LG points
−0.906
−0.538
0
0.538
0.906
3
LGR points
−1
−0.720
−0.167
0.446
0.886
4
LGL points
−1
−0.654
0
0.655
1Figure 6
Legendre-Gauss (LG), Legendre-Gauss-Radau (LGR), and Legendre-Gauss-Lobatto (LGL) with five points.The main characteristic of the Gauss-Lobatto method is the combination of reasonable accuracy with highly sparse constraint Jacobians and Hessians. The pseudospectral method has a more elegant form of Jacobian and it offers spectral accuracy for smooth problems, but the constraint Jacobian is much denser. So there is a balance between accuracy and efficiency. In general, the Gauss-Lobatto method use a limited number of nodes per segment (n
<
4) with many segments, while pseudospectral methods use many nodes per segment (n
>
10) with a limited number of segments.


## 3.4. Control Parameterization Method

The control parameterization method is a further direct transcription technique for solving optimal control problems [38–40]. In this method, only the controls are discretized and represented by a linear combination of basis functions. This approximation scheme yields a suboptimal control for the original problem. The optimal control problem is converted to an approximate nonlinear optimization problem with a finite number of decision variables, which can be solved by nonlinear programming techniques. The advantage of piecewise-constant approximation scheme is due to its simplicity and convergence [41].To solve optimal control problem using the control parameterization method, we approximate the control profileu as (see Figure 7)
(60)
u
(
t
)
≈
u
p
(
t
)
=
σ
k
,
t
∈
[
τ
k
-
1
,
τ
k
)
,
k
=
1
,
…
,
p
,
where p
≥
1 is a given integer, τ
k, k
=
0
,
…
,
p are knot points, and σ
k
∈
R
r, k
=
1
,
…
,
p are vectors containing the approximate control values. The knot points satisfy
(61)
0
=
τ
0
≤
τ
1
≤
τ
2
≤
⋯
≤
τ
p
-
1
≤
τ
p
=
T
.
The approximate control u
p can be written as
(62)
u
p
(
t
)
=
∑
k
=
1
p
-
1
‍
σ
k
χ
[
τ
k
-
1
,
τ
k
)
(
t
)
+
σ
p
χ
[
τ
p
-
1
,
τ
p
)
(
t
)
,
where for a given subinterval I
⊂
[
0
,
T
], the characteristic function χ
I
:
R
→
R is defined by
(63)
χ
I
∶
=
{
1
,
if
t
∈
I
,
0
,
otherwise
.
Note that u
p is a piecewise-constant function with potential discontinuities at the points t
=
τ
k, k
=
1
,
…
,
p
-
1. These points are calledswitching times. Reference [40] surveys the key developments in the control parameterization.Figure 7
Piecewise-constant control approximation with equidistant knot points.


## 4. Computational Issues

Solving optimal trajectory design problems with direct methods yields a large NLP problem. Some computational issues should be implemented to improve efficiency and robustness. Some of these issues are discussed in this section.


### 4.1. Scaling

An important issue that arises in the solution of NLP isscaling. It is known that a poorly scaled problem can lead to either extremely slow convergence or divergence [27]. When the relative size of variables or constraints in a problem is vastly different, the variables should be transformed into a relatively similar scale. In trajectory optimization problems the scaling process can be implemented manually by introducing dimensionless distance unit L, time unit T, velocity unit V, and mass unit M.


### 4.2. Sparse Matrix

A sparse matrix is a matrix filled primarily with zeros; the opposite is referred to as dense matrix. One of the major advantages of using sparse matrix is that finite memory can be saved. In the case of direct collocation methods, it is well-known that the vast majority of the Jacobian and Hessian matrix elements are zero. When a large number of nodes are taken to have a smooth state and control time histories, “out-of-memory” problems can be avoided by using sparse matrices. Another advantage is that the computational time is reduced because only nonzero elements need to be manipulated. Some NLP solvers, such as SNOPT [42] and IPOPT [43], take advantage of sparse techniques. Furthermore IPOPT computes second derivatives with quasi-Newton methods.


### 4.3. Differentiation

Direct methods use gradient-based techniques for solving the NLP, which requires the gradient information of the objective function and constraints with respect to the NLP variables. It is known that providing the analytic derivative results in exact and fast optimization, but sometimes it is impossible or time consuming to compute derivatives analytically. Alternative methods are needed to obtain the necessary gradients. Finite differences are implemented to compute approximate derivatives. Forward and central difference schemes are(64)
d
f
d
x
≃
f
(
x
+
h
)
-
f
(
x
)
h
,
d
f
d
x
≃
f
(
x
+
h
)
-
f
(
x
-
h
)
2
h
,
where h is a properly chosen perturbation, which should be sufficiently small to provide a good approximation to the derivatives but not too small to induce round-off errors.


## 4.1. Scaling

An important issue that arises in the solution of NLP isscaling. It is known that a poorly scaled problem can lead to either extremely slow convergence or divergence [27]. When the relative size of variables or constraints in a problem is vastly different, the variables should be transformed into a relatively similar scale. In trajectory optimization problems the scaling process can be implemented manually by introducing dimensionless distance unit L, time unit T, velocity unit V, and mass unit M.


## 4.2. Sparse Matrix

A sparse matrix is a matrix filled primarily with zeros; the opposite is referred to as dense matrix. One of the major advantages of using sparse matrix is that finite memory can be saved. In the case of direct collocation methods, it is well-known that the vast majority of the Jacobian and Hessian matrix elements are zero. When a large number of nodes are taken to have a smooth state and control time histories, “out-of-memory” problems can be avoided by using sparse matrices. Another advantage is that the computational time is reduced because only nonzero elements need to be manipulated. Some NLP solvers, such as SNOPT [42] and IPOPT [43], take advantage of sparse techniques. Furthermore IPOPT computes second derivatives with quasi-Newton methods.


## 4.3. Differentiation

Direct methods use gradient-based techniques for solving the NLP, which requires the gradient information of the objective function and constraints with respect to the NLP variables. It is known that providing the analytic derivative results in exact and fast optimization, but sometimes it is impossible or time consuming to compute derivatives analytically. Alternative methods are needed to obtain the necessary gradients. Finite differences are implemented to compute approximate derivatives. Forward and central difference schemes are(64)
d
f
d
x
≃
f
(
x
+
h
)
-
f
(
x
)
h
,
d
f
d
x
≃
f
(
x
+
h
)
-
f
(
x
-
h
)
2
h
,
where h is a properly chosen perturbation, which should be sufficiently small to provide a good approximation to the derivatives but not too small to induce round-off errors.


## 5. Practical Examples

### 5.1. A Simple Optimal Control Problem

In this section, a simple optimal control problem [44] is solved with Hermite-Simpson method, and certain issues are set to speed-up the computation time and improving accuracy. Consider the problem of finding u
(
t
) with fixed initial and final time [
0,1
] to minimize the cost
(65)
J
=
x
2
(
1
)
subject to system equations
(66)
x
˙
1
=
0.5
x
1
+
u
,
x
˙
2
=
u
2
+
x
1
u
+
5
4
x
1
2
and to the boundary conditions
(67)
x
1
(
0
)
=
1
,
x
2
(
0
)
=
0
with the following bound box of state and control variables
(68)
-
10
≤
x
1
(
t
)
≤
10
,
-
10
≤
x
2
(
t
)
≤
10
,
-
10
≤
u
(
t
)
≤
10
.
The analytical solution of this problem is
(69)
x
1
a
(
t
)
=
cosh
⁡
(
1
-
t
)
cosh
⁡
(
1
)
,
u
a
(
t
)
=
-
(
tanh
⁡
(
1
-
t
)
+
0.5
)
cosh
⁡
⁡
(
1
-
t
)
cosh
⁡
⁡
(
1
)
.This problem is solved by using Matlab’s “fmincon” with “interior-point” algorithm [45]. The analytic solution and final optimized results are plotted in Figures 8(a) and 8(b), where the green line indicates analytic solutions (i.e., x
1
a and u
a), and the discrete results obtained from NLP solver are marked by blue asterisk and red cross, respectively. Figure 8(c) shows the error of x
1 and u in detail, it can be seen that the error of x
1 is much lower than that of u (1
0
-
10 for x
1 versus 1
0
-
5 for u). This is because the third-order Hermite interpolation polynomial is implemented to approximate the states, while linear interpolation is used to represent the control curve. In order to check the control accuracy, (66) are integrated from the initial condition with the discrete optimal solution u, and then the error of x
1 is plotted in Figure 8(d). It can be found that the error of x
1 in the whole time domain [
0,1
] is below 1
0
-
5.Simulation results of Example 1, 100 points.
(a)
Discrete solution and analytical time history ofx
1
(b)
Discrete solution and analytical time history ofu
(c)
Error ofx
1 and u located at discrete points
(d)
Reintegrate dynamics with discrete control solutionThe “finite differences” method is the default option to obtain the gradient information, however computational efficiency can be remarkably improved by adding analytical Jacobian and Hessian. The comparison of these three different cases is listed in Table3. All of the three cases use 100 discrete points. In case 1, the optimization takes 116 steps and is 17.21 s. In case 2, adding Jacobian only does not reduce the number of iterations but saves almost half of the CPU time (17.21 s versus 8.77 s). In case 3, providing analytical Jacobin and Hessian highly reduces both the number of iterations and the computational time, which drops dramatically down to 0.38 s.Table 3
Improvement with Jacobian and Hessian, 100 points.
Case
Method
Number of iterations
CPU time
1
Finite differences (default)
116
17.21 s
2
Jacobian only
116
8.77 s
3
Jacobian and Hessian
12
0.38 sThe structure of Jacobian is shown in Figure9(a), where the rows are constraints while the columns are NLP decision variables. The diagonal elements are Jacobian of collocation constraints and the last two lines are Jacobian of initial boundary constraints. Figure 9(b) illustrates the structure of Hessian, where both the row and the column represent NLP decision variables. It can be seen that Hermite-Simpson method has highly sparse Jacobian and Hessian structures.Jacobian and Hessian structures.
(a)
Jacobian structure, 100 points
(b)
Hessian structure, 100 pointsBecause only few elements in Figures9(a) and 9(b) are nonzero, sparse matrices can be used. In Figure 8, the accuracy of state and control solution is not so good. A finer mesh with 5000 points is used with “interior-point” method and matrix sparsity, which is impossible with other fmincon algorithms (that use dense matrixes). The results are illustrated in Figure 10. Compared to Figure 8, it is noted that the error of u reduces to around 1
0
-
8. In order to check the accuracy of control solution, the dynamics are reintegrated with the initial condition and optimal control history; the maximal error on x
1 drops down to 1
0
-
8 (see Figure 10(d)). A parametric analysis has been carried out with variable grid points. The outcome is reported in Table 4.Table 4
Comparison of different grid points.
Case
Points
Number of iterations
CPU time
Maximum error ofu
1
100
12
0.38 s
1.07
e
-
5
2
200
12
0.74 s
2.65
e
-
6
3
500
10
1.82 s
4.24
e
-
7
4
1000
13
4.64 s
1.06
e
-
7
5
2000
17
15.55 s
2.82
e
-
8
6
5000
15
65.04 s
9.16
e
-
9Simulation results of Example 1, 5000 points.
(a)
Discrete solution and analytical time history ofx
1
(b)
Discrete solution and analytical time history ofu
(c)
Error ofx
1 and u located at discrete points
(d)
Reintegrate dynamics with discrete control solution


### 5.2. Planar Low-Thrust Orbit Transfer

For a transfer trajectory to be determined, using Cartesian coordinates is the simplest but most disadvantageous choice; this is because a lot of discrete points will be used to catch the rapidly changing position and velocity variables. On the contrary, slowly changing state variables make the NLP problem both efficient and robust. Polar coordinates are then used for the planar two-body dynamics; that is,(70)
r
˙
=
v
r
,
θ
˙
=
v
t
r
,
v
˙
r
=
v
t
2
r
-
μ
r
2
+
u
sin
ϕ
,
v
˙
t
=
-
v
r
v
t
r
+
u
cos
⁡
ϕ
.
In (70), r is the spacecraft radius, θ is phase angle, v
r and v
t are the radial and transversal velocities, respectively, μ is the gravitational constant, u is the propulsive acceleration, and ϕ is thrust angle (see Figure 11). The initial boundary conditions set the spacecraft in its initial circular orbit at t
i:
(71)
r
(
t
i
)
=
1
,
θ
(
t
i
)
=
0
,
v
r
(
t
i
)
=
0
,
v
t
(
t
i
)
=
1
,
while the final boundary conditions place the spacecraft into a target circular orbit at t
f:
(72)
r
(
t
f
)
=
4
,
v
r
(
t
f
)
=
0
,
v
t
(
t
f
)
=
0.5
.
The nondimensional μ and maximum thrust acceleration u
max
⁡ are
(73)
μ
=
1
,
u
max
⁡
=
0.01
.
In the following, the high efficient NLP solver IPOPT [43] is used for large-scale nonlinear optimization.Figure 11
Two-body orbit transfer in polar coordinates.


#### 5.2.1. Time-Optimal Problem

The goal of time-optimal problem is to find the functionsu
(
t
) and ϕ
(
t
) that minimize the performance index
(74)
J
=
t
f
under the dynamic constraints and the boundary conditions above. This problem is solved with N
=
400 uniformly spaced points with Hermite-Simpson method. (That is, 2401 NLP variables: the states and the controls at the nodes plus the final time.) Figure 12 shows the transfer trajectory; the cyan dashed line is the initial guess (although the convergence basin of direct methods is larger than that of indirect methods, a good initial guess is essential; in time-optimal problem, tangential thrust can be used as a simple but suitable approach to guess an initial solution), the thin blue line denotes the final transfer orbit, and the thick red line shows the thrust arc. Figure 13 shows the time history of state variables. The final transfer time is 55.5. Figure 14 shows the profile of the optimal control variables u, ϕ; it can be seen that the engine is on duty along the entire orbit, to shorten the transfer time.Figure 12
Optimal transfer orbit of time-optimal problem.Figure 13
Time history of state variables.Figure 14
Time history of control variables.


#### 5.2.2. Fuel-Optimal Problem

In the fuel-optimal problem, the functionsu
(
t
) and ϕ
(
t
) are sought to minimize the performance index
(75)
J
=
∫
t
i
t
f
‍
u
d
t
.
Figure 15 illustrates the optimal transfer trajectory. This solution has been obtained with N
=
800 mesh points (i.e., 4801 NLP variables). In this problem, tangential thrust with magnitude of 0.5
u
max
⁡ is used to produce the initial guess (cyan dashed line in Figure 15). It can be seen that the thruster is on duty across the periapsis; the only maneuver performed at the apoapsis injects the spacecraft into a nearly circular orbit to acquire the final orbit. The final time of flight is 122.3. Figure 16 illustrates the time history of the state variables, whereas Figure 17 illustrates the time history of control variables. A bang-bang structure is found, indicating the optimality of the solution found.Figure 15
Optimal transfer orbit of fuel-optimal problem.Figure 16
Time history of state variables.Figure 17
Time history of control variables.


## 5.1. A Simple Optimal Control Problem

In this section, a simple optimal control problem [44] is solved with Hermite-Simpson method, and certain issues are set to speed-up the computation time and improving accuracy. Consider the problem of finding u
(
t
) with fixed initial and final time [
0,1
] to minimize the cost
(65)
J
=
x
2
(
1
)
subject to system equations
(66)
x
˙
1
=
0.5
x
1
+
u
,
x
˙
2
=
u
2
+
x
1
u
+
5
4
x
1
2
and to the boundary conditions
(67)
x
1
(
0
)
=
1
,
x
2
(
0
)
=
0
with the following bound box of state and control variables
(68)
-
10
≤
x
1
(
t
)
≤
10
,
-
10
≤
x
2
(
t
)
≤
10
,
-
10
≤
u
(
t
)
≤
10
.
The analytical solution of this problem is
(69)
x
1
a
(
t
)
=
cosh
⁡
(
1
-
t
)
cosh
⁡
(
1
)
,
u
a
(
t
)
=
-
(
tanh
⁡
(
1
-
t
)
+
0.5
)
cosh
⁡
⁡
(
1
-
t
)
cosh
⁡
⁡
(
1
)
.This problem is solved by using Matlab’s “fmincon” with “interior-point” algorithm [45]. The analytic solution and final optimized results are plotted in Figures 8(a) and 8(b), where the green line indicates analytic solutions (i.e., x
1
a and u
a), and the discrete results obtained from NLP solver are marked by blue asterisk and red cross, respectively. Figure 8(c) shows the error of x
1 and u in detail, it can be seen that the error of x
1 is much lower than that of u (1
0
-
10 for x
1 versus 1
0
-
5 for u). This is because the third-order Hermite interpolation polynomial is implemented to approximate the states, while linear interpolation is used to represent the control curve. In order to check the control accuracy, (66) are integrated from the initial condition with the discrete optimal solution u, and then the error of x
1 is plotted in Figure 8(d). It can be found that the error of x
1 in the whole time domain [
0,1
] is below 1
0
-
5.Simulation results of Example 1, 100 points.
(a)
Discrete solution and analytical time history ofx
1
(b)
Discrete solution and analytical time history ofu
(c)
Error ofx
1 and u located at discrete points
(d)
Reintegrate dynamics with discrete control solutionThe “finite differences” method is the default option to obtain the gradient information, however computational efficiency can be remarkably improved by adding analytical Jacobian and Hessian. The comparison of these three different cases is listed in Table3. All of the three cases use 100 discrete points. In case 1, the optimization takes 116 steps and is 17.21 s. In case 2, adding Jacobian only does not reduce the number of iterations but saves almost half of the CPU time (17.21 s versus 8.77 s). In case 3, providing analytical Jacobin and Hessian highly reduces both the number of iterations and the computational time, which drops dramatically down to 0.38 s.Table 3
Improvement with Jacobian and Hessian, 100 points.
Case
Method
Number of iterations
CPU time
1
Finite differences (default)
116
17.21 s
2
Jacobian only
116
8.77 s
3
Jacobian and Hessian
12
0.38 sThe structure of Jacobian is shown in Figure9(a), where the rows are constraints while the columns are NLP decision variables. The diagonal elements are Jacobian of collocation constraints and the last two lines are Jacobian of initial boundary constraints. Figure 9(b) illustrates the structure of Hessian, where both the row and the column represent NLP decision variables. It can be seen that Hermite-Simpson method has highly sparse Jacobian and Hessian structures.Jacobian and Hessian structures.
(a)
Jacobian structure, 100 points
(b)
Hessian structure, 100 pointsBecause only few elements in Figures9(a) and 9(b) are nonzero, sparse matrices can be used. In Figure 8, the accuracy of state and control solution is not so good. A finer mesh with 5000 points is used with “interior-point” method and matrix sparsity, which is impossible with other fmincon algorithms (that use dense matrixes). The results are illustrated in Figure 10. Compared to Figure 8, it is noted that the error of u reduces to around 1
0
-
8. In order to check the accuracy of control solution, the dynamics are reintegrated with the initial condition and optimal control history; the maximal error on x
1 drops down to 1
0
-
8 (see Figure 10(d)). A parametric analysis has been carried out with variable grid points. The outcome is reported in Table 4.Table 4
Comparison of different grid points.
Case
Points
Number of iterations
CPU time
Maximum error ofu
1
100
12
0.38 s
1.07
e
-
5
2
200
12
0.74 s
2.65
e
-
6
3
500
10
1.82 s
4.24
e
-
7
4
1000
13
4.64 s
1.06
e
-
7
5
2000
17
15.55 s
2.82
e
-
8
6
5000
15
65.04 s
9.16
e
-
9Simulation results of Example 1, 5000 points.
(a)
Discrete solution and analytical time history ofx
1
(b)
Discrete solution and analytical time history ofu
(c)
Error ofx
1 and u located at discrete points
(d)
Reintegrate dynamics with discrete control solution


## 5.2. Planar Low-Thrust Orbit Transfer

For a transfer trajectory to be determined, using Cartesian coordinates is the simplest but most disadvantageous choice; this is because a lot of discrete points will be used to catch the rapidly changing position and velocity variables. On the contrary, slowly changing state variables make the NLP problem both efficient and robust. Polar coordinates are then used for the planar two-body dynamics; that is,(70)
r
˙
=
v
r
,
θ
˙
=
v
t
r
,
v
˙
r
=
v
t
2
r
-
μ
r
2
+
u
sin
ϕ
,
v
˙
t
=
-
v
r
v
t
r
+
u
cos
⁡
ϕ
.
In (70), r is the spacecraft radius, θ is phase angle, v
r and v
t are the radial and transversal velocities, respectively, μ is the gravitational constant, u is the propulsive acceleration, and ϕ is thrust angle (see Figure 11). The initial boundary conditions set the spacecraft in its initial circular orbit at t
i:
(71)
r
(
t
i
)
=
1
,
θ
(
t
i
)
=
0
,
v
r
(
t
i
)
=
0
,
v
t
(
t
i
)
=
1
,
while the final boundary conditions place the spacecraft into a target circular orbit at t
f:
(72)
r
(
t
f
)
=
4
,
v
r
(
t
f
)
=
0
,
v
t
(
t
f
)
=
0.5
.
The nondimensional μ and maximum thrust acceleration u
max
⁡ are
(73)
μ
=
1
,
u
max
⁡
=
0.01
.
In the following, the high efficient NLP solver IPOPT [43] is used for large-scale nonlinear optimization.Figure 11
Two-body orbit transfer in polar coordinates.


### 5.2.1. Time-Optimal Problem

The goal of time-optimal problem is to find the functionsu
(
t
) and ϕ
(
t
) that minimize the performance index
(74)
J
=
t
f
under the dynamic constraints and the boundary conditions above. This problem is solved with N
=
400 uniformly spaced points with Hermite-Simpson method. (That is, 2401 NLP variables: the states and the controls at the nodes plus the final time.) Figure 12 shows the transfer trajectory; the cyan dashed line is the initial guess (although the convergence basin of direct methods is larger than that of indirect methods, a good initial guess is essential; in time-optimal problem, tangential thrust can be used as a simple but suitable approach to guess an initial solution), the thin blue line denotes the final transfer orbit, and the thick red line shows the thrust arc. Figure 13 shows the time history of state variables. The final transfer time is 55.5. Figure 14 shows the profile of the optimal control variables u, ϕ; it can be seen that the engine is on duty along the entire orbit, to shorten the transfer time.Figure 12
Optimal transfer orbit of time-optimal problem.Figure 13
Time history of state variables.Figure 14
Time history of control variables.


### 5.2.2. Fuel-Optimal Problem

In the fuel-optimal problem, the functionsu
(
t
) and ϕ
(
t
) are sought to minimize the performance index
(75)
J
=
∫
t
i
t
f
‍
u
d
t
.
Figure 15 illustrates the optimal transfer trajectory. This solution has been obtained with N
=
800 mesh points (i.e., 4801 NLP variables). In this problem, tangential thrust with magnitude of 0.5
u
max
⁡ is used to produce the initial guess (cyan dashed line in Figure 15). It can be seen that the thruster is on duty across the periapsis; the only maneuver performed at the apoapsis injects the spacecraft into a nearly circular orbit to acquire the final orbit. The final time of flight is 122.3. Figure 16 illustrates the time history of the state variables, whereas Figure 17 illustrates the time history of control variables. A bang-bang structure is found, indicating the optimality of the solution found.Figure 15
Optimal transfer orbit of fuel-optimal problem.Figure 16
Time history of state variables.Figure 17
Time history of control variables.


## 5.2.1. Time-Optimal Problem

The goal of time-optimal problem is to find the functionsu
(
t
) and ϕ
(
t
) that minimize the performance index
(74)
J
=
t
f
under the dynamic constraints and the boundary conditions above. This problem is solved with N
=
400 uniformly spaced points with Hermite-Simpson method. (That is, 2401 NLP variables: the states and the controls at the nodes plus the final time.) Figure 12 shows the transfer trajectory; the cyan dashed line is the initial guess (although the convergence basin of direct methods is larger than that of indirect methods, a good initial guess is essential; in time-optimal problem, tangential thrust can be used as a simple but suitable approach to guess an initial solution), the thin blue line denotes the final transfer orbit, and the thick red line shows the thrust arc. Figure 13 shows the time history of state variables. The final transfer time is 55.5. Figure 14 shows the profile of the optimal control variables u, ϕ; it can be seen that the engine is on duty along the entire orbit, to shorten the transfer time.Figure 12
Optimal transfer orbit of time-optimal problem.Figure 13
Time history of state variables.Figure 14
Time history of control variables.


## 5.2.2. Fuel-Optimal Problem

In the fuel-optimal problem, the functionsu
(
t
) and ϕ
(
t
) are sought to minimize the performance index
(75)
J
=
∫
t
i
t
f
‍
u
d
t
.
Figure 15 illustrates the optimal transfer trajectory. This solution has been obtained with N
=
800 mesh points (i.e., 4801 NLP variables). In this problem, tangential thrust with magnitude of 0.5
u
max
⁡ is used to produce the initial guess (cyan dashed line in Figure 15). It can be seen that the thruster is on duty across the periapsis; the only maneuver performed at the apoapsis injects the spacecraft into a nearly circular orbit to acquire the final orbit. The final time of flight is 122.3. Figure 16 illustrates the time history of the state variables, whereas Figure 17 illustrates the time history of control variables. A bang-bang structure is found, indicating the optimality of the solution found.Figure 15
Optimal transfer orbit of fuel-optimal problem.Figure 16
Time history of state variables.Figure 17
Time history of control variables.


## 6. Conclusions

In this survey, some issues related to the direct transcription and collocation for optimal low-thrust space trajectory design are discussed. The principles of direct transcription and collocation are given, and emphasis is put on high-order numerical integration. Hermite-Simpson, high-order Gauss-Lobatto, and the pseudospectral schemes are introduced. For the Gauss-Lobatto method, a scheme to handle arbitrary order is given, together with the procedure to derive analytical gradients. The pseudospectral method has an impressive convergence rate; however, it is more suitable for smooth problems. A simple example (two states, one control) has been solved with the Hermite-Simpson method. This example shows outstanding efficiency when analytical gradient is provided. A second example consists in low-thrust transfers in the two-body model. From this example it can be concluded that compared with pseudospectral method, the Gauss-Lobatto method is more suitable for fuel-optimal problems with discontinued control structures. In both examples, effort is put to give insights on the computational issues for a practical implementation of the method. However, if a long duration problem is encountered (with increasing number of on/off structures) or a spacecraft equipped with very small thrust engines is dealt with (large number of spirals), the NLP solver becomes ill-conditioned. In these cases, indirect approach with homotopic or hybrid methods could be more suitable.



---

*Source: 851720-2014-06-24.xml*