# Optimization and game-theoretical methods for transportation systems

David Rey

► To cite this version:

David Rey. Optimization and game-theoretical methods for transportation systems. Operations Research [math.OC]. Toulouse 3 Paul Sabatier, 2023. tel-04024464

HAL Id: tel-04024464

https://hal.science/tel-04024464

Submitted on 11 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Habilitation à Diriger les Recherches

Université Toulouse III Paul Sabatier

École Doctorale Mathématiques, Informatique, Télécommunications
de Toulouse

## Optimization and game-theoretical methods for
transportation systems

Presentée par

DAVID REY

SKEMA Business School, Sophia Antipolis

Spécialité : RECHERCHE OPÉRATIONNELLE

Soutenue à l'École Nationale de l'Aviation Civile (ENAC) le 6 mars 2023
davant le jury composé de :

LUCE BROTCORNE

GUGLIELMO LULLI

MICHAL TZUR

DIDIER AUSSEL

MARCEL MONGEAU

JEAN-BAPTISTE HIRIART-URRUTY

SONIA CAFIERI

Directrice de Recherche, INRIA Lille

Professeur, Lancaster University

Professeure, Tel Aviv University

Professeur, Université de Perpignan

Professeur, ENAC

Professeur Émérite, Université Toulouse III

Professeure, ENAC

Rapporteure

Rapporteur

Rapporteure

Examenateur

Examenateur

Président

Marraine
Acknowledgements

I would like to warmly thank Luce Brotcorne, Guglielmo Lulli and Michal Tzur for serving as referees for my Habilitation à Diriger les Recherches (HDR). I also thank Didier Aussel, Jean-Baptiste Hiriart-Urruty and Marcel Mongeau for participating in my jury. I am extremely thankful to Sonia Cafieri for her continuous support and dedication throughout this journey.

This manuscript is the fruit of several years of exciting research with my Ph.D. students and my peers. I hope I have not forgetten anyone. First, I would like to thank my former and current Ph.D. students at UNSW Sydney for their efforts and commitment: Mikel Barbara, Shantanu Chakraborty, Fernando H. C. Dias, Ahmed. W. Hammad, Xiaotong (Sharon) Dong, Haoning Xi, Xuefen Liu, Mingye Luan, Zahra Nourmohammadi, Aditya Paul, Esta Qiu, Babak Shahbazi and Xiang Zhang.

I was extremely lucky to spend several years at the Research Centre of Integrated Transport Innovation (rCITI) at UNSW Sydney. This is an outstanding research team where I met many friends and colleagues that have contributed to shaping my academic trajectory. Thank you Sylvia Brohl, Vinayak V. Dixit, Lauren M. Gardner, Hanna Gryzbowska, Sisi Jian, Maria Lee, Wei Liu, Mojtaba Maghrebi, Emily Moylan, Divya J. Nair, Ali Najmi, Ali Nezhad, Taha H. Rashidi, Meead Saberi, Changle Song, Upali Vandebona, S. Travis Waller, Kasun P. Wijayaratna, Zhitao Xiong and Aleksa Zlojutro.

I would also like to thank my colleagues and co-authors: Khaled Almi'ani, Hillel Bar-Gera, Romain Billot, Nour-Eddin El Faouzi, Angelo Furno, Elise Henry, Hassan Hijazi, Michael Levin, Julien Monteil and Jean-Luc Ygnace.

I am grateful for the welcome and the support of SKEMA Business School. Thank you Jan Broekaert, Xavier Brusset, Morteza Davari, Faizal Hafiz, Aïda Jebali, Roland Königsgruber, Xuefei Lu, Philippe Monin, Davide La Torre and Bernardo K. Pagnoncelli.

Last, I would like to thank Rachel, Elise, Helene, Caroline, Lucas and my parents Claire and Philippe.
## Foreword

This *Habilitation à Diriger les Recherches (HDR)* thesis aims to provide an overview my research activities since I defended my Ph.D. in December 2012. I started working in the field of Operations Research (OR) during my Ph.D. which focused on mixed-integer programming methods for air traffic management. After this first research experience, I explored new applications of OR mostly in the field of, but not restricted to, transportation and logistics. This includes notably network design and management, facility location problems, shared and automated mobility, food rescue logistics, incentive mechanisms and pricing.

While diverse in application topics, these studies share common features: they are largely based on optimization and game-theoretical methods ranging from mixed-integer programming to bilevel optimization and mechanism design. This thesis attempts to show selected samples of my contributions in developing such methodologies and their application to targeted transportation systems.
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Bilevel optimization problems in transportation networks</td><td>5</td></tr><tr><td>2.1</td><td>Introduction</td><td>5</td></tr><tr><td>2.2</td><td>Literature review</td><td>6</td></tr><tr><td>2.3</td><td>Computational benchmarking of exact methods for the DNDP</td><td>7</td></tr><tr><td>2.4</td><td>Maintenance scheduling in transportation networks</td><td>13</td></tr><tr><td>2.5</td><td>Scheduling for road network disaster recovery</td><td>29</td></tr><tr><td>2.6</td><td>Conclusion</td><td>35</td></tr><tr><td>3</td><td>Aircraft conflict resolution via mixed-integer nonlinear programming</td><td>37</td></tr><tr><td>3.1</td><td>Introduction</td><td>37</td></tr><tr><td>3.2</td><td>Literature review</td><td>38</td></tr><tr><td>3.3</td><td>Aircraft conflict resolution formulations</td><td>39</td></tr><tr><td>3.4</td><td>Exact solution methods for the aircraft conflict resolution problem</td><td>50</td></tr><tr><td>3.5</td><td>Robust aircraft conflict resolution problem</td><td>56</td></tr><tr><td>3.6</td><td>Numerical experiments</td><td>61</td></tr><tr><td>3.7</td><td>Conclusion</td><td>71</td></tr><tr><td>4</td><td>Modeling approaches for autonomous intersection management</td><td>73</td></tr><tr><td>4.1</td><td>Introduction</td><td>73</td></tr><tr><td>4.2</td><td>Literature review</td><td>74</td></tr><tr><td>4.3</td><td>Conflict point formulation</td><td>77</td></tr><tr><td>4.4</td><td>Hybrid network traffic control policy</td><td>83</td></tr><tr><td>4.5</td><td>Online mechanism design for traffic intersection auctions</td><td>92</td></tr><tr><td>4.6</td><td>Conclusion</td><td>104</td></tr><tr><td>5</td><td>Perspectives</td><td>105</td></tr><tr><td>5.1</td><td>Methodological perspectives</td><td>105</td></tr><tr><td>5.2</td><td>Perspectives in transportation systems</td><td>109</td></tr></table>
# Chapter 1

# Introduction

## Motivation

Transportation systems can be viewed as the blood vessels of a society on the grounds that they provide essential routes for the movements of vital materials that keep every organ alive. They enable society to function by facilitating human interaction and developing economies. Transportation occurs at various levels of human society. It can connect people, provide access to facilities, or ship packages to populations in need. Yet, transportation systems almost always lead to negative externalities. With few exceptions, the act of moving is usually not desirable in itself and it is thus natural to optimize transportation with the goal to minimize such negative externalities. The essence of transportation systems is thus to optimize the movement of such commodities, e.g. people across places or goods between locations. The demand for commodities is linked to human behavior, and it is natural to represent transportation users as utility-maximizing agents. This sets the scene for game-theoretical decision-making frameworks where agents interact with each other through a common medium, i.e. a transportation system.

Several areas of research directly or indirectly connect with transportation systems. For instance, civil infrastructure engineering, urban mobility, logistics, supply chain management may all consider—and even aim to improve—a part of a transportation system relevant to their domain of study. Hence, the optimization of transportation systems encompasses a multitude of efforts across various disciplines including mathematics, economics, computer science, ecology, psychology and humanities. From a methodological standpoint, this thesis does little to address the multiple facets of transportation. The scope of this thesis is largely restricted to operations research (OR) methods for transportation systems with an emphasis on mathematical programming—and more specifically mixed-integer programming—based approaches.

OR has a long history of applications in transportation. Several seminal optimization problems in OR are related to transportation. The assignment and its extension commonly known as the *transportation problem* are typical illustrations of decision problems where a solution is said to be optimal if it minimizes the total cost of connecting two groups of agents. The traveling salesman problem is concerned with the minimization of distance traveled. Facility location aims to identify optimal placements of facilities based on their distance to demand nodes. Network design seeks the optimal configuration of a network so as to minimize the total usage costs—often travel related—induced by commodities. Trajectory optimization
deals with the movements of objects through time and space.

In the OR literature, the traditional treatment of the aforementioned decision problems involves a mix of graph theory, combinatorial optimization, algorithm design and mathematical programming; amongst others. The modeling and the analysis of transportation systems often goes beyond this practice by also embracing a more user-oriented perspective that aims to take into consideration user behavior in the decision-making process. For example, the main difference between a classical network design problem and a transportation network design problem is that the former is purely concerned with cost minimization subject to design constraints while the latter incorporates additional network equilibrium elements to capture the long-term behavior of utility-maximizing—or selfish—users. The incorporation of user behavior in classical OR problems can also be viewed as the problem of finding the optimal incentives to achieve a certain pattern which falls in the realm of mechanism design.

This brief analysis highlights the strong links existing between transportation systems and the fields of optimization and game theory. These links stand at the heart of this thesis which contributes to the development of solutions using OR techniques for decision-making problems that arise in transportation systems. It should be highlighted that while this thesis is not dedicated to the study of classical OR problems, several of the problems discussed therein can be viewed as variants of these seminal decision problems.

## Organization of the thesis

This thesis is organized around three main chapters. Each of these chapters addresses a specific decision problem encountered in transportation systems. The chapters are organized in terms of methodological approaches as well as by the time scale of the decision problems addressed. Chapter 2 summarizes some contributions in the field of bilevel optimization and its applications to transportation network design and maintenance scheduling problems. Such decision-making problems can be viewed as *strategic* in that the decisions involved are based on the long-term behavior of transportation users, which represents a network equilibrium. Chapter 3 focuses on mixed-integer nonlinear programming (MINLP) and its applications to aircraft conflict resolution. This is a *tactical* decision-making problem encountered in air traffic control that seeks short-to-medium term actions to avoid collisions and ensure the safe operations of aircraft in congested airspace. Chapter 4 focuses on modeling approaches for autonomous intersection management in road traffic systems. Autonomous intersection management can be viewed as an *operational* decision-making problem that aims to develop solution approaches for enabling connected and automated mobility at traffic intersections. Various methodologies ranging from mixed-integer linear programming (MILP), decentralized control and mechanism design are explored to enable this future mobility paradigm. An overview of each of these three chapters is provided below, along with some context.

Chapter 2 is concerned with network design problems in transportation. Here the time scale is strategic, that is, long-term decisions are sought. This may include, for example, network capacity addition, link maintenance or repair aspect, or even facility location decisions. In constrat to the classical literature on network design problems, however, network design problems in transportation aim to anticipate the reaction of network users which is typically modeled as a Nash equilibrium problem. This sets the frame for a so-called leader-follower game, also known as a Stackelberg game wherein the leader represents the network planner and the follower represents the collective route choice of network users under some equilibrium conditions. It is critical to emphasize that the strategic nature of the problem is fundamental
since user equilibrium conditions aim to represent users' learning of network conditions over
repeated experiences. This chapter focuses on the analysis and on the development of bilevel
optimization methods for transportation network design problems with a focus on discrete
problems wherein design decisions are represented by integer variables. An initial analysis
of the state-of-the art methods for solving bilevel optimization problems in transportation
networks is conducted. This is followed by the development of a specific method for a multi-
period network maintenance scheduling problem based on a branch-and-price approach. An
application to network disaster recovery is also briefly presented.

Chapter 3 addresses a more tactical decision-making problem in the field of air traffic control. It focuses on the aircraft conflict resolution problem, which was also the topic of my Ph.D., but examines an alternative modeling and optimization framework. Specifically, this chapter summarizes recent efforts to build on and extend state-of-the-art formulations and exact solution methods for the aircraft conflict resolution problem. Given a set of aircraft trajectories that may lead to conflicts (loss of separation) or even collisions, the aircraft conflict resolution problem aims to find a set of conflict-free aircraft trajectories that ensure that minimal separation requirements are maintained at all times. The difficulty of this problems largely depends on its decision space and on the modeling assumptions made. This chapter tackles the continuous form of the problem where neither time nor control actions are discretized. This variant of the aircraft conflict resolution problem is particularly challenging since controlling aircraft headings involves manipulating trigonometric functions that are nonlinear and also involves Euclidean distance-based separation constraints that are nonconvex. Throughout this chapter, a new formulation based on a reformulation technique involving complex numbers is proposed and shown to be empirically competitive compared to state-of-the art algorithms. An extension which aims to incorporate the impact of uncertainty on aircraft predicted trajectories in the conflict resolution problem is also presented. Initial results demonstrate the potential of robust optimization to handle such uncertainty.

Chapter 4 targets an operational problem in the context of autonomous intersection management. Traffic intersections are major bottlenecks of urban transportation networks. Autonomous intersection management is an emerging urban mobility paradigm under which traffic intersections behave autonomously in a signal-free fashion. In this context, users are assumed to be able to reserve space-time trajectories through intersections upon approach. This chapter explores modeling and optimization approaches for autonomous intersection management based on mixed-integer linear programming. A single intersection is first modeled and this formulation is then embedded in a decentralized framework to consider a network of intersections. A stability analysis is conducted to show how a modified max-pressure control policy can accommodate stochastic demands. Autonomous intersection management also opens the door to auction-based mechanisms which hold the potential to improve social welfare by allowing users to bid for priority service. This chapter reports on advances in the application of online mechanism design for developing incentive-compatible mechanisms for traffic intersection auctions. The problem set forth is the design of an online payment mechanism that is able to incentivize utility-maximizing users to reveal their true preferences upon bidding for service. To this end, two Markov chain models are proposed to estimate users’ expected waiting time. These models are incorporated into an online payment mechanism that is shown to be incentive-compatible in the dynamic sense. Numerical experiments are reported to show that, in contrast, static incentive-compatible mechanisms (i.e. that do not account for future arrivals) may fail to ensure truthful user behavior. This motivates the necessity to develop online mechanisms for such time-paced game-theoretic applications in transportation systems.
Perspectives are provided in a Chapter 5. Research directions that may build on and extend this thesis are outlined. In addition, an attempt is made to discuss emerging trends in optimization and game theory and their potential to improve transportation systems.

## Other contributions

Throughout my academic career, I have had the opportunity to work on several interesting topics in the field of transportation and logistics, but also in other fields of research such as construction engineering and epidemic modeling. A brief overview of my contributions other than those presented in the core chapters of this thesis is provided hereafter.

At the time I completed my Ph.D., research on shared mobility systems was in great demand. This lead me to work on one-way carsharing systems based on MILP formulations [Jian et al., 2016, 2019]. I participated in a couple of studies on rideharing systems that involved real-time computing [Najmi et al., 2017], as well as simulation-based optimization to identify user subsidies in a discrete choice modeling context [Song et al., 2021]. I had the insightful opportunity to participate in an experimental economics study on the design of lottery-based incentive mechanisms for promoting offpeak travel [Rey et al., 2016a]. In a more OR context, we explored the design of demand-aware park-and-ride mobility systems [Henry et al., 2022]. Motivated by industry collaborations, I also had the opportunity to work on vehicle routing problems: first in the context of food rescue and logistics where fairness considerations arise [Nair et al., 2016, 2017, Rey et al., 2018], and later in the context of rich, utility-aware dial-a-ride problems [Dong et al., 2020, 2022].

In the same vein as the works summarized in Chapter 2, I participated in several studies on network design problems, some of which involved bilevel optimization formulations and solution methods. In a series of studies, I contributed to study the impact of electric vehicles and the placement of charging infrastructure on network performance [Zhang et al., 2018, 2019a,b]. The impact of facility location decisions was also studied in a sustainability context using bilevel optimization approaches [Hammad et al., 2017b,a]. Other network design studies focused on transit network with elastic demand [Hossein Rashidi et al., 2016], school network capacity optimization [Barbara et al., 2021] and freeway network design in a mixed automated mobility context [Chakraborty et al., 2021].

In other disciplines, I have developed OR techniques to solve decision-making problems in civil engineering. This includes developing new MINLP formulations and algorithms for construction site layout planning problems [Hammad et al., 2016a,b, 2017c], as well as new modeling approaches for job allocation in construction organizations [Shahbazi et al., 2019]. I started working on outbreak modeling and control nearly a decade ago by contributing to risk analysis and modeling [Gardner et al., 2014] and also by proposing combinatorial optimization approaches to solve a variant of the Steiner tree problem [Rey et al., 2016b]. Other efforts in this research direction have been focused on developing multiscale network models to evaluate and optimize the impact of control strategies for mitigating global outbreaks [Chen et al., 2016, 2017, Zlojutro et al., 2019].

The rest of the thesis follows the plan outlined above which includes three core chapters and a final chapter wherein research perspectives are discussed. For conciseness, all mathematical proofs of propositions, lemmas and theorems are omitted, and readers are referred to the associated publications for these details.
# Chapter 2

## Bilevel optimization problems in transportation networks

This chapter is based on Rey et al. [2019], Rey [2020], Rey and Bar-Gera [2020].

### 2.1 Introduction

We first consider the discrete network design problem (DNDP) which is a challenging problem in transportation, introduced by Leblanc [1975]. The DNDP can be formulated as a bilevel optimization problem where the leader problem aims to identify the optimal network design to minimize network travel time and the follower problem represents network users' reaction, typically as a static traffic assignment problem (TAP) under user equilibrium [Wardrop, 1952].

The DNDP can be defined on a network with nodes $N$ and directed links $A$ as a multi-commodity network flow problem with nonlinear link travel time functions. Let $D$ be the set of destination nodes and $d_{is}$ be the demand from node $i \in N$ to destination node $s \in D \subseteq N$. If the pair $(i, s)$ is not an Origin-Destination (OD) pair in the network then $d_{is} = 0$ and to ensure flow conservation we set $d_{ss} = -\sum_{i \in N} d_{is}$. We denote $x_{ij,s}$ the flow of on link $(i, j) \in A$ travelling to destination $s \in D$, and $x_{ij}$ the total flow on link $(i, j) \in A$. Let $t_{ij}$ represent the travel time on link $(i, j) \in A$, typically modelled as a strictly convex function of the total link flow $x_{ij}$ to ensure the uniqueness of the equilibrium link flows. Let $A_1$ be the set of existing links and $A_2$ be the set of candidate links to improve the network, $A = A_1 \cup A_2$. For each link $(i, j) \in A_2$, let $g_{ij}$ be the cost of adding this link to the network and let $y_{ij} \in \{0, 1\}$ be the variable representing this choice. Let $B$ be the available budget for optimization. In the resulting formulation DNDP, $\mathcal{L}$ is the leader problem and $\mathcal{F}$ is the follower problem.

The leader problem $\mathcal{L}$ aims to minimize the total system travel time (TSTT) defined as the sum of $x_{ij}t_{ij}(x_{ij})$ over all links, subject to a budget constraint capturing the cost of link addition decisions $\mathbf{y}$—hereby referred to as the leader variable. The link flow pattern variable $\mathbf{x} = [x_{ij}]_{(i,j) \in A}$ is optimized in the follower problem $\mathcal{F}$, which is the traditional link-based TAP formulation under UE [Beckmann et al., 1956, Leblanc, 1975, Magnanti and Wong, 1984]. The impact of the leader variable $\mathbf{y}$ in the follower is achieved through the constraint $x_{ij} \le y_{ij}M$ wherein $M$ is an upper-bound on the total link flow $x_{ij}$.
$$ (\mathcal{L}) \min_{y} \sum_{(i,j) \in A} x_{ij} t_{ij}(x_{ij}), \\ \text{s.t.} \quad \sum_{(i,j) \in A_2} y_{ij} g_{ij} \leq B, \\ \qquad y_{ij} \in \{0, 1\}, \qquad \forall (i, j) \in A_2, $$

$$ (\mathcal{F}) \quad \begin{array}{ll} \boldsymbol{x} \in \arg \min_{\boldsymbol{x}} & \displaystyle \sum_{(i,j) \in A} \int_{0}^{x_{ij}} t_{ij}(v) dv, \\ \text{s.t.} & \displaystyle \sum_{j \in N:(i,j) \in A} x_{ij,s} - \sum_{j \in N:(j,i) \in A} x_{ji,s} = d_{is}, \quad \forall i \in N, \forall s \in D, \\ & \displaystyle \sum_{s \in D} x_{ij,s} = x_{ij}, \quad \forall (i,j) \in A, \\ & x_{ij} \leq y_{ij} M, \quad \forall (i,j) \in A_2, \\ & x_{ij,s} \geq 0, \quad \forall (i,j) \in A, \forall s \in D. \end{array} \tag{DNDP} $$

The goal of this chapter is threefold. A summary of the state-of-the-art of exact methodologies for the DNDP is first presented and a computational benchmark is conducted to highlight existing computational bottlenecks (Sections 2.2-2.3.3). An application of the DNDP for a network maintenance scheduling problem is then presented (Section 2.4) along with its extension to a network disaster recovery problem (Section 2.5).

## 2.2 Literature review

Several efforts have been proposed to solve the DNDP to global optimality. The seminal work of Leblanc [1975] introduced a Branch-and-Bound (B&B) algorithm for the DNDP which used the system-optimum (SO) relaxation of the TAP to find lower bounds. This relaxation requires fixing all unfixed $y$ variables to 1 to avoid Braess' paradox effects [Braess, 1968], and may thus lead to poor lower bounds. Gao et al. [2005] introduced a mixed-integer nonlinear programming (MINLP) approach based on generalized Benders' decomposition. The authors proposed a benchmark network with a single origin-destination (OD) pair, 12 nodes, 17 existing links and 6 candidate links, commonly referred to as *Gao's instance*. They also provided results for an instance based on the Sioux Falls network which contains 24 zones and nodes, 76 existing links and 5 candidate projects (which may involve one or more links) but computation time is not reported. It was noted by Farvaresh and Sepehri [2013] that despite claims of global optimality, the proposed approach of Gao et al. [2005] may converge to local optimums.

Bilevel optimization problems with a convex follower problem can be reformulated into single-level formulations by representing the follower problem using its Karush-Kuhn-Tucker (KKT) conditions and introducing binary variables to model complementarity slackness conditions [Bard, 2013]. This applies to the DNDP since the follower problem is a TAP which can be represented as a convex nonlinear program (NLP) [Beckmann et al., 1956]. Farvaresh and Sepehri [2011] proposed a mixed-integer linear programming (MILP), single-level reformulation of the DNDP obtained using piecewise linear approximations of link travel time functions. Numerical results are only reported for Gao's instance and an extended network based on Gao's instance which contains 16 nodes, 17 existing links and 25 candidate links, and
2 OD pairs. Farvaresh and Sepehri [2013] proposed a more scalable approach which extends the seminal B&B of Leblanc [1975] by solving the SO relaxation of the DNDP—instead of the TAP—as a mixed-integer nonlinear program (MINLP) to provide tighter lower bounds. This extended B&B algorithm was shown capable to solve instances with up to 100 nodes, 387 links and 15 candidate projects.

Luathep et al. [2011] proposed an SO-relaxation based approach wherein variational in-
equalities (VIs) are iteratively added to ensure UE conditions. Numerical results are re-
ported for Gao’s instance, as well as for an instance on Sioux Falls network with 5 candidate
projects, and a variation with 10 candidate projects for a mixed (discrete and continuous)
case. Fontaine and Minner [2014] proposed a Benders’ decomposition of a single-level MILP
reformulation of a linearized DNDP with piecewise linear approximations of link travel time
functions. Using this decomposition, the authors are able to solve this linearized DNDP
on a network of Berlin’s centre containing 36 zones, 398 nodes, 871 links and 10 candidate
projects.

Wang et al. [2013] studied an extended DNDP where the capacity of candidate links must also be decided. They present two global optimization algorithms which are based on the SO-relaxation of the DNDP, as well as a dynamic outer approximation of link travel time functions to derive lower bounds. Numerical results for a Sioux Falls network instance with up to 10 candidates links and 3 levels of capacity are reported. Wang et al. [2015] extended the DNDP to a variant involving both discrete and continuous decisions variables for adding links and determining their capacity, respectively. The authors used the VI formulation of Luathep et al. [2011] but combine it with an outer approximation of link travel time functions. Numerical results are only reported for Gao’s instance. Bagloee et al. [2017] proposed a B&B algorithm which uses a generalized Benders’ decomposition approach to solve the SO-relaxation of the DNDP at each node of the tree. Results on Sioux Falls and Winninpeg’s network are reported with up to 20 candidate projects.

Although several efforts have been proposed to solve the DNDP or a linear approximation
of this problem, the literature on exact or near-exact approaches remains scarce. Further,
there is no reference datasets for benchmarking solution methods which undermines the
research on this challenging bilevel optimization problem.

## 2.3 Computational benchmarking of exact methods for the DNDP

In this section, we discuss formulations and algorithms to solve the DNDP to optimality. We start by discussing the role of link travel time functions and attempt to categorize existing solution methodologies for the DNDP.

