DISCRETE-TIME OPTIMAL ATTITUDE CONTROL OF
SPACECRAFT WITH MOMENTUM AND CONTROL
CONSTRAINTS

arXiv:1509.04521v2 [cs.SY] 29 May 2016

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR

Abstract. This article solves an optimal control problem arising in attitude
control of a spacecraft under state and control constraints. We first derive the
discrete-time attitude dynamics by employing discrete mechanics. The orientation transfer, with initial and final values of the orientation and momentum
and the time duration being specified, is posed as an energy optimal control
problem in discrete-time subject to momentum and control constraints. Using
variational analysis directly on the Lie group SO(3), we derive first order necessary conditions for optimality that leads to a constrained two point boundary
value problem. This two point boundary value problem is solved via a novel
multiple shooting technique that employs a root finding Newton algorithm.
Robustness of the multiple shooting technique is demonstrated through a few
representative numerical experiments.

1. Introduction
Typical space applications require reorienting a spacecraft or a satellite from
a given initial configuration to a given final configuration over a given duration
of time. Examples of such manoeuvres include positioning star sensors (attitude
estimation sensors) towards deep space, pointing a camera in the desired direction for imaging purposes, positioning solar panels for effective tracking of the sun
for optimal energy harvesting, etc [34]. Such orientation maneuvers are popularly
known as attitude maneuvers, and the class of attitude maneuvers that optimize a
certain performance objective are termed optimal attitude maneuvers. In most applications, attitude control problems are subjected to additional constraints due to
actuator saturation and/or forbidden regions of the state space. Control constraints
are omnipresent — for instance, a momentum delivering device fitted on-board a
spacecraft, e.g., a thruster, or a reaction wheel, has limitations in terms of maximum torque produced and maximum momentum delivered. In case a spacecraft
is fitted with a flexible structure like solar panels, it must inevitably have a restriction of its maximum momentum for the protection of such equipments; such
restrictions translate to state constraints in an optimal control problem. The class
of constrained control problems involving orientation manoeuvres provides a rich
supply of interesting and complex control problems related to space applications.
Few computationally tractable solutions to such problems are known today, and
in this article we propose a new and promising technique to solve a class of such
constrained optimal control problems.
In general, computationally tractable solutions to optimal control problems subject to nonlinear dynamics are difficult to arrive at [35]. Constrained optimal attitude control problems are vastly more challenging than standard optimal control
Key words and phrases. Constrained optimal control, Variational methods, multiple shooting;
(spacecraft) attitude control; state constraints and discrete mechanics.
1

2

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR

problems on Euclidean spaces because (a) the attitude kinematic equations evolve
on the manifold SO(3), and classical techniques developed specifically for dynamics
on Euclidean spaces, therefore, no longer directly apply, and (b) the presence of
state and control constraints in the optimal control problem typically leads to two
point boundary value problems (which represents first order necessary optimality
conditions,) subject to inequality constraints, and these problems cannot be solved
using conventional indirect techniques such as shooting methods [35]. In this article
we provide an algorithmic solution to one such class of constrained optimal attitude
control problems subject to state and control constraints, where the attitude kinematics in discrete time are treated directly on the Lie group SO(3) and are derived
using discrete mechanics [26].
Constrained optimal attitude control problems treated in the literature typically
fall into one of two categories: one that considers only control constraints and the
other that considers state constraints in addition to control constraints. Attitude
control problems with control constraints have been discussed under the framework
of optimal control with various performance indices involving minimization of time,
fuel, and energy. A considerable body of work on the time-optimal attitude control
problems exist [34, 29, 36, 10]. Fuel or energy optimal maneuvers are also of great
interest since on-board energy sources are limited and precious, and there is therefore a natural necessity for using them in an optimal way [3]. Early works on energy
optimal control [2, 17] and fuel optimal control [9, 18] tackled the continuous-time
optimal attitude control problem, arriving at first order necessary conditions for
optimality via the Pontryagin maximum principle. The resulting boundary value
problems were typically solved using shooting or neighborhood extremal methods,
both of which suffer from high sensitivity to initial data. In recent years, attitude
control problems with joint state and control constraints have been attacked from
various directions. Attitude manoeuvres with joint state and control constraints,
considering a linearized dynamics model is addressed in [37, 13]. In these works
the problem is recast as a nonlinear optimization problem, and solved using direct
multiple shooting methods or sequential quadratic programming [6]. Constrained
attitude control problems with attitude kinematics represented in quaternion form
was explored in [31, 24] using pseudospectral methods, and in [38] using particle
swarm optimization. Representation of the attitude kinematics in quaternion form
has a serious disadvantage of non-uniqueness; consequently, boundary conditions
defined for attitude kinematics do not have a unique representation in quaternions,
which presents a problem during computations. Another commonly used attitude
representation is Euler angles, which suffers from the defect of singularity [33]. In
order to avoid issues such as singularity and non-uniqueness, geometric techniques
have been developed to solve constrained attitude control problems [26, 5, 21]. One
such approach to attitude control problems in continuous time with control constraints is discussed in [32]. Projection operators were used there to find a search
direction on the underlying Lie group, while employing Newton’s method to solve
the continuous time two point boundary value problem obtained via the Pontryagin maximum principle on manifolds [1, p. 165]. A different geometric technique
to attitude control problem with state and control constraints is addressed in [14]
where the authors attempted to handle state inequality constraints using penalty
functions. As is well known, penalty functions do not enforce the state inequality
constraints, but add penalties to the cost functions if the constraints are violated.
This constrained optimal control problem leads to a two point boundary value problem that represents first order necessary optimality condition. The boundary value
problem is then solved by employing indirect single shooting method. None of the
these works discussed the issue of digital implementation, a key aspect of which is

ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS

3

the process of discretization of the boundary value problems in time: Euler’s step
and its derivatives are insufficient since these do not account for nor respect the
underlying manifold structure of SO(3).
In this article we tackle the issue of discretization up front by deriving a discrete time model via discrete mechanics [21, 19], and then employing variational
analysis to arrive at first order necessary conditions for optimality. Discrete time
models obtained via discrete mechanics are more accurate than other standard discretization schemes such as Euler’s step because they preserve certain invariance
properties like kinetic energy, momentum, etc, of the system, and the computations
can be done directly on the manifold SO(3), (because this discretization respects
the manifold SO(3),) thereby eliminating the problems associated with parametric
representations of attitude. The presence of state inequality constraints in our attitude control problem makes it challenging because the resulting boundary value
problems obtained using variational analysis are subject to inequality constraints on
the states. Such constrained boundary value problems cannot in general be solved
using classical multiple shooting methods precisely because of those inequality constraints. A non-classical multiple shooting technique has recently been proposed
in [11] for solving constrained boundary value problems arising in optimal control
problems with state constraints. There the complementary slackness conditions
arising due to inequality constraints on the states are represented in the form of
equality constraints using the Fischer-Burmeister function [11, Equation (2.7)] —
a technique that can be highly inefficient in terms of computation because all the
inequality constraints are considered at each iteration irrespective of them being
active or not. In contrast, we propose a multiple shooting algorithm that deals with
state inequality constraints in such a way that the dimension of the boundary value
problem remains the same even when the state inequality constraints are active,
making it more efficient in terms of time and memory complexity.
We reiterate that the thrust of our contribution is towards computational tractability of optimal attitude manoeuvres. The advantages of our technique are five
fold: First, the discrete time model derived via discrete mechanics avoids the need
for discretization at later stages, and eliminates issues associated with parametric
representations of the attitude kinematic. Second, an indirect multiple shooting
method, which provides more accurate solution to the problem unlike direct techniques, is employed to solve the optimal control problem. Integration steps of
the discrete kinematic equations represented by rotation matrices are calculated
at each discrete instant, which provide exact attitude trajectories for manoeuvres
unlike other standard schemes [35]. Third, the multiple shooting method that we
employ here is more robust to initial guesses as compared to indirect single shooting methods that are comparable in terms of accuracy [35]. Fourth, the proposed
algorithm, being a multiple shooting method, can be implemented on a parallel
architecture for fast computation. Fifth, the discrete time model obtained using
discrete mechanics results in a boundary value problem that can be reduced to
difference equations with momentum and co-momentum dynamics only. The dimension of this reduced difference equation model is half of the original one, which
further contributes to savings in terms of time and memory. Moreover, scaling of
the co-state variables such that the co-states are invariant to step length selected
for a given manoeuvre, improves the radius of convergence of the multiple shooting
algorithm.
This article unfolds as follows: In §2 we employ discrete mechanics to derive a
discrete time model of the attitude dynamics. In §3 the energy optimal control

4

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR

