# Computing Bounds on $L_{\infty}$-induced Norm for Linear Time Invariant
Systems Using Homogeneous Lyapunov Functions

Hassan Abdelraouf1, Gidado-Yisa Immanuel2 and Eric Feron3 *This work was
supported by King Abdullah University of Science and Technology (KAUST)1Ph.D.
candidate, Mechanical Engineering, KAUST, Thuwal, Saudi Arabia
<EMAIL_ADDRESS>candidate , School of Aerospace
Engineering, Georgia Institute of Technology, Atlanta, 30332, USA
<EMAIL_ADDRESS>Department of Computer, Electrical and
Mathematical Sciences and Engineering, KAUST, Thuwal, Saudi Arabia.
<EMAIL_ADDRESS>
###### Abstract

Quadratic Lyapunov function has been widely used in the analysis of linear
time invariant (LTI) systems ever since it has shown that the existence of
such quadratic Lyapunov function certifies the stability of the LTI system. In
this work, the problem of finding upper and lower bounds for the
$L_{\infty}$-induced norm of the LTI system is considered. Quadratic Lyapunov
functions are used to find the star norm, the best upper on the
$L_{\infty}$-induced norm, by bounding the unit peak input reachable sets by
inescapable ellipsoids. Instead, a more general class of homogeneous Lyapunov
functions is used to get less conservative upper bounds on the
$L_{\infty}$-induced norm and better conservative approximations for the
reachable sets than those obtained using standard quadratic Lyapunov
functions. The homogeneous Lyapunov function for the LTI system is considered
to be a quadratic Lyapunov function for a higher-order system obtained by
Lifting the LTI system via Kronecker product. Different examples are provided
to show the significant improvements on the bounds obtained by using
Homogeneous Lyapunov functions.

##  I INTRODUCTION

For single-input,single-output (SISO) continuous-time invariant systems, the
${L}_{\infty}$-induced norm, the signal peak-to-peak gain of a transfer
function, is the $\ell_{1}$ norm that can be computed by
$\int_{0}^{\infty}|h(t)|\text{dt}$, where $h(t)$ is the impulse response of
the system. Shamma [1] shows that there does not exist a closed form
expression for the $\ell_{1}$ norm of the LTI system. Given the state space
model of the system, an approximate value of the $\ell_{1}$ norm can be
obtained by only simulation. The problem is to find the final time that
provides an accurate approximation for the integration value. This fact
motivates the search for alternative methods to get reliable upper bounds for
the $\ell_{1}$ norm. Authors in [2] introduce the star-norm to provide a valid
upper bound for the $\ell_{1}$ norm by the construction of inescapable
ellipsoids using quadratic lyapunov functions. But the star-norm calculated in
[2] represents a very conservative upper bound on the $\ell_{1}$ norm for
stiff systems, such as introduced in [3]. This motivated us to introduce new
methods for computing more conservative approximations for the $\ell_{1}$ norm
of linear time invariant systems.

Several methods in the control theory literature are used in computing
$\ell_{1}$ norm including the originating article [4] in which the author
computes the $\ell_{1}$ norm for both discrete and continuous time systems
with some initial results on controller design for disturbance rejection. In
[5], a comparison between the $L_{\infty}$ induced norm of a discrete system
and its RMS gain (the maximum magnitude of the its frequency response) is made
to see how this difference affects $\mathbf{H}_{\infty}$ optimal controller
design. An improved upper and lower bounds for a discrete time system ’s
$\ell_{1}$ norm are introduced in [6] . Based on these bounds, the worst case
$L_{\infty}$ induced norm of the discrete-time systems with diagonal
perturbations can be computed accurately. The $\ell_{1}$ optimal controller
design for discrete-time MIMO systems is introduced in [7]. The controller is
designed to make the system internally stable and optimally track a bounded
input signal. The results in [7] are extended to $\ell_{1}$ optimal controller
design for continuous time systems in [8]. The problem of $\ell_{1}$ optimal
controller with full state feedback is considered in [9] where author showed
that instead of using state feedback dynamic linear controller with arbitrary
higher order, the $\ell_{1}$ optimal controller can be nonlinear memoryless
static feedback.

In this note, we introduce different techniques in approximating the
$\ell_{1}$ norm of LTI systems. The method used in [2] depends on computing
the star-norm by using quadratic Lyapunov function to bound the reachable
sets. This methods produces good bounds for some systems. But for stiff
systems, as we will show later, these bounds become more conservative. Our
approach is based on using higher order homogeneous Lyapunov function in
generating the inescapable ellipsoids to reduce the conservatism in the
computed upper bound of the $\ell_{1}$ norm. This higher order Lyapunov
function can be considered as a quadratic Lyapunov function for a higher order
”lifted” system. The authors of [10] show how the lifted systems are generated
for linear time varying systems using a recursive algorithm based on the
Kronecker product. In addition, they show that the higher order homogeneous
Lyapunov function that certifies stability for linear varying can be
considered as a quadratic Lyapunov function for the lifted system to a higher
degree. The homogeneous Lyapunov functions are also used to obtain better
approximations for some performance metrics for linear time varying systems in
[11]. The authors show that the bounds on system’s peak norms obtained using
higher order homogeneous Lyapunov function are more accurate and less
conservative than the bounds resulted using quadratic Lyapunov functions. In
[11], the effect of control inputs is considered in the lifting process. In
this work, we extend the work introduced in [10] and [11] and use homogeneous
Lyapunov function to compute less conservative upper bounds for the linear
time invariant system’s $\ell_{1}$ norm. These bounds are less than the bounds
introduced in [2]. In this work, the reachable sets for LTI systems using unit
peak input are approximated using higher order homogeneous Lyapunov functions,
then the star norm is computed based on the generated ellipsoid to get better
approximation for the system’s $\ell_{1}$ norm. Additionally, we introduce
another technique to better approximate the integral value
$\int_{0}^{\infty}|Ce^{At}B|\text{dt}$ by calculating the integral to a
specific time $T_{0}$ then approximate an upper bound for the remaining
integral $\int_{T_{0}}^{\infty}|Ce^{At}B|\text{dt}$ by computing the star-norm
for a an equivalent new system. The upper bound for the remaining integral is
also improved by lifting the new system. We show that it is more powerful that
computing the star norm of the original system directly. To illustrate the
results, we compute upper bounds on $\ell_{1}$ norm for different types of
systems like: systems with high damping, systems with low damping and systems
with stiff mass matrix, by using our proposed approaches that gives very
accurate approximations compared with the methods introduced in the
literature.

##  II Notation

Denote the sets of non-negative and positive integers by $\mathbb{Z}_{+}$ and
$\mathbb{Z}_{++}$ respectively. $\mathbb{R}_{+}$ and $\mathbb{R}_{++}$ denotes
the set of non-negative and positive real numbers respectively. The set of
positive definite $n\times n$ matrices is denoted by
$\mathbb{S}_{++}^{n}\subset\mathbb{R}^{n\times n}$. For
$P\in\mathbb{R}^{n\times n}$, $P\succ 0$ means that $P$ is a positive definite
matrix and the function $V(x)=x^{T}Px$ is positive for all non-zero $x\in
R^{n}$. The $n\times n$ identity matrix is denoted by $I_{n}$.

###  II-A Kroncker product

For matrices $A\in\mathbb{R}^{n\times m}$ and $B\in\mathbb{R}^{p\times q}$,
the Kroncker product of $A$ and $B$ is denoted by $A\otimes
B\in\mathbb{R}^{np\times mq}$ and is given by

$A\otimes B:=\left[\begin{array}[]{ccc}a_{11}B&\cdots&a_{1m}B\\\
\vdots&\ddots&\vdots\\\ a_{n1}B&\cdots&a_{nm}B\end{array}\right]$ (1)

where $A$ matrix entries are represented via subscript. We will introduce some
properties of the kroncker product defined in [12] and [13]. For
$A\in\mathbb{R}^{n\times m}$, the $d\textsuperscript{th}$ kroncker power
$A^{\otimes d}$ for all $d\in\mathbb{Z}_{++}$ is defined recursively with a
base $A^{\otimes 0}=1$ by:

$\displaystyle A^{\otimes 1}$ $\displaystyle=A$ (2) $\displaystyle A^{\otimes
d}$ $\displaystyle=A\otimes A^{\otimes(d-1)},\quad d=2,3,\dots$

The following important properties of kroncker product are used in this work:

$\displaystyle A\otimes(B+C)=A\otimes B+A\otimes C$ (3a)
$\displaystyle(A\otimes B)(C\otimes D)=(AC)\otimes(BD)$ (3b)
$\displaystyle(AB)^{\otimes n}=A^{\otimes n}B^{\otimes n}$ (3c)

where all matrices $A$,$B$,$C$ and $D$ are with proper dimensions that permit
the formation of products $AC$ and $BD$. In [12],corollary (7) states that
property (3b) can be generalized as

$A_{1}B_{1}\otimes A_{2}B_{2}\otimes\cdots\otimes A_{n}B_{n}\\\
=\left(A_{1}\otimes A_{2}\otimes\cdots\otimes A_{n}\right)\left(B_{1}\otimes
B_{2}\otimes\cdots\otimes B_{n}\right),$ (4)

such that all matrices $(A_{1},A_{2},\dots A_{n})$ and $(B_{1},B_{2},\dots
B_{n})$ are with proper dimensions that allows all multiplications in the
right hand side.

The $d\textsuperscript{th}$ kroncker sum of a symmetric matrix
$A\in\mathbb{R}^{n\times n}$ is defined by

$A^{\oplus d}=\sum_{k=1}^{d}I_{n}^{\otimes(k-1)}\otimes A\otimes
I_{n}^{(d-k)}.$ (5)

For example,

$\displaystyle A^{\oplus 3}$ $\displaystyle=I_{n}^{\otimes 0}\otimes A\otimes
I_{n}^{\otimes 2}+I_{n}^{\otimes 1}\otimes A\otimes I_{n}^{\otimes
1}+I_{n}^{\otimes 2}\otimes A\otimes I_{n}^{\otimes 0}$ (6)
$\displaystyle=A\otimes I_{n}\otimes I_{n}+I_{n}\otimes A\otimes
I_{n}+I_{n}\otimes I_{n}\otimes A.$

###  II-B Norms

For a vector $v\in\mathbb{R}^{n}$, we define the ${L}_{\infty}$ norm as

$\|v\|_{\infty}=\sup_{i}\left|v_{i}\right|$ (7)

where $v=[v_{1},v_{2},\dots,v_{n}]^{T}$. The function space
$\mathcal{L}_{\infty}$ is defined as the set of functions $f(t)$ for which the
$L_{\infty}$ norm:

$\|f\|_{\infty}=\sup_{t\in\mathbb{R}+}\|f(t)\|_{\infty}$ (8)

is bounded. The $L_{\infty}$-induced norm of an operator $\mathcal{G}$ which
maps $u$ to $y=\mathcal{G}u$ is defined by

$\|\mathcal{G}\|_{\infty}=\sup_{\underset{u\neq
0}{u\in\mathcal{L}_{\infty}}}\frac{\|\mathcal{G}u\|_{\infty}}{\left\|u\right\|_{\infty}}$
(9)

Given that the operator $\mathcal{G}$ is considered as the following single-
input single-output linear time invariant system

