# Conformal Prediction Regions for Time Series
using Linear Complementarity Programming

Matthew Cleaveland† † Department of Electrical and Systems Engineering,
University of Pennsylvania  Insup Lee⋄ † Department of Electrical and Systems
Engineering, University of Pennsylvania  George J. Pappas† † Department of
Electrical and Systems Engineering, University of Pennsylvania  and Lars
Lindemann∗ † Department of Electrical and Systems Engineering, University of
Pennsylvania

(March 2023)

###### Abstract

Conformal prediction is a statistical tool for producing prediction regions of
machine learning models that are valid with high probability. However,
applying conformal prediction to time series data leads to conservative
prediction regions. In fact, to obtain prediction regions over $T$ time steps
with confidence $1-\delta$, previous works require that each individual
prediction region is valid with confidence $1-\delta/T$. We propose an
optimization-based method for reducing this conservatism to enable long
horizon planning and verification when using learning-enabled time series
predictors. Instead of considering prediction errors individually at each time
step, we consider a parameterized prediction error over multiple time steps.
By optimizing the parameters over an additional dataset, we find prediction
regions that are not conservative. We show that this problem can be cast as a
mixed integer linear complementarity program (MILCP), which we then relax into
a linear complementarity program (LCP). Additionally, we prove that the
relaxed LP has the same optimal cost as the original MILCP. Finally, we
demonstrate the efficacy of our method on a case study using pedestrian
trajectory predictors.

##  1 Introduction

Autonomous systems perform safety-critical tasks in dynamic and uncertain
environments where system errors can be dangerous and costly, e.g., a mistake
of a self-driving car navigating in urban traffic. In such scenarios, accurate
uncertainty quantification is vital for ensuring safety of learning-enabled
systems. In this work, we focus on quantifying uncertainty of time series
data. For instance, for a given time series predictor that predicts the future
trajectory of a pedestrian, how can we quantify the uncertainty and accuracy
of this predictor?

Conformal prediction (CP) has emerged as a popular method for statistical
uncertainty quantification [1, 2]. It aims to construct prediction regions
that contain the true quantity of interest with a user-defined probability. CP
does not require any assumptions about the underlying distribution or the
predictor itself. Instead, one only needs calibration data that is
exchangeable or independent and identically distributed. That means that CP
can seamlessly be applied to learning-enabled predictors like neural networks,
without the need to analyze the underlying architecture or change the training
procedure [3].

One challenge of CP is that its standard variant can not directly be used to
construct valid prediction regions for time series. This is because datapoints
at different time steps are not exchangeable. Recently, variants of CP have
been developed for time series, such as adaptive conformal inference [4, 5].
However, their coverage guarantees are weaker and only asymptotic. While these
types of guarantees are useful in many contexts, they are clearly insufficient
in many safety-critical applications such as self-driving cars.

Figure 1: Sample trajectory of a pedestrian (teal stars), the pedestrian
predictions using a social LSTM from [6] (blue stars), the conformal
prediction regions for our approach (green circles), and the conformal
prediction regions for the approach from [7] (red circles). Notably, our
approach is less conservative while still providing valid prediction regions.

In previous works [8, 7, 9], this issue was mitigated by instead considering a
calibration dataset that consists of full rollouts of the time series, which
do not violate the exchangeability assumption. At design time, one runs a
separate instance of CP for each of the $T$ time steps of the time series
predictor (also referred to as a trajectory predictor). Then at runtime, these
CP instances are combined with the predictor’s $T$ predictions to compute
prediction regions. However, to obtain prediction regions that are valid over
all $T$ time steps with confidence $1-\delta$, each individual prediction
region has to be valid with confidence $1-\delta/T$.111A simple union bound or
Bonferroni inequality argument can be used to show this result. This results
in conservative prediction regions.

Contributions:  To address this conservatism, this paper proposes constructing
prediction regions for time series using linear complementary programming. Our
main idea is based on i) mapping trajectory prediction errors to a single
parameterized conformal scoring function to avoid union bounding, and ii)
finding the optimal parameters from an additional calibration dataset. As
opposed to existing works, our method enables non-conservative long horizon
planning and verification. Our contributions are as follows:

  * •

We propose the conformal scoring function
$R:=\max(\alpha_{1}R_{1},\ldots,\alpha_{T}R_{T})$ for time horizon $T$,
prediction errors $R_{t}$, and parameters $\alpha_{t}$. Using $R$, we show
that we can obtain non-conservative prediction regions with $1-\delta$
confidence using conformal prediction.

  * •

We formulate the problem of finding the parameters $\alpha_{t}$ that minimize
$R$ from an additional calibration dataset as a mixed integer linear
complementarity program, and we show how we can derive a linear
complementarity program that has the same optimal value.

  * •

On a pedestrian case study, we demonstrate that our method produces much
smaller valid conformal prediction regions compared to [8, 7, 9], see Figure
1.

###  1.1 Related Work

Conformal prediction was originally proposed by Vovk in [2, 1] to quantify
uncertainty of prediction models. The original method, however, required
training a prediction model for each calibration datapoint, which is
computationally intractable for learning-enabled predictors. Inductive
conformal prediction, also referred to as split conformal prediction,
addresses this issue by splitting the calibration data into two sets, one used
for training the prediction model and one for applying conformal prediction
[10]. Split conformal prediction has been extended to provide conditional
probabilistic guarantees [11], to handle distribution shifts [12, 13], and to
allow for quantile regression [14]. Further, split conformal prediction has
been used to construct probably approximately correct prediction sets for
machine learning models [15, 16], to perform out-of-distribution detection
[17, 18], to guarantee safety in autonomous systems [19], and to quantify
uncertainty for F1/10 car motion predictions [20]. Additionally, in [21] the
authors encode the width of the generated prediction sets directly into the
loss function of a neural network during training.

However, the aforementioned methods cannot directly be applied to time series
data. Multiple works have adapted conformal prediction algorithms for the time
series domain, including enbpi [22], adaptive conformal inference [4], and
fully adaptive conformal inference [23, 24]. However, they only provide
averaged coverage guarantees over long horizons. Regardless, these methods
have been used for analyzing financial market data [25], synthesizing safe
robot controllers [26], and dynamically allocating compute resources [27]. For
a more comprehensive overview of the conformal prediction landscape see [3].

Our paper is motivated by [8, 7, 9]. These works apply inductive conformal
prediction in a time series setting, which requires construction individual
prediction regions with confidence $1-\delta/T$ to obtain $1-\delta$ coverage
for all time steps. This results in unnecessarily conservative prediction
regions that are impractical for downstream tasks such as motion planning. We
overcome this limitation by considering a parameterized conformal scoring
function that is defined over multiple time steps.

##  2 Background

In this section, we present background on time series predictors and split
conformal prediction for time series.

###  2.1 Time Series Predictors

We want to predict future values $Y_{1},\ldots,Y_{T}$ of a time series for a
prediction horizon of $T$ from past observed values $Y_{T_{obs}},\ldots,Y_{0}$
for an observation length of ${T_{obs}}$. Let $\mathcal{D}$ denote an unknown
distribution over a finite-horizon time series
$Y:=(Y_{T_{obs}},\ldots,Y_{0},Y_{1},\ldots,Y_{T})$, i.e., let

$(Y_{T_{obs}},\ldots,Y_{0},Y_{1},\ldots,Y_{T})\sim\mathcal{D}$

denote a random trajectory drawn from $\mathcal{D}$ where
$Y_{t}\in\mathbb{R}^{m}$ is the value at time $t$. We make no assumptions
about the distribution $\mathcal{D}$, but we do assume availability of a
calibration dataset in which each trajectory is drawn independently from
$\mathcal{D}$. Particularly, define

$D_{cal}:=\\{Y^{(1)},\ldots,Y^{(n)}\\}$

where the element
$Y^{(i)}:=(Y_{T_{obs}}^{(i)},\ldots,Y_{0}^{(i)},Y_{1}^{(i)},\ldots,Y_{T}^{(i)})$
is independently drawn from $\mathcal{D}$, i.e., $Y^{(i)}\sim\mathcal{D}$.

A key challenge lies in constructing accurate predictions of future values
$Y_{1},\ldots,Y_{T}$ of the time series. We assume here that we are already
given a time series predictor $h:\mathbb{R}^{m(T_{obs}+1)}\to\mathbb{R}^{mT}$
that maps ($Y_{T_{obs}},\ldots,Y_{0}$) to estimates of the next $T$ values
($Y_{1},\ldots,Y_{T}$), denoted as

$(\hat{Y}_{1},\ldots,\hat{Y}_{T}):=h(Y_{T_{obs}},\ldots,Y_{0}).$

We make no assumptions about $h$. For example, it can be a long short-term
memory network (LSTM) [28] or a sliding linear predictor with extended Kalman
filters [29]. Using the predictor $h$, we can obtain prediction for the
calibration data in $D_{cal}$. Specifically, for each $Y^{(i)}\in D_{cal}$, we
obtain

$\displaystyle(\hat{Y}_{1}^{(i)},\ldots,\hat{Y}_{T}^{(i)}):=h(Y_{T_{obs}}^{(i)},\ldots,Y_{0}^{(i)}).$
(1)

The predictor $h$ may not be exact, and the inaccuracy of its predictions
$(\hat{Y}_{1},\ldots,\hat{Y}_{T})$ is unknown. We will quantify prediction
uncertainty using the calibration dataset $D_{cal}$ and split conformal
prediction, which we introduce next.

###  2.2 Split Conformal Prediction

Conformal prediction was introduced in [2, 1] to obtain valid prediction
regions for complex predictive models such as neural networks without making
assumptions on the distribution of the underlying data. Split conformal
prediction is a computationally tractable variant of conformal prediction [10]
where it is assumed that a calibration dataset is available that has not been
use to train the predictor. Let $R^{(0)},\ldots,R^{(k)}$ be $k+1$ exchangeable
random variables.222Formally, this means that the joint distribution of
$R^{(0)},\ldots,R^{(k)}$ is the same as the joint distribution of
$R^{(\sigma(0))},\ldots,R^{(\sigma(k))}$ for any permutation $\sigma$ on
$\\{0,\ldots k\\}$. Intuitively, exchangeability is a slightly weaker form of
independent and identically distributed (i.i.d.) random variables. The
variable $R^{(i)}$ is usually referred to as the _nonconformity score_. In
supervised learning, it is often defined as
$R^{(i)}:=\|Z^{(i)}-\mu(X^{(i)})\|$ where the predictor $\mu$ attempts to
predict the output $Z^{(i)}$ based on the input $X^{(i)}$. A large
nonconformity score indicates a poor predictive model. Our goal is to obtain a
prediction region for $R^{(0)}$ based on $R^{(1)},\ldots,R^{(k)}$, i.e., the
random variable $R^{(0)}$ should be contained within the prediction region
with high probability.

