<table>
  <thead>
    <tr>
      <th>Step</th>
      <th>Algorithm</th>
      <th>Description</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>1</td>
      <td>Algorithm 4: Solution algorithm for the 2D ACRP</td>
      <td></td>
    </tr>
    <tr>
      <td>2</td>
      <td>Input: A, ̂θ, ̂v, q, q̄, θ, θ̄, ε</td>
      <td></td>
    </tr>
    <tr>
      <td>3</td>
      <td>Output: q*, θ*, LB, UB</td>
      <td></td>
    </tr>
    <tr>
      <td>4</td>
      <td>if MIQP 3.4 is infeasible then</td>
      <td></td>
    </tr>
    <tr>
      <td>5</td>
      <td>Return INFEASIBLE</td>
      <td></td>
    </tr>
    <tr>
      <td>6</td>
      <td>if q is feasible (no speed violation) then</td>
      <td></td>
    </tr>
    <tr>
      <td>7</td>
      <td>UB ← LB, q* ← q, θ* ← θ</td>
      <td></td>
    </tr>
    <tr>
      <td>8</td>
      <td>else</td>
      <td></td>
    </tr>
    <tr>
      <td>9</td>
      <td>CONVERGED ← False</td>
      <td></td>
    </tr>
    <tr>
      <td>10</td>
      <td>while CONVERGED = False do</td>
      <td></td>
    </tr>
    <tr>
      <td>11</td>
      <td>q', θ', UB-NLP ← Solve Model 3.2 as NLP with fixed z</td>
      <td></td>
    </tr>
    <tr>
      <td>12</td>
      <td>if UB-NLP &lt; UB then</td>
      <td></td>
    </tr>
    <tr>
      <td>13</td>
      <td>UB ← UB-NLP, q* ← q', θ* ← θ'</td>
      <td></td>
    </tr>
    <tr>
      <td>14</td>
      <td>for i ∈ A do</td>
      <td></td>
    </tr>
    <tr>
      <td>15</td>
      <td>if δ<sup>2</sup><sub>i,x</sub> + δ<sup>2</sup><sub>i,y</sub> &gt; q̄<sup>2</sup><sub>i</sub> then</td>
      <td></td>
    </tr>
    <tr>
      <td>16</td>
      <td>Add constraint (3.22b)</td>
      <td></td>
    </tr>
    <tr>
      <td>17</td>
      <td>else if δ<sup>2</sup><sub>i,x</sub> + δ<sup>2</sup><sub>i,y</sub> &lt; q̄<sup>2</sup><sub>i</sub> then</td>
      <td></td>
    </tr>
    <tr>
      <td>18</td>
      <td>if δ̃<sub>i,x</sub> &gt; δ<sup>2</sup><sub>i,x</sub> then</td>
      <td></td>
    </tr>
    <tr>
      <td>19</td>
      <td>Add segment to Γ<sub>i,x</sub> at δ<sub>i,x</sub>, variable s<sup>k</sup><sub>i,x</sub>, constraints (3.37a), (3.38a), and update (3.39a)</td>
      <td></td>
    </tr>
    <tr>
      <td>20</td>
      <td>if δ̃<sub>i,y</sub> &gt; δ<sup>2</sup><sub>i,y</sub> then</td>
      <td></td>
    </tr>
    <tr>
      <td>21</td>
      <td>Add segment to Γ<sub>i,y</sub> at δ<sub>i,y</sub>, variable s<sup>k</sup><sub>i,y</sub>, constraints (3.37b), (3.38b), and update (3.39b)</td>
      <td></td>
    </tr>
    <tr>
      <td>22</td>
      <td>q, θ, δ<sub>x</sub>, δ<sub>y</sub>, δ̃<sub>x</sub>, δ̃<sub>y</sub>, z, LB ← Solve MIQCP 3.5 and calculate q and θ from δ<sub>x</sub> and δ<sub>y</sub></td>
      <td></td>
    </tr>
    <tr>
      <td>23</td>
      <td>if MIQCP 3.5 is infeasible then</td>
      <td></td>
    </tr>
    <tr>
      <td>24</td>
      <td>Return INFEASIBLE</td>
      <td></td>
    </tr>
    <tr>
      <td>25</td>
      <td>if q is feasible (no speed violation) then</td>
      <td></td>
    </tr>
    <tr>
      <td>26</td>
      <td>UB ← LB, q* ← q, θ* ← θ</td>
      <td></td>
    </tr>
    <tr>
      <td>27</td>
      <td>CONVERGED ← True</td>
      <td></td>
    </tr>
    <tr>
      <td>28</td>
      <td>if (UB-LB)/UB ≤ ε then</td>
      <td></td>
    </tr>
    <tr>
      <td>29</td>
      <td>CONVERGED ← True</td>
      <td></td>
    </tr>
  </tbody>
</table>

### 3.4.3 Decomposition algorithm for the 2D+FL ACRP

We next introduce a two-step decomposition approach for the nonconvex lexicographic 2D+FL conflict resolution problem represented by Model 3.3. The first objective function (3.31) focuses on minimising the number of FL changes. This objective function is null and minimal if all aircraft can remain at their initial FL. Since minimising aircraft FL re-assignment is the highest priority, we need only to identify combinations of aircraft which are non-separable in 2D and ensure that such combinations are not assigned to the same FL. Let $\Omega_I \subseteq 2^A$ be the
set of aircraft combinations which are 2D non-separable, i.e. for any $\omega \in \Omega_I$, the subset of aircraft $\omega$ cannot be separated in 2D. We require:

$$ \sum_{i \in \omega} \rho_{ik} \leq |\omega| - 1, \quad \forall \omega \in \Omega_I, \forall k \in \bigcap_{i \in \omega} \mathcal{Z}_i. \qquad (3.41) $$

Let $\Delta\rho_i \ge 0$ be a variable representing the absolute FL deviation for aircraft $i \in \mathcal{A}$. Objective function (3.31) can be linearised using traditional mathematical programming techniques. Combining the FL separation constraint (3.41) and the linearised objective function yields a compact FL assignment formulation summarized in Model 3.6 which is a MILP with an exponential number of constraints.

**Model 3.6 (FL assignment formulation).**

$$ \min \sum_{i \in \mathcal{A}} \Delta\rho_i \qquad (3.42a) $$

s.t.

$$ \sum_{k \in \mathcal{Z}_i} \rho_{ik} = 1, \quad \forall i \in \mathcal{A}, \qquad (3.42b) $$

$$ \sum_{i \in \omega} \rho_{ik} \leq |\omega| - 1, \quad \forall \omega \in \Omega_I, \forall k \in \bigcap_{i \in \omega} \mathcal{Z}_i, \qquad (3.42c) $$

$$ \Delta\rho_i \geq \sum_{k \in \mathcal{Z}_i} k\rho_{ik} - \hat{\rho}_i, \quad \forall i \in \mathcal{A}, \qquad (3.42d) $$

$$ \Delta\rho_i \geq \hat{\rho}_i - \sum_{k \in \mathcal{Z}_i} k\rho_{ik}, \quad \forall i \in \mathcal{A}, \qquad (3.42e) $$

$$ \rho_{ik} \in \{0, 1\}, \quad \forall i \in \mathcal{A}, k \in \mathcal{Z}_i, \qquad (3.42f) $$

$$ \Delta\rho_i \geq 0, \quad \forall i \in \mathcal{A}. \qquad (3.42g) $$

We propose to solve Model 3.6 by first restricting Constraint (3.41) to subsets of size two, i.e. 2D non-separable aircraft pairs, which can be efficiently identified using Algorithm 3 in pre-processing. We then decompose the 2D+FL problem into a series of 2D conflict resolution problems, one per FL, based on the optimal solution $\rho^*$ of the relaxed Model 3.6. For each FL, we solve the corresponding 2D problem using the exact constraint generation approach of Algorithm 4 with the aircraft set $\mathcal{A}_k = \{i \in \mathcal{A} : \rho_{ik}^* = 1\}$. If Algorithm 4 returns INFEASIBLE for FL $k$, then the corresponding FL separation constraint (3.41) with $\omega = \mathcal{A}_k$ is generated and Model 3.6 is re-solved with the additional constraint(s). The process is repeated until all 2D problems are feasible. Let $\mathcal{Z}$ be the set of all FLs, i.e. $\mathcal{Z} = \cup_{k \in \mathcal{Z}_i : i \in \mathcal{A}} \mathcal{Z}_i$, we denote $q_k$ and $\theta_k$ the vectors of speed and heading controls for aircraft assigned to FL $k \in \mathcal{Z}$, respectively. The proposed approach for the lexicographic 2D+FL conflict resolution problem is summarized in Algorithm 5. Note that this decomposition approach offers the possibility to solve all 2D problems in parallel via the **for** loop at Line 8.

## 3.5 Robust aircraft conflict resolution problem

In this section, we introduce formulations for the robust aircraft conflict resolution problem. We first define the uncertainty model before discussing how it can be incorporated within the ACRP. We then propose a tractable robust counterpart formulation for the robust ACRP.
<caption>**Algorithm 5:** Solution algorithm for the lexicographic 2D+FL conflict resolution problem</caption>

**Input:** $\mathcal{A}$, $[\mathcal{Z}_i]_{i \in \mathcal{A}}$, $\hat{\theta}$, $\hat{v}$, $\hat{\rho}$, $\underline{q}$, $\underline{q}$, $\underline{\theta}$, $\underline{\theta}$
**Output:** $\underline{q}^*$, $\underline{\theta}^*$, $\underline{\rho}^*$
1 &nbsp;&nbsp; $\mathcal{P}, \mathcal{P}_F, \mathcal{P}_S, \mathcal{P}_I \leftarrow \text{Algorithm 3}$
2 &nbsp;&nbsp; $\Omega_I \leftarrow \{(i, j) \in \mathcal{P}_I : \mathcal{Z}_i \cap \mathcal{Z}_j \neq \emptyset\}$
3 &nbsp;&nbsp; $\mathcal{Z} \leftarrow \bigcup_{k \in \mathcal{Z}_i : i \in \mathcal{A}} \mathcal{Z}_k$
4 &nbsp;&nbsp; $\text{CONVERGED} \leftarrow \text{False}$
5 &nbsp;&nbsp; **while** $\text{CONVERGED} = \text{False}$ **do**
6 &nbsp;&nbsp;&nbsp;&nbsp; $\underline{\rho}^* \leftarrow \text{Solve MILP 3.6}$
7 &nbsp;&nbsp;&nbsp;&nbsp; $\text{CONVERGED} \leftarrow \text{True}$
8 &nbsp;&nbsp;&nbsp;&nbsp; **for** $k \in \mathcal{Z}$ **do**
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{A}_k \leftarrow \{i \in \mathcal{A} : \underline{\rho}_{ik}^* = 1\}$
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\underline{q}_k^*, \underline{\theta}_k^* \leftarrow \text{Algorithm 4 with } \mathcal{A} = \mathcal{A}_k$
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** Algorithm 4 returns **INFEASIBLE** **then**
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\Omega_I \leftarrow \Omega_I \cup \mathcal{A}_k$
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\text{CONVERGED} \leftarrow \text{False}$
14 &nbsp;&nbsp;&nbsp;&nbsp; $\underline{q}^*, \underline{\theta}^* \leftarrow [\underline{q}_k^*]_{k \in \mathcal{Z}}, [\underline{\theta}_k^*]_{k \in \mathcal{Z}}$

### 3.5.1 Uncertainty model

We assume that each aircraft has a source of randomness and this affects its current velocity and position. Let $\epsilon_i = [\epsilon_{i,x}, \epsilon_{i,y}]^\top$ be a vector of random variables representing the uncertainty on the velocity components of aircraft $i \in \mathcal{A}$. We next use this vector of random variables to define aircraft-based uncertainty sets.

Given aircraft $i \in \mathcal{A}$, let $v_{i,x}$ and $v_{i,y}$ its velocity components.

$$v_{i,x} = q_i \hat{v}_i \cos(\hat{\theta}_i + \theta_i), \qquad (3.43)$$

$$v_{i,y} = q_i \hat{v}_i \sin(\hat{\theta}_i + \theta_i). \qquad (3.44)$$

**Definition 3.4 (Uncertainty set of aircraft).** The uncertainty set of aircraft $i \in \mathcal{A}$, is defined as:

$$\mathcal{U}_i = \{ \epsilon_i \in \mathbb{R}^2 \mid -\bar{\epsilon}_{i,x} \le \epsilon_{i,x} \le \bar{\epsilon}_{i,x}, -\bar{\epsilon}_{i,y} \le \epsilon_{i,y} \le \bar{\epsilon}_{i,y} \}, \qquad (3.45)$$

where $\bar{\epsilon}_{i,x} \ge 0$ and $\bar{\epsilon}_{i,y} \ge 0$ represent the maximum perturbations on the velocity components $v_{i,x}$ and $v_{i,y}$, respectively, of aircraft $i$.

We denote $[\tilde{v}_{i,x}, \tilde{v}_{i,y}]^\top$ the vector of random aircraft velocity components where the random variables $\tilde{v}_{i,x}$ and $\tilde{v}_{i,y}$ take values in $\tilde{v}_{i,x} \in [-v_{i,x}(1 + \bar{\epsilon}_{i,x}), v_{i,x}(1 + \bar{\epsilon}_{i,x})]$ and $\tilde{v}_{i,y} \in [-v_{i,y}(1 + \bar{\epsilon}_{i,y}), v_{i,y}(1 + \bar{\epsilon}_{i,y})]$, respectively. Accordingly, for each pair of aircraft $(i, j) \in \mathcal{P}$, the random relative velocity components $\tilde{v}_{ij,x}$ and $\tilde{v}_{ij,y}$ are:

$$\tilde{v}_{ij,x} = v_{i,x}(1 + \epsilon_{i,x}) - v_{j,x}(1 + \epsilon_{j,x}) = v_{ij,x} + v_{i,x}\epsilon_{i,x} - v_{j,x}\epsilon_{j,x}, \qquad (3.46a)$$

$$\tilde{v}_{ij,y} = v_{i,y}(1 + \epsilon_{i,y}) - v_{j,y}(1 + \epsilon_{j,y}) = v_{ij,y} + v_{i,y}\epsilon_{i,y} - v_{j,y}\epsilon_{j,y}. \qquad (3.46b)$$
Let $\tilde{v}_{ij,x}$, $\tilde{v}_{ij,y}$ and $\tilde{v}_{ij,y}$, $\tilde{v}_{ij,x}$ be the lower and upper bounds for $\tilde{v}_{ij,x}$ and $\tilde{v}_{ij,y}$, respectively. These bounds can be determined using Eq. (3.46) and the bounds on speed and heading control provided in Eqs. (3.12) and (3.13). The derived bounds on the random relative velocity components can be used to define the random relative velocity box which characterizes all possible trajectories for the pair $(i, j) \in \mathcal{P}$ under uncertainty.

**Definition 3.5** (Random relative velocity box). Consider a pair of aircraft $(i, j) \in \mathcal{P}$. Let $\mathcal{U}_i$ and $\mathcal{U}_j$ be the uncertainty sets of aircraft i and j, respectively. Let $\tilde{\mathcal{B}}_{ij}(\mathcal{U}_i, \mathcal{U}_j)$ be the parametric subset of $\mathbb{R}^2$ defined as

$$ \tilde{\mathcal{B}}_{ij}(\mathcal{U}_i, \mathcal{U}_j) \equiv \{ (\tilde{v}_{ij,x}, \tilde{v}_{ij,y}) \in \mathbb{R}^2 : \underline{v}_{ij,x} \le \tilde{v}_{ij,x} \le \bar{v}_{ij,x}, \underline{v}_{ij,y} \le \tilde{v}_{ij,y} \le \bar{v}_{ij,y} \}. \quad (3.47) $$

$\tilde{\mathcal{B}}_{ij}(\mathcal{U}_i, \mathcal{U}_j)$ is the random relative velocity box of $(i, j) \in \mathcal{P}$ under the uncertainty sets $\mathcal{U}_i$ and $\mathcal{U}_j$.

<figure><img src="image_7.png" /><figcaption>Figure 3.4: Illustration of a two-aircraft conflict in the plane $\{(v_{ij,x}, v_{ij,y}) \in \mathbb{R}^2\}$. The inner box with black lines corresponds to the velocity bounds $\mathcal{B}_{ij}$ in the deterministic scenario while the box with purple dots $\tilde{\mathcal{B}}_{ij}(\mathcal{U}_i, \mathcal{U}_j)$ is the random relative velocity box. The region is hashed in red corresponds to the conflict region $\mathcal{C}_{ij}$. If the uncertainty sets $\mathcal{U}_i$ and $\mathcal{U}_j$ of aircraft i and j are empty, i.e. aircraft trajectories are deterministic, then aircraft i and j are conflict-free. In turn, if the uncertainty sets $\mathcal{U}_i$ and $\mathcal{U}_j$ are such that the random relative velocity $\tilde{\mathcal{B}}_{ij}$ intersects with the conflict region $\mathcal{C}_{ij}$, then there exists a risk of conflict.</figcaption></figure>

The random relative velocity box is illustrated in Figure 3.4 for a two-aircraft conflict. To ensure that a pair of aircraft $(i, j) \in \mathcal{P}$ is separated for any realization of the random variables $\epsilon_i \in \mathcal{U}_i$ and $\epsilon_j \in \mathcal{U}_j$, we redefine (3.17e)-(3.17h) using the random velocity variables $\tilde{v}_{ij,x}$ and $\tilde{v}_{ij,y}$:

$$ \tilde{v}_{ij,y} \hat{x}_{ij} - \tilde{v}_{ij,x} \hat{y}_{ij} \le 0, \quad \text{if } z_{ij} = 1, \quad \forall (i, j) \in \mathcal{P}, \quad (N_1) $$

$$ \tilde{v}_{ij,y} \hat{x}_{ij} - \tilde{v}_{ij,x} \hat{y}_{ij} \ge 0, \quad \text{if } z_{ij} = 0, \quad \forall (i, j) \in \mathcal{P}, \quad (N_0) $$

$$ \tilde{v}_{ij,y} \gamma_{ij}^l - \tilde{v}_{ij,x} \phi_{ij}^l \le 0, \quad \text{if } z_{ij} = 1, \quad \forall (i, j) \in \mathcal{P}, \quad (S_1) $$

$$ \tilde{v}_{ij,y} \gamma_{ij}^u - \tilde{v}_{ij,x} \phi_{ij}^u \ge 0, \quad \text{if } z_{ij} = 0, \quad \forall (i, j) \in \mathcal{P}. \quad (S_0) $$

To find robust aircraft trajectories under uncertainty sets $\mathcal{U}_i$ and $\mathcal{U}_j$, we require that the pairwise separation constraints $(N_0)$, $(N_1)$, $(S_1)$ and $(S_0)$ hold for any $(\tilde{v}_{ij,x}, \tilde{v}_{ij,y}) \in \tilde{\mathcal{B}}_{ij}(\mathcal{U}_i, \mathcal{U}_j)$. We next use state-of-the-art approaches in robust optimization to integrate these constraints into a robust counterpart formulation.
### 3.5.2 Robust counterpart formulation

The constraints $(N_0)$, $(N_1)$, $(S_1)$ and $(S_0)$ are function of the random variables $\tilde{v}_{ij,x}$ and $\tilde{v}_{ij,y}$. We use the approach of Bertsimas and Sim [2004] to reformulate these robust separation constraints as integer-linear constraints with regards to deterministic relative velocity variables $v_{ij,x}$ and $v_{ij,y}$. The constraints $(N_0)$, $(N_1)$, $(S_1)$ and $(S_0)$ are of the form $a\tilde{v}_{ij,x} + b\tilde{v}_{ij,y} \le 0$ (omitting the disjunction), and can be rearranged by separating deterministic and random elements as follows:

$$a(v_{i,x} - v_{j,x}) - b(v_{i,y} - v_{j,y}) + av_{i,x}\epsilon_{i,x} - av_{j,x}\epsilon_{j,x} - bv_{i,y}\epsilon_{i,y} + bv_{j,y}\epsilon_{j,y} \le 0. \quad (3.49)$$

We introduce new variables $\nu_{i,x} \ge 0$ and $\nu_{i,y} \ge 0$ for each $i \in A$ through constraints (3.50) to impose artificial bounds on aircraft velocity components $v_{i,x}$ and $v_{i,y}$:

$$-\nu_{i,x} \le v_{i,x} \le \nu_{i,x}, \quad \forall i \in A, \quad (3.50a)$$

$$-\nu_{i,y} \le v_{i,y} \le \nu_{i,y}, \quad \forall i \in A. \quad (3.50b)$$

Let $RS \equiv \{N_1, N_0, S_1, S_0\}$ be a set of indices corresponding to constraints $(N_0)$, $(N_1)$, $(S_1)$, $(S_0)$, respectively. Further, let $\alpha_k = 1$ if $k = N_1$ or $k = S_1$, and let $\alpha_k = 0$ if $k = N_0$ or $k = S_0$. Let $\Gamma$ be a real parameter that takes values in the range $[0,4]$, where the upper bound is given by the number of decision variables in constraint (3.49). The parameter $\Gamma$ determines the level of robustness for each robust separation constraint. To link the level of robustness $\Gamma$ with each robust separation constraint $k \in RS$ and each aircraft pair $(i, j) \in P$, we introduce real decision variables $\psi_{ij}^k \ge 0$. Further, each constraint of the form (3.49) involves four decision variables, hence for each constraint $k \in RS$ and $(i, j) \in P$, we introduce associated variables $\rho_{ij,x}^{l,k} \ge 0$ and $\rho_{ij,y}^{l,k} \ge 0$ for $l \in \{i, j\}$. These artificial variables reflect the consumption of robustness resources of each velocity variable in (3.49). The following constraints link variables $\psi_{ij}^k$, $\rho_{ij,x}^{l,k}$ and $\rho_{ij,x}^{l,k}$ with $\nu_i^x$ and $\nu_j^x$:

$$\psi_{ij}^k + \rho_{ij,x}^{l,k} \ge \nu_{i,x}\bar{\epsilon}_{i,x}, \quad \text{if } z_{ij} = \alpha_k, \quad \forall (i, j) \in P, \forall l \in \{i, j\}, \forall k \in RS, \quad (3.51a)$$

$$\psi_{ij}^k + \rho_{ij,y}^{l,k} \ge \nu_{i,y}\bar{\epsilon}_{i,y}, \quad \text{if } z_{ij} = \alpha_k, \quad \forall (i, j) \in P, \forall l \in \{i, j\}, \forall k \in RS. \quad (3.51b)$$

Constraints of the form (3.49) can then be rewritten as:

$$a(v_{i,x} - v_{j,x}) - b(v_{i,y} - v_{j,y}) + \psi_{ij}^k \Gamma + \sum_{l \in \{i,j\}} (\rho_{ij,x}^{l,k} + \rho_{ij,y}^{l,k}) \le 0, \quad \text{if } z_{ij} = \alpha_k, \forall (i, j) \in P, \forall k \in RS. \quad (3.52)$$

Combining the above constraints and artificial real variables, we obtain a tractable formulation of the robust separation constraints.

**Proposition 3.4.** The set of robust separation constraints $(N_0)$, $(N_1)$, $(S_1)$ and $(S_0)$ are
equivalent to the following set of integer-linear constraints and real variables:

$$a(v_{i,x} - v_{j,x}) - b(v_{i,y} - v_{j,y}) + \psi_{ij}^k \Gamma + \sum_{l \in \{i,j\}} (\rho_{ij,x}^{l,k} + \rho_{ij,y}^{l,k}) \le 0, \quad \text{if } z_{ij} = \alpha_k, \quad \forall (i,j) \in \mathcal{P}, \forall k \in \mathcal{R}\mathcal{S}, \quad (3.53\text{a})$$

$$\psi_{ij}^k + \rho_{ij,x}^{l,k} \ge \nu_{i,x} \bar{\epsilon}_{i,x}, \quad \text{if } z_{ij} = \alpha_k, \quad \forall (i,j) \in \mathcal{P}, \forall l \in \{i,j\}, \forall k \in \mathcal{R}\mathcal{S}, \quad (3.53\text{b})$$

$$\psi_{ij}^k + \rho_{ij,y}^{l,k} \ge \nu_{i,y} \bar{\epsilon}_{i,y}, \quad \text{if } z_{ij} = \alpha_k, \quad \forall (i,j) \in \mathcal{P}, \forall l \in \{i,j\}, \forall k \in \mathcal{R}\mathcal{S}, \quad (3.53\text{c})$$

$$\rho_{ij,x}^{l,k}, \rho_{ij,y}^{l,k} \ge 0, \quad \text{if } z_{ij} = \alpha_k, \quad \forall (i,j) \in \mathcal{P}, \forall l \in \{i,j\}, \forall k \in \mathcal{R}\mathcal{S}, \quad (3.53\text{d})$$