A major computational challenge inherent to the DNDP is the nonlinearity of the link travel time functions $t_{ij}(x_{ij})$, typically Bureau of Public Roads (BPR) functions. The consensus on link travel time functions is to use strictly convex functions of the form $t_{ij}(x_{ij}) = T_{ij} + c_{ij}x_{ij}^{e_{ij}}$ where $T_{ij}$ represents link free-flow travel time, $c_{ij}$ is a coefficient which captures link capacity and $e_{ij} \ge 1$ is an exponent; to ensure uniqueness of the link flow solution. The vast majority of existing data for link travel time functions of this form assumes $e_{ij} = 4$ for all $(i, j) \in A$. Link travel time functions are present in the objective function of both leader and follower problems in the form $x_{ij}t_{ij}(x_{ij})$ which involves the nonlinear term $x_{ij}^{e_{ij}+1}$.
This term, and sometimes term $x_{ij}^{e_{ij}}$, have been either approximated using piecewise linear functions [Farvaresh and Sepehri, 2011, Luathep et al., 2011, Fontaine and Minner, 2014] or using linear outer approximation schemes [Wang et al., 2013, 2015].

Piecewise linear approximations require discretizing the domain of link flow variable $x$ into $m$ disjunctive segments and requires $O(m)$ auxiliary variables to adjust the approximated nonlinear term according to the segment activated. The domain of link flows is typically defined as $x_{ij} \in [0, \bar{x}_{ij}]$ where $\bar{x}_{ij}$ is the maximum flow that can travel on link $(i, j)$. Determining tight values for the upper bound $\bar{x}_{ij}$ is not trivial and using a conservative value such as the total demand of the network may lead to poor approximation schemes. However, since the nonlinear terms $x_{ij}^{e_{ij}}$ and $x_{ij}^{e_{ij}+1}$ are convex on $[0, \bar{x}_{ij}]$, piecewise linear approximations of these terms do not require any integer variables, as noted by Farvaresh and Sepehri [2011] and in Fontaine and Minner [2014]. Hence, while piecewise linear approximation may require a significant amount of additional variables and linear constraints to achieve high-quality solutions, their impact on computational performance can be moderate.

Outer approximations of link travel time functions attempt to derive compact convex envelopes of nonlinear terms $x_{ij}^{e_{ij}}$ and $x_{ij}^{e_{ij}+1}$ using integer-linear constraints. These schemes are notoriously computationally challenging but guarantee that the original (non-linearized) DNDP is solved to optimality. Only in rare cases have link travel time functions been incorporated without any direct approximation scheme [Farvaresh and Sepehri, 2013, Bagloee et al., 2017].

We next discuss exact methodologies for the DNDP which have adopted one of the above approaches to handle link travel time functions. Since piecewise linear approximations of link travel time functions cannot be guaranteed to converge to optimal solutions unless $m \to \infty$, such methodologies are referred to as linearized DNDP approaches.

### 2.3.1 SO-relaxation based approaches

Despite its potential weak initial lower bound, the SO-relaxation of the DNDP has emerged as a powerful mechanism to conceive iterative solution methods for the DNDP. The SO-relaxation of the DNDP is a single-level optimization problem which ignores the follower objective function as summarized in SO-DNDP.

$$
\begin{array}{ll}
\min_{\mathbf{x},\mathbf{y}} & \displaystyle\sum_{(i,j) \in A} x_{ij} t_{ij}(x_{ij}), \\
\text{s.t.} & \displaystyle\sum_{(i,j) \in A_2} y_{ij} g_{ij} \le B, \\
& \displaystyle\sum_{j \in N:(i,j) \in A} x_{ij,s} - \sum_{j \in N:(j,i) \in A} x_{ji,s} = d_{is}, \quad \forall i \in N, \forall s \in D, \\
& \displaystyle\sum_{s \in D} x_{ij,s} = x_{ij}, \quad \forall (i,j) \in A, \\
& x_{ij} \le y_{ij} M, \quad \forall (i,j) \in A_2, \\
& y_{ij} \in \{0,1\}, \quad \forall (i,j) \in A_2, \\
& x_{ij,s} \ge 0, \quad \forall (i,j) \in A, \forall s \in D.
\end{array}
\qquad (\text{SO-DNDP})
$$

Solving SO-DNDP yields a lower bound on the optimum of DNDP which can serve as starting point of iterative schemes. Leblanc [1975] was the first to propose a customized
B&B algorithm which branches on unfixed $\mathbf{y}$ variables to yield subproblems that refine the initial lower bound. At each node of the B&B tree, Leblanc [1975] approach further relaxed SO-DNDP by temporarily fixing all unfixed $\mathbf{y}$ variables and solving the resulting SO-TAP using a convex programming algorithm, e.g. Frank-Wolfe. This B&B scheme was refined and extended by Farvaresh and Sepehri [2013] which proposed to solve SO-DNDP at each node of the tree using a global MINLP algorithm. Upper bounds are obtained by solving the TAP under UE conditions after obtaining $\mathbf{y}$ from SO-DNDP.

Bagloee et al. [2017] proposes a B&B algorithm which uses generalized Benders' decomposition approach to solve the SO-relaxation of the DNDP at each node of the tree. The proposed algorithm is parameterized using a value (denoted $\alpha$) which influences the solution of the TAP solved therein. The algorithm is not guaranteed to find global optimal solutions if $\alpha \neq 1$.

Luathep et al. [2011] proposed a cutting-plane algorithm to complement formulation SO-DNDP with VIs that characterizes UE conditions [Dafermos, 1980]. Let $\Omega_F$ be the polyhedron of the feasible region of the follower problem $\mathcal{F}$ and let $\mathbf{x}^*$ be the UE link flow pattern, the following VI holds:

$$ \sum_{(i,j) \in A} t_{ij}(x_{ij}^*)(x_{ij} - x_{ij}^*) \geq 0, \quad \forall \mathbf{x} \in \Omega_F. \qquad (2.1) $$

The approach of Luathep et al. [2011] is based on the observation that VIs corresponding to the set of extreme points of the TAP polyhedron is sufficient and necessary to characterize UE conditions. At each iteration, SO-DNDP is solved with a restricted set of VIs and shortest path problems are solved for each OD pair to identify violated VIs of the form (2.1), which are then added as cuts to SO-DNDP until none can be found.

Wang et al. [2013] observed that it was sufficient to iteratively forbid the last $\mathbf{y}$ solution found to construct a UE solution from the SO-relaxation of the DNDP. At each iteration, the proposed algorithm first solves SO-DNDP before calculating the UE cost of the corresponding $\mathbf{y}$ solution using any TAP algorithm to obtain an upper bound. The process is then repeated by adding an interdiction constraint to identify the next-best SO network design. Let $\mathbf{y}^k$ be the optimal solution of SO-DNDP at iteration $k$. At iteration $k+1$, SO-DNDP is solved with the interdiction constraint:

$$ \sum_{(i,j) \in A_2} (y_{ij}(1 - y_{ij}^n) + (1 - y_{ij})y_{ij}^n) \geq 1, \quad \forall n \in \{1, \dots, k\}. \qquad (2.2) $$

The iterative process is repeated until the resulting objective value is greater than or equal to the upper bound.

### 2.3.2 KKT conditions-based approaches

Since the TAP can be formulated as a convex optimization problem, an intuitive approach to solve the DNDP is to replace the follower problem $\mathcal{F}$ by its KKT conditions. Farvaresh and Sepehri [2011] proposed such a direct approach wherein the KKT conditions of the TAP are represented with auxiliary variables. Let $\pi_{is} \geq 0$ be the travel time or path travel time from node $i \in N$ to destination $s \in D$. The UE conditions of the TAP require:
<table><caption>Table 2.1: Summary of exact methodologies for the DNDP.</caption><thead><tr><th>Paper</th><th>Method</th><th>Link travel time function approximation</th><th>Type(s) of optimization problem solved</th></tr></thead><tbody><tr><td>Leblanc [1975]</td><td>SO-relaxation and B&B</td><td>N/A</td><td>Convex NLP</td></tr><tr><td>Gao et al. [2005]</td><td>Generalized Benders' decomposition</td><td>N/A</td><td>MINLP and Convex NLP</td></tr><tr><td>Farvaresh and Sepehri [2011]</td><td>KKT conditions</td><td>Piecewise linear</td><td>MILP</td></tr><tr><td>Luathep et al. [2011]</td><td>SO-relaxation and VIs</td><td>Piecewise linear</td><td>MILP and LP</td></tr><tr><td>Farvaresh and Sepehri [2013]</td><td>SO-relaxation and B&B</td><td>N/A</td><td>MINLP and Convex NLP</td></tr><tr><td>Wang et al. [2013]</td><td>SO-relaxation and interdiction cuts</td><td>Outer linear</td><td>MILP and Convex NLP</td></tr><tr><td>Fontaine and Minner [2014]</td><td>KKT conditions and Benders' decomposition</td><td>Piecewise linear</td><td>MILP and LP</td></tr><tr><td>Wang et al. [2015]</td><td>SO-relaxation and VIs</td><td>Outer linear</td><td>MILP and LP</td></tr><tr><td>Bagloee et al. [2017]</td><td>B&B and Generalized Benders' decomposition</td><td>N/A</td><td>MINLP and Convex NLP</td></tr></tbody></table>

$$ 
\begin{aligned}
t_{ij}(x_{ij}) - \pi_{is} + \pi_{js} &\ge 0, && \forall (i,j) \in A, \forall s \in D, && (2.3a) \\
x_{ij}(t_{ij}(x_{ij}) - \pi_{is} + \pi_{js}) &= 0, && \forall (i,j) \in A, \forall s \in D. && (2.3b)
\end{aligned}
$$

The complementarity slackness conditions (2.3b) is nonlinear and typical mathematical programming approaches require additional binary variables to obtain an integer-linear form suitable for MILP. Farvaresh and Sepehri [2011] observe that this surplus of binary variables yields considerable computational challenge. Fontaine and Minner [2014] proposed an alternative approach to incorporate the KKT conditions of the TAP which consists of replacing the objective of the follower problem $\mathcal{F}$ with a primal-dual constraint which requires a null duality gap. This method is applied to a linearized DNDP, which results in a single-level MILP composed of leader, primal follower and dual follower constraints and variables, with the addition of the primal-dual constraint. This MILP does not require any additional binary variable and the authors propose a Benders' decomposition approach.

A summary of the solution methods discussed is provided in Table 2.1. Most of these methodologies can be adapted to work with piecewise linear or outer- approximations of link travel time functions. While the former only solves the so-called linearized DNDP, piecewise linear approximations are often capable of identifying global optimal solutions of the original problem and are, at the same time, substantially easier to implement.

### 2.3.3 Numerical experiments for the linearized DNDP