$\mathcal{G}:\left\\{\begin{array}[]{l}\dot{x}=Ax+Bu\\\
y=Cx\end{array}\right.$ (10)

where $x\in\mathbb{R}^{n}$ and ($A\in\mathbb{R}^{n\times
n},B\in\mathbb{R}^{n\times 1}$ and $C\in\mathbb{R}^{1\times n}$ ) are the
system’s matrices, the $L_{\infty}$-induced norm in (9) is the $\ell_{1}$ norm
of the system in (10) which can be computed as

$\|\mathcal{G}\|_{1}=\int_{0}^{\infty}|h(t)|\text{dt},$ (11)

where $h(t)$ is the impulse response of the system and can be computed as:

$h(t)=Ce^{At}B$ (12)

##  III Dynamical system lifting

In this section, we incorporate the results introduced in [10] and [11] on
lifting the system dynamics from the original space to a higher order space.
Additionally, we consider the control input in the lifting procedure.

###  III-A Lifting procedure

First, given the LTI system in (10), the state vector $x$ is lifted to
$x^{\otimes d}$ where $d$ is the lifting order. Then, the chain rule is
applied to obtain the derivative of the lifted state vector $x^{\otimes d}$
with respect to time as

$\displaystyle\frac{d}{dt}x^{\otimes d}$
$\displaystyle=\sum_{k=1}^{d}x^{\otimes{(k-1)}}\otimes\dot{x}\otimes
x^{\otimes(d-k)}$ (13)
$\displaystyle=\sum_{(k=1)}^{d}x^{\otimes{k-1}}\otimes(Ax+Bu)\otimes
x^{\otimes(d-k)}.$

By using property (3a), equation (13) can be written as

$\displaystyle\frac{d}{dt}x^{\otimes d}$
$\displaystyle=\sum_{k=1}^{d}x^{\otimes{(k-1)}}\otimes Ax\otimes
x^{\otimes(d-k)}$ (14) $\displaystyle+\sum_{k=1}^{d}x^{\otimes{(k-1)}}\otimes
Bu\otimes x^{\otimes(d-k)}.$

Since $x=I_{n}x$, property (3c) can be used. Therefore, (14) becomes

$\displaystyle\frac{d}{dt}x^{\otimes d}$
$\displaystyle=\sum_{k=1}^{d}(I_{n}^{\otimes{(k-1)}}x^{\otimes{(k-1)}})\otimes
Ax\otimes(I_{n}^{\otimes(d-k)}x^{\otimes(d-k)})$ (15)
$\displaystyle+\sum_{k=1}^{d}(I_{n}^{\otimes{(k-1)}}x^{\otimes{(k-1)}})\otimes
Bu\otimes(I_{n}^{\otimes(d-k)}x^{\otimes(d-k)}),$

hence property (4) can be applied in (15) to be

$\displaystyle\frac{d}{dt}x^{\otimes d}$
$\displaystyle=\sum_{k=1}^{d}(I_{n}^{\otimes(k-1)}\otimes A\otimes
I_{n}^{\otimes(d-k)})x^{\otimes d}$ (16)
$\displaystyle+\sum_{k=1}^{d}(I_{n}^{\otimes(k-1)}\otimes B\otimes
I_{n}^{\otimes(d-k)})ux^{\otimes(d-1)}.$

By using (6), (16) can be summarized as

$\frac{d}{dt}x^{\otimes d}=A^{\oplus d}x^{\otimes d}+B^{\oplus
d}ux^{\otimes(d-1)}.$ (17)

The output equation in (10) is lifted as

$y^{\otimes d}=C^{\otimes d}x^{\otimes d}.$ (18)

To summarize this section, the LTI system (10) can be lifted to a higher order
space with a degree $d\in\mathbb{Z}_{++}$ as

$\displaystyle\dot{\zeta}$ $\displaystyle=A^{\oplus d}\zeta+B^{\oplus d}w$
(19) $\displaystyle\eta$ $\displaystyle=C^{\otimes d}\zeta,$

where $\zeta=x^{\otimes d}\in\mathbb{R}^{n^{d}}$, is the lifted state vector,
$w=ux^{\otimes(d-1)}\in\mathbb{R}^{n^{(d-1)}}$ is the lifted input vector and
$\eta=y^{\otimes d}\in\mathbb{R}$. $A^{\oplus d}\in\mathbb{R}^{n^{d}\times
n^{d}}$, $B^{\oplus d}\in\mathbb{R}^{n^{d}\times n^{d-1}}$ and $C^{\otimes
d}\in\mathbb{R}^{1\times n^{d}}$ are the lifted matrices characterizing the
lifted system dynamics.

##  IV star norm and inescapable ellipsoids ($d=1$)

The authors in [3] and [2] avoided the complexity on computing the
$L_{\infty}$-induced norm of the LTI systems by computing the star norm, an
upper bound on the $L_{\infty}$-induced norm ($\ell_{1}$ norm). The star norm
is obtained by approximating the reachable sets with unit peak input with
inescapable ellipsoids which are defined as follows.

######  Definition 1.

(Reachable set with unit peak input) [3] is defined as the set of all
reachable states from the origin in finite time by unit peak input
$\|u\|_{\infty}\leq 1$

$\mathcal{R}_{\mathrm{up}}\triangleq\left\\{x(T):\begin{array}[]{lr}\dot{x}=Ax+Bu,\quad
x(0)=0\\\ u^{T}u\leq 1,\quad T\geq 0\end{array}\right\\}$ (20)

######  Definition 2.

(Inescapable set) [2] A set $\mathcal{X}$ is said to be inescapable if (1) the
origin ($x(0)\in\mathcal{X}$) and (2) for $x(0)\in\mathcal{X}$ and $u^{T}u\leq
1$, $x(t)$ will evolve inside $\mathcal{X}$ for all future time $t\geq 0$

######  Theorem 1.

[2] consider the stable LTI system (10). Let $P\in\mathbb{R}^{n\times n}$ be a
positive semi-definite. The closed ellipsoid $\\{x:x^{T}Px\leq 1\\}$ is an
inescapable if and only if there exist $\alpha\in R_{+}$ such that:

$\left[\begin{array}[]{cc}A^{T}P+PA+\alpha P&PB\\\ B^{T}P&-\alpha
I\end{array}\right]\leq 0$ (21)

By using Schur complement, and letting $Q=P^{-1}$, The LMI (21) becomes
$AQ+QA^{T}+\alpha Q+BB/\alpha\leq 0$. So, there exists a unique solution if
$(A+\alpha I_{n}/2)$ is a stable matrix. Therefore, $\alpha\in(0,\kappa)$,
where $\kappa=-2\text{max}(\text{real}(\text{eig}(A)))$. For given $\alpha$,
(21) can be solved to get the equivalent inescapable ellipsoid. The objective
is to minimize the maximum output $\|Cx\|_{\infty}$ inside the ellipsoid
$x^{T}Px\leq 1$, so the procedure is as follows:

  1. 1.

For $\alpha$ sweeps from zero to $\kappa$, Solve the following semi-definite
program.

$\displaystyle\underset{P}{\text{minimize}}$ $\displaystyle CP^{-1}C^{T}$ (22)
subject to $\displaystyle P>0$ $\displaystyle\begin{bmatrix}PA+A^{T}P+\alpha
P&PB\\\ B^{T}P&-\alpha I\end{bmatrix}\leq 0$

  2. 2.

For each $\alpha$, compute the upper bound of the output $Cx$ inside the
corresponding inescapable ellipsoid $x^{T}Px\leq 1$, where $P$ is the solution
of (22).

$N_{\alpha}=\underset{x^{T}Px\leq
1}{\text{sup}}||Cx||=||CP^{-1}C^{T}||^{\frac{1}{2}}$ (23)

  3. 3.

Then, the star norm is the lowest of the upper bounds

$H^{\star}=\underset{\alpha\in(0,k)}{\text{inf}}N_{\alpha}.$ (24)

The star norm is the least conservative upper bound determined by inescapable
ellipsoids.

##  V Star norm and inescapable ellipsoids ($d=2$)

In this section, we introduce the main contribution of our work. The LTI
system (10) is lifted to a higher order space such that $d=2$. From (19), the
lifted system is

$\displaystyle\dot{\zeta}$ $\displaystyle=A^{\oplus 2}\zeta+B^{\oplus 2}w$
(25) $\displaystyle\eta$ $\displaystyle=C^{\otimes 2}\zeta,$

such that $\zeta=x^{\otimes 2}\in\mathbb{R}^{n^{2}}$, $w=ux\in\mathbb{R}^{n}$
and $\eta=y^{\otimes d}\in\mathbb{R}$. The matrices
$A\in\mathbb{R}^{n^{2}\times n^{2}}$, $B\in\mathbb{R}^{n^{2}\times n}$ and
$C\in\mathbb{R}^{1\times n^{2}}$ can be written as

$\displaystyle A^{\oplus 2}$ $\displaystyle=A\otimes I_{n}+I_{n}\otimes A$
(26) $\displaystyle B^{\oplus 2}$ $\displaystyle=B\otimes I_{n}+I_{n}\otimes
B$ $\displaystyle C^{\otimes 2}$ $\displaystyle=C\otimes C.$

######  Lemma 2.

If $A$ is a stable matrix, then $A^{\oplus 2}$ is also a stable matrix.

###### Proof.

By using Theorem 4.4.5 in [14]. Let the eignvalues of $A$ are
$\lambda_{1},\lambda_{2},\dots,\lambda_{n}$ and the corresponding eignvectors
are $v_{1},v_{2},\dots,v_{n}$, then the eignvalues of $A^{\oplus 2}$ are the
sum of each pair of the eignvalues of $A$. So,
$\sigma(A)=\\{\lambda_{i}+\lambda_{j}:i=1,2,\dots,n\text{ and }j=1,2,\dots
n\\}$ and the corresponding eignvector of each sum is $v_{i}\otimes v_{j}$.
since $A$ is a stable matrix, so the real part of each $\lambda$ is negative.
Therefore the real part of the sum of any pairs is also negative. So, the real
parts of all eignvalues of $A^{\oplus 2}$ are negative which implies the
stability of $A^{\oplus 2}$. ∎

To clarify the idea of this section which is the main contribution of this
work, we consider a second-order LTI system. Then, the idea can be generalized
to higher order systems. Let $x\in\mathbb{R}^{2}$, then the lifted dynamics
(25) will be

$\displaystyle\dot{\zeta}$
$\displaystyle=\begin{bmatrix}2a_{11}&a_{12}&a_{12}&0\\\
a_{21}&a_{11}&a_{22}&a_{12}\\\ a_{21}&a_{11}&a_{22}&a_{12}\\\
0&a_{21}&a_{21}&2a_{22}\end{bmatrix}\zeta+\begin{bmatrix}2b_{1}&0\\\
b_{1}&b_{2}\\\ b_{1}&b_{2}\\\ 0&2b_{2}\end{bmatrix}w$ (27) $\displaystyle\eta$
$\displaystyle=\begin{bmatrix}c_{1}^{2}&c_{1}c_{2}&c_{1}c_{2}&c_{2}^{2}\end{bmatrix}\zeta,$

where,
$\zeta=\begin{bmatrix}x_{1}^{2}&x_{1}x_{2}&x_{2}x_{1}&x_{2}^{2}\end{bmatrix}^{T}$,
$w=\begin{bmatrix}ux_{1}&ux_{2}\end{bmatrix}^{T}$, and $\eta=y^{2}$. The
objective to get the inescapable ellipsoids which approximates the reachable
set with unit peak input $u^{2}\leq 1$ and the star norm that is considered as
an upper bound on the $L_{\infty}$\- induced norm.

The condition $u^{2}\leq 1$ imposes new inequality constraints on the lifted
states $\zeta$ and lifted inputs $w$ as follows

$\displaystyle w_{1}^{2}$ $\displaystyle=u^{2}x_{1}^{2}\leq
x_{1}^{2}=\zeta_{1}$ (28) $\displaystyle w_{2}^{2}$
$\displaystyle=u^{2}x_{2}^{2}\leq x_{2}^{2}=\zeta_{4}.$

There also exists an equality constraint, $w_{1}x_{2}-w_{2}x_{1}=0$, but this
equality is a function of $w$ and original state $x$. So, we multiply both
sides by $x_{1}$ and $x_{2}$ to get two different quadratic equality
constraints in $\zeta$ and $w$ given by

$\displaystyle w_{1}\zeta_{2}-w_{2}\zeta_{1}$ $\displaystyle=0$ (29)
$\displaystyle w_{1}\zeta_{4}-w_{2}\zeta_{3}$ $\displaystyle=0.$

Suppose there exist a quadratic function $V(\zeta)=\zeta^{T}P\zeta$ where
$P\in\mathbb{S}_{++}^{4}\succ 0$ and $dV(\zeta)/dt\leq 0$ for all $\zeta$ and
$w$ satisfying (27) whenever $\zeta^{T}P\zeta\geq 1$ and $w,\zeta$ satisfy
(28) and (29). Then, the ellipsoid $\zeta^{T}P\zeta\leq 1$ is an inescapable
ellipsoid and contains the reachable set. Therefore, these conditions can be
written as

$\displaystyle\zeta^{T}(PA^{\oplus 2}+{A^{\oplus
2}}^{T}P)\zeta+\zeta^{T}PB^{\oplus 2}w+w^{T}{B^{\oplus 2}}^{T}P\zeta\leq 0$
(30) $\displaystyle\forall\zeta\text{ and }w\text{ whenever,
}\zeta^{T}P\zeta\leq 1,\text{ and }$
$\displaystyle\qquad\qquad\qquad\qquad\qquad(\ref{inequality constraints
n=2})\text{ and }(\ref{equality w and zeta})\text{ are satisfied}.$

Using $\mathcal{S}$ procedure [3], condition (30) holds if there exists
$\alpha\geq 0$, $\beta_{1}\geq 0$, $\beta_{2}\geq 0$, $\gamma_{1}$ and
$\gamma_{2}$ such that for all $\zeta$ and $w$:

$\begin{bmatrix}PA^{\oplus 2}+{A^{\oplus 2}}^{T}P+\alpha P&PB^{\oplus
2}+\gamma_{1}E_{1}+\gamma_{2}E_{1}&z\\\ {B^{\oplus
2}}^{T}P+\gamma_{1}E_{1}^{T}+\gamma_{2}E_{2}^{T}&W&\mathbf{0}_{2}\\\
z^{T}&\mathbf{0}_{2}^{T}&-\alpha\end{bmatrix}\leq 0$ (31)

where $\mathbf{0}_{2}\in\mathbb{R}^{2}$ is the zero vector in
$\mathbb{R}^{2}$,
$z=\begin{bmatrix}\beta_{1}/2&0&0&\beta_{2}/2\end{bmatrix}^{T}$, and the
matrices $E_{1}$, $E_{2}$ and W are

$\displaystyle E_{1}=\begin{bmatrix}0&-1\\\ 1&0\\\ 0&0\\\
0&0\end{bmatrix},\quad E_{1}=\begin{bmatrix}0&0\\\ 0&0\\\ 0&-1\\\
1&0\end{bmatrix},\quad W=\begin{bmatrix}-\beta_{1}&0\\\
0&-\beta_{2}\end{bmatrix}.$ (32)

The same procedure used in section IV to get the inescapable ellipsoid and the
star norm for the lifted system. First, the following SDP is solved for every
$\alpha\in(0,\kappa)$ , where
$\kappa=-2\text{max}(\text{real}(\text{eig}(A^{\oplus 2})))$

$\displaystyle\underset{P}{\text{minimize}}$ $\displaystyle C^{\otimes
2}P^{-1}{C^{\otimes 2}}^{T}$ (33) subject to $\displaystyle P>0,\quad(\ref{LMI
n=2 d=2}).$

Then, find the upper bound of the output in every inescapable ellipsoid
$x^{\otimes 2}Px^{\otimes 2}\leq 1$ as

$N_{\alpha}=||C^{\otimes 2}P^{-1}{C^{\otimes 2}}||^{\frac{1}{2}}.$ (34)

The star norm for the lifted system is the smallest of these upper bounds.

$H_{2}^{\star}=\underset{\alpha\in(0,k)}{\text{inf}}N_{\alpha}$ (35)

Finally, the star norm of the original system can be calculated as the square
root of the star norm obtained for the higher order system ($d=2$).

The lifting procedure can be done for higher order systems ($n>2$) using the
same procedure, but the number of equality and inequality constraints imposed
in the $\mathcal{S}$ procedure will increase. In this case, we have $n$
inequality constraints $\\{w_{i}^{2}\leq\zeta_{i(n+1)-n}:i=1\dots n\\}$. The
number of equality constraints is $n(n-1)/2$:

$\begin{matrix}w_{1}x_{2}-w_{2}x_{1}=0,\dots,w_{1}x_{n}-w_{n}x_{1}=0,\\\
w_{2}x_{3}-w_{3}x_{2}=0,\dots,w_{2}x_{n}-w_{n}x_{2}=0,\\\ \dots
w_{n-1}x_{n}-w_{n}x_{n-1}=0\end{matrix}$ (36)

As in (29) each of both sides of these equality constraints is multiplied by
$x_{1},x_{2},\dots,x_{n}$ to come up with $n^{2}(n-1)/2$ quadratic equality
constraints in $w$ and $\zeta$ that can be imposed easily in the $\mathcal{S}$
procedure.

##  VI Numerical experiments

In this section, we consider three systems: high damping, low damping and
stiff systems. To show the effect of lifting the dynamical system and using
homogeneous Lyapunov function in approximating the reachable sets and
computing upper bounds on the $\ell_{1}$ norm.

###  VI-A System with high damping

Consider the following LTI system

$\dot{x}=\begin{bmatrix}0&1\\\ -4&-4\end{bmatrix}x+\begin{bmatrix}0\\\
1\end{bmatrix}u,\qquad y=\begin{bmatrix}1&1\end{bmatrix}x$ (37)

The exact value of the $\ell_{1}$ norm is obtained using Mathematica by
computing $\int_{0}^{\infty}|Ce^{At}B|dt=\mathbf{0.3177}$. By using method
introduced in section IV without lifting ($d=1$), the computed star norm is
$\mathbf{0.3536}$. After lifting the system to $d=2$, the star norm computed
using the method in section V is $\mathbf{0.3368}$ which is a better
approximation for the upper bound of the $\ell_{1}$ norm. Additionally, in
Fig. 1, the blue and black ellipsoids represents the approximate of the unit
peak input reachable sets using ($d=1$) and ($d=2$) respectively. Which shows
that using the higher order system produces better approximation for the
reachable set. The red line in Fig. 1 represents an attempt at computing a
worst case trajectory, the trajectory generated using the control input $u$
that maximize the gradient of the Lyapunov function ($dV/dt$) at every time
step. this control input can be written as

$\displaystyle u$ $\displaystyle=\underset{u^{2}\leq 1}{\text{argmax}}\quad
x^{T}(PA+A^{T}P)x+2x^{T}PBu$ (38) $\displaystyle=\text{sign}(x^{T}PB)$

Figure 1: Approximate unit peak peak input reachable sets and worst case
trajectory for the high damping system (37) .

The peak output of the worst case trajectory represents the lower bound for
the $\ell_{1}$ norm Fig. 2 shows that value of this lower bound is
$\mathbf{0.3097}$

Figure 2: Red curve represents the ”Worst case” trajectory and the blue and
black dotted lines are the star norm for $d=1$ and $d=2$ respectively for the
high damping system (37)

###  VI-B System with low damping

Consider the following low damping system

$\dot{x}=\begin{bmatrix}0&1\\\ -0.5&-0.5\end{bmatrix}x+\begin{bmatrix}0\\\
1\end{bmatrix}u,\qquad y=\begin{bmatrix}1&1\end{bmatrix}x$ (39)

The exact value of the $\ell_{1}$ norm is by calculating the integral
$\int_{0}^{\infty}|Ce^{At}B|dt$ is $\mathbf{4.3069}$. The computations shows
that the value of the star norm at ($d=2$), $\mathbf{4.5533}$, is less than
the value of star norm obtained at ($d=1$), $\mathbf{4.63}$. In addition, Fig.
3 shows that the approximation of the unit peak input reachable set is less
conservative when using the lifted system ($d=2$).

Figure 3: Approximate unit peak peak input reachable sets and worst case
trajectory for low damping system (39) .

###  VI-C stiff system

Consider the following system

$\dot{x}=\begin{bmatrix}-1&0\\\ 0&-100\end{bmatrix}x+\begin{bmatrix}1\\\
100\end{bmatrix}u,\qquad y=\begin{bmatrix}1&-2\end{bmatrix}x$ (40)

The exact value of the $\ell_{1}$ norm is $\mathbf{3.0412}$. The star norm
computed when ($d=1$) is $\mathbf{10.4600}$ which is very conservative upper
bound for the $\ell_{1}$ norm. However, the star norm computed using the
lifted system is $\mathbf{5.7680}$ which shows the effectiveness of lifting
the system to obtain better bounds. Also, Fig. 4 shows the significant effect
of lifting in reducing the conservatism in reachable set approximations for
system (40).

Figure 4: Approximate unit peak peak input reachable sets and worst case
trajectory for stiff system (40) .

##  VII Alternative method to approximate $\ell_{1}$ norm

In this section, we introduce another method to get better and more accurate
approximation for the $\ell_{1}$ norm which is
$\int_{0}^{\infty}|Ce^{At}B|dt$. This integral is divided in two parts as:

$\int_{0}^{\infty}|Ce^{At}B|dt=\int_{0}^{T_{0}}|Ce^{At}B|dt+\int_{T_{0}}^{\infty}\left|Ce^{At}B\right|dt$
(41)

where $T_{0}>0$ is known, so the first term can be computed exactly. The
second term is the $\ell_{1}$ norm of the following system:

$\displaystyle\dot{z}$ $\displaystyle=Az+e^{AT_{0}}Bu$ (42) $\displaystyle y$
$\displaystyle=Cz$

An upper bound on the $\ell_{1}$ norm of system (42) can be approximated as
the star norm where the system can be lifted to a higher order to obtain
better approximation. Table I shows the upper bound of the $\ell_{1}$ norm of
the low damping system (39) using different values of $T_{0}$. At each
$T_{0}$, the first term in (41) is computed and the second term is
approximated at $d=1$ and $d=2$ for the system (42).

TABLE I: $\ell_{1}$ norm approximation for system (39) [exact value is $4.3069$] $T_{0}\text{(}Sec)$ | Upper bound ($d=1$) | Upper bound ($d=2$)
---|---|---
2 | 4.5683 | 4.5304
5 | 4.4078 | 4.3981
10 | 4.3376 | 4.3332
20 | 4.3096 | 4.3091

This method provides more accurate approximations even for stiff system (40),
The upper bound obtained by using $T_{0}=0.05\text{sec}$ and $d=2$ is
$\mathbf{3.0424}$ which is very close to the exact value $\mathbf{3.0412}$

##  VIII CONCLUSIONS

This work demonstrates how the polynomial homogeneous Lyapunov function can be
used to obtain better approximations for the LTI system’s star norm which is
an upper bound for the system’s $\ell_{1}$ norm. We also demonstrated that the
inescapable ellipsoids, approximations for the unit peak input reachable sets,
generated using higher order Lyapunov functions is more accurate and less
conservative that those produced by using quadratic Lyapunov functions. We
tested that for three types of systems: systems with high damping, systems
with low damping and stiff systems, to show the improvement of the reachable
sets approximation. We only considered lifting the LTI system $d=2$. For
$d>2$, the main difficulty is to transfer the inequality $u^{2}\leq 1$ in the
original space to quadratic inequalities in the higher order space as in (28)
for $d=2$. As a future work, we are working to develop new techniques to deal
with this difficulty. Also, there are future opportunities to use homogeneous
Lyapunov functions in robustness analysis and design better-performing
controllers.

## References

  * [1] J. S. Shamma, “Nonlinear state feedback for $\ell_{1}$ optimal control,” _Systems & control letters_, vol. 21, no. 4, pp. 265–270, 1993.
  * [2] J. Abedor, K. Nagpal, and K. Poolla, “A linear matrix inequality approach to peak-to-peak gain minimization,” _International Journal of Robust and Nonlinear Control_ , vol. 6, no. 9-10, pp. 899–927, 1996.
  * [3] E. M. Feron, “Linear matrix inequalities for the problem of absolute stability of control systems,” Ph.D. dissertation, Stanford University, 1994.
  * [4] M. Vidyasagar, “Optimal rejection of persistent bounded disturbances,” _IEEE Transactions on Automatic Control_ , vol. 31, no. 6, pp. 527–534, 1986\.
  * [5] S. Boyd and J. Doyle, “Comparison of peak and RMS gains for discrete-time systems,” _Systems & control letters_, vol. 9, no. 1, pp. 1–6, 1987.
  * [6] V. Balakrishnan and S. Boyd, “On computing the worst-case peak gain of linear systems,” _Systems & Control Letters_, vol. 19, no. 4, pp. 265–269, 1992\.
  * [7] M. Dahleh and J. Pearson, “$\ell_{1}$ optimal feedback controllers for $MIMO$ discrete time systems,” _IEEE Transactions on Automatic Control_ , vol. 32, no. 4, pp. 314–322, 1987.
  * [8] ——, “$\ell_{1}$ optimal compensators for continuous-time systems,” _IEEE Transactions on Automatic Control_ , vol. 32, no. 10, pp. 889–895, 1987\.
  * [9] J. S. Shamma, “Optimization of the ${L_{\infty}}$-induced norm under full state feedback,” _IEEE Transactions on Automatic Control_ , vol. 41, no. 4, pp. 533–544, 1996.
  * [10] M. Abate, C. Klett, S. Coogan, and E. Feron, “Lyapunov differential equation hierarchy and polynomial Lyapunov functions for switched linear systems,” in _2020 American Control Conference (ACC)_. IEEE, 2020, pp. 5322–5327.
  * [11] ——, “Pointwise-in-time analysis and non-quadratic Lyapunov functions for linear time-varying systems,” in _2021 American Control Conference (ACC)_. IEEE, 2021, pp. 3550–3555.
  * [12] B. J. Broxson, “The kronecker product,” 2006.
  * [13] C. F. Van Loan, “The ubiquitous kronecker product,” _Journal of computational and applied mathematics_ , vol. 123, no. 1-2, pp. 85–100, 2000.
  * [14] R. A. Horn, R. A. Horn, and C. R. Johnson, _Topics in matrix analysis_. Cambridge university press, 1994\.


21-776

# Sub-Optimal Fast Fourier Series Approximation for Initial Trajectory Design

Caleb Gunsaulus Undergraduate Student, Department of Aerospace Engineering,
Iowa State University, USA.  Carl De Vries Graduate Student, Department of
Aerospace Engineering, Iowa State University, USA.  William Brown
Undergraduate Student, Department of Aerospace Engineering, Iowa State
University, USA.  Youngro Lee Graduate Student, Department of Aerospace
Engineering, Iowa State University, USA.  Madhusudan Vijayakumar Graduate
Student, Department of Aerospace Engineering, Iowa State University, USA.
Ossama Abdelkhalik Associate Professor, Department of Aerospace Engineering,
Iowa State University, Ames, Iowa 50011, AIAA Senior Member.

###### Abstract

The Finite Fourier Series (FFS) Shape-Based (SB) trajectory approximation
method has been used to rapidly generate initial trajectories that satisfy the
dynamics, trajectory boundary conditions, and limitation on maximum thrust
acceleration. The FFS SB approach solves a nonlinear programming problem (NLP)
in searching for feasible trajectories. This paper extends the development of
the FFS SB approach to generate sub optimal solutions. Specifically, the
objective function of the NLP problem is modified to include also a measure
for the time of flight. Numerical results presented in this paper show several
solutions that differ from those of the original FFS SB ones. The sub-optimal
trajectories generated using a time of flight minimization are shown to be
physically feasible trajectories and potential candidates for direct solvers.

##  1 Introduction

Current trajectory optimization methods can minimize objectives such as time
of flight or fuel consumption. Two common methods to solve these problems
include direct and indirect methods [1]. Direct methods are very useful, but
can be computationally expensive, and require initial guesses for the design
parameters to start the solution process. Sub-optimal, yet feasible,
trajectories can be a good option to initialize a direct solver. However, many
sub-optimal methods are useful only in limited scenarios, are subject to
geometric requirements, or are not able to satisfy spacecraft thrust
constraints. One approximate method is the Finite Fourier Series Shape-Based
method [2, 3, 4]. The FFS SB method generates a spacecraft trajectory by
solving a NLP problem for a set of Fourier coefficients which approximate the
spacecraft states as a function of time. Depending on the problem, some of the
Fourier coefficients can be determined from the boundary conditions. The rest
of the unknown coefficients are found using an optimizer. In this case,
MATLAB’s fmincon function is used to find trajectories by minimizing an
objective function. To find these trajectories, the optimizer computes the
error from the FFS SB approximation at multiple points along the trajectory
and then minimizes the sum of the squared error. The FFS SB method has been
developed for two-body and three-body systems as well as two dimensional and
three dimensional spaces [3, 4]. The method has been applied to interplanetary
transfers, spacecraft rendezvous, orbit raising, and phasing maneuvers. While
the FFS SB method doesn’t produce optimal solutions, it was shown to produce
cost values (fuel mass) near the optimal solutions obtained using GPOPS in
some test cases [4].

Currently, the FFS SB method does not attempt to minimize objectives such as
time of flight or fuel mass. In this paper, a time of flight penalty is added
to the FFS SB method objective function. This work solves the multi-objective
problem given that the time of flight objective could be competing with the
objective of satisfying equations of motion. The approach is applied to
spacecraft orbit raising and rendezvous test cases which demonstrate feasible
trajectories with shorter flight times are generated using the FFS SB method.

##  2 Finite Fourier Series Shape-Based Method

Mathematically, any given periodic function can be represented as an infinite
summation of sine and cosine terms. Shape-Based methods belong to a domain of
space trajectory optimization which geometrically approximates the shape of a
trajectory. This trajectory is then used to compute the thrust vector and
hence the trajectory solution of the design/optimization problem. One such
popular and robust representation of the trajectory shape can be obtained by
using the Finite Fourier Series. The FFS SB method approximates the states of
the spacecraft namely the radius $r$ and polar angle $\theta$ as follows:

$r(t)=\frac{a_{0}}{2}+\sum_{n=1}^{n_{r}}[a_{n}cos(\frac{n\pi
t}{T})+b_{n}sin(\frac{n\pi t}{T})]$ (1)

$\theta(t)=\frac{c_{0}}{2}+\sum_{n=1}^{n_{\theta}}[c_{n}cos(\frac{n\pi
t}{T})+d_{n}sin(\frac{n\pi t}{T})]$ (2)

where, $n_{r}$ and $n_{\theta}$ are the number of Fourier coefficients used to
represent the radius $r$ and polar angle $\theta$ of the spacecraft states,
respectively. The biggest advantage of the FFS method is that for every
different selection of the Fourier coefficients, a different trajectory shape
is obtained. Additionally, one should note that the states of the trajectory
are discretized at ’n’ points and are only approximated using the FFS at these
discrete points. If the initial and final boundary conditions (BCs) of the
spacecraft states are available, they can be used to solve for some of the
Fourier coefficients in terms of the others.

The number of Fourier terms to be used is problem dependent and is obtained on
a trial and error basis. Once the number of FFS terms and the number of
discretization points are determined, the sine and cosine terms become fixed
and can be calculated once and used repeatedly to decrease computational cost.
The Fourier approximated function can thus be reduced to the following form:

$r(t)=[A]*X+B$ (3)

where A is a square matrix of order $(2*n_{r}+1)\times(2*n_{r}+1)$ consisting
of the sine and cosine terms, X is a column matrix having ($2*n_{r}+1$) rows
representing the unknown Fourier coefficients that have to be solved and B is
the column matrix having ($2*n_{r}+1$) rows that represents the coefficients
that were solved for using the BCs.

##  3 Equations of Motion

The polar equations of motions are developed in two dimensions where $r$ is
the spacecraft radius and $\theta$ is the polar angle [5]. The spacecraft is
acted on by a gravitational force from a single, central body and an external
thrust input $T_{a}$. A local vertical local horizontal coordinate frame is
fixed to the spacecraft. The angle $\alpha$ is the steering angle measured
from the local horizontal and is the angle at which the thrust acceleration is
applied. The gravitational parameter is given by $\mu$, and the equations of
motion are shown in Eqs. (4) and (5).

$\ddot{r}-r\dot{\theta}^{2}+\frac{\mu}{r^{2}}=T_{a}sin(\alpha)$ (4)
$2\dot{r}\dot{\theta}+r\ddot{\theta}=T_{a}cos(\alpha)$ (5)

The steering angle is commonly equal to the flight path angle, $\gamma$, plus
a multiple of $\pi$ because it’s allowed to fire in both the positive and
negative directions in Eq. (6) [2, 5, 6, 7, 8, 9]. The additional term $\pi n$
is added to give the direction. The thrust acceleration is applied in the
direction of the flight path angle when $n=0$ and opposite to the flight path
angle when $n=1$. The steering angle can also be found from the common flight
path angle equation as the ratio of the velocity components in the local
vertical local horizontal frame shown in Eq. (7).

$\alpha=\gamma+\pi n$ (6)

$tan(\alpha)=tan(\gamma)=\frac{\dot{r}}{r\dot{\theta}}$ (7)

The thrust acceleration shown in Eq. (8) can be solved for from Eq. (5).
Substituting it into Eq. (4) and applying Eq. (7) results in a single equation
of motion shown in Eq. (9).

$T_{a}=\frac{2\dot{r}\dot{\theta}+r\ddot{\theta}}{cos(\alpha)}$ (8)
$F=r^{2}(\dot{\theta}\ddot{r}-\dot{r}\ddot{\theta})+\dot{\theta}(\mu-2r\dot{r}^{2})-(r\dot{\theta})^{3}=0$
(9)

The FFS SB method exploits the fact that the equations of motion can be
formulated in a single expression and be collected easily so all the terms are
set equal to zero. The first and second derivatives of the FFS SB state
approximations can be computed from Eqs. (1) and (2). When the states and
state derivatives are evaluated using specific Fourier coefficients, the
resulting values can be plugged into Eq. (9). The result of evaluating this
equation will be non-zero because they are approximations. This error
quantifies how well the FFS SB approximated values fit the equations of
motion. This calculation is performed at multiple points along the trajectory
with the goal to minimize the summation of the squared residual error.

##  4 Sub-Optimal Solution Formulation

The objective function developed to solve the multi-objective function is
composed of a weighted sum of individual objectives. Two objectives are used
and they are 1) the squared error of the FFS approximation to the equations of
motion and 2) the time of flight penalty. The objective function is shown in
Eq. (10)