$$\psi_{ij}^k \ge 0, \quad \text{if } z_{ij} = \alpha_k, \quad \forall (i,j) \in \mathcal{P}, \forall k \in \mathcal{R}\mathcal{S}, \quad (3.53\text{e})$$

$$- \nu_{i,x} \le v_{i,x} \le \nu_{i,x}, \quad \forall i \in \mathcal{A}, \quad (3.53\text{f})$$

$$- \nu_{i,y} \le v_{i,y} \le \nu_{i,y}, \quad \forall i \in \mathcal{A}, \quad (3.53\text{g})$$

$$\nu_{i,x}, \nu_{i,y} \ge 0, \quad \forall i \in \mathcal{A}. \quad (3.53\text{h})$$

Proposition 3.4 establishes the equivalency between constraints $(N_0)$, $(N_1)$, $(S_1)$ and $(S_0)$, which are expressed in terms of random relative velocity variables $\tilde{v}_{ij,x}$ and $\tilde{v}_{ij,y}$, and a tractable integer-linear reformulation of these constraints which uses additional real variables. Using this reformulation, we propose the following robust counterpart formulation of the ACRP.

**Model 3.7 (Robust counterpart nonconvex formulation of the 2D ACRP).**

$$\min \sum_{i \in \mathcal{A}} (1-w)(1-q_i)^2 + w\theta_i^2, \quad (3.54\text{a})$$

s.t.

Motion Equations (3.2),

Robust Separation Constraints and Variables (3.53),

Speed and Heading Control Constraints (3.12), (3.13),

$$v_{ij,x}, v_{ij,y} \in \mathcal{B}_{ij}, \quad \forall (i,j) \in \mathcal{P}, \quad (3.54\text{b})$$

$$z_{ij} \in \{0, 1\}, \quad \forall (i,j) \in \mathcal{P}, \quad (3.54\text{c})$$

$$q_i \le q_i \le \bar{q}_i, \quad \forall i \in \mathcal{A}, \quad (3.54\text{d})$$

$$\theta_i \le \theta_i \le \bar{\theta}_i, \quad \forall i \in \mathcal{A}. \quad (3.54\text{e})$$

The formulation presented in Model 3.7 is nonconvex due to trigonometric functions and non-linear components. To solve the robust optimization problem represented by Model 3.7, we use the complex number formulation of the ACRP summarized in Model 3.2 and adapt it for the robust ACRP. Specifically, the complex number reformulation of the robust is constructed by substituting the original motion equations (3.2) with the reformulated
aircraft velocities constraints (3.55):

$$v_{i,x} = \delta_{i,x} \hat{v}_i \cos(\hat{\theta}_i) - \delta_{i,y} \hat{v}_i \sin(\hat{\theta}_i), \quad (3.55a)$$

$$v_{i,y} = \delta_{i,y} \hat{v}_i \cos(\hat{\theta}_i) - \delta_{i,x} \hat{v}_i \sin(\hat{\theta}_i). \quad (3.55b)$$

The speed and heading control constraints and variables bounds are replaced with (3.22), (3.23) and (3.21), respectively; and the original objective function (3.17a) is replaced with (3.24). Observe that the robust separation constraints (3.53) remain unchanged since aircraft velocity variables $v_{i,x}$ and $v_{i,y}$ are linked to variables $\delta_{i,x}$ and $\delta_{i,y}$ via constraints (3.55). We refer to the resulting formulation as the robust complex number formulation and use Algorithm 4 to solve this problem to optimality.

## 3.6 Numerical experiments

We first introduce the experimental framework used to test the proposed mixed-integer formulations and algorithms in Section 3.6.1. We then explore the behaviour of the proposed 2D objective function in Section 3.6.2. Numerical results for the 2D ACRP are presented in Section 3.6.3, and results for the 2D+FL problem are presented in Section 3.6.4.

### 3.6.1 Experiments design

We test the performance of the proposed mixed-integer formulations and algorithms using four benchmarking problems from the literature: the Circle Problem (CP), the Flow Problem (FP), the Grid Problem (GP) and the Random Circle Problem (RCP). The four types of benchmarking instances are illustrated in Figure 3.5. The CP consists of a set of aircraft uniformly positioned on the circle heading towards its centre. Aircraft speeds are assumed to be identical, hence the problem is highly symmetric (see Figure 3.5a). The CP is notoriously difficult due to the geometry of aircraft initial configuration and has been widely used for benchmarking CD&R algorithms in the literature [Durand and Alliot, 2009, Rey et al., 2015, Cafieri and Omheni, 2017, Cafieri and Rey, 2017, Rey and Hijazi, 2017]. To break the symmetry of CP benchmarking instances, Vanaret et al. [2012] introduced the RCP which builds on the same framework, but aircraft initial speeds and headings are randomly deviated within specified ranges to create randomised instances with less structure (see Figure 3.5b). CP and RCP instances are named CP-N and RCP-N, respectively, where N is the total number of aircraft. Lehouillier et al. [2017] formally introduced two additional structured problems which aim to represent more realistic air traffic configurations: the FP and the GP. The FP consists of two streams of aircraft separated by an angle $\alpha$ and anchored on the circumference of a circle. In each stream, aircraft are separated by at least 5 NM from each other (see Figure 3.5c). The GP consists of two FP instances separated by 15 NM diagonally (see Figure 3.5d). In our experiments, on each stream of aircraft in FP and GP instances, consecutive aircraft are initially separated by 15 NM. FP and GP instances are named FP-N and GP-N, respectively, where N denotes the number of aircraft per stream.

In all experiments, we use a circle of radius 200 NM. For CP, FP and GP instances, all aircraft have an initial speed of 500 NM/h. For RCP instances, aircraft initial speeds are randomly chosen in the range 486-594 NM/h and their initial headings are deviated from a radial trajectory (i.e. towards the centre of the circle) by adding a randomly chosen angle
<figure><img src="image_8.png" /><figcaption>Figure 3.5: Example of 2D benchmarking instances for the Circle Problem (CP), Random Circle Problem (RCP), Flow Problem (FP) and Grid Problem (GP).</figcaption></figure>

between $-\frac{\pi}{6}$ and $+\frac{\pi}{6}$. For FP and GP instances, we use $\alpha = \frac{\pi}{6}$ and $\alpha = \frac{\pi}{2}$, respectively, as proposed by Lehouillier et al. [2017].

We consider a subliminal speed control range of [-6%, +3%] [Bonini et al., 2009]. We consider two heading control ranges, we first assume that aircraft can modify their heading within the range [-30°, +30°] as commonly used in the literature [Cafieri and Omheni, 2017, Rey and Hijazi, 2017], and we also consider a reduced heading control range of [-15°, +15°]. For conflict resolution problems with altitude control, we randomly assign each aircraft $i \in \mathcal{A}$ to a FL $\hat{\rho}_i \in \mathcal{Z}$, and we assume that only adjacent FLs are available for aircraft, i.e. $\mathcal{Z}_i = \{\hat{\rho}_i - 1, \hat{\rho}_i, \hat{\rho}_i + 1\}$.

The proposed approach is referred to as Disjunctive and is compared to two benchmarks from the literature: the method proposed by Rey and Hijazi [2017] named Disjunctive-2017 and an implementation based on the so-called shadow separation constraints named Shadow. The method Disjunctive corresponds to Algorithm 4 for 2D ACRP instances and to Algorithm 5 for 2D+FL ACRP instances. The method Disjunctive-2017 is based on the same disjunctive linear separation conditions as Disjunctive but uses the algorithm proposed by Rey and Hijazi [2017]. This algorithm has 3 steps: i) solve the MIQP relaxation (Model 3.4), ii) solve the MIQCP relaxation (Model 3.5) without any mixed-integer cuts (3.37), and, if a feasible solution has not been obtained, iii) solve Model 3.2 as NLP for fixed $\mathcal{Z}$. Compared to Disjunctive, the method Disjunctive-2017 is identical in the first step only. In step ii) of Disjunctive-2017, the MIQCP relaxation is solved without any constraint generation, thus all convex quadratic cuts (3.34) are added for all aircraft simultaneously and no mixed-integer cuts are generated. In addition, if the NLP fails to find a feasible solution in step iii), then the method fails to yield a feasible solution whereas the method Disjunctive is guaranteed to converge to a global optimal solution if one exists. In the implementation of the method Shadow, the control variables and the algorithm are identical to that of the method Disjunctive, and the only difference between both methods is the set of separation constraints used, i.e. Constraints (3.10) and (3.11) are replaced with the shadow separation conditions and the number of binary variables required to express these On/Off constraints. Specifically, for the 2D ACRP, Disjunctive only requires a single binary variable per aircraft pair ($z_{ij}$) whereas Shadow requires four binary variables per aircraft pair [Pallottino et al., 2002, Alonso-Ayuso et al., 2011, 2016]. All three methods are implemented using the same pre-processing procedure (Algorithm 3) to eliminate conflict-free aircraft pairs.

All 2D ACRPs are solved by implementing Algorithm 4 with an optimality gap $\epsilon = 0.01$
and a time limit of 10 minutes. All 2D+FL problems are solved by implementing Algorithm 5 which calls Algorithm 4 using the same optimality gap, and we allocate a time limit of 10 minutes per FL. All models are implemented using PYTHON on a personal computer with 16 GB of RAM and an Intel i7 processor at 2.9GHz. The MIQPs and MIQCPs are solved with CPLEX v12.10 [CPLEX, 2009] API for PYTHON using default options.

We next conduct a sensitivity analysis on the preference weight *w* in the 2D objective function to explore its impact on aircraft trajectories in Section 3.6.2. We present results on 2D problems in Section 3.6.3 and results on 2D+FL problems in Section 3.6.4.

### 3.6.2 Sensitivity analysis on the preference weight *w*

