# Satellite Attitude Control Using Analytical Solutions to Approximations of the Hamilton-Jacobi Equation

Authors: Stefan LeBel; Christopher J. Damaren
Journal: ISRN Aerospace Engineering
 (2014)

---

## Abstract

The solution to the Hamilton-Jacobi equation associated with the nonlinearℋ
∞ control problem is approximated using a Taylor series expansion. A recently developed analytical solution method is used for the second-, third-, and fourth-order terms. The proposed controller synthesis method is applied to the problem of satellite attitude control with attitude parameterization accomplished using the modified Rodrigues parameters and their associated shadow set. This leads to kinematical relations that are polynomial in the modified Rodrigues parameters and the angular velocity components. The proposed control method is compared with existing methods from the literature through numerical simulations. Disturbance rejection properties are compared by including the gravity-gradient and geomagnetic disturbance torques. Controller robustness is also compared by including unmodeled first- and second-order actuator dynamics, as well as actuation time delays in the simulation model. Moreover, the gap metric distance induced by the unmodeled actuator dynamics is calculated for the linearized system. The results indicated that a linear controller performs almost as well as those obtained using higher-order solutions for the Hamilton-Jacobi equation and the controller dynamics.


---

## Body

## 1. Introduction

The attitude control problem is critical for most satellite applications and has thus attracted extensive interest. While many control methods have been developed to address this problem, most of them are concerned primarily with the optimality of attitude maneuvers [1–4]. In the present work, we shall focus on robust nonlinear control systems. We note that, throughout this paper, by nonlinear ℋ
∞ we mean the ℒ
2-gain of a nonlinear system.Control laws are generally developed based on mathematical models that are, at best, a close approximation of real-world phenomena. For such control methods to have any real practical value, they must be made robust with regard to unmodeled dynamics and disturbances that may act on the system. The study of robust control is therefore an essential part of the application of control theory to physical systems. In general, the development of an optimal nonlinear state feedback control law is characterized by the solution to a Hamilton-Jacobi partial differential equation (HJE) [5], while a robust nonlinear controller is obtained from the solution of one or more Hamilton-Jacobi equations [6–9]. However, no general analytical solution has yet been obtained to solve this optimization problem. Solutions have thus far only been obtained under certain conditions: in the case of linear systems with a quadratic performance index, the HJE reduces to the well-known algebraic Riccati equation (ARE). It is noted that the concept of dissipativity, which is closely related to optimal and robust control, is characterized by a Hamilton-Jacobi inequality [10–12].Extensive work has been carried out to approximate the solution of Hamilton-Jacobi equations through a Taylor series expansion [13–17]. Although such a series expansion results in an infinite-order polynomial, finite-order approximations can be used to obtain suboptimal solutions to an HJE. We also note the work in [18, 19] which uses series solution methods for nonlinear optimal control problems. It has been shown that a local solution to an HJE can be obtained by solving the ARE for the linear approximation of the system [6, 7, 20]. Methods that have been developed over the past decades to attempt to solve this problem include the Zubov procedure [21, 22], the state-dependent Riccati equation [23, 24], the Galerkin method for the equivalent sequence of first-order partial differential equations [25, 26], and the use of symplectic geometry to examine the associated Hamiltonian system [27]. However, one aspect that is lacking in all the above methods is an analytical solution to the approximate equations.The primary purpose of this paper is to develop robust nonlinear controllers based on analytical expressions for approximate solutions to the Hamilton-Jacobi equation. In particular, we shall provide analytical expressions for the second-, third-, and fourth-order terms of the approximation solution. These controllers are then compared through numerical simulations with existing methods from the literature for spacecraft attitude regulation [1–4]. Our objective is to examine the effects of different disturbances and uncertainties on the performance and robustness of the various controllers. More specifically, we include gravity-gradient and geomagnetic torques, as well as unmodeled actuator dynamics and actuation time delays. Moreover, we make use of the gap metric [28] to characterize the difference in the input-output (IO) map of the system induced by the unmodeled actuator dynamics. However, since we cannot calculate the gap between two nonlinear systems, we calculate the gap metric distance for the linearized system only. In contrast with some of the methods used for comparison, which were developed specifically to address the attitude control problem, the method presented in this paper is a general controller synthesis method and has also been applied to spacecraft formation flying [29].The outline of the paper is as follows. In Section2, a detailed description of the general class of systems is given, along with the controller synthesis method that is proposed. This controller is the solution of an appropriate nonlinear ℋ
∞ problem and is taken from the work of James et al. [9]. Then, the nonlinear equations of motion for the satellite attitude dynamics are given in Section 3. Section 4 presents simulation results using the proposed controller and comparisons are made with existing methods. Finally, some conclusions and suggestions for future work are stated in Section 5. We note here that some of these results also appear in past conference proceedings [30, 31] with the present paper containing improvements to the overall presentation.


## 2. Nonlinear Controller Synthesis Approach

This section provides the main results of our approach to robust nonlinear controller synthesis. We begin by describing the class of nonlinear systems with which we are concerned and define robustness in the gap metric. Then, the nonlinearℋ
∞ control problem is presented. Finally, the analytic solutions for the second-, third-, and fourth-order terms in the Taylor series approximation to the solution of the Hamilton-Jacobi equation (HJE) are stated.


### 2.1. Class of Nonlinear Systems

