# Two-stage robust optimization for orienteering problem with stochastic
weights

Ke Shang  Felix T.S. Chan  Stephen Karungaru<EMAIL_ADDRESS>Kenji Terada  Zuren Feng  Liangjun Ke  Department of Information Science and
Intelligent Systems, The University of Tokushima, Japan  State Key Laboratory
for Manufacturing Systems Engineering, Xi’an Jiaotong University, Xi’an, China
Department of Industrial and Systems Engineering, The Hong Kong Polytechnic
University, Hung Hom, Hong Kong

###### Abstract

In this paper, the two-stage orienteering problem with stochastic weights
(OPSW) is considered, where the first-stage problem is to plan a path under
the uncertain environment and the second-stage problem is recourse action to
make sure that the length constraint is satisfied after the uncertainty is
realized. Two recourse models are introduced based on two different
uncertainty realization ways, one is based on sequentially realized weights
which leads to the recourse model proposed by Evers et al. (2014) and the
other is based on concurrently realized weights which leads to a new recourse
model with less variables and less constraints and is computationally more
efficient. Subsequently two two-stage robust models are introduced for OPSW
based on the two different recourse models, and the relationships between the
two-stage robust models and their corresponding static robsut models are
investigated. Theoretical conclusions are drawn which show that the two-stage
robust models are equivalent to their corresponding static robust models with
the box uncertainty set defined, and the two two-stage robust models are also
equivalent to each other even though they are based on different recourse
models. A case study is presented by comparing the two-stage robust models
with an one-stage robust model for OPSW. The numerical results of the
comparative studies show the effectiveness and superiority of the proposed
two-stage robust models for dealing with the two-stage OPSW.

###### keywords:

two-stage robust optimization , stochastic orienteering problem , integer
recourse , box uncertainty set

††journal:

##  1 Introduction

The orienteering problem (OP) is a routing problem which has been widely
studied in the past few decades. It was first introduced by Golden et al.
(1987) and has been developed in terms of the problem variants, solution
algorithms and applications. The original OP aims at planning a path which
starts and ends at the depot location, and visits a subset of nodes in order
to maximize the total collected score while the length of the path cannot
exceed a predefined budget, and each node can only be visited at most one
time. The OP has a wide practical application background. A few examples such
as Unmanned Aerial Vehicle (UAV) mission planning (Mufalli et al. (2012);
Evers et al. (2014)), tourist trip design problem (Vansteenwegen and Van
Oudheusden (2007); Gavalas et al. (2014)) and mobile crowdsourcing problem
(Howe (2008); Yuen et al. (2011)). More detailed surveys on OP are given in
Vansteenwegen et al. (2011); Gunawan et al. (2016).

The stochastic orienteering problem (SOP) is a variant of OP, which assumes
that some parameters in OP are stochastic and uncertain such as the score
associated with each node and the weight (or distance) associated with each
arc. SOP is more appropriate than OP in practical situations. For example, in
a practical environment traffic congestion affects the travel time between
nodes. Ilhan et al. (2008) first considered uncertainties in the score of
nodes and the resulting SOP is called OP with stochastic profits (OPSP).
Campbell et al. (2011); Evers et al. (2014) considered uncertainties in the
travel and service times and the resulting SOP is called OP with stochastic
travel and service times (OPSTS) or OP with stochastic weights (OPSW). Other
variants include the dynamic stochastic OP (DSOP) with stochastic time-
dependent travel times (Lau et al. (2012); Varakantham and Kumar (2013)) and
the stochastic OPTW (SOPTW) with stochastic waiting time (Zhang et al.
(2014)).

In this paper, we focus on the OPSW where the uncertainties lie in the weights
of the arcs. Some works on OPSW have been done in recent years. Campbell et
al. (2011) considered the OPSTS in which a penalty is incurred if a commitment
to a node is made but not completed. A variant of VNS for the OPSTS is
proposed and its performance is evaluated by comparing with a dynamic
programming (DP) approach. Evers et al. (2014) introduced a two-stage
stochastic programming model for the OPSW. The first-stage problem is to plan
a path. The second-stage problem is a recourse action after the uncertain
weights realized, which aborts the execution of the first-stage path and
enforces a direct return to the depot. The objective is to maximize the first-
stage path score plus the expected loss of the score due to the recourse
action. They presented a Sample Average Approximation (SAA) approach and an
OPSW heuristic to solve the problem and the performance of the two approaches
were evaluated. Evers et al. (2014) applied the robust optimization (RO)
methodology to build robust models for UAV mission planning with uncertain
fuel usage between targets, which is an OPSW in nature. The performance of the
robust models are studied in terms of the different uncertainty sets and the
feasibility of the robust solutions.

Inspired by the recourse model proposed by Evers et al. (2014), we consider
the two-stage OPSW, i.e. OPSW with recourse action. The first-stage problem is
to plan a path with the stochastic weights unrevealed. The second-stage
problem is a recourse action to avoid the violation of the length budget after
the uncertainty realized. The recourse action is to abort the execution of the
first-stage path and enforce a direct return to the depot. This kind of
recourse action is necessary especially in the UAV mission planning. The UAV
has to return to the depot safely in the uncertain environment. We notice that
the realization way of the uncertainty is not unique. For example, the
uncertain weights of the first-stage path can be realized sequentially during
the path execution, or the uncertain weights of the first-stage path can be
realized concurrently at the beginning of the path execution. Different
realization ways for the uncertainty can lead to different recourse models. We
thus define two realization ways of the uncertainty: _Sequential realization_
and _Concurrent realization_. The _Sequential realization_ way leads to the
recourse model proposed by Evers et al. (2014), and the _Concurrent
realization_ way leads to a new model with less variables and less
constraints, which is computationally more attractive.

Two-stage robust optimization (RO), also known as adjustable RO and can be
extended to the multi-stage situation, was initially introduced by Ben-Tal et
al. (2004). Compared with traditional one-stage RO, two-stage RO divides the
decision variables into “here and now” decisions and “wait and see” decisions,
which is more flexible and is suitable for modeling two-stage problems. It has
been successfully applied to different applications such as unit commitment
(An and Zeng (2015); Wang et al. (2016)), network flow (Atamtürk and Zhang
(2007); Ordóñez and Zhao (2007)) and portfolio optimization (Takeda et al.
(2008)). In this paper, we apply the two-stage RO paradigm to the two-stage
OPSW for the first time and introduce two two-stage RO models based on two
different recourse models. The two-stage RO models introduced in this paper
are with binary recourse decisions and this kind of problem has largely
resisted solution so far (Hanasusanto et al. (2015)). Computing an optimal
adjustable robust solution is often intractable since it requires computing a
solution for all possible realizations of the uncertainties (Feige et al.
(2007)). Instead of solving the two-stage RO model directly, Bertsimas et al.
(2015) studied the performance of the static solutions for two-stage
adjustable robust linear optimization problems and presented a tight
characterization of the conditions under which a static robust solution is
optimal for the two-stage robust problem. From this point of view, we
introduce two static robust models for the OPSW which correspond to the two
two-stage robust models respectively, and study their performance and the
relationships with the two-stage robust models. We prove that with the _box
uncertainty set_ defined, the two-stage robust models are equivalent to their
corresponding static robust models, and the two two-stage robust models are
also equivalent to each other even though they are based on different recourse
models. These conclusions we obtained indicate that the two-stage robust
models for OPSW can be solved to optimality by solving their corresponding
static robust models, and also we only need to use one static robust model,
which is based on the second recourse model, to deal with two different
uncertainty realization ways.

The main contributions of this paper are summarized as follows:

  1. 1.

Two recourse models are presented for the two-stage OPSW: one is the recourse
model with _Sequential realization_ and the other is the recourse model with
_Concurrent realization_.

  2. 2.

Two two-stage robust models are presented for the first time for the OPSW
based on the two different recourse models.

  3. 3.

Three theorems are established which show the equivalence between the two-
stage robust models and their corresponding static robust models.

  4. 4.

The two-stage robsut models for OPSW are evaluated numerically by comparing
with one-stage robust model for OPSW.

The remainder of the paper is organized as follows. First the deterministic OP
is described in Section 2. Section 3 describes the two-stage OPSW and
introduces two recourse models with different uncertainty realization ways.
Section 4 introduces two two-stage robust models for OPSW and draws some
theoretical conclusions of the equivalence between the two-stage robust models
and their corresponding static robust models. A case study is presented in
Section 5 and we conclude the whole paper in Section 6.

##  2 The deterministic OP

In a deterministic OP, a set of vertices $N$ is given with $|N|$ as its
cardinality. Each vertex $i\in N$ has a score $s_{i}$ associated with it.
Denote $0$ as the depot location where $0\notin N$ and $N^{+}=N\cup\\{0\\}$.
The goal is to plan a path with length limit $L$, that starts and ends at the
depot and visits some vertices in order to maximize the sum of the collected
scores. Each vertex is visited at most one time.

Suppose all nodes $N^{+}$ are on a complete graph $G=(N^{+},A)$ where $A$ is
the set of arcs connecting the vertices in $N^{+}$. The weight of each arc
$(i,j)\in A$ is $d_{ij}$, representing the Euclidean distance from $i$ to $j$.
Let $x_{ij}$ be a binary decision variable, where $x_{ij}=1$ if and only if
arc $(i,j)$ is visited by the path, otherwise $x_{ij}=0$. An auxiliary
variable $u_{i}$ is used to denote the position of node $i$ in the path. The
formulation of the deterministic OP is as follows:

(DOP): maximize $\displaystyle\sum_{i\in N}s_{i}\sum_{j\in
N^{+}\setminus\\{i\\}}x_{ij}$ (1a) subject to $\displaystyle\sum_{(i,j)\in
A}d_{ij}x_{ij}\leq L$ (1b) $\displaystyle\sum_{i\in N}x_{0i}=\sum_{i\in
N}x_{i0}=1$ (1c) $\displaystyle\sum_{i\in
N^{+}\setminus\\{j\\}}x_{ij}=\sum_{i\in N^{+}\setminus\\{j\\}}x_{ji}\leq
1,\forall j\in N$ (1d) $\displaystyle u_{i}-u_{j}+1\leq(1-x_{ij})|N|,\forall
i,j\in N$ (1e) $\displaystyle 1\leq u_{i}\leq|N|,\forall i\in N$ (1f)
$\displaystyle x_{ij}\in\\{0,1\\},\forall(i,j)\in A$ (1g)

Constraint (1b) is the path length constraint. Constraint (1c) guarantees that
the path starts and ends at the depot. Constraint (1d) is the flow
conservation constraint ensuring that a vertex is visited at most once.
Constraint (1e) ensures the connectivity of the path. Constraint (1f) and (1g)
are the boundary and integrality constraints on the auxiliary variables and
decision variables respectively.

##  3 The two-stage orienteering problem with stochastic weights

Suppose the weight of each arc $(i,j)$ is stochastic and uncertain, denote the
stochastic weight of arc $(i,j)$ as $\tilde{d}_{ij}$. In this paper, we
consider $\tilde{d}_{ij}$ as a symmetrically distributed random variable on
the interval $[\bar{d}_{ij}-\hat{d}_{ij},\bar{d}_{ij}+\hat{d}_{ij}]$, where
$\bar{d}_{ij}$ is the expected value of $\tilde{d}_{ij}$ and $\hat{d}_{ij}$ is
the maximum deviation of $\tilde{d}_{ij}$ from its expected value. For
simplicity and convenience, we use $d_{ij}$ to denote the realizations of
$\tilde{d}_{ij}$.

We consider the two-stage OPSW, i.e. OPSW with recourse action. In the two-
stage OPSW, the first-stage problem is to plan a path with the stochastic
weights unrevealed. Due to the randomicity and uncertainty of the stochastic
weights, the first-stage path may violate constraint (1b) after the
uncertainty is realized. So the second-stage problem is a recourse action to
avoid constraint violation after the uncertainty is realized. The recourse
action is to abort the execution of the first-stage path and enforce a direct
return to the depot.

The ways that the uncertainty realized are not unique. Different realization
ways for the uncertainty will lead to different recourse models. We introduce
two realization ways for the uncertainty in two-stage OPSW: _Sequential
realization_ and _Concurrent realization_. _Sequential realization_ means that
the stochastic weights of the first-stage path are realized sequentially
during the path execution. For example, the first-stage path is executed to
node $i$ and the next node is $j$, then the stochastic weight $\tilde{d}_{ij}$
is realized and the stochastic weights of all other unvisited arcs remain
unrevealed. _Concurrent realization_ means that all the stochastic weights of
the first-stage path are realized concurrently at the beginning of the path
execution and the stochastic weights of all other arcs remain unrevealed.

Based on the above two realization ways for uncertainty, we now present two
recourse models for two-stage OPSW.

###  3.1 Recourse model with Sequential realization

