# Displaced Electric Sail Orbits Design and Transition Trajectory Optimization

Authors: Naiming Qi; Mingying Huo; Qiufan Yuan
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Displaced orbits for spacecraft propelled by electric sails are investigated as an alternative to the use of solar sails. The orbital dynamics of electric sails based spacecraft are studied within a spherical coordinate system, which permits finding the solutions of displaced electric sail orbits and optimize transfer trajectory. Transfer trajectories from Earth's orbit to displaced orbit are also studied in an optimal framework, by using genetic algorithm and Gauss pseudospectral method. The initial guesses for the state and control histories used in the Gauss pseudospectral method are interpolated from the best solution of a genetic algorithm. Numerical simulations show that the electric sail is able to perform the transfer from Earth’s orbit to displaced orbit in acceptable time, and the hybrid optimization method has the capability to search the feasible and optimal solution without any initial value guess.


---

## Body

## 1. Introduction

Displaced orbits are non-Keplerian orbits displaced away from the center of the central body, by applying a continuous thrust to counterbalance the gravity. As a result of the need for continuous thrust to maintain displaced orbits, this mission would be impossible for any conventional space propulsion. Solar sails have often been proposed as spacecraft propulsion system to maintain displaced orbits, because solar sails can produce continuous and propellant-less thrust by using solar radiation pressure. Displaced solar sail orbits have been investigated by various authors for applications to early warning of solar storms and in situ observation [1]. As early as 1929, Oberth mentioned displaced solar sail orbits, which relay on solar radiation pressure to generate a displacement between orbit plane and ecliptic plane. More recently, McInnes and Simmons [2–4] proposed large families of displaced solar sail orbits by considering the dynamics in a rotating frame. The dynamics, stability, and control of different families of displaced orbits were investigated in detail. Gong et al. [5, 6] discussed the coupled attitude-orbit dynamics and control of a solar sail in displaced orbits and analyzed the stability of displaced solar sail orbits with passive control.However, the thrust acceleration level of solar sail is rather low to maintain high displacement orbits [7], because its reflecting membrane surface is not lightweight enough. For this problem, we proposed using electric sails as propulsion system to maintain displaced orbits. The electric sail is an innovative propulsion concept which uses the solar wind dynamic pressure to produce continuous and propellant-less thrust. The electric sail (see Figure 1) consists of many tethers which are connected to a solar-powered electron gun. The electron gun emits electrons to maintain the tethers at a high positive potential. The electric field generated by the tethers reflects the solar wind ions to generate thrust without reaction mass. The expansion and maintenance of the electric sail are implemented by spinning the spacecraft around an axis [8].Figure 1
Conceptual sketch of an electric sail.In recent years, experts have done some research on the theory and application of the electric sail. Janhunen [9–12] proposed the idea of the electric sail and proved that the electric sail is possible to be used in the space exploration through theoretical analysis and experimental research. The research results show that the thrust force of the solar sail decays as (
1
/
r
)
2, and that of the electric sail decays as (
1
/
r
)
7
/
6. Therefore, the thrust force of electric sail decays slower than that of solar sail in the outer solar system exploration mission. Quarta and Mengali [8] applied the electric sail into the outer solar system exploration mission and designed the trajectory of the electric sail through an indirect optimization method. Mengali et al. [13] compared the performance of the electric sail and solar sail. Because electric sails reflect the solar wind protons by electrostatic field, the characteristic acceleration of electric sails can be higher than that of solar sails, which uses reflecting membrane surface to generate thrust. Recent results demonstrate that electric sails can generate 1 N thrust with only 100 kg propulsion system mass.In this paper, displaced orbits using electric sails are presented. For the sake of simplicity of expression, the dynamics of a heliocentric orbiting electric sail are investigated in the context of a spherical coordinate system. In this coordinate system, the solutions of displaced electric sail orbits are discussed. The transfer trajectories from Earth’s orbit to displaced orbit are optimized using a hybrid genetic algorithm Gauss pseudospectral method. This method overcomes the shortcomings of the indirect method, which is sensitive to the initial value of costate variable. The structure of the paper is as follows: dynamical equations are derived in Section2; solutions of displaced electric sail orbits are outlined in Section 3; the implementation of the hybrid method for transfer trajectory optimization is described in Section 4; the numerical results and a conclusion follow in Sections 5 and 6.


## 2. Dynamical Equations

Before the description of the orbital dynamics, two reference frames are defined, the inertial frameo
i
x
i
y
i
z
i and the orbital frame o
o
x
o
y
o
z
o. The origin of inertial frame o
i
x
i
y
i
z
i is at the center-of-mass of the Sun. The z
i-axis is along the normal of the ecliptic plane. The other two axes are in the ecliptic plane and form a right-hand triad. The origin of orbital frame o
o
x
o
y
o
z
o is at the center-of-mass of the sail and the x
o-axis is along the sun-sail direction. The z
o-axis is perpendicular to the x
o-axis and the z
i-axis, and the y
o-axis forms a right-handed triad. As shown in Figure 2, the vector equation of motion for the sail can be written as
(1)
r
¨
+
2
ω
×
r
˙
+
ω
×
(
ω
×
r
)
=
-
(
μ
⊙
r
3
)
r
+
a
s
,
where the position vector r and the angular velocity ω of the sail in frame o
o
x
o
y
o
z
o can be written as [
r
0
0
]
T and [
φ
˙
cos
⁡
θ
φ
˙
sin
θ
θ
˙
]
T, respectively.Figure 2
Heliocentric-ecliptic inertial frame and orbital reference frame.The propelling acceleration vectora
s of the electric sail in frame o
o
x
o
y
o
z
o is given by [8]
(2)
a
s
=
a
c
κ
(
r
⊕
r
)
7
/
6
[
cos
⁡
α
sin
α
cos
⁡
β
sin
α
sin
β
]
T
,
where a
c, referred to as sail characteristic acceleration (in analogy to the solar sail case), is the maximum propelling acceleration at r
=
r
⊕
=
1
AU. Thrust control coefficient κ
∈
[
0,1
], as the thrust of the electric sail can be adjusted by the electron gun. α and β are cone and clock angle, respectively.Substituting (2) in (1), the dynamical equations can be obtained as
(3)
r
¨
-
r
φ
˙
2
si
n
2
θ
-
r
θ
˙
2
=
-
μ
⊙
r
2
+
a
c
κ
(
r
⊕
r
)
7
/
6
cos
⁡
α
,
r
θ
¨
+
2
r
˙
θ
˙
-
r
φ
˙
2
sin
θ
cos
⁡
θ
=
a
c
κ
(
r
⊕
r
)
7
/
6
sin
α
cos
⁡
β
,
r
φ
¨
sin
θ
+
2
r
˙
φ
˙
sin
θ
+
2
r
φ
˙
θ
˙
cos
⁡
θ
=
a
c
κ
(
r
⊕
r
)
7
/
6
sin
α
sin
β
.