Consider the nonlinear system shown in Figure1. The plant is given by
(1)
P
:
{
x
˙
(
t
)
=
a
(
x
)
+
b
(
x
)
u
1
(
t
)
,
y
1
(
t
)
=
c
(
x
)
,
where a
(
x
), b
(
x
), and c
(
x
) are assumed to be smooth (i.e., C
∞) nonlinear functions of the plant states with a
(
0
)
=
0 and c
(
0
)
=
0 (i.e., x
(
t
)
=
0 is an equilibrium). The controller is described by
(2)
K
:
{
x
˙
c
(
t
)
=
a
c
(
x
c
)
+
b
c
(
x
c
)
y
2
(
t
)
,
u
2
(
t
)
=
c
c
(
x
c
)
,
where a
c
(
x
c
), b
c
(
x
c
), and c
c
(
x
c
) are smooth nonlinear functions of the controller states. Additionally, the following relations hold:
(3)
u
1
(
t
)
=
u
0
(
t
)
-
u
2
(
t
)
,
y
2
(
t
)
=
y
0
(
t
)
-
y
1
(
t
)
.
In the above, x
(
t
)
∈
ℝ
N is the plant state vector, x
c
(
t
)
∈
ℝ
N
c is the controller state vector, u
2
(
t
)
∈
ℝ
N
u is the control signal, u
0
(
t
)
∈
ℝ
N
u is the exogenous (disturbance) input, u
1
(
t
)
∈
ℝ
N
u is the actual plant input, y
1
(
t
)
∈
ℝ
N
y is the actual plant output, y
0
(
t
)
∈
ℝ
N
y is the reference and/or sensor noise signal, and y
2
(
t
)
∈
ℝ
N
y is the tracking error.Figure 1
Block diagram of system (1) with controller (2).The plant in (1) can be written as
(4)
P
⋆
:
{
x
˙
(
t
)
=
a
(
x
)
+
[
b
(
x
)
0
]
w
(
t
)
+
[
-
b
(
x
)
]
u
(
t
)
,
[
z
1
(
t
)
z
2
(
t
)
]
=
[
0
-
c
(
x
)
]
+
[
0
0
0
1
]
w
(
t
)
+
[
1
0
]
u
(
t
)
,
y
(
t
)
=
[
-
c
(
x
)
]
+
[
0
1
]
w
(
t
)
,
where
(5)
w
(
t
)
=
[
u
0
(
t
)
y
0
(
t
)
]
,
z
(
t
)
=
[
z
1
(
t
)
z
2
(
t
)
]
=
[
u
2
(
t
)
y
2
(
t
)
]
,
u
(
t
)
=
u
2
(
t
)
,
y
(
t
)
=
y
2
(
t
)
.
The generalized system in (4) is shown in Figure 2. The first equation in (4) defines the plant dynamics with state variable x
(
t
), control input u
(
t
), and subject to a set of exogenous inputs w
(
t
), which includes disturbances (to be rejected), references (to be tracked), and/or noise (to be filtered). The second equation defines the penalty variable z
(
t
) representing the outputs of interest, which may include a tracking error, a function of some of the exogenous variables w
(
t
), and a cost of the input u
(
t
) needed to achieve the prescribed control goal. The third equation defines the set of measured variables y
(
t
), which are functions of the plant state x
(
t
) and the exogenous inputs w
(
t
). As we shall see next, we will be concerned with defining an upper bound on the ℒ
2-gain from the disturbance inputs w
(
t
) to the outputs z
(
t
) of the system in (4).Figure 2
Block diagram of generalized system (4) with controller (2).


### 2.2. Robustness in the Gap Metric

In general, the plant and controller are assumed to be causal mappings from their respective inputs and outputs; that is,P
:
𝒰
→
𝒴 and K
:
𝒴
→
𝒰, which satisfy P
0
=
0 and K
0
=
0, where 𝒰 and 𝒴 are appropriate signal spaces. In particular, these mappings can be represented by (1) and (2). Thus, in the feedback configuration of Figure 1, the signals u
i belong to 𝒰 and the signals y
i belong to 𝒴, where i
∈
{
0,1
,
2
}. The input-output relations describing the plant and controller can be represented by their respective graphs:
(6)
𝒢
𝒫
:
=
{
[
u
P
u
]
:
u
∈
𝒰
,
P
u
∈
𝒴
}
,
𝒢
𝒦
:
=
{
[
K
y
y
]
:
K
y
∈
𝒰
,
y
∈
𝒴
}
,
where 𝒢
𝒫, 𝒢
𝒦
⊂
𝒲
=
𝒰
×
𝒴. For a plant P
0 with ℳ
0
:
=
𝒢
𝒫
0 and a plant P
1 with ℳ
1
:
=
𝒢
𝒫
1, the ρ-gap between these two plants is defined as [9, 28]
(7)
ρ
g
(
P
0
,
P
1
)
=
max
⁡
⁡
{
ρ
g
→
(
P
0
,
P
1
)
,
ρ
g
→
(
P
1
,
P
0
)
}
,
where
(8)
ρ
g
→
(
P
0
,
P
1
)
=
sup
⁡
0
≠
v
0
∈
ℳ
0
⁡
inf
⁡
v
1
∈
ℳ
1
⁡
(
∥
v
0
-
v
1
∥
2
+
∥
P
0
v
0
-
P
1
v
1
∥
2
)
1
/
2
(
∥
v
0
∥
2
+
∥
P
0
v
0
∥
2
)
1
/
2
and ρ
g
→
(
P
1
,
P
0
) is similarly defined. The norm used here, ∥
(
·
)
∥, is defined in [9].It has been shown [9, 32] that, for a stable feedback pair (
P
0
,
K
), if a system P
1 is such that
(9)
ρ
g
(
P
0
,
P
1
)
<
∥
Π
𝒩
∥
ℳ
0
∥
∞
-
1
,
then the feedback system (
P
1
,
K
) is also stable. The symbol Π
𝒩
∥
ℳ defines a parallel projection operator [32] and represents the closed-loop mapping from the disturbances w
(
t
) to the outputs z
(
t
) shown in Figure 2. The ℒ
2-gain of Π
𝒩
∥
ℳ is the induced norm defined by
(10)
∥
Π
𝒩
∥
ℳ
∥
∞
=
sup
⁡
0
≠
w
∈
ℒ
2
⁡
∥
z
∥
2
∥
w
∥
2
,
where the ℒ
2-norm is defined by ∥
z
∥
2
=
∫
0
∞
z
T
z
d
t.Thus, to minimize the effects of the disturbancesw
(
t
) on the outputs z
(
t
), which simultaneously achieves optimal robustness, the objective is to design a controller K such that ∥
Π
𝒩
∥
ℳ
0
∥
∞ is minimized. However, no closed-form solution exists to this optimal control problem. Instead, we shall be concerned with the suboptimal problem: given some constant scalar γ, design K such that ∥
Π
𝒩
∥
ℳ
0
∥
∞
<
γ. Note that a γ-iteration procedure can be applied to obtain the optimal solution within an arbitrary tolerance. It will be shown that the controller that satisfies this objective can be obtained from the solution to a single Hamilton-Jacobi equation.


### 2.3. Nonlinearℋ
∞ Control Problem

As indicated in the previous subsection, we are concerned with rendering∥
Π
𝒩
∥
ℳ
0
∥
∞
<
γ. In other words, we wish to bound the ℒ
2-gain from the disturbance inputs w
(
t
) to the outputs z
(
t
) of the system in (4). Without providing the details, it is noted that the HJE corresponding to (4) is dependent on γ and the quadratic term is sign-indefinite (i.e., is neither positive definite nor negative definite). However, by performing a certain transformation [9, 33], the generalized system can be written in a form where the outputs z
(
t
) are not explicit functions of the disturbances w
(
t
). The resulting HJE is no longer dependent on γ and is sign-definite. The system resulting from this transformation is described by
(11)
P
¯
⋆
:
{
x
˙
(
t
)
=
a
(
x
)
+
[
-
b
(
x
)
0
]
w
(
t
)
+
[
-
b
(
x
)
]
+
u
(
t
)
.
[
z
1
(
t
)
z
2
(
t
)
]
=
[
0
-
β
-
1
c
(
x
)
]
+
[
1
0
]
u
(
t
)
,
y
(
t
)
=
[
-
β
-
1
c
(
x
)
]
+
[
0
-
1
]
w
(
t
)
,
where β
=
1
-
γ
-
2. Now, the HJE corresponding to the system in (11) is given by
(12)
∇
V
(
x
)
A
(
x
)
-
1
2
∇
V
(
x
)
R
(
x
)
∇
V
T
(
x
)
+
1
2
Q
(
x
)
=
0
,
where ∇
V
(
x
) is the Jacobian matrix of the storage function V
(
x
) and is defined as
(13)
∇
V
(
x
)
=
∂
V
∂
x
T
=
row
i
{
∂
V
∂
x
i
}
with ro
w
i
{
·
} denoting a row matrix in the index i
∈
{
1,2
,
…
,
N
}. For the system in (11), we have the parameters [9]
(14)
A
(
x
)
=
a
(
x
)
,
R
(
x
)
=
b
(
x
)
b
T
(
x
)
,
Q
(
x
)
=
c
T
(
x
)
c
(
x
)
.
The name used in the literature for the HJE in (12) is not uniform. In the context of robust control, it is sometimes referred to as the Hamilton-Jacobi-Isaacs equation and when used in nonlinear optimal control it is often referred to as the Hamilton-Jacobi-Bellman equation. In [9] it is referred to as the Hamilton-Jacobi-Bellman-Isaacs equation. We will simply call it the HJE.It is shown in [9] that the controller K that solves the suboptimal ℋ
∞ problem for the plant P
⋆ (i.e., renders ∥
Π
𝒩
∥
ℳ
0
∥
∞
<
γ) is related to that (call it K
¯) which solves the same problem for the modified plant P
¯
⋆ as follows: K
=
K
¯
β. We now construct this K. Denote by V
+
(
x
) the unique smooth solution to the HJE (12) that satisfies V
+
(
x
)
≥
0 and V
+
(
0
)
=
0, with
(15)
+
[
A
(
x
)
-
R
(
x
)
∇
V
+
T
(
x
)
]
asymptotically stable. Similarly, we denote by V
-
(
x
) the unique smooth solution to the HJE that satisfies V
-
(
x
)
≤
0 and V
-
(
0
)
=
0, with
(16)
-
[
A
(
x
)
-
R
(
x
)
∇
V
-
T
(
x
)
]
asymptotically stable.Following the approach of James et al. [9], a local solution to the nonlinear ℋ
∞ control problem for the plant in (4) is obtained if the following conditions are satisfied.(1)
There exists aC
3 positive-definite function V
+
(
x
) defined in a neighbourhood of the origin with V
+
(
0
)
=
0 that satisfies the HJE (12).
(2)
There exists aC
3 negative-definite function V
-
(
x
) defined in a neighbourhood of the origin with V
-
(
0
)
=
0 that satisfies the HJE (12), and additionally satisfies
(17)
∇
2
{
-
∇
r
2
e
(
A
+
γ
-
2
β
-
2
R
∇
V
+
T
)
-
1
2
γ
2
β
-
2
+
1
2
γ
-
2
∇
r
2
e
R
∇
r
2
e
T
+
1
2
β
-
4
∇
V
+
R
∇
V
+
T
-
1
2
γ
2
β
-
2
Q
}
<
0
,
where
(18)
r
2
e
(
x
)
=
β
-
2
V
+
(
x
)
+
γ
2
V
-
(
x
)
.
(3)
There exists aC
2 function G
2
(
x
) such that
(19)
∇
r
2
e
(
x
)
G
2
(
x
)
=
γ
2
β
-
1
c
T
(
x
)
in a neighbourhood of the origin.The resulting nonlinear controller is given by(20)
K
:
{
x
˙
c
=
a
(
x
c
)
-
b
(
x
c
)
b
T
(
x
c
)
∇
V
+
T
(
x
c
)
x
˙
c
=
+
G
2
(
x
c
)
[
β
-
1
c
(
x
c
)
+
β
y
2
]
,
u
2
(
x
c
)
=
β
-
2
b
T
(
x
c
)
∇
V
+
T
(
x
c
)
.In the next subsection, we shall examine how to obtain analytical expressions for the approximate solution to the HJE in (12) required for this control law.


### 2.4. Analytical Solutions to HJE Approximation

Our approach for the synthesis of nonlinear controllers is based on the Taylor series approximation of the solution to the Hamilton-Jacobi equation, where each order of the controller is built using the previous orders. The following notation will be used:ro
w
k
{
·
} denotes a row matrix with index k, col
k
{
·
} denotes a column matrix with index k, and ma
t
m
n
{
·
} denotes a matrix with row index m and column index n. It should be emphasized that the symbols k, m, and n used here are dummy indices. In general, A
(
i
,
j
) refers to the (
i
,
j
)th entry of the matrix A. We will denote the positive-semidefinite solution V
+
(
x
) of the HJE (12) simply by V
(
x
).Consider the nonlinear system in (1). We begin by making the assumptions that b
(
x
)
=
B and c
(
x
)
=
C
x, where B and C are constant matrices. From these assumptions, R
(
x
) is a constant matrix, which we will denote simply by R
=
B
B
T, and Q
(
x
) is a quadratic form, which we will write as x
T
Q
x, where Q
=
C
T
C. Thus, the only nonlinearities present are in the system a
(
x
) matrix. For the purpose of the results to be presented, this nonlinear function will be approximated to fourth order. It should be noted that some of these terms may be zero, depending on the system considered. Therefore, we have the following approximation:
(21)
a
(
x
)
=
a
1
(
x
)
+
a
2
(
x
)
+
a
3
(
x
)
+
a
4
(
x
)
,
where
(22)
a
1
(
x
)
=
A
1
x
=
col
k
{
∑
n
A
1
(
k
,
n
)
x
n
}
,
a
2
(
x
)
=
col
k
{
x
T
A
2
k
x
}
=
col
k
{
∑
i
∑
j
x
i
x
j
A
2
k
(
i
,
j
)
}
,
a
3
(
x
)
=
ma
t
m
n
{
x
T
A
3
m
n
x
}
x
=
col
m
{
∑
n
(
∑
i
∑
j
x
i
x
j
A
3
m
n
(
i
,
j
)
)
x
n
}
,
a
4
(
x
)
=
col
k
{
x
T
ma
t
m
n
{
x
T
A
4
k
m
n
x
}
x
}
=
col
k
{
∑
m
∑
n
x
m
x
n
(
∑
i
∑
j
x
i
x
j
A
4
k
m
n
(
i
,
j
)
)
}
and the summations run from 0 to N. Here, A
1, A
2
k, A
3
m
n, and A
4
k
m
n are families of N
×
N square matrices. We shall also find it useful to define the column matrices a
2
i
j
=
col
k
{
A
2
k
(
i
,
j
)
} and a
3
n
j
k
=
col
m
{
A
3
m
n
(
j
,
k
)
} whose entries can be used to form A
2
k and A
3
m
n, respectively.Additionally, consider a storage functionV
(
x
) for which
(23)
∇
V
(
x
)
=
∇
V
1
(
x
)
+
∇
V
2
(
x
)
+
∇
V
3
(
x
)
+
∇
V
4
(
x
)
,
where
(24)
∇
V
1
(
x
)
=
x
T
P
1
=
ro
w
n
{
∑
m
x
m
P
1
(
m
,
n
)
}
,
∇
V
2
(
x
)
=
ro
w
k
{
x
T
P
2
k
x
}
=
ro
w
k
{
∑
i
∑
j
x
i
x
j
P
2
k
(
i
,
j
)
}
,
∇
V
3
(
x
)
=
x
T
ma
t
m
n
{
x
T
P
3
m
n
x
}
∇
V
3
(
x
)
=
ro
w
n
{
∑
m
‍
x
m
(
∑
i
‍
∑
j
‍
x
i
x
j
P
3
m
n
(
i
,
j
)
)
}
,
∇
V
4
(
x
)
=
ro
w
k
{
x
T
ma
t
m
n
{
x
T
P
4
k
m
n
x
}
x
}
∇
V
4
(
x
)
=
ro
w
k
{
∑
m
‍
∑
n
‍
x
m
x
n
(
∑
i
‍
∑
j
‍
x
i
x
j
P
k
m
n
(
i
,
j
)
)
}
.
Here, P
1, P
2
k, P
3
m
n, and P
4
k
m
n are families of N
×
N square matrices. In general, P
(
i
,
j
) refers to the (
i
,
j
)th entry of the matrix P. We shall also find it useful to define the column matrices p
2
i
j
=
col
k
{
P
2
k
(
i
,
j
)
} and p
3
m
i
j
=
col
n
{
P
3
m
n
(
i
,
j
)
} whose entries can be used to form P
2
k and P
3
m
n, respectively.It is important to recognize that, whilea
i
(
x
) may be zero for some i, the corresponding ∇
V
i
(
x
) is not necessarily zero as well. This is because the present method involves a Taylor series expansion of the nonlinear solution to the HJE. Substituting (21) and (23) into the HJE in (12) and grouping terms of the same order yields
(25)
𝒪
[
∥
x
∥
2
]
:
∇
V
1
(
x
)
a
1
(
x
)
-
1
2
∇
V
1
(
x
)
R
∇
V
1
T
(
x
)
+
1
2
x
T
Q
x
=
0
,
(26)
𝒪
[
∥
x
∥
k
]
:
∇
V
(
k
-
1
)
a
c
1
+
∇
V
1
a
(
k
-
1
)
+
∑
i
=
2
k
-
2
‍
∇
V
i
a
c
(
k
-
i
)
=
0
for k
≥
3, where a
c
1
=
a
1
(
x
)
-
R
∇
V
1
T
(
x
)
=
A
c
1
x with A
c
1
=
A
1
-
B
B
T
P
1 and a
c
i
=
a
i
(
x
)
-
(
1
/
2
)
R
∇
V
i
T
(
x
) for i
≥
2. Thus, at each order k, the objective is to solve for the unknown ∇
V
(
k
-
1
). Note that in (26) the summation term is equal to zero for k
=
3, since i
>
k
-
2.We now present the general expressions to compute the unknowns∇
V
(
k
-
1
) in (23). The first-order solution, P
1
=
P
1
T, is obtained by solving the ARE corresponding to (25), which is given by
(27)
P
1
A
1
+
A
1
T
P
1
-
P
1
B
B
T
P
1
+
C
T
C
=
0
,
where it has been noted that R
=
B
B
T and Q
=
C
T
C. We will assume that (
A
1
,
B
) is controllable and (C
,
A
1) is observable so that P
1 is positive definite. Then, the higher-order solutions are obtained by solving (26) recursively for increasing values of k.The second-order solution is given by(28)
p
2
i
j
T
=
-
a
2
i
j
T
P
1
A
c
1
-
1
,
i
,
j
=
1
,
…
,
N
,
where the column matrices p
2
i
j and a
2
i
j were defined above. We emphasize that the multiplications indicated in (28) are standard matrix multiplications.The third-order solution is given by(29)
ma
t
m
n
{
P
3
m
n
(
i
,
j
)
}
=
-
[
ma
t
m
n
{
p
2
m
i
T
a
c
2
j
n
}
+
P
1
ma
t
m
n
{
A
3
m
n
(
i
,
j
)
}
+
P
1
ma
t
m
n
{
A
3
m
n
(
i
,
j
)
}
]
A
c
1
-
1
,
i
,
j
=
1
,
…
,
N
,
where a
c
2
j
n
=
a
2
j
n
-
(
1
/
2
)
B
B
T
p
2
j
n. Here, ma
t
m
n
{
P
3
m
n
(
i
,
j
)
} consists of a square matrix (with row index m and column index n) containing the (
i
,
j
)th entries of each P
m
n for given i and j. Again, all of the indicated multiplications are standard matrix multiplications.The fourth-order solution is given by(30)
ma
t
m
n
{
P
4
n
m
k
(
i
,
j
)
}
=
-
[
ma
t
m
n
{
p
3
m
i
j
T
a
c
2
k
n
}
ma
t
m
n
{
A
4
n
m
k
(
i
,
j
)
}
P
1
+
ma
t
m
n
{
p
2
m
i
T
a
c
3
n
j
k
}
+
ma
t
m
n
{
A
4
n
m
k
(
i
,
j
)
}
P
1
]
A
c
1
-
1
,
i
,
j
,
k
=
1
,
…
,
N
,
where a
c
3
n
j
k
=
a
3
n
j
k
-
(
1
/
2
)
B
B
T
p
3
n
j
k. Similar comments on the multiplications involved above apply here.The negative-semidefinite solution of the HJE,V
-
(
x
), can be determined using (23)–(30) with the proviso that the matrix P
1 is replaced with the negative-definite solution of the Riccati equation in (27) which we will denote by Q
1 and A
c
1 is replaced by A
u
1
=
A
1
-
B
B
T
Q
1. The matrices corresponding to those in (24) will be denoted by Q
1, Q
2
k, Q
3
m
n, and Q
4
k
m
n and they are obtained by solving (27)–(30) with Q
1 replacing P
1, q
2
i
j replacing p
2
i
j, and so forth. In the next section, the satellite attitude dynamics are presented.


## 2.1. Class of Nonlinear Systems

Consider the nonlinear system shown in Figure1. The plant is given by
(1)
P
:
{
x
˙
(
t
)
=
a
(
x
)
+
b
(
x
)
u
1
(
t
)
,
y
1
(
t
)
=
c
(
x
)
,
where a
(
x
), b
(
x
), and c
(
x
) are assumed to be smooth (i.e., C
∞) nonlinear functions of the plant states with a
(
0
)
=
0 and c
(
0
)
=
0 (i.e., x
(
t
)
=
0 is an equilibrium). The controller is described by
(2)
K
:
{
x
˙
c
(
t
)
=
a
c
(
x
c
)
+
b
c
(
x
c
)
y
2
(
t
)
,
u
2
(
t
)
=
c
c
(
x
c
)
,
where a
c
(
x
c
), b
c
(
x
c
), and c
c
(
x
c
) are smooth nonlinear functions of the controller states. Additionally, the following relations hold:
(3)
u
1
(
t
)
=
u
0
(
t
)
-
u
2
(
t
)
,
y
2
(
t
)
=
y
0
(
t
)
-
y
1
(
t
)
.
In the above, x
(
t
)
∈
ℝ
N is the plant state vector, x
c
(
t
)
∈
ℝ
N
c is the controller state vector, u
2
(
t
)
∈
ℝ
N
u is the control signal, u
0
(
t
)
∈
ℝ
N
u is the exogenous (disturbance) input, u
1
(
t
)
∈
ℝ
N
u is the actual plant input, y
1
(
t
)
∈
ℝ
N
y is the actual plant output, y
0
(
t
)
∈
ℝ
N
y is the reference and/or sensor noise signal, and y
2
(
t
)
∈
ℝ
N
y is the tracking error.Figure 1
Block diagram of system (1) with controller (2).The plant in (1) can be written as
(4)
P
⋆
:
{
x
˙
(
t
)
=
a
(
x
)
+
[
b
(
x
)
0
]
w
(
t
)
+
[
-
b
(
x
)
]
u
(
t
)
,
[
z
1
(
t
)
z
2
(
t
)
]
=
[
0
-
c
(
x
)
]
+
[
0
0
0
1
]
w
(
t
)
+
[
1
0
]
u
(
t
)
,
y
(
t
)
=
[
-
c
(
x
)
]
+
[
0
1
]
w
(
t
)
,
where
(5)
w
(
t
)
=
[
u
0
(
t
)
y
0
(
t
)
]
,
z
(
t
)
=
[
z
1
(
t
)
z
2
(
t
)
]
=
[
u
2
(
t
)
y
2
(
t
)
]
,
u
(
t
)
=
u
2
(
t
)
,
y
(
t
)
=
y
2
(
t
)
.
The generalized system in (4) is shown in Figure 2. The first equation in (4) defines the plant dynamics with state variable x
(
t
), control input u
(
t
), and subject to a set of exogenous inputs w
(
t
), which includes disturbances (to be rejected), references (to be tracked), and/or noise (to be filtered). The second equation defines the penalty variable z
(
t
) representing the outputs of interest, which may include a tracking error, a function of some of the exogenous variables w
(
t
), and a cost of the input u
(
t
) needed to achieve the prescribed control goal. The third equation defines the set of measured variables y
(
t
), which are functions of the plant state x
(
t
) and the exogenous inputs w
(
t
). As we shall see next, we will be concerned with defining an upper bound on the ℒ
2-gain from the disturbance inputs w
(
t
) to the outputs z
(
t
) of the system in (4).Figure 2
Block diagram of generalized system (4) with controller (2).


## 2.2. Robustness in the Gap Metric

In general, the plant and controller are assumed to be causal mappings from their respective inputs and outputs; that is,P
:
𝒰
→
𝒴 and K
:
𝒴
→
𝒰, which satisfy P
0
=
0 and K
0
=
0, where 𝒰 and 𝒴 are appropriate signal spaces. In particular, these mappings can be represented by (1) and (2). Thus, in the feedback configuration of Figure 1, the signals u
i belong to 𝒰 and the signals y
i belong to 𝒴, where i
∈
{
0,1
,
2
}. The input-output relations describing the plant and controller can be represented by their respective graphs:
(6)
𝒢
𝒫
:
=
{
[
u
P
u
]
:
u
∈
𝒰
,
P
u
∈
𝒴
}
,
𝒢
𝒦
:
=
{
[
K
y
y
]
:
K
y
∈
𝒰
,
y
∈
𝒴
}
,
where 𝒢
𝒫, 𝒢
𝒦
⊂
𝒲
=
𝒰
×
𝒴. For a plant P
0 with ℳ
0
:
=
𝒢
𝒫
0 and a plant P
1 with ℳ
1
:
=
𝒢
𝒫
1, the ρ-gap between these two plants is defined as [9, 28]
(7)
ρ
g
(
P
0
,
P
1
)
=
max
⁡
⁡
{
ρ
g
→
(
P
0
,
P
1
)
,
ρ
g
→
(
P
1
,
P
0
)
}
,
where
(8)
ρ
g
→
(
P
0
,
P
1
)
=
sup
⁡
0
≠
v
0
∈
ℳ
0
⁡
inf
⁡
v
1
∈
ℳ
1
⁡
(
∥
v
0
-
v
1
∥
2
+
∥
P
0
v
0
-
P
1
v
1
∥
2
)
1
/
2
(
∥
v
0
∥
2
+
∥
P
0
v
0
∥
2
)
1
/
2
and ρ
g
→
(
P
1
,
P
0
) is similarly defined. The norm used here, ∥
(
·
)
∥, is defined in [9].It has been shown [9, 32] that, for a stable feedback pair (
P
0
,
K
), if a system P
1 is such that
(9)
ρ
g
(
P
0
,
P
1
)
<
∥
Π
𝒩
∥
ℳ
0
∥
∞
-
1
,
then the feedback system (
P
1
,
K
) is also stable. The symbol Π
𝒩
∥
ℳ defines a parallel projection operator [32] and represents the closed-loop mapping from the disturbances w
(
t
) to the outputs z
(
t
) shown in Figure 2. The ℒ
2-gain of Π
𝒩
∥
ℳ is the induced norm defined by
(10)
∥
Π
𝒩
∥
ℳ
∥
∞
=
sup
⁡
0
≠
w
∈
ℒ
2
⁡
∥
z
∥
2
∥
w
∥
2
,
where the ℒ
2-norm is defined by ∥
z
∥
2
=
∫
0
∞
z
T
z
d
t.Thus, to minimize the effects of the disturbancesw
(
t
) on the outputs z
(
t
), which simultaneously achieves optimal robustness, the objective is to design a controller K such that ∥
Π
𝒩
∥
ℳ
0
∥
∞ is minimized. However, no closed-form solution exists to this optimal control problem. Instead, we shall be concerned with the suboptimal problem: given some constant scalar γ, design K such that ∥
Π
𝒩
∥
ℳ
0
∥
∞
<
γ. Note that a γ-iteration procedure can be applied to obtain the optimal solution within an arbitrary tolerance. It will be shown that the controller that satisfies this objective can be obtained from the solution to a single Hamilton-Jacobi equation.


## 2.3. Nonlinearℋ
∞ Control Problem

As indicated in the previous subsection, we are concerned with rendering∥
Π
𝒩
∥
ℳ
0
∥
∞
<
γ. In other words, we wish to bound the ℒ
2-gain from the disturbance inputs w
(
t
) to the outputs z
(
t
) of the system in (4). Without providing the details, it is noted that the HJE corresponding to (4) is dependent on γ and the quadratic term is sign-indefinite (i.e., is neither positive definite nor negative definite). However, by performing a certain transformation [9, 33], the generalized system can be written in a form where the outputs z
(
t
) are not explicit functions of the disturbances w
(
t
). The resulting HJE is no longer dependent on γ and is sign-definite. The system resulting from this transformation is described by
(11)
P
¯
⋆
:
{
x
˙
(
t
)
=
a
(
x
)
+
[
-
b
(
x
)
0
]
w
(
t
)
+
[
-
b
(
x
)
]
+
u
(
t
)
.
[
z
1
(
t
)
z
2
(
t
)
]
=
[
0
-
β
-
1
c
(
x
)
]
+
[
1
0
]
u
(
t
)
,
y
(
t
)
=
[
-
β
-
1
c
(
x
)
]
+
[
0
-
1
]
w
(
t
)
,
where β
=
1
-
γ
-
2. Now, the HJE corresponding to the system in (11) is given by
(12)
∇
V
(
x
)
A
(
x
)
-
1
2
∇
V
(
x
)
R
(
x
)
∇
V
T
(
x
)
+
1
2
Q
(
x
)
=
0
,
where ∇
V
(
x
) is the Jacobian matrix of the storage function V
(
x
) and is defined as
(13)
∇
V
(
x
)
=
∂
V
∂
x
T
=
row
i
{
∂
V
∂
x
i
}
with ro
w
i
{
·
} denoting a row matrix in the index i
∈
{
1,2
,
…
,
N
}. For the system in (11), we have the parameters [9]
(14)
A
(
x
)
=
a
(
x
)
,
R
(
x
)
=
b
(
x
)
b
T
(
x
)
,
Q
(
x
)
=
c
T
(
x
)
c
(
x
)
.
The name used in the literature for the HJE in (12) is not uniform. In the context of robust control, it is sometimes referred to as the Hamilton-Jacobi-Isaacs equation and when used in nonlinear optimal control it is often referred to as the Hamilton-Jacobi-Bellman equation. In [9] it is referred to as the Hamilton-Jacobi-Bellman-Isaacs equation. We will simply call it the HJE.It is shown in [9] that the controller K that solves the suboptimal ℋ
∞ problem for the plant P
⋆ (i.e., renders ∥
Π
𝒩
∥
ℳ
0
∥
∞
<
γ) is related to that (call it K
¯) which solves the same problem for the modified plant P
¯
⋆ as follows: K
=
K
¯
β. We now construct this K. Denote by V
+
(
x
) the unique smooth solution to the HJE (12) that satisfies V
+
(
x
)
≥
0 and V
+
(
0
)
=
0, with
(15)
+
[
A
(
x
)
-
R
(
x
)
∇
V
+
T
(
x
)
]
asymptotically stable. Similarly, we denote by V
-
(
x
) the unique smooth solution to the HJE that satisfies V
-
(
x
)
≤
0 and V
-
(
0
)
=
0, with
(16)
-
[
A
(
x
)
-
R
(
x
)
∇
V
-
T
(
x
)
]
asymptotically stable.Following the approach of James et al. [9], a local solution to the nonlinear ℋ
∞ control problem for the plant in (4) is obtained if the following conditions are satisfied.(1)
There exists aC
3 positive-definite function V
+
(
x
) defined in a neighbourhood of the origin with V
+
(
0
)
=
0 that satisfies the HJE (12).
(2)
There exists aC
3 negative-definite function V
-
(
x
) defined in a neighbourhood of the origin with V
-
(
0
)
=
0 that satisfies the HJE (12), and additionally satisfies
(17)
∇
2
{
-
∇
r
2
e
(
A
+
γ
-
2
β
-
2
R
∇
V
+
T
)
-
1
2
γ
2
β
-
2
+
1
2
γ
-
2
∇
r
2
e
R
∇
r
2
e
T
+
1
2
β
-
4
∇
V
+
R
∇
V
+
T
-
1
2
γ
2
β
-
2
Q
}
<
0
,
where
(18)
r
2
e
(
x
)
=
β
-
2
V
+
(
x
)
+
γ
2
V
-
(
x
)
.
(3)
There exists aC
2 function G
2
(
x
) such that
(19)
∇
r
2
e
(
x
)
G
2
(
x
)
=
γ
2
β
-
1
c
T
(
x
)
in a neighbourhood of the origin.The resulting nonlinear controller is given by(20)
K
:
{
x
˙
c
=
a
(
x
c
)
-
b
(
x
c
)
b
T
(
x
c
)
∇
V
+
T
(
x
c
)
x
˙
c
=
+
G
2
(
x
c
)
[
β
-
1
c
(
x
c
)
+
β
y
2
]
,
u
2
(
x
c
)
=
β
-
2
b
T
(
x
c
)
∇
V
+
T
(
x
c
)
.In the next subsection, we shall examine how to obtain analytical expressions for the approximate solution to the HJE in (12) required for this control law.


## 2.4. Analytical Solutions to HJE Approximation

Our approach for the synthesis of nonlinear controllers is based on the Taylor series approximation of the solution to the Hamilton-Jacobi equation, where each order of the controller is built using the previous orders. The following notation will be used:ro
w
k
{
·
} denotes a row matrix with index k, col
k
{
·
} denotes a column matrix with index k, and ma
t
m
n
{
·
} denotes a matrix with row index m and column index n. It should be emphasized that the symbols k, m, and n used here are dummy indices. In general, A
(
i
,
j
) refers to the (
i
,
j
)th entry of the matrix A. We will denote the positive-semidefinite solution V
+
(
x
) of the HJE (12) simply by V
(
x
).Consider the nonlinear system in (1). We begin by making the assumptions that b
(
x
)
=
B and c
(
x
)
=
C
x, where B and C are constant matrices. From these assumptions, R
(
x
) is a constant matrix, which we will denote simply by R
=
B
B
T, and Q
(
x
) is a quadratic form, which we will write as x
T
Q
x, where Q
=
C
T
C. Thus, the only nonlinearities present are in the system a
(
x
) matrix. For the purpose of the results to be presented, this nonlinear function will be approximated to fourth order. It should be noted that some of these terms may be zero, depending on the system considered. Therefore, we have the following approximation:
(21)
a
(
x
)
=
a
1
(
x
)
+
a
2
(
x
)
+
a
3
(
x
)
+
a
4
(
x
)
,
where
(22)
a
1
(
x
)
=
A
1
x
=
col
k
{
∑
n
A
1
(
k
,
n
)
x
n
}
,
a
2
(
x
)
=
col
k
{
x
T
A
2
k
x
}
=
col
k
{
∑
i
∑
j
x
i
x
j
A
2
k
(
i
,
j
)
}
,
a
3
(
x
)
=
ma
t
m
n
{
x
T
A
3
m
n
x
}
x
=
col
m
{
∑
n
(
∑
i
∑
j
x
i
x
j
A
3
m
n
(
i
,
j
)
)
x
n
}
,
a
4
(
x
)
=
col
k
{
x
T
ma
t
m
n
{
x
T
A
4
k
m
n
x
}
x
}
=
col
k
{
∑
m
∑
n
x
m
x
n
(
∑
i
∑
j
x
i
x
j
A
4
k
m
n
(
i
,
j
)
)
}
and the summations run from 0 to N. Here, A
1, A
2
k, A
3
m
n, and A
4
k
m
n are families of N
×
N square matrices. We shall also find it useful to define the column matrices a
2
i
j
=
col
k
{
A
2
k
(
i
,
j
)
} and a
3
n
j
k
=
col
m
{
A
3
m
n
(
j
,
k
)
} whose entries can be used to form A
2
k and A
3
m
n, respectively.Additionally, consider a storage functionV
(
x
) for which
(23)
∇
V
(
x
)
=
∇
V
1
(
x
)
+
∇
V
2
(
x
)
+
∇
V
3
(
x
)
+
∇
V
4
(
x
)
,
where
(24)
∇
V
1
(
x
)
=
x
T
P
1
=
ro
w
n
{
∑
m
x
m
P
1
(
m
,
n
)
}
,
∇
V
2
(
x
)
=
ro
w
k
{
x
T
P
2
k
x
}
=
ro
w
k
{
∑
i
∑
j
x
i
x
j
P
2
k
(
i
,
j
)
}
,
∇
V
3
(
x
)
=
x
T
ma
t
m
n
{
x
T
P
3
m
n
x
}
∇
V
3
(
x
)
=
ro
w
n
{
∑
m
‍
x
m
(
∑
i
‍
∑
j
‍
x
i
x
j
P
3
m
n
(
i
,
j
)
)
}
,
∇
V
4
(
x
)
=
ro
w
k
{
x
T
ma
t
m
n
{
x
T
P
4
k
m
n
x
}
x
}
∇
V
4
(
x
)
=
ro
w
k
{
∑
m
‍
∑
n
‍
x
m
x
n
(
∑
i
‍
∑
j
‍
x
i
x
j
P
k
m
n
(
i
,
j
)
)
}
.
Here, P
1, P
2
k, P
3
m
n, and P
4
k
m
n are families of N
×
N square matrices. In general, P
(
i
,
j
) refers to the (
i
,
j
)th entry of the matrix P. We shall also find it useful to define the column matrices p
2
i
j
=
col
k
{
P
2
k
(
i
,
j
)
} and p
3
m
i
j
=
col
n
{
P
3
m
n
(
i
,
j
)
} whose entries can be used to form P
2
k and P
3
m
n, respectively.It is important to recognize that, whilea
i
(
x
) may be zero for some i, the corresponding ∇
V
i
(
x
) is not necessarily zero as well. This is because the present method involves a Taylor series expansion of the nonlinear solution to the HJE. Substituting (21) and (23) into the HJE in (12) and grouping terms of the same order yields
(25)
𝒪
[
∥
x
∥
2
]
:
∇
V
1
(
x
)
a
1
(
x
)
-
1
2
∇
V
1
(
x
)
R
∇
V
1
T
(
x
)
+
1
2
x
T
Q
x
=
0
,
(26)
𝒪
[
∥
x
∥
k
]
:
∇
V
(
k
-
1
)
a
c
1
+
∇
V
1
a
(
k
-
1
)
+
∑
i
=
2
k
-
2
‍
∇
V
i
a
c
(
k
-
i
)
=
0
for k
≥
3, where a
c
1
=
a
1
(
x
)
-
R
∇
V
1
T
(
x
)
=
A
c
1
x with A
c
1
=
A
1
-
B
B
T
P
1 and a
c
i
=
a
i
(
x
)
-
(
1
/
2
)
R
∇
V
i
T
(
x
) for i
≥
2. Thus, at each order k, the objective is to solve for the unknown ∇
V
(
k
-
1
). Note that in (26) the summation term is equal to zero for k
=
3, since i
>
k
-
2.We now present the general expressions to compute the unknowns∇
V
(
k
-
1
) in (23). The first-order solution, P
1
=
P
1
T, is obtained by solving the ARE corresponding to (25), which is given by
(27)
P
1
A
1
+
A
1
T
P
1
-
P
1
B
B
T
P
1
+
C
T
C
=
0
,
where it has been noted that R
=
B
B
T and Q
=
C
T
C. We will assume that (
A
1
,
B
) is controllable and (C
,
A
1) is observable so that P
1 is positive definite. Then, the higher-order solutions are obtained by solving (26) recursively for increasing values of k.The second-order solution is given by(28)
p
2
i
j
T
=
-
a
2
i
j
T
P
1
A
c
1
-
1
,
i
,
j
=
1
,
…
,
N
,
where the column matrices p
2
i
j and a
2
i
j were defined above. We emphasize that the multiplications indicated in (28) are standard matrix multiplications.The third-order solution is given by(29)
ma
t
m
n
{
P
3
m
n
(
i
,
j
)
}
=
-
[
ma
t
m
n
{
p
2
m
i
T
a
c
2
j
n
}
+
P
1
ma
t
m
n
{
A
3
m
n
(
i
,
j
)
}
+
P
1
ma
t
m
n
{
A
3
m
n
(
i
,
j
)
}
]
A
c
1
-
1
,
i
,
j
=
1
,
…
,
N
,
where a
c
2
j
n
=
a
2
j
n
-
(
1
/
2
)
B
B
T
p
2
j
n. Here, ma
t
m
n
{
P
3
m
n
(
i
,
j
)
} consists of a square matrix (with row index m and column index n) containing the (
i
,
j
)th entries of each P
m
n for given i and j. Again, all of the indicated multiplications are standard matrix multiplications.The fourth-order solution is given by(30)
ma
t
m
n
{
P
4
n
m
k
(
i
,
j
)
}
=
-
[
ma
t
m
n
{
p
3
m
i
j
T
a
c
2
k
n
}
ma
t
m
n
{
A
4
n
m
k
(
i
,
j
)
}
P
1
+
ma
t
m
n
{
p
2
m
i
T
a
c
3
n
j
k
}
+
ma
t
m
n
{
A
4
n
m
k
(
i
,
j
)
}
P
1
]
A
c
1
-
1
,
i
,
j
,
k
=
1
,
…
,
N
,
where a
c
3
n
j
k
=
a
3
n
j
k
-
(
1
/
2
)
B
B
T
p
3
n
j
k. Similar comments on the multiplications involved above apply here.The negative-semidefinite solution of the HJE,V
-
(
x
), can be determined using (23)–(30) with the proviso that the matrix P
1 is replaced with the negative-definite solution of the Riccati equation in (27) which we will denote by Q
1 and A
c
1 is replaced by A
u
1
=
A
1
-
B
B
T
Q
1. The matrices corresponding to those in (24) will be denoted by Q
1, Q
2
k, Q
3
m
n, and Q
4
k
m
n and they are obtained by solving (27)–(30) with Q
1 replacing P
1, q
2
i
j replacing p
2
i
j, and so forth. In the next section, the satellite attitude dynamics are presented.


## 3. Attitude Dynamics and Control

### 3.1. The Attitude Dynamics and Kinematics

The attitude dynamics of a rigid spacecraft are given by Euler's equation:(31)
ω
˙
=
-
I
-
1
ω
×
I
ω
+
I
-
1
u
,
where ω
=
[
ω
1
ω
2
ω
3
]
T are the body angular velocities, I is the moment of inertia matrix, and u
1
=
[
u
1
u
2
u
3
]
T are the body torques. The notation ω
× is the matrix representation of the cross product and is defined as
(32)
ω
×
=
[
0
-
ω
3
ω
2
ω
3
0
-
ω
1
-
ω
2
ω
1
0
]
.While many representations are possible to define the spacecraft attitude kinematics, the modified Rodrigues parameters (MRPs) are chosen here because they are polynomial in the states, which fit nicely with the present controller synthesis approach, and they possess neither singularities nor norm constraints when used in conjunction with the shadow parameters. The MRP vectorσ
=
[
σ
1
σ
2
σ
3
]
T can be defined in terms of the principal rotation axis e
^
=
[
e
1
e
2
e
3
]
T and principal rotation angle Φ of Euler's theorem according to
(33)
σ
=
e
^
tan
⁡
(
Φ
4
)
.
The attitude kinematics using MRPs are defined by
(34)
σ
˙
=
1
2
[
1
2
(
1
-
σ
T
σ
)
1
+
σ
×
+
σ
σ
T
]
ω
,
where 1 is the identity matrix.Upon closer inspection of (33), it is seen that the MRPs encounter a singularity for rotations of Φ
=
±
2
π rad. This corresponds to a complete rotation in either direction about the principal axis. To circumvent this, another set of MRPs, called the shadow parameters and denoted by σ
S, is used in conjunction with the regular MRPs. By switching from one set to the other at rotations of Φ
=
±
π rad, it is possible to avoid any singularities. The parameter switching occurs on the surface defined by
(35)
σ
T
σ
=
σ
S
T
σ
S
=
1
.
The kinematics are identical for both the regular and the shadow parameters. However, when the switching surface is encountered, both the MRPs and their rates must be converted from one set to the other. This can be accomplished with the following relations:
(36)
σ
˙
S
=
-
σ
˙
(
σ
T
σ
)
+
1
+
(
σ
T
σ
)
2
(
σ
T
σ
)
2
σ
σ
T
ω
,
σ
S
=
-
σ
(
σ
T
σ
)
.
More details can be found in the text by Schaub and Junkins [34].Defining the state vectorx
=
[
ω
1
ω
2
ω
3
σ
1
σ
2
σ
3
]
T and grouping terms of the same order, the attitude dynamics in (31) and kinematics in (34) can be expressed as
(37)
x
˙
(
t
)
=
a
1
(
x
)
+
a
2
(
x
)
+
a
3
(
x
)
+
B
u
1
(
t
)
,
y
1
(
t
)
=
C
x
,
where
(38)
a
1
(
x
)
=
A
1
x
=
[
0
0
1
4
1
0
]
[
ω
σ
]
,
B
=
[
I
-
1
0
]
,
C
=
1
,
and the second- and third-order terms are given by
(39)
a
2
(
x
)
=
[
-
I
-
1
ω
×
I
ω
1
2
σ
×
ω
]
,
a
3
(
x
)
=
[
0
(
1
2
σ
σ
T
-
1
4
σ
T
σ
1
)
ω
]
,
respectively. It should be noted that these third-order attitude dynamics are exact so that we may take a
4
(
x
)
=
0 and hence A
4
k
m
n
=
0.


### 3.2. The Attitude Controller

In this section, the proposed controller synthesis methods will be applied to the satellite attitude dynamics given by (37)–(39). For simplicity, we shall assume that the dynamics are formulated in principal axes so that the inertia matrix is given by I
=
diag
⁡
{
I
1
,
I
2
,
I
3
}. Let us begin by comparing the definitions of a
2 and a
3 in (22) with the specific ones given in (39). From this, we identify the nonzero elements of the matrices A
2
k and A
3
m
n as follows:
(40)
A
2,1
23
=
A
2,1
32
=
k
1
2
,
A
2,2
13
=
A
2,2
31
=
k
2
2
,
A
2,3
12
=
A
2,3
21
=
k
3
2
,
A
2,4
53
=
-
A
2,4
62
=
-
A
2,5
43
=
A
2,5
61
=
A
2,6
43
=
-
A
2,6
51
=
1
2
,
where k
1
=
(
I
2
-
I
3
)
/
I
1, k
2
=
(
I
3
-
I
1
)
/
I
2 and k
3
=
(
I
1
-
I
2
)
/
I
3. In addition,
(41)
2
A
3,41
44
=
A
3,41
55
=
A
3,41
66
=
2
A
3,42
45
2
A
3,41
44
=
2
A
3,42
54
=
2
A
3,43
64
=
2
A
3,43
46
=
1
2
,
2
A
3,51
54
=
2
A
3,51
45
=
A
3,52
44
=
2
A
3,52
55
2
A
3,41
44
=
A
3,52
66
=
2
A
3,53
56
=
2
A
3,53
65
=
1
2
,
2
A
3,61
64
=
2
A
3,61
46
=
2
A
3,62
65
=
2
A
3,56
45
2
A
3,41
44
=
A
3,63
44
=
A
3,63
55
=
2
A
3,63
66
=
1
2
.Given the definitions ofA
1, B, and C in (38), the positive-definite solution of the algebraic Riccati equation in (27) is easily determined. The nonzero elements are given by
(42)
P
1
(
i
,
i
)
=
I
i
1
+
(
I
i
2
)
,
P
1
(
i
,
3
+
i
)
=
P
1
(
3
+
i
,
i
)
=
I
i
,
P
1
(
3
+
i
,
3
+
i
)
=
4
1
+
(
I
i
2
)
,
i
=
1,2
,
3
.
The negative-definite solution (nonzero elements) is given by
(43)
Q
1
(
i
,
i
)
=
-
I
i
1
+
(
I
i
2
)
,
Q
1
(
i
,
3
+
i
)
=
Q
1
(
3
+
i
,
i
)
=
I
i
,
Q
1
(
3
+
i
,
3
+
i
)
=
-
4
1
+
(
I
i
2
)
,
i
=
1,2
,
3
.
The corresponding closed-loop matrices A
c
1
=
A
1
-
B
B
T
P
1 and A
u
1
=
A
1
-
B
B
T
Q
1 are readily determined:
(44)
A
c
1
(
i
,
i
)
=
-
I
i
-
1
1
+
(
I
i
2
)
,
A
c
1
(
i
,
3
+
i
)
=
I
i
-
1
,
A
c
1
(
3
+
i
,
i
)
=
1
4
,
A
c
1
(
3
+
i
,
3
+
i
)
=
0
,
i
=
1,2
,
3
,
A
u
1
(
i
,
i
)
=
I
i
-
1
1
+
(
I
i
2
)
,
A
u
1
(
i
,
3
+
i
)
=
I
i
-
1
,
A
u
1
(
3
+
i
,
i
)
=
1
4
,
A
u
1
(
3
+
i
,
3
+
i
)
=
0
,
i
=
1,2
,
3
.Using the above quantities, the entries inP
2
k (via p
2
i
j), P
3
m
n, and P
4
k
n
m can be calculated using (28), (29), and (30). The same equations can be used to determine Q
2
k (via q
2
i
j), Q
3
m
n, and Q
4
k
n
m with P
1 replaced with Q
1 and A
c
1 replaced with A
u
1. The dynamic controller in (20) can be made specific to the attitude control problem:
(45)
x
˙
c
=
a
(
x
c
)
-
B
B
T
∇
V
+
T
(
x
c
)
+
G
2
(
x
c
)
[
β
-
1
C
x
c
+
β
y
2
]
,
u
2
(
x
c
)
=
β
-
2
B
T
∇
V
+
T
(
x
c
)
,
where a
(
x
c
) is determined using (22) in conjunction with (40) and (41) and ∇
V
+
T
(
x
c
) is determined using (24) and the solutions in (27)–(30). The observer gain G
c
(
x
c
) which is defined by (19) can be determined as follows. Since, r
2
e
(
x
)
=
β
-
2
V
+
(
x
)
+
γ
2
V
-
(
x
), we have used (24):
(46)
∇
r
2
e
(
x
)
=
x
T
R
e
2
(
x
)
=
x
T
[
(
β
-
2
P
1
+
γ
2
Q
1
)
ro
w
k
(
β
-
2
P
2
k
x
+
γ
2
Q
2
k
x
)
+
⋯
=
x
T
+
ro
w
k
(
β
-
2
P
2
k
x
+
γ
2
Q
2
k
x
)
+
⋯
]
.
Using this in (19) yields the following expression for the observer gain:
(47)
G
(
x
c
)
=
γ
2
β
-
1
[
(
β
-
2
P
1
+
γ
2
Q
1
)
+
row
k
(
β
-
2
P
2
k
x
c
+
γ
2
Q
2
k
x
c
)
+
⋯
]
-
1
C
T
.
The condition in (17) needs to be satisfied in a region of the origin. Using the lowest order terms in the expansions for A
=
a
(
x
), ∇
V
+, ∇
V
-, and ∇
r
2
e, the Hessian matrix defined by (17) is given by
(48)
-
(
β
-
2
P
1
+
γ
2
Q
1
)
(
A
1
+
γ
-
2
β
-
2
B
B
T
P
1
)
-
(
A
1
+
γ
-
2
β
-
2
B
B
T
P
1
)
T
(
β
-
2
P
1
+
γ
2
Q
1
)
+
γ
-
2
(
β
-
2
P
1
+
γ
2
Q
1
)
B
B
T
(
β
-
2
P
1
+
γ
2
Q
1
)
+
β
-
4
P
1
B
B
T
P
1
-
γ
2
β
-
2
C
T
C
which must be negative definite. This condition limits the chosen value of γ. In the sequel, we shall refer to the controller in (45) as the ℋ
∞ controller of order k, where k is the order of the approximation adopted for a
(
x
), ∇
V
+, ∇
V
-, and ∇
r
2
e in determining G
2.


## 3.1. The Attitude Dynamics and Kinematics

The attitude dynamics of a rigid spacecraft are given by Euler's equation:(31)
ω
˙
=
-
I
-
1
ω
×
I
ω
+
I
-
1
u
,
where ω
=
[
ω
1
ω
2
ω
3
]
T are the body angular velocities, I is the moment of inertia matrix, and u
1
=
[
u
1
u
2
u
3
]
T are the body torques. The notation ω
× is the matrix representation of the cross product and is defined as
(32)
ω
×
=
[
0
-
ω
3
ω
2
ω
3
0
-
ω
1
-
ω
2
ω
1
0
]
.While many representations are possible to define the spacecraft attitude kinematics, the modified Rodrigues parameters (MRPs) are chosen here because they are polynomial in the states, which fit nicely with the present controller synthesis approach, and they possess neither singularities nor norm constraints when used in conjunction with the shadow parameters. The MRP vectorσ
=
[
σ
1
σ
2
σ
3
]
T can be defined in terms of the principal rotation axis e
^
=
[
e
1
e
2
e
3
]
T and principal rotation angle Φ of Euler's theorem according to
(33)
σ
=
e
^
tan
⁡
(
Φ
4
)
.
The attitude kinematics using MRPs are defined by
(34)
σ
˙
=
1
2
[
1
2
(
1
-
σ
T
σ
)
1
+
σ
×
+
σ
σ
T
]
ω
,
where 1 is the identity matrix.Upon closer inspection of (33), it is seen that the MRPs encounter a singularity for rotations of Φ
=
±
2
π rad. This corresponds to a complete rotation in either direction about the principal axis. To circumvent this, another set of MRPs, called the shadow parameters and denoted by σ
S, is used in conjunction with the regular MRPs. By switching from one set to the other at rotations of Φ
=
±
π rad, it is possible to avoid any singularities. The parameter switching occurs on the surface defined by
(35)
σ
T
σ
=
σ
S
T
σ
S
=
1
.
The kinematics are identical for both the regular and the shadow parameters. However, when the switching surface is encountered, both the MRPs and their rates must be converted from one set to the other. This can be accomplished with the following relations:
(36)
σ
˙
S
=
-
σ
˙
(
σ
T
σ
)
+
1
+
(
σ
T
σ
)
2
(
σ
T
σ
)
2
σ
σ
T
ω
,
σ
S
=
-
σ
(
σ
T
σ
)
.
More details can be found in the text by Schaub and Junkins [34].Defining the state vectorx
=
[
ω
1
ω
2
ω
3
σ
1
σ
2
σ
3
]
T and grouping terms of the same order, the attitude dynamics in (31) and kinematics in (34) can be expressed as
(37)
x
˙
(
t
)
=
a
1
(
x
)
+
a
2
(
x
)
+
a
3
(
x
)
+
B
u
1
(
t
)
,
y
1
(
t
)
=
C
x
,
where
(38)
a
1
(
x
)
=
A
1
x
=
[
0
0
1
4
1
0
]
[
ω
σ
]
,
B
=
[
I
-
1
0
]
,
C
=
1
,
and the second- and third-order terms are given by
(39)
a
2
(
x
)
=
[
-
I
-
1
ω
×
I
ω
1
2
σ
×
ω
]
,
a
3
(
x
)
=
[
0
(
1
2
σ
σ
T
-
1
4
σ
T
σ
1
)
ω
]
,
respectively. It should be noted that these third-order attitude dynamics are exact so that we may take a
4
(
x
)
=
0 and hence A
4
k
m
n
=
0.


## 3.2. The Attitude Controller

In this section, the proposed controller synthesis methods will be applied to the satellite attitude dynamics given by (37)–(39). For simplicity, we shall assume that the dynamics are formulated in principal axes so that the inertia matrix is given by I
=
diag
⁡
{
I
1
,
I
2
,
I
3
}. Let us begin by comparing the definitions of a
2 and a
3 in (22) with the specific ones given in (39). From this, we identify the nonzero elements of the matrices A
2
k and A
3
m
n as follows:
(40)
A
2,1
23
=
A
2,1
32
=
k
1
2
,
A
2,2
13
=
A
2,2
31
=
k
2
2
,
A
2,3
12
=
A
2,3
21
=
k
3
2
,
A
2,4
53
=
-
A
2,4
62
=
-
A
2,5
43
=
A
2,5
61
=
A
2,6
43
=
-
A
2,6
51
=
1
2
,
where k
1
=
(
I
2
-
I
3
)
/
I
1, k
2
=
(
I
3
-
I
1
)
/
I
2 and k
3
=
(
I
1
-
I
2
)
/
I
3. In addition,
(41)
2
A
3,41
44
=
A
3,41
55
=
A
3,41
66
=
2
A
3,42
45
2
A
3,41
44
=
2
A
3,42
54
=
2
A
3,43
64
=
2
A
3,43
46
=
1
2
,
2
A
3,51
54
=
2
A
3,51
45
=
A
3,52
44
=
2
A
3,52
55
2
A
3,41
44
=
A
3,52
66
=
2
A
3,53
56
=
2
A
3,53
65
=
1
2
,
2
A
3,61
64
=
2
A
3,61
46
=
2
A
3,62
65
=
2
A
3,56
45
2
A
3,41
44
=
A
3,63
44
=
A
3,63
55
=
2
A
3,63
66
=
1
2
.Given the definitions ofA
1, B, and C in (38), the positive-definite solution of the algebraic Riccati equation in (27) is easily determined. The nonzero elements are given by
(42)
P
1
(
i
,
i
)
=
I
i
1
+
(
I
i
2
)
,
P
1
(
i
,
3
+
i
)
=
P
1
(
3
+
i
,
i
)
=
I
i
,
P
1
(
3
+
i
,
3
+
i
)
=
4
1
+
(
I
i
2
)
,
i
=
1,2
,
3
.
The negative-definite solution (nonzero elements) is given by
(43)
Q
1
(
i
,
i
)
=
-
I
i
1
+
(
I
i
2
)
,
Q
1
(
i
,
3
+
i
)
=
Q
1
(
3
+
i
,
i
)
=
I
i
,
Q
1
(
3
+
i
,
3
+
i
)
=
-
4
1
+
(
I
i
2
)
,
i
=
1,2
,
3
.
The corresponding closed-loop matrices A
c
1
=
A
1
-
B
B
T
P
1 and A
u
1
=
A
1
-
B
B
T
Q
1 are readily determined:
(44)
A
c
1
(
i
,
i
)
=
-
I
i
-
1
1
+
(
I
i
2
)
,
A
c
1
(
i
,
3
+
i
)
=
I
i
-
1
,
A
c
1
(
3
+
i
,
i
)
=
1
4
,
A
c
1
(
3
+
i
,
3
+
i
)
=
0
,
i
=
1,2
,
3
,
A
u
1
(
i
,
i
)
=
I
i
-
1
1
+
(
I
i
2
)
,
A
u
1
(
i
,
3
+
i
)
=
I
i
-
1
,
A
u
1
(
3
+
i
,
i
)
=
1
4
,
A
u
1
(
3
+
i
,
3
+
i
)
=
0
,
i
=
1,2
,
3
.Using the above quantities, the entries inP
2
k (via p
2
i
j), P
3
m
n, and P
4
k
n
m can be calculated using (28), (29), and (30). The same equations can be used to determine Q
2
k (via q
2
i
j), Q
3
m
n, and Q
4
k
n
m with P
1 replaced with Q
1 and A
c
1 replaced with A
u
1. The dynamic controller in (20) can be made specific to the attitude control problem:
(45)
x
˙
c
=
a
(
x
c
)
-
B
B
T
∇
V
+
T
(
x
c
)
+
G
2
(
x
c
)
[
β
-
1
C
x
c
+
β
y
2
]
,
u
2
(
x
c
)
=
β
-
2
B
T
∇
V
+
T
(
x
c
)
,
where a
(
x
c
) is determined using (22) in conjunction with (40) and (41) and ∇
V
+
T
(
x
c
) is determined using (24) and the solutions in (27)–(30). The observer gain G
c
(
x
c
) which is defined by (19) can be determined as follows. Since, r
2
e
(
x
)
=
β
-
2
V
+
(
x
)
+
γ
2
V
-
(
x
), we have used (24):
(46)
∇
r
2
e
(
x
)
=
x
T
R
e
2
(
x
)
=
x
T
[
(
β
-
2
P
1
+
γ
2
Q
1
)
ro
w
k
(
β
-
2
P
2
k
x
+
γ
2
Q
2
k
x
)
+
⋯
=
x
T
+
ro
w
k
(
β
-
2
P
2
k
x
+
γ
2
Q
2
k
x
)
+
⋯
]
.
Using this in (19) yields the following expression for the observer gain:
(47)
G
(
x
c
)
=
γ
2
β
-
1
[
(
β
-
2
P
1
+
γ
2
Q
1
)
+
row
k
(
β
-
2
P
2
k
x
c
+
γ
2
Q
2
k
x
c
)
+
⋯
]
-
1
C
T
.
The condition in (17) needs to be satisfied in a region of the origin. Using the lowest order terms in the expansions for A
=
a
(
x
), ∇
V
+, ∇
V
-, and ∇
r
2
e, the Hessian matrix defined by (17) is given by
(48)
-
(
β
-
2
P
1
+
γ
2
Q
1
)
(
A
1
+
γ
-
2
β
-
2
B
B
T
P
1
)
-
(
A
1
+
γ
-
2
β
-
2
B
B
T
P
1
)
T
(
β
-
2
P
1
+
γ
2
Q
1
)
+
γ
-
2
(
β
-
2
P
1
+
γ
2
Q
1
)
B
B
T
(
β
-
2
P
1
+
γ
2
Q
1
)
+
β
-
4
P
1
B
B
T
P
1
-
γ
2
β
-
2
C
T
C
which must be negative definite. This condition limits the chosen value of γ. In the sequel, we shall refer to the controller in (45) as the ℋ
∞ controller of order k, where k is the order of the approximation adopted for a
(
x
), ∇
V
+, ∇
V
-, and ∇
r
2
e in determining G
2.


## 4. Numerical Example and Comparisons

In this section, the nonlinear controller presented above is compared through numerical simulations with existing methods from the literature for spacecraft attitude regulation. The purpose of these comparisons is to examine the effects of different disturbances and uncertainties on the performance and robustness of the various controllers. In particular, we include gravity-gradient and geomagnetic torques, as well as unmodeled actuator dynamics and actuation time delays. In addition to these comparisons, we make use of the gap metric [28] to characterize the difference in the input-output (IO) map of the system induced by the unmodeled actuator dynamics.The simulation parameters are as follows. The satellite is in a circular orbit with an inclination ofi
=
87 degrees and a longitude of the ascending node of Ω
=
0. The initial value of the argument of latitude is zero. The altitude is 700 km and we take R
e
=
6378.14 km for the Earth's radius. These orbital parameters will be used to determine the gravity-gradient and geomagnetic disturbance torques acting on the spacecraft. The spacecraft position in the geocentric inertial frame, r
i
(
t
), is determined using a simple Keplerian model. The gravity-gradient torque is then given by u
0
g
=
3
(
μ
/
R
5
)
r
b
×
I
r
b, where μ is the geocentric gravitational parameter, R
=
r
i
T
r
i, and r
b
=
C
b
i
r
i, where the rotation matrix relating the body-fixed frame to the inertial frame may be expressed in terms of the MRPs as C
b
i
=
(
1
-
σ
×
)
2
(
1
+
σ
×
)
2.For the purposes of the geomagnetic disturbance torque, the satellite is assumed to generate a magnetic dipole ofm
=
[
0.1
0.1
0.1
]
T A·m2. The magnetic field model is the tilted dipole model, B
i
(
r
i
), presented in [35], where B
i are the geomagnetic field components expressed in the geocentric inertial frame. The geomagnetic disturbance torque is given by u
0
m
=
m
×
C
b
i
B
i. The satellite inertia matrix is given by I
=
diag
⁡
{
10.0
,
6.3
,
8.5
} kg·m2. In all comparisons, we consider the regulation problem only, hence y
0
=
0 and the input to the controller is y
2
=
-
y
1
=
-
C
x
=
-
x. Hence, we assume perfect measurements of the state.In the case of the disturbance rejection comparisons, the simulations start from the desired attitude and we compare the ability of the different controllers to maintain that attitude. For all other comparisons, the initial states are chosen from Schaub et al. [2] to be ω
(
0
)
=
[
1.4
0.9
0.8
]
T rad/s and σ
(
0
)
=
[
0.87
0
0
]
T. These initial conditions have the satellite oriented almost π rad from the desired attitude with large angular velocities moving it towards this upside-down attitude. All simulations will be performed for one orbit using a 4th-order Runge-Kutta numerical integration method with a step-size Δ
t
=
0.01 s. The ℋ
∞ controllers are designed with γ
=
4.0, which was chosen to satisfy the linear version of the condition in (17); that is, such that the matrix in (48) is negative definite.The methods to be used, in addition to theℋ
∞ controllers of the previous section, are the linear and nonlinear proportional-derivative (PD) laws of Tsiotras [1], the open-loop (OL) optimal control method by Schaub et al. [2], the closed-loop (CL) optimal nonlinear method of Tewari [3], and the sum of squares (SOS) approach of Gollu and Rodrigues [4]. In order to gauge the performance and properly compare the different controllers, we will use the following metrics:
(49)
E
rms
=
[
1
T
∫
0
T
‍
y
2
T
(
t
)
y
2
(
t
)
d
t
]
(
1
/
2
)
,
T
rms
=
[
1
T
∫
0
T
‍
u
2
T
(
t
)
u
2
(
t
)
d
t
]
(
1
/
2
)
,
where y
2
=
y
0
-
y
1
=
-
y
1
=
-
x
=
-
[
ω
1
ω
2
ω
3
σ
1
σ
2
σ
3
]
T, T is the orbital period, E
rms is the RMS tracking error, and T
rms is the RMS control torque (note that -
u
2 are the control torques).


### 4.1. Existing Control Methods

We now provide a brief overview of the controller synthesis methods to be used for comparisons with our method. The interested reader is referred to the appropriate literature for a more detailed exposition of these methods.


#### 4.1.1. Proportional-Derivative Controllers

Tsiotras [1] developed two proportional-derivative (PD) control laws for the attitude control problem. The linear PD controller is given by
(50)
u
2
(
t
)
=
ω
+
2
σ
,
while the nonlinear version is given by
(51)
u
2
(
t
)
=
ω
+
2
(
1
+
σ
T
σ
)
σ
.


#### 4.1.2. Open-Loop Optimal Control

The optimal control method presented by Schaub et al. [2] is designed to minimize the cost function
(52)
J
=
1
2
K
1
g
(
σ
(
t
f
)
)
+
1
2
ω
T
(
t
f
)
K
2
ω
(
t
f
)
︸
ϕ
(
t
f
)
+
∫
t
0
t
f
‍
1
2
[
K
3
g
(
σ
)
+
ω
T
K
4
ω
+
u
2
T
R
u
2
]
︸
p
(
x
,
u
2
,
t
)
d
t
,
where K
1 and K
3 are scalar weights, K
2, K
4, and R are weighting matrices, and
(53)
g
(
σ
)
=
4
σ
T
σ
(
1
+
σ
T
σ
)
2
.
The Hamiltonian relating to this optimal control problem is defined as
(54)
H
=
p
(
x
,
u
,
t
)
+
Λ
T
[
a
(
x
)
-
B
u
2
(
t
)
]
,
where we have used the plant model in (1) and the second relation of (3) with y
1
=
0. The costates, denoted by Λ, have dynamics
(55)
Λ
˙
=
-
∂
H
∂
x
,
Λ
(
t
f
)
=
∂
ϕ
∂
x
|
t
=
t
f
.
Note that the costates are specified at some final time t
f, not the initial time. The optimal control law for this problem is determined from
(56)
∂
H
∂
u
2
=
0
,
which yields
(57)
u
2
(
t
)
=
R
-
1
B
T
Λ
.
The primary disadvantage of this method from a practical point of view is that it requires the solution of a two-point boundary value problem (TPBVP) and results in an open-loop control strategy. For the simulation results presented below, the following weighting parameters are used: K
1
=
5.0, K
2
=
5.0
1,  K
3
=
1.0, K
4
=
1, and R
=
1. Moreover, the maneuver is optimized for a final time t
f
=
60 s.


#### 4.1.3. Closed-Loop Optimal Controller

The optimal control method presented by Tewari [3] is based on obtaining an exact analytical solution to the Hamilton-Jacobi equation. Consider the HJE in (12) with parameters
(58)
A
(
x
)
=
a
(
x
)
,
R
(
x
)
=
B
R
-
1
B
T
,
Q
(
x
)
=
[
x
x
[
2
]
]
T
[
Q
11
Q
12
Q
21
Q
22
]
︸
Q
[
x
x
[
2
]
]
,
where the inertia matrix is defined as I
=
dia
g
i
{
I
i
}, R
=
dia
g
i
{
R
i
} is symmetric and positive definite, Q is symmetric and positive semidefinite, and x
[
2
]
=
[
ω
1
2
ω
2
2
ω
3
2
σ
1
2
σ
2
2
σ
3
2
]
T. It is assumed that V
(
x
) has the same form as Q
(
x
); thus,
(59)
V
(
x
)
=
1
2
[
x
x
[
2
]
]
T
[
P
11
P
12
P
21
P
22
]
︸
P
[
x
x
[
2
]
]
,
where P is symmetric and positive definite. The state feedback controller is given by
(60)
u
2
(
t
)
=
R
-
1
B
T
∇
V
T
(
x
)
,
where
(61)
∇
V
T
(
x
)
=
P
11
x
+
2
dia
g
i
{
x
i
}
P
12
x
+
P
12
x
[
2
]
+
2
dia
g
i
{
x
i
}
P
22
x
[
2
]
.The matrixP is obtained as follows. First, P
11 is calculated from the algebraic Riccati equation in (27) corresponding to the parameters in (58), that is, with R in (27) replaced by B
R
-
1
B
T and Q in (27) replaced by Q
11. Then, the equations in (26) with k
=
3
,
…
,
6, are solved simultaneously for the remaining unknowns. In particular, the nonzero elements of the matrices P
12 and P
22 are
(62)
P
12
11
=
I
1
2
R
1
(
2
Q
12
11
+
k
2
P
11
12
+
k
3
P
11
13
)
(
6
P
11
11
)
,
P
12
16
=
I
1
2
R
1
Q
12
16
(
2
P
11
11
)
,
P
12
22
=
I
2
2
R
2
(
Q
12
12
-
2
k
2
P
11
22
)
(
3
P
11
22
)
,
P
12
25
=
I
2
2
R
2
(
2
Q
12
25
-
3
P
12
22
P
11
25
)
(
2
P
11
22
)
,
P
12
33
=
I
3
2
R
3
Q
12
33
(
3
P
11
33
)
,
P
22
11
=
[
I
1
2
R
1
Q
22
11
-
9
(
P
12
11
)
2
]
(
2
P
11
11
)
,
P
22
12
=
Q
22
12
[
P
11
11
/
(
I
1
2
R
1
)
+
P
11
22
/
(
I
2
2
R
2
)
]
,
P
22
13
=
Q
22
13
[
P
11
11
/
(
I
1
2
R
1
)
+
P
11
33
/
(
I
3
2
R
3
)
]
,
P
22
22
=
[
I
2
2
R
2
Q
22
22
-
9
(
P
12
22
)
2
]
(
2
P
11
22
)
,
P
22
23
=
Q
22
23
[
P
11
33
/
(
I
3
2
R
3
)
+
P
11
22
/
(
I
2
2
R
2
)
]
,
P
22
33
=
[
I
3
2
R
3
Q
22
33
-
9
(
P
12
33
)
2
]
(
2
P
11
33
)
,
where, as before, k
1
=
(
I
2
-
I
3
)
/
I
1, k
2
=
(
I
3
-
I
1
)
/
I
2 and k
3
=
(
I
1
-
I
2
)
/
I
3.For the simulation results presented below, the following weighting parameters are used [3]:
(63)
R
=
dia
g
i
=
{
I
i
-
2
}
,
Q
11
=
(
0.01
)
1
,
Q
12
=
[
(
0.1
)
1
0
0
0
]
,
Q
22
=
[
P
11
11
P
11
12
P
11
13
P
11
21
P
11
22
P
11
23
0
P
11
31
P
11
32
P
11
33
0
0
]
.We now make a few remarks regarding the characteristics of the synthesis method of Tewari [3]. Like the controller presented in this paper, Tewari's closed-loop optimal method results in a polynomial feedback controller. However, his solution is derived specifically for the attitude control problem. In contrast, the method developed in the present paper is applicable to a wider class of systems.


#### 4.1.4. Sum of Squares Controller

A multivariate polynomialp
(
x
) is a sum of squares (SOS) if there exist some polynomials f
i
(
x
), i
=
1
,
…
,
M, such that
(64)
p
(
x
)
=
∑
i
=
1
M
‍
f
i
2
(
x
)
.
The SOS controller synthesis approach relaxes the search for positive- or negative-definite functions to a search for SOS functions. It should be noted, however, that the use of sums of squares is conservative, since p
(
x
) being SOS implies that p
(
x
)
≥
0, while the converse is not true in general. Also, it can be shown that if v
T
P
(
x
)
v is SOS for v
∈
ℝ
n, then P
(
x
)
≥
0 for all x
∈
ℝ
N [36].In applying the SOS controller synthesis method, we first rewrite the system in (1) as
(65)
x
˙
=
A
(
x
)
x
-
B
u
2
(
t
)
,
where we have used the second relation of (3) with y
1
=
0. Note that the matrix A
(
x
) is not unique. The SOS state feedback controller for this problem is given by [4]
(66)
u
2
(
t
)
=
-
K
(
x
)
P
-
1
x
,
where we note that P is now taken to be constant (state independent).Consider the Lyapunov function(67)
V
(
x
)
=
x
T
P
-
1
x
.
Taking the time derivative of this function along the trajectories of the system in (65) with the controller of (66) yields
(68)
V
˙
(
x
)
=
[
A
(
x
)
x
+
B
K
(
x
)
P
-
1
x
]
T
P
-
1
x
+
x
T
P
-
1
[
A
(
x
)
x
+
B
K
(
x
)
P
-
1
x
]
.
Using the change of variables x
=
P
v, this last expression can be written as
(69)
V
˙
(
x
)
=
v
T
[
P
A
T
(
x
)
+
A
(
x
)
P
+
K
T
(
x
)
B
T
+
B
K
(
x
)
]
v
.
The conditions V
(
x
)
>
0 and V
˙
(
x
)
<
0 can be replaced by the conditions that  P is positive definite and -
V
˙
(
x
) is SOS. The second of these conditions will be strengthened to -
[
V
˙
(
x
)
+
ϵ
(
x
)
] being SOS, where ϵ
(
x
) is some SOS function. This semidefinite programming (SDP) problem can then be written as follows:
(70)
find
P
,
K
(
x
)
,
ϵ
(
x
)
s
.
t
.
P
=
P
T
>
0
ϵ
(
x
)
is
SOS
-
v
T
[
P
A
T
(
x
)
+
A
(
x
)
P
+
K
T
(
x
)
B
T
P
A
T
+
B
K
(
x
)
+
ϵ
(
x
)
1
]
v
is
SOS
.
This optimization problem can be solved using the SOSTOOLS software package [37].


### 4.2. Disturbance Rejection

We begin by comparing the controllers presented in this paper with the methods from the literature with regard to disturbance rejection. The two disturbances considered here are the gravity-gradient and geomagnetic torques. For the purposes of this comparison, the simulations start from the desired attitude (i.e.,y
0
=
0) and we compare the ability of the different controllers to maintain that attitude over one complete orbit. Table 1 presents values of the performance metrics for the different controllers.Table 1
Controller disturbance rejection.
ℋ
∞
ℋ
∞
PD
Optimal
SOS
(1st-order)
(2nd-, 3rd-, and 4th-order)
(both)
(CL)
E
rms (×
1
0
-
6)
10.809686
10.809683
1.989095
5.328614
1.319846
T
rms (N·m) (×
1
0
-
6)
3.979559
3.979559
3.977874
3.977763
3.976904The optimal control method of Schaub et al. [2] is not included in the table because it is unable to reject any disturbances, which is entirely due to its open-loop nature. There is no apparent difference in performance between the linear and nonlinear PD controllers. There is also very little difference between any of the nonlinear ℋ
∞ controllers developed using the present method. This is due to the very small magnitude of these disturbance torques, which the linear term in the controller can effectively overcome. With regard to disturbance rejection, the present ℋ
∞ controllers do not perform as well as the existing methods.


### 4.3. Response to Initial Conditions

In this and the following subsections, the disturbance torques are set to zero and the nonzero initial conditions noted at the beginning of this section are applied. All of the resulting controllers described previously are stable.The resulting attitude and control torques when the fourth-orderℋ
∞ controller is applied are given in Figures 3 and 4. The MRP switching can clearly be seen in Figure 3, where we also note that it requires several rotations for the controller to sufficiently slow down the spacecraft. This is due to the small torques applied, as seen in Figure 4. It is important to note that the control torques are continuous, which is a result of the dynamic aspect of the controller developed in this paper. Once again, there is little to distinguish the behaviour of the ℋ
∞ controllers of varying order of approximation. This will be discussed further after the robustness of these controllers has been assessed.Figure 3
MRP trajectories forI
=
diag
⁡
{
10,6.3,8.5
} kg·m2 with 4th-order ℋ
∞ controller.Figure 4
Control torques forI
=
diag
⁡
{
10,6.3,8.5
} kg·m2 with 4th-order ℋ
∞ controller.


### 4.4. Robustness to Actuation Time Delay

The robustness properties of the different controllers are now examined with regard to a time delay in the actuation. Such a delay could represent the finite time required by a satellite on-board computer to take the sensor measurements and calculate the required control signal. The time delay is made equal to an integer number of the numerical integration step-sizeΔ
t. Table 2 indicates the maximum allowable time delay, h
max
⁡, for each controller such that the desired attitude maneuver is achieved within one orbit. As can be seen from these results, the four ℋ
∞ controllers are more robust with regard to this effect than the other control methods. In particular, the present fourth-order ℋ
∞ controller is the most robust. We also note that the open-loop method of Schaub et al. [2] has no robustness to this effect; the time delay results in a nonzero final angular velocity such that the system will rotate endlessly.Table 2
Robustness to actuation time delay.
ℋ
∞
ℋ
∞
ℋ
∞
ℋ
∞
(1st-order)
(2nd-order)
(3rd-order)
(4th-order)
h
max
⁡ (s)
3.93
3.94
3.97
3.98
PD
PD
Optimal
Optimal
SOS
(linear)
(nonlinear)
(OL)
(CL)
h
max
⁡ (s)
1.69
1.09
0.00
1.41
0.81


### 4.5. Robustness to Unmodeled Actuator Dynamics

The robustness properties of the different controllers are now examined with regard to unmodeled actuator dynamics. For the purposes of this study, we make use of first- and second-order actuator models. In practice the actuator dynamics may be far more complex than the ones used here. However, we use these simple models here for the purposes of studying the capabilities of the different control methods. Evidently, as the actuator bandwidth decreases, it becomes harder for the controller to stabilize the system. Thus, we are able to infer the relative robustness properties of the different controllers by examining Figures5–8. In particular, the farther a line reaches towards the left-hand side of the graph, the more robust that controller is with regard to the unmodeled actuator dynamics. As we shall see, the present ℋ
∞ controllers are always more robust than the other methods.Figure 5
RMS tracking error with respect to 1st-order actuator bandwidth.Figure 6
RMS control effort with respect to 1st-order actuator bandwidth.Figure 7
RMS tracking error with respect to 2nd-order actuator bandwidth.Figure 8
RMS control effort with respect to 2nd-order actuator bandwidth.


#### 4.5.1. Unmodeled 1st-Order Actuator Dynamics

We begin by examining the robustness of the controllers with regard to unmodeled first-order actuator dynamics. This shall be accomplished by including the following first-order dynamics model in each component (i
=
1,2
,
3) of the controller output:
(71)
x
˙
=
-
ω
b
x
+
ω
b
u
2
,
i
,
y
i
=
x
,
where ω
b is the actuator bandwidth. Figures 5 and 6 show E
rms and T
rms, respectively, as a function of the actuator bandwidth. It is noted from these figures that all four ℋ
∞ controllers provide nearly the same tracking error and control effort. Moreover, while the closed-loop optimal control law provides better tracking error compared with the ℋ
∞ controllers, the trade-off is that it requires more control effort. Similarly, the SOS controller yields very low tracking error at the expense of greater control effort. The two PD laws, on the other hand, have a higher tracking error and control effort than the other methods.


#### 4.5.2. Unmodeled 2nd-Order Actuator Dynamics

The robustness properties of the different controllers are now examined with regard to unmodeled second-order actuator dynamics. This shall be accomplished by including the following second-order dynamics model in each component (i
=
1,2
,
3) of the controller output:
(72)
[
x
˙
1
x
˙
2
]
=
[
0
1
-
ω
b
2
-
2
ζ
ω
b
]
[
x
1
x
2
]
+
[
0
ω
b
2
]
u
2
,
i
,
y
i
=
x
1
,
where ζ is the actuator damping ratio and ω
b is the actuator bandwidth. All simulations were performed with a damping ratio ζ
=
0.5. Figures 7 and 8 show E
rms and T
rms, respectively, as a function of the actuator bandwidth. It is seen from these figures that the various control methods follow the same trends as in the case of the first-order actuator dynamics.


### 4.6. Robustness in the Gap Metric (Revisited)

We now make use of the gap metric [28] to characterize the difference in the input-output (IO) map of the system induced by the unmodeled actuator dynamics. However, since we cannot calculate the gap between two nonlinear systems, we calculate the gap metric distance for the linearized system only. As the actuator bandwidth and damping ratio change, the value of the gap metric will also vary. Figure 9 shows the gap metric value with respect to the first-order actuator bandwidth. Figure 10 shows the gap metric value with respect to the second-order actuator bandwidth for a damping ratio ζ
=
0.5 as used in the above numerical simulations. As the damping ratio varies, the curve of this figure moves to the left or right slightly, although there does not appear to be any discernible trend. Moreover, it should be noted that as the actuator bandwidth approaches infinity, the effects of the actuator in the controller input-output map become negligible. This is to be expected and is seen in both Figures 9 and 10, where the gap metric approaches zero with increasing actuator bandwidth.Figure 9
Gap metric with respect to 1st-order actuator bandwidth.Figure 10
Gap metric with respect to 2nd-order actuator bandwidth.We now return for a moment to the question of controller robustness in the gap metric. From (9) we have the following small-gain-type stability criterion:
(73)
ρ
g
(
P
0
,
P
1
)
·
∥
Π
𝒩
∥
ℳ
0
∥
∞
︸
≤
γ
<
1
.
In the numerical simulations performed with γ
=
4, it was determined that the closed-loop system is stable in the presence of the unmodeled second-order actuator with bandwidth as low as ω
b
=
0.3 rad/s. With this actuator bandwidth, the calculation of the gap between the linear attitude dynamics with and without the actuator results in a value of ρ
g
=
0.522, which does not satisfy the condition of (73). This could be explained by the conservativeness of the small-gain criterion and the fact that it is a sufficient (but not necessary) condition for stability. However, it can also be attributed to the nonlinear effects of the dynamics not taken into account in the gap calculation here. This emphasizes the need for a method to calculate the gap metric for nonlinear systems.


### 4.7. Overall Assessment of theℋ
∞ Controllers

The previous sections show that the methodology presented here can be used to develop very robust attitude controllers. However, their performance was not as good as some of the other techniques to which it was compared. We have not demonstrated significant benefits for the higher-order terms in the controller which, while disappointing, is an important contribution to the literature. The fact that a linear controller can perform very closely to the higher-order ones is a strong vindication of the most popular approach used for actual attitude control: linear feedback of angular velocity and attitude information. This is consistent with [38] which showed that a linear combination of angular velocity and quaternion (Euler parameter) feedback solves the state feedback nonlinear (suboptimal)-ℋ
∞ control problem for rigid spacecraft attitude control. Our results are entirely consistent with this and we strongly suspect that results analogous to [38] can be obtained for the case of angular velocity and MRPs.


## 4.1. Existing Control Methods

We now provide a brief overview of the controller synthesis methods to be used for comparisons with our method. The interested reader is referred to the appropriate literature for a more detailed exposition of these methods.


### 4.1.1. Proportional-Derivative Controllers

Tsiotras [1] developed two proportional-derivative (PD) control laws for the attitude control problem. The linear PD controller is given by
(50)
u
2
(
t
)
=
ω
+
2
σ
,
while the nonlinear version is given by
(51)
u
2
(
t
)
=
ω
+
2
(
1
+
σ
T
σ
)
σ
.


### 4.1.2. Open-Loop Optimal Control

The optimal control method presented by Schaub et al. [2] is designed to minimize the cost function
(52)
J
=
1
2
K
1
g
(
σ
(
t
f
)
)
+
1
2
ω
T
(
t
f
)
K
2
ω
(
t
f
)
︸
ϕ
(
t
f
)
+
∫
t
0
t
f
‍
1
2
[
K
3
g
(
σ
)
+
ω
T
K
4
ω
+
u
2
T
R
u
2
]
︸
p
(
x
,
u
2
,
t
)
d
t
,
where K
1 and K
3 are scalar weights, K
2, K
4, and R are weighting matrices, and
(53)
g
(
σ
)
=
4
σ
T
σ
(
1
+
σ
T
σ
)
2
.
The Hamiltonian relating to this optimal control problem is defined as
(54)
H
=
p
(
x
,
u
,
t
)
+
Λ
T
[
a
(
x
)
-
B
u
2
(
t
)
]
,
where we have used the plant model in (1) and the second relation of (3) with y
1
=
0. The costates, denoted by Λ, have dynamics
(55)
Λ
˙
=
-
∂
H
∂
x
,
Λ
(
t
f
)
=
∂
ϕ
∂
x
|
t
=
t
f
.
Note that the costates are specified at some final time t
f, not the initial time. The optimal control law for this problem is determined from
(56)
∂
H
∂
u
2
=
0
,
which yields
(57)
u
2
(
t
)
=
R
-
1
B
T
Λ
.
The primary disadvantage of this method from a practical point of view is that it requires the solution of a two-point boundary value problem (TPBVP) and results in an open-loop control strategy. For the simulation results presented below, the following weighting parameters are used: K
1
=
5.0, K
2
=
5.0
1,  K
3
=
1.0, K
4
=
1, and R
=
1. Moreover, the maneuver is optimized for a final time t
f
=
60 s.


### 4.1.3. Closed-Loop Optimal Controller

The optimal control method presented by Tewari [3] is based on obtaining an exact analytical solution to the Hamilton-Jacobi equation. Consider the HJE in (12) with parameters
(58)
A
(
x
)
=
a
(
x
)
,
R
(
x
)
=
B
R
-
1
B
T
,
Q
(
x
)
=
[
x
x
[
2
]
]
T
[
Q
11
Q
12
Q
21
Q
22
]
︸
Q
[
x
x
[
2
]
]
,
where the inertia matrix is defined as I
=
dia
g
i
{
I
i
}, R
=
dia
g
i
{
R
i
} is symmetric and positive definite, Q is symmetric and positive semidefinite, and x
[
2
]
=
[
ω
1
2
ω
2
2
ω
3
2
σ
1
2
σ
2
2
σ
3
2
]
T. It is assumed that V
(
x
) has the same form as Q
(
x
); thus,
(59)
V
(
x
)
=
1
2
[
x
x
[
2
]
]
T
[
P
11
P
12
P
21
P
22
]
︸
P
[
x
x
[
2
]
]
,
where P is symmetric and positive definite. The state feedback controller is given by
(60)
u
2
(
t
)
=
R
-
1
B
T
∇
V
T
(
x
)
,
where
(61)
∇
V
T
(
x
)
=
P
11
x
+
2
dia
g
i
{
x
i
}
P
12
x
+
P
12
x
[
2
]
+
2
dia
g
i
{
x
i
}
P
22
x
[
2
]
.The matrixP is obtained as follows. First, P
11 is calculated from the algebraic Riccati equation in (27) corresponding to the parameters in (58), that is, with R in (27) replaced by B
R
-
1
B
T and Q in (27) replaced by Q
11. Then, the equations in (26) with k
=
3
,
…
,
6, are solved simultaneously for the remaining unknowns. In particular, the nonzero elements of the matrices P
12 and P
22 are
(62)
P
12
11
=
I
1
2
R
1
(
2
Q
12
11
+
k
2
P
11
12
+
k
3
P
11
13
)
(
6
P
11
11
)
,
P
12
16
=
I
1
2
R
1
Q
12
16
(
2
P
11
11
)
,
P
12
22
=
I
2
2
R
2
(
Q
12
12
-
2
k
2
P
11
22
)
(
3
P
11
22
)
,
P
12
25
=
I
2
2
R
2
(
2
Q
12
25
-
3
P
12
22
P
11
25
)
(
2
P
11
22
)
,
P
12
33
=
I
3
2
R
3
Q
12
33
(
3
P
11
33
)
,
P
22
11
=
[
I
1
2
R
1
Q
22
11
-
9
(
P
12
11
)
2
]
(
2
P
11
11
)
,
P
22
12
=
Q
22
12
[
P
11
11
/
(
I
1
2
R
1
)
+
P
11
22
/
(
I
2
2
R
2
)
]
,
P
22
13
=
Q
22
13
[
P
11
11
/
(
I
1
2
R
1
)
+
P
11
33
/
(
I
3
2
R
3
)
]
,
P
22
22
=
[
I
2
2
R
2
Q
22
22
-
9
(
P
12
22
)
2
]
(
2
P
11
22
)
,
P
22
23
=
Q
22
23
[
P
11
33
/
(
I
3
2
R
3
)
+
P
11
22
/
(
I
2
2
R
2
)
]
,
P
22
33
=
[
I
3
2
R
3
Q
22
33
-
9
(
P
12
33
)
2
]
(
2
P
11
33
)
,
where, as before, k
1
=
(
I
2
-
I
3
)
/
I
1, k
2
=
(
I
3
-
I
1
)
/
I
2 and k
3
=
(
I
1
-
I
2
)
/
I
3.For the simulation results presented below, the following weighting parameters are used [3]:
(63)
R
=
dia
g
i
=
{
I
i
-
2
}
,
Q
11
=
(
0.01
)
1
,
Q
12
=
[
(
0.1
)
1
0
0
0
]
,
Q
22
=
[
P
11
11
P
11
12
P
11
13
P
11
21
P
11
22
P
11
23
0
P
11
31
P
11
32
P
11
33
0
0
]
.We now make a few remarks regarding the characteristics of the synthesis method of Tewari [3]. Like the controller presented in this paper, Tewari's closed-loop optimal method results in a polynomial feedback controller. However, his solution is derived specifically for the attitude control problem. In contrast, the method developed in the present paper is applicable to a wider class of systems.


### 4.1.4. Sum of Squares Controller

A multivariate polynomialp
(
x
) is a sum of squares (SOS) if there exist some polynomials f
i
(
x
), i
=
1
,
…
,
M, such that
(64)
p
(
x
)
=
∑
i
=
1
M
‍
f
i
2
(
x
)
.
The SOS controller synthesis approach relaxes the search for positive- or negative-definite functions to a search for SOS functions. It should be noted, however, that the use of sums of squares is conservative, since p
(
x
) being SOS implies that p
(
x
)
≥
0, while the converse is not true in general. Also, it can be shown that if v
T
P
(
x
)
v is SOS for v
∈
ℝ
n, then P
(
x
)
≥
0 for all x
∈
ℝ
N [36].In applying the SOS controller synthesis method, we first rewrite the system in (1) as
(65)
x
˙
=
A
(
x
)
x
-
B
u
2
(
t
)
,
where we have used the second relation of (3) with y
1
=
0. Note that the matrix A
(
x
) is not unique. The SOS state feedback controller for this problem is given by [4]
(66)
u
2
(
t
)
=
-
K
(
x
)
P
-
1
x
,
where we note that P is now taken to be constant (state independent).Consider the Lyapunov function(67)
V
(
x
)
=
x
T
P
-
1
x
.
Taking the time derivative of this function along the trajectories of the system in (65) with the controller of (66) yields
(68)
V
˙
(
x
)
=
[
A
(
x
)
x
+
B
K
(
x
)
P
-
1
x
]
T
P
-
1
x
+
x
T
P
-
1
[
A
(
x
)
x
+
B
K
(
x
)
P
-
1
x
]
.
Using the change of variables x
=
P
v, this last expression can be written as
(69)
V
˙
(
x
)
=
v
T
[
P
A
T
(
x
)
+
A
(
x
)
P
+
K
T
(
x
)
B
T
+
B
K
(
x
)
]
v
.
The conditions V
(
x
)
>
0 and V
˙
(
x
)
<
0 can be replaced by the conditions that  P is positive definite and -
V
˙
(
x
) is SOS. The second of these conditions will be strengthened to -
[
V
˙
(
x
)
+
ϵ
(
x
)
] being SOS, where ϵ
(
x
) is some SOS function. This semidefinite programming (SDP) problem can then be written as follows:
(70)
find
P
,
K
(
x
)
,
ϵ
(
x
)
s
.
t
.
P
=
P
T
>
0
ϵ
(
x
)
is
SOS
-
v
T
[
P
A
T
(
x
)
+
A
(
x
)
P
+
K
T
(
x
)
B
T
P
A
T
+
B
K
(
x
)
+
ϵ
(
x
)
1
]
v
is
SOS
.
This optimization problem can be solved using the SOSTOOLS software package [37].


## 4.1.1. Proportional-Derivative Controllers

Tsiotras [1] developed two proportional-derivative (PD) control laws for the attitude control problem. The linear PD controller is given by
(50)
u
2
(
t
)
=
ω
+
2
σ
,
while the nonlinear version is given by
(51)
u
2
(
t
)
=
ω
+
2
(
1
+
σ
T
σ
)
σ
.


## 4.1.2. Open-Loop Optimal Control

The optimal control method presented by Schaub et al. [2] is designed to minimize the cost function
(52)
J
=
1
2
K
1
g
(
σ
(
t
f
)
)
+
1
2
ω
T
(
t
f
)
K
2
ω
(
t
f
)
︸
ϕ
(
t
f
)
+
∫
t
0
t
f
‍
1
2
[
K
3
g
(
σ
)
+
ω
T
K
4
ω
+
u
2
T
R
u
2
]
︸
p
(
x
,
u
2
,
t
)
d
t
,
where K
1 and K
3 are scalar weights, K
2, K
4, and R are weighting matrices, and
(53)
g
(
σ
)
=
4
σ
T
σ
(
1
+
σ
T
σ
)
2
.
The Hamiltonian relating to this optimal control problem is defined as
(54)
H
=
p
(
x
,
u
,
t
)
+
Λ
T
[
a
(
x
)
-
B
u
2
(
t
)
]
,
where we have used the plant model in (1) and the second relation of (3) with y
1
=
0. The costates, denoted by Λ, have dynamics
(55)
Λ
˙
=
-
∂
H
∂
x
,
Λ
(
t
f
)
=
∂
ϕ
∂
x
|
t
=
t
f
.
Note that the costates are specified at some final time t
f, not the initial time. The optimal control law for this problem is determined from
(56)
∂
H
∂
u
2
=
0
,
which yields
(57)
u
2
(
t
)
=
R
-
1
B
T
Λ
.
The primary disadvantage of this method from a practical point of view is that it requires the solution of a two-point boundary value problem (TPBVP) and results in an open-loop control strategy. For the simulation results presented below, the following weighting parameters are used: K
1
=
5.0, K
2
=
5.0
1,  K
3
=
1.0, K
4
=
1, and R
=
1. Moreover, the maneuver is optimized for a final time t
f
=
60 s.


## 4.1.3. Closed-Loop Optimal Controller

The optimal control method presented by Tewari [3] is based on obtaining an exact analytical solution to the Hamilton-Jacobi equation. Consider the HJE in (12) with parameters
(58)
A
(
x
)
=
a
(
x
)
,
R
(
x
)
=
B
R
-
1
B
T
,
Q
(
x
)
=
[
x
x
[
2
]
]
T
[
Q
11
Q
12
Q
21
Q
22
]
︸
Q
[
x
x
[
2
]
]
,
where the inertia matrix is defined as I
=
dia
g
i
{
I
i
}, R
=
dia
g
i
{
R
i
} is symmetric and positive definite, Q is symmetric and positive semidefinite, and x
[
2
]
=
[
ω
1
2
ω
2
2
ω
3
2
σ
1
2
σ
2
2
σ
3
2
]
T. It is assumed that V
(
x
) has the same form as Q
(
x
); thus,
(59)
V
(
x
)
=
1
2
[
x
x
[
2
]
]
T
[
P
11
P
12
P
21
P
22
]
︸
P
[
x
x
[
2
]
]
,
where P is symmetric and positive definite. The state feedback controller is given by
(60)
u
2
(
t
)
=
R
-
1
B
T
∇
V
T
(
x
)
,
where
(61)
∇
V
T
(
x
)
=
P
11
x
+
2
dia
g
i
{
x
i
}
P
12
x
+
P
12
x
[
2
]
+
2
dia
g
i
{
x
i
}
P
22
x
[
2
]
.The matrixP is obtained as follows. First, P
11 is calculated from the algebraic Riccati equation in (27) corresponding to the parameters in (58), that is, with R in (27) replaced by B
R
-
1
B
T and Q in (27) replaced by Q
11. Then, the equations in (26) with k
=
3
,
…
,
6, are solved simultaneously for the remaining unknowns. In particular, the nonzero elements of the matrices P
12 and P
22 are
(62)
P
12
11
=
I
1
2
R
1
(
2
Q
12
11
+
k
2
P
11
12
+
k
3
P
11
13
)
(
6
P
11
11
)
,
P
12
16
=
I
1
2
R
1
Q
12
16
(
2
P
11
11
)
,
P
12
22
=
I
2
2
R
2
(
Q
12
12
-
2
k
2
P
11
22
)
(
3
P
11
22
)
,
P
12
25
=
I
2
2
R
2
(
2
Q
12
25
-
3
P
12
22
P
11
25
)
(
2
P
11
22
)
,
P
12
33
=
I
3
2
R
3
Q
12
33
(
3
P
11
33
)
,
P
22
11
=
[
I
1
2
R
1
Q
22
11
-
9
(
P
12
11
)
2
]
(
2
P
11
11
)
,
P
22
12
=
Q
22
12
[
P
11
11
/
(
I
1
2
R
1
)
+
P
11
22
/
(
I
2
2
R
2
)
]
,
P
22
13
=
Q
22
13
[
P
11
11
/
(
I
1
2
R
1
)
+
P
11
33
/
(
I
3
2
R
3
)
]
,
P
22
22
=
[
I
2
2
R
2
Q
22
22
-
9
(
P
12
22
)
2
]
(
2
P
11
22
)
,
P
22
23
=
Q
22
23
[
P
11
33
/
(
I
3
2
R
3
)
+
P
11
22
/
(
I
2
2
R
2
)
]
,
P
22
33
=
[
I
3
2
R
3
Q
22
33
-
9
(
P
12
33
)
2
]
(
2
P
11
33
)
,
where, as before, k
1
=
(
I
2
-
I
3
)
/
I
1, k
2
=
(
I
3
-
I
1
)
/
I
2 and k
3
=
(
I
1
-
I
2
)
/
I
3.For the simulation results presented below, the following weighting parameters are used [3]:
(63)
R
=
dia
g
i
=
{
I
i
-
2
}
,
Q
11
=
(
0.01
)
1
,
Q
12
=
[
(
0.1
)
1
0
0
0
]
,
Q
22
=
[
P
11
11
P
11
12
P
11
13
P
11
21
P
11
22
P
11
23
0
P
11
31
P
11
32
P
11
33
0
0
]
.We now make a few remarks regarding the characteristics of the synthesis method of Tewari [3]. Like the controller presented in this paper, Tewari's closed-loop optimal method results in a polynomial feedback controller. However, his solution is derived specifically for the attitude control problem. In contrast, the method developed in the present paper is applicable to a wider class of systems.


## 4.1.4. Sum of Squares Controller

A multivariate polynomialp
(
x
) is a sum of squares (SOS) if there exist some polynomials f
i
(
x
), i
=
1
,
…
,
M, such that
(64)
p
(
x
)
=
∑
i
=
1
M
‍
f
i
2
(
x
)
.
The SOS controller synthesis approach relaxes the search for positive- or negative-definite functions to a search for SOS functions. It should be noted, however, that the use of sums of squares is conservative, since p
(
x
) being SOS implies that p
(
x
)
≥
0, while the converse is not true in general. Also, it can be shown that if v
T
P
(
x
)
v is SOS for v
∈
ℝ
n, then P
(
x
)
≥
0 for all x
∈
ℝ
N [36].In applying the SOS controller synthesis method, we first rewrite the system in (1) as
(65)
x
˙
=
A
(
x
)
x
-
B
u
2
(
t
)
,
where we have used the second relation of (3) with y
1
=
0. Note that the matrix A
(
x
) is not unique. The SOS state feedback controller for this problem is given by [4]
(66)
u
2
(
t
)
=
-
K
(
x
)
P
-
1
x
,
where we note that P is now taken to be constant (state independent).Consider the Lyapunov function(67)
V
(
x
)
=
x
T
P
-
1
x
.
Taking the time derivative of this function along the trajectories of the system in (65) with the controller of (66) yields
(68)
V
˙
(
x
)
=
[
A
(
x
)
x
+
B
K
(
x
)
P
-
1
x
]
T
P
-
1
x
+
x
T
P
-
1
[
A
(
x
)
x
+
B
K
(
x
)
P
-
1
x
]
.
Using the change of variables x
=
P
v, this last expression can be written as
(69)
V
˙
(
x
)
=
v
T
[
P
A
T
(
x
)
+
A
(
x
)
P
+
K
T
(
x
)
B
T
+
B
K
(
x
)
]
v
.
The conditions V
(
x
)
>
0 and V
˙
(
x
)
<
0 can be replaced by the conditions that  P is positive definite and -
V
˙
(
x
) is SOS. The second of these conditions will be strengthened to -
[
V
˙
(
x
)
+
ϵ
(
x
)
] being SOS, where ϵ
(
x
) is some SOS function. This semidefinite programming (SDP) problem can then be written as follows:
(70)
find
P
,
K
(
x
)
,
ϵ
(
x
)
s
.
t
.
P
=
P
T
>
0
ϵ
(
x
)
is
SOS
-
v
T
[
P
A
T
(
x
)
+
A
(
x
)
P
+
K
T
(
x
)
B
T
P
A
T
+
B
K
(
x
)
+
ϵ
(
x
)
1
]
v
is
SOS
.
This optimization problem can be solved using the SOSTOOLS software package [37].


## 4.2. Disturbance Rejection

We begin by comparing the controllers presented in this paper with the methods from the literature with regard to disturbance rejection. The two disturbances considered here are the gravity-gradient and geomagnetic torques. For the purposes of this comparison, the simulations start from the desired attitude (i.e.,y
0
=
0) and we compare the ability of the different controllers to maintain that attitude over one complete orbit. Table 1 presents values of the performance metrics for the different controllers.Table 1
Controller disturbance rejection.
ℋ
∞
ℋ
∞
PD
Optimal
SOS
(1st-order)
(2nd-, 3rd-, and 4th-order)
(both)
(CL)
E
rms (×
1
0
-
6)
10.809686
10.809683
1.989095
5.328614
1.319846
T
rms (N·m) (×
1
0
-
6)
3.979559
3.979559
3.977874
3.977763
3.976904The optimal control method of Schaub et al. [2] is not included in the table because it is unable to reject any disturbances, which is entirely due to its open-loop nature. There is no apparent difference in performance between the linear and nonlinear PD controllers. There is also very little difference between any of the nonlinear ℋ
∞ controllers developed using the present method. This is due to the very small magnitude of these disturbance torques, which the linear term in the controller can effectively overcome. With regard to disturbance rejection, the present ℋ
∞ controllers do not perform as well as the existing methods.


## 4.3. Response to Initial Conditions

In this and the following subsections, the disturbance torques are set to zero and the nonzero initial conditions noted at the beginning of this section are applied. All of the resulting controllers described previously are stable.The resulting attitude and control torques when the fourth-orderℋ
∞ controller is applied are given in Figures 3 and 4. The MRP switching can clearly be seen in Figure 3, where we also note that it requires several rotations for the controller to sufficiently slow down the spacecraft. This is due to the small torques applied, as seen in Figure 4. It is important to note that the control torques are continuous, which is a result of the dynamic aspect of the controller developed in this paper. Once again, there is little to distinguish the behaviour of the ℋ
∞ controllers of varying order of approximation. This will be discussed further after the robustness of these controllers has been assessed.Figure 3
MRP trajectories forI
=
diag
⁡
{
10,6.3,8.5
} kg·m2 with 4th-order ℋ
∞ controller.Figure 4
Control torques forI
=
diag
⁡
{
10,6.3,8.5
} kg·m2 with 4th-order ℋ
∞ controller.


## 4.4. Robustness to Actuation Time Delay

The robustness properties of the different controllers are now examined with regard to a time delay in the actuation. Such a delay could represent the finite time required by a satellite on-board computer to take the sensor measurements and calculate the required control signal. The time delay is made equal to an integer number of the numerical integration step-sizeΔ
t. Table 2 indicates the maximum allowable time delay, h
max
⁡, for each controller such that the desired attitude maneuver is achieved within one orbit. As can be seen from these results, the four ℋ
∞ controllers are more robust with regard to this effect than the other control methods. In particular, the present fourth-order ℋ
∞ controller is the most robust. We also note that the open-loop method of Schaub et al. [2] has no robustness to this effect; the time delay results in a nonzero final angular velocity such that the system will rotate endlessly.Table 2
Robustness to actuation time delay.
ℋ
∞
ℋ
∞
ℋ
∞
ℋ
∞
(1st-order)
(2nd-order)
(3rd-order)
(4th-order)
h
max
⁡ (s)
3.93
3.94
3.97
3.98
PD
PD
Optimal
Optimal
SOS
(linear)
(nonlinear)
(OL)
(CL)
h
max
⁡ (s)
1.69
1.09
0.00
1.41
0.81


## 4.5. Robustness to Unmodeled Actuator Dynamics

The robustness properties of the different controllers are now examined with regard to unmodeled actuator dynamics. For the purposes of this study, we make use of first- and second-order actuator models. In practice the actuator dynamics may be far more complex than the ones used here. However, we use these simple models here for the purposes of studying the capabilities of the different control methods. Evidently, as the actuator bandwidth decreases, it becomes harder for the controller to stabilize the system. Thus, we are able to infer the relative robustness properties of the different controllers by examining Figures5–8. In particular, the farther a line reaches towards the left-hand side of the graph, the more robust that controller is with regard to the unmodeled actuator dynamics. As we shall see, the present ℋ
∞ controllers are always more robust than the other methods.Figure 5
RMS tracking error with respect to 1st-order actuator bandwidth.Figure 6
RMS control effort with respect to 1st-order actuator bandwidth.Figure 7
RMS tracking error with respect to 2nd-order actuator bandwidth.Figure 8
RMS control effort with respect to 2nd-order actuator bandwidth.


### 4.5.1. Unmodeled 1st-Order Actuator Dynamics

We begin by examining the robustness of the controllers with regard to unmodeled first-order actuator dynamics. This shall be accomplished by including the following first-order dynamics model in each component (i
=
1,2
,
3) of the controller output:
(71)
x
˙
=
-
ω
b
x
+
ω
b
u
2
,
i
,
y
i
=
x
,
where ω
b is the actuator bandwidth. Figures 5 and 6 show E
rms and T
rms, respectively, as a function of the actuator bandwidth. It is noted from these figures that all four ℋ
∞ controllers provide nearly the same tracking error and control effort. Moreover, while the closed-loop optimal control law provides better tracking error compared with the ℋ
∞ controllers, the trade-off is that it requires more control effort. Similarly, the SOS controller yields very low tracking error at the expense of greater control effort. The two PD laws, on the other hand, have a higher tracking error and control effort than the other methods.


### 4.5.2. Unmodeled 2nd-Order Actuator Dynamics

The robustness properties of the different controllers are now examined with regard to unmodeled second-order actuator dynamics. This shall be accomplished by including the following second-order dynamics model in each component (i
=
1,2
,
3) of the controller output:
(72)
[
x
˙
1
x
˙
2
]
=
[
0
1
-
ω
b
2
-
2
ζ
ω
b
]
[
x
1
x
2
]
+
[
0
ω
b
2
]
u
2
,
i
,
y
i
=
x
1
,
where ζ is the actuator damping ratio and ω
b is the actuator bandwidth. All simulations were performed with a damping ratio ζ
=
0.5. Figures 7 and 8 show E
rms and T
rms, respectively, as a function of the actuator bandwidth. It is seen from these figures that the various control methods follow the same trends as in the case of the first-order actuator dynamics.


## 4.5.1. Unmodeled 1st-Order Actuator Dynamics

We begin by examining the robustness of the controllers with regard to unmodeled first-order actuator dynamics. This shall be accomplished by including the following first-order dynamics model in each component (i
=
1,2
,
3) of the controller output:
(71)
x
˙
=
-
ω
b
x
+
ω
b
u
2
,
i
,
y
i
=
x
,
where ω
b is the actuator bandwidth. Figures 5 and 6 show E
rms and T
rms, respectively, as a function of the actuator bandwidth. It is noted from these figures that all four ℋ
∞ controllers provide nearly the same tracking error and control effort. Moreover, while the closed-loop optimal control law provides better tracking error compared with the ℋ
∞ controllers, the trade-off is that it requires more control effort. Similarly, the SOS controller yields very low tracking error at the expense of greater control effort. The two PD laws, on the other hand, have a higher tracking error and control effort than the other methods.


## 4.5.2. Unmodeled 2nd-Order Actuator Dynamics

The robustness properties of the different controllers are now examined with regard to unmodeled second-order actuator dynamics. This shall be accomplished by including the following second-order dynamics model in each component (i
=
1,2
,
3) of the controller output:
(72)
[
x
˙
1
x
˙
2
]
=
[
0
1
-
ω
b
2
-
2
ζ
ω
b
]
[
x
1
x
2
]
+
[
0
ω
b
2
]
u
2
,
i
,
y
i
=
x
1
,
where ζ is the actuator damping ratio and ω
b is the actuator bandwidth. All simulations were performed with a damping ratio ζ
=
0.5. Figures 7 and 8 show E
rms and T
rms, respectively, as a function of the actuator bandwidth. It is seen from these figures that the various control methods follow the same trends as in the case of the first-order actuator dynamics.


## 4.6. Robustness in the Gap Metric (Revisited)

We now make use of the gap metric [28] to characterize the difference in the input-output (IO) map of the system induced by the unmodeled actuator dynamics. However, since we cannot calculate the gap between two nonlinear systems, we calculate the gap metric distance for the linearized system only. As the actuator bandwidth and damping ratio change, the value of the gap metric will also vary. Figure 9 shows the gap metric value with respect to the first-order actuator bandwidth. Figure 10 shows the gap metric value with respect to the second-order actuator bandwidth for a damping ratio ζ
=
0.5 as used in the above numerical simulations. As the damping ratio varies, the curve of this figure moves to the left or right slightly, although there does not appear to be any discernible trend. Moreover, it should be noted that as the actuator bandwidth approaches infinity, the effects of the actuator in the controller input-output map become negligible. This is to be expected and is seen in both Figures 9 and 10, where the gap metric approaches zero with increasing actuator bandwidth.Figure 9
Gap metric with respect to 1st-order actuator bandwidth.Figure 10
Gap metric with respect to 2nd-order actuator bandwidth.We now return for a moment to the question of controller robustness in the gap metric. From (9) we have the following small-gain-type stability criterion:
(73)
ρ
g
(
P
0
,
P
1
)
·
∥
Π
𝒩
∥
ℳ
0
∥
∞
︸
≤
γ
<
1
.
In the numerical simulations performed with γ
=
4, it was determined that the closed-loop system is stable in the presence of the unmodeled second-order actuator with bandwidth as low as ω
b
=
0.3 rad/s. With this actuator bandwidth, the calculation of the gap between the linear attitude dynamics with and without the actuator results in a value of ρ
g
=
0.522, which does not satisfy the condition of (73). This could be explained by the conservativeness of the small-gain criterion and the fact that it is a sufficient (but not necessary) condition for stability. However, it can also be attributed to the nonlinear effects of the dynamics not taken into account in the gap calculation here. This emphasizes the need for a method to calculate the gap metric for nonlinear systems.


## 4.7. Overall Assessment of theℋ
∞ Controllers

The previous sections show that the methodology presented here can be used to develop very robust attitude controllers. However, their performance was not as good as some of the other techniques to which it was compared. We have not demonstrated significant benefits for the higher-order terms in the controller which, while disappointing, is an important contribution to the literature. The fact that a linear controller can perform very closely to the higher-order ones is a strong vindication of the most popular approach used for actual attitude control: linear feedback of angular velocity and attitude information. This is consistent with [38] which showed that a linear combination of angular velocity and quaternion (Euler parameter) feedback solves the state feedback nonlinear (suboptimal)-ℋ
∞ control problem for rigid spacecraft attitude control. Our results are entirely consistent with this and we strongly suspect that results analogous to [38] can be obtained for the case of angular velocity and MRPs.


## 5. Conclusions

The results presented here clearly show the trade-off between performance and robustness. Existing methods from the attitude control literature typically focus on performance. In contrast, the method developed in this paper emphasizes robustness. In particular, while the methods from the literature are better at disturbance rejection, the new nonlinear controllers have overall better robustness properties. In general, there is still a need to characterize the trade-off between these two properties. In addition to the results presented in this paper, there is still much room for improvement and many more areas to be explored. Some topics that could be explored in future work include analytical solutions to higher-order terms in the approximation and the effects of varyingγ on the closed-loop response. The need for methods to calculate the gap metric between nonlinear systems was also motivated.



---

*Source: 136315-2014-02-20.xml*