The recourse model with sequentially realized weights was initially introduced
by Evers et al. (2014). In this model, the uncertain weights of the first-
stage path are realized sequentially during the path execution. The
uncertainty realization rule is: suppose the first-stage path is executed to
node $i$ and the next node is $j$, then the stochastic weight $\tilde{d}_{ij}$
is realized and the stochastic weights of all other unvisited arcs remain
unrevealed. Then the recourse action is to abort the execution of the first-
stage path and enforce a direct return to the depot from node $i$ at the
moment that the remaining length budget is insufficient to support a visit to
the next node $j$ plus the expected return length from the next node $j$ to
the depot. Evers et al. (2014) assumed that a certain amount of extra length
budget is available to cover the maximum deviation from the expected length on
any of the arcs to the depot, this safety stock not being part of the length
limit $L$ used in the model.

Denoting the first-stage path as vector $\mathbf{x}$ which contains all
$x_{ij}$, and the weight realizations as vector $\mathbf{d}$ which contains
all $d_{ij}$. Let $x_{ijk}$ be a binary variable, $x_{ijk}=1$ if arc $(i,j)$
is the $k$th arc in the first-stage path, otherwise $x_{ijk}=0$; let $y_{i}$
be a binary variable, $y_{i}=1$ if node $i$ is in the first-stage path but
cannot be reached as a result of the recourse action, otherwise $y_{i}=0$; let
$z_{k}$ be a binary variable, $z_{k}=1$ if the $k$th node in the first-stage
path cannot be reached as a result of the recourse action, otherwise
$z_{k}=0$. With the first-stage path $\mathbf{x}$ and the weight realizations
$\mathbf{d}$, the recourse problem of the two-stage OPSW with sequentially
realized weights is formulated as follows:

(Recourse-Sequential): $\displaystyle RS(\mathbf{x},\mathbf{d})=\max\qquad$
$\displaystyle-\sum_{i\in N}s_{i}y_{i}$ (2a) subject to $\displaystyle
x_{0j1}\geq x_{0j},\forall j\in N$ (2b) $\displaystyle x_{ijk}\geq
x_{ij}+\sum_{l\in N^{+}}x_{li(k-1)}-1,\forall i,j\in N,k=1,...,|N|$ (2c)
$\displaystyle\sum_{k=1}^{K}\sum_{(i,j)\in A}d_{ij}x_{ijk}+\sum_{(i,j)\in
A}\bar{d}_{j0}x_{ijK}\leq L+Mz_{K},\forall K=1,...,|N|$ (2d) $\displaystyle
z_{k}\geq z_{k-1},\forall k=2,...,|N|$ (2e) $\displaystyle y_{j}\geq\sum_{i\in
N^{+}}x_{ijk}+z_{k}-1,\forall j\in N,k=1,...,|N|$ (2f) $\displaystyle
x_{ijk}\in\\{0,1\\},\forall i,j\in N^{+},k=1,...,|N|$ (2g) $\displaystyle
y_{i}\in\\{0,1\\},\forall i\in N$ (2h) $\displaystyle
z_{k}\in\\{0,1\\},\forall k=1,...,|N|$ (2i)

where the objective function (2a) is to minimize the loss in the collected
score as a result of the recourse action. Constraint (2b) identifies the first
arc in the path. Constraint (2c) identifies the order of the other arcs in the
path. Constraint (2d) determines the nodes of the first-stage path
$\mathbf{x}$ that can and cannot be reached based on $\mathbf{d}$, where $M$
is a sufficiently large number. Constraint (2e) makes sure that all nodes in
the path after the first node that cannot be reached, cannot be reached
either. Constraint (2f) identifies the nodes in the first-stage path that
cannot be reached, based on their indexes. A detailed explanation of the model
is given in Evers et al. (2014).

With the first-stage path $\mathbf{x}$ and the weight realizations
$\mathbf{d}$, the objective value of the sequential recourse problem can be
calculated not only by solving the Recourse-Sequential model, but also by an
efficient Forward Checking algorithm which is described in Algorithm 1. With
Algorithm 1, the objective value of the Recourse-Sequential model can be
obtained in time $\mathcal{O}(n)$ where $n$ is the number of nodes in the
first-stage path.

Algorithm 1 Forward Checking algorithm for Recourse-Sequential model

1:The first-stage path $\mathbf{x}$ and the weight realizations $\mathbf{d}$

2:The objective value of $RS(\mathbf{x},\mathbf{d})$

3:Denote the first-stage path $\mathbf{x}$ as a node sequence
$(v_{0},v_{1},...,v_{n},v_{0})$ where $v_{0}$ is the depot and $v_{k}$ is the
$k$th node.

4:$violation=False$

5:for $k=1$ to $n$ do

6: $Length=dist(v_{0},v_{1},...,v_{k})+\bar{d}_{v_{k}v_{0}}$ where
$dist(v_{0},v_{1},...,v_{k})=d_{v_{0}v_{1}}+...+d_{v_{k-1}v_{k}}$.

7: if $Length>L$ then

8: $violation=True$

9: break for loop

10: end if

11:end for

12:if $violation==True$ then

13: $Loss=s_{v_{k}}+s_{v_{k+1}}+...+s_{v_{n}}$

14:else

15: $Loss=0$

16:end if

17:$RS(\mathbf{x},\mathbf{d})=-Loss$

###  3.2 Recourse model with Concurrent realization

We now introduce a recourse model with concurrently realized weights. In this
model, all the stochastic weights of the first-stage path are realized
concurrently at the beginning of the path execution and the stochastic weights
of all other arcs remain unrevealed. Then the recourse action is to find a
node $i$ in the first-stage path and enforce a direct return to the depot from
node $i$ so that the length of the subpath from the depot to node $i$ plus the
expected return length from node $i$ to the depot is within the length limit
$L$ and the loss in the collected score is minimized. Here, we also assume
that a certain amount of extra length budget is available to cover the maximum
deviation from the expected length on any of the arcs to the depot. This
safety stock is not part of the length limit $L$ used in the model.

We introduce a new binary variable $y_{ij}$; $y_{ij}=1$ if arc $(i,j)$ is in
the first-stage path but is cancelled by the recourse action, $y_{ij}=0$ if
arc $(i,j)$ is in the first-stage path and is not cancelled by the recourse
action, or arc $(i,j)$ is not in the first-stage path. Then the recourse
problem of the two-stage OPSW with concurrently realized weights can be
formulated as follows:

(Recourse-Concurrent): $\displaystyle RC(\mathbf{x},\mathbf{d})=\max\qquad$
$\displaystyle-\sum_{j\in N}s_{j}\sum_{i\in N^{+}\setminus\\{j\\}}y_{ij}$ (3a)
subject to $\displaystyle y_{ij}\leq x_{ij},\forall(i,j)\in A$ (3b)
$\displaystyle\sum_{i\in N^{+}\setminus\\{j\\}}y_{ij}\leq\sum_{k\in
N^{+}\setminus\\{j\\}}y_{jk},\forall j\in N$ (3c) $\displaystyle\sum_{(i,j)\in
A}d_{ij}x_{ij}-\sum_{(i,j)\in A}d_{ij}y_{ij}+\sum_{j\in N}\left(\sum_{k\in
N^{+}\setminus\\{j\\}}y_{jk}-\sum_{i\in
N^{+}\setminus\\{j\\}}y_{ij}\right)\bar{d}_{j0}\leq L$ (3d) $\displaystyle
y_{ij}\in\\{0,1\\},\forall(i,j)\in A$ (3e)

where the objective function (3a) is to minimize the loss in the collected
score as a result of the recourse action. Constraint (3b) ensures that the
cancelled arcs are from the first-stage path. Constraint (3c) ensures that the
cancelled arcs compose a subpath of the first-stage path starting at a vertex
of the first-stage path and ending at the depot. Constraint (3d) ensures that
the modified path after the recourse action is within the length limit.

With the first-stage path $\mathbf{x}$ and the weight realizations
$\mathbf{d}$, the objective value of the concurrent recourse problem can be
calculated not only by solving the Recourse-Concurrent model, but also by an
efficient Backward Checking algorithm which is described in Algorithm 2. With
Algorithm 2, the objective value of the Recourse-Concurrent model can be
obtained in time $\mathcal{O}(n)$ where $n$ is the number of nodes in the
first-stage path.

Algorithm 2 Backward Checking algorithm for Recourse-Concurrent model

1:The first-stage path $\mathbf{x}$ and the weight realizations $\mathbf{d}$

2:The objective value of $RC(\mathbf{x},\mathbf{d})$

3:Denote the first-stage path $\mathbf{x}$ as a node sequence
$(v_{0},v_{1},...,v_{n},v_{0})$ where $v_{0}$ is the depot and $v_{k}$ is the
$k$th node.

4:$violation=True$

5:for $k=n$ to $1$ do

6: $Length=dist(v_{0},v_{1},...,v_{k})+\bar{d}_{v_{k}v_{0}}$ where
$dist(v_{0},v_{1},...,v_{k})=d_{v_{0}v_{1}}+...+d_{v_{k-1}v_{k}}$.

7: if $Length\leq L$ then

8: $violation=False$

9: break for loop

10: end if

11:end for

12:if $k\neq n$ and $violation==False$ then

13: $Loss=s_{v_{k+1}}+...+s_{v_{n}}$

14:else if $k\neq n$ and $violation==True$ then

15: $Loss=s_{v_{1}}+...+s_{v_{n}}$

16:else

17: $Loss=0$

18:end if

19:$RC(\mathbf{x},\mathbf{d})=-Loss$

##  4 Two-stage robust optimization for OPSW

In this section, we apply the two-stage RO methodology to model the two-stage
OPSW. In the two-stage OPSW, the first-stage “here and now” decisions are the
binary decision variables $x_{ij}$ described in Section 2. The second-stage
“wait and see” decisions are the binary decision variables $y_{i}$ and $z_{k}$
in the Recourse-Sequential model or the binary decision variables $y_{ij}$ in
the Recourse-Concurrent model.

In the two-stage RO for OPSW, an uncertainty set needs to be defined for the
stochastic weights. We consider the _box uncertainty set_ which is defined by
the $\infty$-norm of the uncertain vector. The reasons that we choose the _box
uncertainty set_ for the two-stage RO for OPSW are as follows:

  1. 1.

It is simple compared with the polyhedral uncertainty set and the ellipsoidal
uncertainty set which is defined by 1-norm and 2-norm respectively (Bertsimas
et al. (2004); Ben-Tal et al. (2009)), and the derived robust counterpart has
the same computational complexity as the original model.

  2. 2.

With the _box uncertainty set_ , we can draw some interesting conclusions in
the following subsections which describe the equivalence between the two-stage
robust models and their corresponding static robust models.

Without loss of generality, the _box uncertainty set_ $\mathcal{U}$ for the
stochastic weights is defined as follows:

$\mathcal{U}=\left\\{\mathbf{d}\in\mathbb{R}^{M}:d_{ij}=\bar{d}_{ij}+\zeta_{ij}\hat{d}_{ij},\forall
i,j\in N^{+},\zeta\in\mathcal{Z}\right\\}$ (4)

where $\mathbf{d}$ is a $M$ dimensional vector with $M=|N^{+}|\times|N^{+}|$,
$\zeta\in\mathbb{R}^{M}$ is the vector of primitive uncertainties, and
$\mathcal{Z}$ is a convex set which is defined as follows:

$\mathcal{Z}=\left\\{\zeta\in\mathbb{R}^{M}:\left\|\zeta\right\|_{\infty}\leq\Theta\right\\}$
(5)

where $\Theta\in[0,1]$ is the parameter controlling the size of $\mathcal{Z}$.

Next, we introduce two two-stage RO models for OPSW based on the Recourse-
Sequential model and the Recourse-Concurrent model respectively.

###  4.1 Two-stage robust model for OPSW with Recourse-Sequential model

Based on the Recourse-Sequential model and the two-stage RO paradigm, we
introduce the following two-stage RO model for OPSW:

(Two-stage-Sequential): maximize $\displaystyle\sum_{i\in N}s_{i}\sum_{j\in
N^{+}\setminus\\{i\\}}x_{ij}+\min_{\mathbf{d}\in\mathcal{U}}RS(\mathbf{x},\mathbf{d})$
(6a) subject to $\displaystyle\sum_{(i,j)\in
A}(\bar{d}_{ij}-\hat{d}_{ij})x_{ij}\leq L$ (6b)
$\displaystyle\eqref{opc:2}-\eqref{opc:6}$ (6c)

where $RS(\mathbf{x},\mathbf{d})$ is the Recourse-Sequential model and
$\mathcal{U}$ is the box uncertainty set. Constraint (6b) is the length limit
on the first-stage path. Without constraint (6b), the first-stage path can be
arbitrarily long providing there exists unvisited nodes, and these nodes can
be included in the first-stage path even some nodes in any particular case
cannot be reached. By adding constraint (6b), we limit the length of the
first-stage path in the most optimistic situation, i.e. all arc weights equal
to their minimum values. With this constraint, the size of the solution space
can be reduced while the problem optimality is maintained.

