# WAVELET CHANNEL ATTENTION MODULE WITH A FUSION NETWORK FOR SINGLE IMAGE
DERAINING

Single image deraining is a crucial problem because rain severely degenerates
the visibility of images and affects the performance of computer vision tasks
like outdoor surveillance systems and intelligent vehicles. In this paper, we
propose the new convolutional neural network (CNN) called the wavelet channel
attention module with a fusion network. Wavelet transform and the inverse
wavelet transform are substituted for down-sampling and up-sampling so feature
maps from the wavelet transform and convolutions contain different frequencies
and scales. Furthermore, feature maps are integrated by channel attention. Our
proposed network learns confidence maps of four sub-band images derived from
the wavelet transform of the original images. Finally, the clear image can be
well restored via the wavelet reconstruction and fusion of the low-frequency
part and high-frequency parts. Several experimental results on synthetic and
real images present that the proposed algorithm outperforms state-of-the-art
methods.

Index Terms—  Wavelet transform, single image deraining, fusion, channel
attention, convolutional neural network

##  1 Introduction

Single image deraining is a crucial problem because rain occludes the
background scene, appears in different locations and decreases the performance
of computer vision tasks like outdoor surveillance systems and intelligent
vehicles [1]. The rain streaks can be seen as linear noises which may vary in
size, direction and density. Moreover, in real cases, when rain accumulation
is dense, the individual streaks cannot be observed clearly. Accumulated rain
streaks reduce the visibility in a manner more similar to fog, and create a
haze-like phenomenon in the background. This foggy phenomenon can be described
as the haze model [2], and the whole model [3] is written as

$J=T\odot(I+\sum_{i}^{n}S_{i})+(1-T)\odot A$ (1)

where $J$ means the observed rain streak image, $I$ means the clear image,
$S_{i}$ means a rain streak layer with the same direction, $n$ means the
maximum number of layers, $A$ means the global atmospheric light, $T$ means
the atmospheric transmission and $\odot$ denotes element-wise multiplication.

Several methods try to analyze visual priors to capture deterministic and
statistical properties of rainy images [4, 5, 6]. However, these methods tend
to introduce undesirable artifacts, since their handcrafted priors from human
observations do not always hold in diverse real-world images. Instead of
applying handcrafted visual priors, recently, deep-learning-based methods [3,
7, 8, 9] are proposed, and these methods usually perform more accurate than
conventional priors-based methods with significant performance gains.

We also consider neural networks for single image deraining since the
deraining model is a crude approximation. CNN-based model can learn and
capture more detailed features from rainy images. Similar to other image
restoration tasks like image deblurring [7], and image dehazing [10], image
deraining can be modeled as an image-to-image mapping problem. From previous
studies [8, 9, 10], low-level features (e.g., edge and frequency) are more
important than high-level features like attribute, texture. To extract these
low-level features, we apply the wavelet transform and propose the wavelet
channel attention module (WCAM) with a fusion network. First, our network
replaces the down-sampling and up-sampling with the discrete wavelet transform
(DWT) and the inverse discrete wavelet transform (IDWT). The network further
captures the various frequency features and bi-orthogonal property of the DWT
for signal recovery. Second, the channel attention [11] module selectively
emphasizes interdependent channel maps by integrating associated features
among all channel maps. Thus, we combine the DWT and the channel attention so
that intermediate feature maps with different frequencies are integrated
effectively. Third, as demonstrated in [8, 12], fusing various levels of
features is beneficial for many computer vision tasks. The DWT is seen as the
fusion of high-frequency and low-frequency images from original images.
Inspired by it, our input is four sub-band images from the DWT, and the output
is four confidence maps that determine the importance of different sub-band
images and fuse them to reconstruct the clear image by the IDWT. In summary,
Our network is encoder-decoder architecture. The proposed WCAMs replace
convolutions in the encoder, and corresponding inverse wavelet channel
attention modules (IWCAMs) also replace convolutions in the decoder.

This paper makes the following contributions: (i) We propose a novel end-to-
end WCAM with a fusion network that captures frequency features and fuses four
sub-band images for single image deraining. (ii) Several experiments show that
the proposed network obtains much better performance than previous state-of-
the-art deraining methods, even in the presence of large rain streaks and rain
streak accumulation.

##  2 Related Work

###  2.1 Single image deraining

Single image deraining is an ill-posed problem and seen as the denoising
problem. Early methods propose hand-craft priors to estimate distributions of
rain streaks and remove them. In [4], Chen and Hsu decompose the background
and rain streak layers based on low-rank priors to restore clear images. Luo
_et al._ [5] use sparse coding with high discriminability so that the derained
image layer and the rain layer can be accurately separated. In [13], patch-
based priors are applied for both the clean background and rain layers in the
form of Gaussian mixture models to remove rain streaks.

Recently, trainable CNNs have been proposed to estimate the clean image from a
rainy input directly. In [7], a dual convolutional neural network for
deraining is proposed. The proposed network consists of two parallel branches,
which respectively recovers the structures and details in an end-to-end
manner. Zhang _et al._ [8] propose a density-aware multi-stream densely
connected convolutional neural network that jointly estimates rain density
estimation and derained image. In [9], the authors introduce the Gaussian-
Laplacian image pyramid decomposition technology to the neural network and
propose a lightweight pyramid network for single image deraining. In [3], the
authors propose the two-stage network. Rain streaks, transmissions, and the
atmospheric light are estimated in the first stage. Derained images are
refined in the second stage by the conditional generative adversarial network.
Different from them, our network is implemented in the wavelet space to
capture detailed frequency information.

###  2.2 Attention mechanisms

Attention plays an important role in human perception and computer vision
tasks. Attention mechanisms give feature maps weights so that features of the
sequence of regions or locations are magnified. Generally, there are two
attention mechanisms: spatial attention and channel attention [14, 15]. Mnih
et al. [16] propose an attention model that spatially selects a sequence of
regions to refine feature maps, and the network not only performs well but is
also robust to noisy inputs. Hu et al. [11] propose the squeeze-and-excitation
module and use global average-pooled features to compute channel-wise
attention. Furthermore, Woo et al.[14] combine the spatial and channel
attention to propose a convolutional block attention module. Their module
sequentially infers attention maps along two separate dimensions, channel and
spatial, then attention maps are multiplied to the input feature map for
adaptive feature refinement, which increases the accuracy of image
recognition. In our work, we integrate channel attention and wavelet transform
so that output feature maps contain frequency features but different magnified
weights.

##  3 PROPOSED METHODS

###  3.1 Wavelet Channel Attention module

We first describe the 2-D discrete wavelet transform (DWT) in our model. We
apply Haar wavelets [17] and it contains four kernels,
$LL^{\top},HL^{\top},LH^{\top},HH^{\top}$, where $L$ and $H$ are the low pass
filter and high pass filter respectively. Both filters are

$L=\frac{1}{\sqrt{2}}[1,1]^{\top},H=\frac{1}{\sqrt{2}}[1,-1]^{\top}$ (2)