## 3. Solutions of Displaced Electric Sail Orbits

Displaced orbits are non-Keplerian orbits displaced away from the center of the central body, by applying a continuous thrust to counterbalance the gravity. Displaced geostationary orbits or displaced non-Keplerian orbits in general can be found by seeking equilibrium solutions to the two-body problem. To obtain the required cone angle, clock angle, and characteristic acceleration, the thrust proportionality coefficientκ is assumed to equal 1 in this section.


### 3.1. General Displaced Orbits

As discussed in [2], the sail orbital period, halo amplitude, and out-of-plane displacement distance of displaced orbits should be constants and can be chosen at will. This also means that the orbital radius, zenithal angle, and azimuthal angle velocity can be chosen to be some particular fixed value. Therefore, the general displaced orbits should have the following characteristics:
(4)
θ
=
θ
d
,
r
=
r
d
,
φ
˙
=
ω
φ
d
,
where θ
d, r
d, and ω
φ
d are constants and are zenithal, orbital radius and azimuthal angular velocity of displaced orbit, respectively.Becauseθ
d, r
d, and ω
φ
d are constant, the time derivative of them should be zero:
(5)
r
˙
d
=
θ
˙
d
=
r
¨
d
=
θ
¨
d
=
φ
¨
d
=
0
.Substituting (4) and (5) in (3), the dynamical equations can be written as
(6)
-
r
d
ω
φ
d
2
si
n
2
θ
d
=
-
μ
⊙
r
d
2
+
a
c
d
(
r
⊕
r
d
)
7
/
6
cos
⁡
α
d
,
-
r
d
ω
φ
d
2
sin
θ
d
cos
⁡
θ
d
=
a
c
d
(
r
⊕
r
d
)
7
/
6
sin
α
d
cos
⁡
β
d
,
0
=
a
c
d
(
r
⊕
r
d
)
7
/
6
sin
α
d
sin
β
d
,
where a
c
d is the required characteristic acceleration of electric sail for the displaced orbit. α
d and β
d are the required cone angle and clock angle, respectively.Letω
*
2
=
μ
⊙
/
r
d
3; the required cone angle, clock angle, and characteristic acceleration of the electric sail to maintain the displaced orbit can be written as
(7)
α
d
=
arctan
(
sin
θ
d
cos
⁡
θ
d
(
ω
*
/
ω
φ
d
)
2
-
si
n
2
θ
d
)
,
β
d
=
arccos
(
sign
⁡
(
-
sin
θ
d
cos
⁡
θ
d
)
)
,
a
c
d
=
r
d
13
/
6
ω
φ
d
2
r
⊕
7
/
6
si
n
2
θ
d
cos
⁡
2
θ
d
+
(
(
ω
*
ω
φ
d
)
2
-
si
n
2
θ
d
)
2
.


### 3.2. Displaced Geostationary Orbits

In this mode of operation, the electric sail orbital period is chosen to be the Earth’s orbital period:(8)
ω
φ
d
=
ω
⊕
=
μ
⊙
r
⊕
3
.Substituting (4) in (6), the required cone angle, clock angle, and characteristic acceleration for displaced geostationary orbits are obtained as
(9)
α
d
=
arctan
(
sin
θ
d
cos
⁡
θ
d
(
r
⊕
/
r
d
)
3
-
si
n
2
θ
d
)
,
(10)
β
d
=
arccos
(
sign
⁡
(
-
sin
θ
d
cos
⁡
θ
d
)
)
,
(11)
a
c
d
=
μ
⊕
r
d
13
/
6
r
⊕
25
/
6
si
n
2
θ
d
cos
⁡
2
θ
d
+
(
(
r
⊕
/
r
d
)
3
-
si
n
2
θ
d
)
2
.As can be seen from (9)–(11), the required cone angle, clock angle, and characteristic acceleration are determined by the displaced orbit parameters θ
d and r
d. The required characteristic acceleration and required cone angle are shown in Figures 3 and 4, respectively. As seen in Figure 3, the required characteristic acceleration increases as displacement increases or the orbital radius decreases.Figure 3
The required characteristic acceleration.Figure 4
The required cone angle.


## 3.1. General Displaced Orbits