To quantify the impact of the preference weight *w* in the proposed 2D objective function (3.24), we conduct numerical experiments on one instance of each of the four types of benchmarking instances for varying values of *w*. For this experiment, we focus on the typical heading control range [-30°, +30°]. We report the total speed deviation $\Sigma_q = \sum_{i \in \mathcal{A}} (1 - q_i)^2$, and the total heading deviation $\Sigma_{\theta} = \sum_{i \in \mathcal{A}} \theta_i^2$. Our goal is to show that by varying the preference weight $w \in ]0, 1[$, the decision-maker can control the desired level of trade-off between total speed deviation and total heading deviation. Recall that in objective function (3.24), *w* is the coefficient of $\delta_{i,y}^2 = (q_i \sin(\theta_i))^2$ which is minimal for $\theta_i = 0$; while $(1-w)$ is the coefficient of $(1 - \delta_{i,x})^2 = (1 - q_i \cos(\theta_i))^2$ which is minimal for $q_i = 1$ and $\theta_i = 0$. Hence, one can expect that increasing (resp. decreasing) *w* will tend to penalise heading (resp. speed) deviations more than speed (resp. heading) deviations.

This behaviour is confirmed in our numerical experiments. Specifically, we solve the 2D instances CP-8, FP-10, GP-10 and one RCP-30 instance for $w = 0.1, \dots, 0.9$ in steps of size 0.1, i.e. for a total of 9 values of *w* per instance. All instances are solved to optimality using Algorithm 4 with no MIQCP iterations, i.e. the MIQP returned a global optimal solution for all tests. The change in the total speed deviation $\Sigma_q$ and in the total heading deviation $\Sigma_{\theta}$ are reported in Figure 3.6. For all four instances tested, we find that increasing *w* monotonically decreases the total heading deviation and monotonically increases the total speed deviation. Further, we observe that in all cases both the total speed and heading deviations are of similar order of magnitudes.

This sensitivity analysis shows that using the proposed 2D objective function, the decision-maker can control which manoeuvre is prioritised by scaling up or down the preference weight *w* accordingly. Higher values of *w* will minimise the total heading deviation while lower values of *w* will minimise the total speed deviation. We use $w = 0.5$ in the numerical experiments presented in the remaining of this study.

### 3.6.3 Results on 2D instances

Results for RCP instances are reported in Tables 3.1 and 3.2 for four instance sizes with 10, 20, 30 and 40 aircraft per group. For each instance size, 100 RCP instances are randomly generated. We compare the performance of the proposed formulations for both the standard and reduced heading control ranges.

Each row in the results tables represents a group of instances (RCP). The header of the results tables is presented from left to right. The left-most column, Instance, identifies
<figure><img src="image_9.png" /><figcaption>Figure 3.6: Sensitivity analysis on the preference weight $w$ of the 2D objective function (3.24). For all figures, $\Sigma_q$ represents the total speed deviation defined as $\sum_{i \in \mathcal{A}} (1 - q_i)^2$ (in red) and $\Sigma_{\theta}$ represents the total heading deviation defined as $\sum_{i \in \mathcal{A}} \theta_i^2$ (in blue).</figcaption></figure>

the instance; $|\mathcal{A}|$ is the number of aircraft and $n_c$ is the number of the conflicts. The next three columns summarize the performance of the pre-processing algorithm: $|\mathcal{P}_F|/|\mathcal{P}|$ is the proportion of conflict-free aircraft pairs; $|\mathcal{P}_I|/|\mathcal{P}|$ is the proportion of non-separable aircraft pairs; and Time is the runtime of Algorithm 3 in seconds. The next six columns summarize the performance of Disjunctive: LB and UB are the lower and upper bound; Gap is the optimality gap in percent calculated using LB and UB; Time is the total runtime in seconds; $n_i$ is the number of iterations of the while loop in Algorithm 4, and $n_t$ represents the proportion of instances that could not be solved within the time limit (10 minutes), i.e. the number of time-outs. The next five columns summarize the performance of Disjunctive-2017: $\Delta$UB is the upper bound of Disjunctive-2017 minus that of Disjunctive; Gap is the optimality gap, Time is the total runtime, $n_t$ is the proportion of time-outs, and Gain is the performance gain in runtime for instances solved within the time limit calculated as the runtime of Disjunctive-2017 minus that of Disjunctive in percentage: a positive value indicates that Disjunctive is faster. The right-most five columns summarize the performance of Shadow relative to Disjunctive, similarly to the five previous columns. For RCP instances, we also report the number of times each method failed to find a feasible solution in columns named $n_f$. This value is always zero for CP, FP and GP instances hence we do not report it in their corresponding results tables.

The experiments performed on RCP instances reveal that the pre-processing procedure (Algorithm 3) can eliminate approximately 8% of aircraft pairs when using a reduced heading control range (Table 3.2), whereas no aircraft pairs are conflict-free using a standard heading control range (Table 3.1). The implementation of Disjunctive on RCP instances (Tables 3.1 and 3.2) reveals that while all 10- and 20-aircraft RCP instances can be solved via the MIQP iteration, 30- and 40-aircraft RCP instances may require additional MIQCP iterations. Since no MIQCP iterations are performed ($n_i = 0$) for RCP-10 and RCP-20, the results obtained
using Disjunctive-2017 are identical to those obtained using Disjunctive. Using the standard heading control range (Table 3.1), on average, RCP-10 instances can be solved in 0.05 s and 0.14 s using Disjunctive and Shadow, respectively. RCP-20 instances require 0.26 s and 2.52 s, on average, using Disjunctive and Shadow, respectively, under the standard heading control range. The performance of the proposed formulations on RCP-10 and RCP-20 using a reduced heading control range (Table 3.2) is of similar order of magnitude due to the relatively low average number of conflicts per instance, i.e. 3.1 and 13.1, respectively. Those values lead to a performance gain of Disjunctive over Shadow of 71.8% and 62.6% for standard and reduced ranges, respectively, for RCP-10 instances; and of 89.6% and 81.3% for RCP-20 instances.

RCP-30 and RCP-40 instances have on average 32.9 and 59.3 conflicts, respectively, and present considerable computing challenges, notably the latter. Overall, we observe that reducing the heading control range tends to improve the performance while retaining comparable optimal solutions, as indicated by the similar UB values obtained. Disjunctive requires an average of 0.4 MIQCP iterations for RCP-30 instances using a standard heading control range (Table 3.1) compared to 1.4 when using a reduced heading control range (Table 3.2). These figures increase to 0.8 and 1.6 in RCP-40 instances. For RCP-30 instances, the average optimality gaps are 0.17%, 4.60% and 2.36% using Disjunctive, Disjunctive-2017 and Shadow, respectively using the standard heading control range (Table 3.1). These figures are reduced to 0.01%, 3.13% and 0.89% using the reduced heading control range (Table 3.2). For RCP-40 instances, the average optimality gaps values are significantly greater: 15%, 51.2% and 29.8% for a standard heading control (Table 3.1); and 13.1%, 52.0% and 25.8% for a reduced heading control range (Table 3.2). This shows that Shadow tends to double the optimality gap compared to Disjunctive, while Disjunctive-2017 tends to triple this figure compared to Disjunctive. Using a standard heading control range (Table 3.1), Disjunctive is able to solve all but 3% of the RCP-30 instances whereas Disjunctive-2017 and Shadow time out on 26% and 20% of these instances, respectively. The performance gain of Disjunctive compared to Disjunctive-2017 and Shadow in terms of runtime are 21.0% and 41.7%, respectively. For RCP-40 instances, Disjunctive can solve 28% of the problems (72% of time-outs), while Disjunctive-2017 and Shadow time out on 73% and 100% of the instances, respectively. We also observe that Disjunctive-2017 and Shadow may occasionally fail to find a UB as competitive as that found by the Disjunctive on RCP-40 instances. The runtime performance gains of Disjunctive on RCP-30 instances are 21.0% and 41.7% compared to Disjunctive-2017 and Shadow, respectively, using standard heading control range; and increase to 57.2% and 60.5% using the reduced heading control range. For RCP-40 instances, the gains of Disjunctive compared to Disjunctive-2017 are 3.73% and 24.6% for standard and reduced heading control ranges, respectively; while Shadow is unable to solve any RCP-40 instances to optimality (100% of time-outs). The method Disjunctive-2017 fails to find a feasible solution for 19% and 48% of the RCP-40 instances using the standard and reduced heading control ranges, respectively.

### 3.6.4 Results on 2D+FL instances

Results on 2D+FL instances are reported using similar tables as in Section 3.6.3. Instead of reporting the number of aircraft and the number of conflicts, we report the average number of aircraft and conflicts per FL, i.e. $\frac{|A|}{|Z|}$, $\frac{n_c}{|Z|}$ respectively. In addition, we include a section corresponding to the FL assignment formulation. In all our numerical experiments on 2D+FL instances, only a single pass through the **while** loop of Algorithm 5 is required. Hence, we only report Obj which is the objective function value of MILP 3.6; and Time which is
<table><thead><tr><th rowspan="2">Instance</th><th rowspan="2">|A|</th><th rowspan="2">n<sub>c</sub></th><th colspan="3">Pre-processing</th><th colspan="7">Disjunctive</th></tr><tr><th>|P<sub>F</sub>|/|P|</th><th>|P<sub>I</sub>|/|P|</th><th>Time</th><th>LB</th><th>UB</th><th>Gap</th><th>Time</th><th>n<sub>i</sub></th><th>n<sub>t</sub></th><th>n<sub>f</sub></th></tr></thead><tbody><tr><td>RCP-10</td><td>10</td><td>3.10</td><td>0</td><td>0</td><td>0.12</td><td>2.2E-4</td><td>2.2E-4</td><td>0.00</td><td>0.05</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-20</td><td>20</td><td>13.1</td><td>0</td><td>0</td><td>0.23</td><td>1.7E-3</td><td>1.7E-3</td><td>0.00</td><td>0.26</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-30</td><td>30</td><td>32.9</td><td>0</td><td>0</td><td>0.45</td><td>7.1E-3</td><td>7.1E-3</td><td>0.17</td><td>135</td><td>0.4</td><td>3</td><td>0</td></tr><tr><td>RCP-40</td><td>40</td><td>59.3</td><td>0</td><td>0</td><td>0.60</td><td>1.8E-2</td><td>2.4E-2</td><td>15.0</td><td>516</td><td>0.8</td><td>72</td><td>0</td></tr></tbody></table>

<table><caption>Table 3.1: Results for 2D RCP instances with a speed control range of [−6%, +3%] and a heading control range of [−30°, +30°]. Times (Time) are reported in seconds. The proportions of conflict-free (|P<sub>F</sub>|/|P|) and non-separable (|P<sub>I</sub>|/|P|) pairs, optimality gaps (Gap), time-outs (n<sub>t</sub>), infeasible solutions (n<sub>f</sub>) and the performance gain (Gain) are reported in %.</caption><thead><tr><th rowspan="2">Instance</th><th colspan="7">Disjunctive-2017</th><th colspan="7">Shadow</th></tr><tr><th>ΔUB</th><th>Gap</th><th>Time</th><th>n<sub>t</sub></th><th>n<sub>f</sub></th><th>Gain</th><th>ΔUB</th><th>Gap</th><th>Time</th><th>n<sub>t</sub></th><th>n<sub>f</sub></th><th>Gain</th></tr></thead><tbody><tr><td>RCP-10</td><td>0.00</td><td>0.00</td><td>0.05</td><td>0</td><td>0</td><td>0.00</td><td>0.00</td><td>0.02</td><td>0.14</td><td>0</td><td>0</td><td>71.8</td></tr><tr><td>RCP-20</td><td>0.00</td><td>0.00</td><td>0.26</td><td>0</td><td>0</td><td>0.00</td><td>0.00</td><td>0.01</td><td>2.52</td><td>0</td><td>0</td><td>89.6</td></tr><tr><td>RCP-30</td><td>0.00</td><td>4.60</td><td>171</td><td>26</td><td>0</td><td>21.0</td><td>0.00</td><td>2.36</td><td>231</td><td>20</td><td>0</td><td>41.7</td></tr><tr><td>RCP-40</td><td>0.01</td><td>51.2</td><td>536</td><td>77</td><td>19</td><td>3.73</td><td>0.02</td><td>29.8</td><td>600</td><td>100</td><td>0</td><td>-</td></tr></tbody></table>

the corresponding computing runtime in seconds. For numerical experiments on 2D+FL problems, we focus on RCP instances named RCP-N-Z where N is the number of aircraft and Z is the number of FLs. We consider three numbers of aircraft: 50, 100 and 150; and two numbers of FLs: 3 and 5. In terms of conflict density, for the same number of aircraft, instances with 3 FLs have a greater number of conflicts compared to instances with 5 FLs and are more computationally challenging. The results are reported in Table 3.3 for the standard heading control range, and in Table 3.4 for the reduced heading control range.

Using the standard heading control range ([-30°, +30°], see Table 3.3), the numerical experiments reveal that among all six groups of 2D+FL instances, only some RCP-150-3 and RCP-150-5 instances may require aircraft to change FLs. That is, for all other 2D+FL instances, all conflicts can be resolved using 2D trajectory control only and thus aircraft do not require performing any FL change. This can be explained by observing that for RCP-150-3 instances, the average number of aircraft per FL is 50 which corresponds to a denser aircraft configuration compared to the 2D RCP-40 instances which all admit feasible solutions. In comparison, only a single RCP-150-5 instance required a FL change for a single aircraft. For RCP-150-5, even though the average density per FL is 30, some instances may have denser FL requiring FL separation. Among all 100 RCP-150-3 instances, 35% of the instances requires a FL change. The maximum number of non-separable pairs (|P<sub>I</sub>|) is 2 and the maximum objective value of Model 3.6 is 1, indicating that only a single aircraft deviated from its initial FL.

Comparing the methods, we find that Disjunctive slightly outperforms the two benchmarks in terms of solution quality by occasionally finding better UBs than Disjunctive-2017 or Shadow. For instances with 3 FLs, the optimality gaps are relatively small, i.e. less than 1.0% for RCP-50-3 and RCP-100-3 using Disjunctive and around 10% for RCP-150-3 However, using Disjunctive-2017 and Shadow, the optimality gaps are considerably larger. Especially
<table><thead><tr><th rowspan="2">Instance</th><th rowspan="2">|A|</th><th rowspan="2">n<sub>c</sub></th><th colspan="2">Pre-processing</th><th colspan="7">Disjunctive</th></tr><tr><th>|P<sub>F</sub>|/|P|</th><th>|P<sub>I</sub>|/|P|</th><th>Time</th><th>LB</th><th>UB</th><th>Gap</th><th>Time</th><th>n<sub>i</sub></th><th>n<sub>t</sub></th><th>n<sub>f</sub></th></tr></thead><tbody><tr><td>RCP-10</td><td>10</td><td>3.10</td><td>8.2</td><td>0</td><td>0.33</td><td>2.2E-4</td><td>2.2E-4</td><td>0.03</td><td>0.04</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-20</td><td>20</td><td>13.1</td><td>7.7</td><td>0</td><td>0.51</td><td>1.7E-3</td><td>1.7E-3</td><td>0.01</td><td>0.24</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-30</td><td>30</td><td>32.9</td><td>7.7</td><td>0</td><td>0.62</td><td>7.2E-3</td><td>7.2E-3</td><td>0.01</td><td>66.3</td><td>1.4</td><td>3</td><td>0</td></tr><tr><td>RCP-40</td><td>40</td><td>59.3</td><td>7.9</td><td>0</td><td>0.75</td><td>1.8E-2</td><td>2.2E-2</td><td>13.1</td><td>389</td><td>1.6</td><td>59</td><td>0</td></tr></tbody></table>

<table><caption>Table 3.2: Results for 2D RCP instances with a speed control range of [−6%, +3%] and a heading control range of [−15°, +15°]. Times (Time) are reported in seconds. The proportions of conflict-free (|P<sub>F</sub>|/|P|) and non-separable (|P<sub>I</sub>|/|P|) pairs, optimality gaps (Gap), time-outs (n<sub>t</sub>), infeasible solutions (n<sub>f</sub>) and the performance gain (Gain) are reported in %.</caption><thead><tr><th rowspan="2">Instance</th><th colspan="7">Disjunctive-2017</th><th colspan="7">Shadow</th></tr><tr><th>ΔUB</th><th>Gap</th><th>Time</th><th>n<sub>t</sub></th><th>n<sub>f</sub></th><th>Gain</th><th>ΔUB</th><th>Gap</th><th>Time</th><th>n<sub>t</sub></th><th>n<sub>f</sub></th><th>Gain</th></tr></thead><tbody><tr><td>RCP-10</td><td>0.00</td><td>0.03</td><td>0.04</td><td>0</td><td>0</td><td>0.00</td><td>0.00</td><td>0.03</td><td>0.11</td><td>0</td><td>0</td><td>62.6</td></tr><tr><td>RCP-20</td><td>0.00</td><td>0.01</td><td>0.24</td><td>0</td><td>0</td><td>0.00</td><td>0.00</td><td>0.01</td><td>1.39</td><td>0</td><td>0</td><td>81.3</td></tr><tr><td>RCP-30</td><td>0.00</td><td>3.13</td><td>155</td><td>22</td><td>0</td><td>57.2</td><td>0.00</td><td>0.89</td><td>167</td><td>10</td><td>0</td><td>60.5</td></tr><tr><td>RCP-40</td><td>0.01</td><td>52.0</td><td>516</td><td>78</td><td>48</td><td>24.6</td><td>0.02</td><td>25.8</td><td>600</td><td>100</td><td>0</td><td>-</td></tr></tbody></table>

for the former which yields an average gap of 15.7% for RCP-50-3 instances, and of 42.0% and 61.6% for RCP-100-3 and RCP-150-3 instances, respectively. Using Shadow, those values are less than 1.0% for RCP-50-3 instances, 53.2% for RCP-100-3 instances and 56.6% for RCP-150-3 instances. From a computational standpoint, the average runtime of Model 3.6 on RCP-150-3 instances is 4.56 s. For RCP-50-3 instances, Disjunctive, Disjunctive-2017 and Shadow can solve all instances in an average time of 0.88 s, 2.81 s and 10.2 s, respectively. For RCP-100-3 instances, which corresponds to an average number of 33.3 aircraft per FL, Disjunctive solves all instances in an average time of 3.3 minutes, whereas using Disjunctive-2017 there is a total of 64% instances that time out, and for the remaining instances that can be solved the runtime is 411 s. Shadow times out on 46% of the instances and requires an average runtime of 9.5 minutes for the instances solved. The denser RCP-150-3 instances, with an average of 50 aircraft per FL, present substantial computational challenges. Disjunctive is able to solve only 40% of these instances (60% of time-outs) and an average of 2.8 MIQCP iterations are required. In comparison, Disjunctive-2017 is unable to solve most of these instances (95% of time-outs) and fails to find a feasible solution for 5% of them. The method Shadow is unable to solve any of these instances within the available time limit. The performance gain using 3 FLs shows that a gain of 87.8% and 51.8% can be obtained using Disjunctive over Disjunctive-2017 and 91.3% and 61.9% using Disjunctive over Shadow for RCP-50-3 and RCP-100-3 instances, respectively.

Increasing the number of FLs from 3 to 5 reduces the density of aircraft per FL which translates into better computational performance for all three methods. The UB value obtained by those instances using Disjunctive-2017 is relatively close to the UB obtained with Disjunctive yielding deviations of 0.08, 0.09 and 0.10 in RCP-50-5, RCP-100-5 and RCP-150-5 instances, respectively (see Table 3.3). Shadow only deviates from the UB found by Disjunctive by 0.08 for RCP-150-5 instances. All RCP-50-5 and RCP-100-5 instances are
solved to optimality within the time limit using all three methods. The optimality gap is relatively small for Disjunctive and Shadow, and does not exceed 0.2% in those instances. Using Disjunctive-2017, it reaches 13.4% for RCP-50-5 instances and 45.5% for RCP-100-5 instances. For RCP-150-5 instances, using Disjunctive solves all problems with an optimality gap of 0.05%. However, Disjunctive-2017 fails to solve all RCP-150-5 instances and yields an average optimality gap of 45.1%. The method Shadow also fails to solve all RCP-150-5 instances and yields an average optimality gap of 1.56%. In terms of runtime, we observe that Disjunctive solves all instances with 5 FLs in less than 142 s (see Table 3.3). The performance of the methods on 2D+FL instances with 5 FLs reveal that gains of 75.6% and 42.7% can be obtained using Disjunctive over Disjunctive-2017, and 81.8% and 98.8% using Disjunctive over Shadow for RCP-50-5 and RCP-100-5 instances, respectively.

Reducing the heading control range (see Table 3.4), yields an overall similar performance. The main differences relative to the results obtained using the standard heading control range are observed in the pre-processing procedure which, as in the 2D RCP instances, eliminates on average 8% of the aircraft pairs. We also find that reducing the heading control range increases the number of non-separable pairs to 1.45% for RCP-150-3 instances and to 0.03% for RCP-150-5 instances. This increase in the proportion of non-separable aircraft pairs is reflected in the solution of the flight assignment formulation (MILP 3.6) which has an average value of 1.62 on RCP-150-3 instances and 0.02 on RCP-150-5 instances, respectively. For RCP-150-3, 67% of the instances required a FL change, the maximum number of non-separable pairs ($|P_1|$) is 6 and the maximum objective value of Model 3.6 is 3, indicating that three aircraft deviated from their initial FL. The methods find the same UBs for RCP-50-5 and RCP-100-5 instances, whereas for RCP-150-5 instances we observe an average deviation of 0.05 using Disjunctive-2017 and of 0.02 using Shadow. The optimality gap is below 1% for all instances using Disjunctive, while it ranges from 13.3% in RCP-50-3 instances to 40.1% in RCP-150-3 instances using Disjunctive-2017, and from 0.01% to 42.1% using Shadow. In terms of runtime, all three methods are able to solve all RCP-50-3 instances in on average 0.88 s, 2.37 s and 11.3 s using Disjunctive, Disjunctive-2017 and Shadow, respectively. For RCP-100-3 instances, we observe that all instances are solved by Disjunctive in 168 s, while 26% and 39% of them time out using Disjunctive-2017 and Shadow, respectively. In addition, Disjunctive-2017 fails to find a feasible solution in 10% of the instances. For those instances that can be solved, the average runtime is 474 s for Disjunctive-2017 and 453 s for Shadow. For RCP-150-3 instances, Disjunctive, Disjunctive-2017 and Shadow time out on 46%, 74% and 100% of the instances, respectively; and Disjunctive-2017 fails to find a feasible solution in 14% of the instances. For the remaining instances, the runtime is 420 s using Disjunctive and 587 s using Disjunctive-2017. The performance gains of Disjunctive compared to both benchmarks are above 62% for RCP-50-3 and RCP-100-3 instances, and near 40% compared to Disjunctive-2017 for RCP-150-3 instances. All RCP-50-5 and RCP-100-5 instances are solved to optimality by all methods. The runtime for those instances is less than 1 s for RCP-50-5 instances and 35.1 s for RCP-100-5 instances using Disjunctive, while Disjunctive-2017 and Shadow yield larger runtimes but of the same order of magnitude. These figures do not carry over to RCP-150-5 instances. Using Disjunctive, all RCP-150-5 instances can be solved with an average runtime of 152 s. However, using Disjunctive-2017 all instances time out and this method fails to find a feasible solution in 15% of the cases. Further, the average UB deviation is 0.12 using Disjunctive-2017 and the average optimality gap is 42.3%. Using Shadow, all RCP-150-5 instances also time out, the average UB deviation is 0.09 and the average optimality gap is 12.3%.
<table><thead><tr><th rowspan="2">Instance</th><th rowspan="2">|A|/|Z|</th><th rowspan="2">n<sub>c</sub>/|Z|</th><th colspan="2">Pre-processing</th><th colspan="2">FL assignment</th><th colspan="8">Disjunctive</th></tr><tr><th>|P<sub>F</sub>|/|P|</th><th>|P<sub>t</sub>|/|P|</th><th>Time</th><th>Obj.</th><th>Time</th><th>LB</th><th>UB</th><th>Gap</th><th>Time</th><th>n<sub>i</sub></th><th>n<sub>t</sub></th><th>n<sub>f</sub></th></tr></thead><tbody><tr><td>RCP-50-3</td><td>15</td><td>3.15</td><td>0</td><td>0.00</td><td>0.42</td><td>0.00</td><td>0.00</td><td>1.2E-2</td><td>1.3E-2</td><td>0.02</td><td>0.88</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-100-3</td><td>33</td><td>35.6</td><td>0</td><td>0.00</td><td>0.86</td><td>0.00</td><td>0.00</td><td>1.4E-1</td><td>1.4E-1</td><td>0.12</td><td>198</td><td>1.2</td><td>0</td><td>0</td></tr><tr><td>RCP-150-3</td><td>50</td><td>64.1</td><td>0</td><td>0.55</td><td>1.02</td><td>0.62</td><td>4.56</td><td>4.5E-2</td><td>5.2E-2</td><td>10.2</td><td>520</td><td>2.8</td><td>60</td><td>0</td></tr><tr><td>RCP-50-5</td><td>10</td><td>3.11</td><td>0</td><td>0.00</td><td>0.62</td><td>0.00</td><td>0.00</td><td>9.4E-3</td><td>9.4E-3</td><td>0.00</td><td>0.26</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-100-5</td><td>20</td><td>14.1</td><td>0</td><td>0.00</td><td>1.25</td><td>0.00</td><td>0.00</td><td>8.3E-3</td><td>8.3E-2</td><td>0.00</td><td>2.51</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-150-5</td><td>30</td><td>34.2</td><td>0</td><td>0.01</td><td>1.92</td><td>0.01</td><td>0.11</td><td>3.3E-1</td><td>3.3E-1</td><td>0.05</td><td>142</td><td>2.3</td><td>0</td><td>0</td></tr></tbody></table>

<table><caption>Table 3.3: Results for 2D+FL RCP instances with a speed control range of [−6%, +3%] and a heading control range of [−30°, +30°]. Times (Time) are reported in seconds. The proportions of conflict-free (|P<sub>F</sub>|/|P|) and non-separable (|P<sub>t</sub>|/|P|) pairs, optimality gaps (Gap), time-outs (n<sub>t</sub>), infeasible solutions (n<sub>f</sub>) and the performance gain (Gain) are reported in %.</caption><thead><tr><th rowspan="2">Instance</th><th colspan="6">Disjunctive-2017</th><th colspan="6">Shadow</th></tr><tr><th>ΔUB</th><th>Gap</th><th>Time</th><th>n<sub>t</sub></th><th>n<sub>f</sub></th><th>Gain</th><th>ΔUB</th><th>Gap</th><th>Time</th><th>n<sub>t</sub></th><th>n<sub>f</sub></th><th>Gain</th></tr></thead><tbody><tr><td>RCP-50-3</td><td>0.02</td><td>15.7</td><td>2.81</td><td>0</td><td>0</td><td>87.8</td><td>0.00</td><td>0.01</td><td>10.2</td><td>0</td><td>0</td><td>91.3</td></tr><tr><td>RCP-100-3</td><td>0.02</td><td>42.0</td><td>411</td><td>64</td><td>0</td><td>51.8</td><td>0.05</td><td>53.2</td><td>572</td><td>46</td><td>0</td><td>61.9</td></tr><tr><td>RCP-150-3</td><td>0.10</td><td>61.6</td><td>600</td><td>95</td><td>5</td><td>-</td><td>0.03</td><td>56.6</td><td>600</td><td>100</td><td>0</td><td>-</td></tr><tr><td>RCP-50-5</td><td>0.08</td><td>13.4</td><td>2.17</td><td>0</td><td>0</td><td>75.6</td><td>0.00</td><td>0.00</td><td>1.43</td><td>0</td><td>0</td><td>81.8</td></tr><tr><td>RCP-100-5</td><td>0.09</td><td>45.5</td><td>67.2</td><td>0</td><td>0</td><td>42.7</td><td>0.00</td><td>0.13</td><td>126</td><td>0</td><td>0</td><td>98.0</td></tr><tr><td>RCP-150-5</td><td>0.10</td><td>45.1</td><td>600</td><td>100</td><td>0</td><td>-</td><td>0.08</td><td>1.56</td><td>600</td><td>100</td><td>0</td><td>-</td></tr></tbody></table>

### 3.6.5 Results on the robust ACRP

In this section, we analyze the performance of the robust complex number formulation for a varying level of robustness (Γ) for a fixed maximum uncertainty of $\bar{\epsilon} = 5\%$. We present the results for RCP instances in Table 3.5 for 3 instance sizes with 10, 20 and 30 aircraft per group. For each instance group, 100 RCP instances are randomly generated and we report the average performance. For all instances, we compare the performance of the proposed formulations for varying level of robustness where $\Gamma = 0$ corresponds to the deterministic case and $\Gamma = 4$ corresponds to the most robust configuration.

Each row in the results tables represents a group of instances (RCP). The header of the results tables is presented from left to right: Γ is the level of robustness; UB is the objective function value obtained after solving Model 3.7; Gap is the optimality gap in percent obtained after solving Model 3.7, it is calculated directly by CPLEX; Time is the total runtime in seconds, n<sub>i</sub> is the number of MIQCP iterations where n<sub>i</sub> = 0 means that an optimal solution was found or the time limit was reached during the initial MIQP solve, n<sub>t</sub> represents the proportion of instances that could not be solved within the time limit (10 minutes), i.e. the number of time-outs and n<sub>0</sub> indicates the number/proportion of infeasible instances.

The results of the experiments on RCP instances are summarized in Table 3.5. Increasing the level of robustness in RCP instances also tends to increase the objective value (UB), and the effect is on average magnified on instances with a large number of aircraft. For RCP-10 instances, the average objective value for $\Gamma = 4$ is two orders of magnitude greater than that obtained using a deterministic configuration. For RCP-20 instances, the average objective values obtained are one order of magnitude greater than those obtained for RCP-10. We observe that RCP-10 instances can be solved within less than a second using for any level of
<table><thead><tr><th rowspan="2">Instance</th><th rowspan="2">|A|/|Z|</th><th rowspan="2">n<sub>c</sub>/|Z|</th><th colspan="2">Pre-processing</th><th colspan="2">FL assignment</th><th colspan="8">Disjunctive</th></tr><tr><th>|P<sub>F</sub>|/|P|</th><th>|P<sub>I</sub>|/|P|</th><th>Time</th><th>Obj.</th><th>Time</th><th>LB</th><th>UB</th><th>Gap</th><th>Time</th><th>n<sub>i</sub></th><th>n<sub>t</sub></th><th>n<sub>f</sub></th></tr></thead><tbody><tr><td>RCP-50-3</td><td>15</td><td>3.15</td><td>7.0</td><td>0.00</td><td>0.52</td><td>0.00</td><td>0.00</td><td>1.2E-2</td><td>1.2E-2</td><td>0.02</td><td>0.88</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-100-3</td><td>33</td><td>35.6</td><td>7.2</td><td>0.00</td><td>0.62</td><td>0.00</td><td>0.00</td><td>1.4E-1</td><td>1.4E-1</td><td>0.15</td><td>168</td><td>2.2</td><td>0</td><td>0</td></tr><tr><td>RCP-150-3</td><td>50</td><td>64.1</td><td>8.1</td><td>1.45</td><td>1.22</td><td>1.62</td><td>5.62</td><td>4.5E-2</td><td>5.2E-2</td><td>8.60</td><td>420</td><td>4.8</td><td>46</td><td>0</td></tr><tr><td>RCP-50-5</td><td>10</td><td>3.11</td><td>8.6</td><td>0.00</td><td>0.52</td><td>0.00</td><td>0.00</td><td>9.4E-3</td><td>9.4E-3</td><td>0.00</td><td>0.26</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-100-5</td><td>20</td><td>14.1</td><td>7.8</td><td>0.00</td><td>1.22</td><td>0.00</td><td>0.00</td><td>8.3E-3</td><td>8.3E-2</td><td>0.00</td><td>35.1</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-150-5</td><td>30</td><td>34.2</td><td>8.0</td><td>0.03</td><td>1.86</td><td>0.02</td><td>0.21</td><td>3.3E-1</td><td>3.3E-1</td><td>0.51</td><td>152</td><td>4.3</td><td>0</td><td>0</td></tr></tbody></table>

<table><caption>Table 3.4: Results for 2D+FL RCP instances with a speed control range of [−6%, +3%] and a heading control range of [−15°, +15°]. Times (Time) are reported in seconds. The proportions of conflict-free (|P<sub>F</sub>|/|P|) and non-separable (|P<sub>I</sub>|/|P|) pairs, optimality gaps (Gap), time-outs (n<sub>t</sub>), infeasible solutions (n<sub>f</sub>) and the performance gain (Gain) are reported in %.</caption><thead><tr><th rowspan="2">Instance</th><th colspan="6">Disjunctive-2017</th><th colspan="6">Shadow</th></tr><tr><th>ΔUB</th><th>Gap</th><th>Time</th><th>n<sub>t</sub></th><th>n<sub>f</sub></th><th>Gain</th><th>ΔUB</th><th>Gap</th><th>Time</th><th>n<sub>t</sub></th><th>n<sub>f</sub></th><th>Gain</th></tr></thead><tbody><tr><td>RCP-50-3</td><td>0.00</td><td>13.3</td><td>2.37</td><td>0</td><td>0</td><td>62.9</td><td>0.00</td><td>0.01</td><td>11.3</td><td>0</td><td>0</td><td>91.3</td></tr><tr><td>RCP-100-3</td><td>0.00</td><td>37.4</td><td>474</td><td>26</td><td>10</td><td>64.6</td><td>0.08</td><td>43.5</td><td>453</td><td>39</td><td>0</td><td>62.4</td></tr><tr><td>RCP-150-3</td><td>0.05</td><td>40.1</td><td>587</td><td>74</td><td>14</td><td>39.7</td><td>0.02</td><td>42.1</td><td>600</td><td>100</td><td>0</td><td>-</td></tr><tr><td>RCP-50-5</td><td>0.00</td><td>11.6</td><td>1.22</td><td>0</td><td>0</td><td>78.6</td><td>0.00</td><td>0.00</td><td>1.43</td><td>0</td><td>0</td><td>81.8</td></tr><tr><td>RCP-100-5</td><td>0.00</td><td>36.3</td><td>54.1</td><td>0</td><td>0</td><td>35.1</td><td>0.00</td><td>0.13</td><td>186</td><td>0</td><td>0</td><td>81.2</td></tr><tr><td>RCP-150-5</td><td>0.12</td><td>42.3</td><td>600</td><td>100</td><td>15</td><td>-</td><td>0.09</td><td>12.3</td><td>600</td><td>100</td><td>0</td><td>-</td></tr></tbody></table>

robustness. RCP-20 instances require less than a minute for Γ ≤ 2, on average, but using Γ equal to 3 and to 4, requires around 2 and 3 minutes, respectively. We also find that 34%, 46% and 67% of RCP-20 instances cannot be solved when Γ is equal to 2, 3 and 4, respectively. The results for RCP-30 instances reveal that all problems with Γ ≥ 2 are infeasible, while only 40% of these instances can be solved with Γ = 1 and all 100 RCP-30 instances are feasible in the deterministic case.

The impact of the maximum uncertainty is next examined. For this experiment, we compare the performance of the proposed formulation for a maximum level of robustness, i.e. Γ = 4 under varying maximum uncertainty $\bar{\epsilon}$: 0% (for the deterministic counterpart), 2.5%, 5%, 7.5% and 10%. The performance is reported similarly as in the previous section.

The results obtained using RCP instances (Table 3.6) reveal that while all problems are feasible under deterministic conditions, the proportion of infeasible problems increase with the maximum uncertainty and the number of aircraft. Notably, all RCP-30-1 instances are found to be infeasible for $\bar{\epsilon} \ge 2.5\%$ (recall that a maximal level of robustness is used). Examining the change in the objective value, for RCP-10 instances we observe that on average the total deviation of aircraft increases by three orders of magnitude. In terms of runtime, RCP-10 instances, are solved within 1 s using any value of $\bar{\epsilon}$. RCP-20 instances, require on average 0.24 s in the deterministic case. Using $\bar{\epsilon} = 5\%$, 93 instances are infeasible, 2 timed-out and among the remaining 5 instances solved to optimality the average runtime is 200 s.
<table><caption>Table 3.5: Summary of results for RCP instances with a speed control range of [−6%, +3%] and a heading control range of [−30°, +30°] using a maximal uncertainty of $\bar{\epsilon} = 5\%$.</caption><thead><tr><th>Γ</th><th>Instance</th><th>UB</th><th>Gap</th><th>Time</th><th>n<sub>i</sub></th><th>n<sub>t</sub></th><th>n<sub>φ</sub></th></tr></thead><tbody><tr><td rowspan="3">0</td><td>RCP-10</td><td>2.2e<sup>-4</sup></td><td>0.00</td><td>0.05</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-20</td><td>1.7e<sup>-3</sup></td><td>0.00</td><td>0.26</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-30</td><td>7.1e<sup>-3</sup></td><td>0.17</td><td>135</td><td>0.4</td><td>3</td><td>0</td></tr><tr><td rowspan="3">1</td><td>RCP-10</td><td>1.0e<sup>-2</sup></td><td>0.00</td><td>0.24</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-20</td><td>9.8e<sup>-2</sup></td><td>0.00</td><td>10.2</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-30</td><td>7.5e<sup>-1</sup></td><td>0.17</td><td>287</td><td>0.4</td><td>10</td><td>60</td></tr><tr><td rowspan="3">2</td><td>RCP-10</td><td>2.7e<sup>-2</sup></td><td>0.00</td><td>0.34</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-20</td><td>2.8e<sup>-2</sup></td><td>0.00</td><td>25.1</td><td>0.0</td><td>8</td><td>34</td></tr><tr><td>RCP-30</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>100</td></tr><tr><td rowspan="3">3</td><td>RCP-10</td><td>4.2e<sup>-2</sup></td><td>0.00</td><td>0.54</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-20</td><td>7.4e<sup>-1</sup></td><td>0.00</td><td>134</td><td>0.0</td><td>10</td><td>46</td></tr><tr><td>RCP-30</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>100</td></tr><tr><td rowspan="3">4</td><td>RCP-10</td><td>5.3e<sup>-2</sup></td><td>0.00</td><td>0.59</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-20</td><td>8.8e<sup>-1</sup></td><td>0.00</td><td>201</td><td>0.0</td><td>15</td><td>67</td></tr><tr><td>RCP-30</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>100</td></tr></tbody></table>

## 3.7 Conclusion

This chapter presented mixed-integer formulations and exact solution methods for aircraft conflict resolution problems (ACRP). We first considered the 2D ACRP with continuous speed and heading control manoeuvres and proposed compact disjunctive separation conditions. The proposed disjunctive separation conditions are linear with regards to aircraft relative velocity variables and only require a single binary variable per pair of aircraft. We introduced a simple pre-processing algorithm to identify aircraft pairs which are conflict-free or non-separable for any combination of controls, which may help in reducing the size of ACRPs by omitting conflict-free pairs. We built on and extended the complex number formulation for the ACRP introduced by Rey and Hijazi [2017] by augmenting its objective function with a preference weight to balance the trade-off between speed and heading deviations. The resulting formulation is a nonconvex mixed-integer program. This 2D formulation is extended to the context of altitude control by flight level (FL) change and we proposed a lexicographic optimization to solve the 2D+FL ACRP which aims to minimise the number of FL changes in priority and resolve outstanding conflicts by 2D trajectory control. The deterministic formulation was also extended to a stochastic variant wherein aircraft velocities are subject to uncertainty and a robust optimization approach was developed.

We presented exact solution algorithms to solve the nonconvex 2D and 2D+FL ACRPs. The proposed algorithms refine the convex relaxations introduced by Rey and Hijazi [2017]. For the 2D problem, the nonconvex formulations are first relaxed to mixed-integer quadratic programs (MIQP) which solution may violate the speed control constraint. If such violations occur, convex quadratic constraints are added together with a constraint generation algorithm that iteratively refines an outer piecewise linear approximation of the speed control constraint by solving a sequence of mixed-integer quadratically constrained programs (MIQCP). The 2D+FL lexicographic optimization problem is solved by decomposing the nonconvex problem into a flight assignment problem and a series of FL-based 2D problems. The proposed flight
<table><caption>Table 3.6: Summary of results for RCP instances with a speed control range of [−6%, +3%] and a heading control range of [−30°, +30°] using a level of robustness of Γ = 4.</caption><thead><tr><th>ε (%)</th><th>Instance</th><th>UB</th><th>Gap (%)</th><th>Time (s)</th><th>n<sub>i</sub></th><th>n<sub>t</sub></th><th>n<sub>θ</sub></th></tr></thead><tbody><tr><td rowspan="3">0.0</td><td>RCP-10</td><td>2.2e<sup>-4</sup></td><td>0.00</td><td>0.05</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-20</td><td>1.7e<sup>-3</sup></td><td>0.01</td><td>0.24</td><td>0.0</td><td>0</td><td>0</td></tr><tr><td>RCP-30</td><td>7.2e<sup>-3</sup></td><td>0.01</td><td>66.3</td><td>1.4</td><td>3</td><td>0</td></tr><tr><td rowspan="3">2.5</td><td>RCP-10</td><td>1.6e<sup>-2</sup></td><td>0.00</td><td>0.47</td><td>0.0</td><td>0</td><td>1</td></tr><tr><td>RCP-20</td><td>8.9e<sup>-1</sup></td><td>0.00</td><td>193</td><td>0.0</td><td>12</td><td>4</td></tr><tr><td>RCP-30</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>100</td></tr><tr><td rowspan="3">5.0</td><td>RCP-10</td><td>5.3e<sup>-2</sup></td><td>0.00</td><td>0.60</td><td>0.0</td><td>0</td><td>1</td></tr><tr><td>RCP-20</td><td>8.8e<sup>-1</sup></td><td>0.00</td><td>200</td><td>0.0</td><td>2</td><td>93</td></tr><tr><td>RCP-30</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>100</td></tr><tr><td rowspan="3">7.5</td><td>RCP-10</td><td>1.3e<sup>-1</sup></td><td>0.00</td><td>0.76</td><td>0.0</td><td>0</td><td>9</td></tr><tr><td>RCP-20</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>100</td></tr><tr><td>RCP-30</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>100</td></tr><tr><td rowspan="3">10.0</td><td>RCP-10</td><td>3.0e<sup>-1</sup></td><td>0.00</td><td>0.89</td><td>0.0</td><td>0</td><td>26</td></tr><tr><td>RCP-20</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>100</td></tr><tr><td>RCP-30</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>100</td></tr></tbody></table>

assignment formulation is based on a reformulation of the FL separation constraint which requires an exponential number of constraints and is embedded into an iterative approach to generate altitude separation constraints as needed.

The performance of the proposed mixed-integer formulations and algorithms was tested on a total of 2072 benchmarking instances. These instances include four types of ACRPs with up to 60 aircraft per instance for 2D problems and 150 aircraft per instance for 2D+FL problems. The performance of the proposed solution algorithms highlights the scalability of the approach compared to existing methods in the literature. Further, we find that the combination of the pre-processing algorithm with the MIQP convex relaxation is sufficient to solve FP and GP instances with up to 60 aircraft, and most RCP instances with up to 30 aircraft. We also find that the number of MIQCP iterations remains low on average when solving larger problems. For 2D+FL lexicographic optimization problems, we find that the pre-processing procedure generates enough altitude separation constraints to solve dense instances with an average of 50 aircraft per FL. Last, the results of the sensitivity analyses conducted highlight that while the robust ACRP can be solved without significantly increasing the level of computational resources required, the impact of on system costs (total deviation) increase rapidly with the level of robustness and/or uncertainty.
# Chapter 4

## Modeling approaches for autonomous intersection management

This chapter is based on Levin and Rey [2017], Rey and Levin [2019] and Rey et al. [2021].

### 4.1 Introduction

Traffic intersections are major bottlenecks of urban transport networks. Signalized traffic intersections are typically controlled so as to improve throughput, minimize vehicle delays or reduce emissions. The emergence of automated mobility in urban networks leads to new operational challenges that have received a growing attention over the past few years. Several new technologies have been proposed for AVs to improve traffic flow Chen and Englund [2016]. The overarching theme consists of replacing signalized traffic intersection with signal-free intersections that are controlled by an intersection manager. This has lead to the development of autonomous intersection management (AIM) wherein the control of signal-free traffic intersections is assumed to be in the hands of an autonomous agent governed by decision-making algorithms that aim to regulate traffic.

Of particular interest is the future paradigm wherein legacy (or human-operated) vehicles have been fully replaced with autonomous vehicles (AV). In this futuristic context, several studies have shown that the management of urban networks may benefit from forecasted technological advancements, such as improved trajectory control, automatic collision avoidance or platooning. However, the intermediate traffic state wherein legacy and autonomous vehicles co-exist has not been nearly as much examined by researchers. In such an intermediate traffic context, urban networks will need to adapt to make the most out of AV technology and allow the transition towards fully autonomous traffic, if this is ever to happen. The first AVs evolving in urban traffic networks are likely to have to abide by the existing infrastructure and legislation. However, this picture may change rapidly. We conjecture that with the increase in AV demand, urban traffic networks will adapt and that AV-specific infrastructure will be available to improve network operations at traffic intersections.

The emergence of automated mobility also offers the opportunity to rethink traffic in-
tersections in the light of users' preferences. Auction-based mechanisms have emerged as promising alternatives to traditional traffic intersection control approaches [Schepperle and Böhm, 2007, Vasirani and Ossowski, 2012, Carlino et al., 2013]. In auction-based mechanisms, users are assumed to be able to declare their preferences, e.g. value of time, to an AIM so as to obtain services commensurate to their need. The intersection can be viewed as a server which role is to process users' service requests. In this context, the intersection manager acts as a controller which decides users' service sequence and users' payments.

The chapter summarizes three approaches that address different facets of AIM. First, a mixed-integer programming formulation is presented to maximize throughput in a fully autonomous mobility context where only AVs share the infrastructure. Second, this formulation is embedded in a stochastic network traffic control with hybrid mobility wherein legacy vehicles (LV) and AVs co-exist. It is shown that in this hybrid mobility context, a max-pressure traffic control policy can stabilize the network. While the design criteria of these approaches have merits they are oblivious to users' preferences. In a last section, a mechanism design perspective is considered wherein users have the possibility to bid for priority service at traffic intersections.

## 4.2 Literature review

We first discuss the literature on autonomous intersection management before discussing studies on network traffic control. We then review efforts that have considered a traffic intersection auction context.

### 4.2.1 Autonomous intersection management

Due to the real-time computing needs of auction-based traffic intersection mechanisms, such approaches are typically conceived in the context of autonomous intersection management (AIM), as proposed in the seminal work of Dresner and Stone [2004, 2008]. In this paradigm, traffic control is assumed to be signal-free and users are assumed to be able reserve space-time trajectories through intersections. Several works have built on and extended the AIM protocol to richer configurations. Fajardo et al. [2011] and Li et al. [2013] proposed AIM protocols based on a First-Come-First-Served (FCFS) policies where vehicles are prioritized based on their arrival time at the intersection. De La Fortelle and Qian [2015] and Altché and De La Fortelle [2016] developed microscopic vehicle trajectory optimization formulations to coordinate vehicles through intersections. In these formulations, the intersection manager decides vehicles' service time and speed while time is discretized. Zhang et al. [2016, 2017] and Malikopoulos et al. [2018] proposed decentralized approaches for traffic control based on FCFS conditions and considered energy consumption as well as vehicle separation and throughput in their formulations. Mirheli et al. [2019] developed a consensus-based approach for cooperative trajectory planning at signal-free intersections and shows that near-optimal solution can be achieved in competitive time. Wu et al. [2019] designed a multi-agent Markov decision process for cooperative trajectory planning in AIM. Levin and Rey [2017] proposed a conflict point model that obviates the need to discretize the intersection space using tiles. This conflict point model was then adapted by Rey and Levin [2019] to accommodate both legacy and automated vehicles with signalized and signal-free traffic phases; and more recently by Chen et al. [2020] to account for pedestrians movements.
### 4.2.2 Network traffic control

The literature on network traffic control spans across several fields such as optimization and control theory, transportation engineering and computer science. The seminal work of Tassiulas and Ephremides [1992] pioneered the research on stability conditions in network traffic control with an application to data packets routing in communication networks. The authors notably introduced the concept of back-pressure algorithm as method for decentralized control of a network of routers. In the context of urban transport networks, there is an extensive body of research on network traffic signal control. Varaiya [2013] proposed a network traffic control policy based on max-pressure, a variant of the back-pressure algorithm wherein the control policy chooses the signal phase that maximizes the pressure at each intersection. Stability was proven assuming each turning movement has a distinct queue. Wongpiromsarn et al. [2012] proposed a pressure-based policy which maximizes throughput for unbounded queues. However, practical limitations such as link length require a careful choice of the pressure function to avoid queue spillback. Building on this effort, Xiao et al. [2014] proposed a pressure-releasing policy that accounts for finite queue capacities. Nonetheless, to more canonically apply the pressure-based routing they assumed that each turning movement has a separate queue, which is often not realistic. Le et al. [2015] proposed a fixed-time policy wherein all phases are assigned a non-zero activation time and proved stability for fixed turning proportions and unbounded queues. Recently, Valls et al. [2016] propose a convex optimization approach to traffic signal control that decouples the stability of the system from the choice of traffic control policy.

Starting from the seminal work of Smith [1979], several efforts have also attempted to model the impact of traffic signal optimization on route choice [Gregoire et al., 2014, Zaidi et al., 2016]. Le et al. [2017] used utility functions in max-pressure control to influence routing. However, while Tassiulas and Ephremides [1992]'s policy is provably throughput-optimal, max-pressure route choice makes no guarantees on the efficiency of the travel times. We assume fixed route choice in our network traffic control model to focus on green and blue phases signal control, and leave route choice for later studies.

### 4.2.3 Traffic intersection auctions

Traffic intersection auctions have received a growing attention over the past few years, notably with the advent of connected and automated vehicle technology. Schepperle and Böhm [2007] proposed a subsidy-based mechanism for slot allocation which aim to balance vehicle waiting time. First-In First-Out (FIFO) constraints are accounted for but incentive-compatibility is not guaranteed. The “effect of starvation” and its impact on fairness are discussed. Vasirani and Ossowski [2012] developed a policy based on combinatorial auctions for the allocation of reservations at traffic intersections. The auction winner pays a price that is exactly the bid that was submitted, which is not incentive-compatible. Carlino et al. [2013] proposed several types of auctions where users can bid on phases or reservations (one vehicle at a time). All drivers can participate in the auction but the candidates only include drivers at the front of their lane. Winners split the cost of the second-highest bid with proportional payment, which yields a static incentive-compatible mechanism and the strategyproofness of the dynamic case is not guaranteed. Lloret-Battle and Jayakrishnan [2016] proposed envy-minimizing mechanisms for traffic intersections with traditional phasing. The authors developed a revenue-neutral, pareto-efficient approach that minimizes envy and delays by considering monetary transfers across users. A formal proof of incentive-compatibility is not provided although sim-
<table><caption>Table 4.1: Summary of the state-of-the-art on traffic intersection auctions. The column **Mechanism** describes the mechanism; **Auction** indicates the type of auction; **IC** indicates whether the mechanism is incentive-compatible; **F** indicates if fairness is considered in the mechanism; **E** indicates if the mechanism is efficient, i.e. social welfare is optimized; **C** indicates if intersection capacity is considered or optimized.</caption><thead><tr><th>Paper</th><th>Mechanism</th><th>Auction</th><th>IC</th><th>F</th><th>E</th><th>C</th></tr></thead><tbody><tr><td>Schepperle and Böhm [2007]</td><td>subsidy-based</td><td>static</td><td>no</td><td>some</td><td>no</td><td>some</td></tr><tr><td>Vasirani and Ossowski [2012]</td><td>combinatorial auction</td><td>static</td><td>no</td><td>no</td><td>yes</td><td>yes</td></tr><tr><td>Carlino et al. [2013]</td><td>marginal cost</td><td>static</td><td>yes</td><td>yes</td><td>yes</td><td>some</td></tr><tr><td>Lloret-Battle and Jayakrishnan [2016]</td><td>envy-minimizing</td><td>static</td><td>no</td><td>yes</td><td>yes</td><td>yes</td></tr><tr><td>Sayin et al. [2018]</td><td>marginal cost</td><td>static</td><td>yes</td><td>no</td><td>yes</td><td>yes</td></tr><tr><td>Censi et al. [2019]</td><td>karma system</td><td>static</td><td>no</td><td>yes</td><td>yes</td><td>some</td></tr><tr><td>Lin and Jabari [2020]</td><td>transferable utility</td><td>static</td><td>no</td><td>some</td><td>yes</td><td>yes</td></tr><tr><td>Rey et al. [2021]</td><td>expected marginal cost</td><td>dynamic</td><td>yes</td><td>no</td><td>yes</td><td>some</td></tr></tbody></table>

ulations show that arbitrage is unlikely. Sayin et al. [2018] proposed an auction protocol and mechanism in which all vehicles within range of the intersection manager communicate their value of time, and the intersection manager assigns reservations to those vehicles via a static auction. The proposed mechanism is incentive-compatible in the static sense, i.e. assuming all participants in the auction are known at the time payments are determined. Censi et al. [2019] introduced a credit-based auction mechanism, i.e. a karma system in which agents pay other agents karma for priority. Agents with a low priority today have an incentive to lose bids so they can achieve higher priority tomorrow by acquiring more karma. The authors attempt to calculate a Nash equilibria of user bids, and they show that a centralized strategy can be more unfair than some of the Nash equilibria identified. Recently, Lin and Jabari [2020] proposed a mechanism for pricing intersection priority based on transferable utility games. In each game, players have the possibility to trade time among themselves and “winners pay losers to gain priority”. The authors provide empirical evidence that their approach is robust against adversarial user behavior, but no formal proof of incentive-compatibility is presented.

A synthesis of the state-of-the-art on traffic intersection auctions is presented in Table 4.1. The column **Mechanism** describes the mechanism of the corresponding paper. The column **Auction** indicates the type of auction: static or dynamic. A static auction is an auction in which the set of participants is assumed known, whereas a dynamic auction allows participants to enter and leave the auction dynamically. Observe that, in a static traffic intersection auction, losers may need participate in multiple auction rounds until being serviced during which new users may arrive and further delay these losers. The column **IC** indicates whether the mechanism is incentive-compatible or not. The column **F** indicates if fairness is considered in the mechanism, either to some extent (some) or via constraints within the mechanism (yes). The column **E** indicates if the mechanism is efficient, i.e. if social welfare is optimized. The column **C** indicates if intersection capacity is considered in the mechanism, either to some extent (some) or explicitly optimized (yes). Table 4.1 highlights that, to the best of our knowledge, existing mechanisms for traffic intersection auctions are restricted to static auctions. This synthesis also emphasizes that fairness considerations and intersection capacity are seldom optimized jointly, thus underlining existing research gaps in the design of traffic intersection auction mechanisms.
<figure><img src="image_10.png" /><figcaption>Figure 4.1: 4-approach, 2-lane intersection. Four vehicle paths are available $WS_R$, $WE_R$, $WE_L$, $WN_L$ for the West approach. All conflict points on each trajectory are represented with a rhombus.</figcaption></figure>

## 4.3 Conflict point formulation

We now develop an MILP to determine the optimal reservations for a single intersection. We replace the tiles in early work on reservations Dresner and Stone [2004] with conflict points, which are the intersections between different turning movement paths. Conflict points are an alternative method of ensuring collision avoidance that are more amenable to mathematical programming formulations. Figure 4.1 illustrates the conflict points for a standard four-approach intersection. Details on the computation of geometrical parameters are provided in Levin and Rey [2017].

Two vehicles are separated if they do not occupy the same conflict point at the same time. This is similar to the separation requirement for tiles [Dresner and Stone, 2004, 2006], but conflict points are only placed where potential conflicts occur. Conflict points have previously been used to model reservations by Zhu and Ukkusuri [2015], although their formulation was intended for use within dynamic traffic assignment rather than as a deployable solution to reservation-based intersection control. Conflict points have also been used for aircraft separation models Rey et al. [2012, 2015, 2016c], but aircraft require different separation constraints.

### 4.3.1 Notation

Consider a single intersection with set of vehicles $\mathcal{V}$. The intersection has sets of entry and exit lanes, $\Gamma^-$ and $\Gamma^+$, respectively. Let $\mathcal{C}$ be the set of all conflict points. Each vehicle requests a turning movement traveling from an entry lane to an exit lane. We assume all specified turning movements are legal and physically feasible. Each vehicle requires a path $\rho_i$ through the intersection, which is an ordered set of conflict points. For each vehicle $i$, we
must determine $t_i(c)$, which is the time that $i$ passes through point $c$, for all $c \in \rho_i$. The first point in $\rho_i$ is the entrance to the intersection, and the last point is the exit point for the intersection. Then $t_i(\gamma_i^-)$ is the time at which $i$ enters the intersection.

We assume that vehicles travel at a uniform speed through the intersection, but that the speed is a decision variable. The reason for this assumption is that varying speeds can greatly complicate determining how long vehicles occupy conflict points because vehicles may occupy multiple conflict points simultaneously. Although the uniform speed assumption may be relaxed in future work, we use it here to focus on the exposition of the conflict avoidance formulation.

Let $\underline{U}_i$ and $\overline{U}_i$ be the minimum and maximum speed for which $i$ can travel through the intersection, respectively. We require a minimum speed to ensure that vehicles to do not spend arbitrarily long times in the intersection. For each vehicle $i$, and conflict point $c \in \rho_i$, $i$ takes some time $\tau_i(c)$ to completely pass through $c$. $\tau_i(c)$ is a function of the vehicle speed during the turning movement. Thus, when $i$ moves through the intersection, $c$ is unusable during the time interval $[t_i(c), t_i(c) + \tau_i(c))$. An appropriate buffer may be added to $\tau_i(c)$ to ensure sufficient separation. Observe that $t_i(\gamma_i^+) + \tau_i(\gamma_i^+)$ is the time at which $i$ exits the intersection.

### 4.3.2 Vehicle ordering

Assume that vehicles do not change lanes while under the control of the intersection manager. Therefore, FIFO holds on all lanes. Let $e_i$ denote the earliest time that $i$ could reach the intersection. Then

$$t_i(\gamma_i^-) \geq e_i, \quad \forall i \in \mathcal{V}. \qquad (4.1)$$

Assume without loss of generality that $e_i < e_j$, so $i$ comes before $j$ in the lane queue. Then $t_i(\gamma_i^-) < t_j(\gamma_j^-)$ must hold as well. Because $i$ occupies $\gamma_i^-$ for $\tau_i(\gamma_i^-)$:

$$t_i(\gamma_i^-) + \tau_i(\gamma_i^-) \leq t_j(\gamma_j^-), \quad \forall i, j \in \mathcal{V}: \gamma_i^- = \gamma_j^-, e_i < e_j. \qquad (4.2)$$

### 4.3.3 Travel time between conflict points

For some vehicle $i$ and any two conflict points $c_1, c_2 \in \rho_i$ with $c_2$ after $c_1$ in $\rho_i$, let $d_i(c_1, c_2)$ be the distance $i$ travels between $c_1$ and $c_2$. Note that $d_i(c_1, c_2)$ may not be the Euclidean distance between $c_1$ and $c_2$ because turning vehicles may travel along an elliptical path (e.g. for left turns). The minimum and maximum vehicle speeds place constraints on the travel time through the intersection:

$$\frac{d_i(\gamma_i^-, \gamma_i^+)}{\overline{U}_i} \leq t_i(\gamma_i^+) - t_i(\gamma_i^-) \leq \frac{d_i(\gamma_i^-, \gamma_i^+)}{\underline{U}_i}, \quad \forall i \in \mathcal{V}. \qquad (4.3)$$

Furthermore, we require that vehicles maintain the same speed throughout the intersection. This is to ensure that vehicles move through conflict points properly. If several conflict points are close together, with distance smaller than vehicle length, speed changes could result in the vehicle blocking a conflict point but not triggering the corresponding constraint.
Therefore:

$$ \frac{t_i(c) - t_i(\gamma_i^-)}{d_i(\gamma_i^-, c)} = \frac{t_i(\gamma_i^+) - t_i(\gamma_i^-)}{d(\gamma_i^-, \gamma_i^+)} \quad \forall i \in \mathcal{V}, \forall c \in \rho_i. \quad (4.4) $$

### 4.3.4 Time spent at conflict points

For each vehicle $i$ traversing the intersection, we can determine the time spent at each conflict point $c \in \rho_i$, $\tau_i(c)$. based on geometry of $i$ and turning movement shape (straight or circular). Specifically, let $D_i(c)$ be the spacing required by $i$ at $c$. This spacing should depend on the vehicle speed in case an emergency braking is required to avoid a collision for unforeseen circumstances (such as a vehicle stall or a pedestrian jaywalking). The spacing for a given speed may be derived from the speed-density relationship, and can be used to specify the desired flow-density relationship for reservation travel. Equation (4.4) gives the speed of vehicle $i$, $u_i$, as:

$$ u_i = \frac{d_i(\gamma_i^-, \gamma_i^+)}{t_i(\gamma_i^+) - t_i(\gamma_i^-)}. \quad (4.5) $$

Using the inverse speed of the vehicle, the time vehicle $i$ spends at conflict point $c$ is:

$$ \tau_i(c) = \frac{D_i(c)}{u_i}, \quad \forall i \in \mathcal{V}, \forall c \in \rho_i. \quad (4.6) $$

For simplicity in the formulation, assume that the fundamental diagram is triangular Newell [1993], Yperman et al. [2005]. The fundamental diagram is defined as:

$$ q(k) = \min \{ \overline{U}_i k, w(K-k) \}, \quad (4.7) $$

where $k$ is the density, $K$ is the jam density, $w$ is the congested wave speed, and $q(k)$ is the flow-density relationship. We assume that $K$ and $w$ are constants.

Speeds are constrained by spacing only in the congested region, and from equation (4.7) the speed satisfies:

$$ u_i = \frac{w}{k} (K - k). \quad (4.8) $$

Let $L_i(c)$ be the distance the vehicle travels around $c$. If the turning movement is a straight line, $L_i(c) = l_i$, where $l_i$ is the length of $i$. For left and right turns, the vehicle will travel for a greater distance while occupying $c$. For circular turns (e.g. left and right turns in a symmetric intersection), $L_i(c) = 4R \sin^{-1}(l_i / (2R))$ where $R$ is the turn radius.

Using the fact that spacing is the inverse of density, i.e. $D_i(c) = \frac{1}{k}$ and $K = \frac{1}{L_i(c)}$, equation (4.8) can be rewritten as

$$ D_i(c) = L_i(c) \frac{u_i + w}{w}. \quad (4.9) $$
From equation (4.9), equation (4.6) can then be rewritten as

$$
\begin{align}
\tau_i(c) &= L_i(c) \frac{u_i + w}{u_i w}, \nonumber \\
&= \frac{L_i(c)}{w} + \frac{L_i(c)}{u_i}, \nonumber \\
&= \frac{L_i(c)}{w} + \frac{L_i(c) (t_i(\gamma_i^+) - t_i(\gamma_i^-))}{d_i (\gamma_i^-, \gamma_i^+)}, \tag{4.10}
\end{align}
$$

which is linear in the decision variables of conflict point arrival times.

4.3.5 Separation at conflict points

For any two vehicles i, j ∈ V, let δ_{ij}(c) ∈ {0, 1} denote whether i and j overlap at c. δ_{ij}(c) is sensitive to the order of vehicles, so each pair of vehicles has two such δ variables per shared conflict point. δ_{ij}(c) = 1 if j arrives at c after i has left. Thus, δ_{ij}(c) = 1 if t_j(c) ≥ t_i(c) + τ_i(c). This is written as

$$
t_i(c) + \tau_i(c) - t_j(c) \le (1 - \delta_{ij}(c)) M_{ij}, \quad \forall i, j \in V, \forall c \in \rho_i \cap \rho_j, \quad (4.11)
$$

where $M_{ij}$ is a large positive constant. If $t_j(c) \ge t_i(c) + \tau_i(c)$ then $\delta_{ij}(c) = 1$; otherwise $\delta_{ij}(c) = 0$. Note that it suffices to have $M_{ij} \ge \max\{t_i(c) + \tau_i(c) - t_j(c)\}$ for this separation constraint to be valid.

Separation is guaranteed for any two pair of vehicles i, j ∈ V and any conflict point c if j enters c after i leaves, or if i enters c after j leaves:

$$
\delta_{ij}(c) + \delta_{ji}(c) = 1, \quad \forall i, j \in \mathcal{V} : i < j, \forall c \in \rho_i \cap \rho_j. \quad (4.12)
$$

Constraint (4.12) does not depend on vehicle ordering, and therefore occurs once per pair of vehicles per shared conflict point.

Note that the number of binary $\delta_{ij}(c)$ variables may be reduced for vehicles with $\gamma_i^- = \gamma_j^-$. For such pairs vehicles, the ordering at the entry point determines the ordering at any conflict points that they share. Constraints (4.11) and (4.12) may therefore be replaced with

$$
t_i(c) + \tau_i(c) \le t_j(c), \quad \forall i, j \in \mathcal{V} : \gamma_i^- = \gamma_j^-, e_i < e_j, \forall c \in \rho_i \cap \rho_j, \quad (4.13)
$$

for pairs of vehicles i, j ∈ V with $\gamma_i^- = \gamma_j^-$.

4.3.6 Objective function

The objective function represents how the intersection manager aims to operate traffic. An intuitive objective consists in maximizing intersection throughput. While this approach may maximize social welfare, such an intersection control policy may also lead to imbalanced vehicle delays, thus potentially impacting acceptability. Alternatively, the intersection man- ager may seek to minimize the maximum delay experienced or maximize fairness. These alternative objective functions aim to incorporate some level of fairness in the allocation of reservation assignments at the vehicular scale. More elaborate fair objective functions can
also be designed to balance reservation assignments across routes or traffic approaches (e.g. North, South, East, West).

Because the constraints guarantee complete separation of vehicles through the intersection, the objective function decides in which order vehicles should move. Denote by $f_i > 0$ the weight of each vehicle, which indicates its priority for movement. The weights can model a variety of situations. For instance, $f_i = 1$ for all vehicles is maximum throughput. $f_i = e_i$ weights vehicles by earliest arrival time, which is similar to a fairness based objective. Note that $f_i = e_i$ is not the same as the FCFS policy used in previous work on reservations, because FCFS results in absolute priority for the vehicle with the earliest arrival time. The objective $f_i = e_i$ may delay early-arriving vehicles to find the combination with the greatest fairness.

The proposed objective function values movements by the time vehicles exit the intersection, $t_i(\gamma_i^+) + \tau_i(\gamma_i^+)$:

$$ \min \sum_{i \in \mathcal{V}} f_i \left( t_i(\gamma_i^+) + \tau_i(\gamma_i^+) \right). \qquad (4.14) $$

Exit times are chosen because vehicle speeds are an implicit decision variable. Both vehicle delays and slower turn speeds will increase the objective value. Note that the objective function (4.14) will never be equal to 0 if $f_i > 0$. To obtain a measure of delay for vehicle $i$, we can subtract the constant term $\frac{d_i(\gamma_i^-, \gamma_i^+)}{\bar{U}_i} + e_i$, which is the earliest time that vehicle $i$ could exit the intersection.

### 4.3.7 Mixed integer linear program

Combining the above constraints and objective function results in the MILP represented by (4.15) through (4.24), henceforth referred to as the *Conflict-Point Intersection Control*
(CPIC) formulation.

$$ \min \sum_{i \in \mathcal{V}} f_i \left( t_i \left( \gamma_i^+ \right) + \tau_i \left( \gamma_i^+ \right) \right), \qquad (4.15) $$

$$ \text{s.t.} \quad t_i(\gamma_i^-) \ge e_i, \qquad \forall i \in \mathcal{V}, \qquad (4.16) $$

$$ t_i(\gamma_i^-) + \tau_i(\gamma_i^-) \le t_j(\gamma_j^-), \qquad \forall i, j \in \mathcal{V}: \gamma_i^- = \gamma_j^-, e_i < e_j, \quad (4.17) $$

$$ \tau_i(c) = \frac{L_i}{w} + \frac{L_i(t_i(\gamma_i^+) - t_i(\gamma_i^-))}{d_i(\gamma_i^-, \gamma_i^+)}, \qquad \forall i \in \mathcal{V}, \forall c \in \rho_i, \qquad (4.18) $$

$$ \frac{d_i(\gamma_i^-, \gamma_i^+)}{\bar{U}_i} \le t_i(\gamma_i^+) - t_i(\gamma_i^-) \le \frac{d_i(\gamma_i^-, \gamma_i^+)}{\underline{U}_i}, \qquad \forall i \in \mathcal{V}, \qquad (4.19) $$

$$ \frac{t_i(c) - t_i(\gamma_i^-)}{d_i(\gamma_i^-, c)} = \frac{t_i(\gamma_i^+) - t_i(\gamma_i^-)}{d(\gamma_i^-, \gamma_i^+)}, \qquad \forall i \in \mathcal{V}, \forall c \in \rho_i, \qquad (4.20) $$

$$ t_i(c) + \tau_i(c) \le t_j(c), \qquad \forall i, j \in \mathcal{V}: \gamma_i^- = \gamma_j^-, e_i < e_j, \forall c \in \rho_i \cap \rho_j, \qquad (4.21) $$

$$ t_i(c) + \tau_i(c) - t_j(c) \le (1 - \delta_{ij}(c))M_{ij}, \qquad \forall i, j \in \mathcal{V}: \gamma_i^- \ne \gamma_j^-, \forall c \in \rho_i \cap \rho_j, \qquad (4.22) $$

$$ \delta_{ij}(c) + \delta_{ji}(c) = 1, \qquad \forall i, j \in \mathcal{V}: \gamma_i^- \neq \gamma_j^-, i < j, \forall c \in \rho_i \cap \rho_j, \qquad (4.23) $$

$$ \delta_{ij}(c) \in \{0, 1\}, \qquad \forall i, j \in \mathcal{V}: \gamma_i^- \neq \gamma_j^-, \forall c \in \rho_i \cap \rho_j, \qquad (4.24) $$

The decision variables are the arrival times at each conflict point for each vehicle $i$, $t_i(c)$, and vehicle ordering at conflict points, determined by $\delta_{ij}(c)$ variables. The times spent at conflict points, $\tau_i(c)$, is defined by the speed at which vehicle $i$ moves through the intersection, which is implicitly defined by vehicle-specific arrival times at conflict points. It is not known whether Model CPIC is NP-hard. However, similar collision avoidance models for aircraft [Chiang et al., 1997] can be generalized as a motion planning problem in the presence of obstacles, which is NP-hard in most cases [Reif and Sharir, 1994].

Model CPIC is always feasible because there is no upper bound on vehicle entry and exit times. If necessary, vehicles could be given arbitrarily long delays to ensure that collision avoidance constraints are satisfied, this is summarized in the following proposition.

**Proposition 4.1.** *Model CPIC has at least one feasible solution.*

Vehicles' waiting time at the intersection entry points is unbounded, hence the values $M_{ij}$ may also be unbounded. Specifically, if $\delta_{ij} = 0$, then Constraint (4.22) becomes $t_i(c) + \tau_i(c) - t_j(c) \le M_{ij}$, which should hold for any values of the decision variables involved therein. Since $\max\{t_i(c) + \tau_i(c) - t_j(c)\} \to \infty$, the big-M values should be chosen carefully to maintain feasibility and avoid cutting out the optimal solution. The proof of Proposition 4.1 provides a method to construct a feasible solution to Model CPIC. We can use this solution to tighten the big-M values $M_{ij}$. Let $i = |\mathcal{V}| - 1$ denote the last vehicle in the sorted set $\mathcal{V}$ and let $\bar{t} = t_{|\mathcal{V}|-1}(\gamma_{|\mathcal{V}|-1}^+) + \tau_{|\mathcal{V}|-1}(\gamma_{|\mathcal{V}|-1}^+)$ denote the latest instant reserved by the intersection. Then $M = \bar{t} - e_0$ is a valid big-M value for any pair of vehicles.
## 4.4 Hybrid network traffic control policy

We now embed the intersection-based conflict point formulation into a stochastic network model wherein LVs and AVs co-exist.

### 4.4.1 Stochastic network model

Consider now a traffic network $\mathcal{G} = (\mathcal{N}, \mathcal{A})$ with a set of intersections $\mathcal{N}$ connected by a set of links $\mathcal{A}$. The set of links is partitioned into three subsets: internal links that connect two intersections, denoted $\mathcal{A}_o \subset \mathcal{A}$, source links at which vehicles enter the network, denoted $\mathcal{A}_r \subset \mathcal{A}$, and sink links at which vehicles exit, denoted $\mathcal{A}_s \subset \mathcal{A}$. We model congestion at network intersections using point-queues of infinite size and we are interested in the evolution of queue lengths over the entire network.

We consider two classes of vehicles and lanes: autonomous vehicles (AVs), denoted $a$, and legacy (or human-driven) vehicles (LVs), denoted $l$. We assume that each link of the network consists of a set of lanes which are either restricted to AVs (AV-lanes) or available to both AVs and LVs (LV-lanes). We use $\mathcal{A}_a$ and $\mathcal{A}_l$ to denote AV-lanes and LV-lanes, respectively, and we assume that vehicle movement between different classes of lanes are forbidden. Although AVs can use LV-lanes, we do not model any type of interaction at the traffic-flow level among AVs and LVs on LV-lanes: if an AV uses an LV-lane we assume that it behaves as an LV. We assume that each class of lanes is served by a color-coded traffic phase. Specifically, we assume that LV-lanes are served exclusively by traditional green signal phases. In turn, we assume that AV-lanes are served exclusively by signal-free blue phases. Blue phases differ from green phases in that they directly control vehicles' trajectory within the intersection. The proposed hybrid network control policy presented hereafter chooses which phase (green or blue) should be activated at each intersection of the network and each time period over a discretized time horizon.

Let $x_i(t) \in \mathbb{R}_+$ be the number of vehicles on link $i \in \mathcal{A}$ seeking to enter the intersection at time $t$. Although we discretize vehicles in our numerical experiments, integer queue lengths are not necessary for the analytical results presented hereafter. Let $x(t)$ be the array of all queue lengths at time $t$. $x(t)$ is the state of the network, and the state space is $\mathcal{X} = \{x_i(t) \ge 0 : i \in \mathcal{A}\}$. We consider discretized time and we assume fixed phase time of length $\Delta t$. Further, we assume that the state of the network $x(t)$ is known at each time period $t = 0, \Delta t, 2\Delta t, \dots$. The goal is to design a throughput-optimal network traffic control policy that optimally selects a traffic signal control at each time period $[t, t + \Delta t]$.

The proposed stochastic network traffic control formulation is based on the concept of vehicle movements which are formally defined below.

**Definition 4.1.** A movement $(i, j) \in \mathcal{A}^2$ is a vehicle trajectory from lane $i$ to lane $j$ across a common intersection $n \in \mathcal{N}$ in the network. We denote $\mathcal{M}$ the set of all movements in the network.

AVs communicate their position with IMs to make use of the AIM protocol, hence we assume that movement-specific queues are known for AVs. In contrast, LV-queues can be detected through loop detectors and flow sensors currently in use for traffic signals but their destination is assumed unknown. Specifically, let $\mathcal{A}_a \subset \mathcal{A}$ be the set of AV-restricted lanes. For these lanes, we assume known movements queues, i.e. if $j \in \mathcal{A}_a$, then $x_i(t) = \sum_{j \in \mathcal{A}_a} x_{ij}(t)$
and $x_{ij}(t)$ is known. In contrast, for other lanes $i \in \mathcal{A}_l = \mathcal{A} \setminus \mathcal{A}_a$, only $x_i(t)$ is known since route choice for LVs is assumed unknown.

For each lane $i \in \mathcal{A}$ we assume that lane capacity $C_i$ is known and determined assuming a triangular fundamental diagram relationship, that is $C_i = \frac{\bar{U}_i w K}{\bar{U} + w}$, where $\bar{U}_i$ is the free-flow speed on lane $i$, $K$ is the jam density and $w$ is the congestion wave speed. We also assume that the maximum, unconditional movement service rate $\bar{s}_{ij}$ is known for each movement and determined based on lost time, specifically: unconditional movement service rates for all movements $(i, j)$ are calculated as $\bar{s}_{ij} = \min\{C_i, C_j\} \frac{\Delta t - L}{\Delta t}$.

Let $P_{ij}(t)$ be a random variable denoting the turning proportion from lane $i$ to $j$ at time $t$ with known mean $p_{ij}$. Let $D_i(t)$ be a random variable denoting the external incoming traffic onto lane $i$ at time $t$ with known mean $d_i$. We denote $\mathbf{p}$, $\mathbf{d}$ and $\bar{\mathbf{s}}$ the vectors of mean turning proportions, mean demands and unconditional movement service rates, respectively. These upper bounds on movement service rates represent the maximum number of vehicles that can be moved from lane $i$ to lane $j$ during time period $t$ when no conflicting movement with $(i, j)$ is activated. From these rates, we can determine maximum, unconditional lane service rates $\bar{s}_i = \sum_{j \in \mathcal{A}:(i,j) \in \mathcal{M}} \bar{s}_{ij}$. For convenience if $i$ and $j$ do not correspond to a possible movement in the network, e.g. they belong to different intersections or they are both entry or exit lanes of the same intersection, we assume that $p_{ij} = 0$. Hence, we can define the sets of AV-movements $\mathcal{M}_a \equiv \{(i,j) \in \mathcal{A}_a^2 : p_{ij} \neq 0\}$ and LV-movements $\mathcal{M}_l \equiv \{(i,j) \in \mathcal{A}_l^2 : p_{ij} \neq 0\}$.

Traffic at network intersections is coordinated by phases which are determined by the selection of the activation matrix. In addition, we also introduce the concept of service matrix which is used in the proposed traffic control formulations.

**Definition 4.2.** An activation matrix $\beta(t)$ is a $|\mathcal{A}| \times |\mathcal{A}|$-matrix wherein all entries take value 0 (inactive) or 1 (active), i.e. $\beta_{ij}(t) \in \{0,1\}$ at time $t$.

**Definition 4.3.** A service matrix $\alpha(t)$ is a $|\mathcal{A}| \times |\mathcal{A}|$-matrix wherein all entries take a value between 0 (not serviced) and 1 (fully serviced), i.e. $\alpha_{ij}(t) \in [0, 1]$ at time $t$.

The entries of an activation matrix characterize the *activeness* of the corresponding phase: $\beta_{ij}(t) = 1$ means that movement $(i,j)$ is active during phase $t$ whereas $\beta_{ij}(t) = 0$ means that movement $(i,j)$ is inactive. The entries of the service matrix characterize the *service level* of active movements during phase $t$: $\alpha_{ij}(t) = 1$ corresponds to a maximal service level for of movement $(i,j)$, whereas $\alpha_{ij}(t) = 0$ means that movement $(i,j)$ cannot serve any vehicles during time period $t$. Fractional service level values model situations where conflicting movements, i.e. posing a safety risk, simultaneously have non-zero activation values. The activation and service matrices are linked through the movement-based constraints $\alpha_{ij}(t) \le \beta_{ij}(t)$. These linking constraints ensure that an intersection can only service vehicles on movement $(i,j)$ if this movement is active $\beta_{ij}(t) = 1$. Further, for priority movements, we set $\alpha_{ij}(t) = \beta_{ij}(t)$ which implies that an activated priority movements has full service level. In the proposed traffic control policy, the selection of the activation matrix requires the solution of two mathematical optimization problems. Green phases serve LVs using a MILP that captures FIFO blocking effects on lanes. Blue phases serve AVs using an modified version of the MILP presented in Section 4.3 wherein the objective function is adapted for network traffic control. These MILPs are referred to as GREEN and BLUE. Further details on these MILPs can be found in Rey and Levin [2019].

Let $Y_i(t)$ be a random variable with mean $y_i(t)$ denoting the number of vehicles serviced in lane $i$ at time $t$. The vector $\mathbf{Y}(t)$ is endogenous to the service matrix $\alpha(t)$ selected by the
control policy. Note that the mean $y_i(t)$ of $Y_i(t)$ is unknown and $y_i(t)$ is modeled as a control variable in the proposed traffic phase optimization formulation. The proposed stochastic network traffic control model is summarized by the lane-queue evolution equation (4.25):

$$x_j(t + \Delta t) = x_j(t) - Y_j(t) + \sum_{i \in A} P_{ij}(t)Y_i(t) + D_j(t), \quad \forall j \in A. \quad (4.25)$$

Note, if $j \notin A_r$ then $D_j(t) = 0$. Conversely, if $j \in A_r$ then $j$ has no predecessor links, thus $\sum_{i \in A} P_{ij}(t) = 0$. Although AV-lanes and other lanes have identical queue evolution equations, the information available is more accurate for AV-lanes, this used in the calculation of network control policies.

In related works, the lane service rate vector $y(t)$ is commonly calculated as the minimum between the supply (lane or movement service capacity) and the demand (lane or movement queue length) [Varaiya, 2013, Le et al., 2015]. However, this modeling approach does not always capture the interdependency between the activation of possibly conflicting movements with lane or movement service capacity and queue length. Precisely, previous efforts have assumed that a set of activation matrices is provided for each intersection and that lane or movement service capacities can be pre-processed accordingly. This overlooks the impact of queue length on intersection capacity. In contrast our proposed integrated approach aims to accurately estimate the expected number of serviced vehicles at each time period by leveraging the available lane-queue length information (note that this information is also assumed available in the aforementioned papers).

### 4.4.2 Stability analysis

We present a network traffic control policy for intersection control combining green (LV-lane restricted) and blue phases (AV-lane restricted) and prove that it maximizes throughput. The proposed network traffic control policy works by repeatedly solving GREEN and BLUE at each time period $t$ based on the network state $x(t)$ and combining local (intersection-level), optimal activation matrices into a network-wide activation matrix $\alpha(t)$.

Let $\mathbb{M}$ be the set of service matrices. A policy is a function $\pi: \mathcal{X} \to \mathbb{M}$ that chooses a service matrix $\alpha(t) \in \mathbb{M}$ for every state $x(t)$. We use the concept of strong stability to characterize the proposed stochastic queuing process [Leonardi et al., 2001, Wongpiromsarn et al., 2012, Zaidi et al., 2016]:

**Definition 4.4.** Let $\bar{t}$ be a time period index. A stochastic queue evolution process is strongly stable under policy $\pi$ if and only if there exists $K < \infty$ such that the network state $x(t)$ verifies:

$$\limsup_{\bar{t} \to \infty} \mathbb{E} \left[ \frac{1}{\bar{t}} \sum_{t=1}^{\bar{t}} |x(t)| \right] < K. \quad (4.26)$$

For brevity, we hereby referred to strong stability as stability.

To show that the proposed policy is stabilizing and define the stability region of the system, we introduce artificial decision variables in GREEN and BLUE MILPs. Let $A_{ro} = A_r \cup A_o$ and let $\gamma_i(t) \in [0, 1]$ be a decision variable representing the lane service supply on
lane $i \in A_{ro}$. Recall that $\bar{s}_i = \sum_{j \in A:(i,j) \in M} \bar{s}_{ij}$ is the unconditional lane supply. Consider the lane supply constraints of the form:

$$y_i(t) \leq \gamma_i(t)\bar{s}_i, \quad \forall i \in A_{ro}, \qquad (4.27a)$$

$$\gamma_i(t) \leq \beta_{ij}(t), \quad \forall (i, j) \in \mathcal{M}. \qquad (4.27b)$$

Constraints (4.27) require that $\gamma_i(t)$ be at least $y_i(t)/\bar{s}_i$ and impose that $\gamma_i(t)$ be null if any movement emanating from lane $i$ is inactive. Observe that since $\beta_{ij}(t) = 0 \Rightarrow y_{ij}(t) = 0$, FIFO conditions on lanes imply that if there exists a lane $j$ such that $p_{ij} > 0$ and $\beta_{ij}(t) = 0$, then $y_i(t) = 0$, which is consistent with Constraints (4.27). Observe that Constraints (4.27a) are equivalent to the class-specific lane supply constraints:

$$x_i(t)p_i(t) \leq \gamma_i(t)\bar{s}_i, \quad \forall i \in A_l, \qquad (4.28a)$$

$$\sum_{v \in V_i^n(t)} z_v \leq \gamma_i(t)\bar{s}_i, \quad \forall i \in A_a. \qquad (4.28b)$$

Lane supply constraints (4.27) can be incorporated in GREEN by adding (4.28a) and (4.27b) to GREEN. Similarly, (4.27) can be incorporated in BLUE by adding (4.28b) and (4.27b) to BLUE. Observe that variables $\gamma_i(t) \in [0, 1]$ do not influence the optimal solutions of GREEN or BLUE MILPs since $y_i(t) \leq \bar{s}_i$ for all $i \in A_{ro}$. We henceforth refer to the MILPs with lane supply constraints (4.28) and variables $\gamma_i(t)$, $\forall i \in A_{ro}$ as the extended GREEN and BLUE MILPs.

Let $\Omega_{GREEN}$ and $\Omega_{BLUE}$ denote the feasible regions of the extended GREEN and BLUE MILPs, respectively. Further, let $\Omega = \Omega_{GREEN} \cup \Omega_{BLUE}$ be the feasible region verifying the constraint sets of both GREEN and BLUE, and let $\text{CONV}(\Omega)$ be the convex hull of $\Omega$. We say that an infinite sequence of controls $\{\gamma(t)\}$ for $t = 1, 2, \dots$ is admissible if $\{\gamma(t) \in \Omega\}$ for all time periods $t$.

**Definition 4.5.** For any lane $i \in A_{ro}$, let $\hat{\gamma}_i = \lim_{t \to \infty} \frac{1}{t} \sum_{t=1}^{\bar{t}} \gamma_i(t)$ be the average lane supply of $i$ for sequence $\{\gamma(t) \in \Omega\}$. An admissible control sequence $\{\gamma(t) \in \Omega\}$ accommodates flow $f$ if

$$f_i < \hat{\gamma}_i \bar{s}_i, \quad \forall i \in A_{ro}. \qquad (4.29)$$

Definition 4.5 states that a flow $f$ can be accommodated if, on average, $\hat{\gamma}_i \bar{s}_i$ vehicles can be served on lane $i \in A_{ro}$. From Proposition 2 of Varaiya [2013], $\hat{\gamma} \in \text{CONV}(\Omega)$ if and only if $\{\gamma(t) \in \Omega\}$ is an admissible control sequence. We define the stability region of the system accordingly.

**Definition 4.6.** Let $\mathcal{D}$ be the set of demand rate vectors such that there exists an admissible control sequence $\{\gamma(t) \in \Omega\}$ with $\hat{\gamma} \in \text{CONV}(\Omega)$ defined as:

$$\mathcal{D} = \{ \mathbf{d} \in \mathbb{R}_+^{|\mathcal{A}_r|} : (f_i \leq \hat{\gamma}_i \bar{s}_i, i \in A_{ro}, \{\gamma(t) \in \Omega\}) \}. \qquad (4.30)$$

Let $\mathcal{R}$ be the interior of $\mathcal{D}$ i.e. $\mathcal{R} = \{\mathbf{d} \in \mathcal{D}^\circ\}$. $\mathcal{R}$ is the stability region of the system.

The proposed hybrid max-pressure policy selects at each time period $t$ the service matrix $\alpha(t)$ which maximizes the network-wide pressure among all possible green and blue phases.
**Definition 4.7.** The network traffic control policy $\pi^*(\boldsymbol{x}(t))$, defined as:

$$ \pi^*(\boldsymbol{x}(t)) = \arg\max \left\{ \sum_{n \in \mathcal{N}} (Z_G^n(\boldsymbol{x}(t)) \vee Z_B^n(\boldsymbol{x}(t))) : \boldsymbol{\alpha}(t) \in \mathbb{M} \right\}. \quad (4.31) $$

is hereby referred to as the hybrid max-pressure network control policy for legacy and autonomous vehicles.

At each intersection $n \in \mathcal{N}$, the hybrid max-pressure policy selects the phase (green or blue) and service matrix ($\boldsymbol{\alpha}(t)$) maximizing the local pressure. The service matrix is either $\boldsymbol{\alpha}_G^n(t)$ if the green phase is activated or $\boldsymbol{\alpha}_B^n(t)$ if the blue phase is activated. The implementation of the hybrid max-pressure policy requires the resolution of the two MILPs GREEN and BLUE at each intersection of the network. We next show that policy $\pi^*$ is stabilizing for any demand rate in the stability region of the system $\mathcal{R}$.

**Theorem 4.1.** The stochastic queue evolution process (4.25) is stable under the hybrid max-pressure policy $\pi^*$ (4.31) for any demand rates vector $\boldsymbol{d} \in \mathcal{R}$.

A natural extension of Theorem 4.2 is that the pure network traffic control policies wherein only GREEN or BLUE is used to coordinate traffic are also stable.

**Corollary 4.1.** The pure pressure-based network traffic control policy consisting of policy $\pi^*$ (4.31) with only green (respectively, blue) phases coordinated by GREEN (respectively, BLUE) are stable for any demand rates in the stability region $\mathcal{R}$.

Theorem 4.2 proves that the proposed hybrid network control policy $\pi^*$ (4.31) stabilizes any demand vector in the stability region of the system $\mathcal{R}$ (4.30). According to Tassiulas and Ephremides [1992], and as discussed in Wongpiromsarn et al. [2012] and Varaiya [2013] for the case of signalized traffic control, this is equivalent to throughput optimality since it shows that the stability region of policy (4.31) is a superset of the stability region of any policy.

Corollary 4.1 establishes that pure policies based on GREEN or BLUE traffic control models also maximize throughput. We note that stability of the pure green network traffic control case is an extension of the work of Varaiya [2013]. Varaiya [2013] proposed a network traffic control policy for a single class of vehicles and assumed that each movement had a dedicated queue. In addition, it was assumed that movement capacities are exogenous to the traffic signal control policy. We have both relaxed this framework by only requiring knowledge of lane-queues, and extended the formulation to two classes of lanes. Further, we introduced pressure-based formulations for both green and blue phases which account for FIFO blocking effects on upstream lanes, and account for the loss of capacity due to conflicting movements being simultaneously activated.

We are now ready to present our decentralized network traffic control algorithm used to implement the proposed hybrid max-pressure network control policy. The pseudo-code of the proposed policy is summarized in Algorithm 6. At each time period $t$, we calculate the optimal GREEN and BLUE phases at each intersection of the network $n \in \mathcal{N}$ based on the current state of the network $\boldsymbol{x}(t)$. The phase with the highest local pressure is selected for each intersection.
<table><caption>**Algorithm 6:** Hybrid max-pressure network control policy</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Input:</td><td><i>G</i> = (<i>N</i>, <i>A</i>), <b><i>d</i></b>, <b><i>p</i></b>, <b><i>s̄</i></b>, <i>t</i>, <b><i>x</i></b>(<i>t</i>)</td></tr><tr><td>Output:</td><td><i>α</i>(<i>t</i>)</td></tr><tr><td>1</td><td>for <i>n</i> ∈ <i>N</i> do</td></tr><tr><td>2</td><td>&nbsp;&nbsp;&nbsp;&nbsp;<i>Z</i><sup><i>n</i></sup><sub><i>G</i></sub>(<b><i>x</i></b>(<i>t</i>)) ← Solve GREEN</td></tr><tr><td>3</td><td>&nbsp;&nbsp;&nbsp;&nbsp;<i>Z</i><sup><i>n</i></sup><sub><i>B</i></sub>(<b><i>x</i></b>(<i>t</i>)) ← Solve BLUE</td></tr><tr><td>4</td><td>&nbsp;&nbsp;&nbsp;&nbsp;<b><i>α</i></b><sup><i>n</i></sup>(<i>t</i>) ← arg max<sub><b><i>α</i></b><sup><i>n</i></sup><sub><i>G</i></sub>(<i>t</i>),<b><i>α</i></b><sup><i>n</i></sup><sub><i>B</i></sub>(<i>t</i>)</sub>{<i>Z</i><sup><i>n</i></sup><sub><i>G</i></sub>(<b><i>x</i></b>(<i>t</i>)), <i>Z</i><sup><i>n</i></sup><sub><i>B</i></sub>(<b><i>x</i></b>(<i>t</i>))}</td></tr><tr><td>5</td><td><b><i>α</i></b>(<i>t</i>) ← [<b><i>α</i></b><sup><i>n</i></sup>(<i>t</i>)]<sub><i>n</i>∈<i>N</i></sub></td></tr></tbody></table>

### 4.4.3 Numerical experiments

We implement the proposed hybrid network control policy on artificial datasets to test computational performance and analyze the algorithm's behavior. We use a synthesized grid network of size 5 × 5, wherein each of the 25 nodes corresponds to a controlled intersection and each edge represents a bidirectional link between adjacent nodes. All intersections have the same topology as that depicted in Figure 4.2a, i.e. each node has four incoming and four outgoing links, each of which has one LV-lane and one AV-lane. Each incoming lane allows three movements: through, left turn and right turn.

<figure><img src="image_11.png" /><figcaption>Figure 4.2: Figure 4.2a shows all possible LV- and AV-movements on a typical intersection. Figure 4.2b depicts a possible green phase involving priority and yield LV-movements and Figure 4.2c illustrates the conflict-point formulation during blue phase.</figcaption></figure>

We assume that vehicles' routes in the network are fixed. In each instance generated, we randomly and uniformly assign an origin and a destination to each vehicle, a route among these nodes and a departure time within the considered time horizon. Origins and destinations are chosen among nodes at the edge of the grid. The level of travel demand is determined by the *departure rate* of vehicles into the network and the impact of travel demand onto network performance is assessed through a sensitivity analysis.

The time period is set to $\Delta t = 10$ s. We assume that green phases have a *lost time* of 2 s to account for vehicle start-up delays and signal clearance intervals and we conduct a sensitivity analysis on this input parameter. In turn, blue phases are assumed to have zero lost time.
<figure><img src="image_12.png" /><figcaption>Figure 4.3: Total system travel time for a varying proportion of AVs and vehicle departure rate. The results illustrate the trend of the mean total system travel time over 40 simulations on a 5 × 5 grid network with each link having one AV and one LV lane. The 2×Green experiment (benchmark) corresponds to the scenario where each link has the capacity of two LV-lanes.</figcaption></figure>

We use point-queues for all links in the network and we assume that vehicles take three time periods to travel from an intersection to the next intersection on their route. Vehicles travel through 3 links, each with a 10 s free flow time, between each intersection. Hence, in this configuration, it takes 30 s for a vehicle to travel between two adjacent intersections at free flow. We set the time horizon to 30 minutes and we execute Algorithm 6 periodically until all vehicles have exited the network. Vehicles’ speed limit through intersections is assumed to be uniform and equal to 44 ft/s and the wave propagation speed is taken as 11 ft/s. We assume that all vehicles have a length of 17.6 ft and that lanes have a width of 12 ft. We use the triangular fundamental diagram to determine lane capacity which results in 1,440 veh/h or 5 vehicles per time period.

We explore the sensitivity of the proposed policy with regards to the proportion of AVs by varying the proportion of AVs from 0% to 100% in increments of 10%. In all our numerical experiments we assume that AVs always choose AV-lanes. To benchmark the performance of the proposed hybrid network control policy (summarized in Algorithm 6), we also simulate network traffic under a traditional traffic signal configuration wherein AV-lanes and blue phases are nonexistent. Under this configuration all AV-lanes are treated as LV-lanes, and we model this by using single-lane links with twice the lane capacity of the LV-lanes in the network configuration with both LV- and AV-lanes. Thus, in this benchmark there are no AV-lanes and since we assume that AVs behave as LVs and LV-lanes, all vehicles are treated as LVs. This benchmark is hereby referred to as 2×Green. Each experiment is simulated 40 times and average performance is reported.

The simulation framework is implemented in Java on a Windows machine with 8 Gb of RAM and a CPU at 3 GHz. All MILPs are solved with CPLEX 12.8 (Java API) with a time limit of 60 s and default options.
<figure><img src="image_13.png" /><figcaption>Figure 4.4: Average vehicle travel time based on AV proportion and vehicle departure rate. The results depict the mean and standard deviation over 40 simulations on a 5×5 grid network with each link having one AV and one LV lane. The 2×Green experiment (benchmark) corresponds to the scenario where each link has the capacity of two LV-lanes.</figcaption></figure>
The evolution of the total system travel time (TSTT) for a varying departure rate is depicted in Figure 4.3. For this experiment, the green phase lost time is set to 2 s and we vary the departure rate from 4,000 veh/h (lowest demand) to 10,000 (highest demand). As expected, we observe that TSTT increases super-linearly with the departure rate, i.e. travel demand. We find that the market penetration of AVs has a significant impact on TSTT. If less than 50% of AVs are present in the system, we find that the use of dedicated AV lanes and the blue phase is not beneficial for the network in terms of TSTT, even at high demands, as noted by the performance of the benchmark which outperforms the hybrid network control policy for these levels of AV market penetration. Recall that the benchmark represents the TSTT when the pure, green network traffic control policy is used with the equivalent of two LV-lanes. At high demands (i.e. 8,000 veh/h and beyond), we find that a market penetration of more than 50% of AVs outperforms the benchmark. At the highest departure rate tested (i.e. 10,000 veh/h), for a market penetration of 70%, the average TSTT is reduced by 1/2 compared to the benchmark. For AV market penetration rates of 80% and higher, we find that the reduction of the average TSTT is lesser than that achieved at 70% of AV market penetration. This can be explained by the fact that network capacity is not fully used at 100% of AVs compared to 70% since in the former configuration LV-lanes remain empty throughout the experiment.

To further investigate the behavior of the proposed hybrid network control policy, we examine average, vehicle-class travel times relative to the benchmark configuration. Figure 4.4 shows the average vehicle travel time for AVs (blue series), LVs (green series) and overall (orange series) in the network based on AV proportion and departure rate. The benchmark is shown as a dashed flat line in red. Three main trends can be identified: first, we find that increasing the departure rate mainly impacts the mean relative travel time of LVs with regards to the benchmark. Second, the proportion of AVs minimizing the mean vehicle travel time in the network is relatively robust to the travel demand and is found to be around 70% of AVs. This suggests that at lower or higher market penetration rates, network capacity is not saturated and AV- or LV-lanes are under-utilized, respectively. This insight can help in identifying the optimal AV market penetration rate to deploy AV lanes and blue phases. Third, for high levels of travel demand, a sufficiently high proportion of AVs improves on the benchmark, i.e. the average vehicle travel time (over both LVs and AVs) obtained using the hybrid network control policy is lower than the average vehicle travel time obtained with the pure, green network control policy.

For a departure rate of 4,000 veh/h, LVs and AVs' average travel time remain similar and the hybrid network control policy performs similarly to the benchmark. Increasing the departure rate from 5,000 to 7,000 veh/h, we observe congestion effects impacting LVs' average travel time, while AVs' average travel time remain only slightly penalized. Further increasing the departure rate to 8,000 veh/h and beyond yields a new pattern: for a proportion of AVs greater or equal to 60%, the aggregate average vehicle travel time improves on the benchmark although LVs' average travel time remain more penalized than that of AVs. At a departure rate of 10,000 veh/h, we find that both LVs and AVs' average travel time are lower than the average travel time in the benchmark configuration and considerable travel time reduction are achieved from 60% of AVs and beyond. Specifically, at this departure rate and with a market penetration of AVs of 70%, we find that the average vehicle travel time is reduced by approximately 1/2 compared to the benchmark. We also observe that the mean relative vehicle (LVs and AVs) travel time exhibits a convex-shaped profile with a minima at 70% of AV market penetration.
## 4.5 Online mechanism design for traffic intersection auctions

We now turn our attention to market-based AIM mechanisms and present two online payment mechanisms for traffic intersection auctions. For this study, we restrict ourselves to a simplified intersection model and focus on the game-theoretical aspects of dynamic user pricing.

### 4.5.1 Dynamic traffic intersection auction framework

We consider a traffic intersection with a finite number of access lanes. We consider discrete time and we assume that at each time period, users arrive on access lanes with a known probability. Upon arriving at the front of their lane, users are requested to declare their delay cost which will be used as their bid in a combinatorial auction. We assume that users seek to minimize a generalized cost function which is a linear combination of their expected waiting time and their payment to the intersection manager. We assume that users are serviced sequentially by a single server, i.e. the traffic intersection and that service is nonpreemptive, i.e. cannot be interrupted once started. This does not prevent several users to simultaneously traverse the intersection. Instead, the proposed traffic intersection auction only aims to determine the sequence in which users are to be serviced and the payments they should be charged based on their declared preferences.

The goal of the intersection manager is to maximize social welfare which is defined as the total generalized user cost. To achieve optimal social welfare, we seek to determine incentive-compatible user payments to ensure the truthful declaration of their delay costs. To ensure truthful user behavior, traffic intersection lanes are serviced in order of decreasing declared user delay costs. The expected waiting time of users is a function of the declared delay costs and the state of the system. In turn, the payment of users is to be determined by the intersection manager so as to maximize social welfare and ensure truthfulness.

To present the proposed online mechanism for determining incentive-compatible payments, we consider the case of a user $i$ arriving at the front of its lane and declaring a delay cost of $v_i$. Let $Q$ be the number of access lanes of the intersection. We denote $\bar{q}_i$ (resp. $\underline{q}_i$) the number of lanes which are occupied by users with greater (resp. lower) delay cost than $i$. Further, we denote $q_0$ the number of empty lanes. Summing $\bar{q}_i$, $\underline{q}_i$ and $q_0$ corresponds to the $Q - 1$ lanes of the intersection which are not occupied by user $i$. Hence, the following equation holds for any user $i$ at any time period:

$$Q = \underline{q}_i + \bar{q}_i + q_0 + 1. \qquad (4.32)$$

Accordingly, we can define the pricing queue of the proposed dynamic traffic intersection auction.

**Definition 4.8 (Pricing queue).** In a traffic intersection with $Q$ access lanes, the pricing queue is a dynamic set of at most $Q$ users which represents the users at the front of their lane at a given time period.

The pricing queue consists of the set of users which have no users in front of them in their access lane to the intersection. In the proposed traffic intersection auction, only users in the pricing queue are participants. This implies that users who are behind other users
are not eligible to participate until they reach the front of their lane. This auction design aims to obviates blocking effects induced by non-overtaking conditions in access lanes of the intersection. Under such typical, non-overtaking conditions users waiting behind other users in their lane cannot bypass users in front of them, thus preventing them to be serviced before reaching the front of their lane. In the proposed approach, only users in the pricing queue participate in the auction. While this dynamic auction framework restricts the number of participants to at most Q users, every user queueing in the intersection will eventually reach the front of its lane and participate in a single instance of the auction. This is in contrast to approaches in the literature where auction “losers” may participate in several auctions. We refer readers to Carlino et al. [2013] who presented pricing mechanisms for traffic intersection auctions in which blocked users are allowed to “vote” for their lane.

The objective of the traffic intersection auction is to determine incentive-compatible user payments upon users joining the pricing queue. Let $W_i(v_i)$ be the expected waiting time of user $i$ declaring a delay cost of $v_i$. Let $v_i^*$ be the true delay cost of user $i$, the expected waiting time cost of user $i$ declaring a delay cost of $v_i$ is $v_i^*W_i(v_i)$. Let $P_i(v_i)$ be the payment of user $i$ declaring $v_i$, which is to be determined by the intersection manager. We denote $C_i(v_i)$ the generalized cost for user $i$ and define the user objective function as:

**Definition 4.9 (User objective function).** The objective function of user $i$ is:

$$ \min_{v_i} C_i(v_i) = \min_{v_i} (v_i^* W_i(v_i) + P_i(v_i)). \qquad (4.33) $$

Following the approach of Dolan [1978], we will show that setting the payment $P_i(v_i)$ equal to the expected marginal delay cost user $i$ imposes to other users is incentive-compatible, i.e. the user objective is minimal for $v_i = v_i^*$.

We assume that the intersection manager has knowledge of the probability distribution of users’ valuation (delay costs), and of the arrival probability of users on access lanes of the intersection. These assumptions are reasonable since one can assume that the intersection manager can observe user behavior (arrival rate and valuations) and learn these probabilities over time. Let $F(v_i)$ be the cumulative distribution function representing the probability that a user declares a delay cost $v \le v_i$. Further, let $\underline{v}$ and $\bar{v}$ be lower and upper bounds on users’ delay costs. For conciseness, we hereby use the term lower-bidding user (resp. higher-bidding user) to refer to a user with a declared delay cost lower (resp. higher) than that of user $i$.

An overview of the proposed dynamic auction process is summarized in Figure 1. Users are assumed to arrive dynamically and submit their bids (delay cost) to the intersection manager. The intersection manager then determines a user payment $P_i(v_i)$ based on the user’s bid $v_i$ using either the proposed queue- or lane-based model. The user is then charged $P_i(v_i)$ and moves through the intersection.

### 4.5.2 Online incentive-compatible mechanisms

We propose two online mechanisms to determine the expected waiting time and payment of users. The proposed mechanisms differ in the way the expected waiting time of users is determined. We first consider a queue-based model which only tracks the number of users bidding higher or lower than the reference user (Section 4.5.2). We next extend this model to a lane-based approach which tracks users’ arrival lane (Section 4.5.2). A generic
<figure><img src="image_14.png" /><figcaption>Figure 4.5: Online traffic intersection auction process</figcaption></figure>

payment mechanism is proposed to determine incentive-compatible payments for both queue-
and lane-based models which result in two online mechanisms (Section 4.5.3). The queue-
based mechanism assumes that the arrival probability of users is uniform across lanes of
the intersection and provides a relatively low-dimensional mathematical framework. This
uniformness assumption is relaxed in the lane-based mechanism which can accommodate
non-uniform lane arrival probabilities at the expense of an extended state space.

Queue-based model

We use a MC to determine the expected waiting time of user *i* under different system states.
The state of the queuing system w.r.t. user *i* is represented by the number of lanes occupied
by lower bidding vehicles lanes *q*<sub>*i*</sub> the number of empty lanes *q*<sub>empty</sub>. We denote *q*<sub>*i*</sub> the vector
*q*<sub>*i*</sub> = (*q*<sub>*i*</sub>, *q*<sub>empty</sub>) characterizing the state of the system w.r.t. user *i*. User *i* is serviced when
all lanes are either occupied by a lower-bidding user or empty. Observe that once a lane is
occupied by a lower-bidding user w.r.t. user *i*, it remains in this state until *i* is serviced.
Hence, from the perspective of user *i*, the state space of the MC can be characterized by the
inequality *q*<sub>*i*</sub> + *q*<sub>empty</sub> ≤ *Q* − 1. This leads to the following characterization of the state space.

**Definition 4.10 (State space in the queue-based model).** In the queue-based model, the state
space of the MC for user i declaring a delay cost of $v_i$ in a pricing queue of size $Q$ is:

$$S_i^q = \{ \mathbf{q}_i = (\underline{q}_i, q_\emptyset) \in \mathbb{N}^2 : \underline{q}_i + q_\emptyset \le Q - 1 \}. \quad (4.34)$$

For each state $\mathbf{q}_i \in S_i^q$, we are interested in determining the cost-to-go from this state through the MC. If $\underline{q}_i + q_\emptyset = Q - 1$, then (4.32) yields $\bar{q}_i = 0$ which is equivalent to require that the pricing queue does not contain any user with declared delay costs higher than that of i. Hence, if $\underline{q}_i + q_\emptyset = Q - 1$, then the MC converges and user i is serviced. Hence, to determine the state transition probabilities from state $\mathbf{q}_i$ to state $\mathbf{q}_i'$, we assume that $\underline{q}_i + q_\emptyset < Q - 1$.

To calculate the expected waiting time of user i in the queue-based model, we enumerate system states with varying number of lower-bidding and empty lanes w.r.t. user i. We wish to define a function $f_q$ which gives the transition between the state at time step $t$ to the state at time $t+1$. The transition function of the MC, $f_q$, depends on the realization of random variables. For this problem, there are two random variables affecting the transition: user delay costs, a continuous random variable; and whether a new vehicle arrives on each lane, a boolean value. Hence, the uncertain data of the problem can be represented by a vector $\xi \in \mathbb{R}^Q \times 2^Q$ where $\mathbb{R}^Q$ represents user delays costs and $2^Q$ represents vehicle arrivals. Since the state space of the queue-based model is $S_i^q$, the transition function maps the current state $\mathbf{q}_i \in S_i^q$ and the uncertain data $\xi \in \mathbb{R}^Q \times 2^Q$ to the next state, i.e. $f_q : S_i^q \times \mathbb{R}^Q \times 2^Q \to S_i^q$. We denote $g_q(\mathbf{q}_i)$ the one-step-cost of the MC representing the service time of the system in state $\mathbf{q}_i$ and we assume that the one-step-cost is deterministic. This assumption is plausible since users' behavior is assumed to be deterministic, i.e. users cannot change their declared delay cost over time.

Let $W_q(v_i, \mathbf{q}_i)$ be the expected waiting time of user i declaring $v_i$ for state $\mathbf{q}_i$, which represents the cost-to-go of the MC representing the state evolution in the queue-based model. Let $\mathcal{T}_i^q = \{\mathbf{q}_i \in S_i^q : q_\emptyset + \underline{q}_i = Q - 1\}$ be the set of terminal states in the queue-based mechanism. The expected waiting time of user i declaring $v_i$ for state $q$ can be determined as:

$$W_q(v_i, \mathbf{q}_i) = \begin{cases} 0 & \text{if } \mathbf{q}_i \in \mathcal{T}_i^q, \\ \mathbb{E}[W_q(v_i, f_q(\mathbf{q}_i, \xi))] + g_q(\mathbf{q}_i) & \text{otherwise.} \end{cases} \quad (4.35)$$

Note that the term $\mathbb{E}[W_q(v_i, f_q(\mathbf{q}_i, \xi))]$ represents the expected value of the expected waiting time $W_q(v_i, f_q(\mathbf{q}_i, \xi))$ after the transition function $f_q$. If $\mathbf{q}_i$ is not a terminal state, i.e. $\underline{q}_i + q_\emptyset < Q - 1$, the expected waiting time of user i can be calculated using the probability $\text{Pr}[\mathbf{q}_i'|\mathbf{q}_i]$ of transitioning from state $\mathbf{q}_i$ to state $\mathbf{q}_i'$. The expected waiting time of user i declaring $v_i$ is:

$$\mathbb{E}[W_q(v_i, f_q(\mathbf{q}_i, \xi))] = \sum_{\mathbf{q}_i' \in S_i^q} \text{Pr}[\mathbf{q}_i'|\mathbf{q}_i] W_q(v_i, \mathbf{q}_i'). \quad (4.36)$$

Combining Eqs. (4.35) and (4.36) yields the following system of linear equations:

$$W_q(v_i, \mathbf{q}_i) = 0, \quad \forall \mathbf{q}_i \in \mathcal{T}_i^q, \quad (4.37a)$$

$$W_q(v_i, \mathbf{q}_i) (1 - \text{Pr}[\mathbf{q}_i|\mathbf{q}_i]) - \sum_{\mathbf{q}_i' \in S_i^q : \mathbf{q}_i' \neq \mathbf{q}_i} W_q(v_i, \mathbf{q}_i') \text{Pr}[\mathbf{q}_i'|\mathbf{q}_i] = g_q(\mathbf{q}_i), \quad \forall \mathbf{q}_i \in S_i^q \setminus \mathcal{T}_i^q. \quad (4.37b)$$
Let $p$ be the probability that a user arrives on a lane of the intersection at the next time period. Recall that $F(v_i)$ is the cumulative distribution function of users bids. On each lane, three events may occur with the following probabilities: i) a lane can become or remain empty with probability $(1-p)$; ii) a lower-bidding user can arrive with probability $pF(v_i)$; and iii) a higher-bidding user can arrive with probability $p(1-F(v_i))$. Observe that there are $(q_0^{q_0+1})$ ways to choose $q'_0$ empty lanes among $q_0+1$ lanes. There are $q'_i - q_i$ lanes which become occupied by lower-bidding users and there are $(q'_i - q'_0)$ ways to choose $q'_i - q_i$ of the remaining $q_0 + 1 - q'_0$ non-empty lanes to become occupied by a lower-bidding user. Accordingly, the transition probability from state $(q_i, q_0)$ to state $(q'_i, q'_0)$ is:

$$ \mathrm{Pr}[\boldsymbol{q}'_i | \boldsymbol{q}_i] = \binom{q_0+1}{q'_0} (1-p)^{q'_0} \binom{q_0+1-q'_0}{q'_i-q_i} (pF(v_i))^{q_0+1-q'_0-q'_i+q_i}. \quad (4.38) $$

The expected waiting time defined in Eq. (4.35) is a recursive function which can be calculated by solving a series of systems of linear equations where each equation corresponds to a possible state $\boldsymbol{q}_i \in S_i^q$. Since all states such that $q_i + q_0 = Q - 1$ are terminal, their corresponding expected waiting time is null and these variables can be eliminated from the system of equations. Hence, we need only to consider the set of pairs of nonnegative integers $(q_i, q_0)$ such that $q_i + q_0 \le Q - 2$. Observe that each variable $q_i$ or $q_0$ can take $Q - 2 + 1$ values, thus the number of equations in the system is $\sum_{k=0}^{Q-1} k = \binom{Q}{2}$.

Further, observe that $q_i$ can never decrease while $q_0$ can increase or decrease. Hence, for any state $\boldsymbol{q}_i = (q_i, q_0)$, we can calculate $W_q(v_i, \boldsymbol{q}_i)$ in a dynamic programming fashion starting with the maximum number of lanes, i.e. $k = Q - 1$, and iterating downwards until $k = q_i$. At each iteration $k$, we solve the system of $Q-k$ linear equations with $Q-k$ variables corresponding to states $\boldsymbol{q}_i = (k, l)$ for $0 \le l \le Q - 1 - k$.

## Lane-based model

We now consider an alternative approach where lane arrival probabilities can be non-uniform across lanes of the intersection. We abuse notation and denote $p_j$ the probability of a user arriving on lane $j$ at the next time period. Let $\mathcal{L}_i = \{\sigma_0, \sigma_{\bar{v}_i}, \sigma_{\bar{v}_i}\}$ be the set of possible lane-based states with regards to user $i$ bidding $v_i$; where $\sigma_0$ represents an empty lane, $\sigma_{\bar{v}_i}$ the arrival of a user bidding lower than $v_i$ and $\sigma_{\bar{v}_i}$ the arrival of a user bidding higher than $v_i$. Let $z_{ij} \in \mathcal{L}_i$ be the state of lane $j$ w.r.t. user $i$. We denote $z_i = (z_{i,1}, \dots, z_{i,Q-1})$ the state w.r.t. user $i$ in the lane-based model.

**Definition 4.11 (State space for lane-based model).** In the lane-based model, the state space of the MC for user $i$ declaring a delay cost of $v_i$ in a pricing queue of size $Q$ is

$$ S_i^z = \{z_i = (z_{i,j_1}, \dots, z_{i,j_{Q-1}}) \in \mathcal{L}_i^{Q-1}\}. \quad (4.39) $$

In the lane-based model, the number of possible states w.r.t. user $i$ in an intersection with $Q$ lanes is $|\mathcal{L}_i|^{Q-1}$, hence a four-lane intersection has $3^3 = 27$ possible states.

Let $f_z : S_i^z \times \mathcal{R}^Q \times 2^Q \to S_i^z$ be the MC transition function of the lane-based model and let $g_z(\boldsymbol{z}_i)$ be the one-step cost. Let $\mathcal{T}_i^z \subset S_i^z$ be the set of terminal states, i.e. $\mathcal{T}_i^z = \{\boldsymbol{z}_i \in$
$S_i^z : z_{ij} = \sigma_0 \lor z_{ij} = \sigma_{v_i}, j = 1, \dots, Q-1 \}$. The expected waiting time $W_z(v_i, z_i)$ of user $i$ declaring $v_i$ for state $z_i$ in the lane-based model is thus:

$$W_z(v_i, z_i) = \begin{cases} 0 & \text{if } z_i \in \mathcal{T}_i^z, \\ \mathbb{E}[W_z(v_i, f_z(z_i, \xi))] + g_z(z_i) & \text{otherwise,} \end{cases} \quad (4.40)$$

with:

$$\mathbb{E}[W_z(v_i, f_z(z_i, \xi))] = \sum_{z'_i \in S_i^z} \mathrm{Pr}[z'_i|z_i] W_z(v_i, z'_i). \quad (4.41)$$

This leads to the following system of linear equations:

$$W_z(v_i, z_i) = 0 \quad \forall z_i \in \mathcal{T}_i^z, \quad (4.42a)$$

$$W_z(v_i, z_i)(1 - \mathrm{Pr}[z_i|z_i]) - \sum_{z'_i \in S_i^z: z'_i \neq z_i} W_z(v_i, z'_i) \mathrm{Pr}[z'_i|z_i] = g_z(z_i) \quad \forall z_i \in S_i^z \setminus \mathcal{T}_i^z. \quad (4.42b)$$

Unlike the systems of equations (4.37), the system (4.42) does not admit an intuitive decomposition which could be used to solve the system of equations via a recursive algorithm. Thus, in our numerical experiments, (4.42) is solved in a single step using traditional linear algebra codes.

To determine transition probabilities in the lane-based mechanism, it is necessary to track which user will be serviced at the next time period. For any non-terminal state, there is at least one higher-bidder in the pricing queue, i.e. $\bar{q}_i \ge 1$, and the highest-bidder in the queue will be serviced next. Let $V_{>i}$ be the set of higher-bidders, i.e. $V_{>i} = \{j_1, j_2, \dots, j_{\bar{q}_i}\}$ with $v_i < v_{j_1} < v_{j_2} < \dots < v_{j_{\bar{q}_i}}$. If there is more than a single user bidding higher than user $i$, i.e. $\bar{q}_i > 1$, it is unknown which of these higher-bidders will be the highest-bidder. Since users' delay costs are assumed to all follow the same probability distribution, the probability that a higher-bidder is the highest-bidder is uniform, i.e. $1/\bar{q}_i$.

Let $z_{ij}(k)$ be the state of lane $j$ with regards to user $i$ assuming lane $k$ is the moving lane, i.e. the lane occupied with the highest-bidder in the pricing queue. Let $\mathrm{Pr}[z'_{ij}|z_{ij}(k)]$ be the transition probability of lane $j$ with regards to user $i$ from state $z_{ij}(k)$ to state $z'_{ij}$. Since the arrival process of a lane is assumed to be independent of that of other lanes, the transition probability from state $z_i$ to $z'_i$ denoted $\mathrm{Pr}[z'_i|z_i]$ is

$$\mathrm{Pr}[z'_i|z_i] = \mathrm{Pr}[(z'_{i,1}, \dots, z'_{i,Q-1})|(z_{i,1}, \dots, z_{i,Q-1})] = \frac{1}{\bar{q}_i} \sum_{k \in V_{>i}} \prod_{j=1}^{Q-1} \mathrm{Pr}[z'_{ij}|z_{ij}(k)]. \quad (4.43)$$

In the lane-based mechanism, transition probabilities are function of current lane state. If lane $j$ is empty or if lane $j$ is occupied by the highest bidder in the pricing queue, then we say that this lane is open. Otherwise, lane $j$ is either occupied by a lower-bidding user or a higher-bidding user who will not be serviced next and we say that his lane is closed. Let $V_{-i}$ be the set of users in the pricing queue other than user $i$. Let $O_{i,j} \subseteq L_i$ and $C_{i,j} \subseteq L_i$ be the
set of open and closed states for lane j with regards to user i, respectively:

$$ \mathcal{O}_{i,j} = \begin{cases} \{\sigma_0, \sigma_{\bar{v}_i}\} & \text{if } j \in \arg\max\{v_k : k \in \mathcal{V}_{-i}\}, \\ \{\sigma_0\} & \text{otherwise,} \end{cases} \qquad (4.44) $$

$$ \mathcal{C}_{i,j} = \begin{cases} \{\sigma_{\underline{v}_i}\} & \text{if } j \in \arg\max\{v_k : k \in \mathcal{V}_{-i}\}, \\ \{\sigma_{\underline{v}_i}, \sigma_{\bar{v}_i}\} & \text{otherwise.} \end{cases} \qquad (4.45) $$

For any lane j in state $\sigma \in \mathcal{O}_{i,j}$, we have the following transition probabilities:

$$ \Pr[\sigma_0|\sigma] = 1 - p_j \qquad \forall j \in \{1, \dots, Q-1\}, \forall \sigma \in \mathcal{O}_{i,j}, \qquad (4.46a) $$

$$ \Pr[\sigma_{\underline{v}_i}|\sigma] = p_j F(v_i) \qquad \forall j \in \{1, \dots, Q-1\}, \forall \sigma \in \mathcal{O}_{i,j}, \qquad (4.46b) $$

$$ \Pr[\sigma_{\bar{v}_i}|\sigma] = p_j(1 - F(v_i)) \qquad \forall j \in \{1, \dots, Q-1\}, \forall \sigma \in \mathcal{O}_{i,j}. \qquad (4.46c) $$

Otherwise if lane j is in state $\sigma \in \mathcal{C}_{i,j}$, this lane remains in its current state with probability one, since the corresponding user cannot be serviced at the next time period.

$$ \Pr[\sigma'|\sigma] = \begin{cases} 1 & \text{if } \sigma' = \sigma, \\ 0 & \text{otherwise,} \end{cases} \qquad \forall j \in \{1, \dots, Q-1\}, \forall \sigma \in \mathcal{C}_{i,j}. \qquad (4.47) $$

Using these lane-based transition probabilities, we can determine intersection-based transition probabilities via Eq. (4.43).

The queue-based model can be viewed as a special case of the lane-based model. Let $\mathbf{1}_{z_{ij}=\sigma}$ be the indicator function taking value 1 if $z_{ij} = \sigma$ and 0 otherwise. We first define the mapping between the state spaces of both queue- and lane-based models.

**Definition 4.12.** Let $h_i : S_i^z \to S_i^q$ be a function mapping the state space of the lane-based model to the state space of the queue-based model: $h_i : \mathbf{z}_i \mapsto \mathbf{q}_i = h_i(\mathbf{z}_i)$ with $\underline{q}_i = \sum_{j=1}^{Q-1} \mathbf{1}_{z_{ij}=\sigma_{\underline{v}_i}}$ and $q_0 = \sum_{j=1}^{Q-1} \mathbf{1}_{z_{ij}=\sigma_0}$.

The proposition below highlights the correspondence between queue- and lane-based models.

**Proposition 4.2.** If lane arrival probabilities are uniform, i.e. $p = p_j$ for all $j = 1, \dots, Q-1$, and if the one-step-costs of the queue- and lane-based models are such that $g_z(\mathbf{z}_i) = g_q(\mathbf{q}_i)$ for all states $\mathbf{z}_i \in S_i^z$ and $\mathbf{q}_i \in S_i^q$ such that $\mathbf{q}_i = h_i(\mathbf{z}_i)$, then $W_z(v_i, \mathbf{z}_i) = W_q(v_i, \mathbf{q}_i)$ for all such states.

Proposition 4.2 establishes a correspondence between the proposed queue- and lane-based models for determining the expected waiting times of users and shows that under uniform lane arrival probabilities and one-step-costs, both models are equivalent.

We next present a payment mechanism that can be used with both of these models to determine incentive-compatible payments.
### 4.5.3 Payment mechanism

The proposed payment mechanism can be used with both queue- and lane-based models for determining the expected waiting time of users. Hence, in this section, we simply denote $W_i(v_i)$ the expected waiting time of user $i$ bidding $v_i$. This expected waiting time can be replaced by $W_q(v_i, q_i)$ for the queue-based model, or by $W_z(v_i, z_i)$ for the lane-based model. To determine incentive-compatible user payments, we calculate the expected marginal cost user $i$ imposes on other users. For this, we first examine the period of time over which an extra user in the lane of user $i$ is expected to have an impact on other users.

The busy period w.r.t. user $i$ is the expected time required to clear the queues if $i$ declares a minimal delay cost, i.e. $v_i = \underline{v}$. Observe that if $v_i = \underline{v}$, then $q_i = 0$, since no user can bid lower than $\underline{v}$ i.e. $F(\underline{v}) = 0$. We measure the impact of user $i$ onto other users throughout the remaining busy period, i.e. the duration between the start of service of user $i$ if $i$ declares $v_i$, and the start of service of $i$ if $i$ had declared a minimal delay cost. Hence, the remaining busy period w.r.t. user $i$ bidding $v_i$ denoted $\Delta W_i(v_i)$ can be defined as:

$$ \Delta W_i(v_i) = W_i(\underline{v}) - W_i(v_i). \qquad (4.48) $$

To determine the marginal impact of user $i$ on other users, we can split the remaining busy period w.r.t. user $i$ into two components: the expected portion of $\Delta W_i(v_i)$ which affects users who arrived before user $i$; and the expected portion of $\Delta W_i(v_i)$ which affects users who are expected to arrive after user $i$. We denote $B_i(v_i)$ and $A_i(v_i)$ these quantities, respectively, and intuitively we define

$$ \Delta W_i(v_i) = B_i(v_i) + A_i(v_i). \qquad (4.49) $$

Calculating $B_i(v_i)$ and $A_i(v_i)$ correctly is not trivial, especially the latter since the exact number of users arriving after $i$ that will be delayed by user $i$ is unknown. We next propose a method to determine $B_i(v_i)$, then use Eq. (4.49) together with Eq. (4.48) to obtain $A_i(v_i)$.

To determine the marginal delay caused by user $i$ onto users already in the pricing queue, we need only to consider those $q_i$ users which bid lower than $i$ since others users are not delayed by user $i$. Let $\mathcal{V}_{<i}$ be the set of all users bidding lower than $i$ ordered by increasing declared delay costs, i.e. $\mathcal{V}_{<i} = \{j_1, j_2, \dots, j_{q_i}\}$ with $v_{j_1} < v_{j_2} < \dots < v_{j_{q_i}} < v_i$. To determine the expected portion of the remaining busy period w.r.t. user $i$ which affects user $j \in \mathcal{V}_{<i}$, we take the perspective of an extra user bidding $v_j$ in the lane of user $i$. The expected period of time during which this extra user affects user $j$ is the difference between the waiting time of the extra user if the extra user is bypassed by user $j$ and that if the extra user bypasses user $j$ in the pricing queue. We abuse notation and denote $W_i(v_j : z_{ij} = \sigma_{\bar{v}_i})$ and $W_i(v_j : z_{ij} = \sigma_{\underline{v}_i})$ these expected waiting times, respectively.

We propose the following formula for $B_i(v_i)$:

$$ B_i(v_i) = \sum_{k=1}^{q_i} W_i(v_{jk} : z_{ij_k} = \sigma_{\bar{v}_i}) - W_i(v_{jk} : z_{ij_k} = \sigma_{\underline{v}_i}). \qquad (4.50) $$

The definitional relationship linking the time periods $B_i(v_i)$ and $\Delta W_i(v_i)$, i.e. Eq. (4.49), implicitly assumes that $B_i(v_i) \le \Delta W_i(v_i)$. However, it is not obvious that the right-hand side of Eq. (4.50) is never greater than $\Delta W_i(v_i)$. Since this is necessary for the correctness of the proposed formula for the time period $B_i(v_i)$, i.e. Eq. (4.50), we prove this relationship in Proposition 4.3.
**Proposition 4.3.** For any user *i* bidding $v_i$,

$$ \sum_{k=1}^{q_i} W_i(v_{jk} : z_{ijk} = \sigma_{\bar{v}_i}) - W_i(v_{jk} : z_{ijk} = \sigma_{\underline{v}_i}) \le \Delta W_i(v_i). \quad (4.51) $$

Proposition 4.3 together with Eq. (4.49) ensure that the expected portion of the remaining busy period corresponding to future arrivals is nonnegative, i.e. $\Delta W_i(v_i) - B_i(v_i) = A_i(v_i) \ge 0$. Using Eq. (4.50), the expected marginal delay cost that user *i* imposes on users who arrived before her can be determined as:

$$ MB_i(v_i) = \sum_{k=1}^{q_i} (W_i(v_{jk} : z_{ijk} = \sigma_{\bar{v}_i}) - W_i(v_{jk} : z_{ijk} = \sigma_{\underline{v}_i}))v_{jk}. \quad (4.52) $$

To determine the expected marginal delay cost that user *i* imposes on future arrivals, denoted $MA_i(v_i)$, we use the result of Dolan [1978] who observed that this marginal delay cost can be calculated by integrating $\frac{dA_i(v)}{dv}v$ over the bid range $[v, v_i)$. Further, the author observed that the state is constant for any bid comprised between two consecutive bids of users in the pricing queue. Thus, for any segment corresponding to a pair of consecutive bids in the sequence $\{\underline{v}, v_{j_1}, \dots, v_{j_{q_i-1}}, v_{j_{q_i}}, v_i\}$, the expected marginal delay cost imposed by user *i* onto users bidding in this segment can be calculated by integrating $\frac{dA_i(v)}{dv}v$ over the corresponding domain. Specifically, let $v_{j_0} = \underline{v}$ and let $v_{j_{q_i+1}} = v_i$. Let $D(k) = (v_{jk}, v_{jk+1})$ be the domain corresponding to the bid segment $(v_{j_m}, v_{j_n})$ for $m < n$ such that $j_m, j_n \in V_{<i}$. Accordingly, the expected marginal delay cost imposed by user *i* on future arrivals can be determined as:

$$ MA_i(v_i) = \sum_{k=0}^{q_i} \int_{D(k)} \frac{dA_i(v)}{dv} v dv = \sum_{k=0}^{q_i} \int_{D(k)} \frac{-dW_i(v)}{dv} v dv. \quad (4.53) $$

The expected marginal delay cost imposed by user *i* declaring $v_i$ on other users is thus:

$$ MC_i(v_i) = MB_i(v_i) + MA_i(v_i). \quad (4.54) $$

We next give the main result.

**Theorem 4.2.** If the user objective function is $\min_{v_i} v_i^* W_i(v_i) + P_i(v_i)$ and users are serviced in order of decreasing declared delay costs, the payment $P_i(v_i) = MC_i(v_i)$ is incentive-compatible.

Theorem 4.2 proves that the expected marginal delay cost as defined by Eq. (4.54) represents an incentive-compatible user payment in the dynamic sense, i.e. by taking into account the expected impact of users on future arrivals. This provides a basis to implement the proposed queue- and lane-based mechanisms in an online fashion. Upon reaching the front of their lane, users are asked to declare their delay cost $v_i$ and receive a corresponding payment $P_i(v_i)$ which ensures strategyproof user behavior in the long run.

The proposed queue- and lane-based traffic intersection mechanisms are computationally efficient. In terms of computational resources, both mechanisms require the solution of Eqs.
(4.52) and (4.53) to determine user payments. Eq. (4.52) requires the solution of multiple
systems of linear equations ((4.37) or (4.42) depending on the mechanism selected) which
can be accomplished using standard linear algebra codes. To compute Eq. (4.53), standard
numerical integration and differentiation techniques can be used in combination with codes for
systems of linear equations. As will be shown in our numerical experiments, both mechanisms
have low average computation time per user and can thus be expected to be executed in real-
time to determine incentive-compatible payments.

4.5.4 Numerical experiments

We implement the proposed online mechanisms by simulating the arrival and departure of users at an intersection using a discrete time process. At every time period, random trials are conducted for each empty lane to simulate the stochastic arrival process of users. For all our experiments, we use a unit one-step cost, i.e. $g_q(q_i) = g_z(z_i) = 1$ for all states $q_i \in S_i^q$ and $z_i \in S_i^z$. We set the time period duration to 1 s, i.e. representing a uniform user service time of 1 s. This value is chosen based on the geometric configuration of typical 4-lane intersection and average vehicle speeds [Levin and Rey, 2017]. Users' delay costs are randomly drawn from a uniform distribution $U(5, 10)$, i.e. $v = \$5/hour$ and $\bar{v} = \$10/hour$. These values are chosen based on a study on value of time for automated vehicles [Wong et al., 2018]. This delay cost range corresponds to user bids between 0.14 c and 0.28 c for a time period duration of 1 s which aims to represent users' willingness to pay for service priority at traffic intersections.

To quantify the value of using a dynamic auction compared to a static auction, we imple-
ment the static mechanism outlined by Dolan [1978] for priority queueing systems. This static
mechanism is a Vickrey-Clarke-Groves auction which is known to be incentive-compatible in
the static sense, i.e. all participating users are assumed to be known when determining pay-
ments [Vickrey, 1961, Clarke, 1971, Groves and Loeb, 1975]. The static mechanism utilizes
the same priority queueing model than that of the proposed online mechanism, i.e. users in
the pricing queue are sorted by decreasing declared delay costs, but differ in the determina-
tion of waiting time and user payments. In the static mechanism, the waiting time of user
i is only function of the position of user i in the pricing queue, e.g. if user i is third in the
pricing queue, her waiting time is the sum of the two one-step costs for transitioning from
the current state to the state where user i served. The payment of user i is sum of user
valuations over lower-bidders in the pricing queue. This payment is incentive-compatible in
the static sense [Dolan, 1978].

In each experiment, we simulate the service of one million users and report the average user behavior by segmenting users' delay cost. Specifically, we segment the delay cost range [5, 10] into 30 uniform bins and group all serviced users throughout the simulation into these bins. We then report average waiting times, payments and generalized user costs based on the average quantities obtained for each of these 30 bins. The simulations are implemented in Python and the linear systems of equations (4.37) and (4.42) are solved using Numpy's linear algebra solver on a Windows 10 machine with 8 Gb of RAM and a CPU of 2.7 Ghz. We compare the behavior of the proposed traffic intersection auctions mechanisms on a four-lane intersection under non-uniform lane arrival probabilities on a four-lane intersection. We consider non-uniform lane arrival probabilities $p_j$ for all lanes $j$ of a four-lane intersection.

We use cardinal directions, i.e. East (E), South (S), West (W) and North (N) to refer
to all $Q = 4$ lanes of the intersection and assign lane arrival probabilities of 0.50, 0.25,
<figure><img src="image_15.png" /><figcaption>(a) Static mechanism.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Queue-based mechanism.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 4.5: Relative user generalized cost ($C_i(v_i) - C_i(v_i^*)$) / $C_i(v_i^*)$. Results of a simulation with 1 million users with lane-based arrival probabilities of 0.50, 0.25, 0.15 and 0.10. Figure 4.6a illustrates the static mechanism; Figure 4.6b illustrates the queue-based mechanism and Figure 4.5c illustrates the lane-based mechanism.</figcaption></figure>
0.15 and 0.10, respectively. For the queue-based mechanism, we use the average lane arrival probability, i.e. $p = \sum_{j=1}^{Q} p_j = 0.25$. We report lane-based heatmaps of the relative difference in terms of relative user generalized cost $(C_i(v_i) - C_i(v_i^*)) / C_i(v_i^*)$ for the static, queue- and lane-based mechanisms in Figure 4.5. Figure 4.6a shows the outcome of the simulation using the static mechanism. We find that the static mechanism does not promote truthful user behavior for all four lanes of the intersection, i.e. regardless of the lane arrival probability the mechanism is not incentive-compatible, and on all lanes users may reduce their generalized cost in the long run by over-reporting their delay cost. The outcome of the same simulation obtained using the queue-based mechanism is summarized in Figure 4.6b. The outcome of these simulation highlight that using the average lane arrival probability in the queue-based mechanism fails to yields incentive-compatible outcome for the East, West and North lanes which have a lane arrival probability different that $p = 0.25$. Only users on the South lane with $p_S = 0.25 = p$ receive incentive-compatible payments. In turn, Figure 4.5c illustrates that the lane-based mechanism is incentive-compatible for each lane of the intersection.

Additional numerical results can be found in Rey et al. [2021].

## 4.6 Conclusion

This chapter first presented a method for optimizing reservation-based intersection controls in AIM before embedding this conflict point formulation in a stochastic network traffic model. A pressure-based network traffic control policy for coordinating legacy vehicles (LV) and autonomous vehicles (AV) was proposed under the assumption that LVs and AVs share the urban network infrastructure. This policy combines green (for LVs) and blue (for AVs) phases to coordinate traffic at network intersections. We characterized the stability region of the proposed queuing system and showed that the proposed decentralized hybrid network control policy is stable, i.e. that it maximizes network throughput, under conventional travel demand conditions.

The last part of this chapter focus on market-based auction mechanisms for AIM. Assuming that users can declare their delay cost privately to the intersection manager, we investigated the computation of incentive-compatible payments for users to minimize users’ generalized cost which is defined as a linear combination of expected waiting time and user payment. We introduced two Markov chain models to determine users’ expected waiting time, and presented a payment mechanism that can be implemented with both models. We showed that the proposed online traffic intersection mechanisms are incentive-compatible in the dynamic sense and thus maximize social welfare.
# Chapter 5

## Perspectives

The goal of this thesis was to provide an overview of some contributions in optimization and game theory, and their application to decision-making problems encountered in transportation systems. These contributions have been organized into three main chapters that differ by the operations research (OR) methodologies used therein, as well as by the decision-making problems selected in each chapter. Chapter 2 focused on contributions in the field of bilevel optimization and its applications to transportation network design problems. Chapter 3 turned to mixed-integer nonlinear programming and its applications to aircraft conflict resolution problems. Chapter 4 embraced a broader modeling perspective by summarizing different optimization and game-theoretical approaches to autonomous intersection management. In each of these three chapters, novel approaches that leveraged state-of-the-art advances in OR techniques were conceived, developed and implemented. The contributions of this thesis are thus both theoretical and numerical, and they lie both in OR as well as in the field of transportation systems science.

Next, an attempt is made to identify research perspectives that could build on these works as well as to provide guidelines for extending the OR methodologies used therein. The rest of this section is organized in two parts: methodological perspectives are examined in Section 5.1, and perspectives in transportation are discussed in Section 5.2.

### 5.1 Methodological perspectives

This section attempts to synthesize the state-of-the-art in the fields of bilevel optimization (Section 5.1.1) and mixed-integer programming (Section 5.1.2) and discusses how methodological developments could be steered to achieve greater efficiency. This section also extends the discussion to machine learning and data-driven optimization approaches (Section 5.1.3), and their potential to boost optimization and game-theoretical methods.

#### 5.1.1 Bilevel optimization

Bilevel optimization is a rapidly growing field of research. Chapter 2 discussed its applications to transportation network design problems, which provide a natural framework to develop bilevel optimization approaches given the existing hierarchical interactions among multiple players such as transportation authorities and users. Of the many challenges faced
by researchers in this field, one of most stringent is the computational scalability and the quality of solution methods for bilevel optimization problems. It can roughly be claimed that bilevel optimization problems whose follower problem(s) enjoy convexity have a substantial advantage over problems that do not enjoy this property. This is largely due to the fact that bilevel optimization problems with convex follower problem can be reformulated into single-level optimization problems by making use of the follower's first-order optimality conditions, such as Karush Kuhn Tucker (KKT) conditions or strong duality. It should be noted that the single-level reformulation of bilevel optimization problems is not trivial and, often, not equivalent. It is non-trivial because after embedding the follower's optimality conditions in a single-level formulation, the latter is a nonconvex optimization problem, e.g. a mathematical program with complementarity constraints (MPCC) that requires further care to be solved efficiently to optimality. In general, the single-level MPCC reformulation is not equivalent to the bilevel optimization problem in the sense that their solution sets may not always coincide. However, it is well-known that, under certain constraint qualification conditions, an optimal (global) solution of the single-level MPCC reformulation corresponds to an optimal (global) solution of the bilevel optimization problem [Dempe and Dutta, 2012]. Nevertheless, bilevel optimization problems with a convex follower problem are globally much more tractable than their counterparts. From a practical standpoint, however, this limitation may tend to incentivize modelers to simplify their representation of the follower problem so as to achieve convexity. This simplification is not always desirable as it prevents, notably, the treatment of bilevel optimization problems with integer follower problems even though substantial progress has been achieved for this class of bilevel optimization problems [Fischetti et al., 2017, Liu et al., 2021]. Nevertheless, these approaches remain quite challenging to implement and may require certain assumptions on the structure of the bilevel optimization problems that can be handled. Even in the case of linear bilevel optimization, computational tractability and scalability remains a concern [Kleinert et al., 2020, 2021]. Further, by large, most efforts in the community have been dedicated to optimistic bilevel optimization problems and few methods are available for solving pessimistic bilevel optimization problems [Aussel and Svensson, 2019, Zeng, 2020].

Transportation systems often aim to improve social welfare using specific performance metrics, e.g. congestion, pollution. Hence, such systems shall be designed to anticipate the reaction of the users of the system whose goals are unlikely to be aligned with system-level objectives. This was the context of the transportation network design problems discussed in Chapter 2 wherein users' collective reaction is modeled as a convex optimization problem that characterizes a Nash equilibrium. Yet, as outlined therein, solving large scale transportation network design problems to optimality is extremely challenging. This is partly due to the fact that traditional branch-and-bound algorithms that branch on complementary slackness conditions cannot scale easily due to the large number of nodes and links—let alone paths—in large-scale transportation networks. There is thus a need to develop more scalable exact algorithms that can handle the large dimensions of realistic city-scale transportation networks. A possible avenue of research is to explore the potential of column generation procedures to iteratively obtain new path variables and embed this approach in a traditional complementary slackness based branch-and-bound algorithm. One caveat is that iteratively generating columns may imply that some follower constraints must also be generated iteratively, thus leading to a row-and-column generation framework wherein newly generated columns also create new follower constraints for which complementary slackness conditions must hold.

Efforts are also needed to go beyond the classical bilevel network design problem under Nash equilibrium. While this model has been extensively studied, it relies on several
fundamental modeling assumptions such as link travel cost convexity. The incorporation of multiple players, as opposed to a system-level approach also holds great potential for development. This is reinforced by emerging trends in on-demand mobility and logistics where several service providers compete for transportation demand [Perboli et al., 2021]. Bilevel optimization offers the possibility to consider games with multiple leaders and/or followers and provide new insights for market design and pricing in transportation systems.

A complementary avenue of research is to attempt to refine the leader's model of the follower by incorporating uncertainty. Stochastic optimization—and stochastic programming especially—have greatly contributed to extending deterministic formulations to account for the uncertainty of model parameters. Much fewer efforts have been devoted to the study of stochastic bilevel optimization problems [Patriksson, 2008, Fampa et al., 2008]. A colleague once told me: “solving a bilevel optimization problem is working for the leader”. This is a general statement about bilevel optimization that recognizes the fact that the follower problem is only the leader's representation of the follower's best response. It is thus natural to question whether the leader has perfect knowledge of the follower's strategy. In the event that the answer to this question is not yes, stochastic bilevel optimization could provide elements of answer to design richer formulations that account for the leader's uncertainty regarding the follower's response. For example, a transportation network user may adopt a mixed strategy based on travel time probability distributions to choose her route. Assuming that the leader has knowledge of such distributions, it could therefore be appropriate for the leader to anticipate this mixed strategy. This highlights that stochastic bilevel optimization holds substantial potential to enrich existing deterministic approaches by introducing uncertainty at either or both levels.

### 5.1.2 Mixed-integer programming and solvers

Mixed-integer programming (MIP) has become a standard approach to solve challenging decision-making problems. Recent advances in off-the-shelf MIP solver capabilities and competition among MIP solvers has lead to substantial computational improvements over the past decades [Achterberg and Wunderling, 2013]. This means that the size of MIP problems that can be solved in reasonable time (and memory) today is much greater than that it was in the 2000s [Lodi, 2010]. This trend has supported the democratization of MIP solvers in the OR community and beyond. This has led to a paradigm shift: a mixed-integer linear programming (MILP) problem can now be viewed as tractable, even though no polynomial-time exact algorithm is available to solve it. Not all MILPs enjoy this treatment. There remains families of MILPs, such as vehicle routing problems, that are much too challenging to solve optimally using off-the-shelf solvers. However, some MILPs have become relatively easy to solve using off-the-shelf MIP solvers. One of the immediate consequences is that it has now become reasonable to devise solution methods wherein a MILP is solved repeatedly under the assumption that a (near-)optimal solution can be found rather quickly. This has also boosted the development of matheuristics wherein mathematical programming problems are solved in the course of a broader heuristic algorithm [Fischetti and Fischetti, 2018]. Other classes of MIP problems have also benefited from this progress, e.g. mixed-integer quadratic programming (MIQP) and mixed-integer quadratically-constrained programming (MIQCP). Chapter 3 illustrates this by showing that a challenging nonconvex mixed-integer nonlinear programming (MINLP) problem can be solved to optimality in reasonable time by an algorithm which repeatedly solves MIQPs and MIQCPs. This suggests that the use of certain MIP formulations within broader algorithmic frameworks holds great potential.
Transportation systems are rich in decision-making problems, many of which can be viewed as variants of classical OR problems. There are therefore strong incentives to develop the practice of using MIP in solution approaches to decision-making problems encountered in transportation systems, even when real-time capabilities are at stake. Chapter 4 advocates that vehicle scheduling within an autonomous intersection management context can be solved efficiently using MILP up to certain levels of demand. This off-the-shelf approach fits well operational problems that can be solved repeatedly using a rolling horizon framework. For instance, in air traffic control, it is difficult to solve conflict resolution problems with a large lookahead since there is often too much uncertainty associated to aircraft's future positions beyond a 30-minute horizon. Hence, from a practical standpoint, it is meaningful to adopt a rolling horizon framework with a balanced lookahead and to repeatedly solve a small-to-medium scale MILP that can provide competitive solutions. The main drawback of this approach is that this may be too myopic in that future events are ignored even though some knowledge of possible scenarios may be available. There are thus opportunities to enrich MIP-based approaches by incorporating discounted lookahead information or even attempting to learn data features over time. This could be achieved by using surrogate objective functions that aim to anticipate the impact of future events.

MIP provides a flexible modeling framework for OR practitioners. Cutting-edge solution methods for challenging OR problems often go beyond the implementation of direct, off-the-shelf MIP solvers. Modern MIP solvers now offer additional features such as the automatic and customizable application of Benders' decomposition to MILPs [Bonami et al., 2020]. MIP modeling environments also provide new tools to implement branch-and-cut-and-price algorithms wherein column generation tasks are made easier to handle for the non-expert, even though the learning of such modeling environments itself may present some challenges [Maher et al., 2017]. This trend supports the development of MIP and also advocates for the use of MIP solvers as a default go-to modeling and solution approach. One notable side-effect is the computational opportunities that arise from the ubiquity of MIP solvers: under the hypothesis that their computational performance will continue to improve, there may be a time when formulating and solving a MIP outperforms dedicated combinatorial algorithms. An immediate benefit is the facilitation of complex solution methodologies such as branch-and-cut-and-price algorithms. In the latter, a pricing algorithm is required to price columns as part of the column generation process. In the context of vehicle routing problems, the pricing subproblem is often a constrained shortest path problem that is NP-hard. Historically, such subproblems have been solved extensively with customized label correcting algorithms that may be difficult to implement efficiently. The uprising of MIP solvers raises the question: *When, if ever, will MIP solvers become competitive with label-correcting algorithms for pricing subproblems that arise when solving vehicle routing problems via branch-and-price algorithms?* This also applies to the context of bilevel optimization where MIP solvers may be used to solve follower problems (e.g. Nash equilibrium problems) faster than dedicated algorithms, thus potentially helping in rapidly obtaining competitive bounds. More broadly, this perspective may someday shift the practice of algorithm benchmarking in certain families of OR problems and this may provide opportunities to develop innovative, MIP-powered solution approaches.

### 5.1.3 Machine learning

Machine learning (ML) has been quite successful in solving challenging optimization problems [Bengio et al., 2021]. Most ML-powered approaches to optimization have been targeted
towards large-scale problem instances and are therefore heuristic by design. More recently, the OR and ML communities started to explore the potential of ML-powered approaches to enhance mathematical programming-based solution algorithms for challenging optimization problems. This may take the form of data-driven heuristics that are triggered within the course of branch-and-cut-and-bound algorithm in an attempt to obtain competitive bounds to prune the tree. One drawback lies in the vast amount of computational resources required for efficient model training in order to obtain a well-trained predictive model. To overcome this challenge, carefully selected samples should be obtained to balance training time and the quality of the predictive model. The convex approximation of congestion effects developed in Chapter 2 to obviate the need to solve repeated follower problems can be viewed as an attempt to integrate predictive modeling—in this case nonlinear regression—within a branch-and-price framework. Other applications of ML are possible: for example, it may be used to learn the likelihood that a variable be part of the optimal solution [Tahir et al., 2021]. Another instance is the use of ML to learn optimal branching rules [Lodi and Zarpellon, 2017]. A compromise between avoiding resource-intensive model training and reaping the benefits of ML-powered approaches may be found in the field of reinforcement learning (RL) which is closely related to Markov decision processes (MDP). RL consists of a training-free approach to ML: in RL, optimal decision rules are learnt *online*, i.e. over repeated learning experiences. Hence, casting a branch-and-bound tree search as a MDP opens the door to a suite of RL-based algorithms that could help in identifying optimal tree search decisions such as branching, node selection rules or the decision of when triggering certain heuristics.

So far the interplay between ML and OR has been examined from the perspective of ML-
powered approaches to enhance solution algorithms for OR problems. The reciprocal, i.e. the
use of OR techniques to enhance solution approaches for ML, may be more challenging since
ML usually targets large scale problem instances. It is unlikely that tools from mathemati-
cal programming or combinatorial optimization could enhance the solving of large-scale ML
problems from a computational efficiency standpoint. Nevertheless, recent works on feature
selection for predictive modeling, a classical problem in the ML community, suggest that
bilevel optimization may bring novel insights in regards to solution quality. Feature selection
can be cast as a bilevel optimization problem where the leader represents the modeler and the
follower represents the training problem [Agor and Özaltın, 2019]. If the target ML model
admits a convex optimization formulation, such as support vector machines, one may refor-
mulate the latter using its first-optimality conditions to obtain a single-level mathematical
program with complementary constraints that is amenable to classical solution approaches of
the bilevel optimization literature. While scaling up may remain a challenge, such OR-driven
approaches may assist in discovering global optimal solutions of certain ML problems. Hence,
this stream of research provides an exciting sandbox to test the potential of OR techniques
to assist solving ML problems.

## 5.2 Perspectives in transportation systems

This section adopts a different stance and attempts to discuss perspectives and trends in
transportations systems. The realm of stochastic, dynamic and online transportation prob-
lems is first discussed (Section 5.2.1). Future transportation systems are then considered
both from the perspective of technological advances and business models (Section 5.2.2).
This section concludes by discussing pressing transportation challenges that are inherently
tied to sustainability (Section 5.2.3).
5.2.1 Stochastic, dynamic and online transportation problems

Transportation systems have evolved considerably over the past few years. The advent of mo-
bile communication technologies and the scaling up of computing capabilities has opened the
door to new frameworks for on-demand mobility and logistics. This has been first leveraged
by transportation service providers, e.g. Uber, Lift, Didi; who have developed competitive
solutions for urban mobility. Over time, public transportation authorities have also embraced
this new paradigm, and several cities worldwide have welcome new mobility solutions that
attempt to complement classical services. This shift from traditional fixed schedule mobil-
ity systems towards demand-responsive transportation is reflected in the scientific literature
where a fast-growing number of studies are devoted to the latter [Vazifeh et al., 2018].

On-demand transportation presents specific challenges: anticipating future demand and adjusting transportation solutions over time call for stochastic and dynamic solution approaches [Soeffker et al., 2021]. Often, dynamic transportation models aim to capture either time-varying traffic conditions or time-dependent pricing policies. It should be emphasized that this is different to multi-period formulations: in the latter, the time scale is larger and at each time period a “copy” of the problem is solved—this is the case of the multi-period DNDP studied in Chapter 2. In contrast, dynamic routing problems target a much finer time resolution, e.g. time-dependent shortest path problems aim to find least-cost paths in time-dependent graphs where link weights vary based on the arrival time at the head of the link. There exists also a distinction between dynamic optimization problems where the data is time-dependent and an a priori solution is sought, and online optimization problems where new data is available over time. In the latter, the solution of the problem is a policy, that is, a decision rule. Online optimization problems are highly relevant to modern transportation given the emergence of on-demand mobility and logistics systems where user requests are not known in advance, such as in same-day delivery systems [Voccia et al., 2019]. In several situations, the modeler may have access to historical data that allows her to obtain probability distributions of, for example, user demands. In this context, the decision-making problem can be cast as a stochastic programming problem where an a priori solution is sought by optimizing the expected value—or other non-risk-neutral metrics—of the objective function under such probability distribution. If the problem is both stochastic and online, i.e. new data is obtained on-the-fly, parameters’ probability distributions can be learnt or refined over time. This provides room to explore related solution concepts such as RL, which is widely popular in the artificial intelligence community yet remains in its infancy among OR and transportation practitioners [Powell, 2019]. RL, much like approximate dynamic programming, can be viewed as a heuristic to solve large scale MDPs which in turn can be used to model several stochastic and online decision-making problems encountered in transportation. There exists thus ample opportunities to develop innovative solution approaches to modern transportation problems, notably in the realm of on-demand mobility and logistics.

5.2.2 Future transportation systems

Recent technological developments in automated and autonomous transportation systems
bring interesting perspectives. While, to date, autonomous cars are not pattroling our streets,
unmanned aerial vehicles—or drones—are increasingly being tracked for their potential to
change urban and rural logistics. In the longer term, urban air mobility may emerge as a
novel mode of transportation [Al Haddad et al., 2020]. The development of autonomous trans-
portation systems offers a rich playground for the OR community. An example is provided in
Chapter 4 which explores how autonomous intersection management could be designed and implemented on top of—or alongside—the existing transportation infrastructure. Unmanned air mobility and logistics brings exciting decision-making problems to the table such as finding optimal location of vertiports [Shin et al., 2022] or the operations of urban airways to route drone and air taxi traffic [Rajendran and Srinivas, 2020]. Future air traffic management systems are also expected to benefit from automation. This is largely discussed in Chapter 3 at the tactical decision-making level in the context of air traffic control. At a broader scale, automation holds great potential to sustainably improve air traffic flow management [Bertsimas et al., 2011]. Several recent efforts have also investigated the potential of market mechanisms to account for users’ preferences in the decision process [Castelli et al., 2011, Dal Sasso et al., 2019]. On the road, autonomous mobility provides numerous opportunities to rethink how we move. Shared autonomous vehicles are increasingly being considered as opposed to private autonomous vehicle ownership. In a shared autonomous vehicle context, these vehicles are as a service and thus are not owned by consumers [Krueger et al., 2016]. Automated mobility also has implications regarding users’ value of time which could be used to dampen congestion effects if travelers are busy working in their vehicles during their commute [Monteil, 2014, Rashidi et al., 2020].

The future of transportation is steered towards new business models. As the num-
ber of mobility or logistics service providers increase, transportation systems are becom-
ing marketplaces where transportation may be priced “sur mesure”. Customized user pric-
ing is at the heart of the mobility-as-a-service (MaaS) and the logistics-as-a-service (LaaS)
paradigms [Polydoropoulou et al., 2020]. Both MaaS and LaaS are receiving an increasing
attention from both practitioners and researchers due to their potential to shift the way
people and goods are moved thanks to path-based pricing and integrated payments. The
endgoal of MaaS and LaaS is to provide users the possibility of making seamless payments
for their transportation needs. This would allow users to make multimodal trips, possibly
via competing service providers, by making a single payment that captures precisely users’
willingness-to-pay for such services. Pricing in MaaS and LaaS systems may require the
coordination of multiple non-cooperative service providers as well as consumers competing
for resources. This has motivated researchers to study transportation systems through the
lens of two-sided markets wherein buyers and sellers interact through a virtual marketplace.
In recent efforts, we have investigated the potential of single-leader multi-follower games for
two-sided platform pricing in MaaS systems [Xi et al., 2022], and other efforts have been
dedicated to study collaboration across MaaS platforms [Pantelidis et al., 2020]. This stream
of research highlights the pivotal role of business models in shaping future MaaS and LaaS
transportation systems.

It is worth emphasizing that, business implications aside, future transportation systems that rely on path-based pricing face considerable computational challenges. Path-based pricing is often synonymous with non-additive pricing [Han and Lo, 2004]. An example of non-additive, path-based pricing consists of determining optimal user subsidies to encourage the use of less-congested or low-emission routes in a congested road network. Analogously, one could seek optimal user subsidies in a multimodal mobility network where the links of the network are operated by multiple, non-cooperative service providers. These concepts are promising yet they remain computationally challenging to implement. This steams from the inherent complexity of path finding with non-additive costs in medium to large transportation networks. Therefore, innovative methodologies are required to develop non-additive path-based pricing models at scale.
### 5.2.3 Sustainable transportation

Climate change, ongoing conflicts and, more recently, the covid-19 pandemic have also contributed to rethinking transportation systems. Remote working initiatives and e-commerce deliveries have become part of our daily routine. While the “pandemic effect” could be ephemeral, it has triggered several research efforts toward studying and developing more resilient transportation systems capable of absorbing non-recurrent as well as recurrent perturbations such as pollution peaks and adverse weather effects [Budd and Ison, 2020]. The electrification of transportation systems also poses a greater threat to energy networks. Several countries have adopted policies that strongly encourage the adoption of electric vehicles and/or aim to forbid the use of fuel vehicles in a not-so-far future. One may thus expect that cities with a high rate of electrified transportation fleets will face considerable challenges in the context of energy demand surges [Muratori, 2018]. This motivates research efforts focused on developing integrated modeling frameworks that are able to capture the interplay between transportation and energy networks.

Advances in optimization and transportation systems modeling are also expected to lift the potential of humanitarian logistics. Transportation is often central to crisis-induced logistical challenges such as post-disaster network recovery, as outlined in Chapter 2. This also concerns more operational decision-making problems such as humanitarian aid and food rescue logistics. Such ressource allocation and/or routing problems share similarities with their profit-driven counterparts but may differ in the choice of the objective function(s) or in the properties of the solutions sought wherein fairness and equity play a central role [Eisenhandler and Tzur, 2019]. At the same time, climate change puts sustainability at the heart of existing and future transportation systems [Zhao et al., 2020]. Although sustainability-driven transportation systems have existed for several years in the literature, there is now more than ever a pressing need to engineer these solutions and develop incentives for practitioners to implement them. From a modeling standpoint, this provides a narrative to support the development of incentive mechanisms for promoting sustainable transportation such as active transportation modes [Kaspi et al., 2014]. In this regard, game theory and mechanism design offer several methodological solutions to develop next generation transportation systems. Pricing problems can directly benefit from game-theoretical enhancements [Brotcorne et al., 2008]. Pricing is central to several transportation systems such as fare structure in public transportation, travel demand balancing, but also last-mile delivery systems. Designing sustainability-driven incentive mechanisms by adjusting service prices is a promising research direction. One of the main challenge is to develop both implementable and scalable pricing mechanisms. Eliciting user preferences requires incentive-compatible mechanisms. Yet, standard truthful mechanisms such as Vickrey-Clarke-Groves auctions may require private data that users themselves may not be accurately aware of. For example, one may not know her value of time with enough precision. Another limitation of such mechanisms arise in the context of combinatorial auctions: it might not be practical for users to report their preferences for all combinations of a mobility bundle since the number of such combinations may be very large. There is thus a need to develop incentive mechanisms that are practical, and that also target goals that are beneficial for the greater good.
Bibliography

Tobias Achterberg and Roland Wunderling. Mixed integer programming: Analyzing 12 years of progress. In *Facets of combinatorial optimization*, pages 449–481. Springer, 2013.

Joseph Agor and Osman Y Özaltın. Feature selection for classification models via bilevel optimization. *Computers & Operations Research*, 106:156–168, 2019.

Christelle Al Haddad, Emmanouil Chaniotakis, Anna Straubinger, Kay Plötner, and Constantinos Antoniou. Factors affecting the adoption and use of urban air mobility. *Transportation research part A: policy and practice*, 132:696–712, 2020.

Antonio Alonso-Ayuso, Laureano F Escudero, and F Javier Martín-Campo. Collision avoidance in air traffic management: A mixed-integer linear optimization approach. *IEEE Transactions on Intelligent Transportation Systems*, 12(1):47–57, 2011.

Antonio Alonso-Ayuso, Laureano F Escudero, and F Javier Martín-Campo. Exact and approximate solving of the aircraft collision resolution problem via turn changes. *Transportation Science*, 50(1):263–274, 2014.

Antonio Alonso-Ayuso, Laureano F Escudero, and F Javier Martín-Campo. An exact multi-objective mixed integer nonlinear optimization approach for aircraft conflict resolution. *TOP*, 24(2):381–408, 2016.

Florent Altché and Arnaud De La Fortelle. Analysis of optimal solutions to robot coordination problems to improve autonomous intersection management policies. In *Intelligent Vehicles Symposium (IV), 2016 IEEE*, pages 86–91. IEEE, 2016.

Didier Aussel and Anton Svensson. Is pessimistic bilevel programming a special case of a mathematical program with complementarity constraints? *Journal of Optimization Theory and Applications*, 181(2):504–520, 2019.

Saeed Asadi Bagloee, Majid Sarvi, and Michael Patriksson. A hybrid branch-and-bound and Benders decomposition algorithm for the network design problem. *Computer-Aided Civil and Infrastructure Engineering*, 32(4):319–343, 2017.

Hillel Bar-Gera. Traffic assignment by paired alternative segments. *Transportation Research Part B: Methodological*, 44(8):1022–1046, 2010.

Mikel Barbara, David Rey, and Ali Akbarnezhad. Optimizing location of new public schools in town planning considering supply and demand. *Journal of Urban Planning and Development*, 147(4):04021057, 2021.

Jonathan F Bard. *Practical bilevel optimization: algorithms and applications*, volume 30. Springer Science & Business Media, 2013.