Pulse-Width Predictive Control for LTV Systems with
Application to Spacecraft Rendezvous

arXiv:1511.00869v1 [math.OC] 3 Nov 2015

R. Vazqueza,∗, F. Gavilana , E. F. Camachob
a Departamento

de Ingenierı́a Aeroespacial. Escuela Superior de Ingenieros. Universidad de
Sevilla, Camino de los Descubrimientos s/n, 41092,Sevilla, Spain
b Departamento de Ingenierı́a de Sistemas y Automática. Escuela Superior de Ingenieros.
Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092,Sevilla, Spain

Abstract
This work presents a model predictive controller (MPC) that is able to handle
linear time-varying (LTV) plants with PWM control. The MPC is based on
a planner that employs a PAM or impulsive approximation as a hot-start and
then uses explicit linearization around successive PWM solutions for rapidly
improving the solution by means of linear programming. As an example, the
problem of rendezvous of spacecraft for eccentric target orbits is considered. The
problem is modeled by the LTV Tschauner-Hempel equations, whose transition
matrix is explicit; this is exploited by the algorithm for rapid convergence. The
efficacy of the method is shown in a simulation study.
Keywords: Spacecraft autonomy, Space robotics, Pulse-width modulation,
Trajectory planning, Optimal trajectory, Linear Time-Varying Systems.

1. Introduction
Aerospace systems often need to be controlled by using pulse-width modulated (PWM) actuators, i.e., actuators whose output level is fixed and can
only be turned on and off, such as spacecraft thrusters. It would be therefore
desirable to use control design methods that directly take into account pulsed
actuators. However, most feedback design and motion planning methods ignore
∗ Corresponding

author, Fax number +34954486041, Email address: rvazquez1@us.es
Email addresses: fgavilan@us.es (F. Gavilan), eduardo@esi.us.es (E. F. Camacho)

Preprint submitted to Control Engineering Practice

February 22, 2018

variable width pulses and approximate the control variables either by impulses
(which produce instantaneous changes in some combination of the states) or
pulse-amplitude modulated (PAM) control. However, neither impulsive actuation nor PAM actuation capture with precision the behavior of pulsed actuators
such as spacecraft thrusters. A more realistic model has to take into account
that, typically, thrusters are ON-OFF actuators, i.e., the thrusters are not able
to produce arbitrary forces, but instead can only be switched on (producing
the maximum amount of force) or off (producing no force). These switching
times are the only signals that can be controlled. This type of control signal
is usually referred to as Pulse-Width Modulated (PWM). Control design with
PWM actuation poses a challenge because the system becomes nonlinear in the
switching times, even if the system is linear.
One can find in the literature several procedures to find an equivalent PWM
solution starting from a PAM solution (for instance in Shieh et al. (1996); Ieko
et al. (1999); Bernelli-Zazzera et al. (1998)). These methods allow to, given the
PAM inputs of a system, compute PWM inputs that produce a system output
optimally approximating the output of the system when driven by the PAM
signals. The results are based on the so-called Principle of Equivalent Areas,
which computes the PWM signal so that it covers the same area as the PAM
signal. However, while these procedures are quite effective in the sense that the
output produced by the approximate PWM signals is very similar to the one
produced by PAM signals, they assume that the plant is linear time-invariant.
In this paper, Model Predictive Control (MPC) is used to directly find PWM
signals to control the system. MPC (see, e.g., Camacho and Bordons (2004))
is a family of methods that originated in the late seventies and has developed
considerably since then. In MPC, the process model is used to predict the future
plant outputs, based on past and current values and on the proposed optimal
future control actions. These actions are calculated by the optimizer taking into
account the cost function as well as the constraints. Since the plant is nonlinear
in the control signals (ON-OFF times), the underlying optimization problem is
nonlinear and possibly non-convex. To solve the problem, the algorithm starts
2

from an initial guess computed by solving an optimal linear program with PAM
or impulsive actuation, approximate the solution with ON-OFF thrusters, and
then iteratively linearize around the obtained solutions to improve the PWM
solution. While the idea of linearization to specifically compute optimal PWM
control signals in the context of MPC is, to the best knowledge of the authors,
original, it must be noted that local linearization techniques have been used for
optimal trajectory problems in other contexts (see e.g. Kim et al. (2002)).
As an application the problem of rendezvous of spacecraft is considered, i.e.,
the controlled close encounter of two space vehicles. Autonomous spacecraft
rendezvous capabilities are becoming a necessity as access to space continues
increasing. The field has become very active in recent years, with a rapidly
growing literature. Among others, approaches based on trajectory planning and
optimization (Breger and How (2008); Arzelier et al. (2013, 2011); Louembet
et al. (2015); Deaconu et al. (2015, 2014); D’Amico et al. (2013); Gaias et al.
(2014)) and predictive control (Richards and How (2003); Rossi and Lovera
(2002); Asawa et al. (2006); Gavilan et al. (2009, 2012); Larsson et al. (2006);
Hartley et al. (2012); Leomanni et al. (2014); Jewison et al. (2015); Weiss et al.
(2012)) are emerging.
Classically, in these approaches the problem of rendezvous is modeled by
using impulsive maneuvers; one computes a sequence of (possibly optimal) impulses (usually referred to as ∆V ’s) to achieve rendezvous.
Recently, Vazquez et al. (2011, 2014) introduced a trajectory planning algorithm algorithm for spacecraft rendezvous that was able to incorporate PWM
control signals. The former considered the linear time-invariant Clohessy-Wiltshire
model (target orbiting in a circular Keplerian orbit, see Clohessy and Wiltshire
(1960)). The latter extended the approach to elliptical target orbits by using
the linear time-varying Tschauner-Hempel model (see Tschauner and Hempel
(1965)). Both methods start from an initial guess computed by solving an optimal linear program with PAM or impulsive actuation, approximate the solution
with ON-OFF thrusters, and then iteratively linearize around the obtained solutions to improve the PWM solution. For both circular and elliptical target
3

orbits the algorithms are simple and reasonably fast, and simulations favorably compare with an impulsive-only approach. These results were extended
in Vazquez et al. (2015) to a decreasing-horizon model predictive controller able
to take into account orbital perturbations, disturbances or model errors.
In this paper, a receding horizon model predictive controller with PWM
inputs is formulated for general LTV plants and both alternatives (PAM or impulsive starting guess) are discussed in detail, with an application to rendezvous
given at the end of the paper.
The structure of the paper is as follows. In Section 2 the plant model is
introduced. Three types of inputs are considered: PWM, PAM and impulsive.
Section 3 follows with a formulation of the underlying optimization problem.
Section 4 describes a method that solves the planning problem using PWM
signals. Section 5 develops the model predictive controller. Next, Section 6 describes the application to spacecraft rendezvous. Section 7 presents a simulation
study of the method applied to spacecraft rendezvous. The paper finishes with
some remarks in Section 8.

2. System Model
Consider a linear time-varying system given as
ẋ = A(t)x + B(t)u,

(1)