As discussed in [2], the sail orbital period, halo amplitude, and out-of-plane displacement distance of displaced orbits should be constants and can be chosen at will. This also means that the orbital radius, zenithal angle, and azimuthal angle velocity can be chosen to be some particular fixed value. Therefore, the general displaced orbits should have the following characteristics:
(4)
θ
=
θ
d
,
r
=
r
d
,
φ
˙
=
ω
φ
d
,
where θ
d, r
d, and ω
φ
d are constants and are zenithal, orbital radius and azimuthal angular velocity of displaced orbit, respectively.Becauseθ
d, r
d, and ω
φ
d are constant, the time derivative of them should be zero:
(5)
r
˙
d
=
θ
˙
d
=
r
¨
d
=
θ
¨
d
=
φ
¨
d
=
0
.Substituting (4) and (5) in (3), the dynamical equations can be written as
(6)
-
r
d
ω
φ
d
2
si
n
2
θ
d
=
-
μ
⊙
r
d
2
+
a
c
d
(
r
⊕
r
d
)
7
/
6
cos
⁡
α
d
,
-
r
d
ω
φ
d
2
sin
θ
d
cos
⁡
θ
d
=
a
c
d
(
r
⊕
r
d
)
7
/
6
sin
α
d
cos
⁡
β
d
,
0
=
a
c
d
(
r
⊕
r
d
)
7
/
6
sin
α
d
sin
β
d
,
where a
c
d is the required characteristic acceleration of electric sail for the displaced orbit. α
d and β
d are the required cone angle and clock angle, respectively.Letω
*
2
=
μ
⊙
/
r
d
3; the required cone angle, clock angle, and characteristic acceleration of the electric sail to maintain the displaced orbit can be written as
(7)
α
d
=
arctan
(
sin
θ
d
cos
⁡
θ
d
(
ω
*
/
ω
φ
d
)
2
-
si
n
2
θ
d
)
,
β
d
=
arccos
(
sign
⁡
(
-
sin
θ
d
cos
⁡
θ
d
)
)
,
a
c
d
=
r
d
13
/
6
ω
φ
d
2
r
⊕
7
/
6
si
n
2
θ
d
cos
⁡
2
θ
d
+
(
(
ω
*
ω
φ
d
)
2
-
si
n
2
θ
d
)
2
.


## 3.2. Displaced Geostationary Orbits

In this mode of operation, the electric sail orbital period is chosen to be the Earth’s orbital period:(8)
ω
φ
d
=
ω
⊕
=
μ
⊙
r
⊕
3
.Substituting (4) in (6), the required cone angle, clock angle, and characteristic acceleration for displaced geostationary orbits are obtained as
(9)
α
d
=
arctan
(
sin
θ
d
cos
⁡
θ
d
(
r
⊕
/
r
d
)
3
-
si
n
2
θ
d
)
,
(10)
β
d
=
arccos
(
sign
⁡
(
-
sin
θ
d
cos
⁡
θ
d
)
)
,
(11)
a
c
d
=
μ
⊕
r
d
13
/
6
r
⊕
25
/
6
si
n
2
θ
d
cos
⁡
2
θ
d
+
(
(
r
⊕
/
r
d
)
3
-
si
n
2
θ
d
)
2
.As can be seen from (9)–(11), the required cone angle, clock angle, and characteristic acceleration are determined by the displaced orbit parameters θ
d and r
d. The required characteristic acceleration and required cone angle are shown in Figures 3 and 4, respectively. As seen in Figure 3, the required characteristic acceleration increases as displacement increases or the orbital radius decreases.Figure 3
The required characteristic acceleration.Figure 4
The required cone angle.


## 4. Transition Trajectory Optimization

The minimum-time electric sail transfer problem from Earth’s orbit to the displaced orbit is studied based on a hybrid genetic algorithm Gauss pseudospectral method. This method overcomes the shortcomings of the indirect method, which is sensitive to the initial value of costate variable.


### 4.1. Problem Description

A spherical coordinate system is established to analyze the heliocentric equations of motion for an electric sail. The dynamical equations of the electric sail take normalized astronomical distance unit (AU) as the unit length and astronomical time unit as the unit time, respectively. The dynamical equations are as follows:(12)
r
˙
=
v
r
,
θ
˙
=
ω
θ
,
φ
˙
=
ω
φ
,
v
˙
r
=
r
ω
φ
2
si
n
2
θ
+
r
ω
θ
2
-
μ
⊙
r
2
+
a
c
κ
(
r
⊕
r
)
7
/
6
cos
⁡
α
,
ω
˙
θ
=
ω
φ
2
sin
θ
cos
⁡
θ
-
2
v
r
ω
θ
r
+
a
c
κ
r
⊕
7
/
6
r
13
/
6
sin
α
cos
⁡
β
,
ω
˙
φ
=
-
2
v
r
ω
φ
r
-
2
ω
θ
ω
φ
cot
θ
+
a
c
κ
r
⊕
7
/
6
sin
α
sin
β
r
13
/
6
sin
θ
.In the trajectory optimization for spacecraft, the objective function to be minimized usually is fuel consumption or transfer time. However, since the electric sail can use the solar wind dynamic pressure to produce thrust without the need for fuel consumption, the minimized transfer time will always be selected as the optimization performance index [8, 13]. In this context of the optimal control problem, the objective function to be minimized is given by
(13)
J
=
t
f
-
t
0
.Assuming that the electric sail is in the Earth’s orbit at the initial timet
0, the initial state of the electric sail should be
(14)
r
(
t
0
)
=
r
⊕
,
θ
(
t
0
)
=
π
2
,
φ
(
t
0
)
=
0
,
v
r
(
t
0
)
=
0
,
ω
θ
(
t
0
)
=
0
,
ω
φ
(
t
0
)
=
ω
⊕
.Assuming that the electric sail arrived at the target displaced orbit at the final timet
f, the final state of the electric sail should be
(15)
r
(
t
f
)
=
r
d
,
θ
(
t
f
)
=
θ
d
,
φ
(
t
f
)
=
ω
⊕
(
t
f
-
t
0
)
,
v
r
(
t
f
)
=
0
,
ω
θ
(
t
f
)
=
0
,
ω
φ
(
t
f
)
=
ω
⊕
,
where r
d and θ
d are parameters of displaced orbit.In the flight of the electric sail, the inequality constraints are as(16)
g
1
=
r
min
⁡
-
r
≤
0
,
g
2
=
α
-
α
max
⁡
≤
0
,
where r
min
⁡ is the minimum allowable distance between the electric sail and the Sun, considering the heat resistance of the electric sail. a
max
⁡ is the maximum allowable thrust angle of the electric sail, which is justified by the necessity of preventing possible mechanical instabilities.


### 4.2. Hybrid Optimization Method