$J=(1-\omega)\textbf{F}^{T}\textbf{F}+\omega\left(\frac{ToF}{ToF_{0}}\right)$
(10)

where, $J$ is the objective function, F is the error in the Fourier series
approximation of the equations of motion evaluated multiple points on the
trajectory, $ToF$ is the time of flight, $ToF_{0}$ is the initial time of
flight guess, and $\omega$ is the weight applied to the ToF.

Now that the time of flight is included in the objective function, the goal is
to select a weight which will minimize the combination of the weighted
equation of motion squared error and the weighted time of flight.
Modifications must be made to ensure the trajectory’s boundary conditions are
satisfied since the time of flight is now a design parameter. The final polar
angle is a free parameter in the orbit raising formulation, however,
spacecraft rendezvous requires a fixed polar angle. The final polar angle at
which rendezvous occurs will change throughout the optimization process as the
time of flight varies.

The steps to compute the new rendezvous polar angle are shown in Eqs. (11),
(12), and (13). The orbital rate for the target orbit can be computed from the
gravitational parameter and the semi-major axis for the final, circular orbit.
Then the initial guess for the time of flight is subtracted from the current
ToF trial value. The change in the time of flight is multiplied by the orbital
rate and yields the change in mean anomaly (polar angle) for the target orbit
across that time period. Adding the change in mean anomaly to a reference mean
anomaly which corresponds to the initial guess for the time of flight results
in the new rendezvous polar angle.