where x ∈ Rn is the state, u ∈ Rm is the input (control) vector, and A(t) and
B(t) are, respectively, n × n and n × m matrices depending on time t ≥ 0.
Considering that, for some time tk ≥ 0, initial conditions x(tk ) ∈ Rn are
given and the input is known, the solution to (1) for t > tk is given by
Z t
x(t) = Φ(t, tk )x(tk ) +
Φ(t, s)B(s)u(s)ds,

(2)

tk

where Φ(t, tk ) is the system transition matrix, see for instance Rugh (1996).
This matrix can be computed numerically (or analytically if possible) as the

4

unique solution to the linear matrix differential equation
Φ̇(t, tk )
Φ(tk , tk )

= A(t)Φ(t, tk ),
=

t > tk

(3)

I.

(4)

To obtain an unified notation in term of the inputs, denote by Bi (t) the i-th
column of B(t), corresponding to the i-th input ui (t), for i = 1, . . . , m. In the
paper, time intervals starting at some initial time tk and ending at tk+1 = tk +T
are considered, where T will be an adequate sample time. Then equation (2)
can be written as
x(tk+1 )

=

Φ(tk+1 , tk )x(tk ) +

m Z
X

tk+1

Φ(tk+1 , s)Bi (s)ui (s)ds,

(5)

tk

i=1

The objective is solving the problem with PWM inputs. In addition, two
other types of inputs are considered; they will be used as an intermediate step
towards computing PWM inputs by the algorithm. All types of input are analyzed in the following sections.
2.1. Pulse width-modulated (PWM) control
In the PWM case, each input ui is a pulse starting at time τk,i (relative to
tk ) with pulse width κk,i , with constant magnitude umax = uW
k,i , as shown in
Fig. 1, i.e.,


0,



ui(t) =
uW
k,i ,



 0,

t ∈ [tk , tk + τk,i ] ,
(6)

t ∈ [tk + τk,i , tk + τk,i + κk,i ] ,
t ∈ [tk + τk,i + κk,i , tk+1 ] ,

with κk,i > 0, τk,i > 0 and τk,i + κk,i < T , where the last constraint prevent the
PWM signal to spill over to the next time interval. Then, substituting ui (t) in
(5) one obtains
x(tk+1 )

=

Φ(tk+1 , tk )x(tk ) +

Z
m
X
i=1

tk +τk,i +κk,i

!

Φ(tk+1 , s)Bi (s)ds uW
k,i ,

(7)

tk +τk,i

and denoting
W
Bk,i
(τk,i , κk,i ) =

Z

tk +τk,i +κk,i

Φ(tk+1 , s)Bi (s)ds,
tk +τk,i

5

(8)

T
umax
·
¿
t
Figure 1: PWM Variables.

one can write the solution as
x(tk+1 )

=

Φ(tk+1 , tk )x(tk ) +

m
X

W
Bk,i
(τk,i , κk,i )uW
k,i .

(9)

i=1

There is an important difference between a PWM input and the PAM or impulsive inputs that will be subsequently introduced. While the latter can in
principle take positive o negative values at different times, the former is fixed
either as positive or negative for all time. Thus, typically a PWM model has
twice number of inputs than a PAM/impulsive model. To make this explicit in
the model (10), denote with a plus or minus super-index the positive or negative
inputs, as follows
m h
i
X
+
W+
−
−
W−
W
W
x(tk+1 ) = Φ(tk+1 , tk )x(tk )+
Bk,i
(τk,i
, κ+
)u
−
B
(τ
,
κ
)u
,
k,i k,i
k,i k,i
k,i k,i
i=1

(10)
+ +
+
W− −
−
with uW
k,i , τk,i , κk,i and uk,i , τk,i , κk,i denoting, respectively, the magnitude,

start, and width of the positive and negative i-th input pulses.
2.2. Pulse amplitud-modulated (PAM) control
In this case, each control ui (t) in (5) is constant inside the interval [tk , tk+1 ],
and equal to uA
k,i . Then, substituting ui (t) in (5) one obtains
x(tk+1 )

=

Φ(tk+1 , tk )x(tk ) +

m Z
X
i=1

tk+1


Φ(tk+1 , s)Bi (s)ds uA
k,i ,

(11)

tk

and denoting by
A
Bk,i
=

Z

tk+1

Φ(tk+1 , s)Bi (s)ds,
tk

6

(12)

the solution can be written as
x(tk+1 )

=

Φ(tk+1 , tk )x(tk ) +

m
X

A A
Bk,i
uk,i .

(13)

i=1

2.3. Impulsive control
In this case, ui (t) = uIk,i δ(t − (tk + τk,i )), where δ(t) is Dirac’s delta function,
tk + τk,i is the instant at which the impulse is given, and uk,i is the magnitude
of the impulse. Then, assuming 0 < τk,i < T for all i (all the impulses are given
inside the considered time interval) and substituting ui (t) in (5) one obtains
x(tk+1 )

=

Φ(tk+1 , tk )x(tk ) +

m
X

Φ(tk+1 , tk + τi )Bi (tk + τi )uIk,i ,

(14)

i=1
I
and denoting by Bk,i
(τk,i ) = Φ(tk+1 , τi )Bi (tk + τi ),

x(tk+1 )

=

Φ(tk+1 , tk )x(tk ) +

m
X

I
Bk,i
(τk,i )uIk,i .

(15)

i=1

2.4. Discretization and compact notation
Consider now a sequence of time instants tk = t0 +kT , k = 0, . . ., and denote
xk = x(tk ). Then, it is possible to write, for both PAM and impulsive control,
xk+1 = Ak xk + Bk Uk ,

(16)

where Ak = Φ(tk+1 , tk ), and Bk and Uk depend on the input type. In the PWM
case, write
xk+1 = Ak xk + Bk+ Uk+ − Bk− Uk− ,

(17)

+
+
+
W
where Bk+ is a matrix whose i-th column is Bk,i
(τk,i
, κ+
k,i ) and Uk , τk , and
+
+
W
κ+
k are column vectors whose i-th entries are, respectively, uk,i , τk,i and κk,i .

The same definitions (with minus super-index) are used for the negative pulses.
Then, to reach model (16), define






+
+
+
h
i
τ
κ
U
Bk = Bk+ −Bk− , Uk =  k  , τk =  k  , κk =  k  .
τk−
κ−
Uk−
k

(18)

The definitions of Bk are simpler for the other types of actuation. In the
A
PAM case, Bk is a matrix whose i-th column is Bk,i
and Uk a column vector

7

whose i-th entry is uA
k,i . In the impulsive case, Bk is a matrix whose i-th column
I
is Bk,i
(τk,i ), and Uk , τk are column vector whose i-th entries are, respectively,

uIk,i and τk,i .
Next a compact formulation is developed to simplify the notation of the
problem. The state at time tk+j+1 , given the state xk at time tk , and the input
signals from tk to time tk+j , is computed by applying recursively, in the PAM
and impulsive cases, by applying Equation (16):
xk+j+1

= Ak+j+1,k xk +

k+j
X

Ak+j,i Bk,i Uk,i ,

(19)

i=k

where the definition Ak,k = I, Ak+1,k = Ak , and if j > 0 then Ak+j+1,k =
Ak+j Ak+j−1 . . . Ak has been used. Define now Xk and Uk as a stack vector of
Np state and input vectors, respectively, spanning from time tk+1 to time tk+Np
for the state and from time tk to time tk+Np −1 for the controls, where Np is the
planning horizon:

Xk



= 


xk+1
..
.









 , Uk = 



xk+Np



Uk
..
.



.


Uk+Np −1

Similarly, for the impulsive and PWM cases, define



κk
τk



..
.



..
Γk = 
 , Λk = 
.



κk+Np −1
τk+Np −1





Γ
k


 , Υk = 

Λk

Then one can write
Xk = Fk xk + Gk Uk ,

(20)

where Gk is a square, block lower triangular matrix of size mNp , defined as


0

···

0



Bk+1
..
.

···
..
.

0
..
.




,




Ak+Np ,k+2 Bk+1

...

Bk+N p−1

Bk



 Ak+2,k+1 Bk
Gk = 
..


.

Ak+Np ,k+1 Bk

8

(21)

this is, its non-null blocks are defined by (Gk )jl = Ak+j,k+l Bk+l−1 , and the
matrix Fk is defined as:


Ak+1,k



 Ak+2,k
Fk = 
..


.

Ak+Np ,k





.




(22)

It is important to note that, in the impulsive case, Gk is a (nonlinear) function of Γk , whereas in the PWM case it is a (nonlinear) function of Υk . To
avoid lengthy expressions this dependence has been omitted. Another important remark is that, in the PWM case, Uk is fixed whereas in the other cases is
the input variable

3. Formulation of the planning problem
Next the planning problem is formulated, introducing the constraints and
the objective function. The formulation is done for the three types of control
signals.
3.1. Constraints on the problem
First constraints on the state and input are introduced. While only inequality are considered, equality constraints would be treated similarly.
3.1.1. Inequality constraints on the state
In this work it is assumed that the state is subject to inequality constraints
along the planning horizon, which can vary as time advances. These can be
formulated in general as Ak Xk ≤ bk , and using (20), one reaches a expression
in term of inputs, namely
Ak Gk Uk ≤ bk − Ak Fk xk .

9

(23)

3.1.2. Input constraints
Input constraints are different depending on the type of input.
In the PWM case, the inputs Uk are fixed, but the start time of impulse,
tk + τk , and its end, tk + τk + κk , must be within the time interval both for
negative and positive pulses. Thus
0

≤ Γk ,

(24)

0

≤ Λk ,

(25)

≤ T,

(26)

Γk + Λ k
which can be summarized as

AW Υk ≤ bW .

(27)

In the PAM case, the inputs are limited above and below. Thus
U P AM ≤ Uk ≤ U P AM

(28)

In the impulsive case, the inputs uk are limited above and below, but also
the times of impulse, tk + τk , must be within the time interval. Thus
U IM P
0

≤

Uk ≤ U IM P ,

≤ Γk ≤ T.

(29)
(30)

3.2. Objective function
The objective function to be minimized in the planning problem is a combination of the 1-norm of the control signal, which is denote das JU , (which
gives an estimation of fuel consumption in case the control signal is thrust, see
Section 6) and a weighted 2-norm of the state, which is denoted as JX , both
taken over the planning horizon. Thus,
Jk = JU ,k + αJX ,k ,

(31)

where α is a positive constant that allows us to give a relative weight between
input cost and state error. JX ,k is computed as
JX ,k = XkT Qk Xk ,
10

(32)

for Qk > 0. Written in terms of the inputs and the starting point xk , JX ,k is
JX ,k = 2xTk FkT Qk Gk Uk + UkT GTk Qk Gk Uk ,

(33)

an expression in which the constant term xTk FkT Qk Fk xk , which does not play a
role in the planning optimization as it is constant for a given xk , is neglected.
The value of JU ,k does, however, depend on the control type.
3.2.1. PWM control inputs
For the case of PWM control inputs, using definition (6) it can be seen that
the objective function JU (k) is given by:
k+Np −1

JU ,k

=

X

=

UkT Λk

=

AJk Υk ,

 + T +

(Uj ) κj + (Uj− )T κ−
j

j=k

(34)

with AJk defined by blocks as

AJk = 

0

0

0

UkT


.

(35)

The times Γk where inputs start does not play a role in the cost function (only
their duration Λk ).
3.2.2. PAM control inputs
For the case of PAM control inputs, it can be seen that the objective function
JU (k) is given by:
k+Np −1

JU ,k

=

X

T kUj k1 = T kUk k1 .

(36)

j=k

3.2.3. Impulsive control inputs
For the case of impulsive control inputs, JU (k) is given by:
k+Np −1

JU ,k

=

X

kUj k1 = kUk k1 ,

(37)

j=k

where it should be noticed that, as in the PWM case, the location τk,i of the
impulses does not play a role in the cost function.
11

3.3. Planning optimization problem
Now, for each of the input types, one can formulate a planning optimization
problem starting from initial condition xk at time tk , with a planning horizon
of Np , as follows.
3.3.1. PWM control inputs
For PWM control inputs, the planning optimization problem is formulated
as

minΥk
s. t.

Ak Gk (Υk )Uk
AW Υ k

2xTk FkT Qk Gk (Υk )Uk + UkT GTk (Υk )Qk Gk (Υk )Uk + αAJk Υk ,
≤ bk − Ak Fk xk .

(38)

≤ bW .

(39)

Notice that Uk is known, and one has to compute the start and width of the
pulses, contained in Υk (start and duration of pulses), which enter nonlinearly
in the optimization problem. The dependence of Gk on Υk has been made
explicit.
3.3.2. PAM control inputs
For PAM control inputs, the planning optimization problem is formulated
as
minUk
s.t.

2xTk FkT Qk Gk Uk + UkT GTk Qk Gk Uk + αT kUk k1 ,

Ak Gk Uk

≤ bk − Ak Fk xk .

U P AM

≤ Uk ≤ U P AM .

12

(40)

3.3.3. Impulsive control inputs
For impulsive control inputs, the planning optimization problem is formulated as
minUk ,Γk
s. t.

2xTk FkT Qk Gk (Γk )Uk + UkT GTk (Γk )Qk Gk (Γk )Uk + αkUk k1 ,

Ak Gk (Γk )Uk

≤ bk − Ak Fk xk .

U IM P

≤ Uk ≤ U IM P ,

0

(41)

≤ Γk ≤ T.

The dependence of Gk on Γk (location of impulses) has been made explicit to
emphasize that the optimization problem is nonlinear.

4. PWM planning algorithm
In this section the subindex k is kept even though it does not play any role.
For a “pure” planning problem, it could be set to zero. However, k will be useful
when defining the MPC algorithm in Section 5.
Consider now the problem (38), given xk and Uk . Since the problem is
nonlinear, one needs to design an algorithm to solve it. The planning algorithm
is based on starting the problem using either the impulsive or PAM model. The
algorithm is composed of the following steps.
Step 1. Solve either the PAM optimization problem (40), or the impulsive
problem (41) with a fixed Γk , to provide an initial guess of the PWM
solution.
Step 2. The PAM or impulsive control inputs resulting from the optimization
algorithm in Step 1 are converted to a sequence of PWM inputs, denote
this initial sequence by Υ0k . Set i = 0.
Step 3. The trajectory of the system with the PWM inputs Υik is computed
analytically (if possible) or numerically by using equation (20). Denote
the trajectory by Xki .