We implement three solution methods for the linearized DNDP: i) the B&B algorithm of Farvaresh and Sepehri [2013] which is an extension of Leblanc [1975]—SOBB; ii) the SO-relaxation based algorithm of Wang et al. [2013] with interdiction cuts—SOIC, and iii) the primal-dual formulation of Fontaine and Minner [2014] as a single-level MILP (without Benders' decomposition)—MKKT. All three solution methods, SOBB, SOIC and MKKT, are implemented using the piecewise linear approximation of link travel time functions proposed
by Farvaresh and Sepehri [2011]. Hence, in all cases, the optimization problems solved are MILPs and LPs (TAPs within SOBB and SOIC are solved in their linearized form). This differs from the original implementation of Farvaresh and Sepehri [2013] and Wang et al. [2013], wherein outer approximation schemes are used. In addition, the method MKKT is implemented as direct MILP approach unlike the Benders' decomposition scheme proposed in Fontaine and Minner [2014]. To measure the quality of the approximated solutions, the flow pattern corresponding to the best (lowest leader objective value) y solution among all three methods is calculated by solving the TAP as a convex problem.

All methods are implemented in Python. All MILPs and LPs are solved using CPLEX 12.8. Convex TAPs are solved using the Pyomo module and IPOPT. All solution methods were tested and implemented on the same Windows 7 machine with 16Gb of RAM and a CPU of 2.7Ghz, in a single-thread mode with a time limit of 10 minutes. The upper bound on link flows $\bar{x}_{ij}$ was set to $1e^5$ and this value is also used for $M$. The number of segments used in the piecewise linear approximations of link travel time functions is $m = 100$. A scaling factor of $1e^{-3}$ is used to scale travel demand and link capacities as it was found to improve computational performance. For reproducibility purposes, all implemented optimization formulations and codes are provided at https://github.com/davidrey123/DNDP.

New instances for the DNDP based on Sioux Falls network<fnref n="1" /> have been designed to test the implemented solution methods. For these new instances, a total of 15 pairs of new links (total of 30 links) have been created, along with their performance characteristics and addition costs. Among these new links, 5 of the 15 pairs are identical to the ones used by Luathep et al. [2011]. A total of 20 DNDP instances have been created, 10 of these contain 10 new links (5 pairs) and are named SF_DNDP_10, and the remaining contain 20 new links (10 pairs) and are named SF_DNDP_20. In all our numerical experiments, link addition variables are not paired, i.e. the size of $A_2$ is equal to the total number of new links which may or may not be added in pairs; hence instances SF_DNDP_10 require 10 binary variables and instances SF_DNDP_20 require 20 binary variables. All instance data is provided in the public repository available at https://github.com/davidrey123/DNDP in the TNTP format. An extra column has been appended to provide link cost data ($g_{ij}$) and a value of 0 indicates that the link is part of the original network (before optimization).

Two numerical experiments are conducted on instances SF_DNDP_10 and SF_DNDP_20: a budget sensitivity analysis wherein the available budget $B$ is equal to 25%, 50% and 75% of the total cost $\sum_{(i,j) \in A_2} g_{ij}$, and a demand sensitivity analysis wherein the budget is fixed to 50% of the total cost and travel demand is set to 50%, 100% and 150% of the base demand. The average performance of all three methods implemented, i.e. SOBB, SOIC and MKKT, is reported in Tables 2.2 and 2.3 for budget and demand sensitivity analysis, respectively.

We find that 10-project instances SF_DNDP_10 can present considerable computational challenges, especially for SO-relaxation based methods when the budget available is relatively high ($B\% = 75\%$). Increasing the number of candidate links to 20 results in substantial difficulties in most cases, regardless of the method used, since most of the time the methods implemented were not able to converge within the 10-minute time limit. The results of the budget sensitivity analysis on SF_DNDP_10 instances show that MKKT is often faster than other methods at medium and high budgets ($B\% = 50\%, 75\%$). In turn, SOIC and SOBB tend to time-out more frequently at these budgets, which is likely due to the more substantial

<footnote n="1">¹These instances are based on Sioux Falls network as available at https://github.com/bstabler/TransportationNetworks</footnote>
<table><caption>Table 2.2: Budget sensitivity experiment. Methods SOBB, SOIC and MKKT are implemented for instances SF_DNDP_10 and SF_DNDP_20 for a budget B equal to B% = 25%, 50% and 75% of the total cost $\sum_{(i,j) \in A_2} g_{ij}$. The time limit is 10 minutes. AT is the average runtime in seconds, AG is the average relative optimality gap upon termination in % and TO is the proportion of time-outs in %.</caption><thead><tr><th rowspan="3">Method</th><th colspan="6">SF_DNDP_10</th><th colspan="6">SF_DNDP_20</th></tr><tr><th colspan="3">B% = 25%</th><th colspan="3">B% = 50%</th><th colspan="3">B% = 75%</th><th colspan="3">B% = 25%</th><th colspan="3">B% = 50%</th><th colspan="3">B% = 75%</th></tr><tr><th>AT</th><th>AG</th><th>TO</th><th>AT</th><th>AG</th><th>TO</th><th>AT</th><th>AG</th><th>TO</th><th>AT</th><th>AG</th><th>TO</th><th>AT</th><th>AG</th><th>TO</th><th>AT</th><th>AG</th><th>TO</th></tr></thead><tbody><tr><td>SOBB</td><td>56.8</td><td>0.00</td><td>0</td><td>333.8</td><td>0.39</td><td>20</td><td>484.9</td><td>1.75</td><td>60</td><td>600</td><td>2.48</td><td>100</td><td>600</td><td>4.76</td><td>100</td><td>600</td><td>6.07</td><td>100</td></tr><tr><td>SOIC</td><td>33.8</td><td>0.00</td><td>0</td><td>283.4</td><td>0.31</td><td>20</td><td>357.5</td><td>0.69</td><td>40</td><td>600</td><td>3.76</td><td>100</td><td>600</td><td>4.44</td><td>100</td><td>600</td><td>5.26</td><td>100</td></tr><tr><td>MKKT</td><td>57.2</td><td>0.00</td><td>0</td><td>170.6</td><td>0.00</td><td>0</td><td>155.8</td><td>0.00</td><td>0</td><td>600</td><td>5.68</td><td>100</td><td>600</td><td>6.46</td><td>100</td><td>600</td><td>6.21</td><td>100</td></tr></tbody></table>

<table><caption>Table 2.3: Demand sensitivity experiment. Methods SOBB, SOIC and MKKT are implemented for instances SF_DNDP_10 and SF_DNDP_20 for a budget B equal to 50% of the total cost $\sum_{(i,j) \in A_2} g_{ij}$ and a demand of D% = 50%, 100% and 150% of the base demand. The time limit is 10 minutes. AT is the average runtime in seconds, AG is the average relative optimality gap upon termination in % and TO is the proportion of time-outs in %.</caption><thead><tr><th rowspan="3">Method</th><th colspan="6">SF_DNDP_10</th><th colspan="6">SF_DNDP_20</th></tr><tr><th colspan="3">D% = 50%</th><th colspan="3">D% = 100%</th><th colspan="3">D% = 150%</th><th colspan="3">D% = 50%</th><th colspan="3">D% = 100%</th><th colspan="3">D% = 150%</th></tr><tr><th>AT</th><th>AG</th><th>TO</th><th>AT</th><th>AG</th><th>TO</th><th>AT</th><th>AG</th><th>TO</th><th>AT</th><th>AG</th><th>TO</th><th>AT</th><th>AG</th><th>TO</th><th>AT</th><th>AG</th><th>TO</th></tr></thead><tbody><tr><td>SOBB</td><td>537.7</td><td>0.54</td><td>80</td><td>336.6</td><td>0.40</td><td>20</td><td>79.2</td><td>0.00</td><td>0</td><td>600</td><td>2.38</td><td>100</td><td>600</td><td>4.75</td><td>100</td><td>432.2</td><td>0.24</td><td>40</td></tr><tr><td>SOIC</td><td>513.8</td><td>0.56</td><td>80</td><td>284.6</td><td>0.34</td><td>20</td><td>14.6</td><td>0.00</td><td>0</td><td>600</td><td>2.14</td><td>100</td><td>600</td><td>4.44</td><td>100</td><td>254.9</td><td>0.12</td><td>30</td></tr><tr><td>MKKT</td><td>87.9</td><td>0.00</td><td>0</td><td>172.1</td><td>0.00</td><td>0</td><td>65.9</td><td>0.00</td><td>0</td><td>600</td><td>2.58</td><td>100</td><td>600</td><td>6.46</td><td>100</td><td>519.8</td><td>1.43</td><td>50</td></tr></tbody></table>

gap between UE and SO flow patterns when a high budget is available. At low budget (B% = 25%), SOIC outperforms other methods. Instances SF_DNDP_20 result systematically in time-outs and the average optimality gap upon termination tends to increase with the budget. The performance of all three methods is of the same order of magnitude, with optimality gaps in the range of 3%-6%. The demand sensitivity analysis highlights the good performance of MKKT at low demand (D% = 50%) compared to other methods on SF_DNDP_10 instances. Increasing the demand to 150% of the base demand yields substantially easier problems, notably for SF_DNDP_20 instances, and SOIC is found to outperform other methods in terms of number of time-outs. To assess the quality of the linear approximation of link travel time functions, we compare the network travel time obtained by solving the TAP in its convex NLP form for the best y solution found among all three methods. The deviation observed is on average of 0.28% with a standard deviation of 0.25%.

Overall, this study shows that even medium-size instances for the linearized DNDP can present considerable computational challenges. The numerical experiments revealed that some 10-project instances and most 20-project instances could not be solved within the imposed time limit (10 minutes). Although this time limit is relatively small, all methods were implemented using a piecewise linear approximation of link travel time functions which is expected to outperform outer approximations approaches.
## 2.4 Maintenance scheduling in transportation networks

Road maintenance and repair operations are conducted regularly in major cities. This is due to the numerous links that connect the transport network as well as to the high travel demand which contributes in deteriorating road quality. In addition, the organization of urban areas often relies on the road network to incorporate utility networks, such as electricity, telecommunications and water networks. As a result, lane or road closures—hereby referred to as maintenance projects—are frequent and can lead to significant delays within the city if they are conducted without any coordination across the network. In this section, we address the problem of determining an optimal schedule for the conduction of maintenance projects in urban transportation networks. The impact of link capacity reductions induced by maintenance projects onto network performance is measured by the TSTT. We propose a novel solution method for the Network Maintenance Scheduling Problem (NMSP) that finds a schedule which minimizes the TSTT over a planning period subject to project schedule constraints.

The first effort to account for traffic dynamics in network maintenance scheduling can be attributed to Chang et al. [2001] who proposed a tabu search approach to solve the optimization problem. The authors implement their approach on a 6-project scenario derived from the Columbus, Ohio network which comprises 5,441 nodes and 12,658 links. All projects are of equal duration and at most 2 projects can be conducted in parallel (i.e. 2 teams are available), thus considerably reducing the number of possible patterns. Subsequent works have mainly focused on meta-heuristics approaches to solve this bilevel scheduling problem. Ng et al. [2009] used a cell transmission model to account for users' reaction and proposed a genetic algorithm-based solution method. The authors implemented the proposed methodology on a test network with 13 nodes and 24 links where each link can be eligible for yearly maintenance over a 3-year planning period. Lee [2009] considered a framework that combines ant colony optimization with a microscopic traffic simulation model to represent traffic delays caused by maintenance projects. The framework is implemented on a realistic case study derived from a sewer system construction project on a network with 292 nodes and 178 links. A total of 35 maintenance projects are considered with a maximum of 5 teams over a planning period of 270 days. A similar approach is used by Xiong et al. [2014] where an ant colony optimization approach is used along with a static traffic assignment model to account for users' reaction. The proposed approach is tested on a 6-project scenario derived from the Chicago Sketch network which includes 933 nodes and 2,950 links over a planning period of 36 months. Chu and Chen [2012] proposed a hybrid dynamic model to simulate threshold-based maintenance based on pavement life-cycle. The bilevel optimization formulation is solved using a modified tabu search algorithm on a network with 70 nodes and 270 links. A 10-year planning period is considered over which all links are eligible for maintenance operations. Zheng et al. [2012] introduced a detailed modeling framework for representing maintenance project features, such as day/night time construction strategies, and used an approximation of the UE travel times to account for traffic delays. Fontaine and Minner [2017] proposed an exact approach for the bilevel network maintenance planning problem which is based on a Benders' decomposition method of a single-level reformulation of the original problem. The optimality conditions of the TAP are used as equilibrium constraints. In their formulation, they assume that a budget is available for maintenance costs. Numerical results on the Sioux Falls network and an extended Sioux Falls network of 100 nodes of 317 links are reported for planning periods of 10 and 8 time steps, respectively.

We propose a new bilevel programming formulation and develop a solution methodology
ogy based on maintenance project patterns enumeration. This representation of the NMSP
involves a number of decision variables which grows exponentially with the number of main-
tenance projects. We propose a new Branch-and-Price (B&P) algorithm to solve this bilevel
problem. Namely, we introduce customized branching and bounding rules, and develop a
Column Generation (CG) procedure to solve the Linear Programming (LP) relaxation of a
restricted problem at each node of the tree. Problem-specific pricing algorithms are proposed
to iteratively introduce profitable project patterns in the formulation. In addition, we intro-
duce a statistical regression model to approximate congestion effects and derive approximated
lower bounds within the optimization framework.

2.4.1 Network maintenance scheduling

The coordination of maintenance projects within a transport network must satisfy operational
constraints (e.g. spatial and temporal restrictions) and should have a minimal impact onto
network performance. Specifically, we focus on a planning period (e.g. one year) during
which all projects must be completed in consecutive time periods (e.g. one week). The
consecutiveness assumption is consistent with the operations of traffic management centers
in which road maintenance projects are conducted without interruption from their starting
date. It is important to stress that we are mostly concerned with maintenance projects
that have an impact during peak-hour traffic. Strictly overnight work and other off-peak
road maintenance projects with no impact on traffic conditions may be scheduled without
accounting for their influence on network performance. Hence, in this section, we focus on
maintenance projects that span one or more time periods and reduce the capacity of the
associated links during the entire duration of the project.

The proposed network maintenance scheduling problem can be defined on a directed graph
G = (N, A) where the set of nodes N represents travel demand centroids or intersections,
and the set of arcs A represents links among nodes. We assume that the set of projects is
known at the beginning of the planning horizon and we assume that the recovery time of each
project is known deterministically. The network travel time at each time period is determined
based on network link states and its associated traffic assignment pattern determined under
the UE principle. We denote Qw the demand for a specific OD pair w ∈ W, Πw its set of
paths and fw the flow on path k ∈ Πw. In addition, we denote δw the fixed binary value
indicating if link a ∈ A belongs to path k ∈ Πw (1) or not (0). Finally, we denote xa the flow
on link a ∈ A. Let ma denote the state of link a ∈ A. The TAP corresponding to a network
state vector m is summarized in (2.4). Note that bold symbols are used to denote vectors
throughout the chapter.

$$
x^*(m) \in \arg \min \quad \sum_{a \in A} \int_{0}^{x_a} t_a(v, m_a) dv, \qquad (2.4a)
$$

$$
\text{s.t.} \quad \sum_{k \in \Pi_w} f_k^w = Q_w, \qquad \forall w \in W, \tag{2.4b}
$$

$$
\sum_{w \in W} \sum_{k \in \Pi_w} f_k^w \delta_{a,k}^w = x_a, \qquad \forall a \in A, \tag{2.4c}
$$

$$
f_k^w \ge 0, \qquad \forall w \in W, k \in \Pi_w, \quad (2.4d)
$$

$$
x_a \ge 0, \qquad \forall a \in A. \tag{2.4e}
$$
Observe that the link flow pattern $\boldsymbol{x}^*(\boldsymbol{m})$ depends on the network state vector $\boldsymbol{m}$. The network travel time under UE, also known as the total system travel time, corresponding to network state $\boldsymbol{m}$, denoted $TSTT(\boldsymbol{m})$, is thus calculated as follows:

$$TSTT(\boldsymbol{m}) = \boldsymbol{x}^*(\boldsymbol{m})^\top \boldsymbol{t}(\boldsymbol{x}^*(\boldsymbol{m}), \boldsymbol{m}) = \sum_{a \in A} x_a^*(\boldsymbol{m}) t_a(x_a^*(\boldsymbol{m}), m_a). \quad (2.5)$$

Let $P$ be the set of network maintenance projects that must be conducted and let $T$ be the set of time periods within the planning period considered. Each maintenance project $p \in P$ is defined by the following attributes: the set of links affected by the project $A_p \subseteq A, \forall p \in P$; link capacity reduction ratios $R_a \in [0, 1], \forall a \in A_p$; and project duration $D_p \in \{1, \dots, |T|\}$. The set of links affected by each project is typically a subset of the network links and we assume that each link can only be affected by at most one project during the planning period. The link capacity reduction $R_a$ is assumed to be a percentage of the maximal capacity of link $a \in A$. A reduction of 0% for all links of the network ($R_a = 0, \forall a \in A$) corresponds to the base scenario in which no project is active and all links operate at maximal capacity. If a project $p$ is active, then all links $a \in A_p$ see their capacity reduced by $100 \times R_a\%$. Each project $p$ must be active for $D_p$ consecutive time periods and completed within the planning period. Hence, for each project $p \in P$, the set of feasible start times $F_p$ can be defined as:

$$F_p = \{0, \dots, |T| - D_p\}. \quad (2.6)$$

We introduce binary decision variables to represent project start times:

$$g_{s,p} = \begin{cases} 1 & \text{if project } p \text{ starts at time period } s, \\ 0 & \text{otherwise,} \end{cases} \quad \forall p \in P, \forall s \in F_p. \quad (2.7)$$

and impose the project schedule constraints:

$$\sum_{s \in F_p} g_{s,p} = 1, \quad \forall p \in P. \quad (2.8)$$

To connect maintenance projects with the network performance metric (TSTT) we propose to enumerate all the capacity patterns $c$ resulting from all possible maintenance project combinations. The number of possible project combination is $2^{|P|}$. For each project combination $\bar{P} \subseteq P$—i.e. all projects in $\bar{P}$ are active—the capacity of the links $\cup_{p \in \bar{P}} A_p$ is reduced. We refer to the corresponding link capacity vector $c$ as a capacity pattern. Since links are associated with maintenance projects, it is more adequate to handle project patterns than capacity patterns. Let $\Sigma = \{0, \dots, 2^{|P|} - 1\}$ be the set of all project patterns, for each pattern $\sigma \in \Sigma$, we denote $TSTT(\sigma)$ its associated TSTT under UE conditions. That is, $TSTT(\sigma)$ represents the total network travel time corresponding to the unique TAP link flow solution under UE if pattern $\sigma$ is active. Observe that each project pattern $\sigma$ corresponds to a unique network state vector $\boldsymbol{m}$ with $TSTT(\sigma) = TSTT(\boldsymbol{m})$. Further, each pattern can be represented as a binary vector where each component indicates if the corresponding project is active in this pattern or not. For instance, consider a scenario with 10 projects: $P = \{1, \dots, 10\}$; the pattern consisting of all odd numbered projects is $\sigma = \{1, 3, 5, 7, 9\}$ and this pattern can be encoded as the binary vector $\varphi = [1, 0, 1, 0, 1, 0, 1, 0, 1, 0]$ where $\varphi_p \in \{0, 1\}$ indicate the activeness of $p \in P$—we henceforth refer to $\varphi$ as the binary representation of $\sigma$. Throughout this chapter we use alternatively $\boldsymbol{m}$, $\sigma$ and $\varphi$ to denote a pattern whenever convenient.
The objective of the NMSP is to minimize the overall TSTT over the planning period T. Hence a possible formulation for the NMSP is to view the decision problem as that of selecting the optimal sequence of patterns subject to project schedule constraints. We introduce binary decision variables to represent the pattern selection at each time interval:

$$k_{\sigma}^{t} = \begin{cases} 1 & \text{if pattern } \sigma \text{ is selected at time period } t, \\ 0 & \text{otherwise,} \end{cases} \quad \forall \sigma \in \Sigma, \forall t \in T. \quad (2.9)$$

and impose the pattern selection constraints:

$$\sum_{\sigma \in \Sigma} k_{\sigma}^{t} = 1, \quad \forall t \in T. \quad (2.10)$$

The objective function of the upper level of the NMSP can then be formulated as:

$$\min \sum_{t \in T} \sum_{\sigma \in \Sigma} TSTT(\sigma) k_{\sigma}^{t}. \quad (2.11)$$

To link the project schedule variables with the pattern selection variables, we introduce two binary matrices $[\gamma_{\sigma,p}]$ and $[\zeta_{s,p}^t]$ defined as:

$$\gamma_{\sigma,p} = \begin{cases} 1 & \text{if project } p \text{ is active in pattern } \sigma, \\ 0 & \text{otherwise,} \end{cases} \quad \forall p \in P, \forall \sigma \in \Sigma, \quad (2.12)$$

$$\zeta_{s,p}^{t} = \begin{cases} 1 & \text{if } s \le t \le s + D_p - 1, \\ 0 & \text{otherwise,} \end{cases} \quad \forall p \in P, \forall s \in F_p, \forall t \in T. \quad (2.13)$$

Matrix $[\gamma_{\sigma,p}]$ maps maintenance projects to patterns while matrix $[\zeta_{s,p}^t]$ maps a project start time $s \in F_p$ to $D_p$ consecutive time periods. The linking constraints are:

$$\sum_{\sigma \in \Sigma} k_{\sigma}^{t} \gamma_{\sigma,p} = \sum_{s \in F_p} g_{s,p} \zeta_{s,p}^{t}, \quad \forall p \in P, \forall t \in T. \quad (2.14)$$

Let $g$ denote the vector of variables $g_{s,p}, \forall p \in P, \forall s \in F_p$ and $k$ denote the vector of variables $k_{\sigma}^t, \forall \sigma \in \Sigma, \forall t \in T$. The space delimited by constraints (2.8), (2.10) and (2.14) represents the feasible region for the binary variables $g$ and $k$. We show that, under these constraints, the integrality constraints on variables $k$ can be relaxed.

**Proposition 2.1.** *The integrality constraints on variables k can be relaxed.*

Proposition 2.1 allows us to relax the integrality constraints on variables $k_{\sigma}^t$ and helps in significantly reducing the number of subproblems in a Branch-and-Bound (B&B) framework. Further, because of the pattern selection constraint (2.10), the upper bound on variables $k_{\sigma}^t$ can be dropped. The resulting formulation for the NMSP is summarized in Model 2.1.
**Model 2.1 (NMSP).**

$$z_{NMSP} = \min \sum_{t \in T} \sum_{\sigma \in \Sigma} TSTT(\sigma) k_{\sigma}^{t}, \qquad (2.15)$$

$$\text{s.t.} \quad \sum_{s \in F_p} g_{s,p} = 1, \qquad \forall p \in P, \qquad (2.16)$$

$$\sum_{\sigma \in \Sigma} k_{\sigma}^{t} = 1, \qquad \forall t \in T, \qquad (2.17)$$

$$\sum_{\sigma \in \Sigma} k_{\sigma}^{t} \gamma_{\sigma,p} = \sum_{s \in F_p} g_{s,p} \zeta_{s,p}^{t}, \qquad \forall p \in P, \forall t \in T, \qquad (2.18)$$

$$k_{\sigma}^{t} \ge 0, \qquad \forall \sigma \in \Sigma, \forall t \in T, \qquad (2.19)$$

$$g_{s,p} \in \{0, 1\}, \qquad \forall p \in P, \forall s \in F_p. \qquad (2.20)$$

Model 2.1 is a bilevel MIP with an exponential number of variables (recall that $|\Sigma| = 2^{|P|}$). Further each weight $TSTT(\sigma)$ in the objective function must be obtained by solving the TAP for the corresponding pattern. Being a combination of MIP and bilevel programming, both of which are difficult optimization problems, the NMSP is intractable for large instances and an efficient solution methodology is required to solve the problem on realistic scenarios.

To solve the NMSP using the representation summarized in Model 2.1 we propose to work with a restricted set of pattern variables and use a Column Generation (CG) algorithm to introduce additional project patterns as required. Let $\lambda^t$ be the dual variable associated to constraint (2.17) and let $\mu_p^t$ be the dual variable associated to constraint (2.18). To price unused pattern variables $k_{\sigma}^t$ we need to determine their reduced-cost $r_{\sigma}^t$, given by:

$$r_{\sigma}^{t} = TSTT(\sigma) - \lambda^{t} - \sum_{p \in P} \mu_{p}^{t} \gamma_{\sigma,p}. \qquad (2.21)$$

Therefore determining the reduced-cost of a pattern variable requires solving the TAP for this particular pattern.

A solution to the NMSP is globally optimal if and only if the reduced cost of all variables is positive. To efficiently explore the feasible region of Model 2.1, we introduce the concept of 1-distance optimality. Let $d_1(\sigma)$ be the 1-distance neighborhood of pattern $\sigma$, i.e. the set of patterns obtained by modifying each bit (one at a time) of the binary representation of $\sigma$ to its complementary value. We say that a solution is 1-distance optimal if the reduced cost of all patterns in the 1-distance neighborhood of the solution is positive. Let $\Sigma^* = \{\sigma \in \Sigma : \exists t \in T, k_{\sigma}^t > 0\}$ be the set of patterns used in the solution; and let $T_{\sigma} = \{t \in T : k_{\sigma}^t > 0\}, \forall \sigma \in \Sigma$ be the set of time periods at which pattern $\sigma$ is used in the solution. For any solution, the condition of 1-distance optimality can be checked efficiently using Proposition 2.2.

**Proposition 2.2.** Let $\sigma \in \Sigma^*$ and let $\sigma' \in d_1(\sigma)$ where $\sigma'$ is obtained from $\sigma$ by modifying the activeness of project $p' \in P$. Define:

$$\mu_{p'}^{*} = \begin{cases} + \max\{\mu_{p'}^{t} : t \in T_{\sigma}\}, & \text{if } \gamma_{\sigma,p'} = 0, \\ - \min\{\mu_{p'}^{t} : t \in T_{\sigma}\}, & \text{otherwise.} \end{cases} \qquad (2.22)$$
There exists $t^* \in T$ such that $r_{\sigma'}^{t^*} < 0$ if and only if $TSTT(\sigma) - TSTT(\sigma') + \mu_{p'}^* > 0$.

Since the NMSP is a mixed-integer program, the CG-based approach must be integrated within a B&P algorithm to guarantee integral solutions. Further, the lower-level must be accounted for within the optimization framework. The proposed B&P methodology is composed of several components. Specifically, the B&P algorithm requires a method to construct an initial feasible solution and initialize the columns pool, branching and node selection rules, a CG procedure to solve the LP relaxation of the Restricted Master Problem (LRMP) at each node of the B&P tree, pricing algorithms to identify negative reduced-cost pattern variables. In addition, a statistical regression model to approximate congestion effects will be integrated within the B&P algorithm to derive approximated lower bounds during the optimization.

We next introduce the concept of first-order TSTT effects of maintenance projects and discuss how this can be used to inform the solution method for the NMSP.

### 2.4.2 First-order TSTT effects of maintenance projects

To measure the level of interaction among projects we examine the variation of the TSTT when projects are active and inactive. This approach is highly combinatorial: for instance, to measure the variation of the TSTT when $n_p$ projects are either active or inactive, we need to consider $2^{|P|-n_p}$ subpatterns and, for each of them, compare the resulting TSTT values when the $n_p$ projects are either active or inactive. To investigate the interaction among maintenance projects, we thus focus on low interaction levels. Measuring project interaction also plays a central role in the proposed solution methodology for the NMSP.

For presentation purposes, we introduce the following notation: we denote $\varphi_p \in \{0, 1\}$ the activeness of project $p \in P$. We introduce the concatenation operator $\oplus$ to represent the concatenation of mutually exclusive projects in a pattern. Further, we abuse notation and use $\tilde{\sigma}$ to denote a subpattern, i.e. a pattern with a subset of projects, whenever necessary. For instance, $\sigma = [\varphi_p = 1, \varphi_q = 0] \oplus \tilde{\sigma}$ is a project pattern where $p$ is active, $q$ is inactive and the activeness of the remaining projects is given by $\tilde{\sigma}$ (which is a subpattern of the subset of projects $P \setminus \{p, q\}$).

Using the above notation, we define the first order TSTT effects of a project $p$ with regards to subpattern $\tilde{\sigma}$, $\Delta_{p,\tilde{\sigma}}$:

$$ \Delta_{p,\tilde{\sigma}} = TSTT([\varphi_p = 1] \oplus \tilde{\sigma}) - TSTT([\varphi_p = 0] \oplus \tilde{\sigma}). \quad (2.23) $$

We denote $\Delta_p$ the vector of first order TSTT effects of project $p$ for all subpatterns of the subset of projects $P \setminus \{p\}$, i.e. $\Delta_p = [\Delta_{p,\tilde{\sigma}}]_{\tilde{\sigma}}$. $\Delta_p$ represents the influence of a single project onto the TSTT experienced in the network.

We next present a tailored experimental design based on first-order TSTT information to construct a statistical regression model that can be used to approximate TSTT values.

### 2.4.3 Statistical approximation of TSTT effects

Throughout this section, we assume that a set of observed TSTT values and their corresponding project patterns is available. This data will be eventually obtained in a pre-processing
step prior to the execution of the B&P algorithm which is presented at the end of the section. We denote $\Sigma_{obs} \subset \Sigma$ the set of observed patterns.

The TSTT corresponding to a project pattern $\sigma$ can be approximated using regression analysis by using project activeness variables $\varphi_p, \forall p \in P$ as explanatory variables. Due to the nonlinear congestion effects inherent to traffic networks, we propose an exponential regression model of the form:

$$TSTT(\varphi) = e^{b_y + \sum_{p \in P} b_p \varphi_p + \epsilon_{\varphi}}. \qquad (E)$$

where $\varphi = [\varphi_p]_{p \in P}$ is the binary representation of a pattern, $b_y, b_p, \forall p \in P$ are regression coefficients for each project and a y-intercept, and $\epsilon_{\varphi}$ is a normally distributed error on the prediction of pattern $\varphi$.

We use an Ordinary Least Squares (OLS) regression to fit Model (E) to observed TSTT data. The resulting statistical model can then provide some information on the TSTT corresponding to a non-observed pattern, i.e. for which the TAP has not been solved. This procedure can be helpful in the pricing subproblem within the CG algorithm since the reduced-cost of a pattern variable depends on its TSTT value. To increase the likelihood that negative reduced-cost pattern variables are correctly identified, we propose to use a scaling factor $\bar{\psi}$ to approximate the prediction error $\epsilon_{\varphi}$ in the statistical model. We define the approximated TSTT of pattern $\varphi$ as:

$$aTSTT(\varphi) = \bar{\psi}e^{b_y + \sum_{p \in P} b_p \varphi_p}. \qquad (2.24)$$

The value of $\bar{\psi}$ should be chosen such that $aTSTT(\varphi)$ is a lower bound on $TSTT(\varphi)$ for most patterns. Assuming that $\bar{\psi}$ is chosen accordingly, this approximation can be used to efficiently determine lower bounds within the resolution of the NMSP. This stems from the functional form of $aTSTT(\varphi)$ which is convex if the integrality constraints on variables $\varphi$ are relaxed.

**Proposition 2.3.** Let $\varphi \in [0, 1]^{|P|}$. The function $aTSTT(\varphi)$ is convex with regards to $\varphi$.

Using approximated TSTT values within the calculation of patterns' reduced-cost gives:

$$\hat{r}^t(\varphi) = aTSTT(\varphi) - \lambda^t - \sum_{p \in P} \mu_p^t \varphi_p, \quad \forall t \in T. \qquad (2.25)$$

If $\hat{r}^t(\varphi) \ge 0$, then there is a high likelihood that the true reduced-cost of the corresponding pattern variable $k_{\sigma}^t$ (i.e. $\varphi$ is the binary representation of $\sigma$) is positive.

To determine an initial value for $\bar{\psi}$, we determine the 5th percentile of the cumulative distribution function of random variable $\frac{TSTT(\varphi)}{e^{b_y + \sum_{p \in P} b_p \varphi_p}}$ over the observed TSTT data $\Sigma_{obs}$ and denote this value $\psi_{5th}$. This initial value for $\bar{\psi}$ provides a conservative lower bound on the TSTT values since it guarantees that the approximated TSTT is a lower bound on TSTT for 95% of the observed patterns. In order to refine this bound, we introduce an ordered set of scaling factors $\Psi$ and initialize this set with $\psi_{5th}$. The set $\Psi$ will be updated during the execution of the B&P algorithm as new TSTT data becomes available in order to provide a less conservative bound. At any point, the scaling factor $\bar{\psi}$ used in (2.24) is the average of all scaling factors in $\Psi$.

We now explain how we select a subset of the patterns and use their TSTT—obtained by solving the TAP—to build a robust statistical model. Project patterns can be characterized
by the number of active projects which follows a binomial distribution, i.e. the number of patterns with $i$ active projects is $\binom{|P|}{i}$. To span the entire spectrum of patterns, we iterate over the number of active projects $i$ and, for each $i$, we randomly select a few patterns. Since the number of patterns with $i$ active projects follows a binomial distribution, we wish to sample patterns accordingly. The ratio $\binom{|P|}{i}/2^{|P|}$ gives the proportion of patterns with $i$ active projects over the total number of patterns. We propose to pick an initial number of patterns $n_S$ with a varying density of active projects and explore their 1-distance neighborhood. Let $S_i = \left\lfloor \frac{n_S}{2^{|P|}} \right\rfloor$ be the number of randomly selected patterns with $i$ active projects and let $\Sigma_S \subset \Sigma$ be the set of randomly selected patterns, i.e. $|\Sigma_S| = \sum_{i=0..|P|} S_i$. We initialize the set of observed patterns as $\Sigma_{obs} = \Sigma_S \cup (\cup_{\sigma \in \Sigma_S} d_1(\sigma))$, where $d_1(\sigma)$ is the 1-distance neighborhood of pattern $\sigma$. We solve the TAP for each pattern in $\Sigma_{obs}$ and use the resulting TSTT values to fit Model (E). We then used the resulting TSTT values to estimate projects' first order TSTT effects based on the data gathered. Since not all TAP problems are expected to be solved during the execution of Algorithm 1, we denote $\tilde{\Delta}_p$ the vector of first order TSTT effects obtained for Project $p \in P$. The objective of this process is twofold: (i) it provides an initial set of patterns of reasonable size to build a robust statistical model; (ii) it provides data on projects' first-order TSTT effects which is shown to have a high impact on branching operations. The pseudo-code of this pre-processing step is summarized in Algorithm 1.

<table><caption>**Algorithm 1:** Pre-processing algorithm</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Input:</td><td>$N, A, P, A_p, \forall p \in P, R_a, \forall a \in A, n_S$</td></tr><tr><td>Output:</td><td>$b_y, b_p, \forall p \in P, \Delta_p, \forall p \in P, \Psi, \Sigma_{obs}$</td></tr><tr><td>1</td><td>$\Sigma_{obs} \leftarrow \emptyset$</td></tr><tr><td>2</td><td>for $i = 0..|P|$ do</td></tr><tr><td>3</td><td>&nbsp;&nbsp;&nbsp;&nbsp;$S_i \leftarrow \left\lfloor \frac{n_S}{2^{|P|}} \right\rfloor$</td></tr><tr><td>4</td><td>&nbsp;&nbsp;&nbsp;&nbsp;for 1..$S_i$ do</td></tr><tr><td>5</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\sigma \leftarrow$ randomly select a pattern with $i$ active projects</td></tr><tr><td>6</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\Sigma_{obs} \leftarrow \Sigma_{obs} \cup (\sigma \cup d_1(\sigma))$</td></tr><tr><td>7</td><td>for $\sigma \in \Sigma_{obs}$ do</td></tr><tr><td>8</td><td>&nbsp;&nbsp;&nbsp;&nbsp;$TSTT(\sigma) \leftarrow$ Solve TAP($\sigma$)</td></tr><tr><td>9</td><td>Fit Model (E) on observed TSTT values</td></tr><tr><td>10</td><td>$\tilde{\Delta}_p, \forall p \in P \leftarrow$ first-order TSTT effects</td></tr><tr><td>11</td><td>$\psi_{5th} \leftarrow$ 5th percentile of $\frac{TSTT(\varphi)}{e^{b_y + \sum_{p \in P} b_p \varphi_p}}$</td></tr><tr><td>12</td><td>$\Psi \leftarrow \{\psi_{5th}\}$</td></tr></tbody></table>

The pre-processing algorithm solves a number of TAP which grows quadratically with $|P|$. The exact number of TAP solved depends on $n_S$ and on the number of times a modified pattern has already been processed by the algorithm. The accuracy of this method is discussed in Section 2.4.7.

### 2.4.4 Optimization framework and numerical experiments

In this section, we present the optimization framework developed to solve the bilevel NMSP. This framework ultimately results in a B&P algorithm that integrates a collection of models and algorithms. B&P algorithms are used to solve mixed integer programs with a large number of variables [Barnhart et al., 1998] and have only recently been used in game-theoretic
contexts such as zero-sum Bayesian games [Halvorson et al., 2009] and Stackelberg security
games [Jain et al., 2010]. In B&P, instead of solving the LP relaxation of the integer program
at each node of an enumeration tree—as in traditional B&B frameworks—a CG algorithm is
used to solve the LRMP and iteratively add columns to the Master Problem.

Implementing a B&P algorithm involves the resolution of a pricing subproblem at each
iteration of the CG algorithm. To guarantee global optimality, this pricing subproblem should
be solved to optimality or at least a valid lower bound should be produced and passed to the
children nodes of the current node in the B&P tree. This is typically achieved by solving an
optimization problem that determines the minimal reduced cost value. Deriving such a valid
lower bound is difficult since the existence of Braess’ Paradox effects in the TAP invalidates a
traditional optimization approach aiming to minimize the reduced cost $r_{\sigma}^{t}$ of pattern variables
$k_{\sigma}^{t}$ as defined by Eq. (2.21).

A global lower bound on $r_{\sigma}^{t}$ can be derived by solving the System Optimum (SO) TAP for
the base case configuration, i.e. all projects are inactive: since the resulting SO-TSTT value
is a lower bound on the TSTT of any pattern $\sigma \in \Sigma$, this can be used to derive a valid lower
bound on the minimal reduced cost of pattern variables. However, this global lower bound is
expected to be very weak since there is often a substantial gap between the TSTT values at
SO and at UE in congested networks. Extensive numerical experiments have shown that this
bound systematically results in negative bounds on reduced costs, thus invalidating its use
for bounding purposes. To overcome this challenge, we revert to an approximation approach.

The convex approximation of TSTT values outlined in Section 2.4.3 provides a method
to determine a lower bound on the approximated TSTT of a non-observed pattern. This
bounding technique can be used in two ways: (i) to determine an approximated lower bound
at each node of the B&P tree by minimizing $aTSTT(\varphi)$ with regards to project schedule
constraints; (ii) to determine an approximated lower bound on the price of pattern variables
within the CG algorithm by minimizing $\hat{r}^{t}(\varphi)$. In both cases the optimization problem can
be solved using a convex algorithm for which efficient implementations are readily available—
recall that we relax the integrality constraints on project activeness variables.

Formally, Model 2.2 represents the Approximated NMSP (ANMSP). Since the objective
function of Model 2.2 is convex with regards to the decision variables and all constraints are
linear, the ANMSP is a convex optimization problem. Observe that the objective function
of Model 2.2 is expressed in terms of variables $y_{p}^{t}$ which indicate project activeness at each
time period. Although it is not guaranteed that $z_{ANMSP} \le z_{NMSP}$, it is highly likely due to
the presence of the scaling factor $\bar{\psi}$.
**Model 2.2 (Approximated NMSP).**

$$z_{ANMSP} = \min \sum_{t \in T} \bar{\psi} e^{b_y + \sum_{p \in P} b_p y_p^t}, \qquad (2.26)$$

$$\text{s.t.} \quad \sum_{s \in F_p} g_{s,p} = 1, \qquad \forall p \in P, \qquad (2.27)$$

$$\sum_{s \in F_p} g_{s,p} \zeta_{s,p}^t = y_p^t, \qquad \forall p \in P, \forall t \in T, \qquad (2.28)$$

$$g_{s,p} \in [0, 1], \qquad \forall p \in P, \forall s \in F_p, \qquad (2.29)$$

$$y_p^t \in [0, 1], \qquad \forall p \in P, \forall t \in T. \qquad (2.30)$$

Similarly, Model 2.3 finds a lower bound on the reduced-cost of pattern variables $k_{\sigma}^t$ by solving a separable convex optimization problem. For each time period $t \in T$, the integrality constraints in the formulation of the approximated reduced cost $\hat{r}^t(\varphi)$ are relaxed and the minimum over all time periods $\hat{r}$ is a likely lower bound on the true minimal reduced-cost of pattern variables $k_{\sigma}^t$.

**Model 2.3 (Lower bound on approximated reduced-cost).**

$$\hat{r} = \min_{t \in T} \min \left\{ \bar{\psi} e^{b_y + \sum_{p \in P} b_p \varphi_p} - \lambda^t - \sum_{p \in P} \mu_p^t \varphi_p : \varphi \in [0, 1]^{|P|} \right\}. \qquad (2.31)$$

If $\hat{r} > 0$, then the lower bound on the approximated reduced-cost is positive and the CG algorithm can be terminated with a statistical guarantee. Otherwise, this bound on the reduced-cost provides a classical lower bound on the objective value of the LRMP which may be used to prematurely terminate the CG algorithm [Lübbecke and Desrosiers, 2005]. Using constraint (2.17), the sum of all pattern variables $k_{\sigma}^t$ is: $\sum_{t \in T} \sum_{\sigma \in \Sigma_R} k_{\sigma}^t = |T|$. Let $z$ be the current value of objective function (2.11) and let $z^*$ be its optimal value:

$$\hat{z} = z + |T| \hat{r} \le z^* \le z. \qquad (2.32)$$

The approximated lower bound $\hat{z}$ provides an approximated optimality gap within the CG algorithm.

### 2.4.5 Overview of the B&P algorithm for the NMSP

An overview of the proposed B&P algorithm is summarized in Algorithm 2.

The proposed B&P algorithm starts with a pre-processing step during which the TAP is solved for a subset of the patterns selected such that projects' first-order TSTT effects can be estimated. A MILP approximating the NMSP is then solved to obtain an upper bound and initialize the set of columns (pattern variables) in the B&P framework. The algorithm then enters a **while** loop that iterates over the nodes in the enumeration tree. At first pass, the
<table><caption>Table 2.4: Network Characteristics.</caption><thead><tr><th>Network</th><th>Nb of Zones</th><th>Nb. of Nodes</th><th>Nb. of Links</th><th>Nb. of Trips</th><th>OD Flow Factor</th></tr></thead><tbody><tr><td>Sioux Falls</td><td>24</td><td>24</td><td>76</td><td>360,000</td><td>1</td></tr><tr><td>BerlinMPF</td><td>98</td><td>975</td><td>2,184</td><td>23,648.5</td><td>10</td></tr><tr><td>Barcelona</td><td>110</td><td>1,020</td><td>2,522</td><td>184,679.5</td><td>2</td></tr></tbody></table>

ANMSP is solved (Model 2.2) at the root node and if its solution is not integral, the LRMP of the NMSP is solved using CG. The best of the two lower bounds obtained is stored. If the solution of the LRMP is also not integral, a rounding heuristic is triggered and the algorithm branches on the solution providing the highest lower bound. This procedure is repeated until the enumeration tree has been fully explored or the approximated optimality gap is closed.

<table><caption>Algorithm 2: Overview of the Branch-and-Price Algorithm</caption><thead><tr><th>Step</th><th>Description</th></tr></thead><tbody><tr><td>1</td><td>Pre-processing step: get observed patterns and first-order TSTT effects</td></tr><tr><td>2</td><td>Solve MILP to obtain an upper bound for the NMSP and an initial set of columns</td></tr><tr><td>3</td><td><strong>while</strong> optimal = False <strong>do</strong></td></tr><tr><td>4</td><td>&nbsp;&nbsp;&nbsp;&nbsp;Select node in the tree</td></tr><tr><td>5</td><td>&nbsp;&nbsp;&nbsp;&nbsp;Solve ANMSP at current node</td></tr><tr><td>6</td><td>&nbsp;&nbsp;&nbsp;&nbsp;<strong>if</strong> solution is not integral <strong>then</strong></td></tr><tr><td>7</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Solve LRMP by CG at current node</td></tr><tr><td>8</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<strong>if</strong> solution is not integral <strong>then</strong></td></tr><tr><td>9</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Select solution with highest lower bound between ANMSP and CG</td></tr><tr><td>10</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Try rounding heuristic to find an improved upper bound</td></tr><tr><td>11</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Branch on solution which provides the highest lower bound</td></tr></tbody></table>

Details on each component of Algorithm 2 can be found in Rey et al. [2019].

### 2.4.6 Numerical experiments

Testing the proposed models and algorithms for the NMSP requires the following input data: a transport network, flow-dependent link travel time functions, maintenance projects with their characteristics (duration, links affected and capacity reduction) and a planning period. Given the many possibilities we focused on specific instances for the NMSP that we find representative of competitive realistic scenarios. All the results are presented for three transport networks, namely Sioux Falls, BerlinMPR and Barcelona, available at https://github.com/bstabler/TransportationNetworks which characteristics are presented in Table 2.4. The Sioux Falls is a classical benchmark network in traffic equilibrium whereas BerlinMPF and Barcelona are larger, more realistic transport networks. We use an OD flow factor on the latter two networks to increase the baseline congestion effects.

For each network, three network maintenance scenarios are used to evaluate the models and the algorithms with respectively 10, 15 and 20 projects. Maintenance projects are composed of spatially clustered links (i.e. all links in $A_p$ are connected) and, when a project is active, each affected link has its capacity reduced by 50%. The list of links affected by each maintenance project in each of the three networks considered can be found at https://github.com/davidrey123/Network_Maintenance_Scheduling.
In the 10-project scenario each project affects 100 links and this scenario induces $2^{10} = 1,024$ patterns on the network. Given that the number of pattern variables in Model 2.1 is $|T|2^{|P|}$, NMSP instances based on 10-project scenarios can be handled by most commercial MILP solvers without resorting to column generation, assuming that all 1,024 TSTTT values can be pre-processed and that $|T|$ is of reasonable size. Hence 10-project scenarios, are used to explore problem behavior and benchmark some components of the B&P algorithm that do not require any column pricing. We implicitly make the assumption that these components will perform similarly when the number of projects increases and a CG algorithm is required to manage the LP relaxation of the master problem.

In the 15- and 20-project scenarios, each project affects 50 links and these scenarios induce 32,768 and 1,048,576 patterns, respectively. Extensive testing on these scenarios indicates that solving Model 2.1 with all pattern variables poses memory management and runtime issues. Further, pre-processing all TSTTT values may be prohibitive, even for medium-size transport networks. Hence these scenarios are used to benchmark the proposed B&P algorithm. All the numerical results are presented in terms of relative TSTTT (rTSTTT) defined as the ratio of the TSTTT by the base TSTTT (when all projects are inactive).

We consider two planning periods: $|T| = 25$ and $|T| = 50$, which can be thought of as half-year-long and year-long planning periods if the time unit is assumed to represent a week. We assume that the duration of projects is normally distributed with a mean of $|T|/3$ and a standard deviation of $|T|/10$; and for all the tests conducted in the optimization framework we randomly generate NMSP instances by sampling project durations from this distribution.

The computing environment is as follows: all the models and algorithms are implemented in Python on a Windows machine with 8Gb of RAM and an Intel i7 processor at 2.9GHz. The linear and integer linear optimization problems are solved with CPLEX v12.7.1 [CPLEX, 2014] in single-threaded mode: CPLEX's primal Simplex algorithm is used for the LP relaxations of Model 2.1 and CPLEX's MIP algorithm is used for solving the proposed MILPs to find an initial feasible solution. In addition, we use IPOPT [Wächter and Biegler, 2006] to solve the convex constrained optimization problems that arise within the CG algorithm. Finally, we use the TAPAS algorithm [Bar-Gera, 2010] with an accuracy corresponding to average excess cost of $1e^{-10}$ to solve all TAPs in our tests.

### 2.4.7 Testing of the statistical regression model

We test the statistical model to approximate TSTTT values on all three networks for three sets of projects of size 10, 15 and 20 using the proposed pre-processing procedure (Algorithm 1) with a varying parameter $n_S = 1, 10$ and $100$. For each combination of network, project size and $n_S$ value, we fitted the proposed statistical model and tested its performance with an out-of-sample validation. For 10-project instances, the number of patterns is relatively small ($2^{10} = 1,024$) and using $n_S = 1$ in Algorithm 1 results in using 10.74% of the patterns. We use all remaining patterns for the out-of-sample validation of the 10-project models. Since 10.74% is already a significant proportion of patterns, higher values of $n_S$ are not tested for 10-project instances. For 15-project instances, 0.75% ($n_S = 1$) to 5.23% ($n_S = 100$) of the patterns are used, whereas for 20-project instances these figures are reduced to 0.04% ($n_S = 1$) to 0.23% ($n_S = 100$). We use 10,000 patterns for the out-of-sample validation of 15- and 20-project models. The R-square ($R^2$) and Root-Mean-Square-Error (RMSE) in terms of relative TSTTT values are summarized in Table 2.5 wherein the values in parenthesis indicate the corresponding proportion of patterns used for fitting the statistical model.
<table><caption>Table 2.5: Validation of the exponential statistical regression model (E). The R² and RMSE values are expressed in terms of relative TSTT and are obtained from an out-of-sample validation of the exponential statistical regression model (E). The patterns used for fitting the models are obtained using Algorithm 1 with different values of the tuning parameter $n_S$. The proportion of patterns used for fitting the models with regards to the entire sample space ($2^{|P|}$) is provided in parenthesis.</caption><thead><tr><th rowspan="2">|P|</th><th rowspan="2">n<sub>S</sub> (% patterns)</th><th colspan="2">Sioux Falls</th><th colspan="2">BerlinMPF</th><th colspan="2">Barcelona</th></tr><tr><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th></tr></thead><tbody><tr><td>10</td><td>1 (10.74)</td><td>0.977</td><td>0.209</td><td>0.961</td><td>0.283</td><td>0.997</td><td>0.051</td></tr><tr><td>15</td><td>1 (0.75)</td><td>0.968</td><td>0.194</td><td>0.974</td><td>0.152</td><td>0.984</td><td>0.035</td></tr><tr><td>15</td><td>10 (0.95)</td><td>0.960</td><td>0.219</td><td>0.978</td><td>0.140</td><td>0.979</td><td>0.041</td></tr><tr><td>15</td><td>100 (5.23)</td><td>0.975</td><td>0.169</td><td>0.983</td><td>0.105</td><td>0.985</td><td>0.034</td></tr><tr><td>20</td><td>1 (0.04)</td><td>0.943</td><td>0.224</td><td>0.946</td><td>0.397</td><td>0.945</td><td>0.528</td></tr><tr><td>20</td><td>10 (0.05)</td><td>0.947</td><td>0.222</td><td>0.963</td><td>0.296</td><td>0.961</td><td>0.423</td></tr><tr><td>20</td><td>100 (0.23)</td><td>0.961</td><td>0.186</td><td>0.978</td><td>0.191</td><td>0.963</td><td>0.404</td></tr></tbody></table>

The results highlight that the fitness of the proposed statistical models is of good quality with an R-square value consistently greater than 0.94, even if only a small fraction of patterns are used for fitting the models (e.g. $n_S = 1$). Increasing the number of samples through the $n_S$ parameter tends to increase R-square values and reduces RMSE values. We find that model fitness is impacted by the network and the associated maintenance projects: larger networks, in terms of number of nodes and links (such as BerlinMPF and Barcelona), appear to be more sensitive to the number of projects compared to smaller networks (Sioux Falls). Although increasing the number of projects appears to reduce R-square values and to increase RMSE values, the fitness of the statistical models remains significant. This suggests that the proposed exponential statistical regression model together with the proposed pre-processing procedure (Algorithm 1) which explore the 1-distance neighborhood of selected patterns, are appropriate to estimate TSTT effects on congested networks. Indeed the proposed exponential approximation cannot capture various effects that occur in transportation networks, such as cases of Braess' paradox. However, recent findings suggest that the existence of Braess' paradox in real transport networks may be scarce [Bekhor and Sorani, 2017], thus not diminishing the usefulness of the proposed approximation too much.

## B&P vs B&B

To show the performance of the proposed B&P algorithm, we compare the solution returned by Algorithm 2 with that obtained when solving the same problem using an exact approach for 10-project scenarios. To solve the latter, all TSTT values of the instance patterns are pre-calculated and Model 2.1 is then solved using the proposed B&B procedure. The relative gap between these two solutions are then calculated and the histograms of their distribution is plotted for each network and number of time periods. For this benchmark, 10.74% of the patterns are used for fitting Model (E). For each network and each number of time periods, we solve 50 randomly generated instances. The results are summarized in Figure 2.1. Over the 300 generated instances, we find that the highest optimality gap observed never exceeds
<figure><img src="image_2.png" /><figcaption>Figure 2.1: Benchmark of the B&B and B&P algorithms: relative gap after solving 50 randomly generated instances from three 10-project scenarios (Sioux Falls, BerlinMPF and Barcelona) with $|T| = 25,50$ and normally distributed project durations $[D_p]_{p \in P} \sim N(\frac{|T|}{3}, (\frac{|T|}{10})^2)$.</figcaption></figure>

5% and the average optimality gap is 0.57%. In addition, the proposed B&P algorithm is able to find the optimal solution in 35% of the cases. Finally, the average proportion of patterns used within the B&P algorithm is 19.9%.

This benchmark shows that the proposed B&P is capable of finding optimal or near-optimal solutions for the bilevel network maintenance scheduling problem using only a fraction of the patterns. Without accounting for the time spent solving the TAPs, the average runtime of the proposed B&B approach is 18.2s and that of the B&P is 56.9s. This gap in runtime is expected for such small instances since the number of columns remains manageable for a direct problem resolution using B&B which obviates the need for column generation. However, extensive implementation reveals that a direct resolution approach by B&B is intractable for larger instances.

## Results on large instances

We now present numerical results for the proposed B&P algorithm on instances derived from the Barcelona network. Results for the Sioux Falls and the Berlin MPF networks are available in Rey et al. [2019]. For each network, we consider four experiments consisting of instances with 15 and 20 maintenance projects and with 25 and for 50 time periods. For each experiment, we generate 20 random instances with project durations following a Normal distribution with a mean of $|T|/3$ and a standard deviation of $|T|/10$. Hence a total of 240 instances are tested regrouped in four experiments per network. The objective of these tests is to assess the performance of the proposed solution algorithm with regards to the input data.
The results are summarized in Table 2.6 which give the average performance of the B&P algorithm for each experiment (i.e. column in the tables). For each experiment, the mean value over the 20 instances is given together with the standard deviation in parenthesis. The header of the rows is described from top to bottom. The first and second rows give the number of projects and the number of patterns in each scenario. The third row gives the number of time periods in each experiment. The next nine rows provide runtime information. Since the runtime for solving TAPs may change considerably depending on the size of the network considered, we distinguish the runtime of the B&P algorithm spent solving the TAPs problems from the rest, hereby referred to as the upper-level problem and denoted U-NMSP. The next seven rows are expressed in % of the total runtime spent solving U-NMSP while the following two rows give the total runtime spent solving the U-NMSP and the TAP (lower-level problem) in hours. FEAS is the total time spent to find an initial and feasible solutions; ANMSP is the total time spent solving Model 2.2 to find lower bounds; LRMP is the total time spent solving the LRMP; PH is the total time consumed by the pricing heuristic; ARC is the total time spent solving convex approximations on reduced-costs (Model 2.3); ES is the total time spent explicitly searching for negative reduced-cost columns; RH is the total time spent in the rounding heuristic. The next six rows give statistics (in %) for the proportion of patterns used within the B&P algorithm, the proportion of TAP solved, the proportion of false negative approximations, the mean value for $\bar{\psi}$ used in the convex approximations, the optimality gap after solving the root node in the B&P tree and the proportion of time spent solving the root node in the B&P algorithm (with regards to the total solve time). The last two rows indicate the number of iterations and the value of the objective function $z_{NMSP}$ in terms of relative TSTT—defined as the ratio of the TSTT by the base TSTT—upon termination of the B&P algorithm.

We find that the mean total runtime of the upper level problem (U-NMSP) increases with the number of projects and the number of time periods but that it is robust across network topologies. Instances with 15 projects require on average 5 ($|T| = 25$) to 20 ($|T| = 50$) minutes; whereas instances with 20 projects require on average 40 ($|T| = 25$) to 90 ($|T| = 50$) minutes for the upper-level problem, i.e. solving the LRMP, the convex approximation problems and the pricing algorithms. Overall, most of the runtime for the upper-level problem is spent solving the LRMP and within the explicit search (ES) pricing algorithm. The pricing heuristic (PH) is seen to be robust to the number of time periods but less so to the number of projects. For 20-project instances, a significant amount of time is also spent exploring the 1-distance neighborhood of the patterns within the PH algorithm. The robustness with regards to the number of time periods is a direct consequence of Proposition 2.2 which obviates the exploration of all time periods to price patterns in the 1-distance neighborhood of the current solution. Although not as robust, the ES procedure is more sensitive to the number of projects than that of time periods.

The time spent solving the TAP varies considerably based on the network size and demand. The total TAP solve time is correlated with the number of TAP solved. For 15-project instances, the implementation on the three networks considered reveals that approximately 15% to 25% of the TAPs are typically solved (i.e. 5,000 to 8,000 TAPs), whereas for 20-project instances these figures change to 1.2% to 1.5% (i.e. 12,000 to 16,000 TAPs). We also find that the number of pattern variables used within the LRMP is reduced by a factor of 2 to 5 compared to the number of TAP solved. This is reflected in the high number of false negative approximations observed (often in the 60% - 80% range) and this can be attributed to the conservative lower bound which is used to initialize $\bar{\psi}$.
<table><caption>Table 2.6: Performance of the B&P Algorithm on the Barcelona network. Average performance of the B&P algorithm on the Barcelona network with 15 and 20 projects and with 25 and 50 time periods. In each experiment, 20 instances with random project durations ($([D_p]_{p \in P} \sim N(\frac{|T|}{3}, (\frac{|T|}{10})^2)$) are solved. The mean of each experiment is given with the standard deviation in parenthesis.</caption><thead><tr><th></th><th></th><th colspan="2">15</th><th colspan="2">20</th></tr><tr><th></th><th>|P|</th><th colspan="2">3e+4</th><th colspan="2">1e+6</th></tr><tr><th></th><th>|Σ|</th><th>25</th><th>50</th><th>25</th><th>50</th></tr><tr><th></th><th>|T|</th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td rowspan="7">Time (% U-NMSP)</td><td>FEAS</td><td>0.38 (0.23)</td><td>0.25 (0.45)</td><td>0.02 (0.04)</td><td>0.03 (0.04)</td></tr><tr><td>ANMSP</td><td>2.30 (1.02)</td><td>2.88 (1.32)</td><td>0.23 (0.07)</td><td>0.81 (0.16)</td></tr><tr><td>LRMP</td><td>50.45 (11.41)</td><td>76.55 (8.88)</td><td>32.54 (6.98)</td><td>38.11 (5.52)</td></tr><tr><td>PH</td><td>5.62 (1.78)</td><td>1.52 (0.49)</td><td>30.06 (8.21)</td><td>28.50 (6.87)</td></tr><tr><td>ARC</td><td>23.59 (5.78)</td><td>11.64 (4.48)</td><td>1.77 (0.58)</td><td>1.66 (0.29)</td></tr><tr><td>ES</td><td>15.07 (4.61)</td><td>6.79 (3.27)</td><td>30.18 (10.70)</td><td>25.66 (6.22)</td></tr><tr><td>RH</td><td>2.37 (9.23)</td><td>0.27 (0.46)</td><td>4.49 (14.85)</td><td>3.10 (5.18)</td></tr><tr><td rowspan="2">Time (h)</td><td>U-NMSP</td><td>0.03 (0.03)</td><td>0.15 (0.10)</td><td>0.64 (0.41)</td><td>1.07 (0.48)</td></tr><tr><td>TAP</td><td>18.25 (7.57)</td><td>18.38 (6.36)</td><td>95.99 (19.98)</td><td>103.91 (21.44)</td></tr><tr><td rowspan="8">(in %)</td><td>Patterns used</td><td>3.02 (1.17)</td><td>3.58 (1.32)</td><td>0.39 (0.09)</td><td>0.60 (0.18)</td></tr><tr><td>TAP solved</td><td>16.56 (5.70)</td><td>17.37 (4.98)</td><td>1.28 (0.21)</td><td>1.33 (0.28)</td></tr><tr><td>False negative</td><td>82.79 (4.34)</td><td>83.14 (5.16)</td><td>69.93 (7.59)</td><td>61.81 (4.08)</td></tr><tr><td>Mean ψ̃</td><td>99.91 (0.77)</td><td>100.39 (0.78)</td><td>101.48 (2.43)</td><td>99.18 (2.22)</td></tr><tr><td>Root node gap</td><td>1.49 (0.54)</td><td>1.53 (0.49)</td><td>2.36 (0.71)</td><td>5.91 (6.35)</td></tr><tr><td>Root node time</td><td>52.96 (17.72)</td><td>50.73 (24.24)</td><td>22.11 (27.28)</td><td>21.58 (22.96)</td></tr><tr><td>Nb. iterations</td><td>3.05 (1.02)</td><td>3.35 (1.01)</td><td>4.10 (1.58)</td><td>3.45 (1.20)</td></tr><tr><td>z<sub>NMSP</sub></td><td>32.64 (1.20)</td><td>65.16 (2.69)</td><td>42.73 (4.19)</td><td>92.88 (7.30)</td></tr></tbody></table>

The root node gap is on average relatively tight, with low values of 1.5% observed for the 15-project instances on the Barcelona network and high values of 5.7 to 6.2% observed for 20-project instances on all three networks tested. The time spent in the B&P algorithm at the root node is often substantial and varies between 18% and 61% of the total solve time for the Sioux Falls network. Since this accounts for the time spent in TAPs, this proportion changes considerably based on the network. For the Barcelona network, on average, half of the total solve time in 15-project instances is spent at the root node whereas this represents only 22% of the total solve time in 20-project. This is due to the greater time required to solve TAPs on the Barcelona network compared to the smaller networks tested.

We find that the number of iterations required by the B&P algorithm only slightly increases with the number of projects and time periods, with an average of 3 to 5 iterations frequently achieved. This low number of iterations together with the analysis of the B&P algorithm at the root node suggests that most of the progress is performed in the CG process at the root node. Subsequent iterations then attempt to close the gap by branching but may also offer the possibility to generate additional patterns and further improve the solution.
## 2.5 Scheduling for road network disaster recovery

The scope of this section is the long-term problem of planning the reconstruction of a road network which was damaged by a disaster. The recovery of the entire road network is assumed to require a substantial amount of resources, e.g. budget, working crew, equipment; which are limited in availability. We assume that the recovery of the road network will take place over a pre-defined planning horizon (e.g. 12-24 months), which is divided into time periods (e.g. 2-4 weeks). We assume that the time periods are long enough so that travelers have the time to adjust their route choice preferences under the UE principle. We assume that the reconstruction of the road network can be organized into recovery projects, each of which represents a set of roads damaged by the disaster. The network recovery problem of interest is that of optimizing the order in which the recovery projects should be completed subject to resource availability constraints and assuming UE conditions. We build on and extend the pattern formulation of the NMSP presented in Section 2.4.1 to account for network recovery effects, and introduce customized scheduling heuristics.

### 2.5.1 Formulation of the network recovery problem

We use the same notation system as in Section 2.4.1. We extend this framework to characterize the evolution of links of the network beyond two states (which is the case for the NMSP). Let $\gamma_{p,m,\sigma}$ be a binary parameter indicating if project $p$ is in state $m$ in pattern $\sigma$ (1) or not (0):

$$ \gamma_{p,m,\sigma} = \begin{cases} 1 & \text{if project } p \text{ is in state } m \text{ in pattern } \sigma, \\ 0 & \text{otherwise,} \end{cases} \quad \forall p \in P, \forall m \in M, \forall \sigma \in \Sigma. \quad (2.33) $$

Let $\omega_{p,m}^{s,t}$ be a binary parameter indicating if project $p$ is in state $m$ at time period $t$, if $p$ started at time period $s$ (1) or not (0):

$$ \omega_{p,m}^{s,t} = \begin{cases} 1 & \text{if project } p \text{ is in state } m \text{ at time period } t, \\ 0 & \text{if } p \text{ started at time period } s, \\ 0 & \text{otherwise,} \end{cases} \quad \forall p \in P, \forall m \in M, \forall s \in F_p, \forall t \in T. \quad (2.34) $$

For each project $p$, for each state $m$ and for each time period $t$, the following linking constraints link variables $g_{s,p}$ and $k_{\sigma}^t$:

$$ \sum_{\sigma \in \Sigma} k_{\sigma}^{t} \gamma_{p,m,\sigma} = \sum_{s \in F_p} g_{s,p} \omega_{p,m}^{s,t}, \quad \forall p \in P, \forall m \in M, \forall t \in T. \quad (2.35) $$

Constraints (2.35) combine two representations of the concept of project activeness, i.e. both its left-hand-side and its right-hand-side take value 1 if project $p$ is in state $m$ at time $t$, and 0 otherwise. To model recovery resources availability, e.g. work crew or equipment, we assume that each project $p \in P$ in state $m \in M$ requires a known amount of recovery resources per time period denoted $R_{p,m}$. The limited availability in recovery resources can be represented by the total amount of recovery resources available at time period $t \in T$ denoted $R_t$. Thus, we impose the following resource availability constraint:

$$ \sum_{p \in P} \sum_{m \in M} \sum_{s \in F_p} g_{s,p} \omega_{p,m}^{s,t} R_{p,m} \leq R_t, \quad \forall t \in T. \quad (2.36) $$
The objective function of the leader problem is to minimize the total network travel time over the planning horizon. For each project pattern $\sigma$ there exists a unique corresponding network link state vector $\mathbf{m}$. Formally, if pattern $\sigma$ is selected at time $t$, i.e. $k_{\sigma}^t = 1$, then the set of links $\{a \in A_p : p \in P\}$ are in state $m$ if $\gamma_{p,m,\sigma} = 1$.

We can now present the complete formulation of the proposed bilevel network recovery problem. The leader problem is summarized in (2.37). The leader variables are project start times $g_{s,p}$ and pattern selection variables $k_{\sigma}^t$. The follower problem is only implicitly represented by the network travel time $TSTT(\sigma)(\sigma)$ in the objective function of the leader.

$$ \min \quad \sum_{t \in T} \sum_{\sigma \in \Sigma} k_{\sigma}^t TSTT(\sigma)(\sigma), \qquad (2.37a) $$

$$ \text{s.t.} \quad \sum_{s \in F_p} g_{s,p} = 1, \qquad \forall p \in P, \qquad (2.37b) $$

$$ \sum_{\sigma \in \Sigma} k_{\sigma}^t = 1, \qquad \forall t \in T, \qquad (2.37c) $$

$$ \sum_{\sigma \in \Sigma} k_{\sigma}^t \gamma_{p,m,\sigma} = \sum_{s \in F_p} g_{s,p} \omega_{p,m}^{s,t}, \qquad \forall p \in P, \forall m \in M, \forall t \in T, \qquad (2.37d) $$

$$ \sum_{p \in P} \sum_{m \in M} \sum_{s \in F_p} g_{s,p} \omega_{p,m}^{s,t} R_{p,m} \leq R_t, \qquad \forall t \in T, \qquad (2.37e) $$

$$ g_{s,p} \in \{0, 1\}, \qquad \forall p \in P, s \in F_p, \qquad (2.37f) $$

$$ k_{\sigma}^t \in \{0, 1\}, \qquad \forall t \in T, \sigma \in \Sigma. \qquad (2.37g) $$

We detail now how the required model parameters, i.e. binary matrices $\gamma_{p,m,\sigma}$ and $\omega_{p,m}^{s,t}$, and recovery resources $R_{p,m}$, can be determined based on the number of states $|M|$. We consider a three-states case with the states damaged, repair and restored. In damaged state, links are deteriorated and perform accordingly. Then, links evolve to the repair state for $D_p$ time periods. Finally, links evolve to the restored state in which they recover their original link performance characteristics. In this three-states case, we assume that $R_{p,damaged} = R_{p,restored} = 0$ and that $R_{p,repair} > 0$, thus ensuring that recovery resources are only consumed during the $D_p$ time periods when $p$ is active. The matrix $\omega_{p,m}^{s,t}$ can then be determined by partitioning the time line in three segments corresponding to the three states of the project lifecycle.

$$ \omega_{p,damaged}^{s,t} \equiv \begin{cases} 1 & \text{if } 0 \le t \le s-1, \\ 0 & \text{otherwise,} \end{cases} \quad \forall p \in P, \forall s \in F_p, \forall t \in T. \qquad (2.38) $$

$$ \omega_{p,repair}^{s,t} \equiv \begin{cases} 1 & \text{if } s \le t \le s + D_p - 1, \\ 0 & \text{otherwise,} \end{cases} \quad \forall p \in P, \forall s \in F_p, \forall t \in T. \qquad (2.39) $$

$$ \omega_{p,restored}^{s,t} \equiv \begin{cases} 1 & \text{if } s + D_p \le t \le |T|, \\ 0 & \text{otherwise,} \end{cases} \quad \forall p \in P, \forall s \in F_p, \forall t \in T. \qquad (2.40) $$

The choice regarding the number of states used in the model depends on the ability to estimate the impacts during each state at the time that the scheduling of recovery projects is determined. In our numerical experiments, we consider a three-state case where $M =$
{damaged, repair, restored} but due to the lack of data on link performance characteristics post-disaster and during reconstruction, we assume that links in either state damaged or repair behave identically. Thus, in our numerical experiments, although $|M| = 3$, only $2^{|P|}$ project patterns are used to model the network recovery problem.

## 2.5.2 Scheduling heuristics and numerical experiments

We propose three heuristic algorithms for the bilevel network recovery problem that build on existing methods in the scheduling literature [Smith, 1956]. All three heuristics are greedy algorithms that differ in the criterion used to sort recovery projects. Once sorted, recovery projects are scheduled as soon as possible while respecting the resource availability constraint (2.36). It should be noted that these heuristics are not guaranteed to find a feasible solution if one exists.

### Shortest Processing Time (SPT)

A natural scheduling heuristic is to sort recovery projects by increasing duration $D_p$. This is known as the Shortest Processing Time (SPT) rule [Smith, 1956]. Using the SPT rule is intuitive in that by completing shortest projects first, the benefits of recovering the associated nominal link characteristics are expected to yield network-wide improvements early in the planning horizon. This heuristic does not require any knowledge of the congestion effects associated with damaged links and their reconstruction.

### Largest Average First-Order (LAFO)

To incorporate congestion effects within the scheduling process, we can estimate the average first-order effects of each recovery project. The first-order effects of a recovery project represent the network travel time reduction observed before and after this project is completed. These effects depend on the completion of other recovery projects, hence calculating exactly the average first-order effect of a recovery project requires comparing network travel time for all $|M|^{P}-1$ combinations of other recovery projects. We denote $\Delta_p$ the distribution of first-order effects for project $p \in P$ and we denote $\bar{\Delta}_p$ the average first-order effect of $p$. The Largest Average First-Order (LAFO) heuristic consists of sorting projects by decreasing average first-order effects $\bar{\Delta}_p$. This heuristic is expected to complete the most impactful (in terms of network travel time reduction) projects first, thus benefiting the entire network for subsequent time periods.

### Largest Approximated Smith's Ratio (LASR)

Smith's ratio, also known as weighted-SPT, consists of sorting projects by decreasing "weight"-to-duration ratio [Smith, 1956]. This rule is optimal for the problem of minimizing the total weighted completion time in a single-machine scheduling problem. In the context of the bilevel network recovery scheduling problem, the "weight" of each project is unknown a priori due to the nonlinear congestion effects inherent to traffic equilibrium. A possible adaptation of Smith's ratio consists of using projects' average first-order network travel time effects to approximate the "weight" of each project. Hence, we sort projects by decreasing $\bar{\Delta}_p/D_p$ values.
The LASR heuristic requires as much computational effort as the LAFO heuristic but captures both projects impact and duration in the decision process. Both LAFO and LASR heuristics require the calculation of project's average first-order effects, which require the solution of $|M|^{|P|-1}$ TAPs. To avoid the possibly resource-intensive computational effort, heuristics algorithms can be used to approximate the exact average first-order effects by only solving a subset of the required TAPs. We use Algorithm 1 to sample subsets of the $|M|^{|P|-1}$ combinations and estimate first-order effects. In our numerical experiments, we compare the behavior of the LAFO and LASR heuristics where all project patterns are used to determine the average first-order effects; and the more practical case where only a subset of project patterns is used for this calculation.

We next report numerical results for the proposed formulation and solution methods. All numerical experiments are conducted on a desktop computer with 16Gb of RAM and a CPU of 3.6GHz. The optimal solution of the bilevel network recovery scheduling problem is obtained by solving the leader problem (2.37) using Mixed-Integer Linear Programming (MILP) commercial software after having solved all $|M|^{|P|}$ TAPs using TAPAS [Bar-Gera, 2010]. The MILPs are solved with CPLEX 12.8 Python API.

### 2.5.3 Numerical results

We conduct numerical experiments on two types of disaster scenario: a **corridor scenario** where all links of a corridor of the network are damaged, and a **crossing scenario** where all links of cross-section of the network are damaged. The motivation behind both disaster scenarios is that a natural disaster such as an earthquake may damage the network within a specific topological region. In the proposed corridor and crossing scenarios, the topology of the damaged region is assumed to follow a fault-line, i.e. longitudinal, pattern. We use the Berlin Mitte Centre network which contains 36 centroids, 398 nodes and 871 links to benchmark the proposed heuristics. This network is available at https://github.com/bstabler/TransportationNetworks.

For the corridor scenario, damaged links are organized into 10 recovery projects and all damaged links are assumed to be affected by a 95% reduction in link capacity. For the crossing scenario, damaged links are organized into 9 recovery projects and all damaged links are assumed to be affected by a 90% reduction in link capacity. In both scenarios, the free-flow travel times of damaged links are ten times greater than their default values (additional link characteristics are provided with the transport network data). Two maps depicting the damaged and recovery projects in the corridor and crossing scenarios are provided in Figure 2.2.

The planning horizon is set to 52 time periods (in weeks) and we assume that project durations are uniformly distributed in the range {1, ..., 26}, i.e. the maximum project duration is half of the planning horizon. We assume that all projects require a unit resource per time period, i.e. $R_p = 1$ for all $p \in P$. The available recovery resources $R_t$ are assumed to remain constant across time periods.

For both corridor and crossing scenarios, we randomly generate 10 instances (total of 20 instances). For each instance, we vary the available amount of recovery resources $R_t$ from 3 to 9 (corridor) or 8 (crossing). Note that in this case, $R_t$ represents the maximum number of projects that can be *under repair* in parallel. The lower bound 3 was determined by observing that the problem is often infeasible if $R_t < 3, \forall t \in T$.
<figure><img src="image_3.png" /><figcaption>Figure 2.2: Map of the Berlin Mitte Centre network with colored recovery projects for the corridor and crossing scenarios.</figcaption></figure>

TAPs on the Berlin Mitte Centre network solved with TAPAS require on average 1s per TAP, hence the pre-processing time required to implement the MILP and to exactly calculate projects' first-order effects is on average $2^{10} = 1,024s$ for the Corridor scenario, and $2^9 = 512s$ for the Crossing scenario. In contrast, implementing heuristics LAFO and LASR using the sampling algorithm of Rey et al. [2019] to approximate project's average first-order effects requires the solution of 113 and 88 TAPs, for the Corridor and Crossing scenarios, respectively. Further, in all our numerical experiments, the run time of the heuristics (without pre-processing) is significantly less than 0.1 s. This can be expected, since the heuristics essentially consists of sorting a list of $|P|$ values by a provided criterion (SPT, LAFO or LASR) and creating a schedule based on this sorted list.

All numerical results are presented in terms of normalized network travel time. Let $TSTT(\sigma)(\sigma_0)$ be the network travel time corresponding to the pattern $\sigma_0$ wherein all recovery projects have been completed, i.e. $\gamma_{p,restored,\sigma_0} = 1$, $\gamma_{p,damaged,\sigma_0} = \gamma_{p,repair,\sigma_0} = 0$, $\forall p \in P$. The normalized network travel time of pattern $\sigma$ is $TSTT(\sigma)(\sigma) = \frac{TSTT(\sigma)(\sigma)}{TSTT(\sigma)(\sigma_0)}$.

The average performance of the heuristics is summarized in Table 2.7. For each type of disaster (Corridor and Crossing) and for each level of recovery resources ($R_t$), we report the minimum (min), average (avg) and maximum (max) optimality gap over all 10 instances of this group. We also report how frequently each heuristic provides the best optimality gap (best), and how frequently each heuristic is unable to produce a feasible solution (inf). Average performance over all levels of recovery resources is reported in the last three rows.

LAFO-all is the heuristic which most frequently provides the best optimality gap for Corridor instances: this heuristic dominates other heuristics 62.86% of the time. For Crossing instances, LASR-all dominates 63.33% of them time. LAFO-all is the heuristic less frequently infeasible in Crossing scenario instances, whereas LASR-all outperforms other heuristics in this regard on Corridor scenario instances. However, we find that, on average, LASR is the most robust heuristic and achieves an overall average optimality gap of 0.52% for the Corridor
<table><caption>Table 2.7: Summary of the performance of heuristics SPT, LAFO and LASR. All data is in percentage. For each type of disaster scenario (Corridor and Crossing) and each level of recovery resources ($R_t$) the minimum (min), average (avg) and maximum (max) optimality gap over all 10 instances of this group is reported. Columns “best” indicate the frequency that this heuristic was the best (least optimality gap), and columns “inf” indicate the frequency that this heuristic did not produce a feasible solution. The bottom five rows indicate average performance over all values of recovery resources $R_t$.</caption><thead><tr><th rowspan="2">R<sub>t</sub></th><th rowspan="2">Heuristic</th><th rowspan="2">Type</th><th colspan="5">Corridor</th><th colspan="5">Crossing</th></tr><tr><th>min</th><th>avg</th><th>max</th><th>best</th><th>inf</th><th>min</th><th>avg</th><th>max</th><th>best</th><th>inf</th></tr></thead><tbody><tr><td rowspan="3">3</td><td>SPT</td><td></td><td>1.51</td><td>81.84</td><td>100.00</td><td>60.00</td><td>80.00</td><td>0.16</td><td>41.16</td><td>100.00</td><td>50.00</td><td>30.00</td></tr><tr><td>LAFO</td><td>All</td><td>0.71</td><td>76.02</td><td>100.00</td><td>60.00</td><td>70.00</td><td>0.83</td><td>59.17</td><td>100.00</td><td>30.00</td><td>20.00</td></tr><tr><td></td><td>(sampling)</td><td>(0.85)</td><td>(81.13)</td><td>(100.00)</td><td>(80.00)</td><td>(80.00)</td><td>(1.23)</td><td>(26.74)</td><td>(100.00)</td><td>(90.00)</td><td>(10.00)</td></tr><tr><td rowspan="3">4</td><td>LASR</td><td>All</td><td>0.43</td><td>60.38</td><td>100.00</td><td>100.00</td><td>60.00</td><td>0.22</td><td>30.61</td><td>100.00</td><td>60.00</td><td>30.00</td></tr><tr><td></td><td>(sampling)</td><td>(0.06)</td><td>(60.29)</td><td>(100.00)</td><td>(70.00)</td><td>(60.00)</td><td>(0.65)</td><td>(38.21)</td><td>(100.00)</td><td>(50.00)</td><td>(30.00)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">5</td><td>SPT</td><td></td><td>0.33</td><td>2.75</td><td>9.65</td><td>30.00</td><td>0.00</td><td>0.11</td><td>8.70</td><td>28.07</td><td>30.00</td><td>0.00</td></tr><tr><td>LAFO</td><td>All</td><td>0.16</td><td>12.33</td><td>36.48</td><td>30.00</td><td>0.00</td><td>0.36</td><td>15.21</td><td>52.42</td><td>0.00</td><td>0.00</td></tr><tr><td></td><td>(sampling)</td><td>(0.07)</td><td>(12.31)</td><td>(36.39)</td><td>(60.00)</td><td>(0.00)</td><td>(0.84)</td><td>(6.83)</td><td>(18.98)</td><td>(50.00)</td><td>(0.00)</td></tr><tr><td rowspan="3">6</td><td>LASR</td><td>All</td><td>0.22</td><td>1.02</td><td>3.23</td><td>40.00</td><td>0.00</td><td>0.00</td><td>0.46</td><td>0.87</td><td>70.00</td><td>0.00</td></tr><tr><td></td><td>(sampling)</td><td>(0.07)</td><td>(1.33)</td><td>(8.00)</td><td>(50.00)</td><td>(0.00)</td><td>(0.32)</td><td>(4.96)</td><td>(24.20)</td><td>(10.00)</td><td>(0.00)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">7</td><td>SPT</td><td></td><td>0.08</td><td>1.68</td><td>5.78</td><td>20.00</td><td>0.00</td><td>0.00</td><td>3.90</td><td>20.13</td><td>50.00</td><td>0.00</td></tr><tr><td>LAFO</td><td>All</td><td>0.26</td><td>1.35</td><td>8.36</td><td>60.00</td><td>0.00</td><td>0.17</td><td>3.40</td><td>30.32</td><td>30.00</td><td>0.00</td></tr><tr><td></td><td>(sampling)</td><td>(0.14)</td><td>(1.32)</td><td>(8.49)</td><td>(60.00)</td><td>(0.00)</td><td>(0.24)</td><td>(0.70)</td><td>(2.07)</td><td>(30.00)</td><td>(0.00)</td></tr><tr><td rowspan="3">8</td><td>LASR</td><td>All</td><td>0.09</td><td>0.57</td><td>1.36</td><td>40.00</td><td>0.00</td><td>0.00</td><td>0.35</td><td>0.73</td><td>40.00</td><td>0.00</td></tr><tr><td></td><td>(sampling)</td><td>(0.14)</td><td>(1.36)</td><td>(8.36)</td><td>(50.00)</td><td>(0.00)</td><td>(0.03)</td><td>(3.52)</td><td>(16.16)</td><td>(50.00)</td><td>(0.00)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">9</td><td>SPT</td><td></td><td>0.17</td><td>0.74</td><td>1.72</td><td>20.00</td><td>0.00</td><td>0.00</td><td>1.79</td><td>16.00</td><td>50.00</td><td>0.00</td></tr><tr><td>LAFO</td><td>All</td><td>0.10</td><td>0.82</td><td>3.85</td><td>70.00</td><td>0.00</td><td>0.00</td><td>3.11</td><td>29.85</td><td>50.00</td><td>0.00</td></tr><tr><td></td><td>(sampling)</td><td>(0.14)</td><td>(0.82)</td><td>(4.18)</td><td>(40.00)</td><td>(0.00)</td><td>(0.09)</td><td>(0.28)</td><td>(0.75)</td><td>(30.00)</td><td>(0.00)</td></tr><tr><td rowspan="3">10</td><td>LASR</td><td>All</td><td>0.10</td><td>0.57</td><td>1.54</td><td>10.00</td><td>0.00</td><td>0.00</td><td>0.17</td><td>0.39</td><td>40.00</td><td>0.00</td></tr><tr><td></td><td>(sampling)</td><td>(0.14)</td><td>(0.85)</td><td>(3.95)</td><td>(60.00)</td><td>(0.00)</td><td>(0.01)</td><td>(0.25)</td><td>(0.57)</td><td>(50.00)</td><td>(0.00)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">11</td><td>SPT</td><td></td><td>0.12</td><td>0.35</td><td>0.82</td><td>30.00</td><td>0.00</td><td>0.00</td><td>0.13</td><td>0.45</td><td>20.00</td><td>0.00</td></tr><tr><td>LAFO</td><td>All</td><td>0.05</td><td>0.43</td><td>1.34</td><td>60.00</td><td>0.00</td><td>0.00</td><td>0.06</td><td>0.36</td><td>60.00</td><td>0.00</td></tr><tr><td></td><td>(sampling)</td><td>(0.15)</td><td>(0.82)</td><td>(3.94)</td><td>(40.00)</td><td>(0.00)</td><td>(0.00)</td><td>(0.06)</td><td>(0.36)</td><td>(100.00)</td><td>(0.00)</td></tr><tr><td rowspan="3">12</td><td>LASR</td><td>All</td><td>0.06</td><td>0.42</td><td>1.34</td><td>20.00</td><td>0.00</td><td>0.00</td><td>0.01</td><td>0.10</td><td>90.00</td><td>0.00</td></tr><tr><td></td><td>(sampling)</td><td>(0.00)</td><td>(0.41)</td><td>(1.40)</td><td>(70.00)</td><td>(0.00)</td><td>(0.00)</td><td>(0.13)</td><td>(0.36)</td><td>(60.00)</td><td>(0.00)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">13</td><td>SPT</td><td></td><td>0.08</td><td>0.23</td><td>0.61</td><td>30.00</td><td>0.00</td><td>0.00</td><td>0.03</td><td>0.17</td><td>30.00</td><td>0.00</td></tr><tr><td>LAFO</td><td>All</td><td>0.02</td><td>0.27</td><td>1.08</td><td>70.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.01</td><td>80.00</td><td>0.00</td></tr><tr><td></td><td>(sampling)</td><td>(0.00)</td><td>(0.48)</td><td>(3.74)</td><td>(70.00)</td><td>(0.00)</td><td>(0.00)</td><td>(0.00)</td><td>(0.01)</td><td>(100.00)</td><td>(0.00)</td></tr><tr><td rowspan="3">14</td><td>LASR</td><td>All</td><td>0.07</td><td>0.31</td><td>1.08</td><td>20.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.01</td><td>80.00</td><td>0.00</td></tr><tr><td></td><td>(sampling)</td><td>(0.00)</td><td>(0.21)</td><td>(0.90)</td><td>(90.00)</td><td>(0.00)</td><td>(0.00)</td><td>(0.00)</td><td>(0.01)</td><td>(100.00)</td><td>(0.00)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">15</td><td>SPT</td><td></td><td>0.02</td><td>0.11</td><td>0.18</td><td>40.00</td><td>0.00</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>LAFO</td><td>All</td><td>0.00</td><td>0.06</td><td>0.28</td><td>90.00</td><td>0.00</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td>(sampling)</td><td>(0.00)</td><td>(0.06)</td><td>(0.28)</td><td>(100.00)</td><td>(0.00)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="3">16</td><td>LASR</td><td>All</td><td>0.00</td><td>0.07</td><td>0.28</td><td>80.00</td><td>0.00</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td>(sampling)</td><td>(0.00)</td><td>(0.07)</td><td>(0.28)</td><td>(90.00)</td><td>(0.00)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">17</td><td>SPT</td><td></td><td><strong>0.02</strong></td><td><strong>1.24</strong></td><td><strong>16.86</strong></td><td><strong>32.86</strong></td><td><strong>11.43</strong></td><td><strong>0.00</strong></td><td><strong>4.51</strong></td><td><strong>43.91</strong></td><td><strong>38.33</strong></td><td><strong>5.00</strong></td></tr><tr><td>LAFO</td><td>All</td><td><strong>0.00</strong></td><td><strong>3.38</strong></td><td><strong>36.48</strong></td><td><strong>62.86</strong></td><td><strong>10.00</strong></td><td><strong>0.00</strong></td><td><strong>10.51</strong></td><td><strong>82.11</strong></td><td><strong>41.67</strong></td><td><strong>3.33</strong></td></tr><tr><td></td><td>(sampling)</td><td>(0.00)</td><td>(2.73)</td><td>(36.39)</td><td>(64.29)</td><td>(11.43)</td><td>(0.00)</td><td>(4.17)</td><td>(65.87)</td><td>(66.67)</td><td>(1.67)</td></tr><tr><td rowspan="3">18</td><td>LASR</td><td>All</td><td><strong>0.00</strong></td><td><strong>0.52</strong></td><td><strong>3.23</strong></td><td><strong>44.29</strong></td><td><strong>8.57</strong></td><td><strong>0.00</strong></td><td><strong>0.28</strong></td><td><strong>1.86</strong></td><td><strong>63.33</strong></td><td><strong>5.00</strong></td></tr><tr><td></td><td>(sampling)</td><td>(0.00)</td><td>(0.71)</td><td>(8.36)</td><td>(68.57)</td><td>(8.57)</td><td>(0.00)</td><td>(3.00)</td><td>(28.21)</td><td>(53.33)</td><td>(5.00)</td></tr></tbody></table>
scenario instances, and of 0.28% for the Crossing scenario instances, followed by SPT. Hence while LAFO may frequently perform best, it also is more prone to poor performance. This is reflected in the average maximum optimality gaps observed which are substantially higher for this heuristic compared to that of SPT and LASR. If only a subset of patterns are used to approximate project's first-order effects, we find that the maximum optimality gaps increase but average performance is marginally impacted for LASR. In contrast, we find that LAFO on average benefits from this approximation, which suggests while first-order effects are critical as highlighted by the good performance of LASR, their calculation and incorporation in the proposed heuristics require further research.

## 2.6 Conclusion

The DNDP is a challenging bilevel optimization problem with critical implications for network design in urban transport systems. Over the past decade, exact methods have become more capable to solve non-trivial decisions problems. We synthesized the literature on exact methodologies and discussed the main approaches developed to solve the DNDP to optimality. In total, nine papers have been examined. Although some of these papers have not addressed the DNDP in its present form, they have presented a methodology for a related problem which can be reduced to the DNDP. The characteristics of the solution methods examined have been categorized in two types of approaches: SO-relaxation based and KKT conditions based methods. A computational benchmarking of the DNDP was conducted by adapting three solution methods and comparing their performance using the same piecewise linear link travel time function approximation. For this benchmark, a total of 20 new instances for the DNDP have been created. For reproducibility purposes, all implementation formulations, codes and benchmarking instances are available at https://github.com/davidrey123/DNDP.

The problem of scheduling maintenance projects in congested transport networks can be viewed as a multi-period application of the DNDP. TSTT can be used as a metric to represent the impact of road capacity reductions onto network performance by accounting for users' reaction in the schedule optimization. A novel pattern-based formulation has been proposed which extends previous work by providing a compact MIP representation based on the enumeration of project patterns. An exponential regression model was introduced to approximate network congestion effects based on maintenance projects' first order TSTT effects. This model is used in combination with a pre-processing algorithm to explore the project pattern space based on the concept of 1-distance neighborhood which proved to be very effective in capturing projects' impact on network congestion. We showed that the integer-relaxation of the proposed statistical model is convex and that this relaxation provides an efficient method to approximate TSTT values. MILP models to find feasible solutions as well as customized branching and node selection rules are proposed and calibrated. A CG algorithm was then presented to solve the NMSP on large instances where pre-processing all TSTT values is intractable. The proposed CG algorithm relies on a fast pricing heuristic based on 1-distance optimality, approximate lower bounds derived from convex approximations of the NMSP and an explicit search procedure. The proposed B&P algorithm was tested on instances derived from 15- and 20-project scenarios on three network topologies with 25 and 50 time periods.

Our implementation demonstrated that the proposed B&P algorithm is able to solve all instances in reasonable time and by using only a fraction of the project patterns during the exploration. This is an important result since it suggests that competitive solutions can be
found efficiently even on large instances. From a runtime perspective, the bottleneck of the
algorithm is the time spent solving the TAP and the time spent in the explicit search pro-
cedure. While the former depends on the network size, the latter depends on the number
of projects considered. Both of these steps may be challenging from a computational stand-
point and additional work is necessary to identify more efficient approaches to manage this
complexity.

The pattern-based formulation for the NMSP was then extended to incorporate network recovery effects. Three heuristics employed from the scheduling literature were adapted for this bilevel network recovery problem. Numerical experiments were designed on a realistic transport network wherein two disaster scenarios were considered. Both disaster scenarios aim to represent the possible damage caused by disasters that exhibit a fault-line topology, such as earthquakes. A total of 20 instances based on these disaster scenarios were randomly generated and used to benchmark the proposed solution methods. The numerical results obtained highlight the performance of the proposed adaptation of Smith's ratio for the bilevel network recovery scheduling problem at hand, but also reveal that the proposed heuristic methods may fall short of optimal solutions when the amount of recovery resources is low compared to the demand of recovery projects. This research provides a basis to support decision-making in the context of disaster recovery. The proposed formulation focused on the reconstruction of the road network under limited recovery resources availability which can help in informing planning policies.
# Chapter 3

## Aircraft conflict resolution via mixed-integer nonlinear programming

This chapter is based on Rey and Hijazi [2017], Dias et al. [2022], Dias and Rey [2022].

### 3.1 Introduction

Air traffic control (ATC) is an extremely dynamic and constrained environment where many decisions need to be taken in a short amount of time. Adopting automation within such an environment can be vital to reduce controller workload and improve airspace capacity [Durand et al., 1997, Barnier and Allignol, 2009, Rey et al., 2016c]. Traditional methods for air traffic control have been exhaustively used and are reaching their limits, hence automated approaches are receiving a significant and growing attention in the field [Vela et al., 2009b]. The International Civil Aviation Organization (ICAO) determines all regulations related to civil aviation [ICAO, 2010]. One of its main roles is to set separation standards for commercial aviation. We focus on aircraft separation for en-route traffic. During cruise stage, separation conditions require a minimum of 5 Nautical Miles (NM) horizontally or 1000 feet (ft) vertically between any pair of aircraft. A conflict between two or more aircraft is a loss of separation among these aircraft. Air traffic networks are organised in flight levels which are separated by at least 1000 ft, hence during cruise stage, most conflicts occur among aircraft flying at the same flight level. Congested air traffic networks can lead to loss of separation between aircraft which impairs flight safety and may result in collisions. The aircraft conflict resolution problem (ACRP) can be formulated as an optimization problem in which the objective is to find least-deviating conflict-free trajectories for a set of aircraft. Different strategies have been used to address this problem based on the type of deconfliction manoeuvres available, namely: speed control (acceleration or deceleration), heading control, vertical control (flight level reassignment) or a combination of these manoeuvres.

This chapter presents the complex number formulation for aircraft conflict resolution introduced by Rey and Hijazi [2017] and further developed by Dias et al. [2022] for the deterministic case, and by Dias and Rey [2022] for the robust case. Nonconvex mixed-integer programming (MIP) formulations for ACRPs and exact solution methods are proposed
for the two-dimensional (2D) conflict resolution problem with speed and heading control represented as continuous variables; as well as for an extended problem wherein flight level (FL) changes are also available. Theoretical and numerical contributions are provided, and extensive numerical results on benchmarking instances are reported. For reproducibility purposes, all formulations and instances are made available online at the public repository https://github.com/acrp-lib/acrp-lib.

The chapter is organised as follows. We review the state-of-the-art on aircraft conflict resolution, highlight existing research gaps and position our contributions to the field in Section 3.2. We next characterise pairwise aircraft separation conditions and present nonconvex aircraft conflict resolution formulations for the 2D problem, as well as for an extension of the 2D problem with FL changes, hereby referred to as 2D+FL in Section 3.3. We present exact solution methods for the 2D and 2D+FL problems that build on existing convex relaxations and propose novel decomposition and constraint generation algorithms in Section 3.4. The robust aircraft conflict resolution problem is discussed in Section 3.5. Numerical results are provided in Section 3.6 and concluding remarks are discussed in Section 3.7.

## 3.2 Literature review

The first exact global optimization approach for aircraft conflict resolution problems was proposed by Pallottino et al. [2002], who introduced two MIP formulations: a first model based on speed control only and a second model based on heading control only which required that all aircraft fly at the same speed. In the proposed MIP formulation for conflict resolution with speed control, the authors derived linear pairwise aircraft separation constraints based on the geometric construction introduced by Bilimoria [2000]. These separation conditions are obtained by projecting the shadow of an intruder aircraft onto the trajectory of a reference aircraft, thus we henceforth refer to these separation conditions as the shadow separation conditions. Frazzoli et al. [2001] were the first to observe that this geometric construction provided a basis to characterise the set of aircraft pairwise conflict-free trajectories via linear half-planes in the relative velocity (speed and heading) plane. The authors introduced a nonconvex formulation for the conflict resolution problem with speed and heading control, and proposed a convex relaxation based on semi-definite programming as well as a heuristic algorithm to find feasible solutions on problems with up to 10 aircraft.

The shadow separation conditions were subsequently used in several formulations. Alonso-Ayuso et al. [2011] proposed a mixed-integer linear programming (MILP) formulation for conflict resolution by speed and altitude control and reported solving instances with up to 50 aircraft in competitive time. Alonso-Ayuso et al. [2014] proposed a two-step approach in which only heading control is available for deconfliction and the available angle changes are discretised. The same group of authors also proposed a nonconvex formulation involving trigonometric functions based on the shadow separation conditions that enables speed, heading and altitude control [Alonso-Ayuso et al., 2016]. The authors used a mixed-integer nonlinear programming (MINLP) solver to solve the resulting nonconvex formulations and reported results for the 2D ACRP with up to 7 aircraft on structured instances and up to 20 aircraft on un-structured instances. Alternative representations of pairwise aircraft separation based on conflict points have been proposed by several authors. In Vela et al. [2009c,b,a], the authors proposed several MILP formulations which aim to minimise fuel consumption, incorporate air traffic controller workload in the objective function, and account for the impact of uncertainty on trajectory prediction due to wind effects. Omer [2015] proposed a
space-discretised MILP formulation involving a finite set of turning angles. In contrast to most other approaches, the heading control manoeuvres consist of two actions: a first heading change for collision avoidance and a second turn to recover aircraft's initial heading. Rey et al. [2012, 2016c] proposed linear upper bounds for the ACRP with speed control only and the resulting MILP formulations can solve realistic large-scale instances to optimality within a few seconds.

More recently, nonlinear global optimization approaches received an increasing attention in the literature. Omer and Farges [2013] proposed a hybrid algorithm which uses the optimal solution of a MILP as the starting point for solving a nonlinear formulation of the same problem. Cafieri [2012] introduced a MINLP approach that provides a compact nonlinear separation condition for conflict resolution. This formulation was then used by Cafieri and Rey [2017] for conflict resolution with speed control only which highlights that subliminal speed control alone may not be sufficient to resolve all conflicts in dense traffic scenarios. Using a similar framework, Cafieri and Omheni [2017] presented a two-step approach where a maximum number of conflicts are first solved using speed control only and outstanding conflicts are solved by heading control. Cerulli et al. [2021] proposed a formulation based on bi-level optimization with multiple follower problems, each of which representing a two-aircraft separation problem. The authors presented two formulations, one using speed control only and another using heading control only. A cut generation algorithm is proposed to solve the corresponding bi-level optimization problems. Recently, Pelegrín and d'Ambrosio [2022] conducted a review of the literature on mathematical programming approaches for aircraft conflict resolution and have shown that the disjunctive linear separation conditions introduced by Rey and Hijazi [2017] and the shadow separation conditions are equivalent. This result combined with the theorem provided in this study – linking the disjunctive separation conditions with the definitional nonlinear separation conditions (see Theorem 3.1) – further motivates the comparison of both separation conditions in terms of computational performance which we address in our numerical experiments.

This review of the literature highlights that despite recent improvements in the development of optimization approaches for aircraft conflict resolution, there remains significant open challenges in the design of scalable and exact global optimization approaches. The rest of the chapter is organized as follows. Formulations based on disjunctive linear separation conditions for the deterministic aircraft conflict resolution problem are first presented (Section 3.3) and exact solution methods are then developed (Section 3.4). An extension to the robust aircraft conflict resolution problem is presented (Section 3.5) before numerical experiments are results are summarized (Section 3.6).

## 3.3 Aircraft conflict resolution formulations

We present mixed-integer formulations for the ACRP with speed, heading and altitude control. We first focus on the 2D problem under velocity, i.e. speed and heading, control which aims to represent a single flight level during cruise stage air traffic conditions. We present disjunctive linear separation conditions in Section 3.3.1. We then characterise 2D separation conditions for a pair of aircraft based on velocity control bounds in Section 3.3.2, before introducing a compact nonconvex formulation in Section 3.3.3. We then extend this non-convex formulation to the case of multiple separated flight levels (FL), denoted 2D+FL in Section 3.3.4.
### 3.3.1 Disjunctive linear separation conditions

Our goal is to find least-deviating conflict-free trajectories for a set of aircraft in cruise stage, i.e. flying at a fixed altitude and at constant speed. We denote $t = 0$ the time instant representing aircraft current positions at the time of decision, also referred to as aircraft initial positions. We assume that all aircraft are separated at $t = 0$ and seek to derive separation conditions to ensure that aircraft trajectories are separated for any time $t \ge 0$. Let $\mathcal{A}$ be the set of aircraft. For each $i \in \mathcal{A}$, $[\hat{x}_i, \hat{y}_i]^\top$ is the aircraft initial position in the 2D plane, $\hat{v}_i$ is its nominal speed (in NM/h) and $\hat{\theta}_i$ is its heading angle. Assuming uniform motion laws, aircraft motion can be described as: $p_i(t) = [x_i(t), y_i(t)]^\top$, where $x_i(t) = \hat{x}_i + q_i \hat{v}_i \cos(\hat{\theta}_i + \theta_i)t$ and $y_i(t) = \hat{y}_i + q_i \hat{v}_i \sin(\hat{\theta}_i + \theta_i)t$. In this model, the decision variables are $q_i$, which is the speed control variable that determines the acceleration or deceleration with regards to the initial speed $\hat{v}_i$ ($q_i$ equals to 1 means no speed variation) and $\theta_i$, which is the heading control variable that determines the deviation with regards to the initial trajectory ($\theta_i$ equal to 0 means no deviation in heading angle).

Let $\mathcal{P} = \{i, j \in \mathcal{A} : i < j\}$ be the set of aircraft pairs, the relative motion of $(i, j) \in \mathcal{P}$ is denoted $p_{ij}(t) = p_i(t) - p_j(t)$. Let $d$ be the minimum separation distance (e.g. 5 NM). We next define the notion of 2D separation for a pair of aircraft.

**Definition 3.1** (2D separated trajectories). The trajectories of a pair of aircraft $(i, j) \in \mathcal{P}$ are said to be 2D separated, i.e. conflict-free, if and only if:

$$ ||p_{ij}(t)|| \ge d, \quad \forall t \ge 0. \qquad (3.1) $$

Let $v_{ij} = v_i - v_j$ be the 2D relative velocity vector of $(i, j) \in \mathcal{P}$, i.e. $v_{ij} = [v_{ij,x}, v_{ij,y}]^\top$:

$$ v_{ij,x} = q_i \hat{v}_i \cos(\hat{\theta}_i + \theta_i) - q_j \hat{v}_j \cos(\hat{\theta}_j + \theta_j), \qquad (3.2a) $$

$$ v_{ij,y} = q_i \hat{v}_i \sin(\hat{\theta}_i + \theta_i) - q_j \hat{v}_j \sin(\hat{\theta}_j + \theta_j). \qquad (3.2b) $$

Aircraft relative velocity equations are linear with regards to speed control variables $q_i$ and $q_j$, but nonlinear with regards to heading control variables $\theta_i$ and $\theta_j$. Expanding the expression in Eq. (3.1) and denoting $\hat{p}_{ij}$ the relative initial position of aircraft pair $(i, j) \in \mathcal{P}$, we obtain a second-order polynomial function:

$$ f_{ij}(t) \equiv ||v_{ij}||^2 t^2 + 2\hat{p}_{ij} \cdot v_{ij}t + ||\hat{p}_{ij}||^2 - d^2 \ge 0. \qquad (3.3) $$

Since the coefficient of the second-order term is positive, $f_{ij}(t)$ admits a minimum which corresponds to the time instant of minimum separation between aircraft $i$ and $j$. Let $t_{ij}^{\min}(v_{ij})$ be the function corresponding to this time instant parameterised by aircraft relative velocity vector, then:

$$ t_{ij}^{\min}(v_{ij}) = \frac{-\hat{p}_{ij} \cdot v_{ij}}{||v_{ij}||^2}. \qquad (3.4) $$

Evaluating Eq. (3.3) at $t_{ij}^{\min}(v_{ij})$ yields a time-independent separation condition [Rey and Hijazi, 2017, Cafieri and Rey, 2017, Cafieri and Omheni, 2017]:

$$ f_{ij}(t_{ij}^{\min}(v_{ij})) = \frac{-(\hat{p}_{ij} \cdot v_{ij})^2}{||v_{ij}||^2} + ||\hat{p}_{ij}||^2 - d^2 \ge 0. \qquad (3.5) $$
Let $g_{ij}(v_{ij}) \equiv ||v_{ij}||^2 f_{ij}(t_{ij}^{\min}(v_{ij}))$. We abuse notation and rewrite $g_{ij}(\cdot)$ and $t_{ij}^{\min}(\cdot)$ as functions of aircraft velocity components $v_{ij,x}$ and $v_{ij,y}$. After expanding and factorising, we obtain:

$$g_{ij}(v_{ij,x}, v_{ij,y}) = v_{ij,x}^2(\hat{y}_{ij}^2 - d^2) + v_{ij,y}^2(\hat{x}_{ij}^2 - d^2) - v_{ij,x}v_{ij,y}(2\hat{x}_{ij}\hat{y}_{ij}) \geq 0. \quad (3.6)$$

Assuming aircraft are initially separated, if $t_{ij}^{\min}(v_{ij,x}, v_{ij,y}) \leq 0$, then they are diverging and do not incur any risk of future conflict. If $t_{ij}^{\min}(v_{ij,x}, v_{ij,y}) > 0$ and $g(v_{ij,x}, v_{ij,y}) \geq 0$, aircraft are converging but separation is ensured. Otherwise, if $t_{ij}^{\min}(v_{ij,x}, v_{ij,y}) > 0$ and $g(v_{ij,x}, v_{ij,y}) < 0$, there is a loss of separation and aircraft trajectories should be adjusted to avoid it. Hence, pairwise aircraft separation conditions for $(i, j) \in \mathcal{P}$ can be written as:

$$g_{ij}(v_{ij,x}, v_{ij,y}) \geq 0 \lor t_{ij}^{\min}(v_{ij,x}, v_{ij,y}) \leq 0. \quad (3.7)$$

To linearise the separation condition (3.6) with regards to variables $v_{ij,x}$ and $v_{ij,y}$, we adopt the approach proposed by Rey and Hijazi [2017] and recall the main steps hereafter. Observe that the solutions of the equation $g_{ij}(v_{ij,x}, v_{ij,y}) = 0$ can be identified by alternatively fixing variable $v_{ij,x}$ and $v_{ij,y}$, and calculating the roots of the resulting single-variable quadratic equations. Isolating each variable, we obtain the discriminants:

$$\left\{ \begin{array}{l} \Delta_{v_{ij,x}} = 4d^2 v_{ij,y}^2 (\hat{x}_{ij}^2 + \hat{y}_{ij}^2 - d^2), \\ \Delta_{v_{ij,y}} = 4d^2 v_{ij,x}^2 (\hat{x}_{ij}^2 + \hat{y}_{ij}^2 - d^2). \end{array} \right. \quad (3.8)$$

Assuming aircraft are initially separated, then $\hat{x}_{ij}^2 + \hat{y}_{ij}^2 - d^2 \geq 0$ and the discriminants are positive, and the roots of equation $g_{ij}(v_{ij,x}, v_{ij,y}) = 0$ are the lines defined by the system of equations:

$$(\hat{y}_{ij}^2 - d^2)v_{ij,x} - (\hat{x}_{ij}\hat{y}_{ij} + d\sqrt{\hat{x}_{ij}^2 + \hat{y}_{ij}^2 - d^2})v_{ij,y} = 0, \quad (3.9a)$$

$$(\hat{y}_{ij}^2 - d^2)v_{ij,x} - (\hat{x}_{ij}\hat{y}_{ij} - d\sqrt{\hat{x}_{ij}^2 + \hat{y}_{ij}^2 - d^2})v_{ij,y} = 0, \quad (3.9b)$$

$$(\hat{x}_{ij}^2 - d^2)v_{ij,y} - (\hat{x}_{ij}\hat{y}_{ij} + d\sqrt{\hat{x}_{ij}^2 + \hat{y}_{ij}^2 - d^2})v_{ij,x} = 0, \quad (3.9c)$$

$$(\hat{x}_{ij}^2 - d^2)v_{ij,y} - (\hat{x}_{ij}\hat{y}_{ij} - d\sqrt{\hat{x}_{ij}^2 + \hat{y}_{ij}^2 - d^2})v_{ij,x} = 0. \quad (3.9d)$$

Let us emphasise that if all coefficients in (3.9a)-(3.9d) are non-zero, then (3.9a) is identical to (3.9c) and (3.9b) is identical to (3.9d). Observe that

$$\hat{x}_{ij}\hat{y}_{ij} \pm d\sqrt{\hat{x}_{ij}^2 + \hat{y}_{ij}^2 - d^2} = 0, \\ d^2(\hat{x}_{ij}^2 + \hat{y}_{ij}^2 - d^2) = \hat{x}_{ij}^2\hat{y}_{ij}^2,$$

⇒

$$\hat{x}_{ij}^2\hat{y}_{ij}^2 - d^2(\hat{x}_{ij}^2 + \hat{y}_{ij}^2 - d^2) = 0,$$

⇔

$$(\hat{x}_{ij}^2 - d^2)(\hat{y}_{ij}^2 - d^2) = 0.$$

⇔

Eqs. (3.9a), (3.9b), (3.9c) and (3.9d) define two lines, denoted $R_1$ and $R_2$, in the plane $\{(v_{ij,x}, v_{ij,y}) \in \mathbb{R}^2\}$ and the sign of $g_{ij}(v_{ij,x}, v_{ij,y})$ can be characterised based on the position of $(v_{ij,x}, v_{ij,y})$ relative to these lines. Recall that according to Eq. (3.4), the sign of the dot
<figure><img src="image_4.png" /><figcaption>Figure 3.1: Illustration of a two-aircraft conflict in the plane $\{(v_{ij,x}, v_{ij,y}) \in \mathbb{R}^2\}$. The red lines represent the lines P and N. The dashed blue lines correspond to the linear equations $R_1$ and $R_2$ that are the roots of $g_{ij}(v_{ij,x}, v_{ij,y}) = 0$. The sign of $g_{ij}(v_{ij,x}, v_{ij,y})$ is shown by the + and - green symbols.</figcaption></figure>

product $\hat{p}_{ij} \cdot v_{ij}$ indicates aircraft convergence or divergence. Let (P) be the equation of the line corresponding to $\hat{p}_{ij} \cdot v_{ij}$:

$$v_{ij,x}\hat{x}_{ij} + v_{ij,y}\hat{y}_{ij} = 0. \qquad (P)$$

The line defined by (P) splits the plane $\{(v_{ij,x}, v_{ij,y}) \in \mathbb{R}^2\}$ in two half-planes, each of which representing converging and diverging trajectories, respectively. This is illustrated in Figure 3.1a which depicts a two-aircraft conflict in the plane $\{(v_{ij,x}, v_{ij,y}) \in \mathbb{R}^2\}$. The sign of $g_{ij}(v_{ij,x}, v_{ij,y})$ is shown by the + and - green symbols and the hashed pink region corresponds to $g_{ij}(v_{ij,x}, v_{ij,y}) \ge 0$. The hashed green half-plane delimited by (P) represents diverging trajectories, i.e. $t_{ij}^{\min}(v_{ij,x}, v_{ij,y}) \le 0$.

Consider the line normal to (P), denoted (N):

$$v_{ij,y}\hat{x}_{ij} - v_{ij,x}\hat{y}_{ij} = 0. \qquad (N)$$

Recall that any point $(v_{ij,x}, v_{ij,y})$ such that $t_{ij}^{\min} \le 0$ or $g_{ij}(v_{ij,x}, v_{ij,y}) \ge 0$ corresponds to a pair of conflict-free trajectories. Hence, the conflict-free region is nonconvex and is represented by the union of the green and pink hashed regions in Figure 3.1a. Reciprocally, the conflict region, corresponding to conflicting trajectories is convex and represented by the non-hashed region in Figure 3.1a. An equivalent expression of Eq. (3.9) was proposed by Frazzoli et al. [2001] which observed that the set of conflict-free trajectories could be characterised by the union of two half-planes. We next show through Lemmas 3.1 and 3.2 that (N) is a bisector of the angle formed by lines $R_1$ and $R_2$ in the conflict region ($g_{ij}(v_{ij,x}, v_{ij,y}) \le 0$) and can be used to generate two disjunctive but convex conflict-free regions.

**Lemma 3.1.** The lines (P) and (N) are bisectors of the angles formed by the two lines $R_1$ and $R_2$ representing the solutions of $g_{ij}(v_{ij,x}, v_{ij,y}) = 0$.

**Lemma 3.2.** $g_{ij}(v_{ij,x}, v_{ij,y}) \le 0$ for all points $(v_{ij,x}, v_{ij,y})$ of the normal line (N).
Lemmas 3.1 and 3.2 assert that (N) can be used to split the conflict-free region into two convex but disjunctive regions. We model this disjunction using variable $z_{ij} \in \{0, 1\}$ defined as:

$$v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \le 0, \quad \text{if } z_{ij} = 1, \quad \forall (i, j) \in \mathcal{P}, \qquad (3.10a)$$

$$v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \ge 0, \quad \text{if } z_{ij} = 0, \quad \forall (i, j) \in \mathcal{P}. \qquad (3.10b)$$

In each convex sub-region, the lines defined by (3.9a)-(3.9d) delineate the conflict-free region. The expressions of these lines depend on aircraft initial positions, i.e. $\hat{x}_{ij}, \hat{y}_{ij}$. Recall that we denote $R_1$ and $R_2$ the equation of these lines. Integer-linear separation conditions with regards to aircraft velocity components can be derived as follows:

$$v_{ij,y} \gamma_{ij}^{l} - v_{ij,x} \phi_{ij}^{l} \le 0, \quad \text{if } z_{ij} = 1, \quad \forall (i, j) \in \mathcal{P}, \qquad (3.11a)$$

$$v_{ij,y} \gamma_{ij}^{u} - v_{ij,x} \phi_{ij}^{u} \ge 0, \quad \text{if } z_{ij} = 0, \quad \forall (i, j) \in \mathcal{P}, \qquad (3.11b)$$

where $\gamma_{ij}^l, \phi_{ij}^l$ and $\gamma_{ij}^u, \phi_{ij}^u$ are coefficients of the lines (3.9a)-(3.9d) corresponding to the roots of $g_{ij}(v_{ij,x}, v_{ij,y}) = 0$. The proposed linear disjunction is illustrated in Figure 3.1b which depicts the resulting convex sub-regions corresponding to the disjunction $z_{ij} \in \{0, 1\}$ for a two-aircraft conflict. This leads to the following result.

**Theorem 3.1.** The disjunctive linear separation conditions (3.10)-(3.11) fully characterises the set of aircraft pairwise conflict-free trajectories as given by Eq. (3.1)

Theorem 3.1 asserts that the disjunctive separation conditions (3.10)-(3.11) are equivalent to the definitional nonlinear separation conditions. Further, these disjunctive separation conditions are linear with regards to aircraft velocity variables $v_{ij,x}$ and $v_{ij,y}$, and only require a single binary variable per pair of aircraft. This is expected to improve on the so-called shadow separation conditions which are also linear with regards to aircraft velocity variables, but require four binary variables per pair of aircraft [Pallottino et al., 2002, Alonso-Ayuso et al., 2011, 2016].

### 3.3.2 Characterisation of conflict-free trajectories based on velocity control bounds

To characterise the set of 2D conflict-free trajectories, we examine the relative velocity vector $v_{ij}$ as a function of aircraft trajectory control bounds. For each aircraft $i \in \mathcal{A}$, we assume that the speed rate variable is lower bounded by $\underline{q}_i$ and upper bounded by $\bar{q}_i$, i.e.

$$\underline{q}_i \le q_i \le \bar{q}_i, \quad \forall i \in \mathcal{A}. \qquad (3.12)$$

We assume that the heading deviation is lower and upper bounded by $\underline{\theta}_i$ and $\bar{\theta}_i$, i.e.

$$\underline{\theta}_i \le \theta_i \le \bar{\theta}_i, \quad \forall i \in \mathcal{A}. \qquad (3.13)$$

To derive lower and upper bounds on relative velocity components $v_{ij,x}$ and $v_{ij,y}$, we re-arrange Eq. (3.2) using trigonometric identities:

$$v_{ij,x} = q_i \hat{v}_i \cos(\hat{\theta}_i) \cos(\theta_i) - q_i \hat{v}_i \sin(\hat{\theta}_i) \sin(\theta_i) - q_j \hat{v}_j \cos(\hat{\theta}_j) \cos(\theta_j) + q_j \hat{v}_j \sin(\hat{\theta}_j) \sin(\theta_j), \qquad (3.14a)$$

$$v_{ij,y} = q_i \hat{v}_i \sin(\hat{\theta}_i) \cos(\theta_i) + q_i \hat{v}_i \cos(\hat{\theta}_i) \sin(\theta_i) - q_j \hat{v}_j \sin(\hat{\theta}_j) \cos(\theta_j) - q_j \hat{v}_j \cos(\hat{\theta}_j) \sin(\theta_j). \qquad (3.14b)$$
Let $v_{ij,x}, \bar{v}_{ij,x}$ and $v_{ij,y}, \bar{v}_{ij,y}$ be the lower and upper bounds for $v_{ij,x}$ and $v_{ij,y}$, respectively. These bounds can be determined using Eq. (3.14) and the bounds on speed and heading control provided in Eqs. (3.12) and (3.13). The derived bounds on the relative velocity components can be used to define a box in the plane $\{(v_{ij,x}, v_{ij,y}) \in \mathbb{R}^2\}$.

**Definition 3.2 (Relative velocity box).** Consider a pair of aircraft $(i, j) \in \mathcal{P}$. Let $\mathcal{B}_{ij}$ be the subset of $\mathbb{R}^2$ defined as

$$ \mathcal{B}_{ij} = \{(v_{ij,x}, v_{ij,y}) \in \mathbb{R}^2 : \underline{v}_{ij,x} \le v_{ij,x} \le \bar{v}_{ij,x}, \underline{v}_{ij,y} \le v_{ij,y} \le \bar{v}_{ij,y}\}. \quad (3.15) $$

$\mathcal{B}_{ij}$ is the relative velocity box of $(i, j) \in \mathcal{P}$.

The relative velocity box $\mathcal{B}_{ij}$ characterises all possible trajectories for the pair $(i, j) \in \mathcal{P}$ based on the available 2D deconfliction resources, i.e. speed and heading controls. To characterise the set of conflict-free trajectories of a pair of aircraft $(i, j) \in \mathcal{P}$, we compare the relative position of the relative velocity box $\mathcal{B}_{ij}$ with the conflict region of this pair of aircraft. Observe that the conflict region is convex and can be defined by reversing the inequalities (3.11) and omitting the disjunction $z_{ij} \in \{0, 1\}$.

**Definition 3.3 (Conflict region).** Consider a pair of aircraft $(i, j) \in \mathcal{P}$. Let $\mathcal{C}_{ij}$ be the subset of $\mathbb{R}^2$ defined as

$$ \mathcal{C}_{ij} = \{(v_{ij,x}, v_{ij,y}) \in \mathbb{R}^2 : v_{ij,y}\gamma_{ij}^l - v_{ij,x}\phi_{ij}^l \ge 0 \wedge v_{ij,y}\gamma_{ij}^u - v_{ij,x}\phi_{ij}^u \le 0\}. \quad (3.16) $$

$\mathcal{C}_{ij}$ is the conflict region of $(i, j) \in \mathcal{P}$.

The conflict region of a pair of aircraft represents the set of relative velocity vectors which corresponds to a conflict. The relative positions of the relative velocity box $\mathcal{B}_{ij}$ and the conflict region $\mathcal{C}_{ij}$ can be examined to determine the existence or not of a potential conflict. For any pair $(i, j) \in \mathcal{P}$, if $\mathcal{B}_{ij} \cap \mathcal{C}_{ij} = \emptyset$, then aircraft i and j are separated for any combination of controls; conversely if $\mathcal{B}_{ij} \subset \mathcal{C}_{ij}$ then i and j cannot be separated via speed or heading control within the assumed control bounds; otherwise, $\mathcal{B}_{ij}$ and $\mathcal{C}_{ij}$ intersect but do not completely overlap. This is illustrated in Figure 3.2 which illustrates the three possible cases. Figure 3.2a illustrates the case where aircraft i and j are separated for any combination of speed and heading control—we say that such pairs are conflict-free. Figure 3.2b depicts the case where $\mathcal{B}_{ij}$ and $\mathcal{C}_{ij}$ only partially intersect—we say that such pairs are separable. Last, Figure 3.2c illustrates the case where $\mathcal{B}_{ij} \subset \mathcal{C}_{ij}$—we say that such pairs are non-separable. The following propositions provide methods to efficiently determine if, given controls bounds on speed and heading, a pair of aircraft is either conflict-free or non-separable.

**Proposition 3.1 (Conflict-free aircraft pair).** Consider a pair of aircraft $(i, j) \in \mathcal{P}$, and let $\text{LP}(i, j)$ be the feasibility problem defined as:

$$ \text{LP}(i, j) : \left\{ \begin{array}{l} v_{ij,y}\gamma_{ij}^l - v_{ij,x}\phi_{ij}^l \ge 0, \\ v_{ij,y}\gamma_{ij}^u - v_{ij,x}\phi_{ij}^u \le 0, \\ (v_{ij,x}, v_{ij,y}) \in \mathcal{B}_{ij}. \end{array} \right. $$

The pair of aircraft $(i, j)$ is conflict-free for any 2D control if and only if the feasibility problem $\text{LP}(i, j)$ is infeasible.
**Proposition 3.2 (Non-separable aircraft pair).** Consider a pair of aircraft $(i, j) \in \mathcal{P}$. The pair $(i, j)$ is non-separable if and only if each of the four extreme points of $B_{ij}$ corresponds to a conflict.

Using Propositions 3.1 and 3.2, we can design an efficient pre-processing algorithm to partition the set of aircraft pairs $\mathcal{P}$ of a 2D ACRP instance into three categories: conflict-free pairs denoted $\mathcal{P}_F$, separable pairs denoted $\mathcal{P}_S$ and non-separable pairs $\mathcal{P}_I$.

**Algorithm 3: Pre-processing of aircraft pairs**

**Input:** $\mathcal{A}, \hat{\theta}, \hat{v}, \underline{q}, \bar{q}, \underline{\theta}, \bar{\theta}$
**Output:** $\mathcal{P}, \mathcal{P}_F, \mathcal{P}_S, \mathcal{P}_I$
1 $\mathcal{P} \leftarrow \{i, j \in \mathcal{A} : i < j\}$
2 $\mathcal{P}_F, \mathcal{P}_S, \mathcal{P}_I \leftarrow \emptyset$
3 **for** $(i, j) \in \mathcal{P}$ **do**
4 &nbsp;&nbsp;&nbsp;&nbsp;Solve $LP(i, j)$
5 &nbsp;&nbsp;&nbsp;&nbsp;**if** $LP(i, j)$ is infeasible **then**
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{P}_F \leftarrow \mathcal{P}_F \cup \{(i, j)\}$
7 &nbsp;&nbsp;&nbsp;&nbsp;**else**
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k \leftarrow 0$
9 &nbsp;&nbsp;&nbsp;&nbsp;**for** $(v_{ij,x}, v_{ij,y}) \in E(B_{ij})$ **do**
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $g_{ij}(v_{ij,x}, v_{ij,y}) < 0 \land t_{ij}^{\min}(v_{ij,x}, v_{ij,y}) > 0$ **then**
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k \leftarrow k + 1$
12 &nbsp;&nbsp;&nbsp;&nbsp;**if** $k = 4$ **then**
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{P}_I \leftarrow \mathcal{P}_I \cup \{(i, j)\}$
14 &nbsp;&nbsp;&nbsp;&nbsp;**else**
15 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{P}_S \leftarrow \mathcal{P}_S \cup \{(i, j)\}$

To identify conflict-free pairs, $LP(i, j)$ is solved and a pair of aircraft is conflict-free if and only if the LP is infeasible. Observe that the feasibility problem $LP(i, j)$ is linear and can be solved by enumerating all extreme points of its feasible region and test if this corresponds to a conflict or not. Since $LP(i, j)$ consists of four bound constraints and two shared constraints, there is a total of 13 extreme points to test (the combinations of the bound constraints of a variable can be excluded). To identify non-separable pairs, we denote $E(B_{ij})$ the set of extreme points of the relative velocity box $B_{ij}$ for any pair $(i, j) \in \mathcal{P}$ and use the separation condition (3.3) to determine if all extreme points are conflicts or not. This procedure is summarized in Algorithm 3 (we use boldface to denote vectors). Observe that pairwise variables and constraints need only to be indexed by the set of separable pairs $\mathcal{P}_S$ since pairs in $\mathcal{P}_F$ are always conflict-free. Further, any 2D conflict resolution problem such that $|\mathcal{P}_I| > 0$ is infeasible.

The nonconvex formulation of the 2D ARCP is summarized in Model 3.1.
<figure><img src="image_5.png" /><figcaption>Figure 3.2: Three aircraft configurations illustrating conflict-free (3.2a), separable (3.2b) and non-separable (3.2c) pairs. Aircraft initial positions and trajectories are depicted at the bottom part of each sub-figure, where black dots correspond to initial positions and red lines represent initial trajectories. The top part of each sub-figure shows a graph of the relative velocity box $B_{ij}$ represented by a black rectangle, and the conflict region $C_{ij}$ represented by a red hashed area. Dashed blue lines represent the solutions of $g(v_{ij,x}, v_{ij,y}) = 0$ ($R_1$ and $R_2$). In all cases, aircraft i and j have nominal speeds $v_i = v_j = 500$ NM/h and are initially separated by 30 NM in the x-axis direction.</figcaption></figure>

**Model 3.1 (Nonconvex 2D Formulation).**

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

s.t.

$$ v_{ij,x} = q_i \hat{v}_i \sin(\hat{\theta}_i + \theta_i) - q_j \hat{v}_j \sin(\hat{\theta}_j + \theta_j), \quad \forall (i,j) \in \mathcal{P}, \qquad (3.17b) $$

$$ v_{ij,y} = q_i \hat{v}_i \cos(\hat{\theta}_i + \theta_i) - q_j \hat{v}_j \cos(\hat{\theta}_j + \theta_j), \quad \forall (i,j) \in \mathcal{P}, \qquad (3.17c) $$

$$ v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \le 0, \quad \text{if } z_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}, \qquad (3.17d) $$

$$ v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \ge 0, \quad \text{if } z_{ij} = 0, \quad \forall (i,j) \in \mathcal{P}, \qquad (3.17e) $$

$$ v_{ij,y} \gamma_{ij}^{l} - v_{ij,x} \phi_{ij}^{l} \le 0, \quad \text{if } z_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}, \qquad (3.17f) $$

$$ v_{ij,x} \gamma_{ij}^{u} - v_{ij,y} \phi_{ij}^{u} \ge 0, \quad \text{if } z_{ij} = 0, \quad \forall (i,j) \in \mathcal{P}, \qquad (3.17g) $$

$$ v_{ij,x}, v_{ij,y} \in B_{ij}, \quad \forall (i,j) \in \mathcal{P}, \qquad (3.17h) $$

$$ z_{ij} \in \{0, 1\}, \quad \forall (i,j) \in \mathcal{P}, \qquad (3.17i) $$

$$ q_i \le q_i \le \bar{q}_i, \quad \forall i \in \mathcal{A}, \qquad (3.17j) $$

$$ \underline{\theta}_i \le \theta_i \le \bar{\theta}_i, \quad \forall i \in \mathcal{A}. \qquad (3.17k) $$

Model 3.1 provides a compact formulation for the ACRP with speed and heading control
which requires a single binary variable per pair of aircraft. Note that coefficients $\gamma_{ij}^l$, $\phi_{ij}^l$, $\gamma_{ij}^u$ and $\phi_{ij}^u$ (present in (3.10)) can be pre-processed based on the sign of $\hat{x}_{ij}$ and $\hat{y}_{ij}$. For implementation details, a fully reproducible formulation can be found at: https://github.com/acrp-lib/acrp-lib.

### 3.3.3 Complex formulation for the 2D aircraft conflict resolution problem

An alternative way to represent aircraft motion equation is via the complex number formulation introduced by Rey and Hijazi [2017]. Let $V_i$ be the complex number defined as:

$$V_i = q_i(\cos(\theta_i) + i \sin(\theta_i)), \quad \forall i \in \mathcal{A}. \qquad (3.18)$$

Representing $V_i$ in its polar form with $\delta_{i,x} = q_i \cos(\theta_i)$ and $\delta_{i,y} = q_i \sin(\theta_i)$, yield:

$$V_i = \delta_{i,x} + i\delta_{i,y}, \quad \forall i \in \mathcal{A}. \qquad (3.19)$$

The magnitude of $V_i$ is $|V_i| = \sqrt{\delta_{i,x}^2 + \delta_{i,y}^2} = q_i$ and its argument is $\arg(V_i) = \arctan \frac{\delta_{i,x}}{\delta_{i,y}} = \theta_i$. This approach is inspired by complex number formulation for the optimal power flow problem in power systems, such as in Hijazi et al. [2017] and Coffrin et al. [2015]. Accordingly, the relative motion equations of a pair of aircraft can be rewritten as:

$$v_{ij,x} = \delta_{i,x} \hat{v}_i \cos(\hat{\theta}_i) - \delta_{i,y} \hat{v}_i \sin(\hat{\theta}_i) - \delta_{j,x} \hat{v}_j \cos(\hat{\theta}_i) + \delta_{j,y} \hat{v}_j \sin(\hat{\theta}_j), \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.20a)$$

$$v_{ij,y} = \delta_{i,y} \hat{v}_i \cos(\hat{\theta}_i) - \delta_{i,x} \hat{v}_i \sin(\hat{\theta}_i) - \delta_{j,y} \hat{v}_j \cos(\hat{\theta}_i) + \delta_{j,x} \hat{v}_j \sin(\hat{\theta}_j), \quad \forall (i,j) \in \mathcal{P}_S. \quad (3.20b)$$

The variables $\delta_{i,x}$ and $\delta_{i,y}$ are used as the main control variables in this formulation and their bounds are derived from the bounds of the original control variables $q_i$ and $\theta_i$:

$$\underline{q}_i \cos(\max\{|\bar{\theta}_i|, |\underline{\theta}_i|\}) \le \delta_{i,x} \le \bar{q}_i, \quad \forall i \in \mathcal{A}, \qquad (3.21a)$$

$$\bar{q}_i \sin(\underline{\theta}_i) \le \delta_{i,y} \le \bar{q}_i \sin(\bar{\theta}_i), \quad \forall i \in \mathcal{A}. \qquad (3.21b)$$

The speed control constraint (3.12) can be reformulated in quadratic form as:

$$\underline{q}_i^2 \le \delta_{i,x}^2 + \delta_{i,y}^2, \quad \forall i \in \mathcal{A}, \qquad (3.22a)$$

$$\bar{q}_i^2 \ge \delta_{i,x}^2 + \delta_{i,y}^2, \quad \forall i \in \mathcal{A}. \qquad (3.22b)$$

The heading control constraint (3.13) can be reformulated in linear form as:

$$\delta_{i,x} \tan(\underline{\theta}_i) \le \delta_{i,y} \le \delta_{i,x} \tan(\bar{\theta}_i), \quad \forall i \in \mathcal{A}. \qquad (3.23a)$$

To design the objective function, we introduce a preference weight $w \in ]0, 1[$ to balance the trade-offs among velocity controls, i.e. speed and heading. We extend the objective function proposed by Rey and Hijazi [2017] as follows:

$$\underset{i \in \mathcal{A}}{\text{minimise}} \sum w \delta_{i,y}^2 + (1-w)(1-\delta_{i,x})^2. \qquad (3.24)$$

We next show that the objective function (3.24) attains its minimum value when aircraft have deviation-free trajectories, i.e. $q_i = 1$ and $\theta_i = 0$ for all $i \in \mathcal{A}$.
**Proposition 3.3.** The 2D objective function (3.24) is convex with regards to decision variables $\delta_{i,x}$ and $\delta_{i,y}$ for any value $w \in ]0, 1[$, and is minimal for $q_i = 1$ and $\theta_i = 0$ for all aircraft $i \in \mathcal{A}$.

Proposition 3.3 shows that the proposed 2D objective function achieves an optimal value for deviation-free aircraft trajectories for any preference weight $w \in ]0, 1[$. We will show in Section 3.6.2 that increasing $w$ increases the penalisation of heading deviations whereas decreasing $w$ increases the penalisation of speed deviations. The complex number nonconvex formulation for the 2D ACRP is summarized in Model 3.2.

**Model 3.2 (Complex number nonconvex 2D formulation).**

$$ \min \sum_{i \in \mathcal{A}} w \delta_{i,y}^2 + (1-w)(1-\delta_{i,x})^2, \qquad (3.25a) $$

s.t.

$$ v_{ij,x} = \delta_{i,x} \hat{v}_i \cos(\hat{\theta}_i) - \delta_{i,y} \hat{v}_i \sin(\hat{\theta}_i) - \delta_{j,x} \hat{v}_j \cos(\hat{\theta}_i) + \delta_{j,y} \hat{v}_j \sin(\hat{\theta}_j), \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.25b) $$