$n=\sqrt{\frac{\mu}{r_{f}^{3}}}$ (11)

$\Delta M=n(ToF-ToF_{0})$ (12)

$\theta_{f}=\theta_{f,0}+\Delta M$ (13)

##  5 Optimal Solution Formulation

The FFS generates the sub-optimal thrust acceleration by penalizing ToF in the
objective function, as shown in the previous design process and results. In
order to examine the similarity between the sub-optimal solutions and optimal
solutions, General Purpose OPtimal control Software (GPOPS-II), which is a
MATLAB software intended to solve general nonlinear optimal control problems
[10], is employed to find optimal solutions. Since GPOPS-II is a direct method
to solve the optimal control problem, the equations of motion in Eqs. (4) and
(5) are used instead of using the combined equation, which is Eq. (9).
Accordingly, it does not need to consider the motion error itself, so the
objective function only includes the ToF term.

The minimum time orbit raising problem can be constructed as shown in Eq.
(14).

$\displaystyle J$ $\displaystyle=ToF$ (14a) subject to $\displaystyle
Eqs.\eqref{eq:r_eom}\text{ and }\eqref{eq:theta_eom}$ (14b) $\displaystyle
r(0)=r_{i},\theta(0)=\theta_{i},\dot{r}(0)=\dot{r}_{i},\dot{\theta}(0)=\dot{\theta}_{i}$
(14c) $\displaystyle
r(ToF)=r_{f},\dot{r}(ToF)=\dot{r}_{f},\dot{\theta}(ToF)=\dot{\theta}_{f}$
(14d) $\displaystyle-T_{a,max}\leq T_{a}\leq T_{a,max}$ (14e) $\displaystyle
r_{i}\leq r\leq r_{f},0\leq\theta\leq\theta_{f,0},0\leq\dot{r}\leq
1,0\leq\dot{\theta}\leq 1$ (14f)