Formally, given a failure probability $\bar{\delta}\in(0,1)$, we want to
construct a valid prediction region $C$ so that333More formally, we would have
to write $C(R^{(1)},\ldots,R^{(k)})$ as the prediction region $C$ is a
function of $R^{(1)},\ldots,R^{(k)}$. For this reason, the probability measure
$P$ is defined over the product measure of $R^{(0)},\ldots,R^{(k)}$.

$\displaystyle\text{Prob}(R^{(0)}\leq C)\geq 1-\bar{\delta}.$ (2)

We pick
$C:=\text{Quantile}(\\{R^{(1)},\ldots,R^{(n_{1})},\infty\\},1-\bar{\delta})$
which is the $(1-\bar{\delta})$th quantile of the empirical distribution of
the values $R^{(1)},\ldots,R^{(k)}$ and $\infty$. Alternatively, by assuming
that $R^{(1)},\ldots,R^{(k)}$ are sorted in non-decreasing order and by adding
$R^{(k+1)}:=\infty$, we can equivalently obtain $C:=R^{(p)}$ where
$p:=\lceil(k+1)(1-\bar{\delta})\rceil$ with $\lceil\cdot\rceil$ being the
ceiling function, i.e., $C$ is the $p$th smallest nonconformity score. By a
quantile argument, see [12, Lemma 1], one can prove that this choice of $C$
satisfies Equation 2. We remark that $k\geq\lceil(k+1)(1-\delta)\rceil$ is
required to hold to obtain meaningful, i.e., bounded, prediction regions.

######  Remark 1.

Note that split conformal prediction assumes that the non-conformity scores
$R^{(0)},\ldots,R^{(k)}$ are exchangeable. This complicates its use for time
series data where $Y_{t+1}$ generally depends on $Y_{t}$. However, with a
calibration dataset $D_{cal}$ of trajectories independently drawn from
$\mathcal{D}$, the application of conformal prediction is possible [8, 7, 9].
In summary, one defines $T$ non-conformity scores
$R^{(i)}_{t}:=\|\hat{Y}^{(i)}_{t}-Y^{(i)}_{t}\|$ for each time
$t\in\\{1,\ldots,T\\}$ and each calibration trajectory $Y^{(i)}\in D_{cal}$;
$R^{(i)}_{t}$ can be interpreted as the $t$-step ahead prediction error. Using
conformal prediction, for each time $t$ we can then construct a prediction
region $C_{t}$ so that $\text{Prob}(\|\hat{Y}_{t}-Y_{t}\|\leq C_{t})\geq
1-\bar{\delta}$. To then obtain a prediction region over all time steps, i.e,
$\text{Prob}(\|\hat{Y}_{t}-Y_{t}\|\leq C_{t},\;\forall
t\in\\{1,\ldots,T\\})\geq 1-\delta$, we simply set $\bar{\delta}:=\delta/T$.
As previously mentioned, this results in conservative prediction regions that
are not practical.

##  3 Problem Formulation and Main Idea

Our goal is to construct valid conformal prediction regions for time series.
The main idea is to use a parameterized non-conformity score that considers
predictions errors over all $T$ time steps. Particularly, in this paper we
consider the non-conformity score

$\displaystyle
R:=\max(\alpha_{1}R_{p}(Y_{1},\hat{Y}_{1}),\ldots,\alpha_{T}R_{p}(Y_{T},\hat{Y}_{T}))$
(3)

where $\alpha_{1},\ldots,\alpha_{T}>0$ are parameters and where $R_{p}$ is the
prediction error for an individual time step. For the purpose of this paper,
we use the Euclidean distance between the actual and the predicted value of
the time series at time $t$ as

$\displaystyle R_{p}(Y_{t},\hat{Y}_{t}):=\|Y_{t}-\hat{Y}_{t}\|$ (4)

We note that $R_{p}$ can in general be any function of $Y_{t}$ and
$\hat{Y}_{t}$. Within the context of this paper, however, we set $R_{p}$ to be
the Euclidean distance between $Y_{t}$ and $\hat{Y}_{t}$ as in (4).

Taking the maximum prediction error over all $T$ time steps as in equation (4)
will allow us to avoid the overly conservative prediction regions from [8, 7,
9] as we, in our proof, do not need to union bound over the individual
prediction errors. However, the choice of the parameters
$\alpha_{1},\ldots,\alpha_{T}$ is now important. If we apply conformal
prediction as described in Section 2.2 to the nonconformity score in (3) with
$\alpha_{1}=\ldots=\alpha_{T}=1$, this would typically result in obtaining
large prediction regions for the first time step since usually prediction
errors get larger with increasing time, i.e., $R_{p}(Y_{T},\hat{Y}_{T})$ in
(3) would dominate so that $R=R_{p}(Y_{T},\hat{Y}_{T})$. Introducing the
parameters $\alpha_{1},\ldots,\alpha_{T}$ allows for weighting the different
time steps, and we generally expect $\alpha_{t}$ for larger times to be
smaller than $\alpha_{t}$ for smaller times.

Given the parameterized function $R$ as in Equation 3 and a failure threshold
of $\delta\in(0,1)$, we are in this paper particularly interested in computing
the parameters $\alpha_{1},\ldots,\alpha_{T}$ that minimize the constant $C$
that satisfies

$\displaystyle\text{Prob}(R\leq C)\geq 1-\delta.$

For given parameters $\alpha_{1},\ldots,\alpha_{T}$, note that the constant
$C$ can be found by applying conformal prediction to the nonconformity score
$R$ using the calibration set $D_{cal}$. When we have found this constant $C$,
we know that

$\displaystyle\text{Prob}(\|Y_{t}-\hat{Y}_{t}\|\leq C/\alpha_{t},\;\;\forall
t\in\\{1,\ldots,T\\})\geq 1-\delta.$

To solve the aforementioned problem, we split the calibration dataset
$D_{cal}$ into two sets $D_{cal,1}$ and $D_{cal,2}$. We will use the first
calibration dataset $D_{cal,1}$ to compute the values of
$\alpha_{1},\ldots,\alpha_{T}$ that give the smallest prediction region $C$
for $R$ (see Section 4). In a next step, we use these parameters that now
fully define the nonconformity score $R$ along with the second calibration
dataset $D_{cal,2}$ to construct valid prediction regions via conformal
prediction (see Section 5).

##  4 Nonconformity Scores for Time-Series via Linear Programming

We first have to set up some notation. Let us denote the elements of the first
calibration dataset as $D_{cal,1}:=\\{Y^{(1)},\ldots,Y^{(n_{1})}\\}$ where
$n_{1}>0$ is the size of the first calibration dataset. Recall that each
element $Y^{(i)}$ of $D_{cal,1}$ is defined as
$Y^{(i)}=(Y_{T_{obs}}^{(i)},\ldots,Y_{0}^{(i)},Y^{(i)}_{1},\ldots,Y^{(i)}_{T})$.
Using the trajectory predictor $h$, we obtain predictions
$(\hat{Y}_{1}^{(i)},\ldots,\hat{Y}_{T}^{(i)})$ for
$(Y^{(i)}_{1},\ldots,Y^{(i)}_{T})$ as per equation (1).

We then compute the prediction error $R_{t}^{(i)}$ according to equation (4)
for each calibration trajectory and for each time step, i.e., for each
$t\in\\{1,\ldots,T\\}$ and $i\in\\{1,\ldots,n_{1}\\}$. More formally, we
compute for each calibration trajectory

$\displaystyle R_{t}^{(i)}:=R_{p}(Y_{t}^{(i)},\hat{Y}_{t}^{(i)}).$ (5)

We can now cast the problem of finding the parameters
$\alpha_{1},\ldots,\alpha_{T}$ as the following optimization problem.

$\displaystyle\min_{\alpha_{1},\ldots,\alpha_{T}}$
$\displaystyle\quad\text{Quantile}(\\{R^{(1)},\ldots,R^{(n_{1})},\infty\\},1-\delta)$
(6a) s.t. $\displaystyle\quad
R^{(i)}=\text{max}(\alpha_{1}R^{(i)}_{1},\ldots,\alpha_{T}R^{(i)}_{T}),i=1,\ldots,n_{1}$
(6b) $\displaystyle\quad\sum_{j=1}^{T}\alpha_{j}=1$ (6c)
$\displaystyle\quad\alpha_{j}>0,j=1,\ldots,T$ (6d)

where we recall that
$\text{Quantile}(\\{R^{(1)},\ldots,R^{(n_{1})},\infty\\},1-\delta)$ denotes
the empirical $1-\delta$ quantile of the set
$\\{R^{(1)},\ldots,R^{(n_{1})},\infty\\}$. Equations 6c and 6d normalize
$\alpha_{1},\ldots,\alpha_{T}$ and ensure them to be positive.

######  Remark 2.

We note that the optimization problem in equation (6) is always feasible. This
follows since any set of parameters $\alpha_{1},\ldots,\alpha_{T}$ that
satisfy Equations 6c and 6d constitute a feasible solution to (6).

To solve the optimization problem in equation (6), we transform it into a
linear complementarity program [30] that has the same optimal value as (6). To
do so, we first reformulate the $\text{max}()$ operator in equation (6b) as a
set of mixed integer linear program (MILP) constraints, which we then
transform into a set of linear program constraints. We show that this
transformation preserves the optimal value of (6) (see Section 4.1). Then, we
reformulate the Quantile function from equation (6a) into a linear program
(LP) that we further reformulate using its KKT conditions (see Section 4.2),
resulting in a linear complementarity program.