problem is posed as a discrete time optimal control problem, and first order necessary conditions are obtained using variational analysis. Then we discuss scaling of
the variables and reduction of the dynamics to momentum and co-momentum variables in §4. §5 contains an introduction to multiple shooting methods, and provides
a solution to the system of difference equations represented in momentum variables
presented in §4. §6 provides the numerical experiments for large angle manoeuvres
with momentum and control constraints. We conclude in §7 with a brief discussion
of future work. The proofs of our results are presented in a consolidated fashion in
the Appendices.
2. Derivation of Discrete time model via discrete mechanics
This section contains the modeling of the attitude dynamics of the spacecraft
using discrete mechanics [26]. In most optimal control problems involving mechanical systems, some sort of discretization is performed in order to employ numerical
techniques. In the case of discrete mechanics, the variational description is directly
discretized, and discrete time equations are obtained. This approach is advantageous in comparison to usual discretizations of continuous time models because it
preserves certain invariants of the system such as momentum and energy.
For ease of understanding, a quick introduction to discrete mechanics is given
here, followed by discrete time modeling of the attitude dynamics of the spacecraft.
2.1. Introduction. Consider a mechanical system with the configuration space Q
as a smooth manifold. Then the velocity vectors lie on the tangent bundle T Q of the
manifold Q and the Lagrangian for the system can be defined as L : T Q → R [25].
In discrete mechanics, the velocity phase space T Q is replaced by Q × Q which is
locally isomorphic to T Q. Let us consider an integral curve q(t) in the configuration
space such that q(0) = q0 and q(h) = q1 , where h represents the integration step.
Then, the discrete Lagrangian Ld : Q × Q → R, which is an approximation of the
action integral along the integral curve segment between q0 and q1 , can be defined
as [27]
Z h
(2.1)
Ld (q0 , q1 ) ≈
L (q(t), q̇(t)) dt.
0

Pick h > 0, (this quantity plays the role of step length) and consider a grid of the
time domain T = N h as {tk = kh|k = 0, 1, . . . , N } and the corresponding discrete
path space Pd (Q) := qd : {tk }N
k=0 → Q . The discrete trajectory qd ∈ Pd (Q) is
such that qd (tk ) = qk . Now, defining the discrete action sum Gd as
Gd (qd ) :=

N
−1
X

Ld (qk , qk+1 ) .

k=0

Assuming qd (0) = q0 and qd (tN ) = qN fixed, define the variations δ(qd ) as δ(qd (tk )) =
δ(qk ) ∈ Tqk Q which vanishes at the end points (i.e., δ(q0 ) = δ(qN ) = 0). A discrete
path qd is a stationary point of the discrete action sum Gd if Dqd Gd (qd )δ(qd ) = 0
for all δ(qd ) [16, p. 21]. This is equivalent to saying that the points {qk } of the path
qd satisfies the discrete Euler-Lagrange equations, i.e.,
(2.2)

D2 Ld (qk−1 , qk ) + D1 Ld (qk , qk+1 ) = 0

for all k = 1, 2, . . . , N − 1,

where Di is the derivative of the function with respect to the ith argument.
Notice that (2.2) involves qk−1 , qk and qk+1 at kth instant of time, which means
that the difference equations obtained are of second order. To arrive at the discrete
time model with first order difference equations, one needs to use the discrete time

ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS

5

analogue of the Hamiltonian formulation. To this end, recall that the continuous
time Legendre transform is a map FL from the Lagrangian state space T Q to the
Hamiltonian phase space T ∗ Q. Similarly, the discrete time Legendre transforms
F+ Ld , F− Ld : Q × Q → T ∗ Q [26] can be defined as
F+ Ld (qk , qk+1 ) 7→ (qk+1 , pk+1 ) = (qk , D2 Ld (qk , qk+1 )) ,
F− Ld (qk , qk+1 ) 7→ (qk , pk ) = (qk , −D1 Ld (qk , qk+1 )) ,
which are maps from the discrete Lagrangian state space Q × Q to the discrete
Hamiltonian phase space T ∗ Q. The map F+ Ld is the forward discrete Legendre
transform which relates (qk , qk+1 ) to Tq∗k+1 Q, and the map F− Ld is the backward
discrete Legendre transform which relates (qk , qk+1 ) to Tq∗k Q. Let the discrete
Lagrangian map FLd : Q × Q → Q × Q be defined as FLd (qk , qk+1 ) = (qk+1 , qk+2 );
it defines the evolution of the dynamics on the discrete state space. Then the
corresponding discrete Hamiltonian map F̃Ld : T ∗ Q → T ∗ Q can be defined in the
following equivalent ways [26]:
−1
−1
;
= F+ L d ◦ F− L d
F̃Ld := F± Ld ◦ FLd ◦ F± Ld

this is clear from the commuting diagram shown in Figure 1. The discrete Hamil(qk , qk+1 )

(qk+1 , qk+2 )

F+ Ld F− Ld

F− Ld

(qk , pk )

FLd

F̃Ld

(qk+1 , pk+1 )

F+ Ld
F̃Ld

(qk+2 , pk+2 )

Figure 1. Flow of the discrete Lagrangian and Hamiltonian map
tonian map is defined in coordinates as follows:
(2.3)
(qk , pk ) 7→ F̃Ld (qk , pk ) := (qk+1 , pk+1 )

where