The two-stage robust model for OPSW introduced above is an 0-1 integer
programming problem with 0-1 integer recourse. Next, we present its
corresponding static robust model in which the second-stage “wait and see”
decisions become “here and now”. The corresponding static robust model of the
Two-stage-Sequential model is formulated as follows:

(Static-Sequential): maximize $\displaystyle\sum_{i\in N}s_{i}\sum_{j\in
N^{+}\setminus\\{i\\}}x_{ij}-\sum_{i\in N}s_{i}y_{i}$ (7a) subject to
$\displaystyle\eqref{opc:2}-\eqref{opc:6},\eqref{2ro:2}$ (7b)
$\displaystyle\eqref{2rm1:1}-\eqref{2rm1:2},\eqref{2rm1:4}-\eqref{2rm1:8}$
(7c) $\displaystyle\sum_{k=1}^{K}\sum_{(i,j)\in A}d_{ij}x_{ijk}+\sum_{(i,j)\in
A}\bar{d}_{j0}x_{ijK}\leq L+Mz_{K},\forall
K=1,...,|N|,\mathbf{d}\in\mathcal{U}$ (7d)

In the above static robust model, the second-stage decision variables $y_{i}$
and $z_{k}$ are “here and now” and do not depend on the realizations of
uncertain $\mathbf{d}$. Both first-stage decisions $x_{ij}$ and second-stage
decisions $y_{i}$ and $z_{k}$ are selected before the uncertain $\mathbf{d}$
is known. An optimal static robust solution to Static-Sequential can be
computed efficiently with the box uncertainty set $\mathcal{U}$. What
interests us is the relationship between the two-stage robust model Two-stage-
Sequential and its corresponding static robust model Static-Sequential.

In the following, we establish a theorem which describes the equivalence of
the Two-stage-Sequential model and the Static-Sequential model.

######  Theorem 1.

The two-stage robust model Two-stage-Sequential and its corresponding static
robust model Static-Sequential are equivalent.

###### Proof.

It is clear that the optimal solution of the static robust model Static-
Sequential is a feasible solution of the two-stage robust model Two-stage-
Sequential. All we need to show is that the optimal solution of the two-stage
robust model Two-stage-Sequential is a feasible solution of the static robust
model Static-Sequential.

We prove by apagoge. Denote $(x_{ij}^{*},x_{ijk}^{*},y_{i}^{*},z_{k}^{*})$ as
the optimal solution of the two-stage robust model Two-stage-Sequential.
Suppose the optimal solution is infeasible for the static robust model Static-
Sequential, which means

$\exists\mathbf{d}^{\prime}\in\mathcal{U}\text{ and
}K^{\prime},\sum_{k=1}^{K^{\prime}}\sum_{(i,j)\in
A}d^{\prime}_{ij}x_{ijk}^{*}+\sum_{(i,j)\in
A}\bar{d}_{j0}x_{ijK^{\prime}}^{*}>L+Mz_{K^{\prime}}^{*}$ (8)

Because $M$ is a sufficiently large number, so the above condition is only
satisfied by $z_{K^{\prime}}^{*}=0$. This means that the $K^{\prime}$th node
in the first-stage path of the two-stage robust model Two-stage-Sequential is
reachable, but this node is unreachable under the context of the static robust
model Static-Sequential.

Now we consider the second-stage problem
$RS(\mathbf{x}^{*},\mathbf{d}^{\prime})$ where $\mathbf{x}^{*}$ is the first-
stage optimal solution, and denote the optimal solution as
$(y^{\prime}_{i},z^{\prime}_{k})$. Then for $K^{\prime}$,

$\sum_{k=1}^{K^{\prime}}\sum_{(i,j)\in
A}d^{\prime}_{ij}x_{ijk}^{*}+\sum_{(i,j)\in
A}\bar{d}_{j0}x_{ijK^{\prime}}^{*}\leq L+Mz^{\prime}_{K^{\prime}}$ (9)

Comparing constraints (8) and (9), it is clear that $z^{\prime}_{K^{\prime}}$
must be 1. This means the $K^{\prime}$th node in the first-stage path is
unreachable with $\mathbf{d}^{\prime}$. Denote the $K^{\prime}$th node in the
first-stage path as node $j$, then the second-stage optimal value
$RS(\mathbf{x}^{*},\mathbf{d}^{\prime})\leq-\sum_{i\in
N}s_{i}y_{i}^{*}-s_{j}$. Because we are optimizing
$\min_{\mathbf{d}\in\mathcal{U}}RS(\mathbf{x}^{*},\mathbf{d})$,
$(y_{i}^{*},z_{k}^{*})$ is not the optimal second-stage solution, so this is a
contradiction. Thus the hypothesis cannot be established, which means
$(x_{ij}^{*},x_{ijk}^{*},y_{i}^{*},z_{k}^{*})$ is feasible for the static
robust model Static-Sequential.

We conclude that the optimal solution of the two-stage robust model Two-stage-
Sequential is a feasible solution of the static robust model Static-
Sequential. Because the two models have the same objective function value with
the same solutions, then the optimal solution of the two-stage robust model
Two-stage-Sequential is also the optimal solution of the static robust model
Static-Sequential, this implies the two models are equivalent. ∎

######  Remark.

The proof of Theorem 1 does not need the support of the _box uncertainty set_.
We can still draw this conclusion even if the uncertainty set $\mathcal{U}$ is
an arbitrary uncertainty set.

Based on Theorem 1, the two-stage robust model Two-stage-Sequential can be
solved to optimality by solving its corresponding static robust model Static-
Sequential. Comparing the static robust model Static-Sequential with the
original deterministic OP model DOP, many new integer variables are added
which makes the static robust model Static-Sequential computationally
expensive. Evers et al. (2014) proved that the relaxation model with $0\leq
x_{ijk}\leq 1$ and $0\leq y_{i}\leq 1$ of the second-stage problem
$RS(\mathbf{x},\mathbf{d})$ is equivalent to the original
$RS(\mathbf{x},\mathbf{d})$, and the resulting relaxation model provides a
substantial decrease in the computation time. This conclusion can be easily
applied to the static robust model Static-Sequential which leads to the
following proposition:

######  Proposition 1.

The relaxation model with $0\leq x_{ijk}\leq 1$ and $0\leq y_{i}\leq 1$ of the
static robust model Static-Sequential is equivalent to the original static
robust model Static-Sequential.

###### Proof.

This conclusion can be drawn by following the proof way of Theorem 1 in Evers
et al. (2014). ∎

###  4.2 Two-stage robust model for OPSW with Recourse-Concurrent model

Based on the Recourse-Concurrent model and the two-stage RO paradigm, we
introduce the following two-stage RO model for OPSW:

(Two-stage-Concurrent): maximize $\displaystyle\sum_{i\in N}s_{i}\sum_{j\in
N^{+}\setminus\\{i\\}}x_{ij}+\min_{\mathbf{d}\in\mathcal{U}}RC(\mathbf{x},\mathbf{d})$
(10a) subject to $\displaystyle\eqref{opc:2}-\eqref{opc:6},\eqref{2ro:2}$
(10b)

where $RC(\mathbf{x},\mathbf{d})$ is the Recourse-Concurrent model and
$\mathcal{U}$ is the box uncertainty set.

We also consider the corresponding static robust model of Two-stage-Concurrent
instead of solving the two-stage robust model directly, the corresponding
static robust model of Two-stage-Concurrent is formulated as follows:

(Static-Concurrent): maximize $\displaystyle\sum_{i\in N}s_{i}\sum_{j\in
N^{+}\setminus\\{i\\}}x_{ij}-\sum_{j\in N}s_{j}\sum_{i\in
N^{+}\setminus\\{j\\}}y_{ij}$ (11a) subject to
$\displaystyle\eqref{opc:2}-\eqref{opc:6},\eqref{2ro:2}$ (11b)
$\displaystyle\eqref{2rm2:2}-\eqref{2rm2:4},\eqref{2rm2:6}$ (11c)
$\displaystyle\sum_{(i,j)\in A}d_{ij}x_{ij}-\sum_{(i,j)\in
A}d_{ij}y_{ij}+\sum_{j\in N}\left(\sum_{k\in
N^{+}\setminus\\{j\\}}y_{jk}-\sum_{i\in
N^{+}\setminus\\{j\\}}y_{ij}\right)\bar{d}_{j0}\leq
L,\forall\mathbf{d}\in\mathcal{U}$ (11d)

We can readily see that the optimal solution of the static robust model
Static-Concurrent is feasible to the two-stage robust model Two-stage-
Concurrent. With the help of the box uncertainty set, the following theorem
can be established which shows that Two-stage-Concurrent and Static-Concurrent
are equivalent.

######  Theorem 2.

The two-stage robust model Two-stage-Concurrent and its corresponding static
robust model Static-Concurrent are equivalent.

###### Proof.

It is clear that the optimal solution of the static robust model Static-
Concurrent is a feasible solution of the two-stage robust model Two-stage-
Concurrent. All we need to show is that the optimal solution of the two-stage
robust model Two-stage-Concurrent is a feasible solution of the static robust
model Static-Concurrent.

We prove by apagoge. Denote $(x_{ij}^{*},y_{ij}^{*})$ as the optimal solution
of the two-stage robust model Two-stage-Concurrent. Suppose the optimal
solution is infeasible for the static robust model Static-Concurrent, which
means

$\exists\mathbf{d}^{\prime}\in\mathcal{U},\sum_{(i,j)\in
A}d^{\prime}_{ij}x^{*}_{ij}-\sum_{(i,j)\in
A}d^{\prime}_{ij}y^{*}_{ij}+\sum_{j\in N}\left(\sum_{k\in
N^{+}\setminus\\{j\\}}y^{*}_{jk}-\sum_{i\in
N^{+}\setminus\\{j\\}}y^{*}_{ij}\right)\bar{d}_{j0}>L$ (12)

Denote $\mathbf{d}^{u}=\bar{\mathbf{d}}+\Theta\hat{\mathbf{d}}$, according to
the definition of the box uncertainty set $\mathcal{U}$, we know that
$\mathbf{d}^{u}\in\mathcal{U}$ and $\mathbf{d}^{u}\succeq\mathbf{d}^{\prime}$
where $\succeq$ is the element-wise inequality. Based on inequality (12), it
is clear that

$\sum_{(i,j)\in A}d^{u}_{ij}x^{*}_{ij}-\sum_{(i,j)\in
A}d^{u}_{ij}y^{*}_{ij}+\sum_{j\in N}\left(\sum_{k\in
N^{+}\setminus\\{j\\}}y^{*}_{jk}-\sum_{i\in
N^{+}\setminus\\{j\\}}y^{*}_{ij}\right)\bar{d}_{j0}>L$ (13)

Denote $\mathbf{d}^{*}\in\mathcal{U}$ as the optimal value of $\mathbf{d}$
that achieves optimal solution $(x_{ij}^{*},y_{ij}^{*})$ in the two-stage
robust model Two-stage-Concurrent, then

$\sum_{(i,j)\in A}d^{*}_{ij}x^{*}_{ij}-\sum_{(i,j)\in
A}d^{*}_{ij}y^{*}_{ij}+\sum_{j\in N}\left(\sum_{k\in
N^{+}\setminus\\{j\\}}y^{*}_{jk}-\sum_{i\in
N^{+}\setminus\\{j\\}}y^{*}_{ij}\right)\bar{d}_{j0}\leq L$ (14)

and

$\sum_{(i,j)\in A}d^{*}_{ij}x^{*}_{ij}-\sum_{(i,j)\in
A}d^{*}_{ij}y^{\prime}_{ij}+\sum_{j\in N}\left(\sum_{k\in
N^{+}\setminus\\{j\\}}y^{\prime}_{jk}-\sum_{i\in
N^{+}\setminus\\{j\\}}y^{\prime}_{ij}\right)\bar{d}_{j0}>L,\forall\mathbf{y}^{\prime}\in\mathcal{Y}\text{
and }\mathbf{y}^{\prime}\neq\mathbf{y}^{*}$ (15)

where
$\mathcal{Y}=\left\\{\mathbf{y}^{\prime}:\begin{matrix}y^{\prime}_{ij}\leq
y^{*}_{ij},\forall i,j\\\ \sum_{i\in
N^{+}\setminus\\{j\\}}y^{\prime}_{ij}\leq\sum_{k\in
N^{+}\setminus\\{j\\}}y^{\prime}_{jk},\forall j\\\
y^{\prime}_{ij}\in\\{0,1\\},\forall i,j\end{matrix}\right\\}$ is the set which
contains all recourse actions with less cancelled arcs comparing with
$\mathbf{y}^{*}$.

Based on the fact that $\mathbf{d}^{u}\succeq\mathbf{d}^{*}$ and using
inequality (15), we have