###  4.1 Reformulating the Max Operator

We first use ideas from [31, 32] for encoding the max operator in (6b) as an
MILP. Therefore, we introduce the binary variables $b_{t}^{(i)}\in\\{0,1\\}$
for each time $t\in\\{1,\ldots,T\\}$ and calibration trajectory
$i\in\\{1,\ldots,n_{1}\\}$. The equality constraint
$R^{(i)}=\text{max}(\alpha_{1}R_{1}^{(i)},\ldots,\alpha_{T}R_{T}^{(i)})$ in
equation (6b) is then equivalent to the following MILP constraints:

$\displaystyle\quad R^{(i)}\geq\alpha_{t}R_{t}^{(i)},\;t=1,\ldots,T$ (7a)
$\displaystyle\quad
R^{(i)}\leq\alpha_{t}R_{t}^{(i)}+(1-b_{t}^{(i)})M,\;t=1,\ldots,T$ (7b)
$\displaystyle\quad\sum_{t=1}^{T}b_{t}^{(i)}=1$ (7c) $\displaystyle\quad
b_{t}^{(i)}\in\\{0,1\\},\;t=1,\ldots,T$ (7d)

where $M>0$ is a sufficiently large and positive constant.444See [31] for the
exact conditions required of $M$.

By replacing the max operator in Equation 6b with the MILP encoding in
Equation 7, we arrive at the following optimization problem:

$\displaystyle\min_{\alpha_{1},\ldots,\alpha_{T}}$
$\displaystyle\quad\text{Quantile}(\\{R^{(1)},\ldots,R^{(n_{1})},\infty\\},1-\delta)$
(8a) s.t. $\displaystyle\quad
R^{(i)}\geq\alpha_{t}R_{t}^{(i)},\;i=1,\ldots,n_{1}\land t=1,\ldots,T$ (8b)
$\displaystyle\quad
R^{(i)}\leq\alpha_{t}R_{t}^{(i)}+(1-b_{t}^{(i)})M,\;i=1,\ldots,n_{1}$ (8c)
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\;\;\land
t=1,\ldots,T$ $\displaystyle\quad\sum_{t=1}^{T}b_{t}^{(i)}=1$ (8d)
$\displaystyle\quad b_{t}^{(i)}\in\\{0,1\\},\;t=1,\ldots,T$ (8e)
$\displaystyle\quad\eqref{eq:highLevelOptimizationc},\;\eqref{eq:highLevelOptimizationd}$

The constraints in (8) are linear in their parameters, i.e., in $\alpha_{t}$,
$R^{(i)}$, and $b_{t}^{(i)}$. The cost function in (8a) can be reformulated as
a linear complementarity program (LCP), as shown in the next section. Thus,
(8) is an MILCP since the $b_{t}^{(i)}$ variables are binary due to equation
(8e). Commercial solvers can solve MILCPs, but they present scalability issues
for larger time horizons $T$ and calibration data $n_{1}$. Therefore, we
present an LCP reformulation of (8) which will achieve the same optimality in
terms of the constraint function (8a) while providing better scalability.

####  4.1.1 Linear complementarity relaxation of the MILCP in (8)

We start by removing the upper bound on $R^{(i)}$ in equation (8c) and the
binary variables $b_{t}^{(i)}$ in equations (8d) and (8e) from the MILCP in
(8). This results in the optimization problem:

$\displaystyle\min_{\alpha_{1},\ldots,\alpha_{T}}$
$\displaystyle\quad\text{Quantile}(\\{R^{(1)},\ldots,R^{(n_{1})},\infty\\},1-\delta)$
(9a) s.t. $\displaystyle\quad
R^{(i)}\geq\alpha_{t}R_{t}^{(i)},\;i=1,\ldots,n_{1}\land t=1,\ldots,T$ (9b)
$\displaystyle\quad\eqref{eq:highLevelOptimizationc},\;\eqref{eq:highLevelOptimizationd}$

Note that the above optimization problem is now a linear complementarity
program due to the removal of the binary variables $b_{t}^{(i)}$. While the
LCP in (9) is clearly computationally more tractable than the MILCP in (8),
one may ask in what way the optimal solutions of these two optimization
problems are related. We next prove that the optimal values (in terms of the
cost functions in (8a) and (9a)) of the two optimization problems (8) and (9)
are the same.

######  Theorem 1.

Let $\delta\in(0,1)$ be a failure rate and $R_{t}^{(i)}$ be the prediction
errors as in (5) for times $t\in\\{1,\ldots,T\\}$ and calibration trajectories
$i\in\\{1,\ldots,n_{1}\\}$ from $D_{cal,1}$. Then, the optimal cost values of
the optimization problems (8) and (9) are equivalent.

###### Proof.

Assume that
$\alpha_{1}^{*},\ldots,\alpha_{T}^{*},R^{(1),*},\ldots,R^{(n_{1}),*}$ is an
optimal solution to the optimization problem in equation (9). We will first
show that we can use
$\alpha_{1}^{*},\ldots,\alpha_{T}^{*},R^{(1),*},\ldots,R^{(n_{1}),*}$ to
construct a feasible solution to the optimization problem in (8), which we
denote by
$\bar{\alpha}_{1}^{*},\ldots,\bar{\alpha}_{T}^{*},\bar{R}^{(1),*},\ldots,\bar{R}^{(n_{1}),*}$,
without altering the optimal value of the objective function. Particularly,
let

$\displaystyle\bar{\alpha}_{t}^{*}$
$\displaystyle:=\alpha_{t}^{*},\;t=1,\ldots,T$ (10a)
$\displaystyle\bar{R}^{(i),*}$
$\displaystyle:=\text{max}(\bar{\alpha}^{*}_{1}R^{(i)}_{1},\ldots,\bar{\alpha}^{*}_{T}R^{(i)}_{T}),\;i=1,\ldots,n_{1}$
(10b)

Note that there exists binary variables $\bar{b}_{t}^{(i)}$ along with which
$\bar{\alpha}_{1}^{*},\ldots,\bar{\alpha}_{T}^{*},\bar{R}^{(1),*},\ldots,\bar{R}^{(n_{1}),*}$
are a feasible solution to the optimization problem in (8). This follows since
(10a) satisfies (6c) and (6d) and (10b) satisfies (6b) and consequently
(8b-8d) for appropriate choices of $\bar{b}_{t}^{(i)}$ (and $M$).

Note specifically that $\bar{\alpha}_{t}^{*}$ and $\alpha_{t}^{*}$ are the
same, while this does not necessarily hold for ${R}^{(i),*}$ and
$\bar{R}^{(i),*}$. However, it is guaranteed that $\bar{R}^{(i),*}\leq
R^{(i),*}\;\forall i=1,\ldots,n_{1}$ as the constraint in (9b) ensures that
$R^{(i),*}$ is lower bounded by
$\text{max}(\alpha^{*}_{1}R^{(i)}_{1},\ldots,\alpha^{*}_{T}R^{(i)}_{T})=\text{max}(\bar{\alpha}^{*}_{1}R^{(i)}_{1},\ldots,\bar{\alpha}^{*}_{T}R^{(i)}_{T})$.
Since the Quantile function is a monotone operator (its value is non-
increasing if inputs are non-increasing), the value of the objective function
in (8) that corresponds to
$\bar{\alpha}_{1}^{*},\ldots,\bar{\alpha}_{T}^{*},\ldots,\bar{R}^{(1),*},\bar{R}^{(n_{1}),*}$
can not be greater than the value of the objective function in (9) that
corresponds to
$\alpha_{1}^{*},\ldots,\alpha_{T}^{*},R^{(1),*},\ldots,R^{(n_{1}),*}$.

At the same time, we note that the optimal value of the optimization problem
in (9) cannot be greater than the optimal value of the optimization problem in
(8) because it has strictly fewer constraints than (9) so that the feasible
set of (9) contains the feasible of (8). As a consequence, the optimal values
of (8) and (9) are equivalent. ∎

Note that Theorem 1 guarantees that the optimal values of the optimization
problems in (8) and (9) are equivalent. However, we do not guarantee that the
optimal parameters $\alpha_{1}^{*},\ldots,\alpha_{T}^{*}$ obtained by solving
(8) and (9) will be the same as the optimization problem in (8) may have
multiple optimal solutions $\alpha_{1}^{*},\ldots,\alpha_{T}^{*}$.
Nonetheless, the linear reformulation in (9) is computationally tractable and
cost optimal. In our experiments in Section 6.3.1, we also show that we obtain
parameters from (8) and (9) that are almost equivalent.

###  4.2 Reformulating the Quantile Function

We note that the Quantile function in (6a) can be written as the following
linear program [33]:

Quantile $\displaystyle(\\{R^{(1)},\ldots,R^{(n_{1})}\\},1-\delta)$
$\displaystyle=\text{argmin}_{q}$
$\displaystyle\quad\sum_{i=1}^{n_{1}}\left((1-\delta)e_{i}^{+}+\delta
e_{i}^{-}\right)$ (11a) s.t. $\displaystyle\quad
e_{i}^{+}-e_{i}^{-}=R^{(i)}-q,\;i=1,\ldots,n_{1}$ (11b) $\displaystyle\quad
e_{i}^{+},e_{i}^{-}\geq 0,\;i=1,\ldots,n_{1}$ (11c)

where the optimal solution $q$ of (11) is equivalent to
$\text{Quantile}(\\{R^{(1)},\ldots,R^{(n_{1})}\\},1-\delta)$. However, we note
that we can not directly replace
$\text{Quantile}(\\{R^{(1)},\ldots,R^{(n_{1})}\\},1-\delta)$ in the
optimization problem in (6a) with the linear program from equation (11) as
that would result in a problem with an objective function that is itself an
optimization problem. To address this issue, we replace the optimization
problem in (11) with its KKT conditions that are as follows:

$\displaystyle(1-\delta)-u_{i}^{+}+v_{i}=0,\;\;i=1,\ldots,n_{1}$ (12a)
$\displaystyle\delta-u_{i}^{-}-v_{i}=0,\;\;i=1,\ldots,n_{1}$ (12b)
$\displaystyle\sum_{i=1}^{n_{1}}v_{i}=0$ (12c) $\displaystyle
u_{i}^{+}e_{i}^{+}=0,\;\;i=1,\ldots,n_{1}$ (12d) $\displaystyle
u_{i}^{-}e_{i}^{-}=0,\;\;i=1,\ldots,n_{1}$ (12e) $\displaystyle e_{i}^{+}\geq
0,\;\;i=1,\ldots,n_{1}$ (12f) $\displaystyle e_{i}^{-}\geq
0,\;\;i=1,\ldots,n_{1}$ (12g) $\displaystyle
e_{i}^{+}+q-e_{i}^{-}-R^{(i)}=0,\;\;i=1,\ldots,n_{1}$ (12h) $\displaystyle
u_{i}^{+}\geq 0,\;\;i=1,\ldots,n_{1}$ (12i) $\displaystyle u_{i}^{-}\geq
0,\;\;i=1,\ldots,n_{1}$ (12j)

Here, the constraints (12a), (12b) and (12c) encode the stationarity
condition, the constraints (12d) and (12e) encode complementary slackness, the
constraints (12f), (12g), and (12h) encode primal feasibility, and the
constraints (12i) and (12j) encode dual feasibility. This is an LCP because
every constraint is either linear (((12a)-(12c)) and ((12f)-(12j))) or an
equality constraint stating the multiple of two positive variables equals $0$
((12d) and (12e)).555These constraints ensure that one of the two variables
equals $0$. Despite the slight non-linearity, LCPs can be efficiently solved
with existing tools [34].

Because feasible linear optimization problems have zero duality gap [35], any
solution to the KKT conditions in (12) is also an optimal solution to Problem
(11).

###  4.3 Summarizing the Optimization Programs

Based on the results from the previous subsections, we can formulate two
optimization problems that solve (6).

First, we can exactly solve (6) by replacing the max operator in (6b) with the
MILP in (7) and the quantile function in (6a) with its KKT conditions in (12).
The following result follows by construction (for sufficiently large $M$).

######  Proposition 1.

The mixed integer linear complementarity program

$\displaystyle\text{min}_{q,\alpha_{1},\ldots,\alpha_{T}}$ $\displaystyle\quad
q$ (13a) s.t.
$\displaystyle\quad\eqref{eq:highLevelOptimizationMaxMILPb},\;\eqref{eq:highLevelOptimizationMaxMILPc},\;\eqref{eq:highLevelOptimizationMaxMILPd},\;\eqref{eq:highLevelOptimizationMaxMILPe}$
$\displaystyle\quad\eqref{eq:KKTcondsa},\;\eqref{eq:KKTcondsb},\;\eqref{eq:KKTcondsc},\;\eqref{eq:KKTcondsd},\;\eqref{eq:KKTcondse}$
$\displaystyle\quad\eqref{eq:KKTcondsf},\;\eqref{eq:KKTcondsg},\;\eqref{eq:KKTcondsh},\;\eqref{eq:KKTcondsi},\;\eqref{eq:KKTcondsj}$
$\displaystyle\quad\eqref{eq:highLevelOptimizationc},\;\eqref{eq:highLevelOptimizationd}$
is equivalent to the optimization problem in (6).

Second, by further removing the upper bound from the max operator as described
in Section 4.1.1, we can solve a linear complementarity program that has the
same optimal cost value as (6). This result follows from Theorem 2.

######  Proposition 2.

The linear complementarity program

$\displaystyle\text{min}_{q,\alpha_{1},\ldots,\alpha_{T}}$ $\displaystyle\quad
q$ (14a) s.t. $\displaystyle\quad\eqref{eq:highLevelOptimizationMaxMILPb}$
$\displaystyle\quad\eqref{eq:KKTcondsa},\;\eqref{eq:KKTcondsb},\;\eqref{eq:KKTcondsc},\;\eqref{eq:KKTcondsd},\;\eqref{eq:KKTcondse}$
$\displaystyle\quad\eqref{eq:KKTcondsf},\;\eqref{eq:KKTcondsg},\;\eqref{eq:KKTcondsh},\;\eqref{eq:KKTcondsi},\;\eqref{eq:KKTcondsj}$
$\displaystyle\quad\eqref{eq:highLevelOptimizationc},\;\eqref{eq:highLevelOptimizationd}$

has the same optimal value as the optimization problem (6).

##  5 Conformal Prediction Regions for Time-Series

Finally, after obtaining the parameters $\alpha_{1},\ldots,\alpha_{T}$ from
the first calibration dataset $D_{cal,1}$ by solving (13) or (14), we can
apply conformal prediction with the non-conformity score $R$ as defined in (3)
to obtain prediction regions for time series. We do so by following Section
2.2 and using the second calibration dataset $D_{cal,2}$.

More formally, let us denote the elements of the second calibration dataset as
$D_{cal,2}:=\\{Y^{(n_{1}+1)},\ldots,Y^{(n)}\\}$ where $n-n_{1}>0$ is the size
of the second calibration dataset.666Recall the $n$ is the size of $D_{cal}$
and $n_{1}$ is the size of $D_{cal,1}$. For each calibration trajectory
$i\in\\{n_{1}+1,\ldots,n\\}$, we compute the nonconformity score as

$\displaystyle
R^{(i)}:=\max\left(\alpha_{1}R_{p}(Y_{1}^{(i)},\hat{Y}_{1}^{(i)}),\ldots,\alpha_{T}R_{p}(Y_{T}^{(i)},\hat{Y}_{T}^{(i)})\right).$
(15)

We are now ready to state the final result of our paper in which we obtain
prediction regions for
$R:=\max\left(\alpha_{1}R_{p}(y_{1},\hat{y}_{1}),\ldots,\alpha_{T}R_{p}(y_{T},\hat{y}_{T})\right)$.

######  Theorem 2.

Let $Y$ be a trajectory drawn from $\mathcal{D}$ and let
$\hat{Y}=h(Y_{T_{obs}},\ldots,Y_{0})$ be predictions from the time series
predictor $h$. Let $\delta\in(0,1)$ be a failure rate and $R_{t}^{(i)}$ be the
prediction error as in (5) for times $t\in\\{1,\ldots,T\\}$ and calibration
trajectories $i\in\\{1,\ldots,n_{1}\\}$ from $D_{cal,1}$. Let
$\alpha_{1},\ldots,\alpha_{T}$ be obtained by solving (13) or (14). Further,
let $R^{(i)}$ be the nonconformity score as in (15) for calibration
trajectories $i\in\\{n_{1}+1,\ldots,n\\}$ from $D_{cal,2}$. Finally, let
$R=\max\left(\alpha_{1}R_{p}(Y_{1},\hat{Y}_{1}),\ldots,\alpha_{T}R_{p}(Y_{T},\hat{Y}_{T})\right)$.
It then holds that

$\displaystyle\text{Prob}\left(R\leq C\right)\geq 1-\delta$ (16)

where $C:=\text{Quantile}(\\{R^{(n_{1}+1)},\ldots,R^{(n)}\\},1-\delta)$.

###### Proof.

Since $\alpha_{1},\ldots,\alpha_{T}$ are obtained from $D_{cal,1}$, the non-
conformity scores $R^{(i)}$ which are computed from $D_{cal,2}$ are
exchangeable. Consequently, we can apply [12, Lemma 1] to get that Equation 16
holds. ∎

We point out the necessity to split the dataset $D_{cal}$ into the datasets
$D_{cal,1}$ and $D_{cal,2}$ to conform to the exchangeability assumption in
conformal prediction. To convert equation (16) into practical conformal
prediction regions, note that

$\displaystyle R=$
$\displaystyle\text{max}\left(\alpha_{1}\|Y_{1}-\hat{Y}_{1}\|,\ldots,\alpha_{T}\|Y_{T}-\hat{Y}_{T}\|\right)\leq
C$ $\displaystyle\Leftrightarrow\;\;\;\;\;\;\;\;$
$\displaystyle\|Y_{t}-\hat{Y}_{t}\|\leq C/\alpha_{t}\;\;\forall
t\in\\{1,\ldots,T\\}$

Intuitively, the conformal prediction regions provide geometric balls of
radius $C/\alpha_{t}$ around each prediction $\hat{y}_{t}$ due to the use of
the Euclidean norm. For an example, we refer to Figure 1. We remark that other
choices of the prediction $R_{p}$ may lead to different shapes. We summarize
this result next.

######  Corollary 1.

Let all conditions of Theorem 2 hold. Then, it follows that

$\displaystyle\text{Prob}(\|Y_{t}-\hat{Y}_{t}\|\leq C/\alpha_{t}\;\;\forall
t\in\\{1,\ldots,T\\})\geq 1-\delta.$

##  6 Case Study

We compare our method to the approach from [8, 7, 9] that is described in
Section 2.2, which requires conservative union bounding to obtain valid
prediction regions. We colloquially refer to this method as the union bound
approach. We use Gurobi [36] to solve the optimization problems on a Windows
machine running an Intel i7-855OU CPU with 4 cores and 16GB of RAM. The code
for this case study can be found on Github.

###  6.1 ORCA and LSTM

For the case study, we analyze pedestrian location predictions using data
generated from the ORCA simulator [37]. We use the social LSTM [6]
architecture to make predictions up to $2.5$ seconds into the future at a rate
of 8Hz (so that $T=20$ predictions in total) with the previous $2.5$ seconds
as input (so that $T_{obs}=20$) using the trajnet++ framework [38]. For the
exact details of training, we refer the reader to our previous work in [7].

We collect a dataset of $1291$ trajectories, each with $20$ future
predictions. We randomly select $646$ of these points for $D_{cal}$ and use
the rest to form $D_{val}$, which we use for validation. We randomly select
$50$ trajectories from $D_{cal}$ for $D_{cal,1}$ (to compute the $\alpha$
values) and the remaining trajectories for $D_{cal,2}$ (to compute the
conformal prediction regions). We use (14) to compute the $\alpha$ values.
Throughout this section, we set $\delta:=0.05$. In Section 6.2 we show the
results for one realization of this process and in Section 6.3 we run $100$
trials to statistically analyze the results777 This is a random process, as
splitting the data into $D_{cal,1}$, $D_{cal,2}$, and $D_{val}$ is random and
the guarantees of conformal prediction are taken over the distribution of the
calibration and validation data. Thus we need to statistically analyze the
results.