13

Step 4. The system with PWM inputs is linearized around Xki , thus obtaining
a linear, explicit plant with respect to increments, denoted as ∆ik , in the
PWM inputs. Then a quadratic program can be posed and solved to find
the increments that improve the cost function.
Step 5. The resulting solution ∆ik is used to improve the approximation
towards the real solution, by setting Υi+1
= Υik + ∆ik . Increase i by one
k
and go back to Step 3. The process is iterated until the solution converges
or time is up.
Next, all the steps in the algorithm are described.
4.1. Step 1. Computation of PAM/impulsive control input
First, one has to choose if to find an initial guess using a PAM approach or
an impulsive approach. The PAM guess is more suitable if one expects wide
pulses, whereas the impulsive guess is best when the pulses are rather short.
If a PAM guess is chosen, it is computed from (40), setting U IM P = T Uk+
and U IM P = T Uk− , so that the solution can always be converted to PWM
following the procedure of Section 4.2. On the other hand, the impulsive guess
is computed from (41), setting U IM P = Uk+ and U IM P = Uk− . The impulsive
guess also requires to set the impulse location Γk to some pre-determined value,
so only the impulse magnitude (which appears linearly in (41)) is unknown.
Typical positions would be the middle of the interval (all entries of Γk equal to
T /2) or start of the interval (Γk = 0).
4.2. Step 2. Initial PWM solution: Adapting the PAM/impulsive solution
The PAM/impulsive solution from Section 4.1, U, is transformed to a PWM
sequence of inputs, as follows:
I
1. From U extract uA
j,i (or uj,i if the initial solution is of impulsive type) for

j = k, . . . , k + Np − 1 and i = 1, . . . , m. Also extract τj,i if the initial
solution is of impulsive type.

14

2. If the initial solution is of PAM type, set
 T uA
 T uA
A
 Wj,i
 − Wj,i
u
>
0,
, uA
± ,
j,i
j,i < 0,
uj,i
uj,i±
+
−
τj,i
(t) =
τj,i
(t) =


0,
uA
0,
uA
j,i ≤ 0,
j,i ≥ 0,
and if the initial solution is of impulsive type,
 uI

I
 Wj,i± , uA
 − uWj,i± ,
>
0,
j,i
uj,i
uj,i
+
−
τj,i =
τj,i =


0,
uIj,i ≤ 0,
0,

uA
j,i < 0,

(42)

(43)

uIj,i ≥ 0,

3. In the PAM case, the PWM input should be centered in the interval: κ+
j,i =
+
T −τj,i
,
2

κ−
j,i =

−
T −τj,i
.
2

In the impulsive case, the PWM input should be

centered around the chosen τj,i (corrected if necessary to avoid spillover),
i.e.

κ+
j,i

κ−
j,i


τ+


0,
τj,i − j,i

2 < 0,

+
τj,i
+
=
T − τj,i ,
τj,i + 2 > T,


+

 τ − τj,i , otherwise,
j,i
2

τ−


0,
τj,i − j,i

2 < 0,

−
τ
−
=
T − τj,i
,
τj,i + j,i
2 > T,


−

 τ − τj,i , otherwise,
j,i
2

(44)

(45)

+
−
−
4. From τj,i
, τj,i
, κ+
j,i , and κj,i , construct Γk and Λk and thus Υk .

The PWM signals Γk , Λk constructed by this method produce a moderately
similar (but not equal) output to the system driven by PAM or impulsive signals,
but as time advances the output might considerably differ. See Shieh et al.
(1996); Ieko et al. (1999); Bernelli-Zazzera et al. (1998) for more details and
other methods. In addition, the PWM results are not optimal (with respect to
the PWM signals) and they might not even verify the constraints. However, this
solution is only used as an initialization for the optimization algorithm proposed
next. Denote as Υ0k the found solution and set i = 0.
4.3. Step 3. Computation of trajectories under PWM inputs
For the current iteration i, apply (20) to compute the states of the system
Xki at all times, with PWM inputs Υik . The matrix Gk might be needed to
15

compute numerically if an explicit solution for the integrals (10) is not known
or possible.
4.4. Steps 4 and 5. Refined PWM solution: An optimization algorithm
To linearize (20) around inputs Υik , notice from (8) that
∂
B W (τk,i , κk,i )
∂τk,i k,i

=
=

∂
∂τk,i

tk +τk,i +κk,i

Z

Φ(tk+1 , s)Bi (s)ds,
tk +τk,i

Φ(tk+1 , tk + τk,i + κi )Bi (tk + τk,i + κk,i )
−Φ(tk+1 , tk + τk,i )Bi (tk + τk,i ),

(46)

and
∂
B W (τk,i , κk,i )
∂κk,i k,i

=
=

∂
∂κk,i

Z

tk +τk,i +κk,i

Φ(tk+1 , s)Bi (s)ds,
tk +τk,i

Φ(tk+1 , tk + τk,i + κi )Bi (tk + τk,i + κk,i ). (47)

−
Thus, (16) can be explicitly linearized around some given τ̄k+ , τ̄k− and κ̄+
k , κ̄k ,

reaching
xk+1 = Ak xk + Bk Uk + Bk∆τ δτk + Bk∆κ δκk ,

(48)

−
where Bk is computed with τ̄k+ , τ̄k− , κ̄+
k , and κ̄k , and define

Bk∆τ
δτ k

=

h

Bkδτ



τ̄k+ − τk+

= 

+

−

Bkδτ


τ̄k− − τk−

,

i

h
i
+
+
Bk∆κ = Bkδκ
,
Bkδκ


κ̄+ − κ+
k 
,
δκk =  k
−
κ̄−
−
κ
k
k
,

(49)
(50)

where the i-th entries of the Bkδ matrices in (50) are given, respectively, by
+

(Bkδτ̄ )i

=

+
+
+
W+
Φ(tk + T, tk + τ̄k,i
+ κ̄+
k,i )Bi (tk + τ̄k,i + κ̄k,i )uk,i
+
+
+
−Φ(tk + T, tk + τ̄k,i
)Bi (tk + τ̄k,i
)uW
k,i

−

(Bkδτ̄ )i

+

(Bkδκ̄ )i
−

(Bkδκ̄ )i

(51)

−
−
−
W−
= −Φ(tk + T, tk + τ̄k,i
+ κ̄−
k,i )Bi (tk + τ̄k,i + κ̄k,i )uk,i
−
−
−
+Φ(tk + T, tk + τ̄k,i
)Bi (tk + τ̄k,i
)uW
k,i

(52)

=

+
+
+
W+
Φ(tk + T, tk + τ̄k,i
+ κ̄+
k,i )Bi (tk + τ̄k,i + κ̄k,i )uk,i ,

(53)

=

−
−
−
W−
−Φ(tk + T, tk + τ̄k,i
+ κ̄−
k,i )Bi (tk + τ̄k,i + κ̄k,i )uk,i

(54)

16