$\sum_{(i,j)\in A}d^{u}_{ij}x^{*}_{ij}-\sum_{(i,j)\in
A}d^{u}_{ij}y^{\prime}_{ij}+\sum_{j\in N}\left(\sum_{k\in
N^{+}\setminus\\{j\\}}y^{\prime}_{jk}-\sum_{i\in
N^{+}\setminus\\{j\\}}y^{\prime}_{ij}\right)\bar{d}_{j0}>L,\forall\mathbf{y}^{\prime}\in\mathcal{Y}\text{
and }\mathbf{y}^{\prime}\neq\mathbf{y}^{*}$ (16)

Combining inequalities (13) and (16), we can observe that: for the second-
stage problem $RC(\mathbf{x}^{*},\mathbf{d}^{u})$, the recourse action needs
to cancel more arcs than $\mathbf{y}^{*}$ to satisfy the length constraint,
which means $RC(\mathbf{x}^{*},\mathbf{d}^{u})<-\sum_{j\in N}s_{j}\sum_{i\in
N^{+}\setminus\\{j\\}}y^{*}_{ij}$. Because we are optimizing
$\min_{\mathbf{d}\in\mathcal{U}}RC(\mathbf{x}^{*},\mathbf{d})$,
$\mathbf{y}^{*}$ is not the optimal second-stage solution, so this is a
contradiction. Hence, the hypothesis cannot be established, which means
$(x_{ij}^{*},y_{ij}^{*})$ is feasible for the static robust model Static-
Concurrent.

We conclude that the optimal solution of the two-stage robust model Two-stage-
Concurrent is a feasible solution of the static robust model Static-
Concurrent. Because the two models have the same objective function value with
the same solutions, then the optimal solution of the two-stage robust model
Two-stage-Concurrent is also the optimal solution of the static robust model
Static-Concurrent, this implies the two models are equivalent. ∎

Based on Theorem 2, the two-stage robust model Two-stage-Concurrent can be
solved to optimality by solving its corresponding static robust model Static-
Concurrent.

###  4.3 The relationship between Two-stage-Sequential and Two-stage-
Concurrent

Until now, we have introduced two two-stage robust models Two-stage-Sequential
and Two-stage-Concurrent, and also proved that these two models are equivalent
to their corresponding static robust models respectively. In this subsection,
we further investigate the relationships between Two-stage-Sequential and Two-
stage-Concurrent.

First, we investigate the static models Static-Concurrent and Static-
Sequential. Comparing static model Static-Concurrent with static model Static-
Sequential, model Static-Concurrent has less decision variables and less
constraints and is computationally more attractive. The two static models are
based on different recourse models. Next, we show that model Static-Sequential
and model Static-Concurrent are equivalent with the support of the box
uncertainty set.

######  Theorem 3.

The static robust models Static-Concurrent and Static-Sequential are
equivalent.

###### Proof.

First, suppose $(x_{ij}^{*},x_{ijk}^{*},y_{i}^{*},z_{k}^{*})$ is the optimal
solution of the static robust model Static-Sequential, and we introduce
decision variable $y_{ij}$ which is described in Recourse-Concurrent model for
model Static-Sequential, then the optimal solution
$(x_{ij}^{*},x_{ijk}^{*},y_{i}^{*},z_{k}^{*})$ can be mapped to an optimal
solution $(x_{ij}^{*},y_{ij}^{*})$ of model Static-Sequential, and it is clear
that $(x_{ij}^{*},y_{ij}^{*})$ is feasible to the static robust model Static-
Concurrent.

Then, suppose $(x_{ij}^{*},y_{ij}^{*})$ is the optimal solution of the static
robust model Static-Concurrent, and we introduce variables
$x_{ijk},y_{i},z_{k}$ which is described in Recourse-Sequential model for
model Static-Concurrent, then the optimal solution $(x_{ij}^{*},y_{ij}^{*})$
can be mapped to an optimal solution
$(x_{ij}^{*},x_{ijk}^{*},y_{i}^{*},z_{k}^{*})$ of model Static-Concurrent, we
now show that $(x_{ij}^{*},x_{ijk}^{*},y_{i}^{*},z_{k}^{*})$ is feasible to
the static robust model Static-Sequential.

Suppose $z^{*}_{K^{\prime}}=0$ and $z^{*}_{K^{\prime}+1}=1$, then this means
the nodes in the first-stage path become unreachable from the
$(K^{\prime}+1)$th node. Then, the length constraint (11d) in the static
robust model Static-Concurrent is equivalent to

$\sum_{k=1}^{K^{\prime}}\sum_{(i,j)\in A}d_{ij}x^{*}_{ijk}+\sum_{(i,j)\in
A}\bar{d}_{j0}x^{*}_{ijK^{\prime}}\leq L,\forall\mathbf{d}\in\mathcal{U}$ (17)

Based on the definition of the box uncertainty set $\mathcal{U}$, inequality
(17) is equivalent to

$\sum_{k=1}^{K^{\prime}}\sum_{(i,j)\in
A}(\bar{d}_{ij}+\Theta\hat{d}_{ij})x^{*}_{ijk}+\sum_{(i,j)\in
A}\bar{d}_{j0}x^{*}_{ijK^{\prime}}\leq L$ (18)

We transform the left hand side of the above inequality as follows:

$\displaystyle\sum_{k=1}^{K^{\prime}}\sum_{(i,j)\in
A}(\bar{d}_{ij}+\Theta\hat{d}_{ij})x^{*}_{ijk}+\sum_{(i,j)\in
A}\bar{d}_{j0}x^{*}_{ijK^{\prime}}$ (19)
$\displaystyle=\sum_{k=1}^{K^{\prime}-1}\sum_{(i,j)\in
A}(\bar{d}_{ij}+\Theta\hat{d}_{ij})x^{*}_{ijk}+\sum_{(i,j)\in
A}(\bar{d}_{ij}+\Theta\hat{d}_{ij})x^{*}_{ijK^{\prime}}+\sum_{(i,j)\in
A}\bar{d}_{j0}x^{*}_{ijK^{\prime}}$
$\displaystyle\geq\sum_{k=1}^{K^{\prime}-1}\sum_{(i,j)\in
A}(\bar{d}_{ij}+\Theta\hat{d}_{ij})x^{*}_{ijk}+\sum_{(i,j)\in
A}(\bar{d}_{ij}+\bar{d}_{j0})x^{*}_{ijK^{\prime}}$
$\displaystyle\overset{1}{>}\sum_{k=1}^{K^{\prime}-1}\sum_{(i,j)\in
A}(\bar{d}_{ij}+\Theta\hat{d}_{ij})x^{*}_{ijk}+\sum_{(i,j)\in
A}\bar{d}_{j0}x^{*}_{ij(K^{\prime}-1)}$

where relation 1 is due to the triangle inequality.

Based on (18) and (19) we have

$\sum_{k=1}^{K^{\prime}-1}\sum_{(i,j)\in
A}(\bar{d}_{ij}+\Theta\hat{d}_{ij})x^{*}_{ijk}+\sum_{(i,j)\in
A}\bar{d}_{j0}x^{*}_{ij(K^{\prime}-1)}\leq L$ (20)

which implies

$\sum_{k=1}^{K^{\prime}-1}\sum_{(i,j)\in A}d_{ij}x^{*}_{ijk}+\sum_{(i,j)\in
A}\bar{d}_{j0}x^{*}_{ij(K^{\prime}-1)}\leq L,\forall\mathbf{d}\in\mathcal{U}$
(21)

Following the above transformation recursively, finally we can get

$\sum_{k=1}^{K}\sum_{(i,j)\in A}d_{ij}x^{*}_{ijk}+\sum_{(i,j)\in
A}\bar{d}_{j0}x^{*}_{ijK}\leq L,\forall
K=1,...,K^{\prime},\mathbf{d}\in\mathcal{U}$ (22)

So the optimal solution $(x_{ij}^{*},x_{ijk}^{*},y_{i}^{*},z_{k}^{*})$ is
feasible to the static robust model Static-Sequential.

We conclude that the optimal solution of the static robust model Static-
Sequential is a feasible solution of the static robust model Static-
Concurrent, and the optimal solution of the static robust model Static-
Concurrent is a feasible solution of the static robust model Static-
Sequential. Then the two models have the same optimal solution, this implies
the two models are equivalent. ∎

The following corollary shows the equivalence between two two-stage robust
models.

######  Corollary 1.

The two-stage robust models Two-stage-Sequential and Two-stage-Concurrent are
equivalent.

###### Proof.

Based on Theorem 1, Theorem 2 and Theorem 3, we can draw this conclusion. ∎

Based on Theorems 1-3 and Corollary 1, we know that the four models Two-stage-
Sequential, Static-Sequential, Two-stage-Concurrent and Static-Concurrent are
equivalent to each other. It is an interesting conclusion that Two-stage-
Sequential and Two-stage-Concurrent are equivalent with the _box uncertainty
set_ defined even though they are based on different recourse models. We can
use the Two-stage-Concurrent model to deal with the two-stage OPSW with
sequentially realized weights, and the Two-stage-Concurrent model is
computationally more efficient than the Two-stage-Sequential model.

##  5 Case study

In this section, a case study is presented to illustrate the effectiveness of
the proposed two-stage robust models for OPSW.

###  5.1 Test instance

The test instance used in our experiments is based on problem set 3 from
Tsiligirides (1984) which was originally used for the deterministic OP.
Problem set 3 contains 20 instances with the same 33 nodes and 20 varying
length limits. We only consider 3 length limits: 80, 90 and 100. In the
instance, the end point is ignored and the start point is kept as the depot
location. The problem set can be found with URL:
http://www.mech.kuleuven.be/en/cib/op.

To generate the uncertain instances for OPSW, we use the Euclidean distances
between nodes as the expected weights $\bar{d}_{ij}$. Two kinds of uncertain
instances are generated with the deviation values $\hat{d}_{ij}$ chosen as
$0.2\bar{d}_{ij}$ and $0.5\bar{d}_{ij}$ respectively. Then, based on different
length limits and different deviation values, we can get a total of 6
uncertain instances for OPSW.

###  5.2 Experiments

In order to evaluate the effectiveness of the proposed two-stage robust
models, we use one-stage robust model for OPSW as a comparison. In the one-
stage robust model, all the decision variables are “here and now” and there
are no recourse decision variables considered in the model. The one-stage
robust model follows the traditional RO paradigm and is formulated as follows:

(One-stage-RO): maximize $\displaystyle\sum_{i\in N}s_{i}\sum_{j\in
N^{+}\setminus\\{i\\}}x_{ij}$ (23a) subject to
$\displaystyle\eqref{opc:2}-\eqref{opc:6}$ (23b) $\displaystyle\sum_{(i,j)\in
A}d_{ij}x_{ij}\leq L,\forall\mathbf{d}\in\mathcal{U}$ (23c)

As proven in Section 4, the Static-Concurrent model is equivalent to the two-
stage robust models Two-stage-Sequential and Two-stage-Concurrent and is
computationally more efficient than the Static-Sequential model. So we use
Static-Concurrent model to solve the two-stage OPSW with sequential or
concurrent realized weights.

We generate 1000 scenarios for the $\hat{d}_{ij}=0.2\bar{d}_{ij}$ and
$\hat{d}_{ij}=0.5\bar{d}_{ij}$ cases respectively for simulation purposes.
Suppose the uncertain weights $\tilde{d}_{ij}$ are uniformly distributed on
interval $[\bar{d}_{ij}-\hat{d}_{ij},\bar{d}_{ij}+\hat{d}_{ij}]$. Then the
realizations $d_{ij}$ are sampled uniformly on the interval
$[\bar{d}_{ij}-\hat{d}_{ij},\bar{d}_{ij}+\hat{d}_{ij}]$.

For each uncertain instance, the Static-Concurrent model and the One-stage-RO
model are solved by CPLEX 12.6 with $\Theta=0,0.1,...,1$ respectively. The
robust solutions obtained by the Static-Concurrent model and the One-stage-RO
model are then simulated with the 1000 scenarios for the two-stage OPSW with
Recourse-Sequential and Recourse-Concurrent actions. The objective values of
the second-stage recourse problems are calculated by Algorithm 1 and Algorithm
2. The mean objective values and the standard deviations of the robust
solutions are statistically summarized.

###  5.3 Numerical results

Tables 1-6 show the numerical results of the 6 instances with different length
limits and different deviation values. The Obj. in the tables represents the
objective value obtained by the One-stage-RO model or the Static-Concurrent
model. First we can observe that the objective values of the robust solutions
are decreasing as parameter $\Theta$ increases for both the one-stage and two-
stage RO models. As $\Theta$ increases, the size of the uncertainty set
$\mathcal{U}$ is increasing which means the protection level is increasing,
and the resulting robust solution is more conservative.