###  6.2 Single Trial Analysis

For one of the realizations, the $\alpha$ values are shown in Table 1. They
decrease as the prediction horizon increases, which is what we expect as the
prediction error generally increases over time. For that same trial, a
trajectory of the actual and predicted pedestrian locations is shown alongside
the conformal prediction regions in Figure 1. Figure 2 shows scatter plots of
the prediction errors over $D_{val}$ and the conformal prediction region radii
for $3$ time steps.

Table 1: Values of $\alpha$ parameters for ORCA case study. Time step | $\alpha$ | Time step | $\alpha$
---|---|---|---
1 | 0.21134 | 11 | 0.0267
2 | 0.1495 | 12 | 0.0229
3 | 0.1057 | 13 | 0.0207
4 | 0.0945 | 14 | 0.0189
5 | 0.0668 | 15 | 0.0171
6 | 0.0555 | 16 | 0.0158
7 | 0.0423 | 17 | 0.0148
8 | 0.0371 | 18 | 0.0138
9 | 0.0330 | 19 | 0.0125
10 | 0.0290 | 20 | 0.0118

(a) 5 step prediction error

(b) 10 step prediction error

(c) 15 step prediction error

Figure 2: Conformal prediction regions and errors for various prediction
horizons. The green circles show the regions produced by our approach and the
red circles show the regions produced by the union bound approach. Notably,
our approach produces much tighter valid predictions regions.

###  6.3 Statistical Analysis

Over $100$ trials, the average (trajectory level) coverage for our approach
was $0.9558$ and the average coverage for the union bound approach was
$0.9945$. Histograms of the validation coverage for the two approaches are
shown in Figure 3. The union bound approach has larger coverage due to the
conservatism introduced by union bounding. However, this highlights that our
method produces valid conformal prediction regions that are not overly
conservative and significantly tighter than previous approaches.

To further quantify the how much less conservative our prediction regions are,
we also computed the average size of the conformal regions for each time
horizon, which are shown in Figure 4. This again demonstrates that our
approach results in significantly smaller prediction regions.

Figure 3: Validation coverage histograms for each of the two methods for $100$
trials.  Figure 4: Average area of the conformal prediction regions for each
of the three approaches over $100$ trials.

####  6.3.1 Original v. Relaxed Problem Difference

For each trial we also computed the objective function and $\alpha$ values
using Equation 13 to compare to those computed using Equation 14. The average
difference between the objective function values is $1.847*10^{-9}$. The
average $l^{2}$ norm difference between the $\alpha$ values is $3.495*10^{-7}$
and the largest $l^{2}$ norm difference is $3.453*10^{-5}$. These differences
are well within the solver’s default tolerance of $10^{-4}$, which empirically
validates the equivalence from Theorem 1. The average runtime for solving the
original problem was $0.525$ seconds and the relaxed problem was $0.207$
seconds. So the relaxation cuts the computation time in half.

##  7 Conclusion

In this paper, we have presented a method for producing conformal prediction
regions for time series data that are significantly tighter than previous
approaches. To do this, we define a parameterized non-conformity score
function and use optimization to fit the function to calibration data. This
allows us to use standard inductive conformal prediction to get valid
prediction regions. For future work, we plan to integrate our conformal
prediction approach into planning and control frameworks and apply
parameterized non-conformity score functions to the non-time series setting.

## Acknowledgements

This work was generously supported by NSF award CPS-2038873.

## References

  * [1] G. Shafer and V. Vovk, “A tutorial on conformal prediction.” _Journal of Machine Learning Research_ , vol. 9, no. 3, 2008.
  * [2] V. Vovk, A. Gammerman, and G. Shafer, _Algorithmic learning in a random world_. Springer Science & Business Media, 2005.
  * [3] A. N. Angelopoulos and S. Bates, “A gentle introduction to conformal prediction and distribution-free uncertainty quantification,” _arXiv preprint arXiv:2107.07511_ , 2021.
  * [4] I. Gibbs and E. Candes, “Adaptive conformal inference under distribution shift,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 1660–1672, 2021.
  * [5] M. Zaffran, O. Féron, Y. Goude, J. Josse, and A. Dieuleveut, “Adaptive conformal predictions for time series,” in _International Conference on Machine Learning_. PMLR, 2022, pp. 25 834–25 866.
  * [6] A. Alahi, K. Goel, V. Ramanathan, A. Robicquet, L. Fei-Fei, and S. Savarese, “Social lstm: Human trajectory prediction in crowded spaces,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 961–971.
  * [7] L. Lindemann, M. Cleaveland, G. Shim, and G. J. Pappas, “Safe planning in dynamic environments using conformal prediction,” _arXiv preprint arXiv:2210.10254_ , 2022.
  * [8] K. Stankeviciute, A. M Alaa, and M. van der Schaar, “Conformal time-series forecasting,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 6216–6228, 2021.
  * [9] L. Lindemann, X. Qin, J. V. Deshmukh, and G. J. Pappas, “Conformal prediction for stl runtime verification,” _arXiv preprint arXiv:2211.01539_ , 2022.
  * [10] H. Papadopoulos, “Inductive conformal prediction: Theory and application to neural networks,” in _Tools in artificial intelligence_. Citeseer, 2008.
  * [11] V. Vovk, “Conditional validity of inductive conformal predictors,” in _Asian conference on machine learning_. PMLR, 2012, pp. 475–490.
  * [12] R. J. Tibshirani, R. Foygel Barber, E. Candes, and A. Ramdas, “Conformal prediction under covariate shift,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [13] C. Fannjiang, S. Bates, A. Angelopoulos, J. Listgarten, and M. I. Jordan, “Conformal prediction for the design problem,” _arXiv preprint arXiv:2202.03613_ , 2022.
  * [14] Y. Romano, E. Patterson, and E. Candes, “Conformalized quantile regression,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [15] S. Park, O. Bastani, N. Matni, and I. Lee, “Pac confidence sets for deep neural networks via calibrated prediction,” _arXiv preprint arXiv:2001.00106_ , 2019.
  * [16] A. N. Angelopoulos, S. Bates, A. Fisch, L. Lei, and T. Schuster, “Conformal risk control,” _arXiv preprint arXiv:2208.02814_ , 2022.
  * [17] R. Kaur, S. Jha, A. Roy, S. Park, E. Dobriban, O. Sokolsky, and I. Lee, “idecode: In-distribution equivariance for conformal out-of-distribution detection,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 36, no. 7, pp. 7104–7114, Jun. 2022. [Online]. Available: https://ojs.aaai.org/index.php/AAAI/article/view/20670
  * [18] R. Kaur, K. Sridhar, S. Park, S. Jha, A. Roy, O. Sokolsky, and I. Lee, “Codit: Conformal out-of-distribution detection in time-series data,” _arXiv preprint arXiv:2207.11769_ , 2022.
  * [19] R. Luo, S. Zhao, J. Kuck, B. Ivanovic, S. Savarese, E. Schmerling, and M. Pavone, “Sample-efficient safety assurances using conformal prediction,” in _Algorithmic Foundations of Robotics XV: Proceedings of the Fifteenth Workshop on the Algorithmic Foundations of Robotics_. Springer, 2022, pp. 149–169.
  * [20] R. Tumu, L. Lindemann, T. Nghiem, and R. Mangharam, “Physics constrained motion prediction with uncertainty quantification,” _arXiv preprint arXiv:2302.01060_ , 2023.
  * [21] D. Stutz, Krishnamurthy, Dvijotham, A. T. Cemgil, and A. Doucet, “Learning optimal conformal classifiers,” 2022.
  * [22] C. Xu and Y. Xie, “Conformal prediction for dynamic time-series,” _arXiv preprint arXiv:2010.09107_ , 2020.
  * [23] I. Gibbs and E. Candès, “Conformal inference for online prediction with arbitrary distribution shifts,” _arXiv preprint arXiv:2208.08401_ , 2022\.
  * [24] A. Bhatnagar, H. Wang, C. Xiong, and Y. Bai, “Improved online conformal prediction via strongly adaptive online learning,” _arXiv preprint arXiv:2302.07869_ , 2023.
  * [25] W. Wisniewski, D. Lindsay, and S. Lindsay, “Application of conformal prediction interval estimations to market makers’ net positions,” in _International Symposium on Conformal and Probabilistic Prediction with Applications_ , 2020.
  * [26] A. Dixit, L. Lindemann, S. Wei, M. Cleaveland, G. J. Pappas, and J. W. Burdick, “Adaptive conformal prediction for motion planning among dynamic agents,” 2022\. [Online]. Available: https://arxiv.org/abs/2212.00278
  * [27] K. M. Cohen, S. Park, O. Simeone, P. Popovski, and S. Shamai, “Guaranteed dynamic scheduling of ultra-reliable low-latency traffic via conformal prediction,” _arXiv preprint arXiv:2302.07675_ , 2023.
  * [28] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [29] S. X. Wei, A. Dixit, S. Tomar, and J. W. Burdick, “Moving obstacle avoidance: A data-driven risk-aware approach,” _IEEE Control Systems Letters_ , vol. 7, pp. 289–294, 2023.
  * [30] R. W. Cottle, J.-S. Pang, and R. E. Stone, _The linear complementarity problem_. SIAM, 2009.
  * [31] A. Bemporad and M. Morari, “Control of systems integrating logic, dynamics, and constraints,” _Automatica_ , vol. 35, no. 3, pp. 407–427, 1999. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0005109898001782
  * [32] V. Raman, A. Donzé, M. Maasoumy, R. M. Murray, A. Sangiovanni-Vincentelli, and S. A. Seshia, “Model predictive control with signal temporal logic specifications,” in _53rd IEEE Conference on Decision and Control_. IEEE, 2014, pp. 81–87.
  * [33] R. Koenker and G. Bassett, “Regression quantiles,” _Econometrica_ , vol. 46, no. 1, pp. 33–50, 1978. [Online]. Available: http://www.jstor.org/stable/1913643
  * [34] B. Yu, J. E. Mitchell, and J.-S. Pang, “Solving linear programs with complementarity constraints using branch-and-cut,” _Mathematical Programming Computation_ , vol. 11, no. 2, pp. 267–310, 2019.
  * [35] S. Bradley, A. Hax, and T. Magnanti, _Applied Mathematical Programming_. Addison-Wesley Publishing Company, 1977. [Online]. Available: https://books.google.com/books?id=MSWdWv3Gn5cC
  * [36] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2023. [Online]. Available: https://www.gurobi.com
  * [37] J. Van den Berg, M. Lin, and D. Manocha, “Reciprocal velocity obstacles for real-time multi-agent navigation,” in _2008 IEEE international conference on robotics and automation_. Ieee, 2008, pp. 1928–1935.
  * [38] P. Kothari, S. Kreiss, and A. Alahi, “Human trajectory forecasting in crowds: A deep learning perspective,” _IEEE Transactions on Intelligent Transportation Systems_ , 2021.