where $\theta_{f,0}$, which is the initial guess of the final polar angle, is
set as $12\pi$, and the initial guess of ToF is given as 120000 seconds.
Furthermore, the boundaries of the time derivatives of the radial distance and
the polar angle are defined as proper values that provide enough space to
compute the optimal trajectory.

GPOPS computes the optimal solution for the rendezvous case as well. Unlike
the orbit raising, the rendezvous case constrains the final polar angle, which
can be calculated based on equations (11), (12), and (13). Therefore, the new
rendezvous polar angle is included in Eq. (15).

$\displaystyle J$ $\displaystyle=ToF$ (15a) subject to $\displaystyle
Eqs.\eqref{eq:r_eom}\text{ and }\eqref{eq:theta_eom}$ (15b) $\displaystyle
r(0)=r_{i},\theta(0)=\theta_{i},\dot{r}(0)=\dot{r}_{i},\dot{\theta}(0)=\dot{\theta}_{i}$
(15c) $\displaystyle
r(ToF)=r_{f},\dot{r}(ToF)=\dot{r}_{f},\dot{\theta}(ToF)=\dot{\theta}_{f}$
(15d)
$\displaystyle\theta(ToF)=\theta_{f,0}+\sqrt{\frac{\mu}{r_{f}^{3}}}\left(ToF-
ToF_{0}\right)$ (15e) $\displaystyle-T_{a,max}\leq T_{a}\leq T_{a,max}$ (15f)
$\displaystyle r_{i}\leq r\leq
r_{f},0\leq\theta\leq\theta_{f,0},0\leq\dot{r}\leq 1,0\leq\dot{\theta}\leq
1.2$ (15g)

The variables $\theta_{f,0}$ and $ToF_{0}$ are the initial guess for the final
polar angle and the time of flight, and their values are given as $13\pi$ and
120000 seconds, respectively.

##  6 Test Cases

###  6.1 Orbit Raising

The first test case considered is an orbit raising scenario in which the
spacecraft moves from LEO to GEO, and the final polar angle is free. For
comparison, the same design parameters and boundary conditions from Taheri and
Abdelkhalik are used [2]. The initial LEO and final GEO orbits are circular,
and the boundary conditions are shown in Table 1. The gravitational parameter
for the Earth is $\mu=398601.2$ km3/s2. The angular rates of the spacecraft in
the LEO and GEO orbits are then computed from $\dot{\theta}=\sqrt{\mu/r^{3}}$.
The boundary conditions and parameters are converted to canonical units prior
to simulation. All of the figures which follow are presented using canonical
units.

The number of Fourier terms used to approximate the radial and angular states
as well as the number of equally spaced discrete points (DP) along the
trajectory are shown in Table 2. The maximum and minimum thrust accelerations
allowed throughout the flight are designated as $T_{a,max}$, and the time of
flight weight is $\omega$.

Table 1: Boundary Conditions: Orbit Raising Parameter | Value
---|---
$r_{i}$ | 6570 km
$\theta_{i}$ | 0 deg.
$\dot{r}_{i}$ | 0 km/s
$r_{f}$ | 42160 km
$\dot{r}_{f}$ | 0 km/s
Table 2: Design Parameters: Orbit Raising Case 1 Parameter | Value
---|---
$n_{r}$ | 2
$n_{\theta}$ | 3
$DP$ | 160
$T_{a,max}$ |  $\pm$ 0.0102 DU/TU2
$\omega$ | 0.5

The weight for the time of flight penalty is 50% so both the equation of
motion error and time of flight objectives are weighted equally. The low-
thrust FFS SB trajectory generated from penalizing the ToF is compared against
previously verified results shown below [2]. The trajectory and thrust
acceleration shown in Figure 2 and Figure 2 show the time of flight decreases
as expected when the time of flight is penalized. The trajectory resulted in a
ToF reduction of $\Delta ToF=-3.95$ hours - an 11.84% reduction. However, the
equation of motion error increased which indicates the FFS approximation does
not satisfy the equations of motion as well.

Figure 1: FFS Trajectory: Orbit Raising

Figure 2: Thrust Acceleration: Orbit Raising

The reduction in time of flight (at the cost of additional equation of motion
error) is an improvement, but the time of flight weights which produce
feasible trajectories must still be determined. The results above are extended
to show the trade off between the equation of motion error and the time of
flight objectives by plotting both as functions of the time of flight weight
in Figure 4 and Figure 4. The figures include plots of the objectives for
trajectories using FFS SB parameters $n_{r}=2$ and $n{\theta}=3$ as well as
$n_{r}=3$ and $n_{\theta}=4$. The trajectories were generated for weights
between 0.01 and 0.99 in steps of 0.01.

Figure 4 shows that error increases as the weight increases, which is
expected. At low weights the error in the EoM is very small but grows
continually as the normalized ToF becomes increasingly more important in the
objective function. Another important observation is that increasing the
number of Fourier coefficients to approximate the trajectory decreases the
total error which is shown in Figure 4, and reduces the total time of flight
which is shown in Figure 4. These plots were used to determine the trade-offs
in the objective function for the purpose of selecting better time of flight
weights which satisfy the multi-objective optimization.

Figure 3: Optimal EoM vs. Weight

Figure 4: Optimal ToF vs Weight

####  6.1.1 Trajectory Feasibility

The FFS trajectory results are compared against results generated from
numerically integrating the equations of motion. A weight of of 0.01 was
selected for the following test cases. It is completely viable to choose a
larger weight, however, the drawback is that the error in the equations of
motion will increase. This will result in larger position errors in the
spacecraft trajectory FFS approximation.

The numerical integration was implemented using the built in MATLAB function
ode45. This function is used to solve the system of differential equations for
the spacecraft shown in Eqs. (4) and (5). The simulation uses an open loop
control approach in which the FFS SB thrust acceleration profile is
interpolated at each time step. The design parameters used to obtain the
following results are shown in Table 3.

Table 3: Design Parameters: Orbit Raising Case 2 Parameter | Value
---|---
$n_{r}$ | 3
$n_{\theta}$ | 4
$DP$ | 160
$T_{a,max}$ |  $\pm$ 0.0102 DU/TU2
$\omega$ | 0.01

In Figure 6 two FFS trajectories are shown. For one trajectory the time of
flight weight is $\omega=0$, meaning the only term contributing to the
objective function is the term containing the error in the EoM,
$\textbf{F}^{T}\textbf{F}$. The other trajectory has a time of flight weight
of $\omega=0.01$. Comparing the two trajectories it can be seen that the
spacecraft spends much more time closer to the Earth on all but the final
revolution of the penalized trajectory. On this final revolution, the radial
position increases quickly to the desired radius rather than slowly growing
over multiple revolutions as it does in the non-penalized trajectory.

The thrust acceleration curves for these two trajectories are shown in Figure
6 and it can be seen that the penalized trajectory takes much less time to
reach the desired orbit. The equation of motion error between the FFS
approximation and the ode45 trajectory can be seen in Figure 7. The
trajectories and final polar angles are slightly different, however, the
trajectories are still nearly a match.

Figure 5: FFS Trajectories: Orbit Raising

Figure 6: Thrust Acceleration: Orbit Raising

Figure 7: Trajectory Comparison: Orbit Raising

One final orbit raising trajectory was generated with the parameters shown in
Table 4. This test case uses more $n_{r}$ and $n_{\theta}$ terms. In addition,
the number of discretization points had to be decreased from 160 points to 40
points in order to generate Fourier coefficients which produced a feasible
trajectory. Results are shown in figures 9 and 9. This test case was
specifically included to demonstrate how close the trajectories generated
using FFS can be to the optimal solutions generated by GPOPS.

Table 4: Design Parameters: Orbit Raising Case 3 Parameter | Value
---|---
$n_{r}$ | 5
$n_{\theta}$ | 7
$DP$ | 40
$T_{a,max}$ |  $\pm$ 0.0102 DU/TU2
$\omega$ | 0.01

In addition to the FFS SB and associated ode45 trajectories, figures 9 and 9
show the optimal trajectory and thrust acceleration profile generated using
GPOPS. The optimal flight time is 17.38 hours, which corresponds to a 47.87 %
reduction compared to the FFS trajectory with no ToF weight for case 3. This
is very close to the FFS trajectory time of flight, $ToF=17.91$ hours, for the
FFS parameters in Table 5. Even though the optimal time of flight and the FFS
time of flight are very close, the trajectories are still quite different due
to variance in their thrust acceleration profiles. Given these differences,
the optimal trajectory orbits almost 4 times while the FFS trajectory orbits 5
times. The most distinctive feature of the optimal thrust acceleration
solution is that it switches its value from one extreme to the other, which is
referred to as a bang-bang solution.

Figure 8: FFS and GPOPS Trajectories: Orbit Raising

Figure 9: Thrust Acceleration: Orbit Raising

All results for the orbit raising cases presented above are shown in Table 5.
The GPOPS optimal ToF is compared to the non-penalized ToF result generated in
case 3, $ToF=33.34$ hours.

Table 5: Time of Flight Results: Orbit Raising Parameter | Case 2 | Case 3 | GPOPS
---|---|---|---
$n_{r}$ | 3 | 5 | -
$n_{\theta}$ | 4 | 7 | -
$DP$ | 160 | 40 | -
$\omega$ | 0.01 | 0.01 | -
ToF w/o Penalty (hr) | 33.24 | 33.34 | -
ToF w/ Penalty (hr) | 21.53 | 17.91 | -
TOF GPOPS (hr) | - | - | 17.38
ToF Reduction (%) | 35.23 | 46.29 | 47.87

###  6.2 Rendezvous

The same process was used to compare the trajectories obtained for the
rendezvous case. In addition to the boundary conditions given in Table 6, the
same design parameters from the orbit raising test case were used to produce
the results for this test case shown in Table 7. It is worth noting that
identical design parameters do not necessarily result in a good spacecraft
trajectory for each maneuver, however, in this case they did. In general,
different scenarios and boundary conditions will require unique FFS design
parameters.

Table 6: Boundary Conditions: Rendezvous Parameter | Value
---|---
$r_{i}$ | 6570 km
$\theta$ | 0 deg.
$\dot{r}_{i}$ | 0 km/s
$r_{f}$ | 42160 km
$\dot{r}_{f}$ | 0 km/s
Table 7: Design Parameters: Rendezvous Case 4 Parameter | Value
---|---
$n_{r}$ | 3
$n_{\theta}$ | 4
$DP$ | 160
$T_{a,max}$ |  $\pm$ 0.0102 DU/TU2
$\omega$ | 0.01

Two FFS trajectories, with and without penalties on the ToF, were generated
and are shown in Figure 11. The thrust profiles for these two trajectories are
shown in Figure 11 and it is clear that the penalized trajectory completes the
rendezvous in much shorter time. A larger thrust acceleration must be applied
to rendezvous in a shorter time period. Figure 12 shows the FFS trajectory
compares well to the trajectory obtained via numerically integrating the
thrust profile from the FFS approximation.

Figure 10: FFS Trajectories: Rendezvous

Figure 11: Thrust Acceleration: Rendezvous

Figure 12: Trajectory Comparison: Rendezvous

####  6.2.1 Trajectory Feasibility

The optimal trajectory and thrust acceleration profile that achieves the
minimum time rendezvous are shown in Figures 14 and 14. Contrary to the
previous case, the optimal thrust acceleration for the rendezvous case has a
bang-off-bang structure. The solution stays near zero until 15 TU,
oscillations occur around 20 TU, and the thrust profile touches the lower
bound for an instant around 25 TU. This oscillation may be attributed to the
fact that the final polar angle is updated every time a solver inside GPOPS is
called. As a result, an optimal time of flight of 22.08 hours is obtained
which corresponds to a 45.94 % reduction compared to the non-weighted FFS
trajectory time of flight.