Table 1: Numerical results of the instance with $L=80$ and $\hat{d}_{ij}=0.2\bar{d}_{ij}$ $\Theta$ | One-stage RO | Two-stage RO
---|---|---
Obj. | Sequential | Concurrent | Obj. | Sequential | Concurrent
Mean | Std. | Mean | Std. | Mean | Std. | Mean | Std.
0.00 | 710.00 | 680.19 | 41.66 | 681.63 | 41.43 | 710.00 | 693.76 | 27.37 | 695.29 | 24.89
0.10 | 690.00 | 676.52 | 28.56 | 677.93 | 27.62 | 700.00 | 691.55 | 19.85 | 692.43 | 18.59
0.20 | 690.00 | 688.48 | 7.08 | 688.82 | 6.36 | 690.00 | 685.64 | 18.87 | 685.70 | 18.79
0.30 | 680.00 | 679.58 | 6.02 | 679.58 | 6.02 | 680.00 | 679.93 | 1.38 | 679.94 | 1.34
0.40 | 670.00 | 670.00 | 0.00 | 670.00 | 0.00 | 670.00 | 679.77 | 1.50 | 679.78 | 1.47
0.50 | 660.00 | 660.00 | 0.00 | 660.00 | 0.00 | 660.00 | 660.00 | 0.00 | 660.00 | 0.00
0.60 | 650.00 | 650.00 | 0.00 | 650.00 | 0.00 | 650.00 | 650.00 | 0.00 | 650.00 | 0.00
0.70 | 640.00 | 640.00 | 0.00 | 640.00 | 0.00 | 640.00 | 640.68 | 5.17 | 640.68 | 5.17
0.80 | 630.00 | 630.00 | 0.00 | 630.00 | 0.00 | 640.00 | 640.00 | 0.00 | 640.00 | 0.00
0.90 | 630.00 | 630.00 | 0.00 | 630.00 | 0.00 | 630.00 | 630.00 | 0.00 | 630.00 | 0.00
1.00 | 620.00 | 620.00 | 0.00 | 620.00 | 0.00 | 630.00 | 630.00 | 0.00 | 630.00 | 0.00
Table 2: Numerical results of the instance with $L=80$ and $\hat{d}_{ij}=0.5\bar{d}_{ij}$ $\Theta$ | One-stage RO | Two-stage RO
---|---|---
Obj. | Sequential | Concurrent | Obj. | Sequential | Concurrent
Mean | Std. | Mean | Std. | Mean | Std. | Mean | Std.
0.00 | 710.00 | 652.47 | 77.39 | 657.90 | 73.04 | 710.00 | 652.47 | 77.39 | 657.90 | 73.04
0.10 | 680.00 | 662.97 | 40.82 | 664.29 | 39.68 | 680.00 | 680.41 | 25.69 | 681.26 | 24.96
0.20 | 660.00 | 658.71 | 8.26 | 658.85 | 7.74 | 660.00 | 655.91 | 17.83 | 656.05 | 17.58
0.30 | 640.00 | 639.83 | 2.07 | 639.86 | 2.00 | 640.00 | 639.48 | 6.43 | 639.54 | 6.15
0.40 | 620.00 | 620.00 | 0.00 | 620.00 | 0.00 | 630.00 | 629.99 | 0.32 | 629.99 | 0.32
0.50 | 610.00 | 610.00 | 0.00 | 610.00 | 0.00 | 610.00 | 631.71 | 19.40 | 633.11 | 19.23
0.60 | 590.00 | 590.00 | 0.00 | 590.00 | 0.00 | 600.00 | 600.00 | 0.00 | 600.00 | 0.00
0.70 | 570.00 | 570.00 | 0.00 | 570.00 | 0.00 | 580.00 | 580.00 | 0.00 | 580.00 | 0.00
0.80 | 570.00 | 570.00 | 0.00 | 570.00 | 0.00 | 570.00 | 600.00 | 0.00 | 600.00 | 0.00
0.90 | 550.00 | 550.00 | 0.00 | 550.00 | 0.00 | 560.00 | 593.93 | 24.90 | 596.25 | 28.57
1.00 | 540.00 | 540.00 | 0.00 | 540.00 | 0.00 | 550.00 | 570.00 | 0.00 | 570.00 | 0.00
Table 3: Numerical results of the instance with $L=90$ and $\hat{d}_{ij}=0.2\bar{d}_{ij}$ $\Theta$ | One-stage RO | Two-stage RO
---|---|---
Obj. | Sequential | Concurrent | Obj. | Sequential | Concurrent
Mean | Std. | Mean | Std. | Mean | Std. | Mean | Std.
0.00 | 770.00 | 744.96 | 35.74 | 748.44 | 33.24 | 770 | 728.16 | 44.97 | 729.68 | 45.17
0.10 | 760.00 | 750.74 | 25.94 | 751.33 | 25.31 | 760 | 754.09 | 17.22 | 754.47 | 16.44
0.20 | 750.00 | 746.09 | 17.02 | 746.51 | 16.21 | 750 | 752.69 | 10.43 | 752.93 | 9.90
0.30 | 740.00 | 739.80 | 2.36 | 739.85 | 2.25 | 740 | 747.12 | 5.07 | 747.26 | 4.81
0.40 | 730.00 | 730.00 | 0.00 | 730.00 | 0.00 | 730 | 730.00 | 0.00 | 730.00 | 0.00
0.50 | 720.00 | 720.00 | 0.00 | 720.00 | 0.00 | 720 | 720.00 | 0.00 | 720.00 | 0.00
0.60 | 710.00 | 710.00 | 0.00 | 710.00 | 0.00 | 710 | 710.00 | 0.00 | 710.00 | 0.00
0.70 | 700.00 | 700.00 | 0.00 | 700.00 | 0.00 | 710 | 736.43 | 4.83 | 736.71 | 4.74
0.80 | 690.00 | 690.00 | 0.00 | 690.00 | 0.00 | 690 | 690.00 | 0.00 | 690.00 | 0.00
0.90 | 680.00 | 680.00 | 0.00 | 680.00 | 0.00 | 690 | 700.00 | 0.00 | 700.00 | 0.00
1.00 | 670.00 | 670.00 | 0.00 | 670.00 | 0.00 | 680 | 680.00 | 0.00 | 680.00 | 0.00
Table 4: Numerical results of the instance with $L=90$ and $\hat{d}_{ij}=0.5\bar{d}_{ij}$ $\Theta$ | One-stage RO | Two-stage RO
---|---|---
Obj. | Sequential | Concurrent | Obj. | Sequential | Concurrent
Mean | Std. | Mean | Std. | Mean | Std. | Mean | Std.
0.00 | 770.00 | 722.34 | 58.05 | 727.46 | 57.50 | 770 | 719.97 | 59.25 | 725.56 | 58.42
0.10 | 740.00 | 718.99 | 48.56 | 721.56 | 43.38 | 740 | 730.86 | 38.69 | 733.30 | 36.30
0.20 | 720.00 | 718.37 | 10.40 | 718.43 | 10.37 | 720 | 716.06 | 16.54 | 716.74 | 15.18
0.30 | 690.00 | 689.98 | 0.45 | 689.98 | 0.45 | 700 | 699.83 | 2.88 | 699.88 | 2.57
0.40 | 670.00 | 670.00 | 0.00 | 670.00 | 0.00 | 680 | 683.83 | 12.27 | 684.13 | 13.68
0.50 | 650.00 | 650.00 | 0.00 | 650.00 | 0.00 | 660 | 679.90 | 1.18 | 679.92 | 1.09
0.60 | 640.00 | 640.00 | 0.00 | 640.00 | 0.00 | 640 | 701.27 | 26.59 | 702.86 | 26.30
0.70 | 620.00 | 620.00 | 0.00 | 620.00 | 0.00 | 630 | 630.00 | 0.00 | 630.00 | 0.00
0.80 | 610.00 | 610.00 | 0.00 | 610.00 | 0.00 | 610 | 650.00 | 13.75 | 650.21 | 14.85
0.90 | 600.00 | 600.00 | 0.00 | 600.00 | 0.00 | 600 | 649.83 | 26.63 | 652.74 | 31.74
1.00 | 580.00 | 580.00 | 0.00 | 580.00 | 0.00 | 590 | 604.88 | 13.09 | 608.52 | 16.38
Table 5: Numerical results of the instance with $L=100$ and $\hat{d}_{ij}=0.2\bar{d}_{ij}$ $\Theta$ | One-stage RO | Two-stage RO
---|---|---
Obj. | Sequential | Concurrent | Obj. | Sequential | Concurrent
Mean | Std. | Mean | Std. | Mean | Std. | Mean | Std.
0.00 | 800 | 795.39 | 7.86 | 795.68 | 7.84 | 800 | 776.88 | 45.66 | 778.99 | 41.13
0.10 | 800 | 798.06 | 5.57 | 798.14 | 5.53 | 800 | 797.02 | 7.11 | 797.20 | 7.05
0.20 | 790 | 789.83 | 1.86 | 789.83 | 1.86 | 790 | 789.77 | 1.86 | 789.78 | 1.83
0.30 | 790 | 789.98 | 0.45 | 789.99 | 0.32 | 790 | 789.89 | 1.70 | 789.89 | 1.70
0.40 | 780 | 780.00 | 0.00 | 780.00 | 0.00 | 790 | 790.00 | 0.00 | 790.00 | 0.00
0.50 | 780 | 780.00 | 0.00 | 780.00 | 0.00 | 780 | 780.00 | 0.00 | 780.00 | 0.00
0.60 | 760 | 760.00 | 0.00 | 760.00 | 0.00 | 770 | 770.00 | 0.00 | 770.00 | 0.00
0.70 | 760 | 760.00 | 0.00 | 760.00 | 0.00 | 760 | 760.00 | 0.00 | 760.00 | 0.00
0.80 | 740 | 740.00 | 0.00 | 740.00 | 0.00 | 750 | 769.39 | 2.39 | 769.48 | 2.22
0.90 | 740 | 740.00 | 0.00 | 740.00 | 0.00 | 740 | 740.00 | 0.00 | 740.00 | 0.00
1.00 | 730 | 730.00 | 0.00 | 730.00 | 0.00 | 730 | 730.00 | 0.00 | 730.00 | 0.00
Table 6: Numerical results of the instance with $L=100$ and $\hat{d}_{ij}=0.5\bar{d}_{ij}$ $\Theta$ | One-stage RO | Two-stage RO
---|---|---
Obj. | Sequential | Concurrent | Obj. | Sequential | Concurrent
Mean | Std. | Mean | Std. | Mean | Std. | Mean | Std.
0.00 | 800 | 778.79 | 43.57 | 779.34 | 43.51 | 800 | 787.07 | 19.79 | 787.63 | 18.57
0.10 | 790 | 774.72 | 40.31 | 776.13 | 38.69 | 790 | 776.04 | 34.31 | 778.05 | 32.25
0.20 | 780 | 777.50 | 14.56 | 778.10 | 12.25 | 780 | 778.71 | 8.64 | 778.81 | 8.18
0.30 | 750 | 749.71 | 4.56 | 749.77 | 4.15 | 760 | 759.72 | 3.15 | 759.79 | 2.84
0.40 | 730 | 730.00 | 0.00 | 730.00 | 0.00 | 730 | 749.27 | 3.57 | 749.33 | 3.44
0.50 | 710 | 710.00 | 0.00 | 710.00 | 0.00 | 710 | 710.00 | 0.00 | 710.00 | 0.00
0.60 | 690 | 690.00 | 0.00 | 690.00 | 0.00 | 690 | 719.75 | 2.29 | 719.81 | 2.02
0.70 | 670 | 670.00 | 0.00 | 670.00 | 0.00 | 670 | 718.05 | 38.05 | 720.08 | 38.58
0.80 | 650 | 650.00 | 0.00 | 650.00 | 0.00 | 660 | 720.00 | 0.00 | 720.00 | 0.00
0.90 | 630 | 630.00 | 0.00 | 630.00 | 0.00 | 640 | 650.00 | 0.00 | 650.00 | 0.00
1.00 | 620 | 620.00 | 0.00 | 620.00 | 0.00 | 630 | 640.00 | 0.00 | 640.00 | 0.00

For both the one-stage and two-stage robust models, the mean objective values
with concurrent recourse are greater than or equal to the corresponding mean
objective values with sequential recourse. The reason is that the concurrent
recourse has more information on the uncertainty realizations than the
sequential recourse, so the concurrent recourse can make a better recourse
decision and achieve a lower loss of the collected score. However, the gaps
between the mean objective values with sequential recourse and concurrent
recourse are very small which means that the difference between the two
recourse actions is small.

Comparing the mean objective values of the one-stage RO and two-stage RO with
the sequential recourse or concurrent recourse, the two-stage RO achieves
better values than the one-stage RO in most cases. This is because that the
two-stage RO considers the recourse decisions into the model but the one-stage
RO only considers the first-stage decisions. We also notice that the one-stage
RO can achieve better mean objective values than the two-stage RO in some
cases. Figures 1-3 show visual comparisions between the one-stage RO and two-
stage RO with sequential recourse. The comparisions between the one-stage RO
and two-stage RO with concurrent recourse are visually similar. The figures
clearly show that the two-stage RO dominates the one-stage RO in most cases,
which show the effectiveness and superiority of the proposed two-stage robust
models for dealing with the two-stage OPSW.