Thus, defining stack vectors with the increments in the PWM variables at step
i as

∆Γik



= 






δτk
..
.



,




∆Λik = 


∆ik

=



∆Γik
∆Λik



,


δκk+Np −1

δτk+Np −1
and grouping all increments as




δκk
..
.





,

Bk∆ = 

Bk∆τ
Bk∆κ


,

and defining G∆
k as in (55), i.e.,


0

···

0



∆
Bk+1
..
.

···
..
.

0
..
.

∆
Ak+Np ,k+2 Bk+1

...

∆
Bk+N
p−1




,




Bk∆



 Ak+2,k+1 Bk∆

G∆
=
k
..


.

Ak+Np ,k+1 Bk∆

(55)

one can write
i
i
Xki ≈ Fk xk + Gk (Υik )Uk + G∆
k (Υk )∆k ,

(56)

The state constraints (85) become
i
i
i
Ak G∆
k (Υk )∆k ≤ bk − Ak Fk xk − Gk (Υk )Uk .

(57)

The constraints on ∆Γik and ∆Λik are as follows:
− AW ∆ik
∆max

≤ bW − AW Υik ,

(58)

≤ ∆ik ≤ ∆max ,

(59)

where the last constraint (59) is used to avoid large variations that might make
the linearization approximation to fail. All these constraints are might be summarized as
A∆ ∆ik ≤ b∆

17

(60)

Finally, the objective function can be rewritten in terms of ∆ik as Jki = Jki (Υik )+
Jk∆ (Υik , ∆ik ), where
Jki

=

2xTk FkT Qk Gk (Υik )Uk + UkT GTk (Υik )Qk Gk (Υik )Uk + αAJk Υik ,

Jk∆

=

i
i
i T ∆T
i
∆
i
i
2(xTk FkT + UkT GTk (Υik ))Qk G∆
k (Υk )∆k + (∆k ) Gk (Υk )Qk Gk (Υk )∆k

+αAJk ∆ik .

(61)

(62)

Noting that Jk∆ is quadratic in ∆ik , a quadratic optimization program with linear
restriction, formulated on the output increments, can be posed as follows:
min
i

Jk∆ (Υik , ∆ik )

s. t.:

i
i
Ak G∆
k (Υk )∆k

∆k

A∆ ∆ik

(63)
≤ bk − Ak Fk xk − Gk (Υik )Uk .,
≤ b∆ .

The solution ∆k is used to recompute new PWM inputs, Υi+1
= Υik + ∆ik .
k
Then the linearization process can be repeated around the new Υi+1
k , refining
the solution in each iteration.

5. Model Predictive Control with PWM inputs
In this section, building upon the trajectory planning algorithm of Section 4,
which is open-loop and has a finite time-horizon, a closed-loop algorithm is
developed based on the ideas of model predictive control (also known as receding
horizon control). Model predictive control closes the loop by simply re-planning
the maneuver at each time step, after applying just the set of control inputs
corresponding to the first time step, and keeps looking ahead Np time steps.
Thus, the algorithm starts at k = 0 and is repeated for each k. The re-planning
is done from the actual position at each time step, which seldom coincides with
the planned position due to disturbances, thus effectively closing the loop.
However, except at the start, it is not necessary to repeat all the steps of
Section 4. Since the new position should be close to the planned one, one can
apply the linearization scheme of the planning algorithm starting from the last
available linearization. The MPC algorithm is summarized next:
18

Step 1. At time step k = 0 and starting from x0 apply the Planning algorithm
of Section 4, obtaining a set of PWM inputs Υ0 that would optimize
the planning problem (38) for the next Np time steps, if there were no
disturbances.
Step 2. Apply impulses corresponding to the first time instant; save the rest
of impulses. Set k = 1
Step 3. One arrives at xk , which probably is not the intended value of the
state at time k but close. Thus re-planning is necessary.
Step 4. For re-planning, apply the planning algorithm of Section 4. However,
to avoid the initial step of having to use a PAM or impulsive model and
compute an initial guess, use instead as an initial guess the impulses of
Υk−1 that were not used (all of them except those corresponding to time
k − 1) and guess the remaining impulses (at the end) as zeros. In this way,
form an initial guess Υ0k .
Step 4. Apply the linearization algorithm of Section 4.4 using Υ0k as initial
guess to obtain, after iterating, a new set of impulses Υk . Apply the set
of impulses corresponding to time k. Save the rest of impulses.
Step 5. Repeat step 3.

6. Example application: Spacecraft Rendezvous
Rendezvous of spacecraft is the controlled close encounter of two (or more)
space vehicles. This work assumes just two vehicles, one of which is the target
vehicle (which is in a known orbit, and considered passive) and the other is the
chaser spacecraft, which begins from a known position and maneuvers until very
close to target. Only close range rendezvous Fehse (2003) is considered, which
starts at hundreds of meters and ends when the chaser is very close to target (a
few meters with speeds of centimeters per second).

19

There are numerous mathematical models for spacecraft rendezvous; which
one should be used depends on the parameters of the scenario. In Carter (1998)
a survey of numerous mathematical models for spacecraft rendezvous can be
found.
For instance, if the target is orbiting in a circular Keplerian orbit, the general
equations of the relative movement between an active chaser spacecraft close to
a passive target vehicle are linear time-invariant Hill-Clohessy-Wiltshire (HCW)
equations (introduced in Hill (1878) and Clohessy and Wiltshire (1960)). While
these equations are frequently used in the literature, it must be noted that,
in many situations, the HCW equations are not accurate. For instance, if the
target vehicle is moving in a Keplerian eccentric orbit (see Inalhan et al. (2002))
or if some orbital perturbations are taken into account (see for example Humi
and Carter (2008)). A more complex model, the Tschauner-Hempel model
(see Tschauner and Hempel (1965) or Carter (1998)) assumes that the target
vehicle is passive and moving along an elliptical orbit with semi-major axis a
and eccentricity e. The system equations are linear time-varying and cannot
be exactly integrated in time to obtain a discrete transition model; however, if
one substitutes the time t by the eccentric anomaly of the target orbit, E, it
is possible to obtain explicit expressions for the system evolution in the PWM,
impulsive, and PAM actuation cases. This will be the model considered in this
work. The model can be expressed in cartesian coordinates, but also in the
so-called relative orbital elements (see, e.g., Gaias et al. (2014) or Sinclair et al.
(2014)). The former has been chosen for simplicity.
Let us first establish some notation. Define the orbital mean motion n =
pµ
a3 , where µ is the gravitational parameter of the central body around which
the target spacecraft is orbiting.
Now, note that t and E are related in a one-to-one fashion by using Kepler’s
equation:
n(t − tp ) = E − e sin E,

(64)

where tp is the time at periapsis used as a starting point to measure the ec-

20

Figure 2: LVLH frame.