Figure 13: GPOPS Trajectory: Rendezvous

Figure 14: GPOPS Thrust Acceleration: Rendezvous

All results for the rendezvous cases presented above are shown in table 8. The
GPOPS optimal ToF is compared to the non-penalized ToF result that is
generated in case 4, $ToF=40.84$ hours.

Table 8: Time of Flight Results: Rendezvous Parameter | Case 4 | GPOPS
---|---|---
$n_{r}$ | 3 | -
$n_{\theta}$ | 4 | -
$DP$ | 160 | -
$\omega$ | 0.01 | -
ToF w/o Penalty (hr) | 40.84 | -
ToF w/ Penalty (hr) | 25.19 | -
ToF GPOPS (hr) | - | 22.08
ToF Reduction (%) | 38.32 | 45.94

##  7 Conclusion

This paper built on the FFS SB method developed to approximate low-thrust
spacecraft trajectories by adding time of flight to the objective function[2].
A weight is applied to the terms in the objective function to balance the
reduction in the time of flight with a trajectory that satisfies the equations
of motion. The time of flight penalty was tested on orbit raising and
rendezvous maneuvers. In both cases the new objective function generated
trajectories which decreased the ToF. The reduced time of flight trajectories
were validated against numerically integrated results using the non-linear
equations of motion. Additionally, the results from FFS were compared to the
results from GPOPS to see how well the FFS trajectories and thrust
accelerations matched the optimal solution. The results showed that the FFS
obtained solutions that were close to the optimal solutions generated by GPOPS
for the test cases investigated.

## References

  * [1] John T. Betts.  Survey of numerical methods for trajectory optimization.  Journal of Guidance, Control, and Dynamics, 21(2):193–207, 1998\.
  * [2] Ehsan Taheri and Ossama Abdelkhalik.  Shape based approximation of constrained low-thrust space trajectories using fourier series.  Journal of Spacecraft and Rockets, 49(3):535–546, 2012.
  * [3] Ehsan Taheri and Ossama Abdelkhalik.  Fast initial trajectory design for low-thrust restricted-three-body problems.  Journal of Guidance, Control, and Dynamics, 38(11):2146–2160, 2015\.
  * [4] Ehsan Taheri and Ossama Abdelkhalik.  Initial three-dimensional low-thrust trajectory design.  Advances in Space Research, 57(3):889–903, 2016.
  * [5] Bradley J. Wall and Bruce A. Conway.  Shape-based approach to low-thrust rendezvous trajectory design.  Journal of Guidance, Control, and Dynamics, 32(1):95–101, 2009\.
  * [6] D. J. BENNEY.  Escape from a circular orbit using tangential thrust.  Journal of Jet Propulsion, 28(3):167–169, 1958.
  * [7] Frederick W. Boltz.  Orbital motion under continuous tangential thrust.  Journal of Guidance, Control, and Dynamics, 15(6):1503–1507, 1992\.
  * [8] Anastassios E. Petropoulos and James M. Longuski.  Shape-based algorithm for the automated design of low-thrust, gravity assist trajectories.  Journal of Spacecraft and Rockets, 41(5):787–796, 2004.
  * [9] A.E. Petropoulos.  A shape-based approach to automated, low-thrust, gravity-assist trajectory design.  Ph.D. Thesis, 2001.
  * [10] Michael A Patterson and Anil V Rao.  Gpops-ii: A matlab software for solving multiple-phase optimal control problems using hp-adaptive gaussian quadrature collocation methods and sparse nonlinear programming.  ACM Transactions on Mathematical Software (TOMS), 41(1):1–37, 2014\.


# A Third Information-Theoretic Approach
to Finite de Finetti Theorems

Mario Berta Department of Computing, Imperial College, London SW7 2AZ, U.K.,
and Institute for Quantum Information, RWTH Aachen University, Aachen,
Germany. Email<EMAIL_ADDRESS>M.B. acknowledges funding by the
European Research Council (ERC Grant Agreement No. 948139).  Lampros Gavalakis
LAMA - Laboratoire d’Analyse et Mathématiques Appliquées, Gustave Eiffel
University, Marne-la-Vallée Cedex 2, France. Email: lampros.gavalakis@univ-
eiffel.fr. L.G. was supported in part by the European Union’s Horizon 2020
research and innovation program under the Marie Sklodowska-Curie grant
agreement No 101034255 .  Ioannis Kontoyiannis Statistical Laboratory, DPMMS,
University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road,
Cambridge CB3 0WB, U.K. Email<EMAIL_ADDRESS>I.K. was supported in
part by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under
the “First Call for H.F.R.I. Research Projects to support Faculty members and
Researchers and the procurement of high-cost research equipment grant,”
project number 1034.

###### Abstract

A new finite form of de Finetti’s representation theorem is established using
elementary information-theoretic tools. The distribution of the first $k$
random variables in an exchangeable vector of $n\geq k$ random variables is
close to a mixture of product distributions. Closeness is measured in terms of
the relative entropy and an explicit bound is provided. This bound is tighter
than those obtained via earlier information-theoretic proofs, and its utility
extends to random variables taking values in general spaces. The core argument
employed has its origins in the quantum information-theoretic literature.

Keywords —  Exchangeability, de Finetti theorem, entropy, mixture, relative
entropy, mutual information

##  1 Introduction: de Finetti’s theorem

A vector of $n$ random variables $(X_{1},X_{2},\ldots,X_{n})$ is exchangeable
if it has the same distribution as $(X_{\pi(1)},X_{\pi(2)},\ldots,X_{\pi(n)})$
for every permutation $\pi$ on $\\{1,2,\ldots,n\\}$. Similarly, we say that a
process $\\{X_{k}\;;\;k\geq 1\\}$ is exchangeable if
$(X_{1},X_{2},\ldots,X_{n})$ is exchangeable for all $n$.

De Finetti’s celebrated representation theorem [12, 13] characterises all
binary exchangeable processes.

######  Theorem 1.1 (de Finetti’s representation theorem)

A binary process $\\{X_{k}\;;\;k\geq 1\\}$ is exchangeable if and only if its
distribution can be uniquely expressed as a mixture of independent and
identically distributed binary processes.

The importance of de Finetti’s theorem (and its numerous extensions) for the
foundations of Bayesian statistics is illustrated, e.g., by the observation
that it justifies the subjective approach of treating any infinite binary
exchangeable sequence as a sequence of independent coin tosses conditional on
the value of the probability of Heads, which is itself distributed according
to some prior distribution – namely, the mixing measure in Theorem 1.1; see
[16, 3] and the references therein for extensive discussions.

In terms of its practical applicability in statistics and elsewhere, it is
natural to ask whether the same representation applies to finite-length
exchangeable random vectors. As it turns out, the condition that the entire
process be exchangeable is necessary. For example, considering a pair of
binary random variables $(X_{1},X_{2})$ with
$\Pr(X_{1}=1,X_{2}=0)=\Pr(X_{1}=0,X_{2}=1)=1/2$, shows that $(X_{1},X_{2})$
are exchangeable but their joint distribution cannot be expressed as a mixture
of independent and identically distributed (i.i.d.) binary pairs [16].

Nevertheless, the distribution of the first $k$ random variables in an
exchangeable vector of length $n$ admits an approximate de Finetti-style
representation if $k$ is small compared to $n$. Quantitative versions of this
statement have been established by Diaconis [16] for binary random variables
and Diaconis and Freedman [18] for random variables with values in abstract
spaces, with sharp error bounds for the total variation distance. We refer to
such results as ‘finite de Finetti’ theorems. Diaconis’ proof in [16] is based
on a geometric interpretation of the set of exchangeable measures as a convex
subset of the probability simplex. The Diaconis and Freedman proof in [18] is
based on the approximation of the hypergeometric by the binomial distribution.
The general version of Theorem 1.1 due to Hewitt and Savage [22] for random
variables with values in compact Hausdorff spaces was also recovered in [18].
Interestingly, one needs additional assumptions in order to deduce the
infinite version from the – seemingly stronger – finite one.

###  1.1 Information-theoretic approaches

Since the early 1950s, there has been a long line of works establishing core
probabilistic results via information-theoretic ideas and techniques; parts of
this rich history are outlined in [11, 2, 23, 21]. In this spirit, two
different information-theoretic proofs of finite de Finetti theorems were
developed recently, one for binary alphabets [19] and one for finite alphabets
[21]. These results are expressed in terms of the relative entropy, which is a
stronger notion of ‘distance’ than total variation, but the resulting total
variation bounds are generally suboptimal. We briefly recall these results to
highlight their differences with the present development; a short review of
these two approaches may be found in [20].

For any $k$, write $X_{1}^{k}$ for a vector of random variables
$(X_{1},X_{2},\ldots,X_{k})$ and similarly write $x_{1}^{k}$ for a specific
realisation $(x_{1},x_{2},\ldots,x_{k})$ of $X_{1}^{k}$. When the $X_{i}$ all
take values in a finite alphabet, we denote their joint probability mass
function by $P_{X_{1}^{k}}$. If $\mu$ is a probability measure on the simplex
${\cal P}$ of probability mass functions $Q$ on $A$, we define the mixture:

$M_{k,\mu}(x_{1}^{k}):=\int_{\cal P}Q^{k}(x_{1}^{k})\,d\mu(Q),\quad
x_{1}^{k}\in A^{k}.$

The main result of the first information-theoretic proof [19] is that, if
$X_{1}^{n}$ is a binary exchangeable random vector, then there is a
probability measure $\mu$ on the space of Bernoulli distributions such that,
for every $k\leq n$,

$D(P_{X_{1}^{k}}\|M_{k,\mu})\leq\frac{5k^{2}\log n}{n-k},$ (1)

where $D(P\|Q)=\sum_{x\in B}P(x)\log[P(x)/Q(x)]$ denotes the relative entropy
between two probability mass functions on the same finite alphabet $B$, and
$\log$ denotes the natural logarithm throughout. Writing
$\|P-Q\|:=\sup_{B}|P(B)-Q(B)|$ for the total variation distance between $P$
and $Q$, Pinsker’s inequality [10, 24] states that
$\|P-Q\|^{2}\leq\frac{1}{2}D(P\|Q)$. Thus, (1) yields,

$\|P_{X_{1}^{k}}-M_{k,\mu}\|\leq
k\Bigl{(}\frac{5\log{n}}{2(n-k)}\Bigr{)}^{\frac{1}{2}}.$ (2)

This bound is suboptimal in that, as shown in [18], the correct rate with
respect to the total variation distance in (2) is $O(k/n)$. The proof in [19]
is based on an estimate of the dependence between the random variables
$X_{1}^{k}$ conditional on the proportion of 1s it contains.

In the second information-theoretic approach to finite de Finetti theorems
[21], a different finite-$n$ bound was obtained for the relative entropy in
(1), leading to an estimate of the form:

$D(P_{X_{1}^{k}}\|M_{k,\mu})=O\left(\Big{(}\frac{k}{\sqrt{n}}\Big{)}^{1/2}\log{\frac{n}{k}}\right).$
(3)

Although (3) is generally weaker than (1), it holds for random vectors
$X_{1}^{n}$ with values in arbitrary finite alphabets. The proof of the
finite-$n$ bound leading to (3) in [21] was based on a connection between the
Gibbs conditioning principle from statistical mechanics, the information-
theoretic ‘method of types’, and de Finetti’s theorem.

###  1.2 A different information-theoretic approach

The main goal of this note is to present and explore yet another information-
theoretic approach to developing finite de Finetti theorems, based only on
elementary properties of the mutual information. To put it in context, we
first remark that de Finetti-style theorems have attracted increasing
attention in the literature on quantum information theory in recent years. For
example, starting from related ideas in [6], the third part of [8, Theorem 1]
may be viewed as a classical de Finetti theorem with additional constraints,
and a classical de Finetti result may also be obtained from the subsequent
work [4, Theorem 2.4]. Additional relevant results can be found in [7, 5, 25]
and the references therein.

One of our results, given as Corollary 2.4 in the following section, states
that, if $X_{1}^{n}$ is an exchangeable vector of random variables taking
values in finite alphabet $A$, then:

$D(P_{X_{1}^{k}}\|M_{k,\mu})\leq\frac{k(k-1)}{2(n-k-1)}\log|A|.$

This bound, which is stronger than those obtained in [19] and in [21], will be
derived as a consequence of our main result, Theorem 2.1.