In Tables 1-6 we also report the standard deviations of the simulated robust
solutions for both the one-stage RO and the two-stage RO. The standard
deviations can reflect the stabilities of the obtained robust solutions. From
the tables we can see that as parameter $\Theta$ increasing, the standard
deviations tend to decrease, which means the robust solutions are more stable
with a larger uncertainty set. We can also observe that the two-stage RO can
mostly achieve better mean objective values and lower or small standard
deviation values at the same time comparing with the one-stage RO. This
further indicates that the proposed two-stage robust models can efficiently
tackle the two-stage OPSW.

Figure 1: Comparision between the one-stage RO and two-stage RO with
sequential recourse, $L=80$, $\hat{d}_{ij}=0.2\bar{d}_{ij}$(left),
$\hat{d}_{ij}=0.5\bar{d}_{ij}$(right) Figure 2: Comparision between the one-
stage RO and two-stage RO with sequential recourse, $L=90$,
$\hat{d}_{ij}=0.2\bar{d}_{ij}$(left), $\hat{d}_{ij}=0.5\bar{d}_{ij}$(right)
Figure 3: Comparision between the one-stage RO and two-stage RO with
sequential recourse, $L=100$, $\hat{d}_{ij}=0.2\bar{d}_{ij}$(left),
$\hat{d}_{ij}=0.5\bar{d}_{ij}$(right)

##  6 Conclusions

In this paper, we considered the OPSW with recourse actions. Based on
different uncertainty realization ways, we presented two recourse models: one
is the Recourse-Sequential model and the other is the Recourse-Concurrent
model. The Recourse-Concurrent model has less decision variables and less
constraints and is computationally more attractive. We applied the two-stage
RO paradigm to the OPSW and introduced two two-stage RO models based on two
recourse models. We theoretically proved that with the _box uncertainty set_
defined, the two-stage robust models are equivalent to their corresponding
static robust models and the two two-stage robust models are also equivalent
to each other. Subsequently, the two-stage robust models for OPSW can be
determined to optimality by solving their corresponding static models.
Comparative studies between the two-stage robust models and one-stage robust
model for OPSW showed the effectiveness and superiority of the proposed two-
stage robust models for tackling the two-stage OPSW.

We provide the following research directions as our future works:

  1. 1.

The two-stage robust models for OPSW proposed in this paper are based on the
_box uncertainty set_ , therefore we can draw theoretical conclusions on the
equivalence between the two-stage robust models and their corresponding static
robust models. Other uncertainty sets (e.g. the polyhedral uncertianty set)
could be defined in the two-stage robust models and the performance of the
corresponding static robust models can be studied.

  2. 2.

The OPSW considered in this paper is with a two-stage setting where the
decision variables are classified into two categories. As the planned path is
executed dynamically and the nodes are visited sequentially, hence the OPSW
can be viewed as a multi-stage decision making problem. So we can apply the
multi-stage RO methodology and build a multi-stage robust model for the OPSW
with a multi-stage setting.

## Conflict of Interests

The authors declare that there is no conflict of interest regarding the
publication of this manuscript.

## Acknowledgments

This work was supported by National Natural Science Foundation of China (No.
61573277, 71471158), the Research Grants Council of the Hong Kong Special
Administrative Region, China (Project No. PolyU 15201414), the Fundamental
Research Funds for the Central Universities, the Open Research Fund of the
State Key Laboratory of Astronautic Dynamics under Grant 2015ADL-DW403, and
the Scientific Research Foundation for the Returned Overseas Chinese Scholars,
State Education Ministry, Natural Science Basic Research Plan in Shaanxi
Province of China (No. 2015JM6316). The authors also would like to thank The
Hong Kong Polytechnic University Research Committee for financial and
technical support.

## References

  * Evers et al. (2014) L. Evers, K. Glorie, S. Van Der Ster, A. I. Barros, H. Monsuur,  A two-stage approach to the orienteering problem with stochastic weights,  Computers & Operations Research 43 (2014) 248–260.
  * Golden et al. (1987) B. L. Golden, L. Levy, R. Vohra,  The orienteering problem,  Naval research logistics 34 (1987) 307–318.
  * Mufalli et al. (2012) F. Mufalli, R. Batta, R. Nagi,  Simultaneous sensor selection and routing of unmanned aerial vehicles for complex mission plans,  Computers & Operations Research 39 (2012) 2787–2799.
  * Evers et al. (2014) L. Evers, T. Dollevoet, A. I. Barros, H. Monsuur,  Robust uav mission planning,  Annals of Operations Research 222 (2014) 293–315.
  * Vansteenwegen and Van Oudheusden (2007) P. Vansteenwegen, D. Van Oudheusden,  The mobile tourist guide: an or opportunity,  OR insight 20 (2007) 21–27.
  * Gavalas et al. (2014) D. Gavalas, C. Konstantopoulos, K. Mastakas, G. Pantziou,  A survey on algorithmic approaches for solving tourist trip design problems,  Journal of Heuristics 20 (2014) 291–328.
  * Howe (2008) J. Howe, Crowdsourcing: How the power of the crowd is driving the future of business, Random House, 2008\.
  * Yuen et al. (2011) M.-C. Yuen, I. King, K.-S. Leung,  A survey of crowdsourcing systems,  in: Privacy, Security, Risk and Trust (PASSAT) and 2011 IEEE Third Inernational Conference on Social Computing (SocialCom), 2011 IEEE Third International Conference on, IEEE, pp. 766–773.
  * Vansteenwegen et al. (2011) P. Vansteenwegen, W. Souffriau, D. Van Oudheusden,  The orienteering problem: A survey,  European Journal of Operational Research 209 (2011) 1–10.
  * Gunawan et al. (2016) A. Gunawan, H. C. Lau, P. Vansteenwegen,  Orienteering problem: A survey of recent variants, solution approaches and applications,  European Journal of Operational Research (2016).
  * Ilhan et al. (2008) T. Ilhan, S. M. Iravani, M. S. Daskin,  The orienteering problem with stochastic profits,  Iie Transactions 40 (2008) 406–421.
  * Campbell et al. (2011) A. M. Campbell, M. Gendreau, B. W. Thomas,  The orienteering problem with stochastic travel and service times,  Annals of Operations Research 186 (2011) 61–81.
  * Lau et al. (2012) H. C. Lau, W. Yeoh, P. Varakantham, D. T. Nguyen, H. Chen,  Dynamic stochastic orienteering problems for risk-aware applications,  arXiv preprint arXiv:1210.4874 (2012).
  * Varakantham and Kumar (2013) P. Varakantham, A. Kumar,  Optimization approaches for solving chance constrained stochastic orienteering problems,  in: International Conference on Algorithmic DecisionTheory, Springer, pp. 387–398.
  * Zhang et al. (2014) S. Zhang, J. W. Ohlmann, B. W. Thomas,  A priori orienteering with time windows and stochastic wait times at customers,  European Journal of Operational Research 239 (2014) 70–79.
  * Ben-Tal et al. (2004) A. Ben-Tal, A. Goryashko, E. Guslitzer, A. Nemirovski,  Adjustable robust solutions of uncertain linear programs,  Mathematical Programming 99 (2004) 351–376.
  * An and Zeng (2015) Y. An, B. Zeng,  Exploring the modeling capacity of two-stage robust optimization: Variants of robust unit commitment model,  IEEE Transactions on Power Systems 30 (2015) 109–122.
  * Wang et al. (2016) B. Wang, S. Wang, X.-z. Zhou, J. Watada,  Two-stage multi-objective unit commitment optimization under hybrid uncertainties,  IEEE Transactions on Power Systems 31 (2016) 2266–2277.
  * Atamtürk and Zhang (2007) A. Atamtürk, M. Zhang,  Two-stage robust network flow and design under demand uncertainty,  Operations Research 55 (2007) 662–673.
  * Ordóñez and Zhao (2007) F. Ordóñez, J. Zhao,  Robust capacity expansion of network flows,  Networks 50 (2007) 136–145.
  * Takeda et al. (2008) A. Takeda, S. Taguchi, R. H. Tütüncü,  Adjustable robust optimization models for a nonlinear two-period system,  Journal of Optimization Theory and Applications 136 (2008) 275–295.
  * Hanasusanto et al. (2015) G. A. Hanasusanto, D. Kuhn, W. Wiesemann,  K-adaptability in two-stage robust binary programming,  Operations Research 63 (2015) 877–891.
  * Feige et al. (2007) U. Feige, K. Jain, M. Mahdian, V. Mirrokni,  Robust combinatorial optimization with exponential scenarios,  in: International Conference on Integer Programming and Combinatorial Optimization, Springer, pp. 439–453.
  * Bertsimas et al. (2015) D. Bertsimas, V. Goyal, B. Y. Lu,  A tight characterization of the performance of static solutions in two-stage adjustable robust linear optimization,  Mathematical Programming 150 (2015) 281–319.
  * Bertsimas et al. (2004) D. Bertsimas, D. Pachamanova, M. Sim,  Robust linear optimization under general norms,  Operations Research Letters 32 (2004) 510–516.
  * Ben-Tal et al. (2009) A. Ben-Tal, L. El Ghaoui, A. Nemirovski, Robust optimization, Princeton University Press, 2009.
  * Tsiligirides (1984) T. Tsiligirides,  Heuristic methods applied to orienteering,  Journal of the Operational Research Society 35 (1984) 797–809.


# Generation and characterization of ultrabroadband polarization–frequency
hyperentangled photons

Hsuan-Hao Lu<EMAIL_ADDRESS>Quantum Information Science Section, Computational
Sciences and Engineering Division, Oak Ridge National Laboratory, Oak Ridge,
Tennessee 37831, USA  Muneer Alshowkan  Quantum Information Science Section,
Computational Sciences and Engineering Division, Oak Ridge National
Laboratory, Oak Ridge, Tennessee 37831, USA  Karthik V. Myilswamy  School of
Electrical and Computer Engineering and Purdue Quantum Science and Engineering
Institute, Purdue University, West Lafayette, Indiana 47907, USA  Andrew M.
Weiner  School of Electrical and Computer Engineering and Purdue Quantum
Science and Engineering Institute, Purdue University, West Lafayette, Indiana
47907, USA  Joseph M. Lukens  Quantum Information Science Group, Oak Ridge
National Laboratory, Oak Ridge, Tennessee 37831, USA  Research Technology
Office and Quantum Collaborative, Arizona State University, Tempe, Arizona
85287, USA  Nicholas A. Peters  Quantum Information Science Section,
Computational Sciences and Engineering Division, Oak Ridge National
Laboratory, Oak Ridge, Tennessee 37831, USA

13.3(1.4,15) This manuscript has been co-authored by UT-Battelle, LLC, under
contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The US
government retains and the publisher, by accepting the article for
publication, acknowledges that the US government retains a nonexclusive, paid-
up, irrevocable, worldwide license to publish or reproduce the published form
of this manuscript, or allow others to do so, for US government purposes. DOE
will provide public access to these results of federally sponsored research in
accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-
public-access-plan).

We generate ultrabroadband photon pairs entangled in both polarization and
frequency bins through an all-waveguided Sagnac source covering the entire
optical C- and L-bands (1530–1625 nm). We perform comprehensive
characterization of high-fidelity states in multiple dense wavelength-division
multiplexed channels, achieving full tomography of effective four-qubit
systems. Additionally, leveraging the inherent high dimensionality of
frequency encoding and our electro-optic measurement approach, we demonstrate
the scalability of our system to higher dimensions, reconstructing states in a
36-dimensional Hilbert space consisting of two polarization qubits and two
frequency-bin qutrits. Our findings hold potential significance for quantum
networking, particularly dense coding and entanglement distillation in
wavelength-multiplexed quantum networks.

Photonic hyperentanglement typically describes two-photon states that exhibit
simultaneous entanglement in multiple degrees of freedom (DoFs), e.g., orbital
angular momentum, spatial mode, time-frequency, and polarization Kwiat (1997);
Barreiro  _et al._ (2005, 2008); Williams  _et al._ (2017); Reimer  _et al._
(2018); Imany  _et al._ (2019); Ecker  _et al._ (2021); Chen  _et al._ (2022);
Francesconi  _et al._ (2023). The expansion of the Hilbert space enables
deterministic controlled operations between two DoFs within a single photon,
showcasing significant potential for quantum communication protocols including
dense coding Barreiro  _et al._ (2008); Williams  _et al._ (2017) and single-
copy entanglement distillation Simon and Pan (2002); Ecker  _et al._ (2021).
Among the various exploitable DoFs, the polarization DoF has historically
received extensive investigation, primarily due to readily available tools for
state manipulation. On the other hand, time-frequency encoding stands out as a
promising candidate due to its compatibility with established fiber-optic
networks. Specifically, discrete frequency bins Kues  _et al._ (2019); Lu  _et
al._ (2019), a special case under the wider time-frequency umbrella, offer
practical advantages such as straightforward multiplexing, parallel processing
of multiple qubits, and the absence of nested interferometers that typically
require active stabilization.