The DWT for image processing means an image $I$ is convolved and then down-
sampled to obtain the four sub-band images $I_{LL},I_{HL},I_{LH}$ and
$I_{HH}$. The low-pass filter captures smooth surface and texture while the
three high-pass filters extract vertical, horizontal, and diagonal edge-like
information. Even though the downsampling operation is employed, due to the
biorthogonal property of DWT, the original image $I$ can be accurately
reconstructed by the inverse discrete wavelet transform (IDWT). Therefore, the
DWT is seen as four $2\times 2$ convolutional kernels whose weights are fixed
and the stride equals two. Similarily, the IDWT is seen as the transpose
convolution operation whose kernels are identical to DWT’s.

Then, we extend convolutions by combining the DWT and the channel attention
mechanism and propose the wavelet channel attention module (WCAM). A WCAM is a
computational unit that is built upon a transformation $F$ mapping an
intermediate feature map $X\in R^{C\times H\times W}$ to a feature maps $Y\in
R^{4C\times\frac{H}{2}\times\frac{W}{2}}$. Given an intermediate feature map
$X\in R^{C\times H\times W}$, the DWT decomposes $X$ into ${\rm
DWT}(X)=[X_{LL},X_{HL},X_{LH},X_{HH}]\in
R^{4C\times\frac{H}{2}\times\frac{H}{2}}$. $3\times 3$ convolutions and leaky
rectified linear units (LeakyReLUs) are applied to extract various frequency
features from DWT$(X)$ and denoted as $Conv({\rm DWT}(X))$, where $Conv$ is an
operator combining the convolution and the LeakyReLU. Furthermore, we propose
improved feature maps with different frequencies contributing different
weights to restore the clear image. We use the channel attention [11] to
control weights of various channel-wise features. The proposed module
calculates the global average pooling of $X$ and uses $1\times 1$ convolutions
to infer a channel attention map $M_{c}\in R^{4C\times 1\times 1}$. Therefore,
the entire result of WCAM is $Conv({\rm DWT}(X))\odot M_{c}$. The detailed
structure of the WCAM is depicted in Fig. 2(b) and formulated as follows:

$\begin{split}M_{c}&=Conv({\rm AP}(X))\\\ Y&=F(X)=Conv({\rm DWT}(X))\odot
M_{c}\vspace{-0.7cm}\end{split}$ (3)

where AP means global average pooling.

The WCAM reduces the size of feature maps but increases the receptive field to
capture multi-frequency and multi-scale features. To magnify the size of the
feature map, similarly, the inverse wavelet channel attention module (IWCAM)
is proposed. An IWCAM is also a computational unit which can be built upon a
transformation $F^{-1}$ mapping an intermediate feature map $X\in
R^{4C\times\frac{H}{2}\times\frac{W}{2}}$ to a feature map $Y\in R^{C\times
H\times W}$. Given an intermediate feature map $X\in
R^{4C\times\frac{H}{2}\times\frac{W}{2}}$, the IDWT merges $X$ into ${\rm
IDWT}(X)\in R^{C\times H\times W}$ and $3\times 3$ convolutions and
LeaklyReLUs are then adopted and denoted as $Conv({\rm IDWT}(X))\in R^{C\times
H\times W}$. Global average pooling of $X$ is calculated and $1\times 1$
convolutions are used to infer a channel attention map $M_{c}\in R^{4C\times
1\times 1}$. The entire result of IWCAM is $Conv({\rm IDWT}(X))\odot M_{c}$.
The structure of the IWCAM is depicted in Fig. 2(c) and formulated as follows:

$\begin{split}M_{c}&=Conv({\rm AP}(X))\\\ Y&=F^{-1}(X)=Conv({\rm
IDWT}(X))\odot M_{c}\end{split}$ (4)

Fig. 1: The proposed WCAM with a fusion network and components in the network.
(a) The entire network, where green arrows mean skipping connections. (b)
WCAM. (c) IWCAM. (d) The residul WCAM.

###  3.2 Network Architecture

Our network is encoder-decoder structure. The encoder consists of three WCAMs.
Once a WCAM is adopted, the sizes of feature maps become quarter and the
number of channels becomes four times, which not only captures multi-scale
features but various frequency information. At the bottom of the network, the
residual block [18] combining wavelet channel attention is used and shown in
Fig. 1(d). This module aggregates features and makes the learning process
effective, especially in the event of deeper networks. Our decoder consists of
three IWCAMs to generate clear images from extracted features. The terminal
outputs are four confidence maps instead of the restored images. Once
confidence maps for the derived wavelet inputs are predicted, they are
multiplied by the four derived inputs to give the final derained image:

$J={\rm IDWT}(C_{LL}\odot I_{LL},C_{HL}\odot I_{HL},C_{LH}\odot
I_{LH},C_{HH}\odot I_{HH})$ (5)

where $C_{LL}$, $C_{HL}$, $C_{LH}$ ,and $C_{HH}$ are confidence maps for
$I_{LL},I_{HL},I_{LH}$ and $I_{HH}$, respectively. The reason for using the
fusion mechanism is that the low-frequency sub-band plays a role in the
objective quality, while the high-frequency sub-bands can affect the
perceptual quality significantly [19]. When low-frequency parts and high-
frequency parts are optimized separately, the increase of objective quality
cannot decrease perceptual quality. Additionally, like U-net [6], we apply
skipping connections to combine the identical size feature maps from WCAMs and
IWCAMs so that the learning process converges quickly. The entire network
architecture is shown in Fig. 1(a).

##  4 EXPERIMENTAL RESULTS

Fig. 2: Derained results on sample images from the synthetic dataset _Outdoor-
Rain_. (Please zoom-in at screen to view details)  Fig. 3: Derained results on
real rainy images. (Please zoom-in at screen to view details)

###  4.1 Datasets and training details

In this work, _Outdoor-Rain_ dataset [3] is adopted to train and test the
network. This dataset contains clear images, and the corresponding rainy
images generated by Eq.(1). There are 7500 training and validation samples and
1,500 samples for evaluation. Both clean and rainy images are $480\times 720$.
During training, images are cropped to $256\times 256$, the wavelet SSIM loss
[20] and the L1 loss are employed, and RAdam [21] is used as an optimization
algorithm with a mini-batch size of 16. The learning rate starts from 0.0001
and is divided by ten after 100 epochs. The models are trained for 300
iterations. The entire experiments are performed by the Pytorch framework.

###  4.2 Image Deraining results

PSNR and SSIM are chosen as objective metrics for quantitative evaluation. We
select four state-of-the-art works [7, 8, 9, 3] as deep learning-based
benchmarks to make fairly comparisons with our purposed method. For the fair
comparison, all methods are retrained on the same dataset. The comparison
results are shown in Table 1. Table 1 presents our method has the largest PSNR
and SSIM values among all deraining networks, which demonstrates our method
has a superior performance of restoring clean images for this dataset and the
frequency features are beneficial for restoring rainy images. Furthermore, we
perform various methods on synthetic and real rainy photos and results are
depicted in Fig. 2 and Fig. 3. As revealed in Fig. 2 and Fig. 3, since
comparative methods tend to miscalculate rainy concentration, restored images
are dark or have remaining rain and mist. In contrast, our purposed method
predicts better-derained results with balanced colors and detailed edges.