The Gauss pseudospectral method is a direct transcription that transcribes the continuous optimal control problem into a discrete nonlinear programming problem (NLP), which can be solved by well-developed algorithms, for example, SQP. The Gauss pseudospectral method is based on using global polynomial approximations to the dynamic equations at a set of Legendre-Gauss (LG) points. The optimality conditions of the NLP have been found to be equivalent to the discretized optimality conditions of the continuous control problem, which is not true for other pseudospectral methods. This result indicates that the method can take advantage of the properties of both direct and indirect formulations [14].In theory, the trajectory optimization of the electric sail can be implemented based on the Gauss pseudospectral method with sequential quadratic programming (SQP). However, there are following difficulties in the practical problems, when the number of selected Legendre-Gauss (LG) points is large.(
1
) Because the number of design variables is large accordingly in the nonlinear programming (NLP) problem, the feasible initial value is difficult to guess. If the initial value is inappropriate, the NLP problem will converge to an infeasible solution. (
2
) Because the number of equality constraints is large, the optimal solution is difficult to obtain directly based on SQP.So, a hybrid genetic algorithm Gauss pseudospectral optimization method (see Figure5) is proposed to deal with the aforementioned problem. The initial guesses for the state and control histories used in the Gauss pseudospectral method are interpolated from the best solution of a genetic algorithm with less LG points. The proposed hybrid approach includes global search ability of genetic algorithm and the high convergence rate of SQP.Figure 5
The flow chart of electric sail trajectory optimization based on the hybrid optimization method.


### 4.3. Discretization via Gauss Pseudospectral Method

The first step of the Gauss pseudospectral method is changing the time interval from arbitrary bounds[
t
0
,
t
f
] to the interval [
-
1,1
]:
(17)
τ
=
2
t
-
t
f
-
t
0
t
f
-
t
0
.The state variables and control variables are discrete by usingN-order Lagrange interpolation polynomial:
(18)
x
(
τ
)
≈
X
(
τ
)
=
∑
k
=
1
N
x
(
τ
k
)
·
L
k
(
τ
)
,
u
(
τ
)
≈
U
(
τ
)
=
∑
k
=
1
N
u
(
τ
k
)
·
L
k
(
τ
)
,
where τ
k
(
k
=
1
,
…
,
N
) are N LG points. τ
0
=
-
1 and τ
N
+
1
=
1.The differential of state variables and control variables can be approximated as follows:(19)
x
˙
(
τ
i
)
≈
X
˙
(
τ
i
)
=
x
(
τ
0
)
·
D
-
i
+
∑
k
=
1
N
x
(
τ
k
)
·
D
i
k
,
u
˙
(
τ
i
)
≈
U
˙
(
τ
i
)
=
u
(
τ
0
)
·
D
-
i
+
∑
k
=
1
N
u
(
τ
k
)
·
D
i
k
,
where D
-
i
=
L
˙
0
(
τ
i
) and D
i
k
=
L
˙
k
(
τ
i
).The dynamic differential equations of the optimal control problem are transcribed into equality constraints that can be obtained as(20)
D
-
i
·
X
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
X
k
=
t
f
-
t
0
2
f
(
X
k
,
U
k
)
,
i
=
1
,
…
,
N
.The boundary constraints are transcribed via Gauss quadratic integral and can be written as(21)
X
(
τ
f
)
=
X
(
τ
0
)
+
t
f
-
t
0
2
∑
k
=
1
N
ω
k
·
f
(
X
k
,
U
k
)
.So, the dynamic differential equations of the optimal control problem equation (12) are transcribed into equality constraints equation (22). The boundary constraints equation (15) is transcribed into (23). There are 6
N
+
6 equality constraints in the NLP problem:
(22)
h
i
=
D
-
i
·
r
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
r
(
τ
k
)
-
t
f
-
t
0
2
v
r
(
τ
i
)
=
0
,
i
=
1
,
…
,
N
,
h
N
+
i
=
D
-
i
·
θ
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
θ
(
τ
k
)
-
t
f
-
t
0
2
ω
θ
(
τ
i
)
=
0
,
i
=
1
,
…
,
N
,
h
2
N
+
i
=
D
-
i
·
φ
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
φ
(
τ
k
)
-
t
f
-
t
0
2
ω
φ
(
τ
i
)
=
0
,
i
=
1
,
…
,
N
,
h
3
N
+
i
=
D
-
i
·
v
r
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
v
r
(
τ
k
)
-
t
f
-
t
0
2
×
(
a
c
r
⊕
7
/
6
cos
⁡
α
(
τ
i
)
r
(
τ
i
)
7
/
6
r
(
τ
i
)
ω
φ
(
τ
i
)
2
si
n
2
θ
(
τ
i
)
+
r
(
τ
i
)
ω
θ
(
τ
i
)
2
-
μ
⊙
r
(
τ
i
)
2
+
a
c
r
⊕
7
/
6
cos
⁡
α
(
τ
i
)
r
(
τ
i
)
7
/
6
)
=
0
,
i
=
1
,
…
,
N
,
h
4
N
+
i
=
D
-
i
·
ω
θ
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
ω
θ
(
τ
k
)
-
t
f
-
t
0
2
×
(
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
cos
⁡
β
(
τ
i
)
r
(
τ
i
)
13
/
6
ω
φ
(
τ
i
)
2
sin
θ
(
τ
i
)
cos
⁡
θ
(
τ
i
)
-
2
v
r
(
τ
i
)
ω
θ
(
τ
i
)
r
(
τ
i
)
+
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
cos
⁡
β
(
τ
i
)
r
(
τ
i
)
13
/
6
)
=
0
,
i
=
1
,
…
,
N
,
h
5
N
+
i
=
D
-
i
·
ω
φ
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
ω
φ
(
τ
k
)
-
t
f
-
t
0
2
×
(
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
sin
β
(
τ
i
)
sin
θ
(
τ
i
)
r
(
τ
i
)
13
/
6
-
2
v
r
(
τ
i
)
ω
φ
(
τ
i
)
r
(
τ
i
)
-
2
ω
θ
(
τ
i
)
ω
φ
(
τ
i
)
cot
(
θ
(
τ
i
)
)
+
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
sin
β
(
τ
i
)
sin
θ
(
τ
i
)
r
(
τ
i
)
13
/
6
)
=
0
,
i
=
1
,
…
,
N
,
(23)
h
6
N
+
1
=
r
(
τ
0
)
-
r
(
τ
f
)
+
t
f
-
t
0
2
∑
k
=
1
N
w
k
·
v
r
(
τ
k
)
=
0
,
h
6
N
+
2
=
θ
(
τ
0
)
-
θ
(
τ
f
)
+
t
f
-
t
0
2
∑
k
=
1
N
w
k
·
ω
θ
(
τ
k
)
=
0
,
h
6
N
+
3
=
φ
(
τ
0
)
-
φ
(
τ
f
)
+
t
f
-
t
0
2
∑
k
=
1
N
w
k
·
ω
φ
(
τ
k
)
=
0
,
h
6
N
+
4
=
v
r
(
τ
0
)
-
v
r
(
τ
f
)
+
t
f
-
t
0
2
×
∑
k
=
1
N
w
k
·
(
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
sin
β
(
τ
i
)
sin
θ
(
τ
i
)
r
(
τ
i
)
13
/
6
r
(
τ
k
)
ω
φ
(
τ
k
)
2
si
n
2
θ
(
τ
k
)
+
r
(
τ
k
)
ω
θ
(
τ
k
)
2
-
μ
⊙
r
(
τ
k
)
2
+
a
c
r
⊕
7
/
6
κ
(
τ
k
)
cos
⁡
α
(
τ
k
)
r
(
τ
k
)
7
/
6
)
=
0
,
h
6
N
+
5
=
ω
θ
(
τ
0
)
-
ω
θ
(
τ
f
)
+
t
f
-
t
0
2
×
∑
k
=
1
N
w
k
·
(
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
sin
β
(
τ
i
)
sin
θ
(
τ
i
)
r
(
τ
i
)
13
/
6
ω
φ
(
τ
k
)
2
sin
θ
(
τ
k
)
cos
⁡
θ
(
τ
k
)
-
2
v
r
(
τ
k
)
ω
θ
(
τ
k
)
r
(
τ
k
)
+
a
c
r
⊕
7
/
6
κ
(
τ
k
)
sin
α
(
τ
k
)
cos
⁡
β
(
τ
k
)
r
(
τ
k
)
13
/
6
)
=
0
,
h
6
N
+
6
=
ω
φ
(
τ
0
)
-
ω
φ
(
τ
f
)
+
t
f
-
t
0
2
×
∑
k
=
1
N
w
k
·
(
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
sin
β
(
τ
i
)
sin
θ
(
τ
i
)
r
(
τ
i
)
13
/
6
-
2
v
r
(
τ
k
)
ω
φ
(
τ
k
)
r
(
τ
k
)
-
2
ω
θ
(
τ
k
)
ω
φ
(
τ
k
)
cot
(
θ
(
τ
k
)
)
+
a
c
r
⊕
7
/
6
κ
(
τ
k
)
sin
α
(
τ
k
)
sin
β
(
τ
k
)
sin
θ
(
τ
k
)
r
(
τ
k
)
13
/
6
)
=
0
.