$$ v_{ij,y} = \delta_{i,y} \hat{v}_i \cos(\hat{\theta}_i) - \delta_{i,x} \hat{v}_i \sin(\hat{\theta}_i) - \delta_{j,y} \hat{v}_j \cos(\hat{\theta}_i) + \delta_{j,x} \hat{v}_j \sin(\hat{\theta}_j), \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.25c) $$

$$ v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \le 0, \quad \text{if } z_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.25d) $$

$$ v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \ge 0, \quad \text{if } z_{ij} = 0, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.25e) $$

$$ v_{ij,y} \gamma_{ij}^l - v_{ij,x} \phi_{ij}^l \le 0, \quad \text{if } z_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.25f) $$

$$ v_{ij,y} \gamma_{ij}^u - v_{ij,x} \phi_{ij}^u \ge 0, \quad \text{if } z_{ij} = 0, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.25g) $$

$$ q_i^2 \le \delta_{i,x}^2 + \delta_{i,y}^2 \le \bar{q}_i^2, \quad \forall i \in \mathcal{A}, \quad (3.25h) $$

$$ \delta_{i,x} \tan(\theta_i) \le \delta_{i,y} \le \delta_{i,x} \tan(\bar{\theta}_i), \quad \forall i \in \mathcal{A}, \quad (3.25i) $$