We analyze how the WCAM and the fusion help to refine derained results with
three experiments. The first experiment uses convolutions and wavelets without
channel attention and fusion. The second experiment uses the proposed modules
without fusion. The third experiment estimates feature maps to fuse sub-band
images without channel attention. Table 2 compares our method against three
baselines and demonstrates wavelet channel attention and fusion contribute the
best results.

Table 1: Quantitative SSIM and PSNR on the synthetic _Outdoor-Rain_ dataset. | [7] | [8] | [9] | [3] | Ours
---|---|---|---|---|---
PSNR | 17.92 | 21.64 | 18.16 | 21.17 | 24.89
SSIM | 0.676 | 0.788 | 0.723 | 0.742 | 0.813
Table 2: Ablation study shows how the WCAM and the fusion help to refine derained results. | PSNR | SSIM
---|---|---
Ours, w/o attention, w/o fusion | 21.44 | 0.784
Ours, w/o fusion | 23.34 | 0.764
Ours, w/o attention | 23.24 | 0.810
Ours | 24.89 | 0.813

##  5 CONCLUSION

In this paper, the wavelet channel attention module with a fusion network is
proposed for single image deraining. The wavelet transform and the inverse
wavelet transform are substituted for down-sampling and up-sampling to extract
various frequency features. The channel attention effectively controls ratios
of feature maps. Furthermore, the proposed network estimates confidence maps
for each derived wavelet input. Confidence maps and derived inputs are fused
to render final derained results. Experiments on synthetic and real images
verify the superiority of our model compared to the state-of-the-art results.

## References

  * [1] Siyuan Li, Iago Breno Araujo, Wenqi Ren, Zhangyang Wang, Eric K. Tokuda, Roberto Hirata Junior, Roberto Cesar-Junior, Jiawan Zhang, Xiaojie Guo, and Xiaochun Cao,  “Single image deraining: A comprehensive benchmark analysis,”  in The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [2] Earl J McCartney,  “Optics of the atmosphere: scattering by molecules and particles,”  New York, John Wiley and Sons, Inc., 1976. 421 p., 1976.
  * [3] Ruoteng Li, Loong-Fah Cheong, and Robby T Tan,  “Heavy rain image restoration: Integrating physics model and conditional adversarial learning,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2019, pp. 1633–1642.
  * [4] Yi-Lei Chen and Chiou-Ting Hsu,  “A generalized low-rank appearance model for spatio-temporally correlated rain streaks,”  in Proceedings of the IEEE International Conference on Computer Vision, 2013, pp. 1968–1975.
  * [5] Yu Luo, Yong Xu, and Hui Ji,  “Removing rain from a single image via discriminative sparse coding,”  in Proceedings of the IEEE International Conference on Computer Vision, 2015, pp. 3397–3405.
  * [6] Olaf Ronneberger, Philipp Fischer, and Thomas Brox,  “U-net: Convolutional networks for biomedical image segmentation,”  in International Conference on Medical image computing and computer-assisted intervention. Springer, 2015, pp. 234–241.
  * [7] Jinshan Pan, Sifei Liu, Deqing Sun, Jiawei Zhang, Yang Liu, Jimmy Ren, Zechao Li, Jinhui Tang, Huchuan Lu, Yu-Wing Tai, et al.,  “Learning dual convolutional neural networks for low-level vision,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2018, pp. 3070–3079.
  * [8] He Zhang and Vishal M. Patel,  “Density-aware single image de-raining using a multi-stream dense network,”  in The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2018.
  * [9] Xueyang Fu, Borong Liang, Yue Huang, Xinghao Ding, and John Paisley,  “Lightweight pyramid networks for image deraining,”  IEEE transactions on neural networks and learning systems, 2019\.
  * [10] Hao-Hsiang Yang and Yanwei Fu,  “Wavelet u-net and the chromatic adaptation transform for single image dehazing,”  in IEEE International Conference on Image Processing (ICIP). IEEE, 2019, pp. 2736–2740.
  * [11] Jie Hu, Li Shen, and Gang Sun,  “Squeeze-and-excitation networks,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2018, pp. 7132–7141.
  * [12] Wenqi Ren, Lin Ma, Jiawei Zhang, Jinshan Pan, Xiaochun Cao, Wei Liu, and Ming-Hsuan Yang,  “Gated fusion network for single image dehazing,”  in The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2018.
  * [13] Yu Li, Robby T Tan, Xiaojie Guo, et al.,  “Rain streak removal using layer priors,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 2736–2744.
  * [14] Sanghyun Woo, Jongchan Park, Joon-Young Lee, and In So Kweon,  “Cbam: Convolutional block attention module,”  in Proceedings of the European Conference on Computer Vision (ECCV), 2018, pp. 3–19.
  * [15] Chao-Han Yang, Jun Qi, Pin-Yu Chen, et al.,  “Characterizing speech adversarial examples using self-attention u-net enhancement,”  in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2020, pp. 3107–3111.
  * [16] Volodymyr Mnih, Nicolas Heess, et al.,  “Recurrent models of visual attention,”  in Advances in neural information processing systems, 2014, pp. 2204–2212.
  * [17] Stéphane Mallat,  A wavelet tour of signal processing,  Elsevier, 1999.
  * [18] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  “Deep residual learning for image recognition,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [19] Xin Deng, Ren Yang, Mai Xu, and Pier Luigi Dragotti,  “Wavelet domain style transfer for an effective perception-distortion tradeoff in single image super-resolution,”  in The IEEE International Conference on Computer Vision (ICCV), October 2019.
  * [20] Hao-Hsiang Yang, Chao-Han Huck Yang, and Yi-Chang James Tsai,  “Y-net: Multi-scale feature aggregation network with wavelet structure similarity loss function for single image dehazing,”  in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2020, pp. 2628–2632.
  * [21] Liyuan Liu, Haoming Jiang, Pengcheng He, Weizhu Chen, Xiaodong Liu, Jianfeng Gao, and Jiawei Han,  “On the variance of the adaptive learning rate and beyond,”  in Proceedings of the Eighth International Conference on Learning Representations (ICLR 2020), April 2020.


#  Footstep Planning with Encoded Linear Temporal Logic Specifications

Vikram Ramanathan

###### Abstract

This article presents an approach to encode Linear Temporal Logic (LTL)
Specifications into a Mixed Integer Quadratically Constrained Quadratic
Program (MIQCQP) footstep planner. We propose that the integration of LTL
specifications into the planner not only facilitates safe and desirable
locomotion between obstacle-free regions, but also provides a rich language
for high-level reasoning in contact planning. Simulations of the footstep
planner in a 2D environment satisfying encoded LTL specifications demonstrate
the results of this research.

##  I Introduction

Humanoid locomotion is accomplished by changing its footholds. In order to
move feasibly and efficiently, a footstep planner must be able to find
sequences of foot positions and orientations that will realize the robot’s
locomotion objectives. One established solution for the footstep planning
problem is to perform some discrete sampling-based search algorithm over
reachable foot positions and orientations. This approach first determines
candidate footsteps through one of several methods, including:

  * •