In this work, we present an all-waveguided, ultrabroadband hyperentangled
source spanning polarization and frequency-bin DoFs. We perform the first full
quantum state tomography (QST) of polarization and frequency-bin
hyperentangled states, covering multiple dense wavelength division
multiplexing channels across the optical C-band (1530–1565 nm) and L-band
(1565–1625 nm) and extend measurements to high-dimensional Hilbert spaces
through frequency-bin _qudit_ encodings. Through the serial application of
polarization projections and electro-optic-based frequency mixing, our scheme
can probe arbitrary bases in the complete two-photon polarization and
frequency Hilbert space, thus facilitating full state reconstruction with no
constraints on the ground truth state. Our procedure is experimentally
demonstrated for up to 36-dimensional hyperentangled systems, yet is in
principle scalable to much higher dimensions. Overall, our source design and
characterization techniques open new avenues for hyperentanglement generation
and manipulation in these two important DoFs.

Figure 1: (a) Experimental setup. (b) Conceptual diagram of the hyperentangled
source. The background traces the SPDC spectrum measured on an optical
spectrum analyzer (shown in linear scale) Alshowkan  _et al._ (2022). The
stars represent the five channel pairs characterized in this study. [CW:
Continuous-wave laser. WP: liquid crystal waveplate. WDM: wavelength division
multiplexer. PBS: polarizing beamsplitter. PPLN: periodically poled lithium
niobate waveguide. EOM: electro-optic phase modulator. WSS: wavelength-
selective switch. SNSPD: superconducting nanowire detector.]

Figure 1(a) depicts the experimental setup, encompassing the hyperentangled
photon source and two stages of state analyzers for each photon: one for
polarization and one for frequency bins. We operate a continuous-wave laser
around 783 nm and pump a 12 mm-long type-0 periodically poled lithium niobate
(PPLN) ridge waveguide (AdvR) in a fiber Sagnac loop Li  _et al._ (2005); Fan
_et al._ (2007); Vergyris  _et al._ (2017); Alshowkan  _et al._ (2022). We
employ a combination of a liquid crystal waveplate (WP; Thorlabs) and a fiber
polarizing beamsplitter (PBS) to split the laser and coherently pump the
waveguide from both directions. Using a 90-degree rotated fiber in one of the
PBS outputs aligns pump photons and generated photon pairs in both directions
to vertical polarization while traversing the PPLN waveguide. Upon
recombination at the PBS and the 780/1560 nm wavelength-division multiplexer,
the generated photon pairs are in the form of a polarization-entangled state
$\ket{\Psi_{P}}\propto\alpha\ket{HH}+\beta\ket{VV}$, supporting a two-photon
bandwidth of $\sim$18 THz, with the signal and idler photon covering more than
the optical C- and L-band, respectively. The present source, which builds upon
the design presented in Ref. Alshowkan  _et al._ (2022), has previously shown
high-fidelity polarization entanglement across 150 pairs of 25 GHz-wide
channels. In this work, we introduce a notable stability improvement by
tapping 1% of the pump power in the Sagnac loop for active control of the WP.
This allows us to maximize the polarization entanglement (i.e., attain
$|\alpha|\approx|\beta|$) and counter real-time polarization drifts in the
pump laser. With this addition, 31 independent polarization QST trials on the
same 25 GHz-wide channel pair were found to yield a mean fidelity of 98.6%
with a standard deviation of only 0.3% over the course of 18 hours.

The broadband spectral coherence and energy correlation between signal and
idler photons also provide a natural resource for investigating frequency-bin
entanglement and higher-dimensional Hilbert spaces Lu  _et al._ (2018, 2019);
Kues  _et al._ (2019); Lu  _et al._ (2022); Chen  _et al._ (2022). Several
techniques have been explored to create discrete frequency bins: one approach
integrates resonant structures into the pair-generation process Reimer  _et
al._ (2018); Imany  _et al._ (2019); Lu  _et al._ (2022); alternatively,
continuous biphoton spectra can be shaped into bins using external cavities Lu
_et al._ (2018, 2020) or programmable frequency filters Sandoval  _et al._
(2019); Chen  _et al._ (2022). Here we adopt the last of these configurations
by introducing a Fourier transform pulse shaper (Finisar Waveshaper 4000B) to
carve out pairs of energy-correlated frequency bins: Fig. 1(b) highlights a
qubit example with two bins in the L-band ($I_{0}$ and $I_{1}$) and two in the
C-band ($S_{0}$ and $S_{1}$). Each bin is designed with a width of 18 GHz and
spaced 25 GHz apart. Additionally, the same pulse shaper splits the signal and
idler photons into separate optical fibers for subsequent state
characterization. Assuming entanglement in both DoFs, the ideal entangled two-
photon state can be expressed as follows:

$\ket{\Psi_{PF}}=\ket{\Psi_{P}}\otimes\ket{\Psi_{F}}\\\
=(\alpha\ket{HH}+\beta\ket{VV})\otimes\sum_{k=0}^{d-1}\gamma_{k}\ket{\omega_{k}^{(I)}\omega_{d-1-k}^{(S)}},\vspace{-0.1in}$
(1)

where $\omega_{k}^{(I)}=\omega_{0}^{(I)}+k\Delta\omega$ and
$\omega_{d-1-k}^{(I)}=\omega_{d-1}^{(S)}-k\Delta\omega$ denote frequencies
carrying photons chosen from the spectrum such that
$\omega_{0}^{(I)}+\omega_{d-1}^{(S)}=\omega_{p}$ (the pump frequency). The
broadband nature resulting from the type-0 phase matching condition ideally
yields $|\gamma_{k}|\approx\frac{1}{\sqrt{d}}$ for all bins of interest. Any
nonuniformity can be rectified via the pulse shaper Bernhard  _et al._ (2013)
(not required in our experiments).

For a hyperentangled photon pair with $N$ DoFs, each having encoding levels
$d_{1},d_{2},...,d_{N}$, a complete state reconstruction typically requires
$\mathcal{O}\left(\prod_{i=1}^{N}d_{i}^{4}\right)$ linearly independent local
projections across these $N$ DoFs, as the minimum number of independent
parameters describing a mixed state scales quadratically with the Hilbert
space dimension Bertlmann and Krammer (2008). For instance, a total of
$2^{4}\times 3^{4}=1296$ measurements were performed in Ref. Barreiro  _et
al._ (2005) to achieve full tomography of a $(d_{1}\otimes
d_{1})\otimes(d_{2}\otimes d_{2})=(2\otimes 2)\otimes(3\otimes 3)$ system,
which involved a pair of polarization qubits and orbital angular momentum
qutrits. An experimentally simpler approach involves QST of each DoF
independently, using semidefinite programming to compute a lower bound on the
global state fidelity from information in reduced density matrices Chen  _et
al._ (2020, 2022). This approach reduces the number of measurements required
for state characterization to
$\mathcal{O}\left(\sum_{i=1}^{N}d_{i}^{4}\right)$ (a sum over DoFs rather than
a product) , particularly valuable when scaling to higher dimensions. However,
it is crucial to note that these measurements alone cannot provide a complete
characterization of hyperentangled states.

We follow the first approach by performing local projections simultaneously in
both DoFs. For polarization tomography, we employ two motorized polarization
analyzers comprising free-space collimators, a quarter-wave plate, a half-wave
plate, and a polarizing beamsplitter, with a total throughput of $\sim$70%.
Frequency measurement is a more intricate process. Previous methods in the
context of polarization–frequency hyperentanglement Chen  _et al._ (2022);
Francesconi  _et al._ (2023) have relied on Hong–Ou–Mandel interference to
verify frequency entanglement and then infer the density matrix under specific
state assumptions Ramelow  _et al._ (2009). However, for slow (integrating)
detectors, this approach is sensitive to two-dimensional frequency-bin
entanglement only, where the signal and idler photons also share identical
spectra Lingaraju  _et al._ (2019). For more general state analysis applicable
to high-dimensional and nondegenerate photons, we utilize an electro-optic
phase modulator (EOM; EOspace) and wavelength-selective switch (WSS; Finisar)
to implement the necessary projective measurements Lu  _et al._ (2018);
Sandoval  _et al._ (2019); Lu  _et al._ (2020, 2022) for the frequency DoF. It
is crucial to execute the frequency-bin projection _after_ the polarization
analyzer due to the polarization sensitivity of the phase modulation; in this
way, each EOM receives a fixed physical polarization state throughout and can
provide consistent modulation and transmission efficiencies.

To achieve full reconstruction of the actual global state $\rho_{PF}$ with
qubit frequency encoding [$d=2$ in Eq. (1)], we employ a total of $16\times 8$
local projections across the two DoFs: i.e., 16 polarization projections James
_et al._ (2001) paired with 8 frequency-resolved measurements (4 projections
in the Pauli $Z\otimes Z$ basis, and 4 in the $X\otimes X$ basis). In the
$Z\otimes Z$ measurement, the EOMs are turned off and the two WSSs—one for the
signal photon (C-band WSS) and one for the idler photon (L-band
WSS)—demultiplex each photon by color, amounting to a $2\times 2$ joint
spectral intensity measurement. The $X\otimes X$ measurement is equivalent to
preceding the $Z\otimes Z$ measurement with two parallel Hadamard operations,
which we implement probabilistically by driving the EOMs at the bin spacing
(25 GHz) with a modulation index of $1.435$ rad Lu  _et al._ (2020). While we
could implement a more traditional tomographically complete measurement set in
frequency bins, the strong correlations observed in these two mutually
unbiased bases (MUBs) are sufficient to provide a clear entanglement witness
Coles  _et al._ (2017) that, when combined with Bayesian inference, can be
used to estimate the full quantum state with low uncertainty Lu  _et al._
(2018). Significantly, this situation is possible with a completely uniform
prior distribution—i.e., with no a priori restrictions on the form of the
state itself—as Bayes’ theorem automatically extracts all the information
available in a given measurement set via its logical framework MacKay (2003).
We consequently leverage this feature to streamline the number of measurement
settings for Bayesian QST Blume-Kohout (2010); Lukens  _et al._ (2020) below.

Figure 2: Bayesian mean density matrices (real parts) of the $(2\otimes
2)_{P}\otimes(2\otimes 2)_{F}$ hyperentangled state [arrows in Fig. 1(b)] and
its reduced states in individual DoFs. The imaginary components (not shown)
are less than 0.005 for all density matrices.

Figure 2(a) presents the Bayesian mean density matrix $\rho_{PF}$ for a
channel pair with signal and idler photons centered at 1554.2 nm and 1577.0
nm, respectively, obtained from 128 projections with 60 s integration time per
point. Computing the fidelity of 1024 density matrix samples with the target
state [Eq. (1) with $\alpha=\beta$ and $\gamma_{0}=\gamma_{1}$], we find
$\operatorname{\mathcal{F}}_{PF}=94.4(6)\%$—a value that simulations suggest
is constrained by the number of coincidence events, rather than the quality of
the states.We also evaluate the reduced states in polarization
($\rho_{P}=\Tr_{F}\rho_{PF}$) and frequency ($\rho_{F}=\Tr_{P}\rho_{PF}$) by
computing the respective partial traces. The density matrices and fidelities
with respect to $\ket{\Psi_{F}}$ and $\ket{\Psi_{P}}$ in Eq. (1) are presented
in Fig. 2(b,c). Quantitatively, we can lower- and upper-bound the distillable
entanglement per DoF with the coherent information $I_{C}$ (maximized over
one-way communication direction) Devetak and Winter (2005) and logarithmic
negativity $E_{\mathcal{N}}$ Vidal and Werner (2002), respectively. From the
QST data, we obtain the intervals $[I_{C},E_{\mathcal{N}}]=[0.69(3),0.936(9)]$
ebits for $\rho_{P}$ and $[I_{C},E_{\mathcal{N}}]=[0.76(2),0.954(5)]$ ebits
for $\rho_{F}$, confirming clear usable entanglement in both DoFs. Replicating
the tomographic procedure for four additional channel pairs spanning the
spectrum, denoted by stars in Fig. 1(b), we measure the following
hyperentangled state fidelities (counting outward from the spectral center),
denoted as ordered pairs $(\mathcal{F}_{PF},\mathcal{F}_{P},\mathcal{F}_{F})$:
$(93.3(7),94.5(6),95.9(4))\%$, $(93.3(7),94.5(6),96.1(4))\%$,
$(93.7(8),94.8(7),96.7(3))\%$, and $(91.3(9),93.1(8),94.8(5))\%$. Such a
sampling indicates that true hyperentanglement persists across the bandwidth
as expected.

