# A new family of binary sequences with a low correlation via elliptic curves

Lingfei Jin, Liming Ma, Chaoping Xing and Runtian Zhu L. Jin and R. Zhu are
with the Shanghai Key Laboratory of Intelligent Information Processing, School
of Computer Science, Fudan University, Shanghai 200433, China, and with the
State Key Laboratory of Cryptology, P. O. Box, 5159, Beijing 100878, China.
(emails<EMAIL_ADDRESS>23110240100@m.fudan.edu.cn).L. Ma is with the
School of Mathematical Sciences, University of Science and Technology of
China, Hefei 230026, China (e-mail: lmma20@ustc.edu.cn).C. Xing is with School
of Electronic Information and Electrical Engineering, Shanghai Jiao Tong
University, Shanghai 200240, China (email: xingcp@sjtu.edu.cn).

###### Abstract

In the realm of modern digital communication, cryptography, and signal
processing, binary sequences with a low correlation properties play a pivotal
role. In the literature, considerable efforts have been dedicated to
constructing good binary sequences of various lengths. As a consequence,
numerous constructions of good binary sequences have been put forward.
However, the majority of known constructions leverage the multiplicative
cyclic group structure of finite fields $\mathbb{F}_{p^{n}}$, where $p$ is a
prime and $n$ is a positive integer. Recently, the authors made use of the
cyclic group structure of all rational places of the rational function field
over the finite field $\mathbb{F}_{p^{n}}$, and firstly constructed good
binary sequences of length $p^{n}+1$ via cyclotomic function fields over
$\mathbb{F}_{p^{n}}$ for any prime $p$ [8, 10]. This approach has paved a new
way for constructing good binary sequences. Motivated by the above
constructions, we exploit the cyclic group structure on rational points of
elliptic curves to design a family of binary sequences of length $2^{n}+1+t$
with a low correlation for many given integers $|t|\leqslant 2^{(n+2)/2}$.
Specifically, for any positive integer $d$ with $\gcd(d,2^{n}+1+t)=1$, we
introduce a novel family of binary sequences of length $2^{n}+1+t$, size
$q^{d-1}-1$, correlation bounded by $(2d+1)\cdot 2^{(n+2)/2}+|t|$, and a large
linear complexity via elliptic curves.

Index Terms—Binary sequences, low correlation, elliptic curves, automorphisms,
rational places.

##  I Introduction

Sequences with a low correlation have found essential applications in various
domains, including code-division multiple access (CDMA), spread spectrum
systems, broadband satellite communications and cryptography [1, 3, 12, 34].
These applications demand sequences with favorable characteristics, such as a
low correlation (autocorrelation and cross-correlation) and a large family
size. Since the later 1960s, many families of binary sequences with good
properties have been constructed via various methods.

Among them, the well-known Gold sequence family is the oldest binary sequence
of length $2^{n}-1$, size $2^{n}+1$ having three-level out-of-phase auto- and
cross-correlation for odd integer $n$. The family of Gold sequences can be
obtained by XOR-ing two m-sequences [5] of the same length with each other
[2]. Later, Kasami designed a sequence family using m-sequences and their
decimations [11], known as the Kasami sequences. Besides, there are some known
families of sequences of length $2^{n}-1$ with good correlation properties,
such as bent function sequences [21], No sequences [19], and Trace Norm
sequences [13]. In 2011, Zhou and Tang generalized the modified Gold sequences
and obtained a family of binary sequences with length $2^{n}-1$ for $n=2m+1$,
size $2^{\ell n}+\cdots+2^{n}-1$, and correlation $2^{m+\ell}+1$ for each
$1\leqslant\ell\leqslant m$ [35]. In addition to the above mentioned sequences
of length $2^{n}-1$, there is some work devoted to the construction of good
binary sequences with length of forms such as $2(2^{n}-1)$ and $(2^{n}-1)^{2}$
[30, 29, 28, 6]. In [30], Udaya and Siddiqi obtained a family of $2^{n-1}$
binary sequences, length $2(2^{n}-1)$ satisfying the Welch bound for odd $n$.
This family was extended to a larger one with $2^{n}$ sequences and the same
correlation in [29]. In [28], the authors presented two optimal families of
sequences of length $2^{n}-1$ and $2(2^{n}-1)$ for odd integer $n$. In 2002,
Gong presented a family of binary sequences with size $2^{n}-1$, length
$(2^{n}-1)^{2}$, and correlation $3+2(2^{n}-1)$ [6]. These constructions
mainly made use of the cyclic structure of multiplicative group
$\mathbb{F}_{2^{n}}^{*}$. In order to design sequences of length $2^{n}+1$,
the cyclic structure of all rational points of projective lines over
$\mathbb{F}_{2^{n}}$ in [9] is employed to construct a family of binary
sequence of length $2^{n}+1$, size $2^{n}-1$ and correlation bounded by
$2^{(n+2)/2}$ via cyclotomic function fields [10]. Using the cyclic structure
of rational points of elliptic curves, Xing et. al. provided a family of
binary sequences with a low correlation and a large linear span in [33].

There are also some constructions of binary sequences which are based on
finite fields of odd characteristics [14, 20, 22, 24, 27, 7]. Paterson
proposed a family of pseudorandom binary sequences based on Hadamard
difference sets and maximum distance separable (MDS) codes for length $p^{2}$
where $p\equiv 3\ (\text{mod }4)$ is a prime [22]. In 2006, a family of binary
sequences of length $p$ which is an odd prime, family size $(p-1)/2$, and
correlation bounded by $5+2\sqrt{p}$ was given in [24], now known as Weil
sequences. The idea of constructing Weil sequences is to derive sequences from
single quadratic residue based on Legendre sequences using a shift-and-add
construction. In 2010, by combining the $p$-periodic Legendre sequence and the
$(q-1)$-periodic Sidelnikov sequence, Su et al. introduced a new sequence of
length $p(q-1)$, called the Lengdre-Sidelnikov sequence [27]. In 2021, Jin et
al. provided several new families of sequences with low correlations of length
$q-1$ for any prime power $q$ by using the cyclic multiplication group
$\mathbb{F}_{q}^{*}$ and sequences of length $p$ for any odd prime $p$ from
the cyclic additive group $\mathbb{F}_{p}$ in [7]. The idea of constructing
binary sequences is using multiplicative quadratic characters over finite
fields of odd characteristic. By exploiting the cyclic group structure of
rational places of the rational function field $\mathbb{F}_{p^{n}}(x)$, a
family of binary sequence of length $p^{n}+1$, size $p^{n}-2$, and correlation
bounded by $4+2p^{n/2}$ was constructed via automorphisms of cyclotomic
function fields with odd characteristics in [8].

While several constructions of binary sequences with favorable parameters have
been introduced, the available choices for sequence lengths remain limited.
For applications in different scenarios, adding or deleting sequence values
usually destroys the good correlation properties. Consequently, there is a
pressing need for constructing binary sequences that exhibit low correlation
while offering greater flexibility in parameter selection. However,
identifying such binary sequences with novel parameters remains an open and
challenging problem in the field.

TABLE I: PARAMETERS OF SEQUENCE FAMILIES Sequence | Length $N$ | Family Size | Bound of Correlation
---|---|---|---
Gold(odd) [2] |  $2^{n}-1$, $n$ odd | $N+2$ | $1+\sqrt{2}\sqrt{N+1}$
Gold(even) [2] | $2^{n}-1,n=4k+2$ | $N+2$ | $1+2\sqrt{N+1}$
Kasami(small) [11] | $2^{n}-1,n\text{ even}$ | $\sqrt{N+1}$ | $1+\sqrt{N+1}$
Kasami(large) [11] | $2^{n}-1,n=4k+2$ | $(N+2)\sqrt{N+1}$ | $1+2\sqrt{N+1}$
Bent [21] | $2^{n}-1,n=4k$ | $\sqrt{N+1}$ | $1+\sqrt{N+1}$
No [19] | $2^{n}-1,n\text{ even}$ | $\sqrt{N+1}$ | $1+\sqrt{N+1}$
Trace Norm [13] | $2^{n}-1,n\text{ even}$ | $\sqrt{N+1}$ | $1+\sqrt{N+1}$
Tang et al. [29] |  $2(2^{n}-1),n$ odd | $\frac{N}{2}+1$ | $2+\sqrt{N+2}$
Gong [6] |  $(2^{n}-1)^{2}$, $2^{n}-1$ prime | $\sqrt{N}$ | $3+2\sqrt{N}$
Paterson [22] |  $p^{2}$, $p$ prime, $p\equiv 3(\text{mod }4)$ | $N$ | $5+4\sqrt{N}$
Paterson [22] |  $p^{2}$, $p$ prime, $p\equiv 3(\text{mod }4)$ | $\sqrt{N}+1$ | $3+2\sqrt{N}$
Weil [24] |  $p$, odd prime | $(N-1)/2$ | $5+2\sqrt{N}$
Jin et al. [7] |  $p-1$, $p\geqslant 17$ odd prime power | $N+3$ | $6+2\sqrt{N+1}$
Jin et al. [7] |  $p-1$, $p\geqslant 11$ odd prime power | $N/2$ | $2+2\sqrt{N+1}$
Jin et al. [7] |  $p$, $p\geqslant 17$ odd prime | $N$ | $5+2\sqrt{N}$
Jin et al. [7] |  $p$, $p\geqslant 11$ odd prime | $(N-1)/2$ | $1+2\sqrt{N}$
HJMX [8] | $p^{n}+1$, $p$ odd prime | $N-3$ | $4+2\sqrt{N-1}$
JMX [10] | $2^{n}+1$ | $N-2$ | $2\sqrt{N-1}$
XKD [33] | $2^{n}$ | $(N^{2}-2N+2)/3$ | $1+14\sqrt{N}$
XKD [33] | $q+1\pm\sqrt{q},q=4^{n}$ | $(q-\sqrt{q})/2$ | $11\sqrt{q}$
XKD [33] | $q+1\pm\sqrt{2q},q=2\times 4^{n}$ | $(q-\sqrt{2q})/2$ | $(10+\sqrt{2})\sqrt{q}$
Our construction | ${2^{n}}$ | ${N^{2}-1}$ | ${1+14\sqrt{N}}$
Our construction | ${q+1\pm\sqrt{q},q=4^{n}}$ | ${q-1}$ | ${11\sqrt{q}}$
Our construction | ${q+1\pm\sqrt{2q},q=2\times 4^{n}}$ | ${q-1}$ | ${(10+\sqrt{2})\sqrt{q}}$

###  I-A Our main result and techniques

In this paper, we provide an explicit construction of binary sequences of
length $2^{n}+1+t$ with a low correlation via cyclic elliptic function fields
over the finite field $\mathbb{F}_{2^{n}}$. Let $d\geqslant 2$ be a positive
integer with $\gcd(d,2^{n}+1+t)=1$. We establish a family of binary sequences
of length $2^{n}+1+t$, size $q^{d-1}-1$ and correlation upper bounded by
$2(2d+1)\sqrt{q}+1$. For the purpose of comparison, we provide a summary of
various families of binary sequences with a low correlation in Table I.

In order to understand the main idea of this article better, we give a high-
level description of our techniques. We use language of curves and function
fields interchangeably. To produce a sequence of length $q+1+t$, we use the
additive cyclic group structures under the group law $\oplus$ of rational
places of cyclic elliptic curves over a finite field $\mathbb{F}_{q}$ [23].
The automorphism group of an elliptic curve $E$ over $\mathbb{F}_{q}$ is in
the form $T_{E}\rtimes\operatorname{{Aut}}(E,O)$, where
$\operatorname{{Aut}}(E,O)$ is the automorphisms of $E$ fixing the infinity
place $O$, and $T_{E}$ corresponds to translation group of $E$ [16]. In
particular, $T_{E}$ is a cyclic group which is generated by a rational place
$P$ in this case, and there exists a generator $\sigma=\tau_{P}\in T_{E}$ that
makes the set $\\{\sigma^{j}(P)\\}_{j=0}^{q+t}$ being exactly the set of
rational places of $E$. Our sequences are constructed via the absolute trace
of evaluations of rational functions on such cyclically ordered rational
places $\\{\sigma^{j}(P)\\}_{j=0}^{q+t}$. To obtain a lower correlation, we
need to choose rational functions $z_{1},z_{2},\cdots,z_{\ell}$ for some
positive integer $\ell$ carefully. The cyclic group structure of $T_{E}$
allows us to convert the process of determining the correlation of this family
of binary sequences to the estimation of number of rational places on the
Artin-Schreier extensions $E_{j}=E(y)$ where
$y^{2}+y=z_{j}+\sigma^{-u}(z_{j})$ with $1\leqslant j\leqslant\ell$,
$1\leqslant u\leqslant q+t$ and Artin-Schreier extensions $E_{j,k}=E(y)$ where
$y^{2}+y=z_{j}+\sigma^{-u}(z_{k})$ with $1\leqslant j\neq k\leqslant\ell$,
$0\leqslant u\leqslant q+t$. In order to obtain Artin-Schreier curves with
small genus, first we need to choose function $z_{i}$ such that $z_{i}$ and
$\sigma^{-t}(z_{j})$ have less poles. Hence, we choose $z_{i}$ from a Riemann-
Roch space associated to a divisor which is just a place $Q$ with degree
$d\geqslant 2$ of $E$. In order to obtain Artin-Schreier curves with less
rational places, $z_{j}+\sigma^{-t}(z_{k})$ cannot be constant, since constant
field extensions may produce more rational points. Thus, $z_{j}$ can be chosen
as nonzero elements from the complementary dual space of $\mathbb{F}_{q}$ in
the Riemann-Roch space $\mathcal{L}(Q)$.

###  I-B Organization of this paper

In Section II, we provide preliminaries on binary sequences, extension theory
of function fields, elliptic function fields and Artin-Schreier extensions. In
Section III, we present a theoretical construction of binary sequences with a
low correlation via elliptic function fields over finite fields with
characteristic two. In Section IV, we give an algorithm to generate such a
family of binary sequences and provide many numerical results with the help of
the Sage software.

##  II Preliminaries

In this section, we present some preliminaries on binary sequences, extension
theory of function fields, elliptic function fields and Artin-Schreier
extensions.