$$ q_i \cos(\max\{|\bar{\theta}_i|, |\theta_i|\}) \le \delta_{i,x} \le \bar{q}_i, \quad \forall i \in \mathcal{A}, \quad (3.25j) $$

$$ \bar{q}_i \sin(\theta_i) \le \delta_{i,y} \le \bar{q}_i \sin(\bar{\theta}_i), \quad \forall i \in \mathcal{A}, \quad (3.25k) $$

$$ \underline{v}_{ij,x} \le v_{ij,x} \le \bar{v}_{ij,x}, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.25l) $$

$$ \underline{v}_{ij,y} \le v_{ij,y} \le \bar{v}_{ij,y}, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.25m) $$

$$ z_{ij} \in \{0, 1\}, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.25n) $$

$$ \delta_{i,x}, \delta_{i,y} \in \mathbb{R}, \quad \forall i \in \mathcal{A}. \quad (3.25o) $$

Model 3.2 provides a compact formulation for the 2D ACRP with continuous speed and heading control variables which requires a single binary variable per pair of aircraft. This formulation is nonconvex due to the speed lower bound constraint (3.22a) which is nonconvex quadratic. Note that we use the convex hull reformulation of the linear On/Off constraints (3.10) and (3.11) as derived in Hijazi et al. [2014]. This reformulation does not require the introduction of auxiliary variables and is proved to provide the tightest continuous relaxation for each On/Off constraint. Coefficients $\gamma_{ij}^l, \phi_{ij}^l$ and $\gamma_{ij}^u, \phi_{ij}^u$ (present in (3.10) and (3.11)) can be pre-processed based on the sign of $\hat{x}_{ij}$ and $\hat{y}_{ij}$.