intersecting valid sample limb configurations with the environment to
determine which limb configurations are close to contact and then projecting
these configurations onto contact using inverse kinematics [1]

  * •

assuming general bounds on the displacement of each foothold based on the
robot structure and kinematic reachability [2] [3] [4] [5]

  * •

defining an action-effect mapping [6]

Then, a search algorithm, such as A* [7] [6] [8] [4] or Rapidly-exploring
Random Trees (RRT) [5] [9] is used to find a sequence of footsteps. However,
the performance of these search algorithms largely depends on the quality of
the heuristic function that guides the search as well as the granularity of
the discretization. As an alternative, Deits and Tedrake [10] proposed a
footstep planner formulated as a mixed integer optimization program with
continuous decision variables for footstep position and orientation, and
integer decision variables to eliminate non-convex constraints. After
employing the IRIS algorithm [11] to compute large obstacle-free convex
regions within the provided environment, their solution solves a Mixed Integer
Quadratically Constrained Quadratic Program (MIQCQP) with a clever
discretization of cos and sin in its reachability constraint. This approach
solves for globally optimal footstep plans at realtime rates while satisfying
reachability constraints and avoiding obstacles.

##  II Problem Description

In this article, we consider the problem of encoding temporal logic
specifications into a MIQCQP footstep planner to improve and extend its
functionality. Continuous footstep planners like the Mixed Integer
Quadratically Constrained Quadratic Program (MIQCQP) footstep planner proposed
by Deits and Tedrake [10] adopt an optimization framework, thereby avoiding
the use of search heuristics as well as the sampling of discretized footstep
configurations. The problem of integrating temporal logic specifications into
a sampling-based search algorithm such as Probabilistic Road Maps (PRMs) and
Randomly-Exploring Random Trees (RRTs) has been studied [12] [13]. However, we
have adopted the continuous footstep planner to take advantage of its ability
to solve for optimal solutions as well as the ease with which its optimization
framework can be expanded upon to include simplified robot dynamics
constraints and compute for reaction forces.

##  III Solution

As mentioned previously, this work builds upon Deits’ [10] footstep planner.
This planner solves an MIQCQP to determine x, y and $\theta$ for N footsteps.
We first will cover the formulation of this footstep planner before discussing
the encoding of LTL specifications.

###  III-A MIQCQP Footstep Planner

The footstep planner seeks to minimize a quadratic cost function (in terms of
the length of each stride and the distance between the terminal footstep and
the goal footstep position and orientation) subject to obstacle avoidance and
reachability constraints. Obstacle avoidance is facilitated by first
decomposing the non-convex obstacle-free configuration space into a set of
convex “safe” regions [11]. These convex regions are then defined as polygons
by the halfspace respresentation and added as convex constraints to the
program. A binary variable is assigned to each possible footstep and region
pair. If such a binary variable is 1, the corresponding footstep is assigned
to the region. More precisely, the matrix, $H\in\\{0,1\\}^{R\times N}$, where
$R$ is the number of regions and $N$ is the number of footsteps, is
constructed such that if $H_{r,j}$ = 1, then the $j^{th}$ footstep,
$j=\\{1,...,N\\}$, is assigned to the region $r$. This implication is defined
as a mixed-integer linear constraint using the big-M formulation [14],

$-M(1-H_{r,j})+A_{r}f_{j}\leq b_{r}$ (1)

where

$f_{j}=\begin{bmatrix}x_{j}\\\ y_{j}\\\ \theta_{j}\end{bmatrix}$

, $M$ is a sufficiently large positive number, $A_{r}$ defines the normal
vectors corresponding to the supporting halfplanes of the polygonal obstacle-
free regions, and $b_{r}$ defines the offsets between the halfplanes and the
origin. However, since the footstep must only belong to one region instead of
the intersection of all these regions, $\sum_{r=1}^{R}H_{r,j}=1$ must be
enforced.

The footstep reachability constraint is defined as the intersection of two
circular regions offset from the previous foothold (as shown in Eq. 2 below).
These offsets are mirrored for left and right feet.

$\left\lVert\begin{bmatrix}x_{j}\\\
y_{j}\end{bmatrix}-\bigg{(}\begin{bmatrix}x_{j-1}\\\
y_{j-1}\end{bmatrix}+\begin{bmatrix}cos(\theta_{j})&-sin(\theta_{j})\\\
sin(\theta_{j})&cos(\theta_{j})\end{bmatrix}p_{i}\bigg{)}\right\rVert\leq
r_{i}$ (2)

where $i\in\\{1,2\\}$, $j\in\\{2,...,N\\}$, $p_{i}$ are the centers of the
circles, and $r_{i}$ are the radii of the circles. By tuning $p_{i}$ and
$r_{i}$, we can obtain a conservative approximation of a biped’s footstep
reachability. In order to eliminate the nonlinearity, $cos$ and $sin$ are
approximated by linear piece-wise functions along with binary variables to
decide which linear segment of the function to use based on $\theta$’s value.
These binary variables are defined as binary matrices,
$S,C\in\\{0,1\\}^{L\times N}$, where L is the number of piece-wise segments
and N is the number of footsteps. Since $\theta$ cannot belong to multiple
segments, $\sum_{l=1}^{L}S_{l,j}=\sum_{l=1}^{L}C_{l,j}=1$ must be satisfied.
Our piece-wise approximation of $cos$ and $sin$ defines 5 linear segments:

$sin(\theta)=\begin{cases}-\theta-\pi\hskip 15.00002pt-\pi\leq\theta<1-\pi\\\
-1\hskip 38.00008pt1\leq\theta<-1\\\ \theta\hskip 45.00006pt-1\leq\theta<1\\\
1\hskip 45.00006pt1\leq\theta<\pi-1\\\ -\theta+\pi\hskip
20.00003pt\pi-1\leq\theta<\pi\\\ \end{cases}$ (3)

$cos(\theta)=\begin{cases}-1\hskip 43.00009pt-\pi\leq\theta<-\pi/2-1\\\
\theta+\pi/2\hskip 23.00006pt-\pi/2-1\leq\theta<1-\pi/2\\\ 1\hskip
53.0001pt1-\pi/2\leq\theta<\pi/2-1\\\ -\theta+\pi/2\hskip
20.00003pt\pi/2-1\leq\theta<\pi/2+1\\\ -1\hskip
47.50006pt\pi/2+1\leq\theta<\pi\\\ \end{cases}$ (4)

(a) $sin(\theta)$

(b) $cos(\theta)$

Figure 1: Piece-wise Approximations of $cos(\theta)$ and $sin(\theta)$

The graphical representation of this approximation is shown in Fig. 1. We also
bound the change in footstep orientation in every step to $\pi/8$ rad.

Putting the quadratic cost function and all these constraints together, we get
the following MIQCQP:

$\displaystyle\underset{f_{1},...,f_{j},S,C,H}{\text{minimize}}$
$\displaystyle(f_{N}-g)^{T}Q(f_{N}-g)+\sum_{j=1}^{N-1}(f_{j+1}-f_{j})^{T}R(f_{j+1}-f_{j})$
(5) subject to $\displaystyle H_{r,j}\implies A_{r}f_{j}\leq b_{r}$
$\displaystyle
S_{l,j}\implies\begin{cases}\phi_{l}\leq\theta_{j}\leq\phi_{l+1}\\\
s_{j}=g_{l}\theta_{j}+h_{l}\\\ \end{cases}$ $\displaystyle
C_{l,j}\implies\begin{cases}\phi_{l}\leq\theta_{j}\leq\phi_{l+1}\\\
c_{j}=g_{l}\theta_{j}+h_{l}\\\ \end{cases}$
$\displaystyle\left\lVert\begin{bmatrix}x_{j}\\\
y_{j}\end{bmatrix}-\bigg{(}\begin{bmatrix}x_{j-1}\\\
y_{j-1}\end{bmatrix}+\begin{bmatrix}c_{j}&-s_{j}\\\
s_{j}&c_{j}\end{bmatrix}p_{i}\bigg{)}\right\rVert\leq r_{i}$
$\displaystyle\sum_{r=1}^{R}H_{r,j}=\sum_{l=1}^{L}S_{l,j}=\sum_{l=1}^{L}C_{l,j}=1$
$\displaystyle\theta_{j}-\theta_{j-1}\leq\pi/8$

where $r=1...R$, , $j=\\{1,...,N\\}$, and $l=1,...,L$.

###  III-B Linear Temporal Logic

####  III-B1 Preliminaries

In this work, the underlying time domain is discrete since each “execution” or
footstep corresponds to the advancement of a single time-unit. We consider
Linear Temporal Logic (LTL) to describe desired planning behavior. LTL
formulae are defined over a set of atomic propositions, $AP$, according to the
following grammar [15]:

$\phi::=\text{true
}|\,a\,|\,\phi_{1}\wedge\phi_{2}\,|\,\lnot\phi\,|\,\bigcirc\phi\,|\,\phi_{1}U\phi_{2}$

where $a\in AP$, $\bigcirc$ and $U$ are the temporal operators for “next”
(signifies that the proposition is true in the next step) and “until”
(signifies that the proposition preceding the operator is true until the
future moment when the proposition succeeding the operator is true), and
$\lnot$ and $\wedge$ are the negation and conjuction Boolean operators
respectively. The other Boolean operators for disjunction, $\vee$, implication
$\implies$, and equivalence $\iff$ can be derived from the previously defined
grammar:

$\phi_{1}\vee\phi_{2}=\lnot(\lnot\phi_{1}\wedge\lnot\phi_{2})$ (6)
$\phi_{1}\implies\phi_{2}=\lnot\phi_{1}\vee\phi_{2}$ (7)
$\phi_{1}\iff\phi_{2}=(\phi_{1}\implies\phi_{2})\wedge(\phi_{2}\implies\phi_{1})$
(8)

In addition, the temporal operators, $\lozenge$ (eventually) and $\square$
(always) are defined by:

$\lozenge\phi=true\,U\,\phi$ (9) $\square\phi=\lnot\lozenge\lnot\phi$ (10)

Let $\sigma$ be an infinite sequence of states,
$\sigma=A_{0},A_{1},...(2^{AP})^{\omega}$. Let $\sigma_{i}$ denote the word
$\sigma$ from position $i$. The semantics of an LTL formula is defined as a
language $Words(\phi)$ that contains all infinite words over the alphabet
$2^{AP}$ that satisfy $\phi$. Every LTL formula is associated with a single
linear time property. The satisfaction relations are defined inductively by:

$\sigma_{i}\models a\;\text{iff}\;a\in A_{0}$ (11)
$\sigma_{i}\models\phi_{1}\wedge\phi_{2}\;\text{iff}\;\sigma_{i}\models\phi_{1}\;\text{and}\;\sigma_{i}\models\phi_{2}$
(12) $\sigma_{i}\models\lnot\phi\;\text{iff}\;\sigma_{i}\nvDash\phi$ (13)
$\sigma_{i}\models\bigcirc\phi\;\text{iff}\;\sigma_{i+1}\models\phi$ (14)
$\sigma_{i}\models\phi_{1}\;U\;\phi_{2}\;\text{iff}\;\exists j\geq
i\;\text{s.t.}\;\sigma_{j}\models\phi_{2}\;\text{and}\;\sigma_{k}\models\phi_{1}\\\
\forall i\leq k<j$ (15)

####  III-B2 Mixed Integer Encoding

In this section, we will detail the mixed integer encoding of Linear Temporal
Logic specifications. Let $p$ be an atomic proposition defined over some
footstep. Then, this proposition has a corresponding binary variable $P^{k}$,
which is defined to be 1 when $p$ is true at the $k^{th}$ footstep and 0 when
$p$ is false at the $k^{th}$ footstep (where $k=0...N-1$).

To appropriately extend the semantics of Linear Temporal Logic to mixed
integer programming, we redefine the original satisfaction relations outlined
in the previous section in terms of integer variables. Let $\sigma_{k}$ is a
finite run that starts at footstep $k$.

Firstly, $\sigma_{k}\models p$ can be encoded by the following integer
constraint,

$P^{k}=1$ (16)

The conjunction satisfiability relation,
$\sigma_{k}\models\bigwedge_{j=1}^{m}\phi_{j}$, where m is any integer greater
than 1, can be encoded by the following integer constraint,

$\sum_{j=1}^{m}P^{k}_{\phi_{j}}=m$ (17)

The negation satisfiability relation, $\sigma_{k}\models\lnot\phi$ can be
encoded by the following integer constraint,

$P^{k}_{\phi}=0$ (18)

Although the disjunction satisfiability relation,
$\sigma_{k}\models\bigvee_{j=1}^{m}\phi_{j}$, where m is any positive integer
greater than 1, can easily be derived from the conjunction (Eq. 17) and
negation (Eq. 18) satisfiability relations by means of Eq. 6, we propose a
more succinct encoding:

$\sum_{j=1}^{m}P^{k}_{\phi_{j}}>=1$ (19)

Finally, $\sigma_{k}\models\bigcirc\phi$ can be encoded by the following
integer constraint,

$P^{k+1}_{\phi}=1$ (20)

The encoding of the until operator, $U$ can be defined in terms of Eqs. 11-20
above. Specifically, the expansion laws of LTL formulae [15] give us the
following recursive identity for the until operator, $U$:

$\phi_{1}U\phi_{2}=\phi_{2}\vee(\phi_{1}\wedge\bigcirc(\phi_{1}U\phi_{2}))$
(21)

Granted that footstep planning involves a finite horizon, any loops in the
plan are finite and would not lead to circular reasoning (as pointed out in
[16]). As a result, we can avoid using an auxiliary encoding of $U$, which
involves under-approximating its functionality [17]. Instead, our proposed
approach formulates a set of nested mixed-integer constraints (using Eqs.
11-20) over a vector of binary variables, in which each variable corresponds
to the satisfiability of $\phi_{1}U\phi_{2}$ at a particular foothold in the
run.