###  II-A Binary sequences

For a positive integer $N$, a binary sequence of length $N$ is a vector
$\mathbf{s}\in\mathbb{F}_{2}^{N}$, where $\mathbb{F}_{2}$ is the finite field
of two elements. Let $\mathcal{S}$ be a family of binary sequences in
$\mathbb{F}_{2}^{N}$. The number of binary sequences in $\mathcal{S}$ is
called the family size of $\mathcal{S}$. For a binary sequence
$\mathbf{s}=(s_{0},s_{1},\dots,s_{N-1})$, the autocorrelation
$A_{t}(\mathbf{s})$ of $\mathbf{s}$ at delay $t$ for $1\leqslant t\leqslant
N-1$ is defined as

$A_{t}(\mathbf{s}):=\sum_{j=0}^{N-1}(-1)^{s_{j}+s_{j+t}},$

where $j+t$ is identified with the least non-negative integer after taking
modulo $N$. For two distinct binary sequences $\mathbf{u}=(u_{0},u_{1}$,
$\dots,u_{N-1})$ and $\mathbf{v}=(v_{0},v_{1},\dots,v_{N-1})$, their cross-
correlation at delay $t$ with $0\leqslant t\leqslant N-1$ is defined as

$C_{t}(\mathbf{u},\mathbf{v}):=\sum_{j=0}^{N-1}(-1)^{u_{j}+v_{j+t}}.$

For the family of binary sequences $\mathcal{S}$, its correlation is defined
by

$\operatorname{Cor}(\mathcal{S}):=\max\left\\{\max_{\begin{subarray}{c}\mathbf{s}\in\mathcal{S}\\\
1\leqslant t\leqslant
N-1\end{subarray}}\\{|A_{t}(\mathbf{s})|\\},\max_{\begin{subarray}{c}\mathbf{u}\neq\mathbf{v}\in\mathcal{S}\\\
0\leqslant t\leqslant
N-1\end{subarray}}\\{|C_{t}(\mathbf{u},\mathbf{v})|\\}\right\\}.$

In addition, we can define the linear complexity of binary sequences as
follows. For any sequence $\mathbf{0}\neq\mathbf{s}\in\mathcal{S}$, the linear
complexity of $\mathbf{s}$, denoted by $\ell(\mathbf{s})$, is the smallest
positive integer $\ell$ such that there exists $\ell+1$ elements
$\lambda_{0},\lambda_{1},\dots,\lambda_{\ell}\in\mathbb{F}_{2}$ with
$\lambda_{0}=\lambda_{\ell}=1$ satisfying

$\sum_{i=0}^{\ell}\lambda_{i}s_{i+u}=0$

for all $0\leqslant u\leqslant N-1$. Define the linear complexity of the
family of sequences $\mathcal{S}$ by

$\operatorname{LC}(\mathcal{S}):=\min\\{\ell(\bf{s}):\bf{s}\in\mathcal{S}\\}.$

###  II-B Extension theory of function fields

Let $F/\mathbb{F}_{q}$ be an algebraic function field with genus $g$ over the
full constant field $\mathbb{F}_{q}$. Let $\mathbb{P}_{F}$ denote the set of
places of $F$. Any place with degree one of $F$ is called rational, and all
rational places of $F$ form a set denoted by $\mathbb{P}_{F}^{1}$. From the
Serre bound [26, Theorem 5.3.1], the number $N(F)$ of rational places of $F$
is bounded by

$|N(F)-q-1|\leqslant g\cdot\lfloor 2\sqrt{q}\rfloor.$

Here $\lfloor x\rfloor$ stands for the integer part of $x\in\mathbb{R}$.

Let $\nu_{P}$ be the normalized discrete valuation of $F$ with respect to the
place $P\in\mathbb{P}_{F}$. The principal divisor of a nonzero element $z\in
F$ is given by $(z):=\sum_{P\in\mathbb{P}_{F}}\nu_{P}(z)P$. For a divisor $G$
of $F/\mathbb{F}_{q}$, the Riemann-Roch space associated to $G$ is defined by

$\mathcal{L}(G):=\\{z\in F^{*}:\;(z)+G\geqslant 0\\}\cup\\{0\\}.$

If $\deg(G)\geqslant 2g-1$, then $\mathcal{L}(G)$ is a vector space over
$\mathbb{F}_{q}$ of dimension $\deg(G)-g+1$ from the Riemann-Roch theorem [26,
Theorem 1.5.17].

Let $\operatorname{{Aut}}(F/\mathbb{F}_{q})$ denote the automorphism group of
$F$ over $\mathbb{F}_{q}$, i.e.,

$\operatorname{{Aut}}(F/\mathbb{F}_{q})=\\{\sigma:F\rightarrow
F\;|\;\sigma\mbox{ is an }\mbox{automorphism of }F\text{ and
}\sigma|_{\mathbb{F}_{q}}=\mathrm{id}\\}.$

Let us consider the group action of automorphism group
$\operatorname{{Aut}}(F/\mathbb{F}_{q})$ on the set of places of $F$. For any
place $P$ of $F$, $\sigma(P)$ is a place of $F$ for any automorphism
$\sigma\in\operatorname{{Aut}}(F/\mathbb{F}_{q})$. From [18, Lemma 1], we have
the following results.

######  Lemma II.1.

For any automorphism $\sigma\in\operatorname{{Aut}}(F/\mathbb{F}_{q})$,
$P\in\mathbb{P}_{F}$ and $f\in F$, we have

  * (1)

$\deg(\sigma(P))=\deg(P)$;

  * (2)

$\nu_{\sigma(P)}(\sigma(f))=\nu_{P}(f)$;

  * (3)

$\sigma(f)(\sigma(P))=f(P)$, provided that $\nu_{P}(f)\geqslant 0$.

Let $E/\mathbb{F}_{q}$ be a finite extension of $F/\mathbb{F}_{q}$. The
Hurwitz genus formula [26, Theorem 3.4.13] yields

$2{g(E)}-2=[E:F]\cdot(2g(F)-2)+\deg\operatorname{Diff}(E/F),$

where $\operatorname{Diff}(E/F)$ stands for the different of $E/F$. For any
$P\in\mathbb{P}_{F}$ and $Q\in\mathbb{P}_{F}$ with $Q|P$, let $d(Q|P),e(Q|P)$
and $f(Q|P)$ be the different exponent, ramification index and relative degree
of $Q|P$, respectively. The different of $E/F$ is defined as
$\text{Diff}(E/F)=\sum_{Q\in\mathbb{P}_{F}}d(Q|P)Q$ and different exponents
can be determined by the Dedekind different theorem and Hilbert’s different
theorem [26, Theorem 3.5.1 and Theorem 3.8.7]. In particular, the
decomposition group of $Q$ over $P$ is defined as $G_{Z}(Q|P)=\\{\sigma\in{\rm
Gal}(E/F)\colon\sigma(Q)=Q\\}$ and its order is $e(Q|P)f(Q|P)$.

###  II-C Elliptic function fields

An elliptic function field $E$ over $\mathbb{F}_{q}$ is an algebraic function
field over $\mathbb{F}_{q}$ of genus one and a rational place $O$ which is
called the infinite place of $E$. The divisor group Div$(E)$ of
$E/\mathbb{F}_{q}$ is defined as the free abelian group generated by
$\mathbb{P}_{E}$. The set of divisors of degree zero forms a subgroup of
Div$(E)$, denoted by $\text{Div}^{0}(E)$. Two divisors $A$ and $B$ in
$\text{Div}(E)$ are called equivalent if there exist $z\in E^{*}$ such that
$A=B+(z)$, which is denote by $A\sim B$. The set of principal divisors
$\text{Princ}(E)=\\{(x)=\sum_{P\in\mathbb{P}_{E}}\nu_{P}(x)P:x\in E^{*}\\}$ is
called the group of principal divisors of $E/\mathbb{F}_{q}$. The factor group
$\text{Cl}^{0}(E)=\text{Div}^{0}(E)/\text{Princ}(E)$ is called the group of
divisor classes of degree zero of $E/\mathbb{F}_{q}$. There is a bijective map