# Fiber buckling in confined viscous flows: an absolute instability described
by the Ginzburg-Landau equation

Jean Cappello1, Olivia du Roure1, François Gallaire2, Camille Duprat3, and
Anke Lindner1 1 PMMH, ESPCI, Paris, France
2 LMFI, EPFL, Lausanne, Switzerland
3 LadHyX, École Polytechnique, Palaiseau, France

###### Abstract

We explore the dynamics of a flexible fiber transported by a viscous flow in a
Hele-Shaw cell of height comparable to the fiber height. We show that long
fibers aligned with the flow experience a buckling instability. Competition
between viscous and elastic forces leads to the deformation of the fiber into
a wavy shape convolved by a Bell-shaped envelope. We characterize the
wavelength and phase velocity of the deformation as well as the growth and
spreading of the envelope. Our study of the spatio-temporal evolution of the
deformation reveals a linear and absolute instability arising from a local
mechanism well described by the Ginzburg-Landau equation.

The dynamics of a flexible fiber transported in a flow results from a coupling
between deformation and transport and exhibits rich behaviours du Roure et al.
(2019) including bending Cappello et al. (2019); Marchetti et al. (2018); Li
et al. (2013), buckling Guglielmini et al. (2012); Li et al. (2013); Young and
Shelley (2007); Wandersman et al. (2010) or coiling Chakrabarti et al. (2020).
Such dynamics have been studied in different types of flows including simple
shear flows Liu et al. (2018); Forgacs and Mason (1959); Harasim et al.
(2013); Nguyen and Fauci (2014); Becker and Shelley (2001), vortex arrays
Young and Shelley (2007); Wandersman et al. (2010); Quennouz et al. (2015);
Manikantan and Saintillan (2013), Poiseuille flows Steinhauser et al. (2012),
extensional flows Guglielmini et al. (2012); Manikantan and Saintillan (2015);
Liu et al. (2020); Chakrabarti et al. (2020) or sedimentation in quiescent
fluids Xu and Nadim (1994); Li et al. (2013); Marchetti et al. (2018). The
transport dynamics of flexible fibers in confined geometries has only recently
been explored D’Angelo et al. (2009a, b); Cappello et al. (2019), and shows
even richer dynamics, as a consequence of the interaction between the fiber of
evolving shape, and the confining walls.

Figure 1: (a) Chronophotographies of a parallel fiber transported by the
external flow in the moving reference frame of the fiber. Time increases from
top to bottom (time interval is 4 sec). The phase velocity is denoted
$\vec{V}$ and the horizontal arrow illustrates the direction of the
perturbation propagation. The flow velocity, from left to right, is
$u_{0}=0.65$ mm$\cdot$s-1. The fiber has a length $\ell=10130\pm 15$ $\upmu$m,
width $w=75\pm 3$ $\upmu$m, height $h=73\pm 3$ $\upmu$m, Young’s modulus
$E=27\pm 7$ kPa. The fluid viscosity is $\mu=85\pm 3$ mPa$\cdot$s. Channel
height and width are respectively $H=85\pm 3$ $\upmu$m and $W=3500$ $\upmu$m.
Scale bar is 1000 $\upmu$m. (b) Geometry of the fiber and the confining
channel. The pressure-driven flow is sketched in blue.

In this study, we focus on the dynamics of a flexible fiber transported by a
pressure-driven flow in a Hele-Shaw cell (see Fig. 1). The fiber is aligned
with the flow and is confined by the channel top and bottom walls. Typical
experimental observations can be seen in Fig. 1 (a) showing the shape of a
flexible fiber during transport along the microchannel. The fiber, initially
straight, deforms into a sinusoidal shape with a very well defined wavelength
modulated by a bell-shaped envelope. The dynamics of the deformation is
manifold: the amplitude of the deformation increases with time, the
perturbation travels along the fiber from back to front and the envelope
broadens along the fiber while remaining centered at a fixed lengthwise
position.

Such dynamical shape changes can be triggered through various mechanisms
amplifying initial perturbations, classified as either global or local
instabilities. Global instabilities span the entire domain (here the fiber)
and couple the two end conditions. Hence the boundary conditions are essential
in the selection of the deformation mode. This is typically the case for
classical buckling instabilities where the perturbation ranges over the entire
fibre and the perturbation wavelength is function of the fiber length du Roure
et al. (2019).

On the contrary, in the presence of a local amplifying mechanism the
deformations result from an initial disturbance invading the unstable system
via propagating fronts. If the perturbation grows but is simultaneously
advected such that the disturbance eventually decays at any fixed point, it
ultimately leaves the fiber and the instability is said to be convective. In
contrast, as is the case here, if the perturbation grows at a fixed spatial
position, it eventually invades the entire fiber leading to a self-sustained
instability and is said to be absolute.

Absolute instabilities have been observed in open flow hydrodynamics, such as
the flow past bluff bodies (Bénard-von Kármán vortex street) Triantafyllou et
al. (1987), parallel shear flows Gondret et al. (1999); Charru (2012); Huerre
et al. (2000), transition from dripping to jetting Utada et al. (2008);
Gordillo et al. (2001); Gañán-Calvo et al. (2006); Guillot et al. (2007) or
the flow of viscous films down vertical fibers Duprat et al. (2007); Gallaire
and Brun (2017). However, absolute instabilities resulting in the deformation
of a flexible object in interaction with viscous flows have not been reported
yet. In this letter we describe a one-dimensional example of such an
instability developing during the transport of a flexible fiber in a viscous
flow.

We first characterize the spatial and temporal evolution of the observed wave
by studying the wavelength and the phase velocity of the perturbation as a
function of the physical parameters of the system. Then, by the means of the
linear Ginzburg-Landau model - the archetypal evolution model leading to an
absolute or a convective instability Huerre and Monkewitz (1990) \- we
describe the spatio-temporal dynamics of the envelope and show the linear and
local nature of the absolute instability as well as a near zero group
velocity.

## Material and methods.

Figure 2: Characterization of the fiber deformation. (a) Blue: Superposition
of the shapes of the fiber shown in Fig. 1 in the reference frame of the
fiber. Yellow: Superposition of the envelope of the deformation. Darker colors
indicate later time points (time interval is 0.4 sec). (b) Temporal evolution
of the deformation amplitude at a fixed $x$-location on the fiber (purple
dashed line in (a)). The period of the oscillation, $T$, can be extracted from
this curve. (c) Superposition of shapes normalized by their maximum
($\bar{y}=y/{\rm{max}}(y)$) and aligned on their maxima ($x=0$). The
wavelength of the deformation, $\lambda$, is constant. (d) Growth of the
envelope amplitude with time.

Fibers are fabricated in a microchannel of Hele-Shaw geometry (see Figure 1
(b)), using the stop-flow microscope-based projection photo-lithography method
Dendukuri et al. (2008, 2007). This method enables the fabrication, directly
inside a microchannel, of polymeric PEG-based hydrogel particles of
rectangular cross section, surrounded by a Newtonian solution of uncured
PEGDA. It allows for an excellent control of the particle geometry and its
mechanical properties Duprat et al. (2014); Cappello et al. (2019, 2020).
Details can be found in references Dendukuri et al. (2008); Wexler et al.
(2013); Berthet et al. (2016); Duprat et al. (2014); Nagel et al. (2018);
Cappello et al. (2019); Bechert et al. (2019); Cappello et al. (2020). An
important feature of this method is the existence of a liquid gap between the
fiber and the channel top and bottom walls of constant height $b=6.0\pm 1.6$
$\upmu$m. We use channels of fixed height $H=85\pm 3$ $\upmu$m and different
widths ($W=1000$ $\upmu$m - 3500 $\upmu$m). The resulting fiber height is
$h=73\pm 3$ $\upmu$m leading to strong confinement by the top and bottom
walls: $h/H=0.86$. Fiber width varies from $w=73$ $\upmu$m to $w=118$ $\upmu$m
while the length varies from $\ell=9174\pm 15$ $\upmu$m to $\ell=15000\pm 15$
$\upmu$m and Young modulus varies from $E=117\pm 64$ kPa to $E=27\pm 7$ kPa.

Once a fiber is fabricated, flow is turned on with a syringe pump (Nemesys,
Cetoni). Because of the Hele-Shaw geometry of the channel, the flow is plug-
like in the ($Oxy$) plane, except in the vicinity of the lateral walls and of
the fiber, and Poiseuille like in the ($Oxz$) plane (see Fig. 1 (b)). The
fiber is imaged using an Hamamatsu Orca-flash 4.0 camera, at 10 fps. The fiber
is kept in the camera’s field of view by manually moving the microscope stage.

The fiber shapes are extracted from images by standard image processing (with
ImageJ Schneider et al. (2012) and Matlab). The wavelength is measured by
spatial Fourier transform while the envelope is obtained through an Hilbert
transform.

## Results