Let $T\in\\{0,1\\}^{1\times(N-k)}$ be a vector of binary variables of length
$N-k+1$. Assuming Eqs. 11-20, $\sigma_{k}\models\phi_{1}U\phi_{2}$ can be
encoded by the following set of integer constraints:

$T^{i}=P_{\phi_{2}}^{i}\vee(P_{\phi_{1}}^{i}\wedge\bigcirc T^{i+1})\hskip
20.00003pti=k,...,N-1$ (22)

with base case,

$T^{N}=P_{\phi_{2}}^{N}$ (23)

and satisfiability constraint (i.e. $\sigma_{k}$ must satisfy the LTL
specification, $\phi_{1}U\phi_{2}$)

$T^{k}=1$ (24)

Instead of deriving the mixed-integer encodings for safety, liveness, and
persistence LTL formulae from the general definitions above, we can more
efficiently encode them as follows:

Safety (“always”): The safety LTL specification,
$\sigma_{k}\models\square\phi$ can be encoded as the following integer
constraint:

$\sum_{i=k}^{N}P_{\phi}^{i}=N-k+1$ (25)

Liveness (“eventually”): The liveness LTL specification,
$\sigma_{k}\models\lozenge\phi$ can be encoded as the following integer
constraint:

$\sum_{i=k}^{N}P_{\phi}^{i}\geq 1$ (26)

Liveness (“repeated eventually”): The liveness LTL specification,
$\sigma_{k}\models\square\lozenge\phi$ (repeated eventually) can be encoded as
the following set of integer constraints:

$\sum_{i=j}^{N}P_{\phi}^{i}\geq 1\hskip 20.00003ptj=k,...,N$ (27)

Persistence (“eventually always”): The persistence LTL specification,
$\sigma_{k}\models\lozenge\square\phi$ (eventually always) is defined with the
help of Eq. 19, and is encoded as follows:

$\bigvee_{j=k}^{N}\sum_{i=j}^{N}P_{\phi}^{i}=N-j+1$ (28)

##  IV Examples

Our mixed integer encodings facilitate the integration of LTL specifications
in the footstep planner. This is particularly powerful in ensuring safe
performance of region-based locomotion tasks. Specifications can detail
desired locomotion behavior such as the requirement that the robot must
eventually enter a set of regions or that it must reach some region within
some number of footsteps. Specifications can also be used to outline
underlying environment (region-specific) characteristics such as the
requirement that if a footstep is placed in a particular region, the following
footsteps must remain in that region or that access to certain regions is only
granted provided particular regions have already been visited.

In this section, we present three LTL-constrained footstep planning scenarios.
These scenarios were solved using the commercial optimization solver, Gurobi
on a laptop with a 2.3 GHz quad-core processor and 16 GB of memory. Source
code and animations of these scenarios can be found on Github at
https://github.com/vikr53/ltl_footstep_planner .

###  IV-A Liveness

We will first consider encoding a region-based liveness specification that
necessitates that certain obstacle-free regions be visited eventually. The
planner then solves for an optimal set of footsteps that fulfil this
requirement while also reaching some goal footstep position and orientation.

Let $p_{R_{3}}$ and $p_{R_{4}}$ be the atomic propositions, “a footstep is in
Region 1” and “a footstep is in Region 2” respectively.

LTL Formula:

$\lozenge(p_{R_{3}}\vee p_{R_{4}})$ (29)

Mixed-Integer Encoding:  Using the liveness (Eq. 26) and disjunction (Eq. 19)
encodings, we can express the above LTL specification as the following linear
mixed-integer constraint:

$\sum^{N}_{j=1}(H_{R_{3},j}+H_{R_{4},j})\geq 1$ (30)

Figure 2: Scenario 1.  Encoded LTL Specification: $\lozenge(p_{R_{3}}\vee
p_{R_{4}})$. Right-leg footsteps are shown in red stars and left-leg footsteps
are shown in blue circles. The orientation of each footstep is depicted by its
protruding black arrow. Each footstep is numbered in brown. Goal footstep
location and orientation: [$x$: 0, $y$: 1.5, $\theta$: $\pi/2$]

The result from the planner is depicted in Figure 2, where the red stars are
the right leg footsteps, the blue circles are the left leg footsteps and the
green boxes are the obstacle-free convex regions. It also must be noted that
the first two footholds are assumed to be fixed (i.e. the humanoid robot is in
some initial double-support stance). In this scenario, the planner must
eventually reach either region 3 (“R3”) or region 4 (“R4”) before reaching its
goal state - ($x$: 0, $y$: 1.5, $\theta$: $\pi/2$). As is evident from the
simulation, the planner is able to reason that the most efficient way to
satisfy this requirement is to visit region 4 through region 5 before
returning to its goal state in region 5. This optimal result is solved for in
around 25s.

###  IV-B Until

We now demonstrate the encoding of the until, $U$, temporal operator in the
footstep planner. The until operator can be used to ensure that the approach
to some goal state remains within certain pre-defined desirable regions. This
has applications in a wide-variety of locomotive tasks such as surveying and
reconnaissance.

Let $p_{R_{1}}$, $p_{R_{2}}$, and $p_{R_{3}}$ be the atomic propositions, “a
footstep is in Region 1”, “a footstep is in Region 2” and “a footstep is in
Region 3” respectively. The length of the plan is $N=13$.

LTL Formula:

$(p_{R_{1}}\vee p_{R_{2}})Up_{R_{3}}$ (31)

Mixed-Integer Encoding:  First, we define a vector of binary variables,
$T\in\\{0,1\\}^{1\times N}$. Then, we specify the base case (as given by Eq.
23):

$T^{N}=H_{R_{3},N}$ (32)

where R3 denotes the index for Region 3 in matrix $H$. Next, following the
formulation defined in Eq. 22 and Eq. 19, we first encode $\phi_{1}^{k}\equiv
p_{R_{1}}^{k}\vee p_{R_{2}}^{k}$, where $k\in{1,...,N-1}$, and $p_{R_{1}}^{k}$
and $p_{R_{2}}^{k}$ are the atomic propositions, “the $k^{th}$ footstep is in
Region 1”, “the $k^{th}$ footstep is in Region 2”.

$H_{R_{1},k}+H_{R_{2},k}\geq 1-M(1-P_{\phi_{1}}^{k})$ (33)
$H_{R_{1},k}+H_{R_{2},k}\leq 1-m+M(P_{\phi_{1}}^{k})$ (34)

where $m$ and $M$ are sufficiently small and large positive integers
respectively. Then, using the big-M formulation and the conjuction encoding
given by Eq. 17, we encode the term,
$\phi_{1}\wedge\bigcirc(\phi_{1}U\phi_{2})$:

$P_{\phi_{1}}^{k}+T^{k+1}\geq 2-M(1-B^{k})$ (35) $P_{\phi_{1}}^{k}+T^{k+1}\leq
2-m+M(B^{k})$ (36)