### 4.4. Initialization via Genetic Algorithm

As a primary tool to optimize a discrete approximation to a continuous optimal control problem, the GA is not so attractive [15]. However, initialization of the NLP requires only that the solution used for initialization is in the domain of convergence [16]. Because of the global random search ability, the GA is suitable for the initialization of the NLP. The GA applied in this work is based upon the procedure outlined in reference [17].In general, the GA is used to solve unconstrained optimization problem. However the NLP problem in this paper is constrained. The objective function used for ranking not only includes the objective function of the minimum-time problem but should also include the penalties for any constraint violations. The new penalized objective functionJ
GA is as
(24)
J
GA
=
t
f
-
t
0
+
M
(
∑
i
=
1
2
max
⁡
⁡
[
0
,
g
i
]
+
∑
i
=
1
6
N
+
6
|
h
i
|
)
,
where M
=
10,000 is used as a penalty coefficient.


### 4.5. Optimal Solution Calculation via SQP

The initial guesses for the state and control histories used in optimal solution calculation are interpolated from the best solution of the genetic algorithm. Our work presented here was implemented within the MATLAB environment usingfmincon, a constrained optimization routine available with the optimization toolbox. The selected optimization algorithm is the sequential quadratic programming (SQP) algorithm. The NLP problem in the optimal solution calculation is as follows:
(25)
min
⁡
.
J
=
t
f
-
t
0
s
.
t
.
{
h
i
=
0
,
i
=
1
,
…
,
6
N
+
6
g
j
≤
0
,
j
=
1,2
.


## 4.1. Problem Description

A spherical coordinate system is established to analyze the heliocentric equations of motion for an electric sail. The dynamical equations of the electric sail take normalized astronomical distance unit (AU) as the unit length and astronomical time unit as the unit time, respectively. The dynamical equations are as follows:(12)
r
˙
=
v
r
,
θ
˙
=
ω
θ
,
φ
˙
=
ω
φ
,
v
˙
r
=
r
ω
φ
2
si
n
2
θ
+
r
ω
θ
2
-
μ
⊙
r
2
+
a
c
κ
(
r
⊕
r
)
7
/
6
cos
⁡
α
,
ω
˙
θ
=
ω
φ
2
sin
θ
cos
⁡
θ
-
2
v
r
ω
θ
r
+
a
c
κ
r
⊕
7
/
6
r
13
/
6
sin
α
cos
⁡
β
,
ω
˙
φ
=
-
2
v
r
ω
φ
r
-
2
ω
θ
ω
φ
cot
θ
+
a
c
κ
r
⊕
7
/
6
sin
α
sin
β
r
13
/
6
sin
θ
.In the trajectory optimization for spacecraft, the objective function to be minimized usually is fuel consumption or transfer time. However, since the electric sail can use the solar wind dynamic pressure to produce thrust without the need for fuel consumption, the minimized transfer time will always be selected as the optimization performance index [8, 13]. In this context of the optimal control problem, the objective function to be minimized is given by
(13)
J
=
t
f
-
t
0
.Assuming that the electric sail is in the Earth’s orbit at the initial timet
0, the initial state of the electric sail should be
(14)
r
(
t
0
)
=
r
⊕
,
θ
(
t
0
)
=
π
2
,
φ
(
t
0
)
=
0
,
v
r
(
t
0
)
=
0
,
ω
θ
(
t
0
)
=
0
,
ω
φ
(
t
0
)
=
ω
⊕
.Assuming that the electric sail arrived at the target displaced orbit at the final timet
f, the final state of the electric sail should be
(15)
r
(
t
f
)
=
r
d
,
θ
(
t
f
)
=
θ
d
,
φ
(
t
f
)
=
ω
⊕
(
t
f
-
t
0
)
,
v
r
(
t
f
)
=
0
,
ω
θ
(
t
f
)
=
0
,
ω
φ
(
t
f
)
=
ω
⊕
,
where r
d and θ
d are parameters of displaced orbit.In the flight of the electric sail, the inequality constraints are as(16)
g
1
=
r
min
⁡
-
r
≤
0
,
g
2
=
α
-
α
max
⁡
≤
0
,
where r
min
⁡ is the minimum allowable distance between the electric sail and the Sun, considering the heat resistance of the electric sail. a
max
⁡ is the maximum allowable thrust angle of the electric sail, which is justified by the necessity of preventing possible mechanical instabilities.


## 4.2. Hybrid Optimization Method

The Gauss pseudospectral method is a direct transcription that transcribes the continuous optimal control problem into a discrete nonlinear programming problem (NLP), which can be solved by well-developed algorithms, for example, SQP. The Gauss pseudospectral method is based on using global polynomial approximations to the dynamic equations at a set of Legendre-Gauss (LG) points. The optimality conditions of the NLP have been found to be equivalent to the discretized optimality conditions of the continuous control problem, which is not true for other pseudospectral methods. This result indicates that the method can take advantage of the properties of both direct and indirect formulations [14].In theory, the trajectory optimization of the electric sail can be implemented based on the Gauss pseudospectral method with sequential quadratic programming (SQP). However, there are following difficulties in the practical problems, when the number of selected Legendre-Gauss (LG) points is large.(
1
) Because the number of design variables is large accordingly in the nonlinear programming (NLP) problem, the feasible initial value is difficult to guess. If the initial value is inappropriate, the NLP problem will converge to an infeasible solution. (
2
) Because the number of equality constraints is large, the optimal solution is difficult to obtain directly based on SQP.So, a hybrid genetic algorithm Gauss pseudospectral optimization method (see Figure5) is proposed to deal with the aforementioned problem. The initial guesses for the state and control histories used in the Gauss pseudospectral method are interpolated from the best solution of a genetic algorithm with less LG points. The proposed hybrid approach includes global search ability of genetic algorithm and the high convergence rate of SQP.Figure 5
The flow chart of electric sail trajectory optimization based on the hybrid optimization method.


## 4.3. Discretization via Gauss Pseudospectral Method

The first step of the Gauss pseudospectral method is changing the time interval from arbitrary bounds[
t
0
,
t
f
] to the interval [
-
1,1
]:
(17)
τ
=
2
t
-
t
f
-
t
0
t
f
-
t
0
.The state variables and control variables are discrete by usingN-order Lagrange interpolation polynomial:
(18)
x
(
τ
)
≈
X
(
τ
)
=
∑
k
=
1
N
x
(
τ
k
)
·
L
k
(
τ
)
,
u
(
τ
)
≈
U
(
τ
)
=
∑
k
=
1
N
u
(
τ
k
)
·
L
k
(
τ
)
,
where τ
k
(
k
=
1
,
…
,
N
) are N LG points. τ
0
=
-
1 and τ
N
+
1
=
1.The differential of state variables and control variables can be approximated as follows:(19)
x
˙
(
τ
i
)
≈
X
˙
(
τ
i
)
=
x
(
τ
0
)
·
D
-
i
+
∑
k
=
1
N
x
(
τ
k
)
·
D
i
k
,
u
˙
(
τ
i
)
≈
U
˙
(
τ
i
)
=
u
(
τ
0
)
·
D
-
i
+
∑
k
=
1
N
u
(
τ
k
)
·
D
i
k
,
where D
-
i
=
L
˙
0
(
τ
i
) and D
i
k
=
L
˙
k
(
τ
i
).The dynamic differential equations of the optimal control problem are transcribed into equality constraints that can be obtained as(20)
D
-
i
·
X
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
X
k
=
t
f
-
t
0
2
f
(
X
k
,
U
k
)
,
i
=
1
,
…
,
N
.The boundary constraints are transcribed via Gauss quadratic integral and can be written as(21)
X
(
τ
f
)
=
X
(
τ
0
)
+
t
f
-
t
0
2
∑
k
=
1
N
ω
k
·
f
(
X
k
,
U
k
)
.So, the dynamic differential equations of the optimal control problem equation (12) are transcribed into equality constraints equation (22). The boundary constraints equation (15) is transcribed into (23). There are 6
N
+
6 equality constraints in the NLP problem:
(22)
h
i
=
D
-
i
·
r
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
r
(
τ
k
)
-
t
f
-
t
0
2
v
r
(
τ
i
)
=
0
,
i
=
1
,
…
,
N
,
h
N
+
i
=
D
-
i
·
θ
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
θ
(
τ
k
)
-
t
f
-
t
0
2
ω
θ
(
τ
i
)
=
0
,
i
=
1
,
…
,
N
,
h
2
N
+
i
=
D
-
i
·
φ
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
φ
(
τ
k
)
-
t
f
-
t
0
2
ω
φ
(
τ
i
)
=
0
,
i
=
1
,
…
,
N
,
h
3
N
+
i
=
D
-
i
·
v
r
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
v
r
(
τ
k
)
-
t
f
-
t
0
2
×
(
a
c
r
⊕
7
/
6
cos
⁡
α
(
τ
i
)
r
(
τ
i
)
7
/
6
r
(
τ
i
)
ω
φ
(
τ
i
)
2
si
n
2
θ
(
τ
i
)
+
r
(
τ
i
)
ω
θ
(
τ
i
)
2
-
μ
⊙
r
(
τ
i
)
2
+
a
c
r
⊕
7
/
6
cos
⁡
α
(
τ
i
)
r
(
τ
i
)
7
/
6
)
=
0
,
i
=
1
,
…
,
N
,
h
4
N
+
i
=
D
-
i
·
ω
θ
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
ω
θ
(
τ
k
)
-
t
f
-
t
0
2
×
(
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
cos
⁡
β
(
τ
i
)
r
(
τ
i
)
13
/
6
ω
φ
(
τ
i
)
2
sin
θ
(
τ
i
)
cos
⁡
θ
(
τ
i
)
-
2
v
r
(
τ
i
)
ω
θ
(
τ
i
)
r
(
τ
i
)
+
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
cos
⁡
β
(
τ
i
)
r
(
τ
i
)
13
/
6
)
=
0
,
i
=
1
,
…
,
N
,
h
5
N
+
i
=
D
-
i
·
ω
φ
(
τ
0
)
+
∑
k
=
1
N
D
i
k
·
ω
φ
(
τ
k
)
-
t
f
-
t
0
2
×
(
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
sin
β
(
τ
i
)
sin
θ
(
τ
i
)
r
(
τ
i
)
13
/
6
-
2
v
r
(
τ
i
)
ω
φ
(
τ
i
)
r
(
τ
i
)
-
2
ω
θ
(
τ
i
)
ω
φ
(
τ
i
)
cot
(
θ
(
τ
i
)
)
+
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
sin
β
(
τ
i
)
sin
θ
(
τ
i
)
r
(
τ
i
)
13
/
6
)
=
0
,
i
=
1
,
…
,
N
,
(23)
h
6
N
+
1
=
r
(
τ
0
)
-
r
(
τ
f
)
+
t
f
-
t
0
2
∑
k
=
1
N
w
k
·
v
r
(
τ
k
)
=
0
,
h
6
N
+
2
=
θ
(
τ
0
)
-
θ
(
τ
f
)
+
t
f
-
t
0
2
∑
k
=
1
N
w
k
·
ω
θ
(
τ
k
)
=
0
,
h
6
N
+
3
=
φ
(
τ
0
)
-
φ
(
τ
f
)
+
t
f
-
t
0
2
∑
k
=
1
N
w
k
·
ω
φ
(
τ
k
)
=
0
,
h
6
N
+
4
=
v
r
(
τ
0
)
-
v
r
(
τ
f
)
+
t
f
-
t
0
2
×
∑
k
=
1
N
w
k
·
(
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
sin
β
(
τ
i
)
sin
θ
(
τ
i
)
r
(
τ
i
)
13
/
6
r
(
τ
k
)
ω
φ
(
τ
k
)
2
si
n
2
θ
(
τ
k
)
+
r
(
τ
k
)
ω
θ
(
τ
k
)
2
-
μ
⊙
r
(
τ
k
)
2
+
a
c
r
⊕
7
/
6
κ
(
τ
k
)
cos
⁡
α
(
τ
k
)
r
(
τ
k
)
7
/
6
)
=
0
,
h
6
N
+
5
=
ω
θ
(
τ
0
)
-
ω
θ
(
τ
f
)
+
t
f
-
t
0
2
×
∑
k
=
1
N
w
k
·
(
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
sin
β
(
τ
i
)
sin
θ
(
τ
i
)
r
(
τ
i
)
13
/
6
ω
φ
(
τ
k
)
2
sin
θ
(
τ
k
)
cos
⁡
θ
(
τ
k
)
-
2
v
r
(
τ
k
)
ω
θ
(
τ
k
)
r
(
τ
k
)
+
a
c
r
⊕
7
/
6
κ
(
τ
k
)
sin
α
(
τ
k
)
cos
⁡
β
(
τ
k
)
r
(
τ
k
)
13
/
6
)
=
0
,
h
6
N
+
6
=
ω
φ
(
τ
0
)
-
ω
φ
(
τ
f
)
+
t
f
-
t
0
2
×
∑
k
=
1
N
w
k
·
(
a
c
r
⊕
7
/
6
κ
(
τ
i
)
sin
α
(
τ
i
)
sin
β
(
τ
i
)
sin
θ
(
τ
i
)
r
(
τ
i
)
13
/
6
-
2
v
r
(
τ
k
)
ω
φ
(
τ
k
)
r
(
τ
k
)
-
2
ω
θ
(
τ
k
)
ω
φ
(
τ
k
)
cot
(
θ
(
τ
k
)
)
+
a
c
r
⊕
7
/
6
κ
(
τ
k
)
sin
α
(
τ
k
)
sin
β
(
τ
k
)
sin
θ
(
τ
k
)
r
(
τ
k
)
13
/
6
)
=
0
.


## 4.4. Initialization via Genetic Algorithm

As a primary tool to optimize a discrete approximation to a continuous optimal control problem, the GA is not so attractive [15]. However, initialization of the NLP requires only that the solution used for initialization is in the domain of convergence [16]. Because of the global random search ability, the GA is suitable for the initialization of the NLP. The GA applied in this work is based upon the procedure outlined in reference [17].In general, the GA is used to solve unconstrained optimization problem. However the NLP problem in this paper is constrained. The objective function used for ranking not only includes the objective function of the minimum-time problem but should also include the penalties for any constraint violations. The new penalized objective functionJ
GA is as
(24)
J
GA
=
t
f
-
t
0
+
M
(
∑
i
=
1
2
max
⁡
⁡
[
0
,
g
i
]
+
∑
i
=
1
6
N
+
6
|
h
i
|
)
,
where M
=
10,000 is used as a penalty coefficient.


## 4.5. Optimal Solution Calculation via SQP

The initial guesses for the state and control histories used in optimal solution calculation are interpolated from the best solution of the genetic algorithm. Our work presented here was implemented within the MATLAB environment usingfmincon, a constrained optimization routine available with the optimization toolbox. The selected optimization algorithm is the sequential quadratic programming (SQP) algorithm. The NLP problem in the optimal solution calculation is as follows:
(25)
min
⁡
.
J
=
t
f
-
t
0
s
.
t
.
{
h
i
=
0
,
i
=
1
,
…
,
6
N
+
6
g
j
≤
0
,
j
=
1,2
.


## 5. Numerical Simulations

In this section, we present a geostationary displaced orbit (r
d
=
0.8
AU, θ
d
=
8
0
∘), which can be applied into early warning of solar storms and interstellar communication missions. By (9) and (11), we can obtain the required characteristic acceleration a
c
d
=
3.65
mm
/
s
2 and required cone angle α
d
=
9.8
7
∘. Assuming using an electric sail of 1 N thrust and 100 kg mass [9], the payload mass can reach 174 kg, and payload mass fraction can reach 63.5
% in the case of meeting the requirements of characteristic acceleration. However, the required characteristic acceleration is difficult to achieve for solar sail. Assuming using a solar sail of (139 m)2 area, 160 mN thrust, 439 kg propulsion system mass [18], and 174 kg payload mass, the characteristic acceleration of this solar sailcraft is 0.261 mm/s2, which cannot meet the requirements of characteristic acceleration.The transition trajectory between displaced orbit and Earth’s orbit is optimized to investigate the orbital transfer capability of the electric sail and to verify the hybrid optimization method. In this paper, “the orbital transfer capability of the electric sail” means the independent ability of the electric sail to transfer from Earth’s orbit to displaced orbit in reasonable transfer time. The simulation parameters of the electric sail are shown in Table1. In the initial solution calculation via GA, the number of LG points is 12. The population size is 100, and the number of generations is 200. In the optimal solution calculation via SQP, the number of LG points is 60. The convergence criteria for the NLP are a maximum constraint violation of 1
×
1
0
-
9 (TolCon) and a reduction in the cost function of less than 1
×
1
0
-
6 (TolFun) for one iteration. The complete framework is implemented in MATLAB (R2010a) on a personal computer (Intel Core i3) with 2 GB RAM running at 2.67 GHz processor speed.Table 1
Simulation parameters of the electric sail.
Parameters
Value
a
⊕
3.65 mm/s2
α
max
⁡
3
5
°
r
min
⁡
0.5 AUThe optimal solution of the transition trajectory is shown in Figure6. The state histories and control histories are shown in Figures 7 and 8, respectively. The numerical results show that the minimum transfer time is 216.06 days using an electric sail, whose characteristic acceleration is 3.65
mm
/
s
2 and maximum allowable thrust angle is 3
5
∘. The optimal solution meets the boundary constraints very well, and the proposed hybrid optimization method is suitable for the electric sail trajectory optimization.Figure 6
Transition trajectory between Earth’s orbit and displaced orbit.State histories.
(a)
Orbital radius history
(b)
Zenithal history
(c)
Azimuthal history
(d)
Radial velocity history
(e)
Zenith angular velocity history
(f)
Azimuth angular velocity historyControl histories.
(a)
Cone angleα history
(b)
Clock angleβ history
(c)
Proportionality coefficientκ historyAs we can see in Figure7(a), the orbital radius is not monotonically decreased from 1 AU to 0.8 AU but firstly increases from 1 AU to about 1.05 AU for about 50 days and then decreases to 0.8 AU. The reason of this phenomenon is that the electric sail should reduce the azimuth angular velocity by increasing the orbit radius firstly, so that it can satisfy the terminal constraints that azimuth angular and azimuth angular velocity of the electric sail need to be consistent with the Earth. If orbital radius monotonically decreases from 1 AU to 0.8 AU, azimuth angular velocity will also monotonically increase, and the electric sail would be difficult to synchronize with the Earth at last. According to simulation result, we can roughly divide the transition process into three stages. Firstly, from 0th to 60th day, the electric sail mainly reduces its azimuth angular velocity. Then, from 60th to 150th day, the main task of the electric sail is to reduce orbital radius. Finally, from 150th day to the end, zenithal angle decreases to 80 deg, and all terminal constraints are satisfied at last.In the transition process, the attitude of the electric sail needs to be adjusted. To this end, a tunable resistor is placed between the spacecraft and each tether. Because the thrust magnitude depends on the tether potential, the latter can be slightly changed through the potentiometer. As a result, the tether spin plane can be turned by modulating the potentiometer setting with a sinusoidal signal synchronized to the sail rotation [13]. The attitude control of electric sails will be studied in future work.


## 6. Conclusion

In this paper, displaced orbits for spacecraft propelled by electric sails are investigated as an alternative to the use of solar sails. In order to find solutions of displaced orbits and optimize transition trajectories, the dynamics of heliocentric orbiting electric sails are investigated in the context of a spherical coordinate system. In this coordinate system, the solutions of general and geostationary displaced orbits are discussed. By comparing characteristic acceleration between electric sail and solar sail, we found that the electric sail is more suited to high displacement orbits. The transfer trajectory from Earth’s orbit to displaced orbit is also studied in an optimal framework, by using genetic algorithm and Gauss pseudospectral method. The initial guesses for the state and control histories used in the Gauss pseudospectral method are interpolated from the best solution of a genetic algorithm. Numerical simulations show that the electric sail is able to perform the transfer from Earth’s orbit to displaced orbit in acceptable time, and the hybrid optimization method has the capability to search the feasible and optimal solution without any initial value guess. Although the electric sail still deserves more refined theoretical and experimental studies, the proposed study represents a first step toward substantiating this propulsion concept in terms of preliminary mission design.



---

*Source: 932190-2014-06-03.xml*