The evolution of the shape of the fiber shown in Fig. 1 (a) is summarized in
Fig. 2 (a). The fiber adopts a wavy shape (blue curves) whose amplitude
decreases toward the edges as emphasized by the shape envelopes (plotted in
yellow). Darker curves correspond to later time points and one can see that
the deformation amplitude increases with time. In Fig. 2 (a) the back ends
have been aligned to highlight that the wave propagates along the fiber from
back to front at the phase velocity $V=\lambda/T$ with $\lambda$ the
wavelength and $T$ the wave period. The wave period is measured by following
the deformation amplitude $y$ on a fixed abscissa as shown in Fig. 2 (b). The
wavelength is shown to be constant by superposing normalized shapes aligned on
their maxima (Fig. 2 (c)). As time increases, the deformation increases and
spreads along the fiber (see Fig. 2 (a) and (d) and Figs. 4 (a) and (c)) until
saturation (not shown). The evolution of the fiber velocity is given in Fig. 1
of the supp. data. After a transient regime the fiber velocity reaches a
plateau value that precedes the mechanical instability.

Figure 3: Wavelength characterization. (a) Evolution of the wavelength
$\lambda$ as a function of the fiber length $\ell$ (all other parameters are
kept constant. Fiber: $E=0.20\pm 0.11$ kPa, $w=65\pm 3\upmu$m, $h=73\pm
3\upmu$m. Channel: $W=3500$ $\upmu$m, $H=85$ $\upmu$m.). (b) Evolution of
$\lambda$ as a function of $u_{0}$. Two sets of data are represented gathering
59 experiments with different Young moduli, fiber lengths and channel widths,
other parameter being constant ($w=75\pm 9$ $\upmu$m, $h=73\pm 3$ $\upmu$m,
$H=85\pm 3$ $\upmu$m). Yellow and purple: $E=117\pm 64$ kPa, $\ell=15.64\pm
0.47$mm. Red, dark blue and light blue: $E=27\pm 7$ kPa, $\ell=10.187\pm
0.025$ mm. Red: $W=3640\pm 10$ $\upmu$m; yellow and dark blue: $W=3500\pm 10$
$\upmu$m; purple and light blue: $W=1000\pm 10$ $\upmu$m. Dotted black lines
correspond to fits by a power law. (c) Evolution of the wavelength as a
function of the elasto-viscous length $\Lambda$. The data set is the same as
in (b) (same color code) and has been complemented by the green points
corresponding to variation of the width (from $73$ $\upmu$m to $118$
$\upmu$m), other parameters being constant ($E=27\pm 7$ kPa, $\ell=10.23\pm
0.08$ mm, $h=73\pm 3$ $\upmu$m and $u_{0}=3.53$ mm/s). In total 72
measurements are shown in the figure. The dotted line corresponds to a linear
fit of equation $\lambda$ = 1.33 $\Lambda$.

To study the instability we first focus on the spatial and temporal evolution
of the signal by characterizing the wavelength and the phase velocity of the
wave. We characterize the dependence of the wavelength on the different
parameters of the problem in Fig. 3. Contrary to classical buckling
instabilities, the wavelength does not depend on the fiber length (Fig. 3 (a))
which indicates that the instability results from a local mechanism. However,
perturbations are only observable for fibers one order of magnitude larger
than $\lambda$ (gray region at low $\ell$ in Fig. 3 (a)).

The deformation results from the competition between the elastic restoring
force ($F_{\rm{elast}}$) and a forcing by the hydrodynamic force
($F_{\rm{visc}}$). By balancing these two forces one can build an elasto-
viscous length $\Lambda$ Gosselin et al. (2014); Coq et al. (2008);
Chakrabarti et al. (2020):

$\displaystyle
F_{\rm{elast}}=F_{\rm{visc}}\,\,\Leftrightarrow\,\,EI/\Lambda^{2}=\mu
u_{0}\Lambda\,\,\Leftrightarrow\,\,\Lambda=\left(\frac{F_{elast}}{F_{visc}}\right)^{1/3}$
(1)

with $u_{0}$ the flow velocity, $\mu$ the fluid viscosity, $E$ the Young
modulus of the fiber and $I=\frac{hw^{3}}{12}$ its second moment of area.

To test this scaling we plotted in Fig. 3 (b) the evolution of the wavelength
as a function of the flow velocity, for two different Young moduli, different
channel widths and different fiber lengths. While the wavelength depends on
flow velocity and Young modulus, it does not depend on channel width.
Increasing $u_{0}$ or decreasing $E$ results in a decrease of $\lambda$. The
dependence of the wavelength on the flow velocity is nicely fitted by a power
law $\lambda\propto u_{0}^{\alpha}$, with $\alpha=-0.2\pm 0.1$ independent of
$E$. These results are in fair agreement with the scaling given by equation
(1) which predicts a decrease of the wavelength with the flow velocity as
$\lambda\propto u_{0}^{1/3}$ and an increase with the Young modulus. Fig. 3
(c) gathers all the measurements shown in Fig. 3 (b) enriched with experiments
where we varied the fiber width keeping the other parameters constant (green
points). All wavelengths collapse onto a unique linear curve when plotted as a
function of the elasto-viscous length $\Lambda$, confirming $\Lambda$ to be
the relevant lengthscale of the system.

We also analyze the phase velocity $V$ and show its proportionality to the
mean flow velocity, with channel and fiber width having only a small impact
(see Fig. 2 in the supp. data.).

Figure 4: (a) Fiber shape as a function of time. Experimental conditions are:
$u_{0}=3.53$ mm/s, $E=27\pm 7$ kPa, $w=83\pm 3$ $\upmu$m, $h=73\pm 3$ $\upmu$m
and $\mu=85\pm 3$ mPa.s. (b) A few (x, y) shape extracted from (a)
illustrating the colour-code used in (a). Black curves are Gaussian fits of
the envelope. (c) Same representation as (a) for a Ginzburg-Laudau solution
using parameters adjusted from the experiment shown in (a) and (b) (see text
for details) $A=9.7\pm 1.9\times 10^{4}$ $\upmu$m2/s, $t_{0}=1.3\pm 0.8$ s,
$x_{0}=4.57\pm 0.12$ mm, $U=-31\pm 29$ $\upmu$m/s, $\sigma_{0}=1.37\pm 0.07$
s-1). $v^{+}$ and $v^{-}$ are the velocity of the leading and trailing edges
of the wave ($v^{\pm}=U\pm\sqrt{4A\sigma_{0}}$). (d), (e) and (f) correspond
respectively to the evolution of $\Delta$, $\Phi$ and $X$ as a function of
time. Dotted black curves represent fits with the equations (5-7). Only blue
dots are taken into account for the fits, yellow dots are excluded because
they derive from measurements of small deformations and are poorly reliable.
(g) Diffusion coefficient of the envelope $A$ as a function of $u_{0}\Lambda$
and (h) growth rate $\sigma_{0}$ as a function of $u_{0}/\Lambda$ fitted from
the different experiments shown in Fig.3 (same color code). The dotted black
lines are of slope unity. (i) Histogram of the group velocity $U$. The set of
experiments is the same as in Fig. 3.

To further describe the instability, including the envelop dynamics, we plot
the spatio-temporal evolution of the fiber shape (see Fig. 4 (a)). Here,
$x$-axis is the position along the fiber and the color code gives the
deformation of the fiber in the perpendicular direction while time is
represented on the vertical axis from bottom to top. The perturbation is
initiated at roughly the center of the fiber - region where potential ends
effect are less likely to occur - and grows and spreads as time increases. The
signal can be described as a traveling wave modulated by an amplitude (see
Fig. 4 (b)):

$\displaystyle
y(x,t)=\operatorname{Re}\left\\{Y(x,t)e^{i\frac{2\pi}{\lambda}(x-Vt)}\right\\},$
(2)

with $Y(x,t)$ the amplitude of the envelope, a function of time and position.

On Figs. 4 (a) and (b) we observe that the perturbation eventually invades all
the fiber which indicates of an absolute instability. Consequently we use the
linear and real Ginzburg-Landau (G-L) amplitude equation Godreche and
Manneville (1998) to describe the instability:

$\displaystyle\frac{\partial Y}{\partial t}=\sigma_{0}Y-U\frac{\partial
Y}{\partial x}+A\frac{\partial^{2}Y}{\partial x^{2}},$ (3)

where $\sigma_{0}$, $U$ and $A$, are the growth rate, the group velocity, and
the diffusion coefficient of the envelope.

The general solution of equation (3), for a Dirac perturbation of amplitude
$Y_{0}$ at time $t_{0}$ and position $x_{0}$, is a Gaussian function Castaing
(2005):

$\displaystyle Y(x,t)=\Phi(t)e^{-\frac{(x-X(t))^{2}}{\Delta(t)}},$ (4)

with the three quantities:

$\displaystyle\Delta(t)$ $\displaystyle=4A(t-t_{0}),$ (5)
$\displaystyle\Phi(t)$
$\displaystyle=-{Y_{0}}/\left({2\pi\sqrt{2A(t-t_{0})}}\right)e^{\sigma_{0}(t-t_{0})}$
(6) $\displaystyle\text{ and }\,\,X(t)$ $\displaystyle=x_{0}+U(t-t_{0})$ (7)

that only depend on time.

We then fit the envelope extracted from the experiments by a Gaussian function
for every instant (see Fig. 4 (b)) to obtain $\Delta$, $\Phi$ and $X$ as a
function of time (see Fig. 4 (d), (e) and (f)). The temporal evolutions of
these parameters are in good agreement with the equations (5-7) as shown by
the dotted curves in Fig. 4 (d-f) that result from the fitting of the data.
The values of $A$, $t_{0}$, $\sigma_{0}$, $U$ and $x_{0}$ extracted from the
fitting procedure are then used to build the solution of the Ginzburg-Landau
equation $Y(x,t)$. Fig. 4 (c) shows the resulting spatio-temporal evolution
which is in excellent agreement with the experimental observations (Fig. 4
(a)) confirming that this new instability is indeed local and linear.

To link the ad-hoc description of the system dynamics with the linear G-L
equation to the physics of the underlying problem we determine the evolution
of $A$, $\sigma_{0}$ and $U$ as a function of $u_{0}$, $w$ and $E$. We observe
that the diffusion coefficient $A$ scales reasonably well as $u_{0}\Lambda$
(see Fig. 4 (g)). Fig. 4 (h) shows the evolution of the growth rate
$\sigma_{0}$ as a function of $u_{0}/\Lambda$. All points collapse fairly well
onto a master curve, confirming that the characteristic length of the system
is $\Lambda$ and that all velocities are proportional to $u_{0}$. Lastly, as
visible in Fig. 4 (i), $U$ is always close or equal to zero, i.e. the center
of the envelope only slightly moves along the fibers.