$\Phi\colon\left\\{\begin{array}[]{rcl}\mathbb{P}_{E}^{1}&\longrightarrow&\operatorname{Cl}^{0}(E)\\\
P&\longmapsto&[P-O]\end{array}\right.$

from the set of rational places of $E$ and the group of divisor classes of
degree zero, where $[\cdot]$ denotes the equivalence class of $\cdot$ in the
quotient group $\operatorname{Cl}^{0}(E)$. Hence, the set of rational places
$\mathbb{P}^{1}_{E}$ can be endowed with a group operation $\oplus$ as

$P\oplus Q=\Phi^{-1}(\Phi(P)+\Phi(Q)).$

It is easy to see that $P\oplus Q=R$ if and only if $P+Q\sim R+O$. Under this
operation, $O$ is the identity of the additive group $\mathbb{P}_{E}^{1}$. For
convenience, we can write $[2]P=P\oplus P$ and inductively define
$[m+1]P=[m]P\oplus P$ for any integer $m\in\mathbb{Z}$.

######  Lemma II.2.

Let $E/\mathbb{F}_{q}$ be an elliptic function field and let $P,Q$ be two
rational places of $E$. Then $P-Q=(z)$ for some element $z\in E^{*}$ if and
only if $P=Q$.

For elliptic curves, we can consider the isogenies between them, which are
non-constant smooth morphisms that preserve the infinity point. It is well
known that elliptic curves over $\mathbb{F}_{q}$ are isogenous if and only if
they have the same number of $\mathbb{F}_{q}$-rational points [32].

######  Lemma II.3.

The isogeny classes of elliptic curves over $\mathbb{F}_{q}$ for $q=p^{n}$ are
in one-to-one correspondence with the rational integers $t$ whose absolute
value not exceed $2\sqrt{q}$ and satisfying some one of the following
conditions:

  1. (i)

$\gcd(t,p)=1$;

  2. (ii)

$t=\pm 2\sqrt{q}$ if $n$ is even;

  3. (iii)

$t=\pm\sqrt{q}$ if $n$ is even and $p\not\equiv 1\pmod{3}$;

  4. (iv)

$t=\pm p^{(n+1)/2}$ if $n$ is odd and $p=2$ or $3$;

  5. (v)

$t=0$ if either (i) $n$ is odd or (ii) $n$ is even and $p\not\equiv
1\pmod{4}$.

Furthermore, an elliptic curve in the isogeny class corresponding to $t$ has
$q+1+t$ rational points.

The additive group structure of rational places $\mathbb{P}^{1}_{E}$ of any
elliptic function field $E/\mathbb{F}_{q}$ can be found from [23, Theorem 3].

######  Lemma II.4.

Let $\mathbb{F}_{q}$ be the finite field with $q=p^{n}$ elements. Let
$h=\prod_{\ell}\ell^{h_{\ell}}$ be a possible number of rational places of an
elliptic function field $E$ of $\mathbb{F}_{q}$. Then all the possible groups
$\mathbb{P}_{E}^{1}$ are the following

$\mathbb{Z}/p^{h_{p}}\mathbb{Z}\times\prod_{\ell\neq
p}\left(\mathbb{Z}/\ell^{a_{\ell}}\mathbb{Z}\times\mathbb{Z}/\ell^{h_{\ell}-a_{\ell}}\mathbb{Z}\right)$

with

  * (a)

in case (ii) of Lemma II.3: each $a_{\ell}$ is equal to $h_{\ell}/2$, i.e,
$\mathbb{P}_{E}^{1}\cong\mathbb{Z}/(\sqrt{q}\pm
1)\mathbb{Z}\times\mathbb{Z}/(\sqrt{q}\pm 1)\mathbb{Z}.$

  * (b)

in other cases of Lemma II.3: $a_{\ell}$ is an arbitrary integer satisfying
$0\leqslant a_{\ell}\leqslant\min\\{\nu_{\ell}(q-1),[h_{\ell}/2]\\}$. In cases
(iii) and (iv) of Lemma II.3: $\mathbb{P}_{E}^{1}\cong\mathbb{Z}/h\mathbb{Z}.$
In cases (v) of Lemma II.3: if $q\not\equiv-1(\text{mod }4)$, then
$\mathbb{P}_{E}^{1}\cong\mathbb{Z}/(q+1)\mathbb{Z}$; otherwise,
$\mathbb{P}_{E}^{1}\cong\mathbb{Z}/(q+1)\mathbb{Z}$ or
$\mathbb{P}_{E}^{1}\cong\mathbb{Z}/2\mathbb{Z}\times\mathbb{Z}/\frac{q+1}{2}\mathbb{Z}$.

If all rational places of an elliptic function field $E/\mathbb{F}_{q}$ form a
cyclic group under the addition $\oplus$, then $E/\mathbb{F}_{q}$ is called a
cyclic elliptic function field. In particular, we will employ this cyclic
structure of the rational places of elliptic function fields with even
characteristic to construct binary sequences as [33, Lemma 5.3].

######  Lemma II.5.

Let $q=2^{n}$. Let $t$ be an integer satisfying one of the following
conditions:

  * (1)

$t$ is an odd integer between $-2\sqrt{q}$ and $2\sqrt{q}$;

  * (2)

$t=0$;

  * (3)

$t=\pm\sqrt{q}$ if $n$ is even, and $t=\pm\sqrt{2q}$ if $n$ is odd.

Then there exists a cyclic elliptic function field $E$ over $\mathbb{F}_{q}$
has $q+1+t$ rational places.

###  II-D Artin-Schreier extensions

In this subsection, we summarize the main properties of Artin-Schreier
extensions from [17, Theorem 3.1.9].

######  Proposition II.6.

Let $q=2^{n}$ be a prime power. Let $F/\mathbb{F}_{q}$ be the function field
over $\mathbb{F}_{q}$ with genus $g$. Let $f$ be an element in $F$ such that
$f\neq z^{2}+z$ for all $z\in F$. Let

$F^{\prime}=F(y)\text{ with }y^{2}+y=f.$

For any place $P\in\mathbb{P}_{F}$, we define the integer $m_{P}$ by

$m_{P}=\begin{cases}\ell&\text{ if }\exists z\in F\text{ satisfying
}\nu_{P}(u-(z^{2}+z))=-\ell\\\ &\text{ with }\ell\in\mathbb{Z}_{+}\text{ and
}\ell\not\equiv 0(\text{mod }2),\\\ -1&\text{ if
}\nu_{P}(u-(z^{2}+z))\geqslant 0\text{ for some }z\in F.\end{cases}$

Then we have the following results.

  * (a)

$F^{\prime}/F$ is a cyclic Galois extension of degree two. The automorphisms
of $F^{\prime}/F$ are given by $\sigma(y)=y+u$ with $u=0,1.$

  * (b)

$P$ is unramified in $F^{\prime}/F$ if and only if $m_{P}=-1$.

  * (c)

$P$ is totally ramified in $F^{\prime}/F$ if and only if $m_{P}>0$. Denote by
$P^{\prime}$ the unique place of $F^{\prime}$ lying over $P$. Then the
different exponent $d(P^{\prime}|P)$ is given by

$d(P^{\prime}|P)=m_{P}+1.$

  * (d)

Assume that there is at least one place $Q$ satisfying $m_{Q}>0$, then the
full constant field of $F^{\prime}$ is $\mathbb{F}_{q}$ and the genus of
$F^{\prime}$ is

$g(F^{\prime})=2g-1+\frac{1}{2}\sum_{P\in\mathbb{P}_{F}}(m_{P}+1)\cdot\deg(P).$

##  III Binary sequences with a low correlation via elliptic curves

Motivated by the construction of families of binary sequences with a low
correlation via automorphism groups of cyclotomic function fields in [8, 10],
we provide a new construction of binary sequences via elliptic curves.
Compared with the result in [33], our construction achieves the same low
correlation and high linear complexity but offers a larger family size.

###  III-A A new construction of binary sequences via elliptic curves

Let $n$ be a positive integer and $q=2^{n}$. Let $\mathbb{F}_{q}$ be the
finite field with $q$ elements. In this subsection, we shall employ the cyclic
group structure of rational points and automorphism groups of elliptic curves
over $\mathbb{F}_{2^{n}}$ to construct a new family of binary sequences with a
low correlation and a large linear complexity.

Let $E/\mathbb{F}_{q}$ be an elliptic function field defined over
$\mathbb{F}_{q}$ with $q+1+t$ rational places. Let $\tau_{P}$ be the
automorphism of $E$ induced by $\tau_{P}(Q)=P\oplus Q$ for any rational place
$P$ and $Q$ of $E$. Let $T_{E}=\\{\tau_{P}:P\in\mathbb{P}^{1}_{E}\\}$ be the
translation group of $E$. Let $O$ be the infinity place of $O$ and
$\operatorname{{Aut}}(E,O)$ be the automorphism group of $E$ fixing the
infinity place $O$. From [16, Theorem 3.1], the automorphism group of
$E/\mathbb{F}_{q}$ is $\operatorname{{Aut}}(E/\mathbb{F}_{q})\cong
T_{E}\rtimes\operatorname{{Aut}}(E,O)$. Let $F$ be the fixed subfield of $E$
with respect to the translation group $T_{E}$, i.e.,

$F=E^{T_{E}}=\\{z\in E:\sigma(z)=z\text{ for any }\sigma\in T_{E}\\}.$

From Galois theory, $E/F$ is a Galois extension with Galois group ${\rm
Gal}(E/F)\cong T_{E}$. Since the number of rational places of $E$ is $q+1+t$,
there exists a rational place of $F$ which splits completely into $q+1+t$
rational places of $E$ from [26, Theorem 3.1.12 and Theorem 3.7.2].

In particular, we assume that $t$ satisfies one of the three conditions given
in Lemma II.5. Hence, there is a cyclic elliptic function field
$E/\mathbb{F}_{q}$ with $q+1+t$ rational places, i.e., there exists a rational
place $P$ such that the least positive integer $k$ satisfying $[k]P=O$ is
$q+1+t$. Hence, the translation group $T_{E}$ of $E$ is a cyclic group of
order $q+1+t$ with generator $\tau_{P}$. All rational places of $E$ inherits
the cyclic group structure of $T_{E}=\langle\tau_{P}\rangle$ as follows. Let
$P_{j}=\tau_{P}^{j}(O)=[j]P$ be the rational place of $E$ for each $0\leqslant
j\leqslant q+t$. By [26, Theorem 3.7.1], the set
$\\{P_{0},P_{1},\cdots,P_{q+t}\\}$ consists of all rational places of $E$.

######  Proposition III.1.

Let $E/\mathbb{F}_{q}$ be a cyclic elliptic function field defined over
$\mathbb{F}_{q}$ with $q+1+t$ rational places. Let
$T_{E}=\langle\tau_{P}\rangle$ be the translation group of $E$ and
$F=E^{T_{E}}$. Let $d$ be a positive integer and $Q$ be a place of $E$ with
degree $d$. If $\gcd(d,q+1+t)=1$, then $\tau_{P}^{j}(Q)$ are pairwise distinct
places with degree $d$ of $E$ for all $0\leqslant j\leqslant q+t$.

###### Proof.

From [33, Lemma 5.1] or [4], for any $j\in\mathbb{N}$, there exists a unique
automorphism $\tau_{P}^{j}=\tau_{[j]P}$ of
$\operatorname{{Aut}}(E/\mathbb{F}_{q})$ such that

$\tau_{[j]P}(Q)+O-Q-[d][j]P$

is a principal divisor of $E$ for any place $Q$ of $E$ with degree $d$. If
$\tau_{[j]P}(Q)=Q$, then we have $[dj]P-O$ is a principal divisor. By Lemma
II.2, we have $(q+1+t)|dj$ for some $j$. Since $\gcd(d,q+1+t)=1$, we have
$(q+1+t)|j$. Hence, the decomposition group of $Q|Q\cap F$ is the identity
automorphism of $E$. Thus, $\tau_{P}^{j}(Q)$ are pairwise distinct places with
degree $d$ of $E$ for all $0\leqslant j\leqslant q+t$ by [26, Corollary
3.7.2]. ∎

######  Lemma III.2.

Let $E/\mathbb{F}_{q}$ be an elliptic function field defined over
$\mathbb{F}_{q}$ with $q+1+t$ rational places. Let $B_{d}$ be the number of
places with degree $d$ of $E$. Let $\mu(\cdot)$ denote the M$\ddot{o}$bius
function. Then one has

$B_{d}=\frac{1}{d}\cdot\sum_{r|d}\mu\left(\frac{d}{r}\right)\cdot\left(q^{r}+1-\sum_{i=0}^{\lfloor
r/2\rfloor}(-1)^{r-i}\frac{(r-i-1)!\cdot r}{(r-2i)!\cdot
i!}t^{r-2i}q^{i}\right).$

In particular, we have

$B_{2}=\frac{q^{2}+q-t^{2}-t}{2}\text{ and
}B_{3}=\frac{q^{3}-q+t^{3}-3qt-t}{3}.$

###### Proof.

By [26, Theorem 5.1.15], the $L$-polynomial of $E$ is given by

$L_{E}(X)=1+[N(E)-(q+1)]X+qX^{2}=1+tX+qX^{2}.$

Assume that $L_{E}(X)$ factorizes in $\mathbb{C}[X]$ as

$L_{E}(X)=(1-\alpha_{1}X)(1-\alpha_{2}X).$

Thus, we obtain $\alpha_{1}+\alpha_{2}=-t$ and $\alpha_{1}\alpha_{2}=q$. Let
$S_{r}$ be the sum $\alpha_{1}^{r}+\alpha_{2}^{r}$. By [26, Proposition
5.2.9], for all $d\geqslant 2$, we have

$B_{d}=\frac{1}{d}\cdot\sum_{r|d}\mu\left(\frac{d}{r}\right)(q^{r}+1-S_{r}).$

From Waring’s Formula [15, Theorem 1.76], we have

$S_{r}=\sum_{i=0}^{\lfloor r/2\rfloor}(-1)^{r-i}\frac{(r-i-1)!\cdot
r}{(r-2i)!\cdot i!}t^{r-2i}q^{i}.$

This completes the proof. ∎

Let $E$ be a cyclic elliptic function field with $q+1+t$ rational places. Let
$Q$ be a place of $E$ with $\deg(Q)=d\geqslant 2$ such that $\gcd(d,q+1+t)=1$.
The existence of such a place $Q$ can be guaranteed by Lemma III.2. Since
$\deg(Q)=d\geqslant 2g(E)-1=1$, the dimension of the Riemann-Roch space
$\mathcal{L}(Q)$ over $\mathbb{F}_{q}$ is
$\dim\mathcal{L}(Q)=\deg(Q)+1-g(E)=\deg(Q)=d$. It is clear that
$\mathbb{F}_{q}=\mathcal{L}(0)\subseteq\mathcal{L}(Q)$. Hence, there exists a
vector space $V$ over $\mathbb{F}_{q}$ of dimension $d-1$ such that
$\mathcal{L}(Q)=\mathbb{F}_{q}\oplus V$, which is the complementary dual space
of $\mathbb{F}_{q}$ in $\mathcal{L}(Q)$.

Let $\operatorname{Tr}$ be the absolute trace map from $\mathbb{F}_{q}$ to
$\mathbb{F}_{2}$ which is defined by
$\text{Tr}(\alpha)=\sum_{i=0}^{n-1}\alpha^{2^{i}}=\alpha+\alpha^{2}+\cdots+\alpha^{2^{n-1}}$
for any $\alpha\in\mathbb{F}_{q}$. Since $\dim_{\mathbb{F}_{q}}(V)=d-1$, the
set $V\setminus\\{0\\}$ can be listed as
$\\{z_{1},z_{2},\dots,z_{q^{d-1}-1}\\}.$ For any nonzero element $z_{i}\in
V\setminus\\{0\\}$ with $1\leqslant i\leqslant q^{d-1}-1$, a binary sequence
${\mathbf{s}}_{i}$ associated to $z_{i}$ can be defined as follows:

${\mathbf{s}}_{i}=(s_{i,0},s_{i,1},\cdots,s_{i,q+t})\text{ with
}s_{i,j}=\text{Tr}(z_{i}(P_{j}))\text{ for }0\leqslant j\leqslant q+t.$

It is clear that each sequence ${\mathbf{s}}_{i}$ has length $q+1+t$ and the
family of binary sequences ${\mathcal{S}}=\\{{\mathbf{s}}_{i}:1\leqslant
i\leqslant q^{d-1}-1\\}$ has size $q^{d-1}-1$. We are going to show that this
family of binary sequences ${\mathcal{S}}$ has a low correlation and a large
linear complexity.

###  III-B Binary sequences with a low correlation via elliptic curves

In this subsection, we shall show that the family of binary sequences
${\mathcal{S}}=\\{{\mathbf{s}}_{i}:1\leqslant i\leqslant q^{d-1}-1\\}$
constructed as above has a low correlation.

######  Proposition III.3.

For any elements $z_{1},z_{2}\in V\setminus\\{0\\}$, one has

  * (1)

$(z_{1}+\tau(z_{2}))_{\infty}=Q+\tau(Q)$ for any automorphism $\tau\in{\rm
Gal}(E/F)\setminus\\{\rm id\\}$;

  * (2)

$(z_{1}+z_{2})_{\infty}=Q$, provided that $z_{1}\neq z_{2}$.

###### Proof.

For any elements $z_{1},z_{2}\in V\setminus\\{0\\}$, the pole divisors of
$z_{1}$ and $z_{2}$ are given by $(z_{1})_{\infty}=(z_{2})_{\infty}=Q.$ By
Lemma II.1, we have $\nu_{\tau(Q)}(\tau(z_{2}))=\nu_{Q}(z_{2})=-1$ and
$(\tau(z_{2}))_{\infty}=\tau(Q).$ It follows that
$z_{1}+\tau(z_{2})\in\mathcal{L}(Q+\tau(Q)).$ For any automorphism
$\tau\in{\rm Gal}(E/F)\setminus\\{\rm id\\}$, we have $\tau(Q)\neq Q$ by
Proposition III.1. Hence, we have

$(z_{1}+\tau(z_{2}))_{\infty}=Q+\tau(Q)$

from the Strict Triangle Inequality [26, Lemma 1.1.11].

Suppose that $z_{1}+z_{2}\in\mathbb{F}_{q}$. It follows that $z_{1}+z_{2}\in
V\cap\mathbb{F}_{q}=\\{0\\}$, i.e., $z_{1}=z_{2}$ which is a contradiction.
Thus, we obtain $(z_{1}+z_{2})_{\infty}=Q$. ∎

######  Theorem III.4.

Let $q=2^{n}$ be a prime power of $2$ for a positive integer $n$. Let
$E/\mathbb{F}_{q}$ be a cyclic elliptic function field defined over
$\mathbb{F}_{q}$ with $q+1+t$ rational places. Let $F$ be the fixed subfield
of $E$ with respect to its translation group $T_{E}$ generated by $\tau_{P}$.
Let $Q$ be a place of $E$ with $\deg(Q)=d\geqslant 2$ such that
$\gcd(d,q+1+t)=1$. Then there exists a family of binary sequences
${\mathcal{S}}=\\{{\mathbf{s}}_{i}:1\leqslant i\leqslant q^{d-1}-1\\}$ with
length $q+1+t$, size $q^{d-1}-1$ and correlation upper bounded by

$\operatorname{Cor}({\mathcal{S}})\leqslant(2d+1)\cdot\lfloor
2\sqrt{q}\rfloor+|t|.$

###### Proof.

Let $\sigma=\tau_{P}$ be the generator of the translation group $T_{E}$. We
begin by computing the autocorrelation of the family of binary sequences
${\mathcal{S}}=\\{{\mathbf{s}}_{i}:1\leqslant i\leqslant q^{d-1}-1\\}$. For
each $1\leqslant i\leqslant q^{d-1}-1$, the autocorrelation of
${\mathbf{s}}_{i}$ at delay $u$ with $1\leqslant u\leqslant q+t$ is given by

$\displaystyle A_{u}({\mathbf{s}}_{i})$
$\displaystyle=\sum_{j=0}^{q+t}(-1)^{s_{i,j}+s_{i,j+u}}=\sum_{j=0}^{q+t}(-1)^{\operatorname{Tr}(z_{i}(P_{j}))+\operatorname{Tr}(z_{i}(P_{j+u}))}$
$\displaystyle=\sum_{j=0}^{q+t}(-1)^{\operatorname{Tr}(z_{i}(P_{j}))+\operatorname{Tr}((\sigma^{-u}(z_{i}))(P_{j}))}=\sum_{j=0}^{q+t}(-1)^{\operatorname{Tr}((z_{i}+\sigma^{-u}(z_{i}))(P_{j}))}.$

The third equality follows from Lemma II.1. Since $z_{i}\in V\setminus\\{0\\}$
and $1\leqslant u\leqslant q+t$, we have
$(z_{i}+\sigma^{-u}(z_{i}))_{\infty}=Q+\sigma^{-u}(Q)$ by Proposition III.3.
Consider the Artin-Schreier curve $E_{i}=E(y)$ defined by

$y^{2}+y=z_{i}+\sigma^{-u}(z_{i}).$

By Proposition II.6, the place $Q$ and $\sigma^{-u}(Q)$ are the only ramified
places in $E_{i}/E$ and their different exponents are two. From the Hurwitz
genus formula, we have

$2g(E_{i})-2=2[2g(E)-2]+2\deg(Q)+2\deg(\sigma^{-u}(Q)).$

Hence, the genus of $E_{i}$ is $g(E_{i})=2\deg(Q)+1=2d+1$ for each $1\leqslant
i\leqslant q+t$. Let $N_{0}$ denote the cardinality of the set
$S_{0}=\\{0\leqslant j\leqslant
q+t:\operatorname{Tr}((z_{i}+\sigma^{-t}(z_{i}))(P_{j}))=0\\}$. Let $N_{1}$
denote the cardinality of the set $S_{1}=\\{0\leqslant j\leqslant
q+t:\operatorname{Tr}((z_{i}+\sigma^{-t}(z_{i}))(P_{j}))=1\\}$. It is clear
that

$N_{0}+N_{1}=q+1+t.$

If $j\in S_{0}$, then $P_{j}$ splits into two rational places in $E_{i}$ by
[15, Theorem 2.25] and [26, Theorem 3.3.7]. If $j\in S_{1}$, then $P_{j}$ is
inert in $E_{i}$ and any place of $E_{i}$ lying over $P_{j}$ has degree two.
Let $N(E_{i})$ be the number of rational points of the Artin-Schreier curve
$E_{i}$. Then we have

$N(E_{i})=2N_{0}.$

From the Serre bound, we have $|N(E_{i})-q-1|\leqslant g(E_{i})\cdot\lfloor
2\sqrt{q}\rfloor=(2d+1)\cdot\lfloor 2\sqrt{q}\rfloor.$ Hence, we obtain

$|A_{u}({\mathbf{s}}_{i})|=|N_{0}-N_{1}|=|2N_{0}-q-1-t|=|N(E_{i})-q-1-t|\leqslant(2d+1)\cdot\lfloor
2\sqrt{q}\rfloor+|t|.$

Now let us consider the cross-correlation of the family of sequences. For two
distinct sequences ${\mathbf{s}}_{i}$ and ${\mathbf{s}}_{j}$ in
${\mathcal{S}}$ for $1\leqslant i\neq j\leqslant q^{d-1}-1$, the cross-
correlation of ${\mathbf{s}}_{i}$ and ${\mathbf{s}}_{j}$ at delay $u$ with
$0\leqslant u\leqslant q+t$ is given by

$C_{u}({\mathbf{s}}_{i},{\mathbf{s}}_{j})=\sum_{k=0}^{q+t}(-1)^{\operatorname{Tr}(z_{i}(P_{k}))+\operatorname{Tr}(z_{j}(P_{k+u}))}=\sum_{k=0}^{q+t}(-1)^{\operatorname{Tr}((z_{i}+\sigma^{-u}(z_{j}))(P_{k}))}.$

Consider the Artin-Schreier curve $E_{i,j}=E(y)$ defined by

$y^{2}+y=z_{i}+\sigma^{-u}(z_{j}).$

If $1\leqslant u\leqslant q+t$, then we have
$(z_{i}+\sigma^{-u}(z_{j}))_{\infty}=Q+\sigma^{-u}(Q)$ by Proposition III.3.
From the Hurwitz genus formula and Proposition II.6, we have

$2g(E_{i,j})-2=2[2g(E)-2]+2\deg(Q)+2\deg(\sigma^{-u}(Q)).$

Hence, the genus of $E_{i,j}$ is $g(E_{i,j})=2\deg(Q)+1=2d+1$ for each
$1\leqslant i\leqslant q-1$. If $u=0$, then we have
$(z_{i}+\sigma^{-u}(z_{j}))_{\infty}=Q$ by Proposition III.3 and the genus of
$E_{i,j}$ is $g(E_{i,j})=\deg(Q)+1=d+1$ by Proposition II.6. Let
$N_{0}^{\prime}$ denote the cardinality of the set $\\{0\leqslant k\leqslant
q+t:\operatorname{Tr}((z_{i}+\sigma^{-u}(z_{j}))(P_{k}))=0\\}$. Let
$N_{1}^{\prime}$ denote the cardinality of the set $\\{0\leqslant k\leqslant
q+t:\operatorname{Tr}((z_{i}+\sigma^{-u}(z_{j}))(P_{k}))=1\\}$. It is clear
that

$N_{0}^{\prime}+N_{1}^{\prime}=q+1+t.$

Let $N(E_{i,j})$ be the number of rational points on the Artin-Schreier curve
$E_{i,j}$. Similarly, we have

$N(E_{i,j})=2N_{0}^{\prime}.$

From the Serre bound, we have

$|C_{u}({\mathbf{s}}_{i},{\mathbf{s}}_{j})|=|N_{0}^{\prime}-N_{1}^{\prime}|=|2N_{0}^{\prime}-q-1-t|=|N(E_{i,j})-q-1-t|\leqslant(2d+1)\cdot\lfloor
2\sqrt{q}\rfloor+|t|.$

The proof is completed by taking the maximum of $|A_{u}(\mathbf{s}_{j})|$ and
$|C_{u}(\mathbf{s}_{j},\mathbf{s}_{k})|$ among all possible $j,k$. ∎

######  Corollary III.5.

Let $q=2^{n}$ be a prime power for a positive integer $n$. Let $t$ be $0$, or
$\pm\sqrt{q}$ if $n$ is even, or $\pm\sqrt{2q}$ if $n$ is odd. Then there
exists a family of binary sequences with length $q+1+t$, size $q-1$ and
correlation upper bounded by

$\operatorname{Cor}({\mathcal{S}})\leqslant 5\cdot\lfloor
2\sqrt{q}\rfloor+|t|.$

###### Proof.

Let $E$ be a cyclic elliptic function field with $q+1+t$ rational places. By
Lemma III.2, there exists a place $Q$ of $E$ of degree $\deg(Q)=2$. Note that
$\gcd(2,q+1+t)=1$. Thus, this corollary follows from Proposition III.1,
Proposition III.3 and Theorem III.4. ∎

######  Corollary III.6.

Let $q=2^{n}$ be a prime power for a positive integer $n$. Let $t$ be an
integer satisfying one of three items in Lemma II.5 and
$\gcd(3,(-1)^{n}+1+t)=1$. Then there exists a family of binary sequences with
length $q+1+t$, size $q^{2}-1$ and correlation upper bounded by

$\operatorname{Cor}({\mathcal{S}})\leqslant 7\cdot\lfloor
2\sqrt{q}\rfloor+|t|.$

###### Proof.

Let $E$ be a cyclic elliptic function field with $q+1+t$ rational places. By
Lemma III.2 there exists a place $Q$ of $E$ of degree $\deg(Q)=3$ with
$\gcd(3,q+1+t)=1$. Then this corollary follows from Proposition III.1,
Proposition III.3 and Theorem III.4. ∎

###  III-C Binary sequences with a high linear complexity

In this subsection, we shall show that the family of binary sequences
constructed above has a high linear complexity similarly as [33, Theorem 3.2].

######  Theorem III.7.

The linear complexity of the family of binary sequence
${\mathcal{S}}=\\{{\mathbf{s}}_{i}:1\leqslant i\leqslant q^{d-1}-1\\}$
satisfies

$LC({\mathcal{S}})\geqslant\frac{q+1+2t-2(d+1)\sqrt{q}}{2d\sqrt{q}}.$

###### Proof.

It remains to prove that, for every $1\leqslant i\leqslant q^{d-1}-1$, the
linear complexity $\ell(\mathbf{s}_{i})$ of binary sequence $\mathbf{s}_{i}$
in $\mathcal{S}$ satisfies

$\ell(\mathbf{s}_{i})\geqslant\frac{q+1+2t-2(d+1)\sqrt{q}}{2d\sqrt{q}}.$

Let $\ell=\ell(\mathbf{s}_{i})$. Then there exist $\ell+1$ binary elements
$\lambda_{0}=1,,\cdots,\lambda_{\ell}=1\in\mathbb{F}_{2}$ such that

$\sum_{j=0}^{\ell}\lambda_{j}\cdot\operatorname{Tr}(z_{i}(P_{j+u}))=0$

for all integers $u\geqslant 0$. Let $\sigma=\tau_{P}$ be the generator of
translation group $T_{E}$. It is clear that

$\displaystyle\sum_{j=0}^{\ell}\lambda_{j}\cdot\operatorname{Tr}(z_{i}(P_{j+u}))$
$\displaystyle=\operatorname{Tr}\left(\sum_{j=0}^{\ell}\lambda_{j}z_{i}(P_{j+u})\right)$
$\displaystyle=\operatorname{Tr}\left(\sum_{j=0}^{\ell}\lambda_{j}z_{i}(\sigma^{j}(P_{u}))\right)$
$\displaystyle=\operatorname{Tr}\left(\sum_{j=0}^{\ell}\lambda_{j}\sigma^{-j}(z_{i})(P_{u})\right)$
$\displaystyle=\operatorname{Tr}(w_{i}(P_{u})),$

provided that $w_{i}=\sum_{j=0}^{\ell}\lambda_{j}\sigma^{-j}(z_{i})$. Hence,
we obtain

$(\operatorname{Tr}(w_{i}(P_{1})),\operatorname{Tr}(w_{i}(P_{2})),\cdots,\operatorname{Tr}(w_{i}(P_{q+1+t})))={\mathbf{0}}.$

Let $L_{i}=E(y)$ be the Artin-Schreier extension over $E$ generated by
$y^{2}+y=w_{i}$. From the Strict Triangle Inequality [26, Lemma 1.1.11], the
principal divisor of $w_{i}$ is given by
$(w_{i})_{\infty}=\sum_{j=0}^{\ell}\lambda_{j}\sigma^{-j}(Q)$. By Proposition
II.6, we have

$2g(L_{i})-2=2\cdot(2g(E)-2)+\deg(\operatorname{Diff}(L_{i}/E))\leqslant
2(\ell+1)\deg(Q).$

Hence, the genus of $L_{i}$ is upper bounded by

$g(L_{i})\leqslant d(\ell+1)+1.$

It follows that

$2(q+1+t)\leqslant N(L_{i})\leqslant q+1+[d(\ell+1)+1]\cdot\lfloor
2\sqrt{q}\rfloor.$

Hence, we obtain

$\ell\geqslant\frac{q+1+2t-(d+1)\cdot\lfloor 2\sqrt{q}\rfloor}{d\cdot\lfloor
2\sqrt{q}\rfloor}.\qed$

###  III-D Comparison

In [33], Xing et al. provided a construction of binary sequences with a low
correlation and a large linear complexity via elliptic curves. The family of
binary sequences given in [33] is different from our family of binary
sequences given in this manuscript, although both families of binary sequences
share the same length, the same upper bound of correlation and the same lower
bound of linear complexity from Theorem III.4, Theorem III.7 and [33, Theorem
5.5].

Let $d$ be a positive divisor which is relatively prime with $q+1+t$. There
are exactly $r=B_{d}/(q+1+t)$ places of $F=E^{T_{E}}$ of degree $d$ which
split completely in the extension $E/F$. Let $Q_{1},Q_{2},\cdots,Q_{r}$ be the
different places of $E$ with degree $d$ such that their restrictions are
pairwise distinct. Xing et al. chosen a unique element in every set
$\mathcal{L}(Q_{i})\setminus\mathbb{F}_{q}$ for $1\leqslant i\leqslant r$.
Hence, the size of their family of binary sequences is $r=B_{d}/(q+1+t)$.
However, we choose all nonzero elements in the complementary dual linear space
$V$ of $\mathbb{F}_{q}$ in the Riemann-Roch space
$\mathcal{L}(Q)=\mathbb{F}_{q}\oplus V$ for any fixed place of $Q$ of $E$ with
degree $d$ in this manuscript. By Lemma III.2, Corollary III.5 and Corollary
III.6, the size of the new family of binary sequences in this manuscript is
larger than the one given in [33].

##  IV An explicit construction of such a family of binary sequences

The previous section provides a theoretical framework for constructing a new
family of binary sequences with a low correlation and a high linear complexity
using elliptic function fields. In this section, we will make this
construction explicitly. We achieve this by determining the automorphisms and
rational places explicitly. Additionally, we present some numerical results
generated by using the software Sage. It turns out that the actual correlation
of our family of binary sequences is significantly lower than the estimated
upper bounds provided in Theorem III.4.

###  IV-A Algorithm

Now we provide an algorithm to generate the family of binary sequences with a
low correlation and a high linear complexity constructed in Section III. Such
an explicit construction of sequences with a low correlation can be given as
follows.

Algorithm: Explicit construction of sequences with a low correlation

  * •

Step 1: Input a prime power $q=2^{n}$ and a fixed integer $t$ given in Lemma
II.5.

  * •

Step 2: Find a cyclic elliptic function field $E/\mathbb{F}_{q}$ with exactly
$q+1+t$ rational places.

  * •

Step 3: Find a rational place $P=(\alpha_{1},\beta_{1})$ of $E$ with order
$q+1+t$ and calculate $P_{j}=[j]P=(\alpha_{j},\beta_{j})$ for $0\leqslant
j\leqslant q+t$ from [25, Group Law Algorithm 2.3].

  * •

Step 4: Let $\lambda=\frac{y-\beta_{1}}{x-\alpha_{1}}$ and
$v=\frac{\beta_{1}x-\alpha_{1}y}{x-\alpha_{1}}$. The automorphism $\tau_{P}$
of $E$ can be determined explicitly as

$\tau_{P}(x)=\lambda^{2}+a_{1}\lambda-
a_{2}-x-\alpha_{1},\quad\tau_{P}(y)=-(\lambda+a_{1})\tau_{P}(x)-v-a_{3}.$

Let $T_{E}$ be the translation group $E$ which is generated by $\tau_{P}$. Let
$F$ be the fixed subfield of $E$ with respect to $T_{E}$.

  * •

Step 5: Find the Riemann-Roch space of $\mathcal{L}(Q)$ for any place $Q$ with
degree $d\geqslant 2$ which is relatively prime to $q+1+t$ and determine the
vector space $V$ such that $\mathcal{L}(Q)=\mathbb{F}_{q}\oplus V$.

  * •

Step 6: Label all elements of $V\setminus\\{0\\}$ with $\\{z_{i}:1\leqslant
i\leqslant q^{d-1}-1\\}$.

  * •

Step 7: Output a family of sequences
${\mathcal{S}}=\\{{\mathbf{s}}_{i}:1\leqslant i\leqslant q^{d-1}-1\\}$ defined
by

${\mathbf{s}}_{i}=(s_{i,0},s_{i,1},\cdots,s_{i,q+t})\text{ with
}s_{i,j}={\operatorname{Tr}(z_{i}(P_{j}))}\text{ for }0\leqslant j\leqslant
q+t.$

  * •

Step 8: Output the correlation
$\operatorname{Cor}({\mathcal{S}})=\max\\{\\{|\sum_{k=0}^{q+t}(-1)^{s_{i,k}}(-1)^{s_{i,k+u}}|:1\leqslant
i\leqslant q^{d-1}-1,1\leqslant u\leqslant
q+t\\}\cup\\{|\sum_{k=0}^{q+t}(-1)^{s_{i,k}}(-1)^{s_{j,k+u}}|:1\leqslant i\neq
j\leqslant q^{d-1}-1,0\leqslant u\leqslant q+t\\}\\}.$

######  Remark IV.1.

Now we estimate the computational complexity on generating the family of
binary sequences from our construction. We consider the case when $d=2$. In
Step 2, to find an elliptic curve with a specific number of points, one would
typically generate random curves and count their points. One of the most
effective algorithms for counting points on elliptic curves over finite fields
is the Schoof-Elkies-Atkin (SEA) algorithm which takes $O((\log q)^{6})$ time
[31]. Thus the worst case complexity for Step 1 is $O(q^{5}(\log q)^{6})$. In
Step 3, for each element, we need to check if it is a generator which takes at
most $O(q)$, thus the total time is $O(q^{2}(\log q)^{2})$. In Step 7, we do
the evaluation which takes $O(q^{2}(\log q)^{3})$. In Step 8, we need to
compute the correlation between any two sequences as well as the cyclic shift,
thus the computational complexity is $O(q^{4})$. The complexity of other steps
are small, so they can be neglected. Thus, the overall complexity of
generating the family of binary sequences is $O(q^{5}(\log q)^{6})$.

###  IV-B Numerical results

In this subsection, we list some numerical results on the parameters of the
family of binary sequences obtained from elliptic curves with the help of Sage
software. We consider the case where $d=2$ in the following.

TABLE II: PARAMETERS OF OUR SEQUENCES WITH LENGTH EQUAL TO $q$ Field Size | Seq. Length | Family Size | Max Correlation
---|---|---|---
64 | 64 | 63 | 38
128 | 128 | 127 | 60
256 | 256 | 255 | 86
512 | 512 | 511 | 124
1024 | 1024 | 31 | 184
2048 | 2048 | 63 | 276
4096 | 4096 | 127 | 416

Let $q=2^{n}$ be a prime power. By Lemma II.5, there exists a cyclic elliptic
function field $E/{\mathbb{F}}_{q}$ with exactly $q$ rational places. By
Corollary III.6, we can construct a family of binary sequences with length
$q$, family size $q^{2}-1$ and correlation upper bounded by $14\sqrt{q}+1$. We
list some numerical results for this case in Table II. It turns out that the
real correlation of our family of binary sequence of length $q$ is much better
than the estimated upper bounds given in Theorem III.4.

Let $q=2^{n}$ be a prime power. Let $t$ be $\sqrt{q}$ if $n$ is even, or
$\sqrt{2q}$ if $n$ is odd. It is trivial that $\gcd(2,q+1+t)=1$. By Corollary
III.5, we can construct a family of binary sequences with family size $q-1$
and odd length $q+1+t$. We list some numerical results for this case in Table
III.

TABLE III: PARAMETERS OF OUR SEQUENCES OF ODD LENGTHS Field Size | Seq. Length | Family Size | Max Correlation
---|---|---|---
64 | 73 | 63 | 39
128 | 145 | 127 | 57
256 | 273 | 255 | 89
512 | 545 | 511 | 137
1024 | 1057 | 1023 | 191

In fact, we are able to generate more families of binary sequences with
flexible lengths $q+1+t$ with $t$ being an integer between $-2\sqrt{q}$ and
$2\sqrt{q}$.

##  V Conclusion

While significant progress has been made in constructing binary sequences with
good parameters, there remains a shortage of options for sequence lengths in
existing results. Therefore, it is imperative to develop sequences with
alternative lengths to accommodate diverse scenarios. Most existing families
of binary sequences use the cyclic structures inherent in the multiplicative
cyclic group of $\mathbb{F}_{p^{n}}$. Notably, recent work by [8, 10]
capitalized on the cyclic structure of $q+1$ rational points on projective
lines over $\mathbb{F}_{q}$ to design a family of binary sequences with a
length of $q+1$.

In this paper, we exploit the cyclic group structure of rational points of
elliptic curves to design a family of binary sequences of length $2^{n}+1+t$
with a low correlation and a high linear complexity. Our family size surpasses
that of [33]. Furthermore, our construction offers the potential to derive
additional sequences with a low correlation and other desirable parameters.

## References

  * [1] A. N. Akansu, R. Poluri, “Walsh-like nonlinear phase orthogonal codes for direct sequence CDMA communications,” IEEE Trans. Signal Processing, vol. 55, no. 7, pp. 3800–3806, Jul. 2007.
  * [2] S. Boztas and P. V. Kumar, “Binary sequences with Gold-like correlation but larger linear span,” IEEE Trans. Inf. Theory, vol. 40, no. 2, pp. 532–537, Feb. 1994.
  * [3] G. Cui, Y. He, P. Li, W. Wang, “Enhanced timing advanced estimation with symmetric Zadoff-Chu sequences for satellite systems,” IEEE Communication Letters, vol. 19, no. 5, pp. 747–750, May 2015.
  * [4] M. Eichler, Introduction to the Theory of Algebraic Numbers and Functions, New York: Academic, 1951.
  * [5] S. W. Golomb, Shift Register Sequences, Aegean Park Press, 1967.
  * [6] G. Gong, “New designs for signal sets with low cross correlation, balance property, and large linear span: $GF(p)$ case,” IEEE Trans. Inf. Theory, vol. 48, no. 11, pp. 2847–2867, Nov. 2002.
  * [7] L. Jin, D. Chen, L. Qian, J. Teng and S. Chen, “Construction of binary sequences with low correlation via multiplicative quadratic character over finite fields of odd characteristics,” IEEE Trans. Inf. Theory, vol. 67, no. 4, pp. 2236–2244, Apr. 2021.
  * [8] X. Hu, L. Jin, L. Ma and C. Xing, “Binary Sequences With a Low Correlation via Cyclotomic Function Fields of Odd Characteristic”, IEEE Trans Inf. Theory, vol. 70, no. 2, pp. 1397–1407, Feb. 2024.
  * [9] L. Jin, L. Ma and C. Xing, “Construction of optimal locally repairable codes via automorphism groups of rational function fields,” IEEE Trans. Inf. Theory, vol. 66, no. 1, pp. 210–221, Jan. 2020.
  * [10] L. Jin, L. Ma and C. Xing, “Binary sequences with a low correlation via cyclotomic function fields,” IEEE Trans. Inf. Theory, vol. 68, no. 5, pp. 3445–3454, May 2022.
  * [11] T. Kasami, “Weight distribution formula for some class of cyclic codes,” Coordinated Science Library R-285, 1996.
  * [12] P. Kim; A. Vanelli-Coralli, M. Villanti, R. Pedone, S. Cioni, M. Neri, C. Palestini, M. Papaleo, H. Lee and G. E. Corazza, “Direct sequence spectrum spreading techniques for next generation mobile broadband satellite services,” International Journal of Satellite Communications and Networking, vol. 28, pp.157–181, May 2010.
  * [13] A. M. Klapper, “D-form sequences: families of sequences with low correlation values and large linear spans,” IEEE Trans. Inf. Theory, vol. 41, no. 2, pp. 423–431, Feb. 1995.
  * [14] A. Lempel, M. Cohn, and W. Eastman, “A class of balanced binary sequences with optimal autocorrelation properties,” IEEE Trans. Inf. Theory, vol. 23, no. 1, pp. 38–42, Jan. 1977.
  * [15] R. Lidl and H. Niederreiter, Finite Fields, Encyclopedia of Mathematics and Its Applications, Cambridge University Press, Cambridge, 1983.
  * [16] L. Ma and C. Xing, “The group structures of automorphism groups of elliptic curves over finite fields and their applications to optimal locally repairable codes,” J. Combin. Theory Ser. A, vol. 193, 105686, Jan 2023.
  * [17] H. Niederreiter and C. Xing, Rational Points on Curves over Finite Fields: Theory and Applications, LMS 285, Cambridge, 2001.
  * [18] H. Niederreiter and C. Xing, “Sequences with high nonlinear complexity,” IEEE Trans. Inf. Theory, vol. 60, no. 10, pp. 6696–6701, Oct. 2014.
  * [19] J. No and P. V. Kumar, “A new family of binary pseudo random sequences having optimal periodic correlation properties and larger linear span,” IEEE Trans. Inf. Theory, vol. 35, no. 2, pp. 371–379, Feb. 1989.
  * [20] J.-S. No, H. Chung, H.-Y. Song, and K. Yang, “Trace representation of Legendre sequences of Mersenne prime period,” IEEE Trans. Inf. Theory, vol. 42, no. 6, pp. 2254–2255, Nov. 1996.
  * [21] J. Olsen, R. Scholtz, and L. Welch, “Bent-function sequences,” IEEE Trans. Inf. Theory, vol. 28, no. 6, pp. 858–864, Jun. 1982.
  * [22] K. G. Paterson, “Binary sequence sets with favorable correlations from difference sets and MDS codes,” IEEE Trans. Inf. Theory, vol. 44, no. 1, pp. 172-180, Jan. 1998.
  * [23] H.G. Ruck, A note on elliptic curves over finite fields, Math. Computation 149(1987), no. 179, pp. 301–304.
  * [24] J. J. Rushanan, “Weil sequences: A family of binary sequences with good correlation properties,” In 2006 IEEE International Symposium on Information Theory, pp. 1648–1652, July 2006.
  * [25] J. H. Silverman, The Arithmetic of Elliptic Curves, (Graduate Texts in Mathematics), vol. 106, Berlin, Germany: Springer Verlag, 1986.
  * [26] H. Stichtenoth, Algebraic Function Fields and Codes (Graduate Texts in Mathematics), vol. 254, Berlin, Germany: Springer Verlag, 2009.
  * [27] M. Su and A. Winterhof, “Autocorrelation of Legendre-Sidelnikov Sequences,” IEEE Trans. Inf. Theory, vol. 56, no. 4, pp. 1714–1718, Apr. 2010.
  * [28] X. Tang, T. Helleseth, L. Hu, and W. Jiang, “Two new families of optimal binary sequences obtained from quaternary sequences,” IEEE Trans. Inf. Theory, vol. 55, no. 4, pp. 1833–1840, Apr. 2009.
  * [29] X. Tang, P. Udaya, and P. Fan, “Generalized binary Udaya-Siddiqi sequences,” IEEE Trans. Inf. Theory, vol. 53, no. 3, pp. 1225–1230, Mar. 2007.
  * [30] P. Udaya and M. U. Siddiqi, “Optimal biphase sequences with large linear complexity derived from sequences over $\mathbb{Z}_{4}$,” IEEE Trans. Inf. Theory, vol. 42, no. 1, pp. 206–216, Jan. 1996.
  * [31] F. Vercauteren, “The SEA algorithm in characteristic $2$,” Preprint, 2000.
  * [32] W. C. Waterhouse, “Abelian varieties over finite fields,” Ann. Sci. Éc. Norm. Supér., série 4, tome 2, no. 4, pp. 521–560, Apr. 1969.
  * [33] C. Xing, P.V. Kumar and C. Ding, “Low-correlation, large linear span sequences from function fields”, IEEE Trans. Inf. Theory, vol. 49, no. 6, pp. 1439–1446, Jun. 2003.
  * [34] W. Zhang, C. Xie, E. Pasalic, “Large sets of orthogonal sequences suitable for applications in CDMA systems,” IEEE Trans. Inf. Theory, vol. 62, no. 6, pp. 3757–3767, Jun. 2016.
  * [35] Z. Zhou and X. Tang, “Generalized modified Gold sequences,” Designs, Codes and Cryptography, vol. 60, pp. 241–253, 2011.


# Convergent iLQR for Safe Trajectory Planning and Control
of Legged Robots

James Zhu, J. Joe Payne, and Aaron M. Johnson This material is based upon work
supported by the National Science Foundation under grant #CMMI-1943900.All
authors are with the Department of Mechanical Engineering, Carnegie Mellon
University, Pittsburgh, PA, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In order to perform highly dynamic and agile maneuvers, legged robots
typically spend time in underactuated domains (e.g. with feet off the ground)
where the system has limited command of its acceleration and a constrained
amount of time before transitioning to a new domain (e.g. foot touchdown).
Meanwhile, these transitions can instantaneously change the system’s state,
possibly causing perturbations to be mapped arbitrarily far away from the
target trajectory. These properties make it difficult for local feedback
controllers to effectively recover from disturbances as the system evolves
through underactuated domains and hybrid impact events. To address this, we
utilize the fundamental solution matrix that characterizes the evolution of
perturbations through a hybrid trajectory and its 2-norm, which represents the
worst-case growth of perturbations. In this paper, the worst-case perturbation
analysis is used to explicitly reason about the tracking performance of a
hybrid trajectory and is incorporated in an iLQR framework to optimize a
trajectory while taking into account the closed-loop convergence of the
trajectory under an LQR tracking controller. The generated convergent
trajectories recover more effectively from perturbations, are more robust to
large disturbances, and use less feedback control effort than trajectories
generated with traditional methods.

###### Index Terms:

Legged Robots, Trajectory Optimization, Robust Control

##  I Introduction

Legged robotics research has increasingly focused on enabling highly dynamic
and agile motions such as jumping, leaping, and landing [1, 2, 3, 4].
Implementing these capabilities reliably would improve legged robot
performance in applications such as extraterrestrial or urban environment
navigation where jumping up on ledges or leaping across chasms may be
necessary. However, jumping and leaping are dangerous maneuvers, with failure
often resulting in catastrophic outcomes for the robot.

What makes these actions challenging is that they induce trajectories that are
both hybrid and underactuated, which doubly contribute to the difficulty in
controlling legged robots. Broadly speaking, a system is hybrid if it
undergoes discrete changes in state and/or dynamics [5, 6], and it is
underactuated if there exists a direction of acceleration in state space that
can not be commanded by any valid input [7, Ch. 1.2]. Even when an
underactuated system is controllable, driving the system to a desired target
state may require significant time and control effort, neither of which may be
readily available. For instance, a robot jumping in the air can not
arbitrarily choose how much time it has until its feet touchdown on the
ground. This means that the controller needs to spend a lot of effort to
correct tracking errors prior to touchdown, or else discontinuous, unbounded
saltation effects [8] can cause arbitrarily large divergence if incoming
errors are not sufficiently mitigated, e.g. with grazing impacts. Increasing
control gains is one possible solution to improve stability, though that
strategy comes at a large drawback of worsening robustness in the face of
modelling errors and uncertainties[9, Ch. 13]. Instead, this work leverages
nonlinearities in continuous and hybrid dynamics that make some trajectories
easier to stabilize than others, even under equivalent feedback controllers.

Figure 1: Similar quadruped gaits tracked with equivalent LQR controllers
display enormous differences in final tracking performance. Results for 50
paired trials of trajectories sampled from an initial error covariance of
$10^{-2}$ in all directions are shown. The top blue trajectory was generated
with a standard (vanilla) iLQR algorithm, while the bottom green trajectory
was generated with our novel $\chi$-iLQR, which improves the average closed-
loop convergence by 92%. The horizontal distance between the displayed frames
is exaggerated for visual clarity. Only the convergence of the position states
is represented, and does not indicate convergence of the velocity states.

This paper presents a novel adaptation of the iLQR trajectory optimization
algorithm that improves closed-loop convergence under an equivalent feedback
controller (i.e. without changing the LQR controller weights), as demonstrated
in Fig. 1. Our simulation results show that this convergent iLQR ($\chi$-iLQR)
achieves three simultaneous improvements over standard iLQR: superior tracking
performance from initial perturbations, reduced feedback control effort over
the trajectory, and improved robustness to large initial errors. Compared to
existing methods, $\chi$-iLQR has two additional key strengths. Firstly, it is
based on an analysis that is simple to compute compared to methods such as
sum-of-squares. Additionally, $\chi$-iLQR captures the local tracking
performance of a closed-loop trajectory, which directly predicts experimental
results.

##  II Related Works

A strategy that has been used to enable dynamic yet precarious behaviors for
legged robots is leveraging highly accurate, complex models and full-body
trajectory optimization to plan precise motions [3, 10]. While these methods
incorporate feedback controllers to stabilize the generated trajectories,
there has been little focus on how these feedback controllers should be
designed to stabilize closed-loop systems under error and uncertainty.

Robust trajectory planning has been successfully implemented for smooth
systems like wheeled robots, with some recent results being adapted to hybrid
systems like legged robots. These works have focused on optimizing over
uncertainties in system dynamics, such as unknown disturbances and modelling
errors [11, 12, 13]. For example, [11] designs robust closed-loop trajectories
for smooth systems by optimizing the volume reduction of an ellipsoidal
disturbance set, but was not applied to hybrid systems. Risk-sensitive
planning and control is an alternate method that optimizes over the variance
of a cost distribution that evolves through the trajectory [14, 15]. Other
approaches present trajectory optimization algorithms for legged robots over
uncertain terrain [16] and compute a forward reachable set to bound closed-
loop errors [17]. Many of these methods require the distribution of errors to
be prespecified, which is not always clear how to tune. Additional actuation,
such as reaction wheels or tails, also relieves the difficulties of
underactuated systems [18, 19]. However, this comes with obvious tradeoffs of
increased cost, size, and weight.

Separately, consider the problem of quantifying the stability or convergence
properties of a system. A very popular method is Lyapunov analysis, where the
existence of a positive definite differentiable scalar function with negative
definite derivatives, called the Lyapunov function, can guarantee asymptotic
stability of the system [20]. Lyapunov functions can be difficult to compute,
particularly for hybrid systems, and can require methods such as sum-of-
squares [21] or machine learning [22] to be tractable. A similar strategy
known as control barrier functions, which restricts the system from entering
some set of undesirable states, has been successfully implemented on legged
robot hardware [23], but has the same drawback as Lyapunov functions.

A different strategy to analyze the stability and convergence of trajectories
is contraction analysis [24], which tracks the distance between two close
trajectories. If this distance monotonically decreases over the trajectory,
then the system is contractive and asymptotic stability can be guaranteed
[24]. Contraction analysis has been incorporated into path planning and
trajectory optimization algorithms on smooth systems [25, 26], but applying
contraction analysis to hybrid systems is difficult because many mechanical
hybrid systems are not contractive at hybrid events [27]. [28] loosened the
contraction criterion and optimized the stability of open-loop periodic orbits
using monodromy matrix analysis. Here, we extend that work by generalizing to
non-periodic trajectories under feedback control.

##  III Analysis and Planning for Hybrid Systems

This section defines a hybrid system and quantifies the performance of a
closed-loop hybrid trajectory using linearized variational equations. With
this analysis, we can generate a scalar measure of a trajectory’s convergence
which is then incorporated into a trajectory optimization framework.

###  III-A Hybrid Systems

Hybrid systems are a class of dynamical systems that consist of continuous
domains connected by hybrid events[5, 6]. Following the notation in [8], we
describe a hybrid system as a set of discrete modes
$\\{\mathrm{I},\mathrm{J},\ldots,\mathrm{K}\\}$, each with a domain
${D}_{\mathrm{I}}$ and a time-varying vector field $F_{\mathrm{I}}$.
$G_{(\mathrm{I},\mathrm{J})}$ is a guard that triggers a transition between
mode $\mathrm{I}$ and mode $\mathrm{J}$ and $R_{(\mathrm{I},\mathrm{J})}$ is
the reset map defining that transition.

An execution of a hybrid system [29] begins at an initial state $x_{0}\in
D_{\mathrm{I}}$. With input $u_{\mathrm{I}}(t,x)$, the system obeys the
dynamics $F_{\mathrm{I}}$ on $D_{\mathrm{I}}$. If the system reaches guard
surface $G_{(\mathrm{I},\mathrm{J})}$, the reset map
$R_{(\mathrm{I},\mathrm{J})}$ is applied and the system continues in domain
$D_{\mathrm{J}}$ under the corresponding dynamics defined by $F_{\mathrm{J}}$.
The flow $\phi(t,t_{0},x_{0},U)$ describes how the hybrid system evolves from
some initial time $t_{0}$ and state $x_{0}$ until some final time $t$ under
input sequence $U$.

###  III-B Linearized Variational Equations

For both continuous domains and hybrid transitions, linearized variational
equations can be constructed to characterize the evolution of perturbations
$\delta x$ [30]. In each continuous domain, the linearized variational
equation is discretized from timestep $i$ to $i+1$ and is $\delta
x_{i+1}\approx(A_{\mathrm{I}}-B_{\mathrm{I}}K_{\mathrm{I}})\delta x_{i}$ [30]
with $A_{I}$ and $B_{I}$ being the derivatives of the discretized dynamics in
mode $\mathrm{I}$ w.r.t. state $x_{i}$ and control inputs $u_{i}$,
respectively, and $K_{\mathrm{I}}$ are linear feedback gains. The feedback
term drops out for open-loop systems. For hybrid events, the analogous
variational equation is the saltation matrix $\Xi_{(\mathrm{I},\mathrm{J})}$,
which describes the transition between modes $\mathrm{I}$ and $\mathrm{J}$.
The saltation matrix is the first-order approximation of the change in state
perturbations from before the hybrid event at $\delta x(t^{-})$ to
perturbations after $\delta x(t^{+})$ [27], such that $\delta
x(t^{+})\approx\Xi_{(\mathrm{I},\mathrm{J})}\delta x(t^{-})$. This linear
approximation assumes that nearby trajectories undergo the same mode
transition. Computing the saltation matrix relies on the derivatives of the
reset and guard of the transition along with the dynamics in each mode, and is
detailed in [8].

###  III-C Convergence Measure

Consider a trajectory that begins at state $x_{0}=x(t_{0})$ at time $t_{0}$
and is executed until time $t_{f}$ where it arrives at state
$x_{f}=\phi(t_{f},t_{0},x_{0},U)$. The control objective is to bring any
nearby initial state $\bar{x}_{0}=x_{0}+\delta x_{0}$ towards the nominal
trajectory so that at time $t_{f}$,
$\bar{x}_{f}=\phi(t_{f},t_{0},\bar{x}_{0},\bar{U})=x_{f}+\delta x_{f}$ is
closer to $x_{f}$. To characterize the closeness of $\bar{x}_{f}$ and $x_{f}$,
we utilize the fundamental solution matrix, $\Phi$. Following [31], the
fundamental solution matrix is the linearized approximation $\delta
x_{f}\approx\Phi\delta x_{0}$ and represents the transformation of error from
the initial state to final state.

The fundamental solution matrix can be computed by sequentially composing the
linearized variational equations in each continuous domain ($\tilde{A}:=A-BK$)
and the saltation matrices ($\Xi$) at each hybrid event [28]. For a hybrid
trajectory with $N$ domains, the fundamental solution matrix is:

$\displaystyle\Phi$
$\displaystyle=\tilde{A}_{N}\Xi_{(N-1,N)}\dots\Xi_{(2,3)}\tilde{A}_{2}\Xi_{(1,2)}\tilde{A}_{1}$
(1)

Since the fundamental solution matrix captures the change in errors across a
trajectory, the singular values of $\Phi$ characterize error change along
principle axes of state space. The largest singular value, which is equivalent
to the induced 2-norm of $\Phi$, describes the evolution of the most divergent
direction of initial error $\delta x_{0}$. We define the convergence measure,
$\chi$ to be exactly this worst-case value:

$\displaystyle\chi=||\Phi||_{2}$ (2)

$\chi$ is a continuous measure of local convergence, where smaller values of
$\chi$ indicate stronger reduction of worst-case final errors. A value of
$\chi<1$ indicates errors in all directions will shrink.

###  III-D iLQR for Hybrid Systems

The iterative linear quadratic regulator (iLQR) is a trajectory optimization
method that also computes LQR feedback gains over the generated trajectory
[32]. iLQR is convenient because compared to other trajectory optimization
methods like direct collocation [33], it is less computationally intensive and
guarantees a feasible trajectory. We draw from recent work that adapts the
iLQR algorithm for use on hybrid dynamical systems [34, 35].

In brief, iLQR solves the optimal control problem over $N$ discretized
timesteps:

$\displaystyle\min_{U}\quad$
$\displaystyle\ell_{N}(x_{N})+\sum_{i=0}^{N-1}\ell_{i}(x_{i},u_{i})$ (3) where
$\displaystyle x_{0}=x(0)$ (4) $\displaystyle
x_{i+1}=\phi(t_{i+1},t_{i},x_{i},u_{i})$ $\displaystyle\forall i$ (5)

where $\ell_{i}(x_{i},u_{i})$ and $\ell_{N}(x_{N})$ represent the nonlinear
stage cost and terminal cost, respectively, $X:=\\{x_{0},x_{1},...,x_{N}\\}$
is a sequence of states with $x_{i}\in\mathbb{R}^{n}$ the system state at
timestep $i$ and $U:=\\{u_{0},u_{1},...,u_{N-1}\\}$ is a sequence of control
inputs with $u_{i}\in\mathbb{R}^{m}$ the control input at timestep $i$. We
also record the sequence of domains $M:=\\{D_{0},D_{1},...,D_{N}\\}$ with
$D_{i}$ the domain at timestep $i$ such that $x_{i}\in D_{i}$. $\phi$ is the
aforementioned flow of the trajectory.

iLQR computes gradient and Hessian information of the cost, which results in a
quadratic approximation of the cost function. As such, the state and terminal
costs can equivalently be simplified as quadratic functions such that the cost
function is simplified to:

$\displaystyle
J=x_{N}^{T}Q_{N}x_{N}+\sum_{i=0}^{N-1}x_{i}^{T}Q_{i}x_{i}+u_{i}^{T}R_{i}u_{i}$
(6)

with $Q_{i},Q_{N}\in\mathbb{R}^{n\times n}$ and $R_{i}\in\mathbb{R}^{m\times
m}$ all positive definite.

iLQR solves the optimal control problem by alternating between forward passes
that simulate the system under a given control input sequence, and backward
passes that solve for a new locally optimal control sequence. In the backward
pass, the value function, which is the optimal cost to go at any timestep, is
propagated through the trajectory in reverse, and gives locally optimal
feedforward inputs and feedback gains at each timestep. Computing the value
function relies on gradient and Hessian computations of the cost function and
Jacobians of the dynamics, which equates to computing the linearized
variational equations discussed in Sec. III-B. For much greater detail of iLQR
for hybrid systems, see [34, 35].

##  IV Convergent iLQR

Here we present a novel trajectory optimization algorithm called convergent
iLQR or $\chi$-iLQR, summarized in Algorithm 1. In this method, the
convergence measure is added to the cost function from (6) such that the
algorithm minimizes:

$\displaystyle J_{\chi}=\quad$ $\displaystyle
Q_{\chi}\chi+x_{N}^{T}Q_{N}x_{N}+\sum_{i=0}^{N-1}x_{i}^{T}Q_{i}x_{i}+u_{i}^{T}R_{i}u_{i}$
(7)

where $Q_{\chi}$ is a scalar weighting parameter. Since $\chi$ is solely a
function of states and inputs, iLQR uses gradient and Hessian information to
make a quadratic approximation compatible with the other cost terms.

Typically in iLQR, the cost function $J$ is evaluated after each forward pass,
since it is only dependent on the states and inputs of the most recent
trajectory. However, in this case the convergence measure portion of the cost
function is dependent on the feedback gains generated by the algorithm. This
means that the gradient and Hessian terms of the cost function rely on the
feedback gains that are being updated at every timestep in the backward pass.
Due to this, the cost function derivatives are highly coupled with the gains
and become convoluted to compute.

To resolve this, we propose executing two separate backward passes that each
compute a different set of gains. We do this to preserve the convergence
properties of iLQR, though other choices might be possible such as borrowing
the previous set of gains under the assumption that the gains do not change
significantly over iterations. First, the tracking backward pass computes the
feedback gains that will be used as the LQR tracking controller gains and to
compute the convergence measure. It is equivalent to the backward pass in
standard iLQR using the cost function $J$ (6), which solves the Riccati
equation for the most recent trajectory. With the gains generated in the
tracking backward pass $K_{t}$, the convergent cost function $J_{\chi}$ (7)
can be computed. The search backward pass takes $J_{\chi}$ from the tracking
backward pass and computes the gradients of the convergent cost function with
controller gains $K_{t}$. The feedforward inputs $k_{s}$ and the feedback
gains $K_{s}$ from this pass are used to search for an improved trajectory in
the forward pass.

Since $J_{\chi}$ is returned by the tracking backward pass, a line search is
performed after this function call to guarantee the reduction of the cost
function $J_{\chi}$. If the line search condition is not satisfied, the
forward pass and tracking backward pass are looped until the line search
condition is passed.

Within the search backward pass, iLQR requires computation of the gradient and
Hessian of $\chi$. The derivatives of $\chi$ can be computed by leveraging the
singular value decomposition of $\Phi=USV^{T}$ where $S$ is a diagonal matrix
of singular values and the columns of $U$ and $V$ are the left and right
singular vectors, respectively. $\chi$ is the largest singular value of $\Phi$
and let $u_{\chi}$ and $v_{\chi}$ be its corresponding left and right singular
vectors. Following [36], the derivative of $\chi$ with respect to the state at
timestep $i$ is:

$\displaystyle\frac{\partial\chi}{\partial x_{i}}$
$\displaystyle=u_{\chi}^{T}\frac{\partial\Phi}{\partial x_{i}}v_{\chi}$ (8)

$\frac{\partial\Phi}{\partial x_{i}}$ in turn can be computed by using the
product rule along with leveraging the fact that only $\tilde{A}_{i}$ and
$\Xi_{(i,i+1)}$ are functions of $x_{i}$, and all other $\tilde{A}$ and $\Xi$
terms have zero derivatives with respect to $x_{i}$. For notational brevity,
let:

$\displaystyle P_{i}$
$\displaystyle=\tilde{A}_{N}\cdots\Xi_{(i+1,i+2)}\tilde{A}_{i+1}$ (9)
$\displaystyle O_{i}$
$\displaystyle=\Xi_{(i-1,i)}\tilde{A}_{i-1}\cdots\Xi_{(1,2)}\tilde{A}_{1}$
(10)

such that $\Phi=P_{i}\Xi_{(i,i+1)}\tilde{A}_{i}O_{i}$ and $\frac{\partial
P_{i}}{\partial x_{i}}=\frac{\partial O_{i}}{\partial x_{i}}=0$. Thus:

$\displaystyle\frac{\partial\Phi}{\partial x_{i}}$
$\displaystyle=P_{i}\frac{\partial\Xi_{(i,i+1)}}{\partial
x_{i}}\tilde{A}_{i}O_{i}+P_{i}\Xi_{(i,i+1)}\frac{\partial\tilde{A}_{i}}{\partial
x_{i}}O_{i}$ (11)

Derivatives with respect to the input $u_{i}$ follow equivalently. To improve
computational efficiency, $O_{i}$ at each timestep can be computed recursively
during the forward pass and each $P_{i}$ can be computed recursively in the
tracking backward pass. Since the rollout does not yet have feedback gain
information, the initial $O$ values must be computed separately.

Algorithm 1 Convergent iLQR Algorithm

Initialize $U$, $Q_{\chi}$, $Q_{N}$, $Q_{i}$, $R_{i}$, $n_{\text{iterations}}$

$X$, $U$, $M$, $J\leftarrow$ Rollout($U$)

$K_{t}$, $J_{\chi}$, $P\leftarrow$ TrackingBP($X$, $U$, $M$, $J$)

$O\leftarrow$ ComputeO($X$, $U$, $M$, $K_{t}$)

for $i\leftarrow 1\text{ to }n_{\text{iterations}}$ do

$k_{s}$, $K_{s}\leftarrow$ SearchBP($X$, $U$, $M$, $J_{\chi}$, $O$)

repeat

$X$, $U$, $M$, $J$, $O\leftarrow$ ForwardPass($X$, $U$, $M$, $k_{s}$, $K_{s}$)

$K_{t}$, $J_{\chi}$, $P\leftarrow$ TrackingBP($X$, $U$, $M$, $J$)

until LineSearchIsSatisfied($J_{\chi}$) return $X$, $U$, $M$, $K_{t}$

Because the scalar $\chi$ is derived from the norm of the matrix $\Phi$, the
gradient of $\chi$ relies on computing a 3-dimensional tensor of
$\tilde{A}_{i}$ and $\Xi$ derivatives, and the Hessian of $\chi$ is computed
from a 4-dimensional tensor of matrix second derivatives. While recent work
has enabled faster computation of second derivatives of dynamics [37] which
can aid in the computation of the 3-D tensor derivatives, computing 4-D tensor
derivatives is generally untenable. Instead, numerical methods like finite
differences for gradients and BFGS [38] for Hessians can perform at reasonable
speed. In order to approach real-time computation, it is likely that the full
Hessian of $\chi$ is not necessary to find an appropriate search direction and
that a partial computation or even leaving out the Hessian completely is
sufficient to compute optimal trajectories. Future work will address this gap.
Nonetheless, the algorithm in its current form can still be useful for offline
planning for trajectories that are expected to have a high degree of risk,
such as leaping across ledges or traversing narrow beams. In real-world
applications, it can be acceptable for a robot to pause and plan a safe
trajectory before executing these dangerous maneuvers.

TABLE I: Rocket Hopper Convergence Measure and Simulation Results Trial | LQR Parameters | Convergence Measure | Mean Simulated Error Ratio | Mean Simulated Feedback Effort
---|---|---|---|---
$Q_{\chi}$ | $Q_{N}$ | $R_{i_{\text{air}}}$ | $R_{i_{\text{stance}}}$ | Vanilla |  $\chi$-iLQR | %Difference | Vanilla |  $\chi$-iLQR | %Difference | Vanilla |  $\chi$-iLQR | %Difference
1 | 50 | 500$I$ | $0.01I$ | $0.1I$ | 1.01 | 0.71 | -29.70% | 0.42 | 0.33 | -21.72% | $1.74\cdot 10^{-5}$ | $1.6\cdot 10^{-5}$ | -7.16%
2 | 50 | 800$I$ | $0.005I$ | $0.01I$ | 0.78 | 0.51 | -34.50% | 0.32 | 0.24 | -25.74% | $4.66\cdot 10^{-5}$ | $3.25\cdot 10^{-5}$ | -30.31%
3 | 50 | 250$I$ | $0.02I$ | $0.05I$ | 1.14 | 0.94 | -17.70% | 0.49 | 0.45 | -8.00% | $7.12\cdot 10^{-6}$ | $7.05\cdot 10^{-6}$ | -1.01%
4 | 75 | 500$I$ | $0.01I$ | $0.01I$ | 1.01 | 0.68 | -33.24% | 0.41 | 0.32 | -21.73% | $1.89\cdot 10^{-5}$ | $1.62\cdot 10^{-5}$ | -14.20%

##  V Examples and Results

In this section, we demonstrate the convergence improvements of our method on
a spring hopper system and a planar quadruped robot model. Simulation results
show that the improved convergence measure correlates with an improvement in
average tracking performance, robustness to large disturbances, and feedback
control effort. Both examples were implemented in MATLAB, with forward
simulations using the `ode113` function. Cost function gradients were computed
using (11), derivatives of $\tilde{A}_{i}$ and $\Xi_{(i,i+1)}$ were computed
with finite differences, and Hessians were computed with BFGS.

###  V-A Rocket Hopper

####  V-A1 Rocket Hopper Model

This system is made up of a point mass body with a single massless spring leg.
The state of the hopper is characterized by the positions $x_{B}$, $y_{B}$ of
the body, the angle $\theta$ of the leg and their derivatives $\dot{x}_{B}$,
$\dot{y}_{B}$, $\dot{\theta}$ such that the full state is a $6\times 1$
vector. The system has two domains: an aerial phase $D_{1}$ and a stance phase
$D_{2}$. Taking a constant ground height at zero gives a touchdown guard
function $g_{(1,2)}$ that is the height of the foot and a liftoff guard
function $g_{(2,1)}$ that is the ground reaction force applied by the spring
leg. Both reset maps $R_{(1,2)}$ and $R_{(2,1)}$ are identity since position
and velocity are continuous.

The system has two inputs: a hip actuation and an actuation in the direction
of the leg. In the air, this allows the hopper to rotate the leg around the
body and exert a propulsion in the direction of the leg, somewhat akin to a
rocket, though this force can approximate forces from other legs or actuators.
A small rotor inertia in the air ensures the dynamics are well-conditioned
when controlling the massless leg. In stance, the hip torque and rocket force
exert ground reaction forces on the body. The body mass of the hopper was
chosen as 1 kg, spring constant as 250 N/m, and resting leg length as 0.75 m.

####  V-A2 Rocket Hopper Results

The objective for this system is to begin in the air at rest with a height of
2 m and end in the air at rest with the same height displaced 0.2 m
horizontally. The system is given 1.5 s for this trajectory.

We generated four trials of paired trajectories with varied weighting
parameters, shown in Table I, and compared the performance of the standard
(vanilla) iLQR method (where $Q_{\chi}=0$) to $\chi$-iLQR. There is no
reference trajectory to track, so $Q_{i}$ is zero for all trials.

For 3 of these trials, vanilla iLQR generated a trajectory with $\chi>1$,
meaning the worst-case error direction was expansive, see Table I. $\chi$-iLQR
decreases every convergence measure to below 1 so that all error directions
are reduced over the trajectory. On average, $\chi$-iLQR decreased $\chi$ by
28.79% compared to the vanilla method.

To validate these trajectories, each closed-loop trajectory was simulated 100
times with equivalent small random initial perturbations in both positions and
velocities with covariance matrix $\text{cov}(X_{0})=10^{-4}I$. A small
covariance was chosen so that the linearizations assumed in the convergence
measure and LQR control are valid. For each simulation run, the initial error
$\delta x_{0}$ and the final error $\delta x_{f}$ were recorded, along with
the sequence of control inputs $V:=\\{v_{0},v_{1},...,v_{N-1}\\}$. Note that
these inputs are distinct from the nominal feedforward inputs to the system
$U$ because there is additional feedback effort exerted by the actuators.

Two values were recorded during each simulation run to characterize the
convergence properties of the trajectories. The first is the error ratio,
$E=\frac{||\delta x_{f}||_{2}}{||\delta x_{0}||_{2}}$ defined as the ratio of
the final error 2-norm to the initial error 2-norm. A lower error ratio means
better tracking performance, and $E<1$ indicates a net reduction in error on
average. The second value is the feedback effort,
$F=\sum_{i=0}^{N-1}{(v_{i}-u_{i})^{2}}$ which is the sum of squares of the
difference between $V$ and $U$.

Table I shows that the simulation results support our assertion that an
improved convergence measure correlates with an improved mean tracking
performance and feedback effort. The mean error ratio and feedback effort over
the 100 simulations were both lower for trajectories generated with
$\chi$-iLQR. The average improvement over the four trials was 19.30% for mean
error ratio and 13.17% for feedback effort.

None of the simulated runs had an error ratio greater than one, which is
sensible since the worst-case direction occurs with probability zero. However,
even if none of the sampled initial errors aligned exactly with the worst-case
direction predicted by the fundamental solution matrix, nearby initial error
directions still see improvement in convergence, which explains the
improvement in mean simulated error ratio.

###  V-B Planar Quadruped

Here we demonstrate the improvements of $\chi$-iLQR on a more complex robot
model akin a standard quadruped robot. The model is simplified as a planar
quadruped, meaning that all movement occurs in the sagittal plane and the
left-right pairs of legs are constrained to move identically.

####  V-B1 Planar Quadruped Model

In the sagittal plane, we can model the robot with 7 positional states.
$x_{B}$, $y_{B}$, $\theta_{B}$ are the position and orientation of the body.
The front and back sets of legs each have two states for the hip angle
$\alpha_{f},\alpha_{b}$ and knee angle $\beta_{f},\beta_{b}$. Thus the full
state is dimension 14.

This system has four domains: the aerial domain $D_{1}$, front stance domain
$D_{2}$, back stance domain $D_{3}$, and full stance domain $D_{4}$. The
impact guard function is the height of the foot and the guard function for
liftoff is the vertical ground reaction force. The dynamics of the robot body
in the aerial phase follow ballistic motion, while the legs are simplified to
be massless while including the aforementioned rotor inertia. The impact reset
map for each foot consists of a discrete update to the hip and knee
velocities, while the body states are unchanged due to the massless legs. The
liftoff reset map is identity. The input vector for this system is
4-dimensional to actuate the hip and knee joints.

In this model, parallel torsion springs are added to the knee joints. Parallel
joint springs have been utilized to mimic tendons found in animals [39] that
increase the energy efficiency of legged locomotion [40, 41]. Due to the
resonance of the natural spring dynamics, controlling these systems requires
special care [42, 43]. For example, [44] solved for optimal gait timings to
leverage resonant spring frequencies. These spring models of legged robots are
good candidates for $\chi$-iLQR because the dynamics of the stance phase
depend strongly on the leg configuration at touchdown. Thus, a small error in
leg states at touchdown can have a large effect on tracking performance.

The inertial and dimensional properties were chosen to match the Ghost
Robotics Spirit 40 quadruped. The added torsional knee spring has a spring
constant $75\text{ N}\cdot\text{m}\cdot\text{rad}^{-1}$ and rest angle 1.2
rad.

####  V-B2 Planar Quadruped Results

The trajectory optimization task for the planar quadruped is to generate a
gait with a forward velocity of 0.25 m/s. The robot begins in the air with a
body height of 0.3 m. The hip joints begin at an angle of 0.6 rad and the knee
joints begin at 1.2 rad. The terminal target state is translated 0.0875 m in
the x-direction from the initial state. The trajectory is given 0.35 s to
execute. We choose to set a constant input weight of $R_{i}=5\cdot 10^{-4}I$.
The terminal weight is $Q_{N}=500I$ and the convergence weight for $\chi$-iLQR
is $Q_{\chi}=1$.

We set up a similar experiment to the prior example, with the addition of
simulating over a range of covariance magnitudes. This is done to evaluate the
basin of attraction of each trajectory over larger initial errors that
introduce greater nonlinear effects. The two trajectories were evaluated with
6 sets of 100 paired simulation runs with random initial error covariance
magnitudes of $10^{-4}$, $5\cdot 10^{-4}$, $10^{-3}$, $5\cdot 10^{-3}$,
$10^{-2}$, and $5\cdot 10^{-2}$ in each direction. The lowest covariance
magnitude of $10^{-4}$ approximates local linear behavior well, while $5\cdot
10^{-2}$ is the maximum magnitude before some trials begin with the robot’s
feet below the ground.

Table II shows the vanilla cost (6), convergence measure, mean simulated error
ratio, and mean simulated feedback effort of the two trajectories at the
covariance magnitude $10^{-4}$. As expected, the vanilla cost of the
convergent trajectory increases since its optimizing for a different cost
function, while the convergence measure and simulation values improve. We
argue that in dynamic legged locomotion, a costlier nominal trajectory can
often be worth an improvement in the trajectory’s robustness.

TABLE II: Mean convergence results for the quadruped model for vanilla and $\chi$-iLQR with covariance magnitude $10^{-4}$ | Vanilla |  $\chi$-iLQR | %Difference
---|---|---|---
Vanilla Cost | 65.58 | 74.39 | +13.43%
Convergence Measure | 60.52 | 41.35 | -31.68%
Mean Simulated Error Ratio | 6.66 | 4.78 | -28.23%
Mean Simulated Feedback Effort | 0.016 | 0.013 | -16.56%

The mean simulated error ratio for the convergent trajectory at this small
covariance magnitude was 28.23% less and the mean simulated feedback effort
was 16.56% less. Fig. 2 displays a histogram of the error ratio for each of
the trials, with the convergent trajectory having improved performance.

As the magnitude of initial errors grows, the performance of the LQR tracking
controller becomes worse due to the increase in nonlinear effects. Fig. 3
shows the simulation results for each trajectory over a range of initial error
covariance magnitudes. Each pair of lines indicates the success rate of the
respective closed-loop trajectories at maintaining error ratios of less than
50, 10, and 5 respectively. An error ratio of greater than 50 is
representative of a catastrophic failure, which the vanilla trajectory
encounters at a covariance magnitude of $5\cdot 10^{-4}$, while the convergent
trajectory first experiences a failure at covariance magnitude $10^{-2}$. This
difference in performance suggests the convergent trajectory is more robust to
larger initial errors and nonlinearities.

Figure 2: Histogram of error ratio for 100 paired simulated trials of the
quadruped model with small initial perturbations. Error ratio is the 2-norm of
final errors divided by the 2-norm of initial errors.

Even with the $\chi$-iLQR convergence improvements, the controller does not
reduce errors in all directions. The simulation results show there was usually
some error growth, which is reasonable since the body dynamics are fully
unactuated in the aerial phase and the system undergoes multiple hybrid
events. A combination of higher feedback gains and a global footstep planner
could be able to grant this system full convergence. Even so, this work can be
valuable to ensure that the system does not diverge too far from its target
trajectory between iterations of a global planner.

Figure 3: Plots showing success of LQR controllers at tracking vanilla ($V$)
and convergent ($\chi$) trajectories over various initial perturbation
covariances. $V_{50}$ and $\chi_{50}$ indicate proportion of trials where
error ratio was below 50, $V_{10}$ and $\chi_{10}$ below 10, and $V_{5}$ and
$\chi_{5}$ below 5.

##  VI Conclusion

In this work, we present a novel trajectory optimization method, $\chi$-iLQR
that optimizes over the worst-case error growth of a hybrid trajectory. This
method is based on the fundamental solution matrix, which maps the evolution
of perturbations through a trajectory. Incorporating the saltation matrix into
the fundamental solution matrix allows for straightforward handling of hybrid
events. The simulation results presented on two legged robot models
demonstrate that this method produces trajectories with improved tracking
performance, decreased feedback actuation effort, and improved robustness to
large perturbations. Even for a quadrupedal trajectory that was very difficult
to track, $\chi$-iLQR produced a trajectory that was superior at avoiding
failures.

Following this work, we aim to apply these principles to robot hardware and
demonstrate robust performance for maneuvers such as leaping and flipping. We
also aim to apply this work to other hybrid systems like robots that undergo
stick-slip transitions. This work and its extensions will further enable
robots to navigate complex, uncertain environments and unlock worlds for
robots to explore.

## References

  * [1] A. M. Johnson and D. E. Koditschek, “Toward a vocabulary of legged leaping,” in _IEEE Intl. Conference on Robotics and Automation_ , Karlsruhe, Germany, May 2013, pp. 2553–2560.
  * [2] H. Kolvenbach, E. Hampp, P. Barton _et al._ , “Towards jumping locomotion for quadruped robots on the moon,” _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2019.
  * [3] C. Nguyen and Q. Nguyen, “Contact-timing and trajectory optimization for 3d jumping on quadruped robots,” _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2022.
  * [4] Z. Li, X. B. Peng, P. Abbeel _et al._ , “Robust and versatile bipedal jumping control through multi-task reinforcement learning,” _arXiv preprint arXiv:2302.09450_ , 2023.
  * [5] A. van der Schaft and J. Schumacher, “Complementarity modeling of hybrid systems,” _IEEE Transactions on Automatic Control_ , 1998.
  * [6] J. Lygeros, K. H. Johansson, S. Simic _et al._ , “Dynamical properties of hybrid automata,” _IEEE Transactions on Automatic Control_ , 2003.
  * [7] R. Tedrake, _Underactuated Robotics_ , 2022. [Online]. Available: http://underactuated.mit.edu
  * [8] N. J. Kong, J. J. Payne, J. Zhu, and A. M. Johnson, “Saltation matrices: The essential tool for linearizing hybrid dynamical systems,” _arXiv preprint arXiv:2306.06862_ , 2023.
  * [9] M. A. Haidekker, _Linear Feedback Controls_ , M. A. Haidekker, Ed. Elsevier, 2020.
  * [10] J. Norby, A. Tajbakhsh, Y. Yang, and A. M. Johnson, “Adaptive complexity model predictive control,” _arXiv preprint arXiv:2209.02849_ , 2022.
  * [11] Z. Manchester and S. Kuindersma, “Robust direct trajectory optimization using approximate invariant funnels,” _Auton. Robots_ , 2019.
  * [12] J. Morimoto, G. Zeglin, and C. Atkeson, “Minimax differential dynamic programming: application to a biped walking robot,” _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2003.
  * [13] T. Lew, R. Bonalli, and M. Pavone, “Chance-constrained sequential convex programming for robust trajectory optimization,” _European Control Conference_ , 2020.
  * [14] B. Hammoud, M. Khadiv, and L. Righetti, “Impedance optimization for uncertain contact interactions through risk sensitive optimal control,” _IEEE Robotics and Automation Letters_ , 2021.
  * [15] M. Ahmadi, X. Xiong, and A. D. Ames, “Risk-sensitive path planning via CVaR barrier functions: Application to bipedal locomotion,” _IEEE Control Systems Letters_ , 2021.
  * [16] L. Drnach and Y. Zhao, “Robust trajectory optimization over uncertain terrain with stochastic complementarity,” _IEEE Robotics and Automation Letters_ , 2021.
  * [17] S. Kousik, S. Vaskov, F. Bu _et al._ , “Bridging the gap between safety and real-time performance in receding-horizon trajectory design for mobile robots,” _The International Journal of Robotics Research_ , 2020.
  * [18] C.-Y. Lee, S. Yang, B. Bokser, and Z. Manchester, “Enhanced balance for legged robots using reaction wheels,” _IEEE International Conference on Robotics and Automation_ , 2023.
  * [19] Y. Yang, J. Norby, J. K. Yim, and A. M. Johnson, “Proprioception and tail control enable extreme terrain traversal by quadruped robots,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2023\.
  * [20] H. Khalil, _Nonlinear Systems_. Pearson, 2002.
  * [21] A. Papachristodoulou and S. Prajna, “On the construction of Lyapunov functions using the sum of squares decomposition,” _IEEE Conference on Decision and Control_ , 2002.
  * [22] S. Chen, M. Fazlyab, M. Morari _et al._ , “Learning Lyapunov functions for hybrid systems,” _International Conference on Hybrid Systems: Computation and Control_ , 2021.
  * [23] R. Grandia, A. J. Taylor, A. D. Ames, and M. Hutter, “Multi-layered safety for legged robots via control barrier functions and model predictive control,” _IEEE International Conference on Robotics and Automation_ , 2021.
  * [24] W. Lohmiller and J. Slotine, “On contraction analysis for non-linear systems,” _Automatica_ , 1998.
  * [25] A. M. Johnson, J. E. King, and S. Srinivasa, “Convergent planning,” _IEEE Robotics and Automation Letters_ , 2016.
  * [26] N. Kong and A. M. Johnson, “Optimally convergent trajectories for navigation,” in _International Symposium on Robotics Research_ , October 2019\.
  * [27] S. A. Burden, T. Libby, and S. D. Coogan, “On contraction analysis for hybrid systems,” _arXiv preprint arXiv:1811.03956_ , 2018.
  * [28] J. Zhu, N. Kong, G. Council, and A. Johnson, “Hybrid event shaping to stabilize periodic hybrid orbits,” _IEEE International Conference on Robotics and Automation_ , 2022.
  * [29] A. M. Johnson, S. A. Burden, and D. E. Koditschek, “A hybrid systems model for simple manipulation and self-manipulation systems,” _The International Journal of Robotics Research_ , 2016.
  * [30] M. Hirsch, S. Smale, and R. Devaney, _Differential Equations, Dynamical Systems, and an Introduction to Chaos_. Academic Press, 2012.
  * [31] R. Leine and H. Nijmeijer, _Dynamics and Bifurcations of Non-Smooth Mechanical Systems_. Springer-Verlag Berlin Heidelberg, 2004.
  * [32] W. Li and E. Todorov, “Iterative linear quadratic regulator design for nonlinear biological movement systems,” _International Conference on Informatics in Control, Automation and Robotics_ , 2004.
  * [33] M. Kelly, “An introduction to trajectory optimization: How to do your own direct collocation,” _SIAM Review_ , 2017.
  * [34] N. J. Kong, G. Council, and A. M. Johnson, “iLQR for piecewise-smooth hybrid dynamical systems,” _IEEE Conference on Decision and Control_ , 2021.
  * [35] N. Kong, C. Li, G. Council, and A. M. Johnson, “Hybrid iLQR model predictive control for contact implicit stabilization on legged robots,” _IEEE Transactions on Robotics_ , 2023, to appear. Also available at arXiv:2207.04591 [cs.RO].
  * [36] J. Townsend, “Differentiating the singular value decomposition,” Tech. Rep., 2016\. [Online]. Available: https://j-towns.github.io/papers/svd-derivative.pdf
  * [37] S. Singh, R. P. Russell, and P. M. Wensing, “On second-order derivatives of rigid-body dynamics: Theory & implementation,” _arXiv preprint arXiv:2302.06001_ , 2023.
  * [38] C. G. Broyden, “The convergence of a class of double-rank minimization algorithms 1. General considerations,” _IMA Journal of Applied Mathematics_ , 1970.
  * [39] K. Endo and H. Herr, “A model of muscle-tendon function in human walking at self-selected speed,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , 2014.
  * [40] M. Scheint, M. Sobotka, and M. Buss, “Optimized parallel joint springs in dynamic motion: Comparison of simulation and experiment,” _IEEE RAS & EMBS International Conference on Biomedical Robotics and Biomechatronics_, 2010\.
  * [41] K. Mombaur, “Using optimization to create self-stable human-like running,” _Robotica_ , 2009.
  * [42] D. Lakatos, K. Ploeger, F. Loeffl _et al._ , “Dynamic locomotion gaits of a compliantly actuated quadruped with slip-like articulated legs embodied in the mechanical design,” _IEEE Robotics and Automation Letters_ , 2018.
  * [43] J. Ackerman and J. Seipel, “Energy efficiency of legged robot locomotion with elastically suspended loads,” _IEEE Transactions on Robotics_ , 2013.
  * [44] M. J. Pollayil, C. D. Santina, G. Mesesan _et al._ , “Planning natural locomotion for articulated soft quadrupeds,” _IEEE International Conference on Robotics and Automation_ , 2022.