### 3.3.4 Complex number formulation for the 2D+FL aircraft conflict resolution problem

To model FL changes, we assume that each aircraft $i \in \mathcal{A}$ is initially assigned to a base FL denoted $\hat{\rho}_i$. We assume that adjacent FLs are vertically separated (e.g. by 1000 ft). Thus, we
only need to impose separation constraints on pairs of aircraft which share the same FL. Let $\mathcal{Z}_i$ denote the set of available FLs for each aircraft $i \in \mathcal{A}$, and consider the binary variable $\rho_{ik}$ defined as:

$$ \rho_{ik} = \begin{cases} 1 & \text{if aircraft } i \in \mathcal{A} \text{ is assigned to FL } k \in \mathcal{Z}_i, \\ 0 & \text{otherwise.} \end{cases} \qquad (3.26) $$

By design, $\hat{\rho}_i \in \mathcal{Z}_i$ and we require that each aircraft $i \in \mathcal{A}$ be assigned to exactly one FL in its reachable set $\mathcal{Z}_i$ via constraint (3.27):

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

Let $\mathcal{P}_S^3$ be the set of aircraft pairs which are not conflict-free and may share the same FL, i.e. $\mathcal{P}_S^3 = \{(i, j) \in \mathcal{P}_S \cup \mathcal{P}_I : \mathcal{Z}_i \cap \mathcal{Z}_j \neq \emptyset\}$. Let $\varphi_{ij}$ be the binary variable defined as:

$$ \varphi_{ij} = \begin{cases} 1 & \text{if aircraft } i \text{ and } j \text{ are assigned to the same FL,} \\ 0 & \text{otherwise.} \end{cases} \qquad (3.28) $$

Variable $\varphi_{ij}$ can be linked to binary variables $\rho_{ik}$ and $\rho_{jk}$ via the constraint:

$$ \rho_{ik} + \rho_{jk} \leq \varphi_{ij} + 1, \quad \forall (i, j) \in \mathcal{P}_S^3, k \in \mathcal{Z}_i \cap \mathcal{Z}_j. \qquad (3.29) $$

The separation conditions (3.10)-(3.11) are rewritten to account for altitude separation as follows:

$$ v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \leq 0, \quad \text{if } z_{ij} = 1 \text{ and } \varphi_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}_S^3, \qquad (3.30a) $$

$$ v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \geq 0, \quad \text{if } z_{ij} = 0 \text{ and } \varphi_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}_S^3, \qquad (3.30b) $$

$$ v_{ij,y} \gamma_{ij}^l - v_{ij,x} \phi_{ij}^l \leq 0, \quad \text{if } z_{ij} = 1 \text{ and } \varphi_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}_S^3, \qquad (3.30c) $$

$$ v_{ij,y} \gamma_{ij}^u - v_{ij,x} \phi_{ij}^u \geq 0, \quad \text{if } z_{ij} = 0 \text{ and } \varphi_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}_S^3. \qquad (3.30d) $$