(

pk = −D1 Ld (qk , qk+1 ) ,
pk+1 = D2 Ld (qk , qk+1 ) .

2.2. Attitude dynamics in discrete time. We now apply the ideas introduced
in §2.1 to obtain the discrete equations of the attitude dynamics of a spacecraft.
First, the Lagrangian in continuous time is described, and then an approximation of the continuous time Lagrangian is taken to define the discrete Lagrangian
(2.1). Thereafter, the discrete time attitude dynamics is obtained using discrete
Hamiltonian formulation (2.3).
Consider a rigid body with a point, typically chosen to be center of mass, fixed
on it. In order to define the orientation of a rigid body, two coordinate systems are
considered with the origin at that fixed point. One frame fixed to the rigid body
is known as the body frame, and the other is a frame fixed in space, known as the
spatial frame. Let X be the position of the mass element in the body frame. Then
the position of the mass element in reference frame x is related to the body frame
coordinates X by the rotation matrix R(t) ∈ SO(3) as x(t) = R(t)X. Let B be the
region occupied by the body in its reference frame. Let ρ(X) be the density of the
rigid body in the body coordinates at point X. Then the kinetic energy of the rigid
body is [16, p. 243]:
Z
1
ρ(X) kẋk2 d3 X,
K=
2 B

6

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR

which can be rewritten, in view of the left-invariance of the kinetic energy [7, p.
275], as
Z
Z
2
2
1
1
3
K=
ρ(X) ṘX d X =
ρ(X) R−1 ṘX d3 X.
2 B
2 B
We know that the spatial angular velocity vector Ω can be represented in terms of
the body angular velocity vector ω as Ω = R−1 ω. Then
(2.4)

Ω × X = R−1 ω × R−1 x = R−1 (ω × x) = R−1 ẋ = R−1 ṘX.

b ∈ so(3) be a vector space homeomorphism. Then, from (2.4) we
Let R3 ∋ Ω 7→ Ω
b = R−1 Ṙ ∈ so(3), the Lie algebra of SO(3). So,
arrive at the kinematic relation Ω
the kinetic energy can be represented in terms of the spatial frame angular velocity
as
Z




1
b
b ⊤ d3 X = 1 tr ΩJ
b dΩ
b⊤ ,
ρ(X) tr (ΩX)(
ΩX)
K=
2 B
2
where
Z
1
Jd =
ρ(X)XX ⊤d3 X.
2 B
R
b ⊤ Xd
b 3 X is related to Jd by
The body moment of inertia matrix J := 21 B ρ(X)X
the following equation [23]:
J = tr (Jd ) I3×3 − Jd .

If the dissipative and potential forces are absent, then the Lagrangian L : T SO(3) →
R for the system is given by [16, p. 245]
1  b b ⊤
L(R, Ω) := K = tr ΩJ
.
dΩ
2
b given above, we know that Ω
b = R⊤ Ṙ. So, the
By the kinematic relation Ṙ = RΩ
Lagrangian can be written as

1 
L(R, Ṙ) = tr R⊤ ṘJd Ṙ⊤ R .
(2.5)
2
We now proceed to discretize the Lagrangian (2.5). Considering discrete time instants tk = kh for k = 0, 1, . . . , such that R(tk ) = Rk and the approximation
Ṙ(tk ) ≈ (Rk+1h−Rk ) for t ∈ [tk , tk+1 ] [23], the discrete Lagrangian Ld : SO(3) ×
SO(3) → R is defined as:


(Rk+1 − Rk )
Ld (Rk , Rk+1 ) ≈ hL Rk ,
h

 ⊤
Rk (Rk+1 − Rk ) (Rk+1 − Rk )⊤ Rk
h
Jd
= tr
2
h
h

1
=
tr (Rk⊤ Rk+1 − I3×3 )Jd (Rk⊤ Rk+1 − I3×3 )
2h
1
= tr ((I3×3 − Fk ) Jd ) ,
(2.6)
h
where Fk = Rk⊤ Rk+1 . Note that under the discretization technique employed, the
discrete Lagrangian, like its continuous counterpart, is invariant under the action of
the SO(3) group. This property will be useful later when momentum equations will
be derived and the rotation sequence constructed based on the momentum history.
Our objective is to come up with first order difference equations describing the
attitude dynamics of the spacecraft. By the left trivialization of the cotangent
bundle of a Lie group, T ∗ SO(3) can be represented as SO(3) × so(3)∗ , where so(3)∗
denotes the dual of the Lie algebra so(3) [1, p. 254]. We now proceed to find the

ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS

7

discrete time Hamiltonian map (2.3) F̃Ld : SO(3) × so(3)∗ → SO(3) × so(3)∗ such
that

b

Πk = −D1 Ld (Rk , Rk+1 ) ,

 

b
b
b k+1 = D2 Ld (Rk , Rk+1 ) ,
F̃Ld Rk , Πk = Rk+1 , Πk+1 , where Π


Πk ∈ R3 .

b k and Π
b k+1 , the variations in Rk are defined in terms of ηbk ∈ so(3),
In order to find Π
and the expressions D1 Ld (Rk , Rk+1 ) and D2 Ld (Rk , Rk+1 ) are evaluated using the
duality product on so(3) [25, p. 290].
For a given ǫ ∈ R and ηk ∈ R3 , the variation in Rk can be defined as
(2.7)

δRk :=

d
dǫ

ǫ=0

Rk e(ǫbηk ) = Rk ηbk , where ηbk ∈ so(3);

b k ∈ so(3)∗ and ηbk ∈ so(3) is defined as
then the duality product of Π
E
D
 
1  b ⊤
1
b k , ηbk = − d
tr Πk ηbk =: Π
Ld (Rkǫ , Rk+1 ) = tr δRk⊤ Rk+1 Jd
2
dǫ ǫ=0
h

 1
1
= tr ηbk⊤ Rk⊤ Rk+1 Jd = tr Fk Jd ηbk⊤ .
h
h
Hence,


 
 1

b k − 1 Fk Jd ηb⊤ 
tr  Π
bk ∈ so(3),
k  = 0 for all η
 2
h

|
{z
}
Ck

which means that Ck is a symmetric matrix. Its skew-symmetric part is, therefore
zero, and this leads to
c k = Fk Jd − Jd F ⊤ .
hΠ
k

(2.8)

b k+1 ∈ so(3)∗ and ηbk+1 ∈ so(3) gives,
Similarly, the duality product of Π
E

D
 
1
1 b
⊤
ǫ
b k+1 , ηbk+1 = d
=: Π
tr Πk+1 ηbk+1
Ld (Rk , Rk+1
) = tr − Rk⊤ δRk+1 Jd
2
dǫ ǫ=0
h
 1

1
⊤
= − tr Rk⊤ Rk+1 ηbk+1 Jd = tr Jd Fk ηbk+1
,
h
h
and




 1
1

⊤ 
b
tr  Πk+1 − Jd Fk ηbk+1
 = 0 for all ηbk+1 ∈ so(3)
h

 2
|
{z
}
Dk

which means that Dk is a symmetric matrix with its skew-symmetric part equal to
zero. Therefore,
⊤
⊤
b k Fk = F\
b k+1 = Jd Fk − Fk Jd = Fk⊤ Π
Π
k Πk ,
h
leading to the following update equation for the momentum:

(2.9)

(2.10)

Πk+1 = Fk⊤ Πk .

In the presence of control, (2.10) modifies to
(2.11)

Πk+1 = Fk⊤ Πk + huk ,

8

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR

where uk is the control input at kth instant of time. The rigid body equations in
discrete time are finally obtained as:


Rk+1 = Rk Fk ,
Rigid Body
Πk+1 = Fk⊤ Πk + huk ,
Dynamics 
d
hΠk = Fk Jd − Jd Fk⊤ .
3. Optimal control of discrete time attitude dynamics

We state the optimal control problem arising in executing energy optimal attitude manoeuvres of a spacecraft. The spacecraft is assumed to have three actuators,
aligned along the three principal moment of inertia axes. Each actuator has its individual saturation limits. The objective is to find the energy optimal control profile
for orienting the spacecraft from an initial configuration to a desired configuration
in a given duration of time, while obeying pre-specified momentum bounds. First
we pose this requirement as an optimal control problem, and then derive the first
order necessary optimality conditions using variational analysis. Later §5.2 the
boundary value problem obtained as the first order necessary conditions will be
solved using a novel multiple shooting method.
3.1. Problem description. Our objective is to find the energy optimal control
law to manoeuvre a spacecraft from the initial configuration (Ri , Πi ) to the final
configuration (Rf , Πf ) in N discrete time steps satisfying the following constraints:
(1) uik ≤ ci
(2) Πik ≤ di

k = 0, 1, . . . , N − 1, and i = 1, 2, 3,
k = 1, 2, . . . , N − 1, and i = 1, 2, 3.

This problem can be posed as an optimal control problem in discrete time as follows:
(3.1)

minimize J :=
{uk }

subject to
(3.2)

N
−1
X
k=0



Rk+1 = Rk Fk
system of equations Πk+1 = Fk⊤ Πk + huk

d
hΠk = Fk Jd − Jd Fk⊤

1
kuk k22
2

with

(3.3)
boundary conditions (R0 , Π0 ) = (Ri , Πi ), (RN , ΠN ) = (Rf , Πf ), and
(3.4)
constraints

(

2
2
uik ≤ ci
for all k = 0, 1, . . . , N − 1, and i = 1, 2, 3,

2
2
Πik ≤ di
for all k = 1, 2, . . . , N − 1, and i = 1, 2, 3.

Note that the optimal control problem (3.1) has both control and state inequality
constraints. While the individual control inputs are constrained in magnitude, the
performance measure reflects a 2-norm on the control action at each stage.
3.2. Necessary optimality conditions. We represent the variations of Fk in
terms of variations in Πk and then the first order necessary conditions are derived.
• Representation of the variations:
Using (2.7) the variations for the matrix Rk⊤ Rk+1 is defined as

(3.5)
δ Rk⊤ Rk+1 = δRk⊤ Rk+1 + Rk⊤ δRk+1 = −ηbk Rk + Rk ηbk+1 .

ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS

(3.6)

9

[
⊤x = F ⊤x
Using the property F
bF , (3.5) simplifies to

∧
δ Rk⊤ Rk+1 = Fk −Fk⊤ ηk + ηk+1 ,

where ( · )∧ : R3 → so(3). Similarly, for a given ξk ∈ R3 , we define the
variation in Fk ∈ SO(3) as
δFk = Fk ξbk .

(3.7)

dk = Fk Jd − Jd F ⊤ in (3.2) gives the relation beThe implicit equation hΠ
k
tween momentum and change in orientation at kth time instant. So, the
relation between the variations in momentum δΠk and δFk can be obtained
from the implicit equation as [21]
⊤[
[k = δFk Jd − Jd δFk⊤ = F
[
hδΠ
k ξk Fk Jd + Jd Fk Fk ξk ,

which can be further simplified using x
bA + A⊤ x
b = ({tr[A]I3×3 − A} x)∧ to
[k = ((tr (Fk Jd ) I3×3 − Fk Jd ) Fk ξk )∧ .
hδΠ

Lemma 3.1. The matrix (tr (Fk Jd ) I3×3 − Fk Jd ) is invertible if
 s

kξk k
2d3 + d2 − d1
<
,
cos
2
2(d3 + d2 )
b

where Fk = eξk ,
d2 ≤ d3 .

ξk ∈ R3 , and Jd = diag(d1 , d2 , d3 ) such that 0 < d1 ≤

We present a proof of Lemma 3.1 in Appendix A. Armed with Lemma
3.1, we represent the vector ξk in terms of the variations in momentum i.e.
δΠk [21] as
ξk = Bk δΠk ,

(3.8)

−1

where Bk = hFk⊤ (tr (Fk Jd ) I3×3 − Fk Jd ) .
• Necessary optimality conditions:
Let χk ∈ R3 and λk ∈ R3 be the Lagrange multipliers corresponding to the
equality constraints Rk+1 − Rk Fk = 0 and Πk+1 − Fk⊤ Πk − huk = 0 respectively. Similarly let 0 ≤ αik ∈ R and 0 ≤ βki ∈ R be the Lagrange multipliers
2
2
2
corresponding to the inequality constraints uik ≤ ci
and Πik ≤
2
di . Let us justify why the Lagrange multiplier χk ∈ R3 is chosen corresponding to the rotational kinematics Rk+1 − Rk Fk = 0. Rotational
kinematics can be rewritten as Rk⊤ Rk+1 − Fk = 0, where Fk can be identified by its skew symmetric part, which in turn can identified by a vector in
R3 [22].
Claim 3.2. Consider the equality Rk⊤ Rk+1 = Fk . If we assume that the
step length h is small enough such that the relative orientation Rk⊤ Rk+1
between two adjacent time instances tk and tk+1 is less than π2 , i.e.,
kξk k , kζk k <

π
b
b
where Rk⊤ Rk+1 = eζk , Fk = eξk and ζk , ξk ∈ R3 ,
2

then the aforementioned equality is satisfied if and only if the skew symmetric parts of both sides are identical.

10

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR

The augmented performance index can be defined as
minimize Ja :=
|uik |≤ci

N
−1
X
k=0

1
huk , uk i + λk , −Πk+1 + Fk⊤ Πk + huk
2



∨ 1
∨
1
⊤
Fk − Fk⊤ −
Rk⊤ Rk+1 − Rk+1
Rk
+ χk ,
2
2

+

N
−1
X
1
1
hαk , uk ⊙ uk − c ⊙ ci +
hβl , Πl ⊙ Πl − d ⊙ di ,
2
2
l=1

where
c = c1

c2

c3

⊤

,

d = d1

d2

d3

⊤

,

x ⊙ y = x1 y 1

x2 y 2

x3 y 3

⊤

.

Using (3.6), (3.7) and (3.8), the infinitesimal variation of the augmented
performance index is defined as
δJa :=

N
−1
X

{hδuk , uk i +

k=0

Dχ

k

2

, δ Fk − Fk⊤

∨

⊤
− δ Rk⊤ Rk+1 − Rk+1
Rk

∨ E

+ λk , −δΠk+1 + δFk⊤ Πk + Fk⊤ δΠk + hδuk + hδuk , αk ⊙ uk i}

+

N
−1
X

hδΠl , βl ⊙ Πl i ,

l=1

where

αik , βli ≥ 0

(3.9)

for all i = 1, 2, 3,

together with the complementary slackness conditions


(3.10)
βli (Πil )2 − (di )2 = 0 and αik (uik )2 − (ci )2 = 0.

∧

Employing the property that x
bA + A⊤ x
b = {(tr[A]I3×3 − A) x} and that
the variation on the boundary is zero, i.e., η0 = ηN = 0, δΠ0 = δΠN = 0, we
rearrange the terms of the expression δJa and get the following expression
after standard algebraic manipulations:

δJa =

N
−1
X

hδuk , uk + hλk + αk ⊙ uk i

k=0

+

N
−1 D
X
k=1

+

(3.11)

+

 E

⊤
δΠk , −λk−1 + Fk − Bk⊤ F\
k Πk λk

N
−1 D
X
k=1

E
ηk
, Fk (tr (Fk ) I3×3 − Fk ) χk − (tr (Fk−1 ) I3×3 − Fk−1 ) χk−1
2

N
−1 
X
k=1


1
δΠk , βk ⊙ Πk + Bk⊤ (tr (Fk ) I3×3 − Fk ) χk .
2

By first order necessary condition of optimality, δJa = 0 along all possible
variations δuk , ηk , δΠk , and we obtain the co-state equations from (3.11)
as:
(3.12)

(
(Fk ) I3×3 − Fk ) χk = 0,
Co-state (tr (Fk−1 ) I3×3 − Fk−1 ) χk−1 − Fk (tr

1 ⊤
⊤Π
Dynamics βk ⊙ Πk − λk−1 + Fk − Bk⊤ F\
k λk + 2 Bk (tr (Fk ) I3×3 − Fk ) χk = 0,
k

ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS

11

and the optimality condition for the control as
ůk + hλk + αk ⊙ ůk = 0.

(3.13)

From 
(3.10) we know that if ůik lies in the interior of the feasible region
i
C := ui ∈ R (ui )2 < (ci )2 , then αik = 0. So, by (3.13) we have hλik =
−ůk . On the other hand if ůik = ci then by (3.13) and (3.9) we have hλik =
−ci (1 + αik ) ≤ −ci . Similarly for ůik = −ci , we have hλik = ci (1 + αik ) ≥ ci .
Hence, the optimal control can be written in a compact form as

ůik = − min ci , hλik sgn(λik ).
(3.14)
We represent the state and co-state dynamics in terms on momentum and costate corresponding to momentum variables. This technique will reduce the model
of the system and hence the algorithm will perform better in terms of memory and
time requirement.
4. Scaling and model reduction
The orientation boundary constraints (R0 , RN ) = (Ri , Rf ) can be represented
in terms of Fk , which in turn can be computed for a given Πk using the implicit
form
dk = Fk Jd − Jd F ⊤ .
(4.1)
hΠ
k

We represent orientation constraints first in terms of Fk and then we discuss a
technique for computing Fk for a given value of Πk .

4.1. Representing orientation constraints in terms of momentum. We represent the boundary constraints on orientation (R0 , RN ) = (Ri , Rf ) (3.3) in terms
of Fk . Then using (3.2) we reconstruct the orientation Rk variables, we see that
Rf can be represented in terms of Fk for k = 1, 2, . . . , N − 1, as
Rf = RN = Ri F0 F1 F2 . . . FN −1 .
The boundary constraints of orientation can be rewritten as:
Ri⊤ Rf = F0 F1 F2 . . . FN −1 .

(4.2)

To represent the constraints (4.2) in vector form, we have nine equations. From
(3.12) we know that the number of free variables corresponding to orientation kinematics are actually three: χ0 ∈ R3 . So, the boundary conditions (4.2) have to
be represented by three independent constraints. Let us define the maps SO(3) ∋
M 7→ logm(M ) ∈ so(3) and so(3) ∋ x 7→ (x)∨ ∈ R3 . Then the boundary condition
(4.2) is satisfied if
∨
(4.3) Orientation constraints: Cornt := logm Rf⊤ Ri F0 F1 F2 . . . FN −1
= 0.
We need to compute the gradient of the orientation constraints (4.3) w.r.t. momentum Πk . This can be done as follows: From (4.3) we know that
∧

eCornt = Rf⊤ Ri F0 F1 . . . Fk . . . FN −1 ,
leading to

∧

∧
= Rf⊤ Ri F0 F1 . . . DΠik (Fk ) . . . FN −1 .
eCornt DΠik Cornt
After algebraic manipulations we get
∧
= FN⊤−1 FN⊤−2 . . . Fk⊤ DΠik (Fk )Fk+1 . . . FN −1 ,
DΠik Cornt
⊤ x = A⊤ x̂A we conclude that
and using the property Ad
∨

DΠik Cornt = FN⊤−1 FN⊤−2 . . . Fk⊤ DΠik Fk .

12

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR

Therefore, the gradient of the orientation constraints (4.3) w.r.t. momentum Πk is

(4.4)
DΠk Cornt = DΠ1k Cornt DΠ2k Cornt DΠ3k Cornt .

We notice that the orientation constraints (4.3) can be computed only when we
construct the matrix Fk for a given value of the momentum vector Πk . Similarly,
the gradient of the orientation constraints (4.4) can be computed only when we find
the derivative of Fk w.r.t. Πik for k = 1, 2, . . . , N − 1, and i = 1, 2, 3.
4.2. Determining Fk and DΠik Fk in terms of momentum Πk . In order to
calculate the orientation constraints (4.3), Fk is obtained in terms of Πk . We know
that Fk can be obtained from Πk by solving the implicit equation (4.1). To solve
this implicit form, we choose quaternions to parameterize the matrix Fk . Let
 2

2q1 q3 + 2q0 q2
q0 + q12 − q22 − q32
2q1 q2 − 2q0 q3
q02 − q12 + q22 − q32
2q2 q3 − 2q0 q1  ,
Fk :=  2q1 q2 + 2q0 q3
2q1 q3 − 2q0 q2
2q2 q3 + 2q0 q1
q02 − q12 − q22 + q32

and the inertia matrix J and Jd are defined as:
 x

x

−I + Iy + Iz
0
0
I
0 0
1
.
0
Ix − Iy + Iz
0
J :=  0 Iy 0  ,
Jd := 
2
x
y
z
z
0
0
I +I −I
0 0 I

⊤
be the momentum of the body at the kth instant, and
Let Πk := Π1k Π2k Π3k
⊤
q = (q0 , q1 , q2 , q3 ) ; then (4.1) can be represented in the form of nonlinear algebraic
equations as follows:


2q2 q3 (Iz − Iy ) + 2q0 q1 Ix − hΠ1k
2q1 q3 (Ix − Iz ) + 2q0 q2 Iy − hΠ2 
k
g (q(Πk ), Πk ) := 
(4.5)
2q1 q2 (Iy − Ix ) + 2q0 q3 Iz − hΠ3k  = 0.
q02 + q12 + q22 + q32 − 1
For a fixed value of momentum Πk , the system of equation (4.5) has quaternions
as unknown parameters which can be numerically found using Newton’s method
at each instant of time. Since Fk is represented in terms of quaternions, we first
need to find the variations of the quaternions in terms of the momentum vector Πk ;
these can be obtained by taking derivative of (4.5) w.r.t. Πk as
Dq g (q(Πk ), Πk ) DΠk q(Πk ) + DΠk g (q(Πk ), Πk ) = 0,
resulting in

Dq g

(4.6)

where

 ∂q0

∂Π1
 ∂q1k
 ∂Π1
k
(q(Πk ), Πk ) 
 ∂q21
 ∂Πk
∂q3
∂Π1k

|


2q1 Ix
2q2 Iy
Dq g (q(Πk ), Πk ) = 
2q3 Iz
2q0

∂q0
∂Π2k
∂q1
∂Π2k
∂q2
∂Π2k
∂q3
∂Π2k

{z

DΠk q(Πk )

∂q0 
∂Π3k
∂q1 

∂Π3k 
∂q2 
∂Π3k 
∂q3
∂Π3k

}



h
0
=
0
0

2q0 Ix
2q3 (Iz − Iy )
x
z
2q3 (I − I )
2q0 Iy
y
x
2q2 (I − I ) 2q1 (Iy − Ix )
2q1
2q2

0
h
0
0


0
0
,
h
0


2q2 (Iz − Iy )
2q1 (Ix − Iz )
.

2q0 Iz
2q3

ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS

13

The matrix DΠk q can be obtained by solving the linear system (4.6). The derivative
of the matrix Fk w.r.t. Πik is obtained using the chain rule as
DΠik Fk =

(4.7)

3
X

Dqn Fk

n=0

∂qn
.
∂Πik

We now discuss the reduction of the difference equation model (3.2), (3.12) to the
momentum and comomentum dynamics respectively and later the reduced model
is scaled by change of variables so as to make the model invariant under change
in the step length h. Invariance of the difference equation model means that for
a particular manoeuvre the optimal trajectory and corresponding Lagrange multipliers remain identical for different step lengths h. This matter is quite essential
because it largely affects the region of convergence and order of convergence of the
algorithm [12].
4.3. Scaling and model reduction. First we discuss about model reduction and
later scaling of the reduced model by appropriate change of variable. Let us define
a new variable ζk := 2h−2 (tr (Fk ) I3×3 − Fk ) χk then (3.12) can be written as


Nk
⊤Π
λk +
(4.8)
ζk = 0
βk ⊙ Πk − λk−1 + Fk − hNk F\
k
k
h
(4.9)
ζk−1 − Fk ζk = 0
B⊤

where Nk = hk . From (4.9) we conclude ζk = Q⊤
k ζ0 such that Qk = F1 F2 . . . Fk .
So (4.8) and (3.2) can be further reduced to the system of difference equations
(4.10)
⊤
 Πk+1 − Fk Πk − hůk 
Q⊤
\
⊤
βk+1 ⊙ Πk+1 − λk + Fk+1 − hNk+1 Fk+1
Πk+1 λk+1 + Nk+1 k+1
h ζ0

!

= 0.

Assume for the moment that for a manoeuvre the ith element of the control vector
ůk saturates (i.e., ůik = ci ) at the kth instant of time. Clearly, from (3.14)
i
λik ≥ ch . If the step length h is small, λik will take very large values, which makes
the difference equations (4.10) stiff. To avoid this situation, we define a variable
γk := hλk , and modify the difference equation (4.10) to
(4.11)
⊤
 Πk+1 − Fk Πk − hůk 
\
⊤
hβk+1 ⊙ Πk+1 − γk + Fk+1 − hNk+1 Fk+1
Πk+1 γk+1 + Nk+1 Q⊤
k+1 ζ0

where

(4.12)


ůik = − min ci , γki

sgn(γki ),

Nk

!

= 0,

−1

Fk .
= tr Jd Fk⊤ I3×3 − Jd Fk⊤

Note that Fk and its gradient with respect to momentum can be obtained as discussed in §4.2. In the following section we employ multiple shooting method to
solve the system of difference equations (4.11) with boundary conditions (3.3) and
constraints (3.4).
5. Multiple Shooting Methods
Shooting methods were mainly developed for solving ordinary differential equations or difference equations with given boundary conditions. An initial guess is
taken for the unknown initial values of the differential or difference equation variables. Then the variables are computed at the terminal time and compared with a
known value of the variables at the boundaries. Then the initial guess is improved
at each iteration to match with the known boundary values. The multiple shooting

14

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR

methods, the time domain is divided into sub-intervals (time domain decomposition), and the boundary value problems are solved for each sub-interval with the
condition that the boundary values at the common points of the adjacent intervals
are the same. Multiple shooting methods are generalizations of the single shooting
method in the sense that multiple shooting with a single time interval is equivalent
to the shooting method [15].
Multiple shooting methods have many advantages over single shooting methods. The former are more stable, and hence can be applied to stiff problems, and
the time domain decomposition allows one to introduce the initial guess to the
problem with prior knowledge. Multiple shooting methods allow one to compute
the solution of the differential equation at individual intervals, which can be very
efficient in computation using parallel architecture [15].

5.1. A quick introduction to multiple shooting methods. Multiple shooting
methods constitute generalizations of the single shooting method, in which two
point boundary value problems are solved at each iteration for the subintervals of
time domain simultaneously. Let
(5.1)
(5.2)

x ∈ Rn

ẋ = f (x, t)
n

n

B :R ×R →R

n

with boundary conditions

such that B(x(a), x(b)) = 0.

Let the time domain be decomposed into N sub-interval as
t0 = τ0 < τ1 . . . τN −1 < τN = tf ,
and let us consider (N + 1) variables s0 , s1 , . . . , sn , known as the multiple shooting
variables. These multiple shooting variables are the guessed initial values of the
dependent variable x defined in (5.1) at the specified time instants. We now define
initial value problems for each sub-interval as
(5.3)

ẋk = f (xk , t) such that xk : [τk , τk+1 ] → Rn with

(5.4)

xk (τk ) = sk for k = 0, 1, . . . , (N − 1).


Note that the solution xk ( · , sk )|k = 0, 1, . . . , (N − 1) to the initial value problems (5.3) can be a solution to the boundary value problem (5.1) only if the solution
xk of the interval [τk , τk+1 ] matches with the initial condition for the next interval i.e. xk (τk+1 , sk ) = sk+1 . This condition is known as the matching condition.
These matching conditions for each interval can be combined together and can be
represented in stacked form as

(5.5)



s1 − x0 (τ1 , s0 )
s2 − x1 (τ2 , s1 )
..
.









F (s) := 
 = 0,


sN − xN −1 (τN , sN −1 )
B(x(a), x(b))


with xk ( · , sk ) the solution of initial value problem (5.3), and s := s0 , s1 , . . . , sN .
The algebraic equations (5.5) can be solved using Newton type algorithms for multivariable functions is described in Algorithm 1.

ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS

15

Algorithm 1 Multiple shooting method
1:
2:
3:
4:
5:
6:
7:
8:
9:

procedure Root of the function F
Choose appropriate initial guess s and the error tolerance bound δ > 0.
Stopping Criterion:
if kF (s)k∞ < δ then Stop.
end if
Calculate ∆s by solving the linear system Ds F (s)∆s = −F (s).
Update the value of s = s + ∆s.
goto Stopping Criterion.
end procedure

5.2. Multiple shooting method for attitude dynamics. The system dynamics (4.11) with the boundary conditions (4.3) and (Π0 , ΠN ) = (Πi , Πf ) have to be
solved with additional inequality constraints (Πik )2 ≤ (di )2 along with the complimentary slackness conditions (3.10) defined by
(5.6)
(
Πk+1 − Fk⊤ Πk − hůk = 0,


\
⊤
Πk+1 γk+1 + Nk+1 Q⊤
hβk+1 ⊙ Πk+1 − γk + Fk+1 − hNk+1 Fk+1
k+1 ζ0 = 0,
(5.7)

n
boundary conditions ΠN − Πf = 0,

Π0 − Πi = 0,

Cornt = 0,

(5.8)

n

slackness conditions βki (Πik )2 − (di )2 = 0,

where

(Πik )2 ≤ (di )2 ,

−1

Fk ,
Nk = tr Jd Fk⊤ I3×3 − Jd Fk⊤

∨
.
:= logm Rf⊤ Ri F0 F1 F2 . . . FN −1


ůik = − min ci , γki

Cornt

βki ≥ 0,

sgn(γki ),

First we consider the case in which state constraints are not active (i.e., βk = 0
for all k). Then (5.8) is trivially satisfied. If state constraints are not active, then
the necessary conditions are defined by (5.6) with βk = 0 for all k, and (5.7). We
now represent the matching and boundary condition for the system of difference
equations (5.6) and (5.7) in terms of Πk , γk and ζ0 and solve the system of nonlinear
algebraic equations comprises of matching conditions.
5.2.1. Matching conditions for multiple shooting methods. The system of difference
equations (5.6) along with the boundary conditions (5.7) will result in the following
set of matching conditions, which can be solved using Newton’s root finding algorithm with quadratic convergence rate [28].
Solve:
⊤
⊤
⊤
Cornt
(5.9) M(X) := Σ1 Ξ1 . . . Σk Ξk . . . ΣN ΞN Cmtm
= 0,
where

⊤
⊤
γk⊤ . . . Π⊤
γN
ζ0⊤ ,
N
⊤
Σk+1 := Πk+1 − Fk⊤ Πk − hůk ,
⊤



\
⊤
,
Ξk+1 := Nk+1 Q⊤
k+1 ζ0 + Fk+1 − hNk+1 Fk+1 Πk+1 γk+1 − γk


∨
ΠN − Πf
Cmtm :=
, Cornt := logm Rf⊤ Ri F0 F1 F2 . . . FN −1
.
Π0 − Πi
X := Π⊤
0

γ0⊤

. . . Π⊤
k

16

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR

Start

Initialize X and the
Error Tol δ > 0

update X = X + ∆X

Evaluate ∆X by solving
the Linear System:
DX M(X)∆X = −M(X)

Evaluate the Matching
Conditions M(X) and
its gradient DX M(X)

No

Is
kM(X)k∞ ≤
δ?

Yes
Stop
Figure 2. Newton’s root finding algorithm
Let DX M be the Jacobian matrix of the matching conditions M (5.9). Assuming
Xn to be the solution of the system M at the nth iteration of the Newton’s root
finding algorithm, the (n + 1)th iteration is given by
Xn+1 = Xn + ∆Xn ,
where ∆Xn is a solution of the linear system DX M(Xn )∆Xn = −M(Xn ).
It is important to note that Newton’s update ∆Xn is not necessarily unique;
∆Xn is unique if and only if DX M(Xn ) is invertible. Invertibility of the matrix
DX M is proved for a special case (when only momentum dynamics are considered)
in Appendix C.
Notice that ůk defined in (4.12) is not differentiable, but it is Lipschitz continuous. So, we take its generalized gradient [8], defined by

i
i

−1 if γk ≤ c ,
for i = 1, 2, 3.
(Dγk ůk )ii = −1 if − γki ≥ −ci ,


0
otherwise

The resulting nonlinear algebraic equations (5.9) can be solved using the nonsmooth version of Newton’s method in [30], and this is illustrated with the help of
the flowchart shown in Figure 2.

5.2.2. Multiple shooting method with state constraints. Here we discuss the general
case in which momentum constraints are active. By the complementary slackness
conditions (5.8), it is clear that if the momentum constraints are not active, then
βki = 0 for all i = 1, 2, 3, and k = 1, 2, . . . , N − 1. On the other hand, if the

ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS

17

slack variable βki > 0, then the momentum variable Πik corresponding to the slack
variable βki lies on the boundary (i.e., Πik = ±di ) [6].
The original Newton’s root finding algorithm cannot handle momentum inequality constraints; we propose a modified strategy for handling the momentum inequality constraints. The algorithm has three phases:
(1) Projection of the states to the feasible region: After each Newton’s
iteration, the momentum variables {Πk }N
k=0 are projected onto the feasible
domain.
(2) Identification of active constraints: Active constraints can be identified
by computing βki from (5.6). If βki > 0, then the inequality constraint
corresponding to the state Πik is active; else it is inactive.
(3) Enforcing active constraints: When βki > 0, then the equation with βki
in (5.6) is trivially satisfied. So, the matching condition Ξik will be replaced
with the corresponding active constraint Πik = di for the next iteration.
Let
X := Π⊤
0

γ0⊤

. . . Π⊤
k

γk⊤

. . . Π⊤
N

⊤
γN

ζ0⊤

⊤

∈ R6N +9

be a vector in the augmented space; the feasible region Λ is defined as

Λ := X ∈ R6N +9

Πik ≤ di .


It is clear that
Λ is the intersection of the two half spaces X + := X Πik ≤ di

and X − := X Πik ≥ −di , hence a convex set. Let PΛ : R6N +9 → Λ be the
projection onto the feasible region Λ, and M̃ represents the active constraints,
which is a collection of equality state constraints, active inequality state constraints
and co-state constraints corresponding to inactive inequality state constraints. Our
modified multiple shooting method is described in Algorithm 2, and is illustrated
with the help of the flowchart given in Figure 3.

6. Numerical Experiments
The following data have been considered for the simulations: Principal Moment of inertia of the satellite (Ix , Iy , Iz )=(800, 1200, 1000) kgm2 , Sampling time
(T ) = 0.1 s, Range of angles (θ) can vary between [10◦ , 90◦ ] about any axis, Maximum torque or control bound (c) = (20, 20, 20) Nm, maximum momentum (d) =
(70, 70, 70) Nms, time duration (tmax ) can range from 5 s to 30 s.
An open loop optimal control profile
 is obtained
 for a maneouvre of reorienting the satellite for 90◦ about axis √13 , √13 , √13 from initial momentum Πi =
(30, −10, 10) Nms to desired momentum Πf = (0, 0, 0) Nms in 19 s. The optimal
control profile along with the momentum and co-state vectors is shown in Figure 4.

18

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR

Algorithm 2 Modified multiple shooting method
1:
2:
3:
4:
5:

6:
7:

8:
9:

procedure Root of the function M̃ such that X ∈ Λ
Choose appropriate initial guess X and the error tolerance bound δ > 0.
Projection:
⊤
⊤
Compute X 7→ PΛ (X) := Π̃⊤
γ0⊤ . . . Π̃⊤
γk⊤ . . . Π̃⊤
γN
ζ0⊤
0
N
k

where Π̃ik = sgn(Πik ) min di , Πik .
Using (5.6), update co-states


(
⊤
Π̃k − Fk−1
Π̃k−1
i
i
i
<
c
,
−ů
if
ů
i
k
k
.
γ̃k−1
=
where ůk =
i
h
γk−1
elsewhere,
if Π̃ik = di then
Compute the slack variable using (5.6) as

i


Z(Π̃k )Ξk−1 (Π̃k , γ̃k , γ̃k−1 , ζ0 )
Fk
i
\
⊤
βk =
, where Z(Πk ) =
− Nk Fk Πk .
h
Π̃ik
end if
Define active constraints
M̃(X) := Σ1 Ξ̃1 . . . Σk Ξ̃k . . . ΣN
where
(
Πik − di if βki > 0,
Ξ̃ik :=
Ξik
otherwise.

10:

11:
12:
13:
14:
15:
16:

⊤
Cmtm

Ξ̃N

⊤
Cornt

⊤

,

Evaluate the matching conditions M̃(X̃) and its gradient DX M̃(X̃),
⊤
⊤
where X̃ := Π̃⊤
γ̃0⊤ . . . Π̃⊤
γ̃k⊤ . . . Π̃⊤
γ̃N
ζ0⊤ .
0
k
N
if M̃(X̃)

∞

< δ then Stop.

end if
Calculate ∆X̃ by solving the linear system DX M̃(X̃)∆X̃ = −M̃(X̃).
Update the value of X = X̃ + ∆X̃.
goto Projection.
end procedure
Optimal Control

Momentum

20
ux
uy

10

Πx

60

Πy

uz

z

Π

40
0
20
−10
0
−20
0

5

10

15

0

5

CoMomentum
ζx

0

γy

ζy
−50

z

γ

10

15

CoRotation

γx

20

10

z

ζ

−100

0

−150
−10
−200
−20
0

5

10

15

Figure 4. Rotation of 90◦ about

0



5

10

√1 , √1 , √1
3
3
3



15

axis with ini-

tial momentum Πi = (30, −10, 10) Nms, desired momentum Πf
= (0, 0, 0) Nms, time duration = 19 s.

ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS

19

Start

Initialize X and the Error Tol δ > 0

Evaluate PΛ (X), γ̃k , βk

update X = X̃ + ∆X̃

Define active constraints M̃ and X̃

Evaluate ∆X̃ by solving
the Linear System:
DX M̃(X̃)∆X̃ = −M̃(X̃)

Evaluate the matching conditions
M̃(X̃) and its gradient DX M̃(X̃)

No

Is
M̃(X̃)
δ?

∞

≤

Yes
Stop
Figure 3. Modified multiple shooting algorithm for problems
with state constraints
It is clear from the Figure 4 that the control saturates when the absolute values
of the co-state vectors corresponding to the momentum are more than the control
bounds. If the momentum saturates, then the control along that axis goes to zero.
7. Conclusion and Future Direction
A tractable solution to a class of optimal control of spacecraft attitude maneuvers under control and momentum constraints has been presented. We used
discrete mechanics to discretize the continuous time model, which has many advantages over the conventional discretization schemes like Euler’s steps. In particular,
the model so obtained preserves conserved quantities of the body like momentum
and energy. It is efficient in terms of numerics because the model reduces to the
momentum dynamics only. A new multiple shooting algorithm has been proposed
to solve the system of difference equations obtained as first order necessary conditions using variational analysis. This algorithm can be used to solve optimal
control problems with state inequality constraints. In future, convergence results
for the modified multiple shooting methods will be explored. The proposed multiple shooting method uses Newton’s root finding algorithm at each step for finding
the roots of the matching conditions. Newtons’s method has a quadratic rate of

20

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR

convergence in a neighborhood of the true solution, and therefore it is worth exploring the convergence rates of our algorithm in a neighborhood of the matching
conditions. It is important to note that our algorithm can handle box constraints
only; an extension to more general constraints is under development.
Appendix A. Proof of Lemma 3.1
Lemma 3.1. The matrix (tr (Fk Jd ) I3×3 − Fk Jd ) is invertible if
 s

2d3 + d2 − d1
kξk k
<
,
cos
2
2(d3 + d2 )
b

where Fk = eξk ,

ξk ∈ R3 , and Jd = diag(d1 , d2 , d3 ) such that 0 < d1 ≤ d2 ≤ d3 .

Proof. If tr (Fk Jd ) = 0 then the matrix (tr (Fk Jd ) I3×3 − Fk Jd ) is invertible. If
tr (F
k Jd ) 6= 0 then (tr (Fk Jd ) I3×3 − Fk Jd ) is invertible if and only if the maFk Jd
trix I3×3 − tr(F
is invertible. Using Banach lemma [20, p. 193], the matrix
k Jd )

Fk Jd
Fk Jd
is invertible if tr(F
< 1.
I3×3 − tr(F
k Jd )
k Jd )
kJd k
kFk Jd k
Fk Jd
≤
=
tr (Fk Jd )
|tr (Fk Jd )|
|tr (Fk Jd )|

Then
q

Fk Jd
tr(Fk Jd )

2d3 +d2 −d1
2(d3 +d2 )

<1

holds if

kJd k < |tr (Fk Jd )|. Now we show that if cos



kξk k
2



then condition kJd k < |tr (Fk Jd )| is true. Let us choose the quater-

nions (q0, q1 , q2, q3 ) ∈ R4 as a parametrization

for the rotation matrix Fk where
⊤
kξk k
kξk k
e
and (q1 , q2 , q3 ) = ξk sin 2 , the vector ξke represent the unit
q0 = cos 2
vector corresponding to vector ξk . Then
|tr (Fk Jd )| = d1 (q02 + q12 − q22 − q32 ) + d2 (q02 − q12 + q22 − q32 ) + d3 (q02 − q12 − q22 + q32 )

(A.1)

≥ d1 (q02 + q12 − q22 − q32 ) + d2 (q02 − q12 + q22 − q32 ) + d3 (q02 − q12 − q22 + q32 )

= d1 + d2 + d3 − 2 (d1 + d2 )q12 + (d2 + d3 )q22 + (d3 + d4 )q32
≥ d1 + d2 + d3 − 2f˜,

where
f˜ = maximize (d1 + d2 )x + (d2 + d3 )y + (d3 + d4 )z

(A.2)

x,y,z

subject to
x + y + z = 1 − q02 ,
x ≥ 0, y ≥ 0, z ≥ 0.
The optimization problem defined in (A.2) is a linear programming problem and
hence the optimum value will be attained on the vertices of the feasible region.
Without loss of generality,
assume that 0 < d1 ≤ d2 ≤ d3 . Then
f˜ = (d3 +d2 )(1−q02 ).
 q

kξk k
2

and q0 = cos kξ2k k , one can conclude that
3 +d2 −d1
q02 < 2d2(d
. Substituting the value of f˜ to (A.1) gives
3 +d2 )

Given that cos

<

2d3 +d2 −d1
2(d3 +d2 )

|tr (Fk Jd )| ≥ d1 + d2 + d3 − 2(d3 + d2 )(1 − q02 ) > d3 = kJd k .


<

ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS

21

Appendix B. Proof of Claim 3.2
Claim 3.2. Consider the equality Rk⊤ Rk+1 = Fk . If we assume that the step length
h is small enough such that the relative orientation Rk⊤ Rk+1 between two adjacent
time instances tk and tk+1 is less than π2 , i.e.,
π
b
b
kξk k , kζk k < where Rk⊤ Rk+1 = eζk , Fk = eξk and ζk , ξk ∈ R3 ,
2
then the aforementioned equality is satisfied if and only if the skew symmetric parts
of both sides are identical.
Proof. This claim can be proved using the Rodrigues’s formula. Let A := Rk⊤ Rk+1 ∈
SO(3) and B := Fk ∈ SO(3) then there exist vectors a, b ∈ R3 such that
A = e(kakbae )

(B.1)

b

and B = e(kbkbe )

where xe is a unit vector corresponding to the vector x ∈ R3 . Using the Rodrigues’s
formula [4]
A = e(kakbae ) := I3×3 + sin(kak)b
ae + b
ae b
a⊤
e (cos(kak) − 1)

we obtain, the skew-symmetric parts of A and B as

A − A⊤
B − B⊤
= sin(kak)b
ae and S(B) :=
= sin(kbk)bbe .
2
2
Given that S(A) = S(B) and since we know that so(3) ∼ R3 , from (B.2) we have
(B.2)

S(A) :=

(B.3)

sin(kak)ae = sin(kbk)be

case 1:
If kbk = 0 then from (B.3) sin(kak)ae = 0 iff a = 0 which means a = b because
kak ≤ π2 . Using the exponential map (B.1) one can conclude that A = B.
case 2:
If kbk 6= 0 then from (B.3)
sin(kak)
ae = b e
sin(kbk)
which is true iff
(B.4)

b e = ae

and

(B.5)

sin(kak)
= ±1.
sin(kbk)

The map sin :] − π2 , π2 [→] − 1, 1[ is bijective. Given that kbk <
(B.5) that sin(kak) = sin(kbk). Hence

π
2,

we conclude from

kbk = kak

(B.6)

From (B.4) and (B.6) it is clear that a = b. Using the exponential map (B.1), one
can conclude that A = B.

Appendix C. Proof of Theorem C.1
Let us discuss a case in which, only momentum dynamics is considered . In this
case, matching conditions (5.9) modifies to
⊤
⊤
(C.1)
= 0,
C(Y ) := Σ1 Ξ̃1 . . . Σk Ξ̃k . . . ΣN Ξ̃N Cmtm
where

Y := Π⊤
0

γ0⊤

. . . Π⊤
k

γk⊤

. . . Π⊤
N

⊤
γN

⊤

, Σk+1 := Πk+1 − Fk⊤ Πk − hůk

⊤

,

22

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR



⊤



ΠN − Πf
\
⊤
.
, Cmtm :=
Ξ̃k+1 := −γk + Fk+1 − hNk+1 Fk+1
Πk+1 γk+1
Π0 − Πi

We define the gradient of the matching

−J1 K1
 0
−J2

 ..
..

.
DY C(Y ) :=  .
(C.2)
 0
0

 0
0
Bi
0
where

conditions (C.1) as
0
K2
..
.
0
0
0

...
...
..
.

0
0
..
.

0
0
..
.

. . . −JN −1
...
0
...
0

KN −1
−JN
0







,
0 

KN 
Bf




I3×3
0
−Dγk−1 Σ⊤
k
,
, Kk :=
I3×3
DΠk Ξ̃⊤
Dγk Ξ̃⊤
k
k




0 0
I3×3 0
f
=
, B := DΠN Cmtm =
.
I3×3 0
0 0


−DΠk−1 Σ⊤
k
Jk :=
0
Bi := DΠ0 Cmtm

0
0
..
.

Theorem C.1. The Jacobian matrix DY C is non-singular if there exist a time
instant k such that the control ůk is not saturated.
Proof. The matrix DY C is similar to the upper triangular matrix
obtained using the elementary row transformations.

−J1 K1
0 ...
0
0
 0
−J
K
.
.
.
0
0
2
2

 ..
..
..
..
..
..

.
.
.
.
.
e ) :=  .
(C.3)
DY C(Y
 0
0
0
.
.
.
−J
K
N
−1
N
−1

 0
0
0 ...
0
−JN
0
0
0 ...
0
0

DY Ce which is
0
0
..
.







,
0 

KN 
X

e ) is invertible if and only if Jk is invertible for k =
Note that the matrix DY C(Y
1, 2, . . . , N , and X is invertible. First we prove that Jk is invertible,
(C.4)


−DΠk−1 Σ⊤
k
Jk =
0

−Dγk−1 Σ⊤
k
I3×3



=

⊤

\
⊤
Fk−1 − hNk−1 Fk−1
Πk−1
0

hDγk−1 ůk−1
I3×3

!



⊤Π
Lemma C.2. The matrix Fk − hNk F\
k is invertible if and only if the matrix
k
(tr (Fk Jd ) I3×3 − Fk Jd ) is invertible.

D. By Lemma C.2 and Lemma 3.1,
A proof of Lemma C.2 isprovided in Appendix

⊤Π
we know that the matrix Fk − hNk F\
is
invertible.
Hence, from (C.4) it is
k
k
easy to conclude that the matrix Jk is invertible. Let us define the inverse of the
matrix Jk as

 −1


−1
I
0
ak−1 −ak−1
bk−1
,
Jk−1 :=
, Kk := 3×3
ck
a⊤
0
I3×3
k

where


⊤
⊤Π
,
ak = Fk − hNk F\
k
k

bk = hDγk ůk ,

ck = DΠk Ξ̃k .

Taking the first order approximation in h, the matrix X can be approximated as
X = B f + B i J1−1 K1 J2−1 K2 . . . JN−1 KN

.

ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS

23


I3×3
∗


0
,
s

P
N −1
⊤
where s = −r0−1
r
b
r
i+1
i
i+1 such that ri = aN aN −1 . . . ai . Now it is clear
i=0
(C.5)

≈

from (C.5) that the matrix X is invertible if the matrix s is invertible. As the matrix
ai is invertible, so the matrix ri is invertible for all i = 0,1,. . . ,N. So, it is enough to
P −1
⊤
prove that the matrix N
i=0 ri+1 bi ri+1 is invertible for s being invertible. Note that
the matrix bi is a negative semi-definite matrix for i = 1, 2, . . . , N − 1. Assuming
that there exist a time instant k such that the control ůk is not
saturated means
PN
−1
⊤
that bk is negative definite. So, it is concluded that the matrix i=0 ri+1 bi ri+1
is
negative definite and hence invertible.

Appendix D. Proof of Lemma C.2


⊤
Lemma C.2. The matrix Fk − hNk F\
k Πk is invertible if and only if the matrix
(tr (Fk Jd ) I3×3 − Fk Jd ) is invertible.

⊤d
⊤
Proof. We know that F\
k hΠk = Fk hΠk Fk . Then
 using the implicit equation
T
\
⊤
d
hΠk = Fk Jd − Jd Fk , the matrix Fk − hNk Fk Πk can be rewritten as

 

⊤d
⊤Π
hΠ
F
=
F
−
N
F
Fk − hNk F\
k
k
k
k
k
k
k


−1
⊤
dk Fk
= Nk Fk Fk (Nk ) − hΠ



dk Fk
= Nk Fk⊤ tr Jd Fk⊤ I3×3 − Jd Fk⊤ − hΠ

(D.1)

= Nk Fk⊤ (tr (Fk Jd ) I3×3 − Fk Jd ) Fk

−1
= tr Jd Fk⊤ I3×3 − Jd Fk⊤
(tr (Fk Jd ) I3×3 − Fk Jd ) Fk .

By (D.1) it is obvious to conclude
 (tr (Fk Jd ) I3×3 − Fk Jd ) is invert that the matrix
\
⊤

ible if and only if the matrix F − hN F Π is invertible.
k

k

k

k

References
[1] A. A. Agrachev and Y. L. Sachkov, Control Theory from the Geometric Viewpoint,
Springer-Verlag New York, 2003.
[2] M. Athans, P. Falb, and R. Lacoss, Time-, fuel-, and energy-optimal control of nonlinear
norm-invariant systems, Automatic Control, IEEE Transactions on, 8 (1963), pp. 196–202.
[3] S. Bhatt, N. Bedrossian, K. Longacre, and L. Nguyen, Optimal propellant maneuver
flight demonstrations on ISS, Guidance, Navigation, and Control (GNC) Conference AIAA,
(2013).
[4] A. M. Bloch, P. E. Crouch, J. E. Marsden, and T. S. Ratiu, The symmetric representation of the rigid body equations and their discretization, Nonlinearity, 15 (2002), p. 1309.
[5] A. M. Bloch, I. Hussein, M. Leok, and A. Sanyal, Geometric structure-preserving optimal
control of a rigid body, Journal of Dynamical and Control Systems, Springer, 15 (2009),
pp. 307–330.
[6] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004.
[7] F. Bullo and A. D. Lewis, Geometric control of mechanical systems - modeling, analysis,
and design for simple mechanical control systems, Springer-Verlag New York, 2005.
[8] F. H. Clarke, Generalized gradients of Lipschitz functionals, Advances in Mathematics, 40
(1981), pp. 52 – 67.
[9] M. V. Dixon, T. N. Edelbaum, J. Potter, and W. E. Vandervelde, Fuel optimal reorientation of axisymmetric spacecraft, Journal of Guidance, Control, and Dynamics, AIAA, 7
(1970), pp. 1345–1351.
[10] A. Fleming, P. Sekhavat, and I. M. Ross, Minimum-time reorientation of a rigid body,
Journal of Guidance, Control, and Dynamics, AIAA, 33 (2010), pp. 160–170.

24

KARMVIR SINGH PHOGAT, DEBASISH CHATTERJEE, RAVI N. BANAVAR

[11] M. Gerdts, Global convergence of a nonsmooth Newton method for control-state constrained
optimal control problems, SIAM Journal on Optimization, 19 (2008), pp. 326–350.
[12] P. E. Gill, W. Murray, and M. H. Wright, Practical Optimization, Academic Press, 1997.
[13] A. Guiggiani, I. Kolmanovsky, P. Patrinos, and A. Bemporad, Fixed-point constrained
model predictive control of spacecraft attitude, American Control Conference (ACC), (2015),
pp. 2317–2322.
[14] R. Gupta, U. V. Kalabić, S. D. Cairano, A. M. Bloch, and I. V. Kolmanovsky, Constrained spacecraft attitude control on SO(3) using fast nonlinear model predictive control,
2015 American Control Conference (ACC), (2015), pp. 2980–2986.
[15] H. K. Hesse, Multiple shooting and mesh adaptation for PDE constrained optimization problems, PhD thesis, Ruprecht-Karls-University, Heidelberg, June 2008.
[16] D. Holm, T. Schmah, C. Stoica, and D. Ellis, Geometric Mechanics and Symmetry: From
Finite to Infinite Dimensions, OUP Oxford, 2009.
[17] I. Ioslovich, Arbitrary fuel-optimal attitude maneuvering of a non-symmetric space vehicle
in a vehicle-fixed coordinate frame, Automatica, (2003), pp. 557 – 567.
[18] D. Kim and J. D. Turner, Suboptimal asymmetric spacecraft maneuvers using two controls,
Journal of Guidance, Control, and Dynamics, AIAA, 36 (1995), pp. 1855–1860.
[19] M. B. Kobilarov and J. E. Marsden, Discrete geometric optimal control on lie groups,
Robotics, IEEE Transactions on, 27 (2011), pp. 641–655.
[20] P. D. Lax, Functional Analysis, John Wiley & Sons, 2002.
[21] T. Lee, M. Leok, and N. McClamroch, Optimal attitude control of a rigid body using
geometrically exact computations on SO(3), Journal of Dynamical and Control Systems, 14
(2008), pp. 465–487.
[22] T. Lee, M. Leok, and N. McClamroch, Time optimal attitude control for a rigid body,
American Control Conference, (2008), pp. 5210–5215.
[23] T. Lee, N. H. McClamroch, and M. Leok, A lie group variational integrator for the
attitude dynamics of a rigid body with applications to the 3d pendulum, IEEE Conference on
Control Applications, (2005), pp. 962–967.
[24] U. Lee and M. Mesbahi, Quaternion-based optimal spacecraft reorientation under complex attitude constrained zones, AAS/AIAA Astrodynamics Specialist Conference, (2013),
pp. 2013–836.
[25] J. E. Marsden and T. S. Ratiu, Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems, Springer-Verlag New York, 2010.
[26] J. E. Marsden and M. West, Discrete mechanics and variational integrators, Acta Numerica, 10 (2001), pp. 357–514.
[27] S. Ober-Blöbaum, O. Junge, and J. E. Marsden, Discrete mechanics and optimal control:
An analysis, ESAIM: Control, Optimisation and Calculus of Variations, 17 (2011), pp. 322–
352.
[28] B. T. Polyak, Newton’s method and its use in optimization, European Journal of Operational
Research, 181 (2007), pp. 1086–1096.
[29] R. Proulx and I. Ross, Time-optimal reorientation of asymmetric rigid bodies, Advances
in the Astronautical Sciences, AAS, (2002), pp. 1207–1227.
[30] L. Qi and J. Sun, A nonsmooth version of Newton’s method, Mathematical Programming,
58 (1993), pp. 353–367.
[31] I. M. Ross and M. Karpenko, A review of pseudospectral optimal control: from theory to
flight, Annual Reviews in Control, 36 (2012), pp. 182 – 197.
[32] A. Saccon, J. Hauser, and A. P. Aguiar, Optimal control on Lie groups: The projection
operator approach, IEEE Transactions on Automatic Control, 58 (2013).
[33] H. Schaub and J. L. Junkins, Analytical Mechanics of Space Systems, American Institute
of Aeronautics and Astronautics, 2003.
[34] R. C. Scrivener, Sandra L.and Thompson, Survey of time-optimal attitude maneuvers,
Journal of Guidance, Control, and Dynamics, AIAA, 17 (1994), pp. 225–233.
[35] E. Trélat, Optimal control and applications to aerospace: some results and challenges,
Journal of Optimization Theory and Applications, 154 (2012), pp. 713–758.
[36] P. Tsiotras, M. Corless, and J. M. Longuski, A novel approach to the attitude control of
axisymmetric spacecraft, Automatica, IFAC, 31 (1995), pp. 1099–1112.
[37] J. Vandersteen, M. Diehl, C. Aerts, and J. Swevers, Spacecraft attitude estimation
and sensor calibration using moving horizon estimation, Journal of Guidance, Control, and
Dynamics, 36 (2013), pp. 734–742.
[38] Y. Zhuang and H. Huang, Time-optimal trajectory planning for underactuated spacecraft
using a hybrid particle swarm optimization algorithm, Acta Astronautica, 94 (2014), pp. 690–
698.

ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS

Indian Institute of Technology Bombay, Mumbai, India.

25