Figure 3: Bayesian mean density matrices (real parts) of the $(2\otimes
2)_{P}\otimes(3\otimes 3)_{F}$ hyperentangled state and its reduced states in
individual DoFs. The imaginary components (not shown) are less than 0.05 for
all density matrices.

Given the broadband nature of the generated photon pairs, expanding the
frequency dimensions is straightforward. For instance, to create frequency
qutrits ($d=3$), we can simply consider three pairs of 25 GHz-spaced, 18 GHz-
wide bins, ideally resulting in the state
$\ket{\Psi_{PF}}=\ket{\Psi_{P}}\otimes\frac{1}{\sqrt{3}}\left(\ket{\omega_{0}^{(I)}\omega_{2}^{(S)}}+\ket{\omega_{1}^{(I)}\omega_{1}^{(S)}}+\ket{\omega_{2}^{(I)}\omega_{0}^{(S)}}\right)$.
However, as the system’s dimensionality increases, so does the number of
measurements required for QST. In the interests of speed, we take advantage of
the high degree of polarization entanglement and consider measurements in the
$Z\otimes Z$ and $X\otimes X$ MUBs only (the same pair explored for frequency
bins in the qubit example), thus reducing the number of polarization
measurements from 16 to 8. To characterize high-dimensional frequency-bin
entanglement, we leverage a novel method based on random measurements. This
involves applying random phases with the pulse shaper and random frequency
mixing with the EOMs, followed by computational-basis measurements Lu  _et
al._ (2022). We consider a total of 720 measurements (60 s each): 8
polarization projections, 10 different random EOM and shaper settings
(motivated by findings in Lu  _et al._ (2022)), and $d\times d=9$ signal-idler
frequency-bin combinations. Specifically, for each setting, the pulse shaper
applies $2d=6$ random spectral phases between $0$ and $2\pi$ to the
aforementioned frequency bins, and both EOMs receive a sinusoidal voltage with
amplitude $\delta$ chosen randomly between $0$ and $2.32$ radians Lu  _et al._
(2022). The resulting mean density matrix [Fig. 3(a)] shows high-fidelity
hyperentanglement: $\operatorname{\mathcal{F}}_{PF}=90.8(7)\%$ in the
$(2\otimes 2)_{P}\otimes(3\otimes 3)_{F}$ system. We again compute the reduced
states in both DoFs [Fig. 3(b,c)], with measured distillable entanglement
intervals of $[I_{C},E_{\mathcal{N}}]=[0.62(1),0.915(3)]$ and
$[I_{C},E_{\mathcal{N}}]=[1.04(4),1.48(1)]$ ebits for the polarization and
frequency DoFs, respectively, to be compared with the maximum qubit limit of 1
ebit and qutrit limit of 1.58 ebits.

Looking ahead, we see no immediate obstacles to generating even higher-
dimensional hyperentangled states by expanding the frequency dimension, which
in principle is capped only by the ratio of total bandwidth to the bin
spacing. Nevertheless, introducing more spectral content may degrade
polarization entanglement due to increased sensitivity to polarization-mode
dispersion (PMD) Gordon and Kogelnik (2000). This could introduce undesired
polarization–frequency correlations that will ultimately depend on the
specific fiber channel. To augment dimensionality without increasing bandwidth
(and hence protecting the state from PMD impairments), narrower frequency
spacings can be pursued instead, yet it is important to note that commercial
diffractive pulse shapers and WSSs usually have resolutions $\gtrsim$10 GHz.
Ultimately, the optimal solution could involve a fully integrated version of
the Sagnac source Unt , supporting polarization diversity and allowing direct
definition of frequency bins through the optical resonances of microrings.
Importantly, in any case where PMD limits the usable frequency-bin dimension
for an _individual_ state $\ket{\Psi_{PF}}$, the remaining bandwidth can still
be leveraged for _parallelization_ , in which the output is sliced into
subbands that are each sufficiently narrow to evade PMD degradation but
collectively utilize the entire band. Such a source could be used for
wavelength-multiplexed entanglement distribution, but where the intra-channel
frequency-bin entanglement carries explicit quantum information as well as the
polarization DoF.

Finally, the ability to _manipulate_ the expanded Hilbert space is the key to
fully harnessing the potential of hyperentangled states. For example,
controlled unitaries between polarization and frequency DoFs Xu  _et al._
(2023) will be valuable for implementing hyperentanglement-based versions of
protocols such as dense coding Barreiro  _et al._ (2008), superdense
teleportation Graham  _et al._ (2015); Chapman  _et al._ (2020), and
entanglement distillation Simon and Pan (2002); Ecker  _et al._ (2021).
However, the EOMs utilized in our measurement scheme are polarization-
sensitive, which limits their suitability for certain applications in this
context. While we can evade this restriction in QST by placing the EOMs after
polarization projections, such a simplification will not be feasible for
general quantum operations in the joint Hilbert space. Further advances in
polarization-diverse/insensitive frequency modulation techniques Ikuta  _et
al._ (2018); Sandoval  _et al._ (2019) will therefore prove valuable to fully
utilize such states in multidimensional quantum information processing.

## Acknowledgments

The authors would like to thank B. T. Kirby for valuable discussions. This
work was performed in part at Oak Ridge National Laboratory, operated by UT-
Battelle for the U.S. Department of energy under contract no. DE-
AC05-00OR22725. Funding was provided by the U.S. Department of Energy, Office
of Science, Advanced Scientific Computing Research (Field work Proposals
ERKJ378, ERKJ353), and the National Science Foundation (2034019-ECCS).

## References

  * Kwiat (1997) P. G. Kwiat, J. Mod. Opt. 44, 2173 (1997).
  * Barreiro  _et al._ (2005) J. T. Barreiro, N. K. Langford, N. A. Peters, and P. G. Kwiat, Phys. Rev. Lett. 95, 260501 (2005).
  * Barreiro  _et al._ (2008) J. T. Barreiro, T.-C. Wei, and P. G. Kwiat, Nat. Phys. 4, 282 (2008).
  * Williams  _et al._ (2017) B. P. Williams, R. J. Sadlier, and T. S. Humble, Phys. Rev. Lett. 118, 050501 (2017).
  * Reimer  _et al._ (2018) C. Reimer, S. Sciara, P. Roztocki, M. Islam, L. R. Cortés, Y. Zhang, B. Fischer, S. Loranger, R. Kashyap, A. Cino, S. T. Chu, B. E. Little, D. J. Moss, L. Caspani, W. J. Munro, J. Azaña, M. Kues, and R. Morandotti, Nat. Phys. 15, 148 (2018).
  * Imany  _et al._ (2019) P. Imany, J. A. Jaramillo-Villegas, M. S. Alshaykh, J. M. Lukens, O. D. Odele, A. J. Moore, D. E. Leaird, M. Qi, and A. M. Weiner, npj Quantum Inf. 5, 59 (2019).
  * Ecker  _et al._ (2021) S. Ecker, P. Sohr, L. Bulla, M. Huber, M. Bohmann, and R. Ursin, Phys. Rev. Lett. 127, 040506 (2021).
  * Chen  _et al._ (2022) C. Chen, C. Xu, A. Riazi, E. Y. Zhu, A. C. Greenwood, A. V. Gladyshev, P. G. Kazansky, B. T. Kirby, and L. Qian, Phys. Rev. A 105, 043702 (2022).
  * Francesconi  _et al._ (2023) S. Francesconi, A. Raymond, R. Duhamel, P. Filloux, A. Lemaître, P. Milman, M. Amanti, F. Baboux, and S. Ducci, Photon. Res. 11, 270 (2023).
  * Simon and Pan (2002) C. Simon and J.-W. Pan, Phys. Rev. Lett. 89, 257901 (2002).
  * Kues  _et al._ (2019) M. Kues, C. Reimer, J. M. Lukens, W. J. Munro, A. M. Weiner, D. J. Moss, and R. Morandotti, Nat. Photon. 13, 170 (2019).
  * Lu  _et al._ (2019) H.-H. Lu, A. M. Weiner, P. Lougovski, and J. M. Lukens, IEEE Photon. Technol. Lett. 31, 1858 (2019).
  * Alshowkan  _et al._ (2022) M. Alshowkan, J. M. Lukens, H.-H. Lu, B. T. Kirby, B. P. Williams, W. P. Grice, and N. A. Peters, Opt. Lett. 47, 6480 (2022).
  * Li  _et al._ (2005) X. Li, P. L. Voss, J. E. Sharping, and P. Kumar, Phys. Rev. Lett. 94, 053601 (2005).
  * Fan  _et al._ (2007) J. Fan, M. Eisaman, and A. Migdall, Phys. Rev. A 76, 043836 (2007).
  * Vergyris  _et al._ (2017) P. Vergyris, F. Kaiser, E. Gouzien, G. Sauder, T. Lunghi, and S. Tanzilli, Quantum Sci. Technol. 2, 024007 (2017).
  * Lu  _et al._ (2018) H.-H. Lu, J. M. Lukens, N. A. Peters, B. P. Williams, A. M. Weiner, and P. Lougovski, Optica 5, 1455 (2018).
  * Lu  _et al._ (2022) H.-H. Lu, K. V. Myilswamy, R. S. Bennink, S. Seshadri, M. S. Alshaykh, J. Liu, T. J. Kippenberg, D. E. Leaird, A. M. Weiner, and J. M. Lukens, Nat. Commun. 13, 4338 (2022).
  * Lu  _et al._ (2020) H.-H. Lu, E. M. Simmerman, P. Lougovski, A. M. Weiner, and J. M. Lukens, Phys. Rev. Lett. 125, 120503 (2020).
  * Sandoval  _et al._ (2019) O. E. Sandoval, N. B. Lingaraju, P. Imany, D. E. Leaird, M. Brodsky, and A. M. Weiner, Opt. Lett. 44, 1674 (2019).
  * Bernhard  _et al._ (2013) C. Bernhard, B. Bessire, T. Feurer, and A. Stefanov, Phys. Rev. A 88, 032322 (2013).
  * Bertlmann and Krammer (2008) R. A. Bertlmann and P. Krammer, J. Phys. A: Math. Theor. 41, 235303 (2008).
  * Chen  _et al._ (2020) Y. Chen, S. Ecker, J. Bavaresco, T. Scheidl, L. Chen, F. Steinlechner, M. Huber, and R. Ursin, Phys. Rev. A 101, 032302 (2020).
  * Ramelow  _et al._ (2009) S. Ramelow, L. Ratschbacher, A. Fedrizzi, N. Langford, and A. Zeilinger, Phys. Rev. Lett. 103, 253601 (2009).
  * Lingaraju  _et al._ (2019) N. B. Lingaraju, H.-H. Lu, S. Seshadri, P. Imany, D. E. Leaird, J. M. Lukens, and A. M. Weiner, Opt. Express 27, 38683 (2019).
  * James  _et al._ (2001) D. F. James, P. G. Kwiat, W. J. Munro, and A. G. White, Phys. Rev. A 64, 052312 (2001).
  * Coles  _et al._ (2017) P. J. Coles, M. Berta, M. Tomamichel, and S. Wehner, Rev. Mod. Phys. 89, 015002 (2017).
  * MacKay (2003) D. J. C. MacKay, “Probabilities and inference,” in  _Information Theory, Inference, and Learning Algorithms_ (Cambridge University Press, Cambridge, UK, 2003).
  * Blume-Kohout (2010) R. Blume-Kohout, New J. Phys. 12, 043034 (2010).
  * Lukens  _et al._ (2020) J. M. Lukens, K. J. H. Law, A. Jasra, and P. Lougovski, New J. Phys. 22, 063038 (2020).
  * Devetak and Winter (2005) I. Devetak and A. Winter, Proc. R. Soc. A 461, 207 (2005).
  * Vidal and Werner (2002) G. Vidal and R. F. Werner, Phys. Rev. A 65, 032314 (2002).
  * Gordon and Kogelnik (2000) J. Gordon and H. Kogelnik, Proc. Natl. Acad. Sci. 97, 4541 (2000).
  * (34) Manuscript in preparation.
  * Xu  _et al._ (2023) D. Xu, C. Chen, B. T. Kirby, and L. Qian, arXiv preprint arXiv:2304.01470 (2023).
  * Graham  _et al._ (2015) T. M. Graham, H. J. Bernstein, T.-C. Wei, M. Junge, and P. G. Kwiat, Nat. Commun. 6, 7185 (2015).
  * Chapman  _et al._ (2020) J. C. Chapman, T. M. Graham, C. K. Zeitler, H. J. Bernstein, and P. G. Kwiat, Phys. Rev. Appl. 14, 014044 (2020).
  * Ikuta  _et al._ (2018) R. Ikuta, T. Kobayashi, T. Kawakami, S. Miki, M. Yabuno, T. Yamashita, H. Terai, M. Koashi, T. Mukai, T. Yamamoto, and N. Imoto, Nat. Commun. 9, 1997 (2018).