centric anomaly E. The time tp is chosen such that it is equal or less than the
starting time which is denoted as t0 (subtracting, if necessary, any number of
orbital periods). Kepler’s equation is not analytically invertible, but its inverse
can be found numerically with any desired degree of precision (see any Orbital
Mechanics reference, such as Wie (1998)). Denote its inverse by the function
K, i.e. E = K(t). Denote by E0 the true anomaly corresponding to t0 , this is,
E0 = K(t0 ). Then, Ek = K(tk ) = K(t0 + kT ), where T is the sampling time
(not to be confused with the orbital period). Call as rx,k , ry,k , and rz,k the position of the chaser in a local–vertical/local–horizontal (LVLH) frame of reference
fixed on the center of gravity of the target vehicle at time tk . In the (elliptical)
LVLH frame, x refers to the radial position, z to the out-of-plane position (in
the direction of the orbital angular momentum), and y is perpendicular to these
coordinates (not necessarily aligned with the target velocity given that its orbit
is not circular). The velocity and inputs of the chaser in the LVLH frame at
time tk are denoted, respectively, by vx,k , vy,k , and vz,k , and by ux,k , uy,k , and
uz,k .
If there is no actuation (i.e. ux,k = uy,k = uz,k = 0), the resulting transition

21

equation was obtained in a simple form in Yamanaka and Ankersen (2002) as
follows:
xk+1 = Φ(tk+1 , tk )xk

(65)

where
xk

=

T

[rx,k ry,k rz,k vx,k vy,k vz,k ] ,

(66)

and where
−1
Φ(tk+1 , tk ) = YK(tk+1 ) YK(t
,
k)

(67)

with Ytk being the fundamental matrix solution of the Tschauner-Hempel model,
which are expressed in Yamanaka and Ankersen (2002) as a function of true
anomaly θ. However there is a one-to-one relation between E and θ given by
r
θ
1+e
E
tan =
tan ,
(68)
2
1−e
2
which is exploited in the sequel. The explicit expression of the matrices1 is
found in (69) and (70),




YE = 




1 These

0

2/ρ − 3esJ

−c

1/ρ

0

−3ρJ

s(1 + 1/ρ)

0

0

c/ρ

0

0

s/ρ

0

0

α(−es − 3eρ2 Jc)

αρ2 s

0

αs(−1 − ρ )

αes

0

αρ(3esρJ − 3)

α(c + e + cρ2 )

0

0

0

−sα

0

0

(c + e)α

s

0

c(1 + 1/ρ)
0
αρ2 c
2

0





 ,(69)




expressions slightly differ from Yamanaka and Ankersen (2002) because the two

transformation matrices that appear in that paper have been pre-multiplied; also, the reference
axes are not the same.

22


YE−1





×







3J

=
2
(1 − e ) 


2

2

eρ2 (1 + ρ)

−e2 ρ2 s

0

e2 s/α

eρ/α

0

ρ2 (1 + ρ)

−eρ2 s

0

es/α

ρ/α

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

2

−s(ρ + 2ρ + e )

es (1 + ρ)

0

−es(1 + ρ)2

ρ2 (1 − ce) + e2 s2

0





1
+
 (1 − e2 )



0
c−2e/ρ
α
ec−2/ρ
α

s(ρ+1)
ρα
es(ρ+1)
− ρα

−

0



0





(70)




0

0

(c + e)(1 − e )

0

0

−s(1−e2 )
αρ

ρ2 (1 + ρ)

−esρ2

0

es
α

3ρ(c + e) − eρs2

−esc(1 + ρ) − e2 s

0

s
α

ρ
α
c(ρ+1)+e
αρ

0

0

0

s(1 − e2 )

0

0

2

where the following symbols are used (expressed in terms of E):
√
1 − e2
1 − e2 sin E
cos E − e
ρ =
,s=
,c =
,
1 − e cos E
1 − e cos E
1 − e cos E
E − Ê − e(sin E − sin Ê)
n
J = α
,α =
,
2
3/2
(1 − e )
(1 − e2 )3/2

0
c(1−e2 )
αρ

(71)
(72)

where Ê in (72) can be substituted by zero or any other desired reference value
of E. For instance, if when evaluating (67) one chooses Ê = Ek = K(tk ), then
−1
for YK(t
one gets J = 0 and the first matrix in (70) becomes zero.
k)

Using (67), one gets Ak in (16) explicitly, as well as Bk in the impulsive
case (explicitly defined in in terms of Φ(tk+1 , tk )). To obtain the Bk matrix in
the PWM and PAM cases, one needs to solve (8) or (12), respectively, which
involves an integral. Defining
Z

r2

Φ(r3 , s)Bi (s)ds,

bi (r1 , r2 , r3 ) =

(73)

r1

one has that, from (8),
W
Bk,i
(τk,i , κk,i ) = bi (tk + τk,i , tk + τk,i + κk,i , tk+1 ),

(74)

and, from (12),
A
Bk,i
= bi (tk , tk+1 , tk+1 ).

To compute the bi ’s, the following integral is needed
Z
−1
YK(t)
Ci+3 dt,
23

(75)

(76)

where Ci is a 6-element column vector of zeros with a value of one at row i.
For the computation, define the functions fi (t), for i = 1, 2, 3, as the indefinite
integrals of (76), in terms of eccentric anomaly
Z
1 − e cos E
fi (E) = YE−1 Ci+3
dE.
n

(77)

Once the fi ’s are computed, one finds the bi ’s as
bi (r1 , r2 , r3 ) = YK(r3 ) (fi (K(r2 )) − fi (K(r1 )))

(78)

Inserting the expression of (70) in (77) and integrating, one obtains


2
2
2
3
2(1 + 6e )S − 3e(2 + e )E + e Ch + e Ŝ/2

f1 =

−7/2
1 − e2
2α2









2e(8 − e2 )S − (4 + 7e2 − 2e4 )E + eCh + e2 Ŝ/2
0
3/2
−2eC 1 − e2

2 3/2
−2C 1 − e




,




(79)

0



f2 =


 
2 −3 

1−e


2
2α




2

C(4(1 + e ) − eC) − eEh − 3eE 2
eC(10 − 2e2 − eC) − Eh − 3E 2

0
√

2E 1 − e2

3/2

1 − e2 (4S − e(3E + Ŝ/2))






,





(80)

0



2 −5/2

f3 =

1−e
4α2









0
0
√
2 1 − e2 C(2 − eC)
0
0

4 e2 + 1 S − e(6E + Ŝ)





,




(81)

where S = sin(E), Ŝ = sin(2E), C = cos E, h = 6α(Ê − E − e sin(Ê)). Similar
expressions for the B matrices can be found in Ankersen (2010), however using
a slightly different definition of reference axes.
Note that, using these formulas, it is possible to express (16) explicitly for
all actuation types. This greatly speeds up the algorithms.
6.1. Constraints for the rendezvous problem
Besides the input constraints (which were given in Section 3.1.2), the inequality state constraints which were generically specified in Section 3.1.1 are,
24

in general, related to safety and sensing considerations (see e.g. Breger and How
(2008)). In this work, it is considered that during rendezvous the chaser vehicle
has to remain inside a line of sight (LOS) area. To simplify the constraint2 ,
in this work a 2-D LOS area is used as shown in Figure 3. This LOS region
is the intersection of a cone, given by the equations ry ≥ cLOS (rx − rx0 ) and
ry ≥ −cLOS (rx + rx0 ), and the region ry ≥ 0.

Figure 3: Line of Sight region.