## Conclusion

When sufficiently long i.e. $\ell\gg\Lambda$, a parallel fiber transported by
an external flow in a confining geometry - the fiber height being similar to
the channel height - experiences an instability. The fiber deforms in a wavy
shape modulated by an envelope which spreads along the fiber and whose
amplitude increases with time. We show that the wavelength is constant with
time and does not depend on the fiber length. It is set instead by the elasto-
viscous length $\Lambda$ on which elastic forces balance hydrodynamic forces.
By investigating the spatio-temporal evolution of the signal we confirm that
we witness an absolute instability arising from a local amplification
mechanism. We show that the dynamics of the amplitude and the spreading of the
envelope are well described by the linear Ginzburg-Landau equation and the
overall envelope dynamics is completely captured by three parameters: a group
velocity $U$, a diffusion coefficient $A$ and a growth rate $\sigma_{0}$. In
all our experiments the group velocity is close to zero, meaning that the
envelope only slightly moves along the fibers. On the contrary, the
wavelength, the phase velocity, the growth rate and the diffusion coefficient
all depend on $\Lambda$ and/or on $u_{0}$. All velocities are proportional to
the flow velocity as expected in viscous flows. As $\Lambda$ takes into
account all the physical ingredients of the system we can claim that the
instability is fully described by a competition between two forces: a
destabilising hydrodynamic force and a restoring elastic force. In light of
these conclusions we attempted several approaches to derive a dispersion
relation (see reference Cappello (2020)), however further efforts are
necessary for a full understanding of the detailed mechanisms.

## References

  * du Roure et al. (2019) O. du Roure, A. Lindner, E. Nazockdast, and M. Shelley, Annual Review of Fluid Mechanics (2019).
  * Cappello et al. (2019) J. Cappello, M. Bechert, C. Duprat, O. du Roure, F. Gallaire, and A. Lindner, Physical Review Fluids 4, 034202 (2019).
  * Marchetti et al. (2018) B. Marchetti, V. Raspa, A. Lindner, O. du Roure, L. Bergougnoux, É. Guazzelli, and C. Duprat, to be submitted pp. 1–23 (2018).
  * Li et al. (2013) L. Li, H. Manikantan, D. Saintillan, and S. E. Spagnolie, Journal of Fluid Mechanics 735, 705 (2013), ISSN 0022-1120, eprint 1306.4692, URL http://www.journals.cambridge.org/abstract{_}S0022112013005120.
  * Guglielmini et al. (2012) L. Guglielmini, A. Kushwaha, E. S. Shaqfeh, and H. A. Stone, Physics of Fluids 24, 123601 (2012).
  * Young and Shelley (2007) Y.-N. Young and M. J. Shelley, Physical review letters 99, 058303 (2007).
  * Wandersman et al. (2010) E. Wandersman, N. Quennouz, M. Fermigier, A. Lindner, and O. Du Roure, Soft matter 6, 5715 (2010).
  * Chakrabarti et al. (2020) B. Chakrabarti, Y. Liu, J. LaGrone, R. Cortez, L. Fauci, O. du Roure, D. Saintillan, and A. Lindner, Nature Physics 16, 689 (2020).
  * Liu et al. (2018) Y. Liu, B. Chakrabarti, D. Saintillan, A. Lindner, and O. Du Roure, Proceedings of the National Academy of Sciences 115, 9438 (2018).
  * Forgacs and Mason (1959) O. Forgacs and S. Mason, Journal of Colloid Science 14, 473 (1959).
  * Harasim et al. (2013) M. Harasim, B. Wunderlich, O. Peleg, M. Kröger, and A. R. Bausch, Physical review letters 110, 108302 (2013).
  * Nguyen and Fauci (2014) H. Nguyen and L. Fauci, Journal of The Royal Society Interface 11, 20140314 (2014).
  * Becker and Shelley (2001) L. E. Becker and M. J. Shelley, Physical Review Letters 87, 198301 (2001).
  * Quennouz et al. (2015) N. Quennouz, M. Shelley, O. Du Roure, and A. Lindner, J Fluid Mech 769, 387 (2015).
  * Manikantan and Saintillan (2013) H. Manikantan and D. Saintillan, Physics of Fluids 25, 073603 (2013).
  * Steinhauser et al. (2012) D. Steinhauser, S. Köster, and T. Pfohl, ACS Macro Letters 1, 541 (2012).
  * Manikantan and Saintillan (2015) H. Manikantan and D. Saintillan, Physical Review E 92, 041002 (2015).
  * Liu et al. (2020) Y. Liu, K. Zografos, J. Fidalgo, C. Duchêne, C. Quintard, T. Darnige, F. Vasco, S. Huille, O. du Roure, M. S. Oliveira, et al., Soft Matter (2020).
  * Xu and Nadim (1994) X. Xu and A. Nadim, Physics of Fluids 6, 2889 (1994), ISSN 1070-6631, URL http://scitation.aip.org/content/aip/journal/pof2/6/9/10.1063/1.868116.
  * D’Angelo et al. (2009a) M. V. D’Angelo, B. Semin, G. Picard, M. E. Poitzsch, J. P. Hulin, and H. Auradou, Trans. Porous Med. 84, 389 (2009a), ISSN 0169-3913, URL http://link.springer.com/10.1007/s11242-009-9507-x.
  * D’Angelo et al. (2009b) M. D’Angelo, H. Auradou, G. Picard, M. Poitzsch, and J. Hulin, in _Journal of Physics: Conference Series_ (IOP Publishing, 2009b), vol. 166, p. 012001.
  * Triantafyllou et al. (1987) G. S. Triantafyllou, K. Kupfer, and A. Bers, Physical review letters 59, 1914 (1987).
  * Gondret et al. (1999) P. Gondret, P. Ern, L. Meignin, and M. Rabaud, Physical review letters 82, 1442 (1999).
  * Charru (2012) F. Charru, _Instabilités hydrodynamiques_ (EDP Sciences, 2012).
  * Huerre et al. (2000) P. Huerre, G. Batchelor, H. Moffatt, and M. Worster, Perspectives in fluid dynamics pp. 159–229 (2000).
  * Utada et al. (2008) A. S. Utada, A. Fernandez-Nieves, J. M. Gordillo, and D. A. Weitz, Physical review letters 100, 014502 (2008).
  * Gordillo et al. (2001) J. M. Gordillo, A. M. Gañán-Calvo, and M. Pérez-Saborid, Physics of fluids 13, 3839 (2001).
  * Gañán-Calvo et al. (2006) A. M. Gañán-Calvo et al., Journal of Fluid Mechanics 553, 75 (2006).
  * Guillot et al. (2007) P. Guillot, A. Colin, A. S. Utada, and A. Ajdari, Physical review letters 99, 104502 (2007).
  * Duprat et al. (2007) C. Duprat, C. Ruyer-Quil, S. Kalliadasis, and F. Giorgiutti-Dauphiné, Physical review letters 98, 244502 (2007).
  * Gallaire and Brun (2017) F. Gallaire and P.-T. Brun, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375, 20160155 (2017).
  * Huerre and Monkewitz (1990) P. Huerre and P. A. Monkewitz, Annual review of fluid mechanics 22, 473 (1990).
  * Dendukuri et al. (2008) D. Dendukuri, P. Panda, R. Haghgooie, J. M. Kim, T. A. Hatton, and P. S. Doyle, Macromolecules 41, 8547 (2008), ISSN 00249297.
  * Dendukuri et al. (2007) D. Dendukuri, S. S. Gu, D. C. Pregibon, T. A. Hatton, and P. S. Doyle, Lab on a Chip 7, 818 (2007), ISSN 1473-0197, URL http://xlink.rsc.org/?DOI=b703457a.
  * Duprat et al. (2014) C. Duprat, H. Berthet, J. S. Wexler, O. du Roure, and A. Lindner, Lab on a chip 15, 244 (2014), ISSN 1473-0189, URL http://www.ncbi.nlm.nih.gov/pubmed/25360871.
  * Cappello et al. (2020) J. Cappello, V. D’herbemont, A. Lindner, and O. Du Roure, Micromachines 11, 318 (2020).
  * Wexler et al. (2013) J. S. Wexler, P. H. Trinh, H. Berthet, N. Quennouz, O. du Roure, H. E. Huppert, A. Linder, and H. A. Stone, Journal of Fluid Mechanics 720, 517 (2013), ISSN 0022-1120, URL http://www.journals.cambridge.org/abstract{_}S0022112013000499.
  * Berthet et al. (2016) H. Berthet, O. du Roure, and A. Lindner, Applied Sciences 6, 385 (2016), ISSN 2076-3417, URL http://www.mdpi.com/2076-3417/6/12/385.
  * Nagel et al. (2018) M. Nagel, P.-T. Brun, H. Berthet, A. Lindner, F. Gallaire, and C. Duprat, Journal of Fluid Mechanics 835, 444 (2018).
  * Bechert et al. (2019) M. Bechert, J. Cappello, M. Daïeff, F. Gallaire, A. Lindner, and C. Duprat, EPL (Europhysics Letters) 126, 44001 (2019).
  * Schneider et al. (2012) C. A. Schneider, W. S. Rasband, and K. W. Eliceiri, Nature Methods 9, 671 (2012), ISSN 1548-7091, URL http://dx.doi.org/10.1038/nmeth.2089.
  * Gosselin et al. (2014) F. Gosselin, P. Neetzow, and M. Paak, Physical Review E 90, 052718 (2014).
  * Coq et al. (2008) N. Coq, O. Du Roure, J. Marthelot, D. Bartolo, and M. Fermigier, Physics of Fluids 20, 051703 (2008).
  * Godreche and Manneville (1998) C. Godreche and P. Manneville, Oceanographic Literature Review 9, 1723 (1998).
  * Castaing (2005) B. Castaing, _Hydrodynamics and nonlinear instabilities_ , vol. 3 (Cambridge University Press, 2005).
  * Cappello (2020) J. Cappello, Phd thesis, Université de Paris (2020), URL https://tel.archives-ouvertes.fr/tel-03260710.