The main idea of the proof of the general bound in Theorem 2.1 has its origins
in quantum information theory. It is based on an argument that was implicitly
used in the proof of the quantum information-theoretic result in [8, Theorem
1]. The fact that this argument leads to interesting classical de Finetti-
style results appears to have been noted by several authors, and it is even
mentioned as an exercise in the recent text [27].

Our main purpose here is to explore and adapt this argument in order to show
that it can be employed to obtain finite-$n$ bounds for de Finetti’s theorem
in terms of relative entropy, which are both stronger and more general than
those obtained via the earlier information-theoretic proofs [19, 21].

Compared to the Diaconis and Freedman [18] variational distance bound
$k(k-1)/(2n)$, our information-theoretic corollary is weaker in terms of the
dependence on $n$ and $k$, and additionally features the logarithm of the
alphabet size. However, our bound holds in the stronger relative entropy
distance and it is then an interesting question if the logarithm of the
alphabet size is needed for that distance measure. Moreover, following [4, 8],
we note that our bound also holds for de Finetti theorems with linear
constraints. In the same way that Diaconis and Freedman’s result [17] gives a
polynomial time approximation scheme for the minimization of polynomials of
fixed degree over the simplex [14, 15], the versions with linear constraints
from [4] – that feature the logarithm of the alphabet size – give rise to
quasi-polynomial time approximation schemes for the minimization of
polynomials of fixed degree over the simplex intersected with affine
constraints [4]. Additionally, for this quasi-polynomial time approximation
scheme there also exists a matching hardness-of-approximation result [1]. As
such, we conclude that at least under certain complexity-theoretic
assumptions, any finite de Finetti proof strategy that can incorporate linear
constraints will likely need to feature the logarithm of the alphabet size in
the approximation error.

##  2 Finite de Finetti theorems

###  2.1 A general bound under exchangeability

We begin with some preliminary definitions and assumptions.

Let $(A,{\cal A})$ be a measurable space. Write ${\cal M}_{1}(A)$ for the
space of probability measures on $(A,{\cal A})$. For $P,Q\in{\cal M}_{1}(A)$,
the relative entropy between $P$ and $Q$ is,

$D(P\|Q)=\begin{cases}\int{\log{\frac{dP}{dQ}\;dP}},&\text{if }P\ll Q,\\\
+\infty,&\text{otherwise},\end{cases}$

where $dP/dQ$ denotes the Radon-Nikodým derivative of $P$ with respect to $Q$.
For a pair of random variables $(X,Y)$ with values in the measurable spaces
$(A_{1},{\cal A}_{1})$ and $(A_{2},{\cal A}_{2})$, respectively, and with
joint law $P_{X,Y}\in{\cal M}_{1}(A_{1}\times A_{2})$, the mutual information
between $X$ and $Y$ is,

$I(X;Y)=D(P_{X,Y}\|P_{X}\times P_{Y}),$

where $P_{X}\in{\cal M}_{1}(A_{1})$ and $P_{Y}\in{\cal M}_{1}(A_{2})$ denote
the first and second marginal of $P_{XY}$, respectively. In order to consider
mixture measures on given measurable space $(A,{\cal A})$, we first need to
equip ${\cal M}_{1}(A)$ with its own $\sigma$-algebra. For every $F\in{\cal
A}$, define the map $\pi_{F}:{\cal M}_{1}(A)\to[0,1]$ by,

$\pi_{F}(P)=P(F),\quad P\in{\cal M}_{1}(A),$

and let ${\cal F}$ denote the smallest $\sigma$-algebra of subsets of ${\cal
M}_{1}(A)$ that makes the maps $\\{\pi_{F}\;;\;F\in{\cal A}\\}$ measurable.
Now consider a specific measure $\mu\in{\cal M}_{1}({\cal M}_{1}(A))$ and a
given $k\geq 1$. For any $F\in{\cal A}^{k}$ the map $Q\mapsto Q^{k}(F)$ is
${\cal F}$-measurable, so we can define the mixture $M_{k,\mu}\in{\cal
M}_{1}(A)$ as,

$M_{k,\mu}(F):=\int_{{\cal M}_{1}(A)}Q^{k}(F)\,d\mu(Q),\quad F\in{\cal
A}^{k}.$

Throughout, we assume that $(A,{\cal A})$ is standard a Borel space, namely,
that there exists a complete separable metric space whose Borel
$\sigma$-algebra is isomorphic to $\mathcal{A}$. This ensures that regular
conditional probabilities exist [26]. Although in many cases, especially in
potential applications of our results, this assumption may be unnecessary, it
ensures that, in the general case, we will not be committing any measure-
theoretic faux pas.

We can now state our main result.

######  Theorem 2.1 (Exchangeability and information)

Suppose $X_{1}^{n}$ is an exchangeable vector of random variables $X_{i}$ with
values in a standard Borel space $(A,{\cal A})$, and with joint law
$P_{X_{1}^{n}}$. For every $1\leq k\leq n-1$ there exists a probability
measure $\mu=\mu_{k,n}$ on ${\cal M}_{1}(A)$, such that:

$D\big{(}P_{X_{1}^{k}}\big{\|}M_{k,\mu}\big{)}\leq\frac{1}{n-k+1}\sum_{i=1}^{k}{I(X_{1}^{i-1};X_{k}^{n}).}$
(4)

For the proof we will need the following two lemmas. The first one is a simple
identity for the relative entropy, stated without proof. The second one
provides a lower bound on the mutual informations between $X_{1}^{i}$ and
$X_{k}^{n}$, for $k>i$, which will be a key step in the proof of Theorem 2.1.

###### Lemma 2.2

For any random vector $Z_{1}^{k}$ of length $k$, the relative entropy between
the joint law $P_{Z_{1}^{k}}$ of $Z_{1}^{k}$ and the product
$\prod_{i=1}^{k}P_{Z_{i}}$ of its marginals can be expressed as:

$D\Big{(}P_{Z_{1}^{k}}\Big{\|}\prod_{i=1}^{k}P_{Z_{i}}\Big{)}=\sum_{i=1}^{k}I(Z_{1}^{i-1};Z_{i}).$

Recall that conditional mutual information between $X$ and $Y$ given $Z$ is
defined as,

$I(X;Y|Z)=\int{D\Bigl{(}P_{X,Y|Z}(\cdot|z)\Big{\|}P_{X|Z}(\cdot|z)\times
P_{Y|Z}(\cdot|z)\Bigr{)}\,dP_{Z}(z)},$

where $P_{X,Y|Z}(\cdot|z),P_{X|Z}(\cdot|z)$, $P_{Y|Z}(\cdot|z)$, denote the
conditional laws of $(X,Y)$, $X$, $Y$, given $Z$, respectively, whenever these
exist.

###### Lemma 2.3

Suppose $X_{1}^{n}$ are as in Theorem 2.1. Then, for each $1\leq i\leq k\leq
n-1$:

$\sum_{m=k}^{n}{I(X_{1}^{i-1};X_{i}|X_{k+1}^{m})}\leq
I(X_{1}^{i-1};X_{k}^{n}).$

Proof.  Let $1\leq i\leq k\leq n-1$ be fixed. For any $k\leq m\leq n$, by
exchangeability we have,

$I(X_{1}^{i-1};X_{i}|X_{k+1}^{m})=I(X_{1}^{i-1};X_{m}|X_{k}^{m-1}),$

with the obvious convention that $X_{k+1}^{k}$ is just a constant. Summing
over all $k\leq m\leq n$ and using the chain rule for mutual information
yields,

$\sum_{m=k}^{n}{I(X_{1}^{i-1};X_{i}|X_{k+1}^{m})}=\sum_{m=k}^{n}{I(X_{1}^{i-1};X_{m}|X_{k}^{m-1})}=I(X_{1}^{i-1};X_{k}^{n}),$

as claimed.  $\Box$

Proof of Theorem 2.1.  Adding the identities of Lemma 2.3 over all $1\leq
i\leq k$ and dividing by $(n-k+1)$,

$\frac{1}{n-k+1}\sum_{m=k}^{n}{\sum_{i=1}^{k}{I(X_{1}^{i-1};X_{i}|X_{k+1}^{m})}}\leq\frac{1}{n-k+1}\sum_{i=1}^{k}{I(X_{1}^{i-1};X_{k}^{n}).}$

Since the average over $m$ satisfies this upper bound, there is an
$m^{*}\in\\{k,k+1,\ldots,n\\}$ such that,

$\sum_{i=1}^{k}{I(X_{1}^{i-1};X_{i}|X_{k+1}^{m^{*}})}\leq\frac{1}{n-k+1}\sum_{i=1}^{k}{I(X_{1}^{i-1};X_{k}^{n}).}$
(5)

Let $Q^{*}=Q^{*}_{k,n}\in{\cal M}_{1}(A^{m^{*}-k})$ denote the joint law of
$X_{k+1}^{m^{*}}$. Taking $\mu=\mu_{k,n}$ to be the law of the regular
conditional probability $P_{X_{1}|X_{k+1}^{m^{*}}}(\cdot|X_{k+1}^{m^{*}})$,
the mixture $M_{k,\mu}$ can be expressed:

$M_{k,\mu}(F)=\int_{A^{m^{*}-k}}P_{X_{1}|X_{k+1}^{m^{*}}}^{k}(F|x_{k+1}^{m^{*}})\,dQ^{*}(x_{k+1}^{m^{*}}),\quad
F\in{\cal A}^{k}.$

Using the joint convexity of relative entropy,

$\displaystyle D\big{(}P_{X_{1}^{k}}\big{\|}M_{k,\mu}\big{)}$
$\displaystyle=D\Big{(}\int_{A^{m^{*}-k}}P_{X_{1}^{k}|X_{k+1}^{m^{*}}}(\cdot|x_{k+1}^{m^{*}})\,dQ^{*}(x_{k+1}^{m^{*}})\Big{\|}\int_{A^{m^{*}-k}}P_{X_{1}|X_{k+1}^{m^{*}}}^{k}(\cdot|x_{k+1}^{m^{*}})\,dQ^{*}(x_{k+1}^{m^{*}})\Big{)}$
$\displaystyle\leq\int
D\Big{(}P_{X_{1}^{k}|X_{k+1}^{m^{*}}}(\cdot|x_{k+1}^{m^{*}})\Big{\|}P_{X_{1}|X_{k+1}^{m^{*}}}^{k}(\cdot|x_{k+1}^{m^{*}})\Big{)}\,dQ^{*}(x_{k+1}^{m^{*}})$
$\displaystyle=\int
D\Big{(}P_{X_{1}^{k}|X_{k+1}^{m^{*}}}(\cdot|x_{k+1}^{m^{*}})\Big{\|}\prod_{i=1}^{k}P_{X_{i}|X_{k+1}^{m^{*}}}(\cdot|x_{k+1}^{m^{*}})\Big{)}\,dQ^{*}(x_{k+1}^{m^{*}}),$

where the last step follows from the fact that, by exchangeability,
$P_{X_{i}|X_{k+1}^{m^{*}}}$ is the same for all $i\leq k$. Finally, applying
Lemma 2.2 to the integrand, yields,

$D\big{(}P_{X_{1}^{k}}\big{\|}M_{k,\mu}\big{)}\leq\sum_{i=1}^{k}{I(X_{1}^{i-1};X_{i}|X_{k+1}^{m^{*}})},$

and the claimed bound follows from this combined with (5).  $\Box$

###  2.2 Explicit de Finetti-style theorems

For random variables $X_{1}^{k}$ on a discrete (finite or countably infinite)
alphabet $A$, the mutual information always satisfies [9],

$I(X_{1}^{k};Y)=H(X_{1}^{k})-H(X_{1}^{k}|Y)\leq H(X_{1}^{k})\leq k\log|A|,$

where $H(X)=-\sum_{x\in B}P(x)\log P(x)$ denotes the entropy of a random
variable with probability mass function $P$ on a discrete alphabet $B$.
Therefore, Theorem 2.1 immediately yields:

######  Corollary 2.4 (Finite de Finetti theorem for discrete random
variables)

Suppose $X_{1}^{n}$ is an exchangeable vector of random variables $X_{i}$
taking values in a discrete alphabet $A$. For every $1\leq k\leq n-1$ there
exists a probability measure $\mu=\mu_{k,n}$ on ${\cal M}_{1}(A)$, such that:

$D\big{(}P_{X_{1}^{k}}\big{\|}M_{k,\mu}\big{)}\leq\frac{k(k-1)}{2(n-k+1)}H(X_{1})\leq\frac{k(k-1)}{2(n-k+1)}\log|A|.$

Compared with the earlier information-theoretic results (1) and (3), the bound
in Corollary 2.4 is both more general and stronger. Moreover, it can be used
to recover the classical infinite version of de Finetti’s theorem for compact
spaces, under some conditions.

######  Corollary 2.5 (Classical de Finetti theorem for compact spaces)

Let $G$ be a compact metrisable space equipped with its Baire $\sigma$-algebra
${\cal G}$. Suppose the process $\\{X_{k}\;;\;k\geq 1\\}$ is exchangeable and
the $X_{k}$ take values in $G$ and are ${\cal G}$-measurable. If for every $k$
we have $I(X_{1}^{k-1};X_{k}^{n})=o(n)$ as $n\to\infty$, then there is a
probability measure $\mu$ on ${\cal M}_{1}(G)$ such that:

$P_{X_{1}^{k}}=M_{k,\mu},\quad\mbox{for every}\;k\geq 1.$

Note that Corollary 2.5 applies to all finite-alphabet exchangeable processes:
Any finite set $G=\\{a_{0},a_{1},\ldots,a_{m-1}\\}$ is compact with respect to
the metric $d(a_{i},a_{j})=|i-j|$ (mod $m$). And since
$I(X_{1}^{k-1};X_{k}^{n})\leq H(X_{1}^{k-1})\leq(k-1)\log m$, the condition
$I(X_{1}^{k-1};X_{k}^{n})=o(n)$ is always satisfied.

Proof.  Choose and fix $k\geq 1$. Since $G$ is compact and metrisable, it is
also complete and separable and thus standard Borel. And since $\\{X_{k}\\}$
is exchangeable, we can apply Theorem 2.1 in combination with Pinsker’s
inequality to obtain that, for any $n\geq k$,

$\|P_{X_{1}^{k}}-M_{k,\mu_{k,n}}\|\leq\left[\frac{1}{2(n-k+1)}\sum_{i=1}^{k}I(X_{1}^{i-1};X_{k}^{n})\right]^{1/2}.$
(6)

The condition $I(X_{1}^{k-1};X_{k}^{n})=o(n)$ implies that the right-hand side
vanishes as $n\to\infty$ and we are almost done, except that the mixing
measure $\mu_{k,n}$ depends on both $n$ and $k$.

Since $G$ is compact and metrisable, it is also Hausdorff. Following [18], we
note that ${\cal M}_{1}({\cal M}_{1}(A))$ is compact in the weak* topology, so
that there is a $\mu_{k}\in{\cal M}_{1}({\cal M}_{1}(A))$ and a subsequence
$\\{n_{j}\\}$ increasing to infinity such that $\mu_{k,n_{j}}\to\mu_{k}$ as
$j\to\infty$. Then, since the map $\mu\mapsto M_{k,\mu}$ is weak*-continuous,
we have that $M_{k,\mu_{k,n_{j}}}\to M_{k,\mu_{k}}$ in the weak* topology as
$j\to\infty$, and, in view of (6), we also have that $M_{k,\mu_{k,n_{j}}}\to
P_{X_{1}^{k}}$ in total variation. Therefore, we must have that,

$P_{X_{1}^{k}}=M_{k,\mu_{k}}.$

Since $k$ was arbitrary, by the fact that the marginals of the process
$\\{X_{k}\\}$ are necessarily consistent we must also have that
$P_{X_{1}^{k}}=M_{k,\mu_{n}}$, for each $k$ and all $n\geq k$. Using
weak*-compactness again, we can find a subsequence $\\{\mu_{k_{\ell}}\\}$ that
converges to some $\mu\in{\cal M}_{1}({\cal M}_{1}(A))$ in the weak* topology
as $\ell\to\infty$. By the consistency of the marginals and using again the
weak*-continuity of $M_{k,\mu}$ in $\mu$, we then have that for each $k$, as
$\ell\to\infty$, $P_{X_{1}^{k}}=M_{k,\mu_{k_{\ell}}}\to M_{k,\mu}$ in the
weak* topology, completing the proof.  $\Box$

###  2.3 Examples

Except for finite-valued processes (noted after Corollary 2.5) the condition
$I(X_{1}^{k-1};X_{k}^{n})=o(n)$ appears rather technical and is generally not
easy to verify. The following example describes a class of real-valued
exchangeable processes that satisfy the condition
$I(X_{1}^{k-1};X_{k}^{n})=o(n)$. Although these are not covered by Corollary
2.5 or the classical infinite de Finetti theorems (since ${\mathbb{R}}$ is not
compact), our Theorem 2.1 does provide useful bounds.

Example 1.  Consider a finite collection of densities
$\\{f_{1},f_{2},\ldots,f_{m}\\}$ on ${\mathbb{R}}$, such that the differential
entropy $h(f_{i})=-\int f_{i}\log f_{i}$ exists and is finite for each $i$.
Let $\theta$ be an arbitrary random variable with values in
$\\{1,2,\ldots,m\\}$. Define an exchangeable real-valued process
$\\{X_{k}\;;\;k\geq 1\\}$ as follows. Conditional on $\theta=i$, let
$\\{X_{k}\\}$ be i.i.d. with each $X_{k}\sim f_{i}$. Then, since conditioning
reduces the differential entropy [9],

$I(X_{1}^{k-1};X_{k}^{n})=h(X_{1}^{k-1})-h(X_{1}^{k-1}|X_{k}^{n})\leq
h(X_{1}^{k-1})-h(X_{1}^{k-1}|X_{k}^{n},\theta),$

and since $X_{1}^{k-1}$ is conditionally independent of $X_{k}^{n}$ given
$\theta$,

$I(X_{1}^{k-1};X_{k}^{n})=h(X_{1}^{k-1})-h(X_{1}^{k-1}|\theta)\leq(k-1)\bigl{[}h(X_{1})-\min_{1\leq
i\leq m}{h(f_{i})}\bigr{]}.$

Therefore, in this case Theorem 2.1 provides a useful bound of $O(k/n)$ for
the relative entropy.

On the other hand, the following is an example of an exchangeable process
$\\{X_{n}\\}$ on a compact space, for which the condition
$I(X_{1}^{k-1};X_{k}^{n})=o(n)$ fails.

Example 2.  Define a process $\\{X_{k}\;;\;k\geq 1\\}$ as follows. Let
$\theta\sim U[0,1]$ and, conditional on $\theta=x$, let the $X_{k}$ be i.i.d.
with each $X_{k}$ taking the values $x/2$ and $x$ with probability 1/2 each.
Then $\\{X_{k}\\}$ is a mixture of i.i.d. processes and hence exchangeable.
But it is easy to see that the joint law $P_{X_{1}^{n}}$ is never absolutely
continuous with respect to $P_{X_{1}^{k-1}}\times P_{X_{k}^{n}}$. Therefore,
not only do we not have $I(X_{1}^{k-1};X_{k}^{n})=o(n)$, but
$I(X_{1}^{k-1};X_{k}^{n})$ is equal to $+\infty$ and the result of Theorem 2.1
is trivial. One possible explanation for this is that, in some cases, the
mixing measure $\mu$ obtained in the proof of Theorem 2.1 is not the “right”
one.

## References

  * [1] S. Aaronson, R. Impagliazzo, and D. Moshkovitz.  AM with multiple Merlins.  In 2014 IEEE 29th Conference on Computational Complexity (CCC), pages 44–55, Vancouver, BC, June 2014.
  * [2] A.R. Barron.  Information theory in probability, statistics, learning, and neural nets.  In Y. Freundand and R.E. Schapire, editors, Proceedings of the Tenth Annual Conference on Computational Learning Theory (COLT), Nashville, Tennessee, July 1997.  Available at www.stat.yale.edu/~arb4/publications_files/COLT97.pdf.
  * [3] M. Bayarri and J.O. Berger.  The interplay of Bayesian and frequentist analysis.  Statist. Sci., 19(1):58–80, February 2004.
  * [4] M. Berta, F. Borderi, O. Fawzi, and V.B. Scholz.  Semidefinite programming hierarchies for constrained bilinear optimization.  Math. Program., 194(1-2):781–829, July 2022.
  * [5] F. Borderi.  De Finetti methods in quantum information.  PhD thesis, Department of Computing, Imperial College, University of London, London, U.K., April 2022.
  * [6] F.G.S.L. Brandão, M. Christandl, and J. Yard.  Faithful squashed entanglement.  Comm. Math. Phys., 306:805–830, September 2011.
  * [7] F.G.S.L. Brandão and A.W. Harrow.  Product-state approximations to quantum states.  Comm. Math. Phys., 342(1):47–80, 2016.
  * [8] F.G.S.L. Brandão and A.W. Harrow.  Quantum de Finetti theorems under local measurements with applications.  Comm. Math. Phys., 353(2):469–506, 4 2017.
  * [9] T.M. Cover and J.A. Thomas.  Elements of information theory.  John Wiley & Sons, New York, NY, second edition, 2012.
  * [10] I. Csiszár.  Information-type measures of difference of probability distributions and indirect observations.  Studia Sci. Math. Hungar., 2:299–318, 1967.
  * [11] I Csiszár.  Information theoretic methods in probability and statistics.  In 1997 IEEE International Symposium on Information Theory (ISIT), Ulm, Germany, June 1997.
  * [12] B. De Finetti.  Sul significato soggettivo della probabilita.  Fundam. Math., 17(1):298–329, 1931.
  * [13] B. De Finetti.  La prévision: Ses lois logiques, ses sources subjectives.  Ann. Inst. Henri Poincaré, 7(1):1–68, 1937.
  * [14] E. de Klerk, M. Laurent, and P.A. Parrilo.  A PTAS for the minimization of polynomials of fixed degree over the simplex.  Theor. Comput. Sci., 361(2):210–225, 2006.
  * [15] E. de Klerk, M. Laurent, and Z. Sun.  An alternative proof of a PTAS for fixed-degree polynomial optimization over the simplex.  Math. Program., 151:433–457, 2015.
  * [16] P. Diaconis.  Finite forms of de Finetti’s theorem on exchangeability.  Synthese, 36(2):271–281, 1977.
  * [17] P. Diaconis and D.A. Freedman.  De Finetti’s theorem for Markov chains.  Ann. Probab., 8(1):115–130, 1980.
  * [18] P. Diaconis and D.A. Freedman.  Finite exchangeable sequences.  Ann. Probab., 8(4):745–764, 1980.
  * [19] L. Gavalakis and I. Kontoyiannis.  An information-theoretic proof of a finite de Finetti theorem.  Electron. Comm. Probab., 26:1–5, 2021.
  * [20] L. Gavalakis and I. Kontoyiannis.  Information-theoretic de Finetti-style theorems.  In 2022 IEEE Workshop on Information Theory (ITW), pages 71–76, Mumbai, India, 2022.
  * [21] L. Gavalakis and I. Kontoyiannis.  Information in probability: Another information-theoretic proof of a finite de Finetti theorem.  In J.-M. Morel and B. Teissier, editors, Mathematics Going Forward: Collected Mathematical Brushstrokes, volume LNM 2313 of Lecture Notes in Mathematics. Springer, May 2023.
  * [22] E. Hewitt and L.J. Savage.  Symmetric measures on Cartesian products.  Trans. Amer. Math. Soc., 80(2):470–501, 1955.
  * [23] O. Johnson.  Information theory and the central limit theorem.  Imperial College Press, London, U.K., 2004.
  * [24] S. Kullback.  A lower bound for discrimination information in terms of variation.  IEEE Trans. Inform. Theory, 13(1):126–127, January 1967.
  * [25] K. Li and G. Smith.  Quantum de Finetti theorem under fully-one-way adaptive measurements.  Phys. Rev. Lett., 114:160503, April 2015.
  * [26] K.R. Parthasarathy.  Probability measures on metric spaces.  Academic Press, New York, NY, 1967.
  * [27] Y. Polyanskiy and Y. Wu.  Information theory: From coding to learning.  Book preprint, 2023.  Available at people.lids.mit.edu/yp/homepage.