The LOS constraint

0


ALOS =  cLOS

−cLOS

is ALOS xk ≤ bLOS , where


−1 0 0 0 0
0




−1 0 0 0 0  , bLOS =  cLOS rx0


−1 0 0 0 0
cLOS rx0




.


(82)

Using the compact formulation that was developed in Section 2, the con2 More

complicated constraints could be considered, see Gavilan et al. (2012) for examples

including a rotating LOS constraint.

25

straints equations for the state can be rewritten as:
Ac X ≤ bc ,
where Ac and bc are given by:

A


A

Ac = 




(83)





..

.

bLOS





 bLOS

 , bc = 
 ..

 .



bLOS
A





.




(84)

Then, using equation (20), one can reformulate the LOS constraints as constraints for the control signals, starting at time step tk , in the following way:
Ac Gk Uk ≤ bc − Ac Fk x0 .

(85)

7. Simulation Results
For simulations the following values have been used: Np = 50 as planning
horizon, T = 60 s, and ū = 10−1 N/kg. The target orbit has e = 0.7 and
perigee altitude hp = 500 km. Initial conditions were θ0 = 45o , r0 = [0.25 0.4 −
0.2]T km, v0 = [0.005 −0.005 −0.005]T km/s. The LOS constraint (see Vazquez
et al. (2011)) is defined by x0 = 0.001 km and CLOS = tan 30o . For the cost
function, α has been set to 103 and Qk as

Rk+1


..
Qk = 
.





,


(86)

Rk+Np
where Rk is defined as

Rk = h(k − ka ) 

Id3×3

Θ3×3

Θ3×3

Θ3×3


.

(87)

In (87), h is the step function, ka is the desired arrival time for rendezvous, and
Id3×3 , Θ3×3 are respectively the identity matrix and a matrix full of zeros, both
of order 3 by 3. The reason for choosing (87) is that it is desired to arrive at
26

0.5
Starting Point
PWM-MPC
Impulsive-MPC
Planning
LOS constraint

y [km]

0.4

0.3

Safe area
0.2

Constraint
violation
0.1

Restricted area
0
0

0.1

0.2

0.3

0.4

0.5

x [km]
Figure 4: System trajectories in the target orbital plane: open-loop PWM inputs computed
from impulsive solution (dashed), closed-loop Model Predictive Control with PWM inputs
using impulsive model (dot-dashed), and closed-loop Model Predictive Control with PWM
inputs using the PWM planning algorithm (solid).

the origin at time ka (and remain there) and at the same time minimize the
control effort.
In the simulations three algorithms were considered: first, an impulsive openloop trajectory planner, as described in Section 4.1. Next, closed-loop simulations using MPC but considering impulsive instead of PWM actuation in the
model (this algorithm is denoted as impulsive MPC). Finally, closed-loop simulations using MPC, based on the PWM algorithms as explained in Section 5.
The impulses produced by the first and second methods are subsequently transformed to PWM inputs using the algorithm of Section 4.2.
Compare first the algorithms without disturbances. The trajectories (projected on the target orbital plane) are shown in Fig. 4. The open-loop impulsive

27

0.5
Starting Point
PWM-MPC
Impulsive-MPC
Planning
LOS constraint

y [km]

0.4

0.3

Safe area
0.2

Constraint
violation

0.1

Restricted area
0
0

0.1

0.2

0.3

0.4

0.5

x [km]
Figure 5: System trajectories in the target orbital plane, with inexact orbit model: open-loop
PWM inputs computed with the planning algorithm (dashed), closed-loop Model Predictive
Control with PWM inputs using impulsive model (dot-dashed), and closed-loop Model Predictive Control with PWM inputs using the PWM planning algorithm (solid).

solution does not achieve rendezvous and drifts away, whereas the other solutions successfully reach the origin (the simulation is stopped when the chaser
vehicle was 5 meters or less away from the target). The impulsive MPC is able
to mostly compensate its imperfect thruster model. The PWM MPC algorithm
had a cost of 15.0 m/s and the impulsive MPC had a cost of 15.8 m/s. Thus,
while a basic MPC is able to rendezvous, the use of an imperfect model has
some fuel costs. In addition, the impulsive MPC does not satisfy the line-ofsight constraints for a period of time.
Next, Fig. 5 shows a simulation where the real orbit is different from the
reference orbit used in the model (the real eccentricity is e = 0.83, the real
perigeee altitude is hp = 525 km, and the real θ0 = 60o ). Both MPC algorithms

28

reach the origin (as in the previous scenario, the simulation is stopped when the
chaser vehicle was 5 meters or less away from the target). The impulsive MPC
again exits the line-of-sight region. The cost for the PWM MPC algorith was
15.3 m/s, whereas the impulsive MPC had a cost of 15.8 m/s.
Each iteration took less than half a second on a conventional computer, using
MATLAB and the Gurobi optimization package (see Gurobi Optimization, Inc.
(2014)). With a maximum number of iterations of 6, the computation time
remained well below the interval sampling time.

8. Concluding Remarks
This paper has presented a MPC algorithm that computes optimal PWM
inputs for LTV systems. The algorithm is based on an initial approximation
with either PAM or impulsive inputs, followed by iterative explicit linearization.
As an application, the problem of rendezvous in elliptical orbits has been considered. In particular, the algorithm might be particularly useful for satellites with
small specific thrust. The algorithm improves the fuel cost of an impulsive-only
MPC (with the impulses posteriorly transformed to PWM inputs), and is able
to satisfy safety constraints and handle disturbances such as imperfect knowledge of the target’s orbit. This algorithm would help avoiding having to include
a PWM approximation term in the “uncertainty budget” and therefore save
costs. However, inclusion of real-life constraints and more realistic simulations
are needed to validate the method.
Possible future lines of research include studying the convergence of the planning algorithm, guaranteeing constraint satisfaction by including an estimate of
linearization error in the model, or analyzing the stability guarantees of the
MPC design.

Acknowledgments
The authors acknowledge financial support of the Spanish Ministry of Science
and Innovation under grant DPI2008-05818.
29

References
References
Ankersen, F., 2010. Guidance, navigation, control and relative dynamics for
spacecraft proximity maneuvers. Ph.D. thesis.
Arzelier, D., Kara-Zaitri, M., Louembet, C., Delibasi, A., 2011. Using polynomial optimization to solve the fuel-optimal linear impulsive rendezvous problem. J. Guid. Contr. Dynam. 34, 1567–1572.
Arzelier, D., Louembet, C., Rondepierre, A., Kara-Zaitri, M., 2013. A new
mixed iterative algorithm to solve the fuel-optimal linear impulsive rendezvous
problem. J. Opt. Theor. Appl. 159, 210–230.
Asawa, S., Nagashio, T., Kida, T., 2006. Formation flight of spacecraft in earth
orbit via MPC. In: SICE-ICASE International Join Conference.
Bernelli-Zazzera, F., Mantegazza, P., Nurzia, V., 1998. Multi-pulse-width modulated control of linear systems. J. Guid. Contr. Dynam. 21 (1), 64–70.
Breger, L., How, J. P., 2008. Safe trajectories for autonomous rendezvous of
spacecraft. J. Guid. Contr. Dynam. 31 (5), 1478–1489.
Camacho, E., Bordons, C., 2004. Model Predictive Control, 2nd Edition.
Springer-Verlag, pp. 131–205.
Carter, T. E., 1998. State transition matrices for terminal rendezvous studies:
Brief survey and new example. J. Guid. Contr. Dynam. 21 (1), 148–155.
Clohessy, W. H., Wiltshire, R. S., 1960. Terminal guidance systems for satellite
rendezvous. J. Aerosp. Sc. 27 (9), 653–658.
D’Amico, S., Ardaens, J.-S., Gaias, G., Benninghoff, H., Schlepp, B., Jörgensen,
J. L., 2013. Noncooperative rendezvous using angles-only optical navigation:
System design and flight results. Journal of Guidance, Control, and Dynamics
36 (6), 1576–1595.
30