where $B$ is a vector of binary variables, i.e.
$B\in{0,1}^{1\times\\{N-1\\}}$. Finally, we have the following two constraints
to finish encoding the entire specification:

$H_{R_{3},k}+B^{k}\geq 1-M(1-T^{k})$ (37) $H_{R_{3},k}+B^{k}\leq 1-m+M(T^{k})$
(38)

All that remains is simply adding the satisfiability constraint (Eq. 24):

$T^{1}=1$ (39)

Figure 3: Scenario 2.  Encoded LTL Specification: $(p_{R_{1}}\vee
p_{R_{2}})Up_{R_{3}}$. Right-leg footsteps are shown in red stars and left-leg
footsteps are shown in blue circles. The orientation of each footstep is
depicted by its protruding black arrow. Each footstep is numbered in brown.
Goal footstep location and orientation: [$x$: 2, $y$: 1.5, $\theta$: $\pi/2$]

Fig. 3 displays the solved LTL-constrained footstep plan. It is evident that
the 13 planned footsteps abide by the given specification, “until Region 3 is
reached, footsteps must either be in Region 1 or Region 2”, while reaching its
goal state - ($x$: 2, $y$: 1.5, $\theta$: $\pi/2$). This example was solved in
around 2s.

###  IV-C Timed Specifications

In this example, we consider the encoding of timed LTL constraints, where each
footstep is assumed to be a single time unit. Defined over a discrete time
domain, these specifications are relatively simple to encode but allow for the
expression of interesting constraints like “always be in a particular region
or set of regions for the first $p$ steps” or “eventually reach a region
within the $p^{th}$ and $q^{th}$ footstep”, where $p,q$ are arbitrary positive
integers. For the sake of this example, we consider the following
specification:

LTL Formula:

$\square^{7\leq k\leq 15}p_{R_{2}}$ (40)

Mixed-Integer Encoding:  Using Eq. 25, we can write the above specification as
the following constraint:

$\sum_{j=7}^{15}H_{R2,j}\geq 15-7+1$ (41)

Figure 4: Scenario 3.  Encoded LTL Specification: $\square^{7\leq k\leq
15}p_{R_{2}}$. Right-leg footsteps are shown in red stars and left-leg
footsteps are shown in blue circles. The orientation of each footstep is
depicted by its protruding black arrow. Each footstep is numbered in brown.
Goal footstep location and orientation: [$x$: 1.5, $y$: 2.2, $\theta$:
$3\pi/4$]

The footstep plan generated for this scenario is shown in Fig. 4. The plan’s
length is $N=18$ steps. It is evident that the plan satisfies the discrete-
timed specification, “the $7^{th}$ to the $15^{th}$ footstep must be in Region
2”, while finally reaching its goal state - ($x$: 1.5, $y$: 2.2, $\theta$:
$3\pi/4$). The solve time for this example was around 2s.

Aside from encoding surveying behaviors, an interesting application of such
timed specifications is assessing the accessibility of a region. While taking
into account the particular formulation of the footstep planner (the planning
horizon, reachability constraints, the obstacle-free regions and other encoded
specifications), this specification can determine whether a specific region or
set of regions can be reached within some number of steps. For example, in
order to determine whether Region 2 can be reached within 5 footsteps for the
formulation used for this example, we can simply encode the specification,
$\lozenge^{\leq 5}p_{R_{2}}$, into the optimization program. By nature of the
encoding, the planner attempts to satisfy this specification. Thus, if the
MIQCQP is rendered infeasible, we know definitively that it is not possible to
reach Regon 2 in 5 footsteps, i.e. Region 2 is not accessible in 5 steps with
the current formulation. The opposite is true when the MIQCQP is solvable,
i.e. Region 2 is accessible in 5 steps.

##  V Ongoing/Future Work

###  V-A Stride Adjustment

We have also considered the encoding of region-based stride adjustments. Since
the traversibility of terrain can vary significantly over a footstep plan, the
planner may need to tread certain regions more carefully than others. Often,
the first step in reducing the risk of loosing dynamic stability on account of
the terrain is reducing the stride length. This can be easily implemented by
changing the reachability constraints for regions with difficult terrain:

From Eq. 2, $p_{i}$ and $r_{i}$, $i\in\\{1,2\\}$, are the centers and radii of
the circles whose intersection defines the footstep reachability of the robot.
To define the reduced reachability constraint, we merely have to intersect two
relatively smaller and less offset circles with centers and radii, $p_{i}^{s}$
and $r_{i}^{s}$ $i\in\\{1,2\\}$ respectively. Then, we use the big-M
formulation to encode the implication that smaller strides only occur in
specific regions. For example, smaller strides in Region 2 can be encoded as:

$\left\lVert\begin{bmatrix}x_{j}\\\
y_{j}\end{bmatrix}-\bigg{(}\begin{bmatrix}x_{j-1}\\\
y_{j-1}\end{bmatrix}+\begin{bmatrix}cos(\theta_{j})&-sin(\theta_{j})\\\
sin(\theta_{j})&cos(\theta_{j})\end{bmatrix}p_{i}^{s}\bigg{)}\right\rVert-M(1-H_{R_{2},j})\\\
\leq r_{i}^{s}$ (42)

Figure 5: Stride Adjustment Example. Planner makes shorter strides in Region 2
since it is assumed to have difficult terrain. Goal footstep location and
orientation: [$x$: 3.5, $y$: 0.5, $\theta$: 0]

Fig. 5 depicts the result of encoding the above specification in the footstep
planner. It is evident that the planner is able to shorten its stride when the
robot enters Region 2. As mentioned previously, this has several applications,
especially in dealing with difficult terrain. For example, given that a
certain region is known to have more slippery terrain, the robot can shorten
its stride in this region to reduce the risk of slipping. However, it must be
noted that solely planning shorter strides does not guarantee safer
locomotion. Appropriate contact forces must also be planned. We hope to
investigate the integration of our LTL-constrained planner with trajectory
optimization (for planning through contact) [18].

###  V-B Multi-Contact Planning

In addition to specifying region-based locomotion behavior, encoded LTL
specifications can potentially be extended to the more difficult multi-contact
problem, allowing us to encode higher-level reasoning about the safety and
performance of these plans. We propose that this is possible by first defining
the planner in three dimensions, then approximating the COM position as the
center of the support polygon and adding approximate arm reachability
constraints. Lastly, LTL constraints can be defined over the ordering of
contacts. Given the combinatorial complexity of the multi-contact problem,
this formulation would be able to not only make the problem more tractable but
also guarantee the safety and performance of multi-contact plans. We are
actively researching this idea and have done some preliminary tests on the
simpler footstep planner (which has only 2 end-effectors).

Preliminary Results:  To demonstrate that the ordering of end-effectors in
contact need not be defined explicitly (as was done in the formulation of the
original MIQCQP footstep planner [10]), we encode ordering-related LTL
specifications into the planner and let the planner formulate the appropriate
end-effector contact sequence. In the case of the footstep planner, the safe
and desirable ordering of end-effectors in contact is trivial - the right foot
must follow the left foot, and a footstep can be made by either the left or
right foot, not both. This is encoded as the linear mixed-integer constraints:

$LL^{j}+RL^{j}=1$ (43)

$-M(1-LL^{j-1})+RL^{j}\leq 1$ (44)

$M(1-LL^{j-1})+RL^{j}\geq 1$ (45)

where $LL^{j}$ and $RL^{j}$, for $j\in\\{2,...,N\\}$, are binary variables
that determine whether the left leg and right leg make the $j^{th}$ footstep
respectively. In addition, the reachability constraint in Eq. 2 must be
changed to adjust to whether a particular footstep was chosen to be made with
the left or right foot:

$\left\lVert\begin{bmatrix}x_{j}\\\
y_{j}\end{bmatrix}-\bigg{(}\begin{bmatrix}x_{j-1}\\\
y_{j-1}\end{bmatrix}+\begin{bmatrix}cos(\theta_{j})&-sin(\theta_{j})\\\
sin(\theta_{j})&cos(\theta_{j})\end{bmatrix}p_{1}\bigg{)}\right\rVert-M(1-LL^{j})\\\
\leq r_{1}$ (46)

$\left\lVert\begin{bmatrix}x_{j}\\\
y_{j}\end{bmatrix}-\bigg{(}\begin{bmatrix}x_{j-1}\\\
y_{j-1}\end{bmatrix}+\begin{bmatrix}cos(\theta_{j})&-sin(\theta_{j})\\\
sin(\theta_{j})&cos(\theta_{j})\end{bmatrix}p_{2}\bigg{)}\right\rVert-M(1-RL^{j})\\\
\leq r_{2}$ (47)

This encoding was successful and generated similar results to the original
formulation. It must be noted that while these specifications are in effect
equivalent to original implicit encoding for the case of the footstep planner,
it is neither preferable nor feasible to add ordering constraints on a
contact-by-contact basis for the more complicated multi-contact problem.
Instead, encoding such LTL specifications allow the planner to not only
determine the optimal sequence of contacts but also guarantee that this
sequence follows some safety constraints. In the multi-contact planning
problem, these constraints could include simple requirements like
$\square(p_{lleg}\vee p_{rleg})$ (one foot should always be on the ground),
and even criteria for the necessity of multi-contact plans (such as
uncertainty in the stability of certain contacts). This idea seems promising,
and we hope to make progress on its development soon. The code for this
example is available as the script, footstep_planner_contact_ordering_specs.py
in the Github repository linked earlier in this paper.

###  V-C Exploiting Sparse Constraint Matrices

While this planner can solve our footstep planning problem for 10 steps, 5
regions and a few LTL constraints in less than a second, increasing the
planning horizon (i.e. planning for more footsteps), introducing new regions
and encoding additional constraints increase the complexity of solving the
mixed-integer program. Hence, a limitation of this planer is its poor
scalability.

We propose that the sparsity of the $H$, $S$ and $C$ matrices can possibly be
exploited to increase the solve times of plans with large horizons (greater
than 20 footsteps). However, this is merely an idea at the moment, and we hope
to explore this in more detail soon.

## References

  * [1] S. Tonneau, A. Del Prete, J. Pettré, C. Park, D. Manocha, and N. Mansard, “An efficient acyclic contact planner for multiped robots,” _IEEE Transactions on Robotics_ , vol. 34, no. 3, pp. 586–601, 2018.
  * [2] A. Hornung, A. Dornbush, M. Likhachev, and M. Bennewitz, “Anytime search-based footstep planning with suboptimality bounds,” in _2012 12th IEEE-RAS International Conference on Humanoid Robots (Humanoids 2012)_. IEEE, pp. 674–679.
  * [3] P. Karkowski, S. Oßwald, and M. Bennewitz, “Real-time footstep planning in 3d environments,” in _2016 IEEE-RAS 16th International Conference on Humanoid Robots (Humanoids)_. IEEE, 2016, pp. 69–74.
  * [4] J. Kuffner, K. Nishiwaki, S. Kagami, M. Inaba, and H. Inoue, “Motion planning for humanoid robots,” in _Robotics Research. The Eleventh International Symposium_. Springer, 2005, pp. 365–374.
  * [5] L. Baudouin, N. Perrin, T. Moulard, F. Lamiraux, O. Stasse, and E. Yoshida, “Real-time replanning using 3d environment for humanoid robot,” in _2011 11th IEEE-RAS International Conference on Humanoid Robots_. IEEE, 2011, pp. 584–589.
  * [6] J. Chestnutt, M. Lau, G. Cheung, J. Kuffner, J. Hodgins, and T. Kanade, “Footstep planning for the honda asimo humanoid,” in _Proceedings of the 2005 IEEE International Conference on Robotics and Automation_ , April 2005, pp. 629–634.
  * [7] P. Michel, J. Chestnutt, J. Kuffner, and T. Kanade, “Vision-guided humanoid footstep planning for dynamic environments,” in _5th IEEE-RAS International Conference on Humanoid Robots, 2005._ IEEE, 2005, pp. 13–18.
  * [8] Y.-C. Lin, B. Ponton, L. Righetti, and D. Berenson, “Efficient humanoid contact planning using learned centroidal dynamics prediction,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 5280–5286.
  * [9] Z. Xia, J. Xiong, and K. Chen, “Global navigation for humanoid robots using sampling-based footstep planners,” _IEEE/ASME Transactions on mechatronics_ , vol. 16, no. 4, pp. 716–723, 2010.
  * [10] R. Deits and R. Tedrake, “Footstep planning on uneven terrain with mixed-integer convex optimization,” in _2014 IEEE-RAS International Conference on Humanoid Robots_. IEEE, 2014, pp. 279–286.
  * [11] ——, _Computing Large Convex Regions of Obstacle-Free Space Through Semidefinite Programming_. Cham: Springer International Publishing, 2015, pp. 109–124.
  * [12] C. I. Vasile and C. Belta, “Sampling-based temporal logic path planning,” in _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2013, pp. 4817–4822.
  * [13] S. Karaman and E. Frazzoli, “Sampling-based motion planning with deterministic $\mu$-calculus specifications,” in _Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference_. IEEE, 2009, pp. 2222–2229.
  * [14] J. P. Vielma, “Mixed integer linear programming formulation techniques,” _Siam Review_ , vol. 57, no. 1, pp. 3–57, 2015.
  * [15] C. Baier and J.-P. Katoen, _Principles of model checking_ , 2008.
  * [16] A. Biere, K. Heljanko, T. Junttila, T. Latvala, and V. Schuppan, “Linear encodings of bounded ltl model checking,” _Logical Methods in Computer Science_ , vol. 2, pp. 1–64, 11 2006.
  * [17] E. M. Wolff, U. Topcu, and R. M. Murray, “Optimization-based trajectory generation with linear temporal logic specifications,” in _2014 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2014, pp. 5319–5325.
  * [18] M. Posa, C. Cantu, and R. Tedrake, “A direct method for trajectory optimization of rigid bodies through contact,” _The International Journal of Robotics Research_ , vol. 33, no. 1, pp. 69–81, 2014. [Online]. Available: https://doi.org/10.1177/0278364913506757