FL changes are typically less desirable compared to other deconfliction manoeuvres such as speed or heading control [Bilimoria et al., 1996, Hu et al., 2002, Alonso-Ayuso et al., 2011]. This is due to several practical considerations including an increase in fuel consumption, passenger discomfort due to climbing or descending and the need for extended monitoring. Under these considerations, we propose a lexicographic optimization approach for the 2D+FL ACRP. The number of FL changes is first minimised; before the total 2D deviation of flights is minimised. The proposed objective function for minimising the number of FL changes is:

$$ \text{minimise} \sum_{i \in \mathcal{A}} \left| \sum_{k \in \mathcal{Z}_i} k \rho_{ik} - \hat{\rho}_i \right|, \qquad (3.31) $$

The resulting nonconvex lexicographic 2D+FL formulation is summarized in Model 3.3.

**Model 3.3 (Complex number nonconvex lexicographic 2D+FL formulation).**

$$ 1. \min \sum_{i \in \mathcal{A}} \left| \sum_{k \in \mathcal{Z}_i} k \rho_{ik} - \hat{\rho}_i \right|, \qquad (3.32a) $$
$$2. \min \sum_{i \in \mathcal{A}} w \delta_{i,y}^2 + (1-w)(1-\delta_{i,x})^2, \qquad (3.32b)$$

s.t.

$$v_{ij,x} = \delta_{i,x} \hat{v}_i \cos(\hat{\theta}_i) - \delta_{i,y} \hat{v}_i \sin(\hat{\theta}_i) - \delta_{j,x} \hat{v}_j \cos(\hat{\theta}_j) + \delta_{j,y} \hat{v}_j \sin(\hat{\theta}_j), \quad \forall (i,j) \in \mathcal{P}_S^3, \quad (3.32c)$$

$$v_{ij,y} = \delta_{i,y} \hat{v}_i \cos(\hat{\theta}_i) - \delta_{i,x} \hat{v}_i \sin(\hat{\theta}_i) - \delta_{j,y} \hat{v}_j \cos(\hat{\theta}_j) + \delta_{j,x} \hat{v}_j \sin(\hat{\theta}_j), \quad \forall (i,j) \in \mathcal{P}_S^3, \quad (3.32d)$$

$$v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \le 0, \quad \text{if } z_{ij} = 1 \quad \text{and} \quad \varphi_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}_S^3, \quad (3.32e)$$

$$v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \ge 0, \quad \text{if } z_{ij} = 0 \quad \text{and} \quad \varphi_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}_S^3, \quad (3.32f)$$

$$v_{ij,y} \gamma_{ij}^l - v_{ij,x} \phi_{ij}^l \le 0, \quad \text{if } z_{ij} = 1 \quad \text{and} \quad \varphi_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}_S^3, \quad (3.32g)$$

$$v_{ij,y} \gamma_{ij}^u - v_{ij,x} \phi_{ij}^u \ge 0, \quad \text{if } z_{ij} = 0 \quad \text{and} \quad \varphi_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}_S^3, \quad (3.32h)$$

$$\underline{q}_i^2 \le \delta_{i,x}^2 + \delta_{i,y}^2 \le \bar{q}_i^2, \quad \forall i \in \mathcal{A}, \quad (3.32i)$$

$$\delta_{i,x} \tan(\underline{\theta}_i) \le \delta_{i,y} \le \delta_{i,x} \tan(\bar{\theta}_i), \quad \forall i \in \mathcal{A}, \quad (3.32j)$$

$$\underline{q}_i \cos(\max\{|\bar{\theta}_i|, |\underline{\theta}_i|\}) \le \delta_{i,x} \le \bar{q}_i, \quad \forall i \in \mathcal{A}, \quad (3.32k)$$

$$\bar{q}_i \sin(\underline{\theta}_i) \le \delta_{i,y} \le \bar{q}_i \sin(\bar{\theta}_i), \quad \forall i \in \mathcal{A}, \quad (3.32l)$$

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

$$\rho_{ik} + \rho_{jk} \le \varphi_{ij} + 1, \quad \forall (i,j) \in \mathcal{P}_S^3, k \in \mathcal{Z}_i \cap \mathcal{Z}_j, \quad (3.32n)$$

$$\underline{v}_{ij,x} \le v_{ij,x} \le \bar{v}_{ij,x}, \quad \forall (i,j) \in \mathcal{P}_S^3, \quad (3.32o)$$

$$\underline{v}_{ij,y} \le v_{ij,y} \le \bar{v}_{ij,y}, \quad \forall (i,j) \in \mathcal{P}_S^3, \quad (3.32p)$$

$$z_{ij} \in \{0, 1\}, \quad \forall (i,j) \in \mathcal{P}_S^3, \quad (3.32q)$$

$$\delta_{i,x}, \delta_{i,y} \in \mathbb{R}, \quad \forall i \in \mathcal{A}, \quad (3.32r)$$

$$\rho_{ik} \in \{0, 1\}, \quad \forall i \in \mathcal{A}, k \in \mathcal{Z}_i. \quad (3.32s)$$

Compared to Model 3.2, Model 3.3 requires additional binary decision variables $\rho_{ik}$ and $\varphi_{ij}$. The former is used to assign aircraft to separated FLs and the latter ensures that aircraft sharing the same FL are separated via the 2D separation conditions.

## 3.4 Exact solution methods for the aircraft conflict resolution problem

We present exact solution methods for the 2D and the 2D+FL ACRPs that build on and extend the convex relaxations presented by Rey and Hijazi [2017]. We first present a convex relaxation of the 2D ACRP that fully relaxes the speed control constraint in the complex number formulation (Section 3.4.1). This relaxation yields a Mixed-Integer Quadratic Program (MIQP) which solution may violate the speed control bounds of the problem. To eliminate these potential violations, we incorporate convex quadratic constraints and piecewise linear outer approximations in a Mixed-Integer Quadratically Constrained Program (MIQCP). We then propose a constraint generation algorithm to iteratively refine the piece-wise linear approximation and show that our approach converges to optimal solutions of the 2D ACRP (Section 3.4.2).

To solve the lexicographic optimization formulation for the 2D+FL ACRP, we propose a
two-step decomposition approach. We first solve a restricted flight assignment problem that only implicitly accounts for aircraft trajectories and yields an optimal solution with regards to the first objective function (total FL deviation). We then use the optimal solution of the flight assignment formulation to assign aircraft to FLs and solve a series of 2D problems (one per FL) to construct an optimal solution with regards to the second objective function (total 2D deviation). Both steps are iterated until a global solution is found (Section 3.4.3). We next present in detail the proposed exact solution methods.

### 3.4.1 Mixed-integer quadratic relaxation for the 2D ACRP

An initial convex relaxation of Model 3.2 was proposed by Rey and Hijazi [2017] by relaxing the speed control constraint (3.22). The resulting formulation is a MIQP summarized in Model 3.4.

Model 3.4 (MIQP 2D formulation).

$$ \min \sum_{i \in \mathcal{A}} w \delta_{i,y}^2 + (1-w)(1-\delta_{i,x})^2, \qquad (3.33a) $$

s.t.

$$ v_{ij,x} = \delta_{i,x} \hat{v}_i \cos(\hat{\theta}_i) - \delta_{i,y} \hat{v}_i \sin(\hat{\theta}_i) - \delta_{j,x} \hat{v}_j \cos(\hat{\theta}_i) + \delta_{j,y} \hat{v}_j \sin(\hat{\theta}_j), \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.33b) $$

$$ v_{ij,y} = \delta_{i,y} \hat{v}_i \cos(\hat{\theta}_i) - \delta_{i,x} \hat{v}_i \sin(\hat{\theta}_i) - \delta_{j,y} \hat{v}_j \cos(\hat{\theta}_i) + \delta_{j,x} \hat{v}_j \sin(\hat{\theta}_j), \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.33c) $$

$$ v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \le 0, \quad \text{if } z_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.33d) $$

$$ v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \ge 0, \quad \text{if } z_{ij} = 0, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.33e) $$

$$ v_{ij,y} \gamma_{ij}^l - v_{ij,x} \phi_{ij}^l \le 0, \quad \text{if } z_{ij} = 1, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.33f) $$

$$ v_{ij,y} \gamma_{ij}^u - v_{ij,x} \phi_{ij}^u \ge 0, \quad \text{if } z_{ij} = 0, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.33g) $$

$$ \delta_{i,x} \tan(\underline{\theta}_i) \le \delta_{i,y} \le \delta_{i,x} \tan(\bar{\theta}_i), \quad \forall i \in \mathcal{A}, \quad (3.33h) $$

$$ \underline{q}_i \cos(\max\{|\bar{\theta}_i|, |\underline{\theta}_i|\}) \le \delta_{i,x} \le \bar{q}_i, \quad \forall i \in \mathcal{A}, \quad (3.33i) $$

$$ \bar{q}_i \sin(\underline{\theta}_i) \le \delta_{i,y} \le \bar{q}_i \sin(\bar{\theta}_i), \quad \forall i \in \mathcal{A}, \quad (3.33j) $$

$$ \underline{v}_{ij,x} \le v_{ij,x} \le \bar{v}_{ij,x}, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.33k) $$

$$ \underline{v}_{ij,y} \le v_{ij,y} \le \bar{v}_{ij,y}, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.33l) $$

$$ z_{ij} \in \{0, 1\}, \quad \forall (i,j) \in \mathcal{P}_S, \quad (3.33m) $$

$$ \delta_{i,x}, \delta_{i,y} \in \mathbb{R}, \quad \forall i \in \mathcal{A}. \quad (3.33n) $$

Model 3.4 yields a lower bound on the optimal objective value of Model 3.2 and a solution which is a global optimum if the relaxed constraint (3.22) is not violated.

### 3.4.2 Mixed-integer quadratically constrained relaxation and constraint generation algorithm

To tighten the MIQP relaxation given in Section 3.4.1, we build on and extend the MIQCP relaxation proposed by Rey and Hijazi [2017] by incorporating the speed control constraint (3.22) using convex quadratic and piecewise linear constraints. Observe that the speed upper bound constraint (3.22b) is convex quadratic, hence it can be incorporated directly in the MIQCP formulation.
To incorporate the speed lower bound constraint (3.22a), for each aircraft $i \in \mathcal{A}$, we introduce real variables $\tilde{\delta}_{i,x}$ and $\tilde{\delta}_{i,y}$ to approximate $\delta_{i,x}^2$ and $\delta_{i,y}^2$, respectively, via convex quadratic constraints:

$$ \tilde{\delta}_{i,x} \geq \delta_{i,x}^2, \quad \forall i \in \mathcal{A}, \qquad (3.34a) $$

$$ \tilde{\delta}_{i,y} \geq \delta_{i,y}^2, \quad \forall i \in \mathcal{A}. \qquad (3.34b) $$

Accordingly, we require:

$$ q_i^2 \leq \tilde{\delta}_{i,x} + \tilde{\delta}_{i,y}, \quad \forall i \in \mathcal{A}. \qquad (3.35) $$

To impose upper bounds on $\tilde{\delta}_{i,x}$ and $\tilde{\delta}_{i,y}$, we use their McCormick envelopes [McCormick, 1976]:

$$ \tilde{\delta}_{i,x} \leq (\bar{q}_i + \underline{q}_i \cos(\max\{|\theta_i|, |\bar{\theta}_i|\}))\delta_{i,x} - \bar{q}_i \underline{q}_i \cos(\max\{|\theta_i|, |\bar{\theta}_i|\}), \quad \forall i \in \mathcal{A}. \qquad (3.36a) $$

$$ \tilde{\delta}_{i,y} \leq \bar{q}_i(\sin(\theta_i) + \sin(\bar{\theta}_i))\delta_{i,y} - \bar{q}_i^2 \sin(\bar{\theta}_i) \sin(\theta_i), \quad \forall i \in \mathcal{A}. \qquad (3.36b) $$

Constraints (3.34)-(3.36) restrict variables $\tilde{\delta}_{i,x}$ and $\tilde{\delta}_{i,y}$ to convex regions thus providing an initial relaxation of the speed lower bound constraint (3.22a). This initial relaxation is illustrated in Figure 3.3 which depicts the variation of $\delta_{i,x}^2$ and $\delta_{i,y}^2$ (in red) over the domain of $\delta_{i,x}$ and $\delta_{i,y}$ for realistic speed and heading control bounds. The green lines in Figures 3.3a and 3.3b represent the initial McCormick envelopes. This relaxation may still yield infeasible aircraft speeds. To refine this convex relaxation, we introduce mixed-integer cuts that can be generated on-the-fly in a constraint generation algorithm.

We next present the general structure of the proposed mixed-integer cuts before discussing how these cuts are generated iteratively. Let $\Gamma_{i,x}$ and $\Gamma_{i,y}$ be partitions of the domain of variables $\delta_{i,x}$ and $\delta_{i,y}$, respectively. Since $\delta_{i,x}^2$ and $\delta_{i,y}^2$ are convex, any line joining two extremities of a segment in $\Gamma_{i,x}$ and $\Gamma_{i,y}$ refines the initial McCormick upper envelopes. Let $\alpha_{i,x}^k$ (resp. $\alpha_{i,y}^k$) and $\beta_{i,x}^k$ (resp. $\beta_{i,y}^k$) be the slope and the intercept corresponding to segment $k \in \Gamma_{i,x}$ (resp. $k \in \Gamma_{i,y}$). Further, let $s_{i,x}^k$ (resp. $s_{i,y}^k$) be a binary variable taking value 1 if $\delta_{i,x}$ (resp. $\delta_{i,y}$) belongs to segment $k$ of partition $\Gamma_{i,x}$ (resp. $\Gamma_{i,y}$). The proposed mixed-integer cuts take the form of:

$$ \tilde{\delta}_{i,x} \leq \alpha_{i,x}^k \delta_{i,x} + \beta_{i,x}^k, \quad \text{if } s_{i,x}^k = 1, \quad \forall i \in \mathcal{A}, k \in \Gamma_{i,x}, \qquad (3.37a) $$

$$ \tilde{\delta}_{i,y} \leq \alpha_{i,y}^k \delta_{i,y} + \beta_{i,y}^k, \quad \text{if } s_{i,y}^k = 1, \quad \forall i \in \mathcal{A}, k \in \Gamma_{i,y}. \qquad (3.37b) $$

Let $\delta_{i,x}^k$ and $\bar{\delta}_{i,x}^k$ (resp. $\underline{\delta}_{i,y}^k$ and $\bar{\delta}_{i,y}^k$) be the extremities of segment $k \in \Gamma_{i,x}$ (resp. $k \in \Gamma_{i,y}$). Binary variables $s_{i,x}^k$ and $s_{i,y}^k$ are defined as:

$$ s_{i,x}^k \equiv \begin{cases} 1 & \text{if } \delta_{i,x} \in [\underline{\delta}_{i,x}^k, \bar{\delta}_{i,x}^k], \\ 0 & \text{otherwise,} \end{cases} \quad \forall i \in \mathcal{A}, k \in \Gamma_{i,x}, \qquad (3.38a) $$

$$ s_{i,y}^k \equiv \begin{cases} 1 & \text{if } \delta_{i,y} \in [\underline{\delta}_{i,y}^k, \bar{\delta}_{i,y}^k], \\ 0 & \text{otherwise.} \end{cases} \quad \forall i \in \mathcal{A}, k \in \Gamma_{i,y}. \qquad (3.38b) $$
We require the following cut selection constraints:

$$ \sum_{k \in \Gamma_{i,x}} s_{i,x}^k = 1 \quad \forall i \in \mathcal{A}, \qquad (3.39a) $$

$$ \sum_{k \in \Gamma_{i,y}} s_{i,y}^k = 1 \quad \forall i \in \mathcal{A}. \qquad (3.39b) $$

The proposed mixed-integer cuts are illustrated in Figures 3.3a and 3.3b for the case of $|\Gamma_{i,x}| = 2$ and $|\Gamma_{i,y}| = 4$ segments, respectively. The shaded purple region represents the feasible region of $\tilde{\delta}_{i,x}$ and $\tilde{\delta}_{i,y}$ after imposing the mixed-integer cuts. The resulting formulation is a MIQCP with mixed-integer cuts summarized in Model 3.5. This formulation can be solved by off-the-shelf commercial optimization software and provides a relaxation of the Model 3.2 which can be tightened as desired by refining the partitions $\Gamma_{i,x}$ and $\Gamma_{i,y}$ of the domain of variables $\delta_{i,x}$ and $\delta_{i,y}$, respectively.

### Model 3.5 (MIQCP 2D formulation).

$$ \min \sum_{i \in \mathcal{A}} w \delta_{i,y}^2 + (1-w)(1-\delta_{i,x})^2, \qquad (3.40a) $$

s.t.

$$ v_{ij,x} = \delta_{i,x} \hat{v}_i \cos(\hat{\theta}_i) - \delta_{i,y} \hat{v}_i \sin(\hat{\theta}_i) - \delta_{j,x} \hat{v}_j \cos(\hat{\theta}_i) + \delta_{j,y} \hat{v}_j \sin(\hat{\theta}_j), \quad \forall (i,j) \in \mathcal{P}_S, \qquad (3.40b) $$

$$ v_{ij,y} = \delta_{i,y} \hat{v}_i \cos(\hat{\theta}_i) - \delta_{i,x} \hat{v}_i \sin(\hat{\theta}_i) - \delta_{j,y} \hat{v}_j \cos(\hat{\theta}_i) + \delta_{j,x} \hat{v}_j \sin(\hat{\theta}_j), \quad \forall (i,j) \in \mathcal{P}_S, \qquad (3.40c) $$

$$ v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \le 0, \quad \text{if } z_{ij} = 1, \qquad \forall (i,j) \in \mathcal{P}_S, \qquad (3.40d) $$

$$ v_{ij,y} \hat{x}_{ij} - v_{ij,x} \hat{y}_{ij} \ge 0, \quad \text{if } z_{ij} = 0, \qquad \forall (i,j) \in \mathcal{P}_S, \qquad (3.40e) $$

$$ v_{ij,y} \gamma_{ij}^l - v_{ij,x} \phi_{ij}^l \le 0, \quad \text{if } z_{ij} = 1, \qquad \forall (i,j) \in \mathcal{P}_S, \qquad (3.40f) $$

$$ v_{ij,y} \gamma_{ij}^u - v_{ij,x} \phi_{ij}^u \ge 0, \quad \text{if } z_{ij} = 0, \qquad \forall (i,j) \in \mathcal{P}_S, \qquad (3.40g) $$

$$ \delta_{i,x} \tan(\theta_i) \le \delta_{i,y} \le \delta_{i,x} \tan(\bar{\theta}_i), \qquad \forall i \in \mathcal{A}, \qquad (3.40h) $$

$$ \tilde{\delta}_{i,x} \ge \delta_{i,x}^2, \qquad \forall i \in \mathcal{A}, \qquad (3.40i) $$

$$ \tilde{\delta}_{i,y} \ge \delta_{i,y}^2, \qquad \forall i \in \mathcal{A}, \qquad (3.40j) $$

$$ q_i^2 \le \tilde{\delta}_{i,x} + \tilde{\delta}_{i,y}, \qquad \forall i \in \mathcal{A}, \qquad (3.40k) $$

$$ \tilde{\delta}_{i,x} \le (\bar{q}_i + q_i \cos(\max\{|\theta_i|, |\bar{\theta}_i|\}))\delta_{i,x} - \bar{q}_i q_i \cos(\max\{|\theta_i|, |\bar{\theta}_i|\}), \qquad \forall i \in \mathcal{A}, \qquad (3.40l) $$

$$ \tilde{\delta}_{i,y} \le \bar{q}_i(\sin(\theta_i) + \sin(\bar{\theta}_i))\delta_{i,y} - \bar{q}_i^2 \sin(\bar{\theta}_i) \sin(\theta_i), \qquad \forall i \in \mathcal{A}, \qquad (3.40m) $$

$$ \tilde{\delta}_{i,x} \le \alpha_{i,x}^k \delta_{i,x} + \beta_{i,x}^k, \quad \text{if } s_{i,x}^k = 1, \qquad \forall i \in \mathcal{A}, k \in \Gamma_{i,x}, \qquad (3.40n) $$

$$ \tilde{\delta}_{i,y} \le \alpha_{i,y}^k \delta_{i,y} + \beta_{i,y}^k, \quad \text{if } s_{i,y}^k = 1, \qquad \forall i \in \mathcal{A}, k \in \Gamma_{i,y}, \qquad (3.40o) $$

$$ s_{i,x}^k = 1 \quad \text{if } \delta_{i,x} \in [\bar{\delta}_{i,x}^k, \bar{\delta}_{i,x}^k], \qquad \forall i \in \mathcal{A}, k \in \Gamma_{i,x}, \qquad (3.40p) $$

$$ s_{i,y}^k = 1 \quad \text{if } \delta_{i,y} \in [\bar{\delta}_{i,y}^k, \bar{\delta}_{i,y}^k], \qquad \forall i \in \mathcal{A}, k \in \Gamma_{i,y}, \qquad (3.40q) $$

$$ \sum_{k \in \Gamma_{i,x}} s_{i,x}^k = 1 \qquad \forall i \in \mathcal{A}, \qquad (3.40r) $$

$$ \sum_{k \in \Gamma_{i,y}} s_{i,y}^k = 1 \qquad \forall i \in \mathcal{A}, \qquad (3.40s) $$

$$ q_i^2 \le \delta_{i,x}^2 + \delta_{i,y}^2, \qquad \forall i \in \mathcal{A}, \qquad (3.40t) $$

$$ q_i \cos(\max\{|\bar{\theta}_i|, |\theta_i|\}) \le \delta_{i,x} \le \bar{q}_i, \qquad \forall i \in \mathcal{A}, \qquad (3.40u) $$
$$ \bar{q}_i \sin(\theta_i) \le \delta_{i,y} \le \bar{q}_i \sin(\bar{\theta}_i), \quad \forall i \in A, \quad (3.40v) $$

$$ \underline{v}_{ij,x} \le v_{ij,x} \le \bar{v}_{ij,x}, \quad \forall (i, j) \in \mathcal{P}_S, \quad (3.40w) $$

$$ \underline{v}_{ij,y} \le v_{ij,y} \le \bar{v}_{ij,y}, \quad \forall (i, j) \in \mathcal{P}_S, \quad (3.40x) $$

$$ z_{ij} \in \{0, 1\}, \quad \forall (i, j) \in \mathcal{P}_S, \quad (3.40y) $$

$$ \delta_{i,x}, \delta_{i,y} \in \mathbb{R}, \quad \forall i \in A, \quad (3.40z) $$

$$ \tilde{\delta}_{i,x}, \tilde{\delta}_{i,y} \in \mathbb{R}, \quad \forall i \in A, \quad (3.40aa) $$

$$ s_{i,x}^k \in \{0, 1\}, \quad \forall i \in A, k \in \Gamma_{i,x}, \quad (3.40ab) $$

$$ s_{i,y}^k \in \{0, 1\}, \quad \forall i \in A, k \in \Gamma_{i,y}. \quad (3.40ac) $$

<figure><img src="image_6.png" /><figcaption>Figure 3.3: Piecewise linear approximation of $\delta_{i,x}^2$ and $\delta_{i,y}^2$ (in red). The green lines represent the initial McCormick relaxation given by (3.36). The purple shaded regions represent the refined feasible region of $\tilde{\delta}_{i,x}$ and $\tilde{\delta}_{i,y}$ after imposing the mixed-integer cuts (3.37). Control bounds are $q_i = 0.94$, $\bar{q}_i = 1.03$, $\underline{\theta}_i = -\frac{\pi}{6}$ and $\bar{\theta}_i = +\frac{\pi}{6}$.</figcaption></figure>

We propose a constraint generation algorithm which starts from a relaxed formulation and iteratively refines the piecewise linear outer approximation of the quadratic terms $\delta_{i,x}^2$ and $\delta_{i,y}^2$ via mixed-integer cuts (3.37). At each iteration, the proposed constraint generation algorithm examines the solution $(\delta_{i,x}, \delta_{i,y})$ of each aircraft $i \in A$ for speed lower bound violations. If $\delta_{i,x}^2 + \delta_{i,y}^2 < q_i^2$, then at least one of the relaxed auxiliary variables $\tilde{\delta}_{i,x}$ or $\tilde{\delta}_{i,y}$ must be such that $\tilde{\delta}_{i,x} > \delta_{i,x}^2$ or $\tilde{\delta}_{i,y} > \delta_{i,y}^2$. To eliminate the current infeasible solution, the partition(s) $\Gamma_{i,x}$ or $\Gamma_{i,y}$ of the violating variable(s) is augmented by dividing the segment(s) corresponding to $\delta_{i,x}$ or $\delta_{i,y}$ into two segments which meet at $\delta_{i,x}$ or $\delta_{i,y}$ and corresponding binary variable(s) $s_{i,x}^k$ or $s_{i,y}^k$ are added to the formulation. The process is repeated until all aircraft have feasible speed profiles which corresponds to a global optimum of Model 3.2. To further improve the convergence of the solution algorithm, after each solve of the relaxed formulation 3.5, Model 3.2 is solved as a nonlinear program (NLP) by fixing variable $z_{ij}$. The pseudo-code of the resulting algorithm is summarized in Algorithm 4.