Deaconu, G., Louembet, C., Théron, A., 2014. Minimizing the effects of the
navigation uncertainties on the spacecraft rendezvous precision. Journal of
Guidance, Control, and Dynamics 37 (2), 695–700.
Deaconu, G., Louembet, C., Théron, A., 2015. Designing continuously constrained spacecraft relative trajectories for proximity operations. Journal of
Guidance, Control, and Dynamics 38 (7), 1208–1217.
Fehse, W., 2003. Automated Rendezvous and Docking of Spacecraft. Cambridge
University Press.
Gaias, G., D’Amico, S., Ardaens, J.-S., 2014. Angles-only navigation to a noncooperative satellite using relative orbital elements. Journal of Guidance, Control, and Dynamics 37 (2), 439–451.
Gavilan, F., Vazquez, R., Camacho, E. F., 2009. Robust model predictive control
for spacecraft rendezvous with online prediction of disturbance bound. In:
Proceedings of AGNFCS’09, Samara, Russia,.
Gavilan, F., Vazquez, R., Camacho, E. F., 2012. Chance-constrained model predictive control for spacecraft rendezvous with disturbance estimation. Contr.
Eng. Pract. 20 (2), 111–122.
Gurobi Optimization, Inc., 2014. Gurobi optimizer reference manual.
URL http://www.gurobi.com
Hartley, E. N., Trodden, P. A., Richards, A. G., Maciejowski, J. M., 2012.
Model predictive control system design and implementation for spacecraft
rendezvous. Control Engineering Practice 20 (7), 695 – 713.
Hill, G., 1878. Researches in lunar theory. American Journal of Mathematics
1 (3), 5–26, 129–147, 245–260.
Humi, M., Carter, T., 2008. Orbits and relative motion in the gravitational field
of an oblate body. J. Guid. Contr. Dynam. 31 (3), 522–532.

31

Ieko, T., Ochi, Y., Kanai, K., 1999. New design method for pulse-width modulation control systems via digital redesign. J. Guid. Contr. Dynam. 22 (1),
123–128.
Inalhan, G., Tillerson, M., How, J. P., 2002. Relative dynamics and control
of spacecraft formations in eccentric orbits. J. Guid. Contr. Dynam. 25 (1),
48–59.
Jewison, C., Erwin, R. S., Saenz-Otero, A., 2015. Model predictive control with
ellipsoid obstacle constraints for spacecraft rendezvous. In: ACNAAV 2015
IFAC workshop.
Kim, H. J., Shim, D. H., Sastry, S., 2002. Nonlinear model predictive tracking
control for rotorcraft-based unmaned aerial vehicles. In: Proceedings of ACC
2002.
Larsson, R., Berge, S., Bodin, P., Jönsson, U., 2006. Fuel efficient relative orbit control strategies for formation flying and rendezvous within prisma. In:
Proceedings of the 29th AAS guidance and control conference.
Leomanni, M., Rogers, E., Gabriel, S. B., 2014. Explicit model predictive control
approach for low-thrust spacecraft proximity operations. Journal of Guidance,
Control, and Dynamics 37 (6), 1780–1790.
Louembet, C., Arzelier, A., Deaconu, G., 2015. Robust rendezvous planning under maneuvering errors. Journal of Guidance, Control, and Dynamics 38 (1),
76–93.
Richards, A. G., How, J., 2003. Performance evaluation of rendezvous using
model predictive control. AIAA Paper 2003-5507.
Rossi, M., Lovera, M., 2002. A multirate predictive approach to orbit control of
small spacecraft. In: Proceedings of ACC 2002.
Rugh, W. J., 1996. Linear System Theory (2Nd Ed.). Prentice-Hall, Inc., Upper
Saddle River, NJ, USA.
32

Shieh, L.-S., Wang, W.-M., Sunkel, J., 1996. Design of PAM and PWM controllers for sampled-data interval systems. J Dyn Syst Meas Contr. 118 (4),
673–681.
Sinclair, A. J., Sherrill, R. E., Lovell, T. A., 2014. Calibration of linearized
solutions for satellite relative motion. Journal of Guidance, Control, and Dynamics 37 (4), 1362–1367.
Tschauner, J., Hempel, P., 1965. Rendevous zu einem in elliptischer bahn umlaufenden. Ziel. Acta Astronaut. II (2), 104–109.
Vazquez, R., Gavilan, F., Camacho, E. F., 2011. Trajectory planning for spacecraft rendezvous with on/off thrusters. In: Proc. of IFAC World Congress
2011.
Vazquez, R., Gavilan, F., Camacho, E. F., 2014. Trajectory planning for spacecraft rendezvous in elliptical orbits with On/Off thrusters. In: IFAC World
Congress, Cape Town.
Vazquez, R., Gavilan, F., Camacho, E. F., 2015. Model predictive control for
spacecraft rendezvous in elliptical orbits with On/Off thrusters. In: ACNAAV
2015 IFAC workshop.
Weiss, A., Kolmanovsky, I., Baldwin, M., Erwin, R. S., 2012. Model predictive
control of three dimensional spacecraft relative motion. In: American Control
Conference (ACC), 2012. IEEE, pp. 173–178.
Wie, B., 1998. Space vehicle dynamics and control. AIAA.
Yamanaka, K., Ankersen, F., 2002. New state transition matrix for relative
motion on an arbitrary elliptical orbit. J. Guid. Contr. Dynam. 25 (1), 60–66.

33

List of Figure Captions
1

PWM Variables. . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2

LVLH frame. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

3

Line of Sight region. . . . . . . . . . . . . . . . . . . . . . . . . .

25

4

System trajectories in the target orbital plane: open-loop PWM
inputs computed from impulsive solution (dashed), closed-loop
Model Predictive Control with PWM inputs using impulsive model
(dot-dashed), and closed-loop Model Predictive Control with PWM
inputs using the PWM planning algorithm (solid). . . . . . . . .

5

27

System trajectories in the target orbital plane, with inexact orbit model: open-loop PWM inputs computed with the planning
algorithm (dashed), closed-loop Model Predictive Control with
PWM inputs using impulsive model (dot-dashed), and closedloop Model Predictive Control with PWM inputs using the PWM
planning algorithm (solid). . . . . . . . . . . . . . . . . . . . . . .

34

28

