of the origin stay is observed. By definition, the observed departure time precedes the actual one. Thus, $\epsilon_d^i$ represents the positive bias between these two values.

* Reciprocally, the observed arrival time of trip i is defined as the time when the first static event of the destination stay is observed. By definition, the observed arrival time follows the actual one. Thus, $\epsilon_a^i$ represents the positive bias between these two values.

* The observed travel time $T_{obs}^i$ of trip i is defined as the time elapsed between two consecutive stays, i.e., between the observed departure and the observed arrival times. It is an overestimate of the actual travel time $T^i$.

Based on these definitions, we thus have:

$$ \epsilon^i = T_{obs}^i - T^i = \epsilon_d^i + \epsilon_a^i \quad (3.3) $$

## The Framework: Main Steps

In the following, we provide a summary of the speed estimatino framework (see Fig. 3.1) by focusing on the procedures for travel time bias removal and for speed estimation. The interested reader can refer to [132] for the complete description of the approach.

The framework assumes the availability as inputs of: (i) a network partitioning into a set $\mathcal{R}$ of homogeneous regions; (ii) a set of relevant regional trip length information. The estimation of trip length information must be performed beforehand from external data sources, on the same network partitioning that defines the adopted regional scale and paths. We relied on previous work for the definition of the network partitioning [147] as well as for the estimation of the regional trip lengths [146]. Specifically, this input corresponds to the average regional trip length in each region along each possible regional path. We assume this information to be available as a distance matrix $\tilde{L}$, whose rows (i) are the different possible regional paths, and columns (j) are the different regions resulting from the spatial tessellation. The cell value at (i, j) corresponds to the average distance travelled in the $j^{th}$ region, when travelling along the $i^{th}$ regional path $P$. It is equal to zero if the path $P$ does not cross the $j^{th}$ region. This matrix was assumed to be constant over time in our framework, but a time-dependent generalisation could be considered as well if such information can be derived independently from another dataset [148].

## Travel Time Bias Modelling and Removal

We modelled the bias as an additive component on the exact travel time of any trip observed via UAPD<fnref n="7" />. Let $P$ be a regional path and $i$ be an individual trip observed along $P$. We thus have:

$$ T_{P,obs}^i = T_P^i + \epsilon^i \quad (3.4) $$

where $T_{P,obs}^i$, $T_P^i$ and $\epsilon^i$ are, respectively, the observed travel time of trip $i$ along $P$, its exact travel time, and the corresponding travel time bias. The latter can be expressed as:

$$ \epsilon^i = \epsilon_d^i + \epsilon_a^i \quad (3.5) $$

Estimating this bias on an individual trip basis is difficult. However, the estimation of its average is less challenging and can allow de-skewing, on average, the observed travel times. To compute the average bias, we proposed merging overlapping trips and averaging their observed travel times to build a unique aggregated biased travel time information per path and per time slot.

<footnote n="7">⁷Non-additive and more complex forms of bias could have been considered, but we opted for the most natural and simple representation in this work.</footnote>
<figure><img src="image_38.png" /><figcaption>**Figure 3.1:** Methodological framework for regional speed reconstruction.</figcaption></figure>

Let $t$ represent a generic time slot, and let $I_P^t$ be the set of overlapping trips along $P$ that reach destination during $t$, with $n_{t,P} = |I_P^t|$. Averaging Eq. 3.4 over $I_P^t$ gives:

$$ \bar{T}_P^t = \bar{T}_{P,obs}^t - \bar{\epsilon}_P^t \qquad (3.6) $$

where $\bar{T}_P^t$, $\bar{T}_{P,obs}^t$ and $\bar{\epsilon}_P^t$ are, respectively, the average actual travel time, the average observed travel time, and the average bias of trips from set $I_P^t$.

Eq. 3.6 can be simplified by assuming that the bias is independent of the trip path and time<fnref n="8" />. Thus, the distribution of individual biases $\epsilon^i$ can be modelled via a unique random variable $X$. The construction of such a model, and the estimation of its first moment

<footnote n="8">⁸The frequency at which mobile phone users are observed through the mobile network, i.e., the inter-event time, is known to depend on the hour of the day [149]. In particular, longer inter-event times are observed during nighttime and early morning. This could make the bias, in turn, dependent on time. However, these results normally account for all individuals, including the ones that are static and sleeping, while, in this study, we were exclusively interested in the moving ones, for which it might be reasonable to make the simplistic hypothesis of time independence. The spatial independence of the bias is similarly debatable and</footnote>
$\mu_X \equiv E(X)$, can offer an approximation of $\bar{\epsilon}_P^t$ allowing the de-skewing of $\bar{T}_{P,obs}^t$, provided that the sample of individuals associated with this time slot and path is large enough:

$$ \bar{T}_P^t \approx \bar{T}_{P,obs}^t - \mu_X \qquad (3.7) $$

Hence, if the departure and arrival temporal offsets $\epsilon_d^i$ and $\epsilon_a^i$ are themselves modelled by the same random variable Y, Eq. 3.5 gives $X = 2Y$. Now, as discussed in Sec. 2.2.2 but with a finer modelling taking into account way-points, the individual's departure time from a stay position can occur with a uniform probability between the pre-departure communication event and the post-departure communication event. The delay between these two events follows the distribution of the user's inter-event times, which can be assimilated to the population's inter-event times distribution for the sake of simplicity.

Mathematically, this means that departure bias follows a uniform distribution law bounded by the population's inter-event time distribution. Symmetrically, the same reasoning applies to the arrival bias. This modelling, focusing on the features of the whole population in terms of communication patterns instead of considering specific properties of the individual trip, represents a significant improvement with respect to the simple approach adopted in Sec. 2.2.2 to estimate the duration of trips from individual passively collected traces.

We modelled the population's inter-event time by considering an exponential law Z of parameter $\lambda$. While in the literature the inter-event time distribution is often modelled as a truncated power law distribution, we selected an exponential distribution for a matter of simplicity. It requires a single parameter $\lambda$ directly linked to the distribution average. Enhancing this modelling choice represents a possible future research direction.

The considerations above leads to:

$$ Z \sim \text{Exp}(\lambda) \qquad (3.8) $$

$$ Y|Z \sim U(0, z) \qquad (3.9) $$

Hence, the probability density function of Z, and the conditional probability density function of Y given the occurrence of the value z of Z can be written as:

$$ f_Z(z) = \lambda e^{-\lambda z} \qquad (3.10) $$

$$ \text{and } f_{Y|Z}(y|z) = \begin{cases} \frac{1}{z} & 0 \le y \le z, \\ 0 & \text{otherwise,} \end{cases} \qquad (3.11) $$

It can be proved (see [132]) that the probability density function of Y is:

$$ f_Y(y) = \lambda \int_0^{+\infty} \frac{e^{-\lambda(y+z)}}{y+z} dz \qquad (3.12) $$

and that the first two moments of Y are:

$$ E(Y) = \frac{1}{2\lambda} \quad \text{and} \quad V(Y) = \frac{5}{12} \frac{1}{\lambda^2} \qquad (3.13) $$

Those results characterise the random variable Y which models the departure and the arrival offsets. This gives for $X = 2Y$:

$$ \mu_X \equiv E(X) = 2E(Y) = \frac{1}{\lambda} \quad \text{and} \quad V(X) = 4V(Y) = \frac{5}{3} \frac{1}{\lambda^2} \qquad (3.14) $$

It is important to highlight here the importance of the size of $I_P^t$ to ensure that $\mu_X$ is representative of the average bias related to regional path P at time t: the larger the sample,

<footnote>would deserve further study, since mobile phone or social network usages are known to be correlated with socio-demographic characteristics.</footnote>
the more the theoretical average bias $\mu_X$ is representative of the sample's average bias. This is all the more important as with our bias model, as the variance of $X$ increases with mean inter-event time $E(Z) = \frac{1}{\lambda}$:

$$V(X) = \frac{5}{3}E(Z)^2 \qquad (3.15)$$

Eq. 3.15 shows that the larger the mean inter-event time is, the more scattered the trip bias distribution will be, and the more data per time slot and per regional path will be needed to ensure a reliable de-biasing process.

## Regional Speed Estimation

The estimation of regional speed dynamics relies on the computation of mean regional path travel times. Specifically, the individual regional travel time from Eq. 3.4 can be expressed as the fraction of the distance travelled on trip $i$ in $r$ (i.e., $L_{P,r}^i$) over the mean spatial speed of $i$ in region $r$ (i.e., $V_r^i$), as described in Eq. 3.17.

$$\forall i \in I_P^i, \quad T_P^i = \sum_{r \in P} T_{P,r}^i \qquad (3.16)$$

$$T_P^i = \sum_{r \in P} \frac{L_{P,r}^i}{V_r^i} \qquad (3.17)$$

Due to the temporal sparsity of the individual UAPD trajectories, the terms $L_{P,r}^i$, and $V_r^i$ are considered unknown.

Although vehicles may experience different local and instantaneous speeds over an area, their average speeds depend mostly on overall traffic conditions, and mainly on the accumulation (i.e., number of vehicles in the region). These average speeds show little scatter among individuals, and can be approximated by the mean spatial speed of all individuals travelling in the region. This observation has sustained the development of the MFD theory [144, 150] and motivates the choice of a network partitioning into sub-regions of consistent traffic dynamics. On this basis, it is possible to consider this equivalence for each trip $i$, being $t$ the time slot trip $i$ refers to, i.e.,

$$V_r^i = V_r^t, \quad \forall i \qquad (3.18)$$

where $V_r^t$ is the regional spatial mean speed at time slot $t$.

In Eq. 3.17, after summing on the $I_{t,P}$ trips, this gives:

$$\sum_{i=1}^{n_{t,P}} T_P^i = \sum_{i=1}^{n_{t,P}} \sum_{r \in P} \frac{L_{P,r}^i}{V_r^i} = \sum_{r \in P} \sum_{i=1}^{n_{t,P}} \frac{L_{P,r}^i}{V_r^i} \qquad (3.19)$$

Eq. 3.19 can easily be rewritten as follows:

$$\sum_{i=1}^{n_{t,P}} T_P^i = \sum_{r \in P} \frac{1}{V_r^t} \sum_{i=1}^{n_{t,P}} L_{P,r}^i \qquad (3.20)$$

$$n_{t,P} \bar{T}_P^t = \sum_{r \in P} n_{t,P} \frac{\bar{L}_{P,r}^t}{V_r^t} \qquad (3.21)$$

$$\bar{T}_P^t = \sum_{r \in P} \frac{\bar{L}_{P,r}^t}{V_r^t} \qquad (3.22)$$
Again, a significant advantage of this averaging process over the sample $I_P^t$ is that the characterisation of individual regional trip lengths $L_{P,r}^i$ for any individual $i$ becomes unnecessary. Instead, the sample mean value $\bar{L}_{P,r}^t$ turns out to be sufficient. On the condition that the sampling size is large enough, this can be replaced by its static estimate $\hat{L}_{P,r}$, drawn from the exogenous trip length matrix $\hat{L}$ described above:

$$ \bar{T}_P^t \approx \sum_{r \in P} \frac{\hat{L}_{P,r}}{V_r^t} \qquad (3.23) $$

At this stage, the computed distance matrix is used to express, through Eq. 3.23, a relationship between the average travel time along regional path $P$ at $t$, and the underlying, unknown mean traffic spatial speeds of the regions along $P$.

Although the actual average trip duration $\bar{T}_P^t$ is unknown, the knowledge of the average time bias $\mu_X$ allows to estimate it from the average trip duration, i.e., $\bar{T}_{P,obs}^t$, observed from the UAPD data at each time slot $t$ and for each path $P$. Based on Eq. 3.7, we thus get:

$$ \bar{T}_{P,obs}^t - \mu_X \approx \sum_{r \in P} \frac{\hat{L}_{P,r}}{V_r^t} \qquad (3.24) $$

Conversely, the constant distance parameters $\hat{L}_{P,r}$ can be drawn from the aforementioned estimated trip length matrix $\hat{L}$. $\mu_X$ is assumed known as well. $V_r^t$ are the only unknowns of the system. When applying in Eq. 3.24 the change of variable $x_r^t = 1/V_r^t$, we finally get the unbiased system:

$$ \forall t, \quad S^t = \{ \bar{T}_{P,obs}^t - \mu_X = \sum_{r \in P} \hat{L}_{P,r} x_r^t, \quad \forall P \}. \qquad (3.25) $$

In Eq. 3.25, we name $S^t$ the linear system composed of $|R|$ unknowns ($x_r^t, r \in R$) and as many equations as the number of regional paths observed during the reference time slot $t$. The UAPD data analysis and the parameters extracted from the trip length matrix allow to fully characterise the system, which can be rewritten in matrix notation as follows:

$$ \forall t, \quad S^t = \{ \mathbf{T}_{obs}^t - \mu_X = \hat{L}^t |x^t| \} \qquad (3.26) $$

where $\mathbf{T}_{obs}^t$ is the average observed travel time vector and $\hat{L}^t$ is the sub-matrix of $\hat{L}$ restricted to the regional paths observed at time slot $t$.

Given that the number of regional paths will generally exceed the number of regions of the adopted partitioning, $S^t$ is very likely over-determined. Consequently, the system will probably have no exact solution, but an approximated one can be calculated using regression analysis. To this purpose, we applied a non-negative least squares regression method to solve the system. In our case, the non-negative constraint allows for taking into account the non-negative nature of the zonal traffic speed.

Taking the reciprocal values of the solution vector $x_0^t$ gives the optimal speed vector $v_0^t$. This resolution process can be iterated throughout the whole studied time span to estimate the complete temporal speed trends.

## Empirical Evaluation

We performed an extensive evaluation of the proposed framework by leveraging a cleaned and map-matched dataset of GPS trajectories that we down-sampled and artificially biased to mimic real-world mobile phone data. The GPS data were provided by a leader European navigation system provider<fnref n="9" /> from October 2017 and September 2018 over the Greater Lyon

<footnote n="9">⁹https://be-mobile.com</footnote>
<figure><img src="image_39.png" /><figcaption>**Figure 3.2:** Maps of the regions partitioning the city of Lyon, France.</figcaption></figure>

area. The study area focuses on the urban area of Lyon and the neighbouring municipality of Villeurbanne, located inside Lyon's ring road. We partitioned this territory into 16 distinct regions, as displayed in Figure 3.2. The traffic variables were verified to be relatively homogeneous in each region [147]. The trips used in this study were extracted from five typical weekdays, i.e., from Monday, February 12, to Friday, February 16, 2018. Being GPS data rather limited in sample size, we artificially extended the size of the dataset by duplicating each trip 100 times. Moreover, as few trips were observed at night-time in our dataset, the time span selected for our evaluation is restrained to day-time hours, i.e., in-between 5 am and 8 pm. The data from the full month of February 2018 were used for the offline calculation of the trip lengths matrix $\tilde{L}$.

Firstly, each GPS trajectory was transformed into our trip representation (i.e., regional path, actual arrival time, and actual travel time), along with the trip id. This approach corresponds to introducing a spatial down-sampling to replace the precise track information with the regional path feature.

To reproduce the temporal sparsity of mobile phone data, we artificially introduced temporal biases based on a given value of the average inter-event time, as from our bias model. The value of the average inter-event time depends on the population observed or on the type of data chosen: for example, NSD datasets will display shorter inter-event times than CDR and LBSN. Therefore, we generated, for each day of data, five different down-sampled datasets, one per inter-event time value. The selected average inter-event time values were 4, 8, 12, 16, and 20 minutes, to cover a large range of average communication rates. The introduction of the bias was performed on a per-trip basis, by sampling the departure and arrival biases according to the probability density function defined via Eq. 3.12 using the average inter-event time values above. The actual travel time of each trip is thus biased with the sum of the sampled departure and arrival biases to obtain the observed travel time. Finally, the biased trips were grouped by regional paths and 15-minutes slots, averaging the observed travel times in each of the resulting groups thus obtaining the $\tilde{T}_P^t$ values. By using our methodology as from Eq. 3.25, we obtained a speed profile in kilometres per hour, per region, and per 15-minutes slots for each day of the evaluation.

To evaluate the reliability of the results, we compared the estimated regional speed profile to the spatial mean speed $V_r^t$ in region $r$ over $t$. The latter was computed from the original GPS trajectories as the ratio of the total travelled distance $TTD_r^t$ and the total travel time
<table><caption>**Table 3.1:** Average speed errors for each mean inter-event time (without de-biasing).</caption><thead><tr><th>Avg IET</th><th>MAE (km/h)</th><th>RMSAE (km/h)</th><th>MAPE (%)</th><th>RMSAPE (%)</th></tr></thead><tbody><tr><td>4</td><td>7.894771</td><td>18.576767</td><td>10.773218</td><td>21.648403</td></tr><tr><td>8</td><td>12.970168</td><td>31.347989</td><td>16.233049</td><td>33.484246</td></tr><tr><td>12</td><td>16.918677</td><td>41.133133</td><td>20.600466</td><td>42.693440</td></tr><tr><td>16</td><td>19.920772</td><td>48.522253</td><td>24.008746</td><td>49.748655</td></tr><tr><td>20</td><td>22.204453</td><td>54.189049</td><td>26.608480</td><td>55.214964</td></tr></tbody></table>

<table><caption>**Table 3.2:** Daily speed errors with de-biasing in the worst-case scenario (20-minutes IET)</caption><thead><tr><th></th><th>MAE (km/h)</th><th>RMSAE (km/h)</th><th>MAPE (%)</th><th>RMSAPE (%)</th></tr></thead><tbody><tr><td>Day 1</td><td>4.911907</td><td>7.174379</td><td>12.124498</td><td>15.264945</td></tr><tr><td>Day 2</td><td>5.057217</td><td>7.564531</td><td>12.024699</td><td>15.306312</td></tr><tr><td>Day 3</td><td>5.067316</td><td>7.505059</td><td>12.617445</td><td>15.839542</td></tr><tr><td>Day 4</td><td>4.848885</td><td>7.216612</td><td>13.006559</td><td>18.000039</td></tr><tr><td>Day 5</td><td>4.678821</td><td>6.591109</td><td>12.331655</td><td>15.714796</td></tr></tbody></table>

$TTT_r^t$ in region $r$ during $t$, as from Eq. 3.2.

Tab. 3.1 reports the performance metrics (i.e., Mean Absolute Error (MAE), Root Mean Absolute Error (RMSAE), Mean Absolute Percentage Error (MAPE), Root Mean Square Absolute Percentage Error (RMSAPE)) for each day as a function of the average inter-event time, prior to applying the de-biasing procedure. The errors are quite high and increase significantly with the average inter-event time, highlighting the importance of the de-biasing method. This table can be contrasted with the results in Tab. 3.2, where the worst-case scenario, corresponding to an average inter-event time of 20 minutes, is analysed after applying our bias removal procedure. For all metrics and all days, there is a noticeable improvement in performance indicators due to the de-biasing procedure. The errors remain below a 20% limit when considering the daily RMSAPE, which is a significant improvement compared to the average value of 55.2% observed in Tab.3.1.

Fig. 3.3 displays the smoothed speed estimation results for the five different average inter-event time values, again after bias removal. The speed curves for different average inter-event time values almost overlap, confirming the importance and effectiveness of our bias removal approach. The results also show that the framework is generally capable of reproducing speed trends and dynamics (i.e., good match with the blue curves in the figure), with particularly good results for many of the considered urban areas and almost all of the ring-road regions.

The less satisfactory results for some regions could be partially explained by the low quantity or quality of the original GPS data in those areas. This issue may be related to the presence of tunnels or a low penetration rate of the data provider in those areas, which inevitably reduces the representativeness of the average travel times, despite the artificial expansion of the dataset. In situations where the input data do not adequately describe traffic conditions, the bias removal mechanism is inevitably insufficient for producing accurate speed estimations. This highlights the need for more advanced processing solutions, such as those that will be described in Sec. 3.2. It is also worth noting that validating the obtained results on real NSD and CDR datasets remains a task for future work.

## 3.2 Trajectory Inference at Scale with Network Signalling Data

To address the challenge of accurately reconstructing mobility trajectories from real-world, massive individual NSD and fully exploit their potential for urban mobility, we proposed a framework named TRANSIT, extensively described in [151].
<figure><img src="image_40.png" /><figcaption>**Figure 3.3:** Speed dynamics after average bias removal. Ground truth speeds from the raw GPS traces are reported as blue curves.</figcaption></figure>

**TRANSIT** was conceived to explicitly leverage the high sampling rate of NSD and the repetitive nature of human mobility, by building on interesting approaches from the literature recently proposed to overcome location-related limitations of CDR and CDR+ datasets: the *Detect, Expand, Check and REmove* (DECRE) [152] and the *Cumulative Weighted Moving Average* (CWMA) [153, 154] frameworks.

DECRE is built on the principle of removing oscillations from CDR in order to reduce spatial uncertainty for enhanced human mobility modelling. However, although such an approach could effectively improve spatial accuracy by removing the noise deriving from the oscillations, the resulting trajectory is still bounded to the original location information from the cellular network (i.e., antennas coordinates), thus exhibiting large spatial uncertainty. As for CWMA, the approach leverages oscillations to infer with increased accuracy user locations. However, CWMA does not take into account the existence of high regularity in human movements, and consequently in mobile phone events, which TRANSIT exploits to enhance the spatio-temporal accuracy of the inferred mobility trajectories.

**TRANSIT** receives as input the set of NSD events of a mobile device $i^{10}$, denoted by

<footnote n="10">¹⁰As in Sec. 2.2, we make the simplifying assumption in the following that a device corresponds to a given</footnote>
<figure><img src="image_41.png" /><figcaption>**Figure 3.4:** Flowchart of TRANSIT.</figcaption></figure>

$\mathcal{T}^i = \{e_1^i, \dots, e_n^i, \dots, e_{N_i}^i\}$, where $e_n^i$ is the $n^{th}$ NSD event recorded for device $i$. Each NSD event is the result of a communication activity between a mobile device and a base station antenna of the telecommunication network, across all 2G, 3G and 4G technologies; it is defined as a tuple $e_n^i = (c_n^i, t_n^i)$, where $c_n^i$ is the antenna at location $l_n^i$ that handled the network event, and $t_n^i$ is the timestamp of the instant at which the event was recorded. The NSD events in a mobile phone trajectory $\mathcal{T}^i$ are ordered by their timestamps $t_n^i$, and $N_i$ denotes the number of events for device $i$. Then, TRANSIT processes $\mathcal{T}^i$ in two successive phases:

* *Trajectory identification.* The framework labels each NSD event $e_n^i \in \mathcal{T}^i$ as either static, if the user $i$ is deemed to be engaged in an activity at a same location at the event time $t_n^i$, or mobile, if $i$ is performing a movement at $t_n^i$. The labelling factually allows telling apart the continuous time intervals during which an individual is moving or not, and building a set $\mathcal{A}^i$ of static activity sessions and a set $\mathcal{M}^i$ of mobile sessions. As a result, the set $\mathcal{M}^i$ also identifies all the trajectories, i.e., continued sequences of movement in time, of user $i$.

* *Trajectory augmentation.* The framework enhances the trajectories associated to mobile sessions in $\mathcal{M}^i$, by exploiting the fact that the same individual typically performs many trips between two given locations over time, generally following very similar paths. This creates redundancy in the mobility information that can be used to in-

<footnote>individual.</footnote>
crease the spatio-temporal accuracy of the trajectories. The resulting set of mobile sessions possibly augmented trajectories is denoted as $\hat{\mathcal{M}}^i$.

Ultimately, the output of TRANSIT are the set $\mathcal{A}^i$ of static activity sessions of user $i$, and the set $\hat{\mathcal{M}}^i$ of mobile sessions with augmented trajectories. Fig. 3.4 presents a flowchart of the stages of TRANSIT.

### 3.2.1 Trajectory identification

As anticipated, the trajectory segmentation step is applied to the individual set of NSD events $\mathcal{T}^i$ recorded for device $i$, and returns a subset of $\mathcal{T}^i$ where each event is labelled as static or mobile and detected oscillations are removed.

Fig. 3.5 illustrates the process of trajectory identification using TRANSIT. The interpolation of NSD events is portrayed as the black solid line. Fig. 3.5(a) refers to static antennas with daily accumulated association time above $T_w$. Fig. 3.5(b) identifies static activity sessions as obtained from consecutive sequences of static antennas only, and detected oscillations. Fig. 3.5(c) exhibits the final static activity sessions upon removal of oscillations, as well as the consequent detected mobile sessions.

We start by assuming that the time spent by user $i$ at the antenna $c_n^i$ associated to event $e_n^i$ is $t_{n+1}^i - t_n^i$, i.e., the temporal span to the subsequent event $e_{n+1}^i$. Given the high temporal resolution of NSD, this simple approach already provides a very good estimation of the time the user is associated to a given antenna, at a low computational cost. Then, a preliminary labelling is performed to trim down candidate static events. To this end, we calculate the cumulated time spent by user $i$ at each antenna $c_n^i$, on a daily basis. As devices stay connected to a limited set of antennas while still, we expect such antennas to yield a non-negligible cumulated time during the target day. We thus tag as static antennas for user $i$ those antennas with a daily cumulated time above a threshold $T_w$. In our experiments, we set $T_w$ to 20 minutes, which falls within the range of commonly accepted values for the typical minimum duration of a significant activity carried out by an individual at a same location [131, 155], and is employed also with high-frequency longitudinal (e.g., GPS) data [156]. An example is provided in Fig. 3.5(a).

A continued time interval where all events of device $i$ are only associated with static antennas is then denoted as a static activity session $a_k^i$. Typically, during one day, a user can have several static sessions, and each can be composed of one or multiple antennas. The set of all such sessions across the whole observation period is $\mathcal{A}^i = \{a_1^i, \dots, a_{K_i}^i\}$. It is worth highlighting here an important aspect concerning the spatial extent of the static activities detected by TRANSIT: while the approach proposed in Sec. 2.2 for OD matrix estimation detects stationary activities at the level of EDR sectors, TRANSIT detects static sessions at the level of the coverage area of one (or small groups of) antenna(s), which, in urban areas, corresponds to an extent from one to three orders of magnitudes smaller than the average area covered by EDR sectors. This represents a significant advancement in terms of the spatial resolution adopted with TRANSIT. This also represents an improvement with respect to the mean speed estimation framework described in Sec. 3.1. That framework relies on a regional partitioning of the road network, where regions are defined by homogeneous traffic analysis zones. These zones typically cover a significantly larger spatial extent than the coverage area of antennas, especially in dense urban environments.

TRANSIT includes an additional step, addressing the removal of oscillations from the mobile phone traces, to enhance the estimation of the duration and location of the detected static sessions. After the stage above, only part of the antennas are labelled. Unlabelled antennas are either encountered during movements, or the result of oscillations that are known to characterise mobile device association to the radio access infrastructure [156]. Oscillations can in fact affect both static and mobile users. In the former case, they can cause the separation of continuous static activities into different static sessions in $\mathcal{A}^i$ interleaved
<figure><img src="image_42.png" /><figcaption>**Figure 3.5:** Main steps of the trajectory identification via TRANSIT.</figcaption></figure>
by non-static antennas. In order to address the issue, and remove oscillations from $\mathcal{A}^i$, TRANSIT adopts the following heuristic. If (i) two consecutive static sessions $a_k^i$ and $a_{k+1}^i$ present at least one common (static) antenna, and (ii) the number of unique antennas associated to events observed after $a_k^i$ and before $a_{k+1}^i$ is below a threshold $N_o$, we merge all the events in $a_k^i$ and $a_{k+1}^i$ into a new, single static session. The new sessions replaces the former pair in $\mathcal{A}^i$. An example of oscillation detection and static sessions before the merging process is shown in Fig 3.5(b).

The single events identified as oscillations in the previous stage are in fact removed from $\mathcal{T}^i$ entirely, so as to limit uninformative noise in the data. The revised static sessions in $\mathcal{A}^i$ are further filtered based on their total duration, and only those with a time span higher than a threshold $T_s$ are retained. The value of $T_s$ corresponds to the assumed minimum duration of a static activity, so that we do not include, e.g., waiting periods at red traffic lights for pedestrian or vehicular trips, or dwell times at stops for bus trips. For the same reasons explained above in relation to threshold $T_w$, used to identify static antennas, the value of 20 minutes has been adopted for $T_s$ as well.

TRANSIT also enforces consistency in the locations of events associated to static activity sessions, as follows. First, we compute the centroid of the locations $l_n^i$ of all events in each session $a_k^i$; then, the well-known DBSCAN clustering algorithm<fnref n="11" /> is run on the centroids of all $a_k^i \in \mathcal{A}^i$. This lets us group together all static sessions related to a same activity, and compute a consolidated location for the activity as the barycentre of all centroids in a same cluster. The locations $l_n^i$ of all events in each session $a_k^i$ are then replaced with the barycentre of the corresponding cluster. Note that the position of the static activity sessions that are labelled as outliers by the DBSCAN algorithm are left unchanged. An example of the resulting $\mathcal{A}^i$ is in Fig 3.5(c).

Finally, all events that have not been labelled as static are labelled as mobile. This directly identifies the mobile sessions $m_h^i$ of user $i$, as the time-continuous sequences of mobile events; an important remark is that the two static events immediately preceding and following the mobile session are also integrated into $m_h^i$ to get a more precise indication of the start and the end locations of the corresponding trip. As a result, the set of mobile sessions is $\mathcal{M}^i = \{m_1^i, \dots, m_{H_i}^i\}$. Each $m_h^i$ corresponds to one trajectory of user $i$ identified by TRANSIT. An example with two detected mobile sessions is reported in Fig 3.5(c).

### 3.2.2 Trajectory Augmentation

The sequences of NSD events in $\mathcal{T}^i$ that correspond to the single trajectories $m_h^i$ of user $i$ are still affected by the limited spatial accuracy that characterises NSD events as explained in Sec. 1.1.2. In its second phase, TRANSIT thus aims at improving the geographical quality of the movement information. As anticipated, the framework relies on the regularity of human mobility; more precisely, we use the information from multiple similar trajectories identified for a same user to mutually improve their accuracy.

As a first step, a similarity measure is computed for all pairs of mobile session $m_h^i \in \mathcal{M}^i$. We employ the Hausdorff distance [157], which is defined as:

$$d_H(m_{h_1}^i, m_{h_2}^i) = \max\{D(m_{h_1}^i, m_{h_2}^i), D(m_{h_2}^i, m_{h_1}^i)\}, \quad (3.27)$$

$$\text{where } D(m_{h_1}^i, m_{h_2}^i) = \sup_{l_{n_1}^i \in m_{h_1}^i, l_{n_2}^i \in m_{h_2}^i} \inf d(l_{n_1}^i, l_{n_2}^i), \quad (3.28)$$

where $m_{h_1}^i$ and $m_{h_2}^i$ are the two mobile sessions to be compared and $d(\cdot, \cdot)$ is the geodesic distance between the two argument locations. This results in a matrix of pairwise distances between all mobile sessions of a same user $i$.

<footnote n="11">¹¹Details on the parameters used with the DBSCAN algorithm for static session clustering can be found in [151].</footnote>
<figure><img src="image_43.png" /><figcaption>**Figure 3.6:** Trajectory reconstruction for a given pair of Origin-Destination path.</figcaption></figure>

Then, DBSCAN is applied to the distance matrix, in order to group trajectories that have similar spatial geometries, and correspond to diverse trips of the user between the same two static activity locations. Fig. 3.6(a) shows an example of a set of mobile sessions, i.e., trajectories, grouped together in the same cluster by DBSCAN, for the origin-destination activity locations in Fig. 1.1. Based on the result of DBSCAN, we can tell apart the mobile sessions in $M^i$ into two subsets: (i) trajectories that fall into a cluster, i.e., which refer to a path that is recurrent in the mobility of user $i$, and which we denote as the set $M_R^i$; and, (ii) outlier trajectories that represent unique movements of $i$, which are grouped in set $M_O^i = M^i \setminus M_R^i$.

For trajectories in $M_R^i$, TRANSIT operates a spatial augmentation, as follows. First, the average duration is computed for all trajectories assigned to the same spatial cluster by DBSCAN above; this corresponds to the expected time that user $i$ takes to travel between the same origin-destination activity locations. The time information is used to filter out trajectories whose duration deviates from the median by 50% or more: these mobile sessions are considered not representative of the routine mobility patterns along the target path. The retained trajectories in the same cluster are then temporally scaled (i.e., stretched or compressed) in time so as to match the average travel duration for the cluster. Finally, the scaled trajectories are temporarily binned according to a fixed time period of one minute, and the spatial coordinates of all different events that fall in the same time bin are averaged.

The previous steps lead to a set of positions, one per minute, which represent the reconstructed itinerary. If there is no event within a particular time slot, the resulting enhanced trajectory will have missing positions. All trajectories in the cluster are then matched to the reconstructed one and become thus identical in the space dimension. However, they are re-conducted to their original duration (i.e., via compression or stretching) so as to keep them faithful to their recorded travel time in the NSD.

As a result, each original mobile session in $M_R^i$ is replaced by a set of reconstructed positions without any temporal deformation and is enriched with information derived from multiple similar trajectories travelled by the same user. This set of enhanced mobile sessions is referred to as $\tilde{M}_R^i$. Fig. 3.6(b) shows the final spatial trajectory inferred from the cluster in Fig. 3.6(a). Trajectories in $M_O^i$ stay instead unchanged, corresponding to those obtained from the simple interpolation of NSD data. The final set of mobile sessions is $\tilde{M}^i = M_R^i \cup M_O^i$.

### 3.2.3 Performance Evaluation

The performance of TRANSIT has been evaluated by using ground-truth information on the trajectories of a small set of four volunteers. More precisely, we collected a total of 310 high-resolution trajectories using a high-frequency GPS logger running on their smartphones; also,
we recovered the NSD data generated by the volunteers' devices (which used Orange as their network provider) during the same observation period. We stress the fact that, although reduced in size, this validation dataset is among the few allowing a direct comparison of NSD and GPS data. The details of the validation process are reported in the reference paper [151], jointly with extensive sensitive analyses of the selected parameters as well as a thorough investigation of the impact of mobility recurrence on the accuracy of the results.

Here, we report the key figures obtained with TRANSIT on the validation dataset, relative to the two identified state-of-the-art solutions (i.e., DECRE and CWMA). We used the classical precision, recall and F1 metrics to evaluate the performance of the trajectory segmentation approaches. Formally:

$$ \text{Precision} = \frac{TP}{(TP + FP)}, \quad \text{Recall} = \frac{TP}{(TP + FN)} \quad \text{and} \quad F1 = 2 \cdot \frac{\text{Precision} \cdot \text{Recall}}{(\text{Precision} + \text{Recall})} \qquad (3.29) $$

where (i) the number of true positives TP is the number of NSD events labelled as static when the user is also labelled as static in GPS data; (ii) the number of false positives FP represents the number of NSD events labelled as static while the user is in fact mobile according to the ground truth; (iii) the number of false negatives FN maps to the number of NSD events labelled as mobile while the user is static in the GPS data. Specifically, TRANSIT yields a 5% relative improvement, in terms of the F1 metric, with respect to CWMA in the trajectory identification task (we remind that DECRE does not include any classification procedure), with nearly perfect recall of static sessions.

More importantly, concerning trajectory enhancement, TRANSIT attains an average geographical error as low as 190 meters, bringing it down to 130 meters for some users on average, with global improvements in the order of 30% and 41% with respect to CWMA and DECRE, respectively, when used on the same NSD datasets. The improvement raises to 60% when CWMA is used on the CDR+ version of the dataset, and 72% when DECRE is used on the CDR one. These results were obtained using the average of two distance metrics computed on a per-trajectory and per-user basis as follows:

$$ D_{GPS} = \frac{1}{|m_{GPS}|} \sum_{e_{n'} \in m_{GPS}} \min_{e_n \in m_{NSD}} d(l_{n'}, l_n) \qquad (3.30) $$

$$ D_{NSD} = \frac{1}{|m_{NSD}|} \sum_{e_n \in m_{NSD}} \min_{e_{n'} \in m_{GPS}} d(l_n, l_{n'}) \qquad (3.31) $$

where $m_{GPS}$ and $m_{NSD}$ are two corresponding trajectories inferred for a given user from GPS and mobile network data, respectively. The operator $| \cdot |$ denotes the cardinality of the argument set, i.e., the number of samples in the case of a trajectory, and $d(\cdot, \cdot)$ the geodesic distance. We use both metrics as they are complementary: while $D_{GPS}$ is representative of the error observed for continuously tracked user, $D_{NSD}$ measures the error specific to events recorded by the mobile phone network.

The takeaway is that the relative improvement of TRANSIT does not simply derive from the increased temporal resolution of NSD, but also by its capacity to enhance recurrent trajectories via historical data.

### 3.2.4 Large-Scale Applications of TRANSIT

While the validation results were related to a reduced number of users, the interest of TRANSIT reveals at city-wide scales, where it can enable a number of mobility-related applications. In the following, we thus focus on four case studies of TRANSIT related to urban mobility that fully leverage the large-scale NSD datasets $\mathcal{D}_P$ and $\mathcal{D}_L$ provided by Orange Innovation for the cities of Paris and Lyon, respectively, and introduced in Sec. 1.1.2. We stress the fact that such results have been obtained by means of specific optimisations
<figure><img src="image_44.png" /><figcaption>**Figure 3.7:** Average weekly profiles of the number of concurrent trips in (a) Paris and (b) Lyon, as inferred from TRANSIT and smart card data. Normalised versions of the same profiles are presented in (c) and (d).</figcaption></figure>

required in order to process the three months of NSD from the large-scale datasets $D_P$ and $D_L$. In particular, TRANSIT has been implemented in Spark by exploiting the Map-Reduce distributed programming paradigm.

## Urban Mobility and Public Transport

By counting the number of concurrent active trips inferred via TRANSIT over time, we were able to reconstruct accurate temporal profiles of the travel demand in urban regions. For such profiles to be dimensionally correct, a re-scaling procedure was needed to account for the penetration rate of the technology (close to 100% in developed countries like France) and the market share of Orange (at 37% over the French territory). The resulting average weekly demand profiles computed in Paris and Lyon are depicted in blue in Fig. 3.7(a) and Fig. 3.7(b), respectively. Our estimates are that around 1,300,000 individual trips occur at the same time in Paris during commuting peaks, while the figure is at 180,000 for Lyon.

We compare the profiles obtained with TRANSIT with equivalent ones from smart card data, which capture mobility via public transportation systems. For Paris, data were provided by the transportation company Ile-de-France-Mobilité. Concerning Lyon, data were shared by the transportation company Keolis-Lyon. For both cities, public transport data were provided in the same period of the year of NSD, and all smart-card transactions were anonymised in the form of aggregate measures at the scale of the whole agglomeration.

Also in this case, a re-scaling was required: while the TRANSIT trajectories refer to the resident population, the smart card data include both residents and non-residents. In order to make the numbers comparable, we applied a scaling factor of 0.81 to the smart card temporal profile; the factor has been calculated from the raw network signalling data, by computing the average instantaneous fraction of resident subscribers present in the target cities, over the total number of observed users. The weekly profiles from smart cards are su-
perposed to the TRANSIT-inferred ones, as the orange curves in Fig. 3.7(a) and Fig. 3.7(b).

The comparison of the profiles reveals interesting facets of mobility in Paris and Lyon. Clearly, the volume of trips identified by TRANSIT is higher than that reconstructed with smart card data: NSD allows monitoring virtually all transport modes, including those beyond public means, e.g., private vehicles, biking, or walking. This lets us quantify which proportion of trips is performed by underground, buses or tramways, and which using personal means. We find that a significant fraction of trips is performed using public transport in both cities: we estimate the percentages of movements captured by smart card data to 66% and 39% of the total, in Paris and Lyon, respectively. The difference between these values is explained by the more developed multi-modal transit network available in Paris, as required to support mass mobility in such a large metropolis.

In addition, we can investigate the temporal incidence of public transport by looking at versions of the same profiles that are normalised so that the integral of all curves is one. Fig. 3.7(c) and Fig. 3.7(d) show the result. This perspective lets us appreciate how in Paris public transport is especially important during commuting hours, but relatively less used during the lunch break or weekends. A slightly different pattern emerges in Lyon, where public transport is also intensely used around midday, but has a lower incidence during evenings and Sunday mornings. We highlight that obtaining this type of insight is hardly achievable by solely relying on surveying, which demonstrates the value of NSD and a method like TRANSIT that can exploit them.

As a final remark, we highlight that the results in Fig. 3.7 can also be considered as a partial validation of the trajectories inference mechanism performed by TRANSIT in large-scale settings. Indeed, the near-perfect match of the timing of commuting peaks or overnight low mobility among curves proves the capability of our trajectory segmentation approach to identify trips that are very consistent with data collected in the field over time. It is also worth underlining that, differently from the results reported in Sec. 2.2.2 and Sec. 3.1, no de-biasing procedure was applied in relation to the start time of the mobility session or the land use. The fact that the profiles from the two different sources align so well, without any de-biasing mechanism, further proves the relevance of 4G communication events in the dataset and the fundamental contribution of the processing taking place in the trajectory identification and enhancement phases of TRANSIT.

## Popular Paths of Commuting Trips

As a second application, we focused on inferring popular commuting trips within a city. The knowledge of such trips is an extremely precious source of information for transport authorities and city planners as: (i) they represent the largest share of the daily urban traffic demand of a city; (ii) they identify the typical commuting behaviours of travellers which regularly stress the transport network infrastructure, especially during peak hours; (iii) they are hard-to-quantify and characterise at city-scale because of the absence of dedicated sensors or probes that can precisely capture the multi-modal, diverse and time-varying nature of such trips.

By applying our framework to the large-scale datasets $D_P$ and $D_L$, we associated to each user *i* of the two analysed cities a set of trips $\tilde{M}^i$ for the whole period of 3 months. As explained above, $\tilde{M}^i$ can be divided into two subsets: $\tilde{M}_r^i$, a subset of recurrent trips enhanced by TRANSIT, and $M_o^i$, a set of non-recurrent trips of user *u*. Considering that commuting trips are, by definition, recurrent, only subset $\tilde{M}_r^i$ was retained in this analysis. Furthermore, to extract commuting trips from $\tilde{M}_r^i$, we filtered only those trips associated with the two most popular locations of each user, under the constraint that at least 10 trips were present between these two locations. The underlying assumption was that the remaining set should mostly contain the two most popular trips performed by users in their daily routine, i.e., home-to-work and work-to-home trips (commuting trips).

The spatial density ( heatmap ) of the reconstructed trips is represented in Fig. 3.8(a)
for Lyon and Fig. 3.8(b) for Paris. As a first consideration, the recurrent trips appear to
have overall a good match with the multi-modal urban transportation network, graphically
overlapped to the heatmap in Fig. 3.8(c) for Lyon and in Fig. 3.8(d) for Paris. A more
in-depth inspection of the figures highlights that, for both cities, the subway network, the
tramway lines and most important urban roads clearly show up among the commuting trips
reconstructed via TRANSIT. In the case of Paris, NSD trips appear to have a near-perfect
match to the underlying multi-modal transport network. The less evident match for the case
of Lyon, especially characterising some peripheral roads (however present in the heatmap),
can be explained by the lower number of available trips and the lower density of the cellular
network of Lyon in these areas, compared to those from the capital city.

Of course, the fact that the majority of commuting trips maps to the public transporta-
tion network is not unexpected. However, TRANSIT opens the door to a detailed analysis
of these trips: the obtained trips can be map matched to the different transportation lines
and modes (see Sec. 3.3), showing their share of trips, on different days of the week and
at different times of the day. Such information would be highly valuable for any public
transportation company or municipality to perform, e.g., informed planning and efficient
resource deployment. Similarly, from a resilience perspective, highly congested routes or
transit lines could be easily spotted, even in real-time settings, for network monitoring and
vulnerability assessment.

<figure><img src="image_45.png" /><figcaption>**Figure 3.8:** Heatmap of commuting trips in Lyon and Paris.</figcaption></figure>
Human Mobility Analysis during Abnormal Events

As a third application, we used TRANSIT to detect abnormal mobility situations that
can occur in the city. In the perspective of resilient management of urban mobility such
knowledge could prove to be essential to learn how people react to specific kinds of abnormal
urban situations, such as accidents, sports events, and celebrations. For this, we segmented
the city of Paris into a grid with squares of dimension 800m × 800m, with a temporal
bin size of one hour. This spatio-temporal granularity makes it possible to analyse human
mobility at a fine-grained scale. For each zone, we computed the attraction demand profile,
which corresponds to the number of trips having as destination the studied zone at any given
hourly time slot. These profiles were obtained by retaining such trips from the whole set of
trajectories $\bar{M}^i$ computed via TRANSIT on $D_P$ for each user $i$. This allowed us to build a
typical weekly attraction profile for each zone and, at the same time, to distinguish abnormal
patterns during certain events. We considered three such abnormal mobility situations,
reported as examples below.

<figure><img src="image_46.png" /><figcaption>**Figure 3.9:** Typical/atypical weekly temporal demand profile during atypical events.</figcaption></figure>

First of all, on Wednesday, the 1st of May 2019, a bank holiday, the Labour Day march took place near Place d'Italie in Paris. Fig. 3.9(a) shows in blue the typical attraction profile of this zone and in red the attraction profile of the week that includes the demonstration. Whereas the attraction profile was similar to the typical profile for most days, we can see that, on Labour Day, the attraction of the studied zone presents an unexpectedly high peak after midday.

As a second event, we studied the fire of Notre Dame de Paris cathedral, on Monday the
15th of April 2019. Fig. 3.9(b) shows in blue the typical attraction profile of this zone, and in red the attraction profile of the week that includes the abnormal event. We can see a high peak in the attraction profile right after the beginning of the fire on Monday 15th (around 6:30 pm). Contrary to the previous example, this event also affected mobility the following days, when an attraction demand higher than usual is observed in the corresponding area. The attraction profile progressively decreases after the event, but an upsurge on Sunday, the Easter holiday, can be noticed, probably explainable by religious activities and nearby gatherings of tourists and worshippers visiting the area surrounding the cathedral after the fire on this special day.

Finally, we studied another special event, the Marathon of Paris, on Sunday the 14th of April 2019, with its start and end in the proximity of the Arc de Triomphe. A high peak on the attraction profile can be observed at the departure time of the marathon, at 9 am, as shown in Fig. 3.9(c). A second peak is observed a few hours later, more spread over time and lower in magnitude compared to the first one, corresponding to the marathon arrival.

These three examples are representative of the vast potential of TRANSIT towards building mobility profiles of the typical demand attracted by (or emitted from) a given zone, as well as detecting and characterising mobility patterns during abnormal or special events. These ideas were at the core of the real-time anomaly detection approach developed in the context of the ANR project DISCRET and discussed in [158].

## Ring Road Trajectory Analysis

As a fourth application, we leveraged TRANSIT to perform a fine-grained trajectory analysis focused on the Paris *péripérique* (ring road). The mobility flow on this urban highway is usually very high, often leading to heavy congestion, especially during peak hours. Transport authorities are traditionally very interested in the possibility of tracing and quantifying the flows of people moving along city major road axes. Such studies are necessary for urban planning purposes, infrastructure renewal and road maintenance, and can be extremely cost-demanding. They are normally based on travel diaries or GPS trajectory collection and generally end up capturing only a small sample of the flow actually traversing the major axis, with resulting limited accuracy. TRANSIT allows leveraging NSD to access a much larger and more representative sample of this specific population.

In our case study related to the Paris *péripérique*, we considered four different zones of interest: the east, west, north and south entries. The idea was to select a spatial zone and study all the trajectories passing by the respective zone. The enhanced trajectories $\tilde{M}^i$ produced by TRANSIT on $D_P$ allowed us to capture at scale the origin, the destination, and the paths taken by the users passing by the studied zone, the kind of information usually expected in the aforementioned studies. The result for the four zones of the *péripérique* (north, west, south and east) are reported in Fig. 3.10(a), Fig. 3.10(b), Fig. 3.10(c) and Fig. 3.10(d).

The obtained maps underline the major role of the *péripérique* in Paris, allowing people to travel across the city and reach any area of interest. Some interesting patterns can be distinguished as well. For example, the trips coming from the west side of the city show a strikingly different pattern than the three other maps, with trajectories that appear less scattered towards the centre of Paris and concentrated along a few major axis. This can be explained by the fact that the west side of Paris is the richest area of the city, with inhabitants who have a lifestyle involving shorter commuting trips. Moreover, the west side of Paris is also the area with the highest density of offices, including the La Defense and Boulogne neighbourhoods. This could explain why this area attracts a large amount of trips, ever from faraway zones.

These results hint at the numerous perspectives brought by TRANSIT in the study of major road arteries. These include fine-grained temporal analysis, the detection of usage and attraction patterns, origin and destination profiling, etc. Generally speaking, having access to
<figure><img src="image_47.png" /><figcaption>**Figure 3.10:** Heatmap of recurrent trips for the Paris ring road (the black square shows the catchment area)</figcaption></figure>

detailed human mobility trajectories at scale, such as those produced by TRANSIT, enables
the in-depth study of any part of the transportation network.

## 3.3 Multi-modal Path Reconstruction via Hidden Markov Model

The high spatial accuracy of the outputs produced by TRANSIT naturally raises the re-
search question of inferring the precise route taken by a mobile subscriber across the multi-
modal transportation network of a city. In other words, the goal was to reconstruct the
sequence of road segments, public transport stops or railroad paths followed by the observed
user during a trip. This challenge, known as multi-modal map-matching of sparse spatio-
temporal tracks, allows for achieving an extremely high level of detail in reconstructing
<figure><img src="image_48.png" /><figcaption>**Figure 3.11:** Illustration of Hidden Markov Model based map-matching.</figcaption></figure>

mobility trajectories paving the way to precise multi-modal pattern detection and traffic variable estimation. In our specific case of NSD trajectories, we need to project the mobile phone trace, potentially improved through TRANSIT, onto a multi-layer graph used to represent the multi-modal transportation network of the target city. The problem is exemplified in Fig. 3.11(a) from [159]: an ordered sequence of timestamped points ($o_1, o_2, o_3$) corresponding to observations of a moving user has to be matched to a sequence of nodes, e.g., ($v_1, v_5, v_4, v_6$), defining a path on the underlying graph-based representation of the transportation network. We addressed such research problem in [151, 160] by relying on Hidden Markov Model (HMM)-based map-matching.

HMM-based map-matching represents the state-of-the-art for noisy and sparse location data and, specifically, for mobile phone data traces [161]. A HMM can be defined as a 5-tuple $\langle S, O, I, T, R \rangle$, with $S = \{s_0, \dots, s_{N-1}\}$ representing a finite set of states; $O = \{o_0, \dots, o_{M-1}\}$ corresponding to a finite set of observations; $I$ being the probability distribution of the initial states; $T$ representing a set of transition probability. The probability to transit from hidden state $s_m$ to hidden state $s_n$ is denoted as $t(s_m, s_n)$. Finally, $R$ is a set of emission probability. The probability to emit observation $o_k$ from hidden state $s_m$ is denoted as $r(s_m, o_k)$. The matrices $I, T$ and $R$ represent the main hyper-parameters of the HMM. Fig. 3.11(b) shows a possible HMM modelling of the toy example reported in Fig. 3.11(a).

In this context, we modelled the multi-modal transport network under study as a weighted directed multi-layer graph [56], $G = (V, E, L)$, with $L = \{l_q\}_{q=1}^Q$ the set of elementary layers, each representing a specific transport mode $q$ (among a total of $Q$ different options), $V$ the set of nodes and $E$ the set of edges. Each layer $l_q \in L$ contains a node subset $V_q \subset V$ and an edge one $E_q \subset E$. Each $V_q$ corresponds to the set of road intersections, transit/train stops, etc., associated to mode $q$. Each $E_q$ represents the set of connections (road segments, public transport lines, etc.) between the nodes of the same transport mode $q$. Edges are directed according to the allowed direction of movement on the corresponding segment and weights can represent, e.g., travel time, length, capacity or speed information associated to
the segment. In the general model, the edge set $E$ also includes a set of inter-layer edges $E(l_{q_i}, l_{q_j})$ that allows to connect pairs of layers $l_{q_i}, l_{q_j}$, i.e., to switch between the two different transport modes $q_i$ and $q_j$. In our study, we assumed graph $G$ to be composed only of layers related to road, train and public transport, denoted as $G_{road}$, $G_{train}$ and $G_{pt}$ in the following. It is worth noting that, with respect to the general network modelling approach, TRANSIT allows the simplifying assumption to consider layers not connected to each other, i.e., $E_{l_{q_i}, l_{q_j}} = \emptyset$. The rationale behind this assumption is the following: in order to move from one transport mode to another the user will need to spend a period of immobility at a given location (e.g., in the proximity of a bus stop, i.e., under the coverage of a limited subset of cellular antennas). If long enough, TRANSIT detects this waiting time as a static activity of the user. Therefore, a multi-modal trip over these three networks will actually be split in multiple sub-trips by TRANSIT, that can be each separately matched to the specific sub-graph. Inter-modal trips are therefore not possible across the three main sub-graphs $G_{road}$, $G_{pt}$ and $G_{train}$, but can occur within the public transport network, i.e., $G_{pt}$, which covers indeed three different transportation modes (i.e., subway, tramway and bus) and is therefore a multi-layer graph itself.

In [160], we extended the TRANSIT framework with map-matching capabilities to further enhance the spatial accuracy of its output: map-matching is performed on the whole set of mobile sessions produced by TRANSIT, i.e., $\tilde{M}^i$, which includes both raw signalling trajectories ($M_O^i$) and TRANSIT-enhanced recurrent trajectories ($\tilde{M}_R^i$). Specifically, we used the following HMM formulation. Hidden states are modelled as the set of vertices (nodes) of the generic transportation sub-network $G_j$. Emissions are modelled as the unique set of x-y coordinates in $\tilde{M}^i$. This set is composed of antennas coordinates from cellular network in $M_O^i$ and reconstructed positions in $\tilde{M}_R^i$. Ultimately, given a sequence of observations, i.e., sequence of antennas for $M_O^i$ and reconstructed positions for $\tilde{M}_R^i$, the HMM is used to retrieve the most likely sequence of hidden states, i.e., sequence of nodes on the transportation sub-network $G_j$.

In the related literature, a similar problem has been solved in [162, 163, 161]. However, these approaches have two main drawbacks. They leverage supervised machine-learning solutions that require a large amount of labelled cellular trajectories for training the parameters of the models. Ground-truth mode-labelled trajectories are very hard to obtain, especially when dealing with highly dynamic and complex environments, such as urban areas. Moreover, most of the approaches match the cellular trajectories only to road networks, without considering other sub-networks corresponding to alternative transportation modes, such as tramway, subway, bus, etc.

Among the very few exceptions, it is worth mentioning the methodology proposed by Asgari et al. [164]. Their solution, namely CT-Mapper, relies on an unsupervised HMM model, which aims at mapping sparse cellular trajectories<fnref n="12" /> to a multi-layer transportation network. Similarly, in our previous work [151], we studied the application of unsupervised HMM-based map-matching for solving the same problem in more realistic settings, characterised by the usage of real-world raw NSD and the analysis of a dense urban case study.

Our approach builds on such previous work and solves the HMM-based map-matching problem according to a two-steps map-matching procedure. The first phase consists in an optimised Viterbi algorithm [165]. The inputs of the Viterbi procedure are the following: the generic transportation sub-network modelled as a (multi-layer) graph $G_j$, the possible hidden states (set of the nodes of $G_j$), the emissions (the unique set of x-y coordinates in $\tilde{M}^i$), the HMM hyper-parameters (matrices $I, T$ and $R$) and the input trajectory from $\tilde{M}^i$. By calculating all possible paths given the input trajectory, the output of the Viterbi algorithm is the most likely sequence of graph nodes, one for each time instant in the input. For real-time application, due to a large number of states and emissions, execution time is

<footnote n="12">¹²In their work, Asgari et al. used a variant of CDR data, characterised by localisation information produced at a frequency, much higher than that observed with traditional CDR.</footnote>
critical. To improve performance, similarly to [163], we implemented an optimised version of the Viterbi algorithm that eliminates all multiplications by zero thus reducing the search space by keeping only states that can be emitted from each observable state. After inferring the most likely states sequence using the optimized Viterbi implementation discussed above, the output sequence of hidden states (nodes on a given sub-network $G_j$) does not necessarily produce a connected path on the specific transport sub-network. Therefore, as the second phase of the map-matching procedure, the final trajectory is further completed by applying a traditional weighted shortest path detection algorithm on the underlying transportation graph between any two consecutive nodes of the most likely sequence. The final completed sequence of nodes on sub-network $G_j$ represents the map-matched trajectory for the processed trace from $\hat{M}^i$ for user $i$.

### 3.3.1 Empirical Evaluation in the Lyon Case Study

In the following, we briefly describe the parametrisation of the HMM, as well as the main results of the proposed solution that has been applied to the $D_L$ dataset of the city of Lyon. More information about the pre-processing steps, parameter choices and detailed results can be found in [160].

#### Initial Probability

All the nodes of the transportation network are initially equally assigned with a probability of $1/N$ with $N$ representing the total number of nodes in the transportation (sub-)network:

$$ \pi(s_m) = \frac{1}{N} \qquad (3.32) $$

#### Transition Probability

The transition probability corresponds to the probability that a mobile phone user moves on the underlying transportation (sub-)network from hidden state $s_m$ at time $t-1$ to hidden state $s_n$ at time $t$. As in [166], the transition probability is assumed to depend on the travel time over an edge. For the public transport and railway sub-networks, the travel time of each edge is calculated by multiplying the reference free-flow speed associated to each mode<fnref n="13" /> and the edge length<fnref n="14" />. Additionally, for public transport, cross-layers edges connecting the different lines and modes are associated to a typical waiting time of 5 minutes.

Finally, the transition probability $t(s_m, s_n)$ between the generic pair of nodes $s_m$ and $s_n$ is defined to be exponentially decreasing according to the travel-time weighted shortest path between the two nodes $s_m$ and $s_n$. Formally:

$$ t(s_m, s_n) = e^{-\beta \cdot tt_{s_m, s_n}}, \quad tt_{s_m, s_n} = \sum_{\forall (s_u, s_v) \in SP_{mn}} tt_{s_u, s_v} \qquad (3.33) $$

where $(s_u, s_v)$ is the generic edge on the travel-time weighted shortest path $SP_{mn}$ connecting the two nodes $s_m$ and $s_n$ in sub-graph $G_j$, computed via the Dijkstra algorithm. The length of the weighted shortest path $SP_{mn}$ corresponds to the sum of the travel time over each edge $(s_u, s_v)$ belonging to $SP_{mn}$, while $tt_{s_u, s_v}$ denotes the travel time between each pair of nodes $s_u$ and $s_v$. Finally, $\beta$ is a damping factor controlling the impact of the travel time on the transition probability.

<footnote n="13">¹³Speeds on the road network depend on the OSM type of route, ranging from 30 km/h to 90km/h. For the subway, the tramway and the bus the reference speed has been estimated respectively to 30 km/h, 15 km/h and 15 km/h.</footnote>

<footnote n="14">¹⁴The length of an edge has been computed as the geodesic distance between the two adjacent nodes of the edge.</footnote>
## Emission Probability

Similarly to [167], we modelled the emission probability as a Gaussian spatial noise centred around the hidden state $s_m$ and an empirically estimated standard deviation of the distance error between hidden states and observations:

$$r(s_m, o_k) = \frac{1}{\sqrt{2\pi}\alpha} e^{-0.5 \left( \frac{d_{s_m, o_k}}{\alpha} \right)^2} \quad (3.34)$$

where $d_{s_m, o_k}$ is the geodesic distance between the generic observation $o_k$ and the generic node $s_m$, while $\alpha$ is the standard deviation of a Gaussian random variable associated to the error distance between the reconstructed and the real position of the mobile. Both $\alpha$ from Eq. 3.34 and $\beta$ from Eq. 3.33 are hyper-parameters that need to be properly estimated via a sensitivity analysis.

## Results

To evaluate the performance of the proposed HMM map-matching solution, we used a refined version of the same small dataset introduced in Sec. 3.2.3 for the performance evaluation of TRANSIT. This dataset was also used to select the values of the $\alpha$ and $\beta$ parameters of Eq. 3.33 and Eq. 3.34 via a sensitivity analysis fully detailed in [160]. Concerning the GPS trajectories, we manually labelled the transport mode of all trajectories by associating one sub-graph $G_j$ of $G$ for each trajectory. The final ground-truth data contain 111 trajectories related to public transport and 72 to car, for a total of 183 trajectories. TRANSIT was applied to the NSD data available for the same set of volunteers to retrieve the mobile sessions with augmented trajectories. The latter were used as the input of our map-matching approach.

The results of the map-matching procedure were produced with and without prior knowledge on the transportation mode. In the case of map-matching performed without any prior knowledge on transportation mode, trajectories were matched to each of the considered sub-graphs of $G$. The final reconstructed output was chosen as the one with the highest likelihood from the Viterbi algorithm. In the other case, we considered that one could obtain a rough estimation of the transport mode, such as the average speed of the user from the observed NSD trajectory. With this information, and by associating a reference speed to each sub-graph, it becomes possible to identify the most likely candidate sub-graph for the reconstruction of the path for the given NSD trajectory.

To evaluate the map-matching performance, we used, among others, the distance metric defined in Eq. 3.31 to measure the geographical error, denoted as $G_e$ in the following, and the matching rate, denoted as MR and defined as the percentage of correctly map-matched edges by our approach. Formally:

$$MR = \frac{TP}{TP + FN + FP} \quad (3.35)$$

where: (i) the number of true positives TP is the number of edges in common between the ground-truth GPS and NSD map-matched trajectories; (ii) the number of false positives FP represents the number of edges in the NSD map-matched trajectory that do not belong to the corresponding ground-truth GPS map-matched trajectory; (iii) the number of false negatives FN represents the number of edges from the ground-truth GPS map-matched trajectory that do not belong to the NSD map-matched one.

The figures associated to these metrics indicate that the combination of TRANSIT with the map-matching step significantly increases even more the accuracy of the mobility information contained in the NSD data. Specifically, when no prior knowledge about the transport mode was assumed to be available, we obtained an overall matching rate (MR) of 63%, and the geographical error ($G_e$) was 110 meters on average for the complete set
<figure><img src="image_49.png" /><figcaption>**Figure 3.12:** Comparison between popular paths reconstructed by our approach and ground-truth popular paths for three case studies: $C_1$, $C_2$ and $C_3$.</figcaption></figure>

of mobile sessions. This level of spatial precision represents a first enhancement compared to the results in Sec. 3.2.3. The results improve significantly when the mode of transport of each trip is provided to select the sub-graph for the map-matching step, with a nearly halved $G_e$ of 60 meters and a 10% improvement in the MR, reaching a satisfactory value of 77%. Furthermore, the map-matching approach achieved superior results in public transport scenarios than on roads, with values as low as 30 meters and as high as 86% for the metrics $G_e$ and MR, respectively, when the prior knowledge of the transport mode is used. The reason for this outcome is the more complex topology of road networks compared to public
transport, making the map-matching problem more challenging in the former case.

We conducted additional validation by processing the entire dataset $D_L$ from Sec. 1.1.2. We used TRANSIT to extract all (possibly enhanced) trips for a limited subset of OD pairs. For each selected OD pair, we identified the most common paths from the map-matched TRANSIT-reconstructed trajectories, by isolating the transportation network edges with the highest number of occurrences from the total set of map-matched trajectories. We also retrieved typical multi-modal itineraries for the same OD pairs by exploiting a variety of route planners<fnref n="15" />, using the latter as ground truth. The results are graphically presented in Fig. 3.12, indicating an overall good match for most of the OD pairs considered.

In conclusion, the analysis demonstrated the potential for obtaining accurate route-level mobility information by transportation mode via network signalling data using simple, coarse transportation mode inference. Our approach can be applied at the scale of an entire country, and with the possibility of frequently recomputing mobility patterns over time or even in real-time. However, the map-matching approach still has many limitations, especially in terms of its limited ability to automatically perform transport mode detection. Several features, such as the probability that the Viterbi algorithm generates as output, the start time/duration of the trip, and more realistic transition probabilities based on time-dependent travel time information on the links of the transportation networks, could be utilised to improve the proposed approach and achieve even higher accuracy.

## 3.4 Towards Privacy-aware Mobility Inference

TRANSIT provided clear evidence of the necessity for special care in the handling of mobility data retrieved from massive NSD, as they might contain sensitive personal information. We take as a reference the GDPR European regulation, which defines as personal any piece of information pertaining to a particular individual. Based on this definition, mobility trajectories are very sensitive because of their frequent uniqueness. Individuals tend to have very different trajectories over the course of a day, which makes it difficult to extract trends or group them in order to provide more privacy-preserving mobility reports [137]. This problem applies anyway even to sparser and less pervasive sources of data such as CDR, LBSN, GPS and even surveys, which can often provide unique trajectories or indicators of individual mobility.

Let us consider for instance the case of a dynamic OD-matrix, a crucial indicator of mobility describing the flows between origins o and destinations d over a time window $[t, t + \Delta t)$. Typical values for $\Delta t$ range from 15 to 60 minutes, depending on the specific transport application and the data sources used to produce them [168]. A set of OD-matrices is then a regular, relational dataset with three attributes for each flow: an origin o, a destination d, and a time t usually denoting the time step of the beginning of the flow. Although they represent a dramatic simplification compared to trajectories, they still preserve high levels of uniqueness. Flows between specific pairs of zones can be in fact very small and isolated from others, and this problem is amplified when the adopted spatio-temporal resolution is finer<fnref n="16" />. In such scenario, OD-matrices are characterised by a high number of modalities<fnref n="17" />, i.e., the origin and the destination can be selected from sets of up to thousands of areas. These features make OD-matrices thus harder to anonymise than regular relational data.

To address this challenge, we proposed in [169, 170] a methodology to efficiently make dynamic OD-matrices meet GDPR's definition of anonymous data, i.e., the individual is

<footnote n="15">¹⁵https://www.google.fr/maps and https://www.viamichelin.fr/web/Itineraires</footnote>

<footnote n="16">¹⁶As discussed in Sec. 3.2.4, TRANSIT can be used to generate such indicators at very detailed spatio-temporal resolutions. We tested TRANSIT to reconstruct OD flows at spatial resolutions lower than 1 km² and with temporal binning in the order of minutes.</footnote>

<footnote n="17">¹⁷We call modality each pair of OD zones deriving from the adopted segmentation of the territory under analysis.</footnote>
not or no longer identifiable (as from Recital 26 of the GDPR). The widely used criterion of *k-anonymity* [171] was adopted to achieve this objective. An OD-matrix is said to be *k-anonymous* if no flow represents less than *k* individuals. By *k-anonymising* an OD-matrix, we can improve its safety against various categories of attacks [172], namely:

*   **Record linkage**, also called re-identification attack. If a flow in the OD-matrix has only one individual, an attacker can pinpoint a target individual in the data.

*   **Attribute linkage**, also called homogeneity attack. If all the individuals from a particular origin go to the same destination at a given time, then an attacker knowing that a target left the origin area can infer where they went, without needing to pinpoint the target first.

*   **Probabilistic attack.** If an attacker knows anything about a target, then accessing the OD-matrix can improve their knowledge about the target's whereabouts.

Each one of these three attacks is a relaxation of the previous one, the goal being less ambitious for the attacker but the success more likely. Note that probabilistic attacks may be successful even if the target user is not actually in the data, as long as we consider the OD-matrix to be representative of the population flows in the study area.

Technically, *k*-anonymity with *k* as low as 2 is secure against record linkage attacks, although this would understandably not be a satisfying anonymity condition. In particular, it would offer very limited protection against attribute linkage and probabilistic attacks. State-of-the-art usually aims at *k* between 3 and 5 in the case of hard-to-anonymise datasets, and up to 200 for simple datasets [173]<fnref n="18" />.

We used an approach to *k*-anonymisation based on **generalisation and suppression** [171], i.e., reducing the granularity (i.e., precision) of the data until the flows are big enough (i.e., larger than *k*) and suppressing the ones that are not. A toy example of 10-anonymisation generalisation and suppression is illustrated in Fig. 3.13: flows $A \to A$, $A \to B$ and $A \to C$ are aggregated together in order to reach a volume above 10, and similar aggregations are performed on the flows originating from *C* and *D*, and directed towards destinations *B*, *C* and *D*. Flow $B \to A$ is suppressed, which can be preferable when a very coarse aggregation becomes necessary to hide this flow among others.

Finding such a solution that minimises the loss in precision is known to be a NP-hard problem [176]. Historically, the first *k*-anonymisation algorithm was Datafly [171], which relies on a generalisation hierarchy describing how modalities should be merged together. An example of such a generalisation hierarchy is illustrated in Fig. 3.14: for a dataset giving the position of individuals in a study zone, the initial modalities are represented by the leaves of the tree. If we choose to generalise the whole dataset one level higher, then the possible modalities are the parents of the leaves. Datafly finds the best horizontal cut in the hierarchy, meaning all individuals (in our case, the OD flows) in the data are generalised to the same level, independently on the area they are localised into. This **uniform generalisation** approach has the advantage of being scalable to huge volumes of data, as well as data with numerous attributes. Within this category, the OIGH algorithm [177] represents the latest and best adapted solution to an OD-matrix. In a bid to find a finer-grained result, some approaches look for a generalisation at the individual level, which gives a solution akin to a clustering [173]. In the specific field of mobility data, this approach is better

<footnote n="18">¹⁸A complementary property for data anonymisation is *l*-diversity [174], which holds when each generalised value (e.g., the aggregate flow encompassing multiple modalities) covers at least *l* distinct modalities. Together, *k*-anonymity combined with *l*-diversity ensures protection against attribute linkage attacks, and significantly reduces the potential of probabilistic attacks. For even further protection, it is possible to implement *t*-closeness [175], which holds when the distribution of attributes in each group of *k* is no farther than a threshold *t* from the distribution of the overall population. However for OD-matrices, it can be argued that *k*-anonymity alone, given a sufficiently high *k*, is enough to make the individuals no longer identifiable: the areas of origin and destination can be large enough to avoid attribute linkage, and the low number of attributes renders probabilistic attacks rather vacuous.</footnote>
<figure><img src="image_50.png" /><figcaption>Initial OD-matrix</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>10-anonymous OD-matrix</figcaption></figure>

**Figure 3.13:** Example of generalisation and suppression of a simple OD matrix. Note that the flows here have been clustered without any kind of constraint. Some approaches seek particular solutions, notably ones where the generalisations form a partitioning of the domain.

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

**Figure 3.14:** Example of a spatial generalisation hierarchy. The root represents the whole study area, and the children of a node form a partitioning of the parents. An individual present in area A in the data can be generalised to be shown as present in area B, C or D depending on what is necessary in order to hide them in a group of *k* individuals.

represented by Glove [137], which generalises points in a dataset of trajectories. However, these approaches lack scalability for datasets characterised by a huge number of flows, as is the case for data derived from mobile phone CDR and NSD: for example, Liang and Samavi [173] report computing times in the order of hours for low values of *k*, and their approach can be expected to take days for *k* ≥ 10 given its exponential time dependency. Other solutions for general-purpose anonymisation rely on the notion of differential privacy, a robust privacy principle introduced by Dwork et al. [178]. Differential privacy does not apply to a dataset in itself but rather to a randomised algorithm taking the dataset as input and returning a series of query results. An algorithm is said to satisfy differential privacy if the probability of obtaining any output is not significantly different regardless of whether
an individual's data is included or not in the dataset.

To achieve rapid and accurate anonymization of OD-matrices, we developed an adaptive approach, formulated as optimisation problem, that utilises generalisation and suppression over a tree-based hierarchy. The problem aims at finding optimal partitionings (i.e., groupings) of the origin and destination areas that minimise the amount of introduced generalisation, with the constraint of satisfying the generalisation hierarchy and keeping low the number of suppressed records to a fixed maximum threshold. This last constraint is critical to guarantee the representativeness of the produced OD matrices for useful transportation applications. A comprehensive description of the optimisation problem can be found in [170] for interested readers.

We employed state-of-the-art algorithms for dependency-constraint knapsack problems [179, 180] to solve the optimisation problem and developed an adaptive generalisation technique that provides more refined results than uniform generalisation. Using this formulation, we proposed the Adaptive Tree Generalisation (ATG) approach, a lightweight algorithm that efficiently achieves *k*-anonymity. The ATG algorithm was developed in two versions, namely, ATG-Dual and ATG-Soft, which correspond to two variations of the problem we solved.

We evaluated our approaches against an extensive benchmark of anonymisation methods from the state of the art: general-purpose solutions from the field of anonymisation (i.e., OIGH [177] with uniform generalisation and Mondrian [181]), clustering from mobility data anonymisation (i.e., Glove [137] and our custom Glove-sk simplified implementation of Glove for large-scale datasets), and differential privacy (i.e., Laplace Noise). We also considered in the benchmark the results obtained by performing suppression alone of all flows that do not respect the *k*-anonymity constraint.

The approaches have been compared on a variety of datasets: New York City taxis (nyc), available as open-data<fnref n="19" />, and the Senegal (senegal) and Cote d'Ivoire (civ) datasets, available in the scope of the Data for Development (D4D) challenge [104, 103]. Several variations on the senegal dataset were also considered: senegal_crop with only the 599 eastern-most base stations, which account for most of the activity; senegal_big which increases the volumes in the data by summing together flows related to periods pertaining to different, but possibly overlapping, sets of individuals; and senegal_split which artificially increases the number of initial tiles by dividing the base stations in four and distributing the flows non-uniformly between them.

<table><caption>**Table 3.3:** Descriptive statistics of the datasets used for the experiments (#matrices: number of matrices in the whole dataset, where each matrix represents the flows over a time step; #tiles: number of initial tiles over which the matrices are set; density: average graph density of the matrices; avg. #flows: average number of flows among the matrices in the dataset; avg. vol.: average sum of the flows; %anon. flows: fraction of flows that are above *k* = 10; %anon. vol: fraction of individuals that are in flows above *k* = 10.)</caption><thead><tr><th>name</th><th>#matrices</th><th>#tiles</th><th>average #flows</th><th>density</th><th>average volume</th><th>%anon. flows</th><th>%anon. vol</th></tr></thead><tbody><tr><td>nyc</td><td>7242</td><td>263</td><td>3058</td><td>4.4%</td><td>15009</td><td>15.8%</td><td>59.0%</td></tr><tr><td>civ</td><td>632</td><td>1221</td><td>3523</td><td>0.2%</td><td>3117</td><td>0.3%</td><td>2.6%</td></tr><tr><td>senegal</td><td>6752</td><td>1666</td><td>18276</td><td>0.7%</td><td>41027</td><td>5.7%</td><td>36.8%</td></tr><tr><td>senegal_crop</td><td>6528</td><td>599</td><td>13710</td><td>3.8%</td><td>29614</td><td>5.8%</td><td>37.8%</td></tr><tr><td>senegal_split</td><td>360</td><td>6664</td><td>305194</td><td>0.7%</td><td>956742</td><td>6.5%</td><td>46.0%</td></tr><tr><td>senegal_big</td><td>360</td><td>1666</td><td>100322</td><td>3.6%</td><td>956742</td><td>12.3%</td><td>80.7%</td></tr></tbody></table>

Table 3.3 summarises the main characteristics of the datasets: we detail in particular the *density* of the matrix, computed as the number of non-null flows over the total possible number of distinct flows (i.e., the squared number of tiles), averaged over the whole period of observation, and the *percentage of flows that are 10-anonymous* in the original dataset.

<footnote n="19">¹⁹publicly available at: https://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page</footnote>
The relatively low amount of these 10-anonymous flows corresponds to a significant share of the total volumes of the data, as they are naturally associated with the OD pairs traversed by the most significant amount of traffic. The part of the anonymous volume in the original data is given under the column %anon\_vol.

For the choice of performance indicators, we assumed that the generalised (i.e., aggregate) areas provided as output by the anonymisation approach will most likely not be of interest for the final users of the data. Therefore, it appeared reasonable to assume that data users would rather be interested in drawing their own areas of interest and querying the matrix for an estimation of the flow between these areas. In the absence of additional information, the best estimation for such a custom flow is obtained by considering the volumes to be uniformly distributed in their generalised areas, and, as a result, the flow between the arbitrary areas can be derived proportionally to the overlaps with the intersected generalised areas. We call this estimation process **reconstruction**.

It is then relevant to actually reconstruct the anonymised OD matrix over the original segmentation and evaluate the difference with respect to the original OD matrix. The error metric E reported in the results is thus the absolute difference across all flows of the reconstructed anonymised matrix and the original one, i.e., formally:

$$E = \frac{1}{V} \sum_{o,d \in \text{leaves}(T)} |\tilde{v}_{o \to d} - v_{o \to d}| \quad (3.36)$$

where $\tilde{v}_{o \to d}$ denotes the reconstructed volume over the initial tiles (i.e., the leaves of the generalisation hierarchy), and $V = \sum_{o,d \in \text{leaves}(T)} v_{o \to d}$ normalises the error by the total volume of observed flows. Note that the sum is over the leaves of T, corresponding to the initial tiles over which the OD matrix is defined.

We also evaluate the approaches with respect to their total generalisation error $\bar{G}$, defined as follows:

$$\bar{G} = \frac{1}{V^+} \sum_{o \to d \in \mathcal{F}^+} (|o| + |d|)v_{o \to d} \quad (3.37)$$

where $\mathcal{F}^+$ represents the set of anonymous flows $o \to d$ such that $v_{o \to d} \ge k$, and $V^+ = \sum_{o \to d \in \mathcal{F}^+} v_{o \to d}$ the total volume of anonymised flows. When the origins and destinations are aggregated to roughly the same level as is normally the case, $\bar{G}$ represents roughly twice the number of tiles in the origin or destination of the average generalised flow. A value of $\bar{G} = 2$ then means that no generalisation was performed. We also consider the fraction of suppressed volumes S as an indicator of accuracy.

Finally, as a complementary metric to measure the distortion induced by the anonymisation, we evaluate the distribution distance D between each OD matrix and its anonymised versions. This metric is slightly different than the reconstruction loss E as it considers OD-matrices as normalised distributions, not penalising suppression. The distribution distance is given by:

$$D = \sum_{o,d \in \text{leaves}(T)} \left| \frac{\tilde{v}_{o \to d}}{V^+} - \frac{v_{o \to d}}{V} \right| \quad (3.38)$$

with the same notation as the previous definitions.

Table 3.4 summarises the average performance over all matrices in the small datasets nyc, civ, senegal and senegal\_crop, for which we were able to run Glove. We report the mean value of $\bar{G}$, E, D and S, and the total computing time to anonymise all the matrices. Firstly, Glove stands out for its lack of scalability, making it impractical in huge-volume cases.

Among the other solutions, the difference of a few minutes in the computing time is relevant but not decisive. They must instead be compared based on $\bar{G}$, for which ATG-Dual
data resolution. Our approaches facilitate the development of aggregate mobility indicators that enable fine-grained studies on human mobility while maintaining fundamental privacy requirements. However, it is important to note that *k*-anonymity alone can be insufficient to protect against attribute linkage and probabilistic attacks. To achieve complete privacy security, the method should also guarantee *l*-diversity and *t*-closeness. Nevertheless, the *k*-anonymity provided by our approach is a necessary initial step towards a more comprehensive solution and meets the current requirements of the French regulator CNIL, responsible for enforcing GDPR in France. Presently, mobility data are underutilised despite their enormous potential, largely discussed in this chapter, due to their size and personal nature, which makes their handling expensive and legally risky. Implementing an affordable, quick, and reliable anonymisation process should thus go in the direction of facilitating their widespread public use and promoting the extraction of useful insights to optimise and make more resilient transportation networks, as detailed in the next chapter.
Chapter 4

Assessing and Managing
Vulnerabilities in Transportation
Networks

Chapters 2 and 3 reported our findings on how massive data can provide accurate insights into the dynamic performance and actual demand of a large-scale, multi-modal transport system in both recurrent and non-recurrent situations. These insights contribute to *learn- ing from the past*, a crucial aspect of resilience as discussed in Sec. 1.2. Specifically, we demonstrated the efficacy of using passively-collected mobile phone data, particularly NSD, to reliably reconstruct both aggregate and fine-grained spatio-temporal estimations of travel demand and traffic variables in the form of OD matrices, regional speed dynamics, travel mode choices, popular paths, etc.

Although having a proper understanding of travel demand and traffic performance based on past observations is crucial for planning robust and efficient transportation systems, such knowledge alone is insufficient for effectively protecting these systems from disruptive events, for which data may be scarce or completely unavailable. Such events can lead to severe declines in the quality of mobility services and impact vulnerable segments of the transportation system and other urban infrastructures, thus potentially causing cascading failures and significant economic and social losses.

To reduce the consequences of disturbances and infrequent events, vulnerability assess-
ment enables the analysis of how engineered systems, particularly cyber-physical ones, re-
spond to specific kinds of disruptions. As detailed in Sec. 1.2.2, vulnerability assessment
corresponds to the anticipating (knowing what to expect) cornerstone of resilience. It aims
to identify risk situations before they occur and spot the system’s components that might be
unable to cope with such events, thus requiring special attention and enhancements that can
increase the system’s robustness and recovery ability. State-of-the-art studies on transport
vulnerability traditionally neglect the complexity of multi-modal transport topologies and
the impact that traffic conditions and travel demand dynamics can have on the system’s
response to a disruptive event [66, 38].

Additionally, extreme weather and human-related events represent a fundamental risk
category that increasingly challenges transport systems and whose impact on such networks
can be hard to model and anticipate. Nevertheless, with proper metrics and real-time data
supporting continuous monitoring of the entire transport network, its surroundings and its
most critical components, detecting such events rapidly and deploying appropriate strategies
to mitigate their consequences may become possible. These aspects are related to enhancing
the monitoring (knowing what to look for) cornerstone of transport resilience.

We also denoted a scarcity of operational decision-making tools to support resilient-driven
planning and management of multi-modal transport networks by considering the system's
exposure to both recurrent and unexpected disruptive events. These tools correspond to
the *knowing what to do* cornerstone of resilience and are essential for properly *responding* to
frequent disruptions and unlikely events, thereby enabling transport networks to maintain
adequate operations even during abnormal situations.

Based on these observations, this chapter summarises the main contribution of my re-
search on transport vulnerability assessment, monitoring and management, aimed at an-
swering the following questions: (i) What are the most appropriate metrics to identify
vulnerabilities of a transport network? (ii) What approaches can be leveraged to model
the occurrence of disruptions, evaluate the response of the transport network, and measure
the consequences of such disruptions in terms of performance? (iii) How can data-driven,
fine-grained, and time-varying knowledge of travel demand and traffic-related variables be
integrated into the identified approaches for more accurate and rapid detection of disrup-
tions and network vulnerabilities on a large scale? (iv) How can new on-demand mobility
services be optimally designed and integrated into an existing transport network to reduce
vulnerability to disruptive events?

The contribution to these research questions is detailed in this chapter as follows. Sec. 4.1 outlines the findings of a comparative study of approaches and metrics for vulnerability assessment. This study demonstrates the main limitations of traditional system-based and graph-based topological approaches, highlighting the need for novel combined strategies that can enact real-time monitoring of transport networks. Based on these findings, Sec. 4.2 proposes an original approach combining big data processing and complex network metrics for promptly detecting and possibly anticipating disruption occurrence, thus being helpful in implementing real-time monitoring of complex and dynamic networks. This study laid the ground for the research on real-time monitoring of large-scale transport networks, which will be detailed in Chapter 5. Finally, Sec. 4.3 discusses our solution based on multi-modality and stochastic optimisation to support the resilience-oriented design of transport networks.

The reported research activities have been performed in the context of the ANR JCJC PROMENADE project and have been at the core of Elise Henry's PhD thesis on resilience modelling of urban multi-modal transport networks, which I co-directed.

This chapter includes content from the following papers:

1. P. Gauthier, A. Furno, and N.-E. El Faouzi, “Road network resilience: how to identify critical links subject to day-to-day disruptions”, Transportation Research Record, vol. 2672, no. 1, pp. 54–65, 2018.

2. E. Henry, A. Furno, and N.-E. E. Faouzi, “Approach to quantify the impact of dis-
ruptions on traffic conditions using dynamic weighted resilience metrics of transport
networks”, Transportation Research Record, vol. 2675, no. 4, pp. 61–78, 2021.

3. E. Henry, L. Bonnetain, A. Furno, N.-E. El Faouzi, and E. Zimeo, “Spatio-temporal correlations of betweenness centrality and traffic metrics”, in 2019 6th International Conference on Models and Technologies for Intelligent Transportation Systems (MT-ITS). IEEE, 2019, pp. 1–10.

4. E. Henry, A. Furno, N.-E. El Faouzi, and D. Rey “Locating park-and-ride facilities for resilient on-demand urban mobility”, Transportation Research Part E: Logistics and Transportation Review, vol. 158, 102557, 2022

The research reported in this chapter is strongly related to the activities discussed in
Chapter 5. In particular, Chapter 5 describes the development of the PROMENADE ar-
chitectural and technological platform, which leverages cloud, IoT, micro-services and big
data solutions to provide the necessary support for enhancing resilience through real-time
monitoring and control in intelligent large-scale urban environments.
## 4.1 Vulnerability Analysis: a Comparative Study

In Sec. 1.2, we introduced the limitations of topological and system-based approaches to the vulnerability assessment of transport networks. These limitations were empirically confirmed in our work [55], in which we proposed a novel system-based methodology, termed the stress-test approach. This approach is based on dynamic mesoscopic simulation and aims to identify the most critical links within a road network. The methodology was conceived to model both day-to-day recurrent disruptions, such as accidents and mild weather-induced traffic degradation affecting road link capacity, and extreme events leading to sudden and significant variations in travel demand. The results obtained using this methodology were compared with those garnered using topology-based static metrics. This comparison highlighted the limitations of the latter from a functional perspective and those of the former from a performance perspective.

From a general point of view, stress testing consists of evaluating a system's performance under extreme conditions that go beyond normal usage, usually by relying on simulation. Stress testing helps identify a system's weaknesses or failure points and determine how it can be made more resilient. It has been widely used in finance [183, 184], as well as in the medical [185, 186] and the hydro-geology domains [187]. For instance, banks might conduct stress tests to evaluate the resilience of lending portfolios by modelling the impact of severe economic junctures on loan defaults. In the field of transportation, stress testing has been traditionally performed via approaches that remove edges or nodes from a network to model the occurrence of a (spatial) disruption, using performance indicators, e.g., based on total travel time increase, to evaluate the impact of the modelled disruption [46, 188, 189, 82].

The large majority of research works on stress testing of road networks focuses on major disruptive events [189, 190, 191], ignoring day-to-day disruptions that negatively and very frequently affect the performance of traffic networks. It is also worth noting the lack of studies combining system-based with topological approaches, as already mentioned in Sec. 1.2.2. Even though weighted graph-based models and metrics that include traffic flow data exist in the literature, their applicability to system-wide, dynamic traffic models is still neither widespread nor fully understood, being these metrics often inaccurate in highly dynamic environments as well as prohibitive to compute on large scale networks [61, 192].

The research described in this section was motivated by these considerations and represented our first step in the direction of joining solutions from the two categories of vulnerability assessment approaches (i.e., topological and system-based stress-testing), as further detailed in Sec. 4.2.

### 4.1.1 Graph-theory metrics

The topological metrics used in this study were based on Betweenness Centrality (BC), originally proposed in [193]. BC measures how central a link in a graph is by considering the number of the shortest paths that pass through the link, relative to the total number of shortest paths connecting any given pair of nodes. As highlighted in Sec. 1.2.2, BC represents one of the most widely used metrics in the literature to perform transport network vulnerability analysis. Edge BC is formally defined as follows for the generic link $l \in E$ of a graph $G(V, E)$:

$$BC(l) = \sum_{i \neq l \neq j} \frac{\sigma_{ij}(l)}{\sigma_{ij}}, \qquad (4.1)$$

where:

* $\sigma_{ij}(l)$ is the number of shortest paths from node $i$ to node $j$ that traverse link $l$;

* $\sigma_{ij}$ is the total number of shortest paths from node $i$ to node $j$.
An equivalent formula, replacing the generic link *l* with the generic node *n*, corresponds to the formal definition of *node* betweenness centrality, *BC*(*n*). In shortest path computation, links can be unweighted or weighted, as for example by the associated estimated travel cost (*e.g.*, travel time). In this study, we tested and compared both cases.

In order to model different aspects of a road traffic network and accelerate the computation of BC on larger networks, we also considered multiple variants of the BC, mostly consisting in calculating the shortest paths from travel demand entry to travel demand exit nodes only<fnref n="1" />. This definition introduces two advantages: computation time is reduced; the metric becomes more realistic from a demand-aware perspective since individuals tend to start and finish their trips over a subset of road intersections. The formula is the same as Eq. 4.1 with the following exceptions:

* *i* is selected from the entry-nodes subset, *i.e.*, a limited number of intersections where vehicles enter the road network;

* *j* is selected from the exit-nodes subset, *i.e.*, a limited number of intersections where vehicles leave the road network.

In conclusion, we considered four different formulations of the edge BC:

* Unweighted Betweenness Centrality (UBC)

* Travel-Time weighted Betweenness Centrality (TTBC)

* UBC on entry/exit nodes only (UBC entries-exits)

* TTBC from entry to exit nodes only (TTBC entries-exits)

It is worth remarking that we used static travel time information, *i.e.*, travel time in free flow conditions, to compute both *TTBC* and *TTBC* entries-exits.

### 4.1.2 Demand-sensitive Vulnerability Metric

Jenelius *et al.* introduced in [78] the demand-aware metric of link *Importance* (I) to characterize transportation vulnerability. This metric allows measuring the network performance loss induced by a disrupted link *l*, by using a generic metric of travel cost weighted by the traffic demand. Differently from [78] that uses a binary approach in modelling link disruptions (*i.e.*, a link is either fully operational or completely closed), we considered the possibility that the capacity of a link (*i.e.*, the maximum number of vehicles that can traverse the link) might be partially reduced of a percentage *δ* of the original capacity, due to accidents or lane closure. This technique allows considering different possible intensities of the disruptive event on a given link, which is especially helpful to model day-to-day disruptions. Such a metric appeared adequate in the scope of our methodological approach as it captures travel demand and the dynamic phenomenon of congestion (*i.e.*, travel costs increase when the network is congested).

In our study, we considered travel cost on a given origin-destination pair as the travel time to move from the origin to the destination divided by the associated travel distance. It is therefore measured in seconds/kilometres. The importance of a link *l* is thus defined in our work by the following equation:

$$I(l, \delta) = \frac{\sum_{i} \sum_{j \neq i} x_{ij} (c_{ij}^{\delta}(l) - c_{ij}^{0})}{\sum_{i} \sum_{j \neq i} x_{ij}}, \quad (4.2)$$

<footnote n="1">¹In a transportation network, not all nodes are responsible for generating and absorbing travel demand. Entry and exit points serve as specialised nodes that naturally associate with some level of travel demand. Their identification is generally tied to the definition of OD matrix zones and their urban topology, pinpointing specific nodes within each zone where trips can either commence or conclude.</footnote>
where:

* $x_{ij}$ is the demand from origin node i to destination node j (measured as number of vehicles);

* $c_{ij}^{\delta}(l)$ is the mean travel cost (measured in seconds/kilometres) from origin node i to destination node j when link l is disrupted at level $\delta$;

* $c_{ij}^{0}$ is the mean travel cost (measured in seconds/kilometres) from node i to node j in the base case, i.e., no disruption present in the network.

To model traffic dynamics and calculate travel costs on a vehicle basis, we used a dynamic mesoscopic simulator from our research team [194], based on the Lighthill-Whitham-Richards model [195, 196] and a dynamic traffic assignment procedure for distributing vehicles along the alternative paths in the network, according to the traffic conditions at the moment the vehicle is generated.

### 4.1.3 Stress Testing Methodology

Our methodology for road network stress testing aims to identify the most critical links in the road network that are considered vulnerable to different kinds of disruptions. The methodology is composed of the following steps:

#### 1. Simulating disruptive road events:

We proposed two strategies to perform this step. In the first, similarly to [82], we simulated day-to-day road disruptions as link capacity drops. The capacity-disruption level, denoted as $\delta$ in the following, is defined as the reduction in link capacity, expressed as a fraction of the original one. As anticipated, we considered the possibility of modelling partial capacity drops and analysing the evolution of the road network performance depending on the capacity-disruption level. For simplicity, we only allowed 5 possible capacity-disruption levels on each examined link, i.e., $\delta \in \{0\%, 20\%, 40\%, 60\%, 80\%\}$. Ultimately, the following equation allows computing the maximum capacity of the generic link $l$ in the presence of a capacity-disruption level $\delta$:

$$ \delta = 100 \cdot \left( 1 - \frac{q_{max}^{\delta}(l)}{q_{max}^{0}(l)} \right), \qquad (4.3) $$

where:

* $\delta$ is the capacity-disruption level applied to link $l$ (percentage) with $\delta \in \{0\%, 20\%, 40\%, 60\%, 80\%\}$;

* $q_{max}^{\delta}(l)$ is the capacity of link $l$ (measured in vehicles/hours) when it is disrupted at level $\delta$;

* $q_{max}^{0}(l)$ is the capacity of link $l$ in the base case (measured in vehicles/hours);

As a second strategy to simulate disruptive road events, we considered increases in the travel demand on specific entry/exit pairs of nodes of the network. By this approach, it was possible to simulate exceptional situations like city evacuations following extreme events (e.g., flooding, attacks, etc.) that typically put significant strain on the road infrastructure and result in total congestion of the network. This strategy consists in changing the OD matrix, i.e., increasing the traffic flow from given entries to specific exit points, and comparing the stress testing results with the modified demand level.

Based on the selected strategy for disruptive road events, we set the parameters of our mesoscopic simulator (e.g., link capacity, traffic demand) and we simulated the network in the specific setting. For both strategies, travel costs were collected for all vehicles in order to compute the performance metrics described in the following point.
## 2. Computing overall performance loss:

This step aims to quantify the consequences of the simulated disruptive road event on the operation of the whole road network. To this purpose, we used the notion of road network performance, measured via the Importance metric from Sec. 4.1.2. The overall Performance Loss (PL), when link $l$ is disrupted, was defined as an adaption of Eq. 4.2 to the situation where costs are known per vehicle:

$$PL(l, \delta) = \sum_{v=1}^{n} \frac{c_v^{\delta}(l) - c_v^{0}}{n}, \qquad (4.4)$$

where:

* $PL(l, \delta)$ is the overall performance loss when link $l$ is disrupted at level $\delta$ (seconds/kilometres);

* $c_v^{\delta}(l)$ is the travel cost of vehicle $v$ when link $l$ is disrupted at level $\delta$ (seconds/kilometres);

* $c_v^0$ is the travel cost of vehicle $v$ in the base case (seconds/kilometres);

* $n$ is the number of vehicles in the network.

In order to compare links and identify the most critical ones, a unique value of criticality must be associated with each link. To this purpose, we defined a global metric, namely the Stress Test Criticality (STC), to aggregate the multiple loss values derived from Eq. 4.4 for each possible capacity-disruption level (i.e., $\delta \in \{0\%, 20\%, 40\%, 60\%, 80\%\}$) applied to the given link $l$. The metric definition for the generic link $l$ follows:

$$STC(l) = \int_{\delta} PL(l, \delta), \qquad (4.5)$$

where:

* $STC(l)$ is the Stress Test Criticality when link $l$ is stress-tested (seconds/kilometres)

* $\delta$ is the capacity-disruption level of link $l$ (percentage)

* $PL(l, \delta)$ is the overall performance loss (seconds/kilometres) derived from one simulation.

## 3. Analysing the results:

At the end of our simulations, we obtained the values of the STC metric for all or a subset of the network links. We then computed the graph-based static indicators described in Sec. 4.1.1. This step involves analysing and comparing the rankings of network links produced by both the graph-based metrics and the overall performance loss indicator STC. Links with higher values of these metrics were given higher rankings.

### 4.1.4 Evaluation

In this section, we present the results of our stress-testing methodology and discuss the link ranking derived from the different selected metrics. We show that the ranking of critical links on the same network can significantly vary when different indicators are used, thus proving that simple modifications of one centrality indicator can have a relevant impact on the capacity of the metric to capture different facets of network vulnerability. Moreover, we discuss the advantages and drawbacks of each different approach in assessing road network vulnerability, towards the definition of new enhanced metrics and approaches for vulnerability identification in large-scale networks.

Our analysis has been performed in two different case studies: a first one related to a simple virtual network, used as a basic test-bed for our approach; and a second one
<figure><img src="image_53.png" /><figcaption>**Figure 4.1:** Simple virtual road network.</figcaption></figure>

corresponding to a real road network in France, which we used to confirm the validity of our results in a realistic scenario and to support the discussion on advantages and drawbacks of both simulation-based stress testing and topological metrics.

## Toy Network Case Study

The simple toy network is composed of 9 links (roads) and 8 nodes. Nodes correspond to 4 road intersections, 2 entry points and 2 exits. The duration of each simulation was fixed at 10 minutes. Fig. 4.1 depicts the network with numbered links and flow directions.

<table><thead><tr><th colspan="2">Destination</th></tr><tr><th>Origin</th><th>Exit 1</th><th>Exit 2</th></tr></thead><tbody><tr><th>Entry 1</th><td>375</td><td>375</td></tr><tr><th>Entry 2</th><td>250</td><td>1000</td></tr></tbody></table>

<table><caption>**Table 4.1:** Origin-destination matrices for the simple virtual network with different demand levels (A and B). Values are expressed as vehicles per hour.</caption><thead><tr><th colspan="2">Destination</th></tr><tr><th>Origin</th><th>Exit 1</th><th>Exit 2</th></tr></thead><tbody><tr><th>Entry 1</th><td>500</td><td>500</td></tr><tr><th>Entry 2</th><td>1200</td><td>300</td></tr></tbody></table>

We were able to measure the stress test criticality and all of the proposed topological metrics on all the links of the simple network, due to its limited size. To perform stress testing, we used both strategies for modelling disruptions described in Sec. 4.1.3, i.e., link capacity drop and traffic demand increase. Therefore, two different demand levels, reported as two different origin-destination matrices (A and B) in Tab. 4.1, were defined. The measures of stress-test criticality corresponding to the two different OD matrices are distinguished as STC A and STC B, respectively. It is worth noting that STC A and B are calculated with the same formula (Eq. 4.5) and the same simulation parameters, with the only exception of the travel demand. In our simulations, we applied sequentially five capacity-disruption levels (i.e., 0%, 20%, 40%, 60%, 80%) to each link. Then, we measured the network-wide performance loss (i.e., Eq. 4.4) consequent to the disruptions applied to the link.

The overall performance loss when using travel demand A is reported on the y-axis of Fig. 4.2(a), while the corresponding capacity-disruption levels (δ) are reported on the x-axis. Results for different links are depicted with different colours and markers in the figure, using linear interpolation. For the sake of readability, the figure only reports the five most critical links (i.e., those with the highest overall performance loss). Intuitively, the overall performance loss grows as the capacity-disruption level increases. In other words, a link capacity drop translates into an increase in network-wide travel cost.

According to the definition provided in Eq. 4.5, the stress test criticality of a given link corresponds to the area below the overall performance loss curve. Thus, in scenario A, *link 9* is the most critical to the operation of the whole road network, followed by *links 5, 1, 2* and *8*. By using link ranking from stress test criticality as a baseline, it is possible to compare the other link rankings as derived from the different topological metrics. Tab. 4.2
<figure><img src="image_54.png" /><figcaption>**Figure 4.2:** Stress testing on a simple test virtual network in two scenarios: performance loss for the top-5 most critical links with different capacity-disruption levels</figcaption></figure>

reports on such link rankings for both stress test criticality and the whole set of topological metrics.

As a preliminary consideration, it can be observed that the rankings of critical links on the same network may dramatically change depending on the selected metric, due to the different properties of the network captured by each of them. As an example, *link 5* is on top of all the topological rankings whereas, in terms of stress-test criticality, it is ranked second, below *link 9*. The top rank of *link 5* with all the topological metrics can be motivated considering the large number of shortest paths traversing this link: *e.g.*, paths (4, 5), (3, 4, 5), (1, 3, 4, 5), (6, 4, 5), (7, 5, 8) are all shortest paths.

The different rankings issued by STC A can be easily explained. If *link 5* is disrupted, the alternative paths (1, 2, 8) and (1, 2, 9) exist for all individuals departing from *entry 1*. Conversely, when *link 9* is disrupted, no alternative path exists for users willing to travel to *exit 2* from both *entry 1* and *entry 2*, thus resulting in significant congestion and consequent travel time increase for all individuals heading to *exit 2*. Additionally, traffic demand for *exit 2* is very high (see *OD matrix* in Tab. 4.1). That explains why *link 5* is more critical than *link 9* in terms of topology, but less critical than *link 9* when considering demand data, as made possible by our stress-testing methodology (based on dynamic simulations) and captured by the related criticality metric. This simple test clarifies how traditional demand-agnostic approaches may fail in properly ranking edge criticality. It also shows that alternative paths may become the shortest paths of the network when links are disrupted by adverse events, thus attracting traffic flow previously directed through the disrupted links. This is a critical factor that static graph-based approaches fail to capture. However, it is essential to clarify that this does not mean topological metrics are poor indicators of vulnerability. Instead, graph-based models of transportation networks should incorporate dynamic elements derived from real-time data (e.g., traffic sensors, cellular networks, IoT devices, etc.) and captured, e.g., via edge weights. Furthermore, metrics like betweenness centrality should be recalculated almost continuously to quickly identify significant disruptions in the network. These aspects are treated extensively in Sec. 4.2.

Another striking difference worth analysing regards *link 2*: it is considered one of the most critical ones according to the **UBC** metric from entries to exits, whereas it is the least critical one for the **TTBC** and the **TTBC** from entries to exits. The peculiarity of *link 2* is its length: it is the longest in the network, thus demanding more time to travel than the other links. Metrics like **UBC** and **UBC** from entries to exits are not weighted, i.e., all links are valued equally, and are consequently unable to grasp this important aspect. Differently,
<table><caption>**Table 4.2:** Simple network link rankings generated by the different metrics of criticality</caption><thead><tr><th>Link</th><th>STC A</th><th>STC B</th><th>UBC</th><th>UBC entries-exits</th><th>TTBC</th><th>TTBC entries-exits</th></tr></thead><tbody><tr><td>9</td><td>1<sup>st</sup></td><td>4<sup>th</sup></td><td>2<sup>nd</sup></td><td>1<sup>st</sup></td><td>4<sup>th</sup></td><td>3<sup>rd</sup></td></tr><tr><td>5</td><td>2<sup>nd</sup></td><td>2<sup>nd</sup></td><td>1<sup>st</sup></td><td>1<sup>st</sup></td><td>1<sup>st</sup></td><td>1<sup>st</sup></td></tr><tr><td>1</td><td>3<sup>rd</sup></td><td>3<sup>rd</sup></td><td>2<sup>nd</sup></td><td>1<sup>st</sup></td><td>4<sup>th</sup></td><td>3<sup>rd</sup></td></tr><tr><td>2</td><td>4<sup>th</sup></td><td>5<sup>th</sup></td><td>2<sup>nd</sup></td><td>1<sup>st</sup></td><td>7<sup>th</sup></td><td>5<sup>th</sup></td></tr><tr><td>8</td><td>5<sup>th</sup></td><td>1<sup>st</sup></td><td>2<sup>nd</sup></td><td>1<sup>st</sup></td><td>4<sup>th</sup></td><td>3<sup>rd</sup></td></tr><tr><td>6</td><td>5<sup>th</sup></td><td>6<sup>th</sup></td><td>4<sup>th</sup></td><td>2<sup>th</sup></td><td>5<sup>th</sup></td><td>4<sup>th</sup></td></tr><tr><td>4</td><td>5<sup>th</sup></td><td>7<sup>th</sup></td><td>2<sup>nd</sup></td><td>2<sup>th</sup></td><td>2<sup>nd</sup></td><td>2<sup>nd</sup></td></tr><tr><td>7</td><td>5<sup>th</sup></td><td>8<sup>th</sup></td><td>3<sup>th</sup></td><td>1<sup>st</sup></td><td>6<sup>th</sup></td><td>4<sup>th</sup></td></tr><tr><td>3</td><td>5<sup>th</sup></td><td>8<sup>th</sup></td><td>3<sup>th</sup></td><td>2<sup>th</sup></td><td>3<sup>rd</sup></td><td>3<sup>rd</sup></td></tr></tbody></table>

links with high travel times are not considered critical by the analysed weighted approaches, because they are not often part of the shortest paths. The same consideration applies to *link 7*, which is the second-longest link of the network.

Finally, it is worth noting that UBC values are often clustered. In particular, the UBC from entries to exits has many equal values and only two different link ranks. Traditional topological metrics appear to have minimal capability to discriminate link criticality at a fine level. In this case, stress test criticality does not differentiate all links either, but this is due to capacity-disruption levels. For *links 8, 6, 4, 7* and *3*, which are all ranked at the 5th place with the same value, demand levels are not high enough, relatively to the other links, to observe a significant performance loss compared to the base case. For example, when capacity disruption is applied to *link 8*, the overall network performance does not change at all (see Fig. 4.2(a)).

To further investigate this aspect, we stress-tested the same simple network with the second travel demand OD matrix, reported as *B* in Tab. 4.1. Results are shown in Fig. 4.2(b) and Tab. 4.2. Link ranking changes significantly when different demand levels are used. As an example, *Link 8* becomes the most critical link, whereas in the previous case, it involved no performance loss. This is due to the large increase in demand level associated with *exit 1*, which is directly connected to *link 8* (see Fig. 4.1).

**Takeaways:** Critical link ranking is highly variable as different approaches are used. Vulnerability assessment via topological metrics is limited in the sense that such metrics do not usually take into account traffic demand and network re-configurations following disruptive events. Conversely, the simulation-based stress testing approach is able to capture these aspects thus providing more realistic rankings via the proposed performance loss metric. Stress testing can also be used to compare different road networks and sub-networks, by analysing their response to similar stresses. *TTBC* produces better estimations of link criticality with respect to *UBC*, which treats all links equally.

## Real Road Network Case Study

To confirm the results of our previous analysis in a realistic scenario, we considered a real-world network, the DIRIF, situated in the South of Paris, France. It has 868 links and 827 nodes (657 intersections, 86 entries and 84 exits). Its roads are mostly highways. Each simulation duration was fixed at 15 minutes. As the network is much larger than the toy one, and as traffic flow can be very low on some links, we specified a longer simulation duration to ensure that enough vehicles could travel through the whole network and that we had a proper number of travel cost observations. The simulation was performed with realistic demand data from 9:00 AM to 9:15 AM, corresponding to the morning peak time, in order to increase the probability of observing performance loss as a consequence of the higher demand. The network is graphically presented in Fig. 4.3.

Given the large size of this network and the high computation time associated with
<figure><img src="image_55.png" /><figcaption>Figure 4.3: The DIRIF road network in Paris agglomeration</figcaption></figure>

<table><caption>Table 4.3: DIRIF network link values generated by the different metrics of criticality</caption><thead><tr><th></th><th>Link</th><th>STC</th><th>UBC</th><th>UBC entries-exits</th><th>TTBC</th><th>TTBC entries-exits</th></tr></thead><tbody><tr><td rowspan="3">Highest-flow links</td><td><b>95</b></td><td>117</td><td>2628</td><td>83</td><td>2628</td><td>83</td></tr><tr><td><b>93</b></td><td>102</td><td>660</td><td>83</td><td>660</td><td>83</td></tr><tr><td><b>94</b></td><td>101</td><td>1974</td><td>83</td><td>1974</td><td>83</td></tr><tr><td rowspan="4">Highest-BC links</td><td><b>802</b></td><td>42</td><td>192497</td><td>3029</td><td>192717</td><td>3039</td></tr><tr><td><b>803</b></td><td>27</td><td>192521</td><td>3029</td><td>192741</td><td>3039</td></tr><tr><td><b>607</b></td><td>27</td><td>192509</td><td>3029</td><td>192729</td><td>3039</td></tr><tr><td><b>608</b></td><td>44.6</td><td>192449</td><td>3029</td><td>192669</td><td>3039</td></tr><tr><td></td><td><b>397</b></td><td>15.8</td><td>83164</td><td>1139</td><td>83164</td><td>1139</td></tr><tr><td></td><td><b>672</b></td><td>14.5</td><td>10</td><td>1</td><td>10</td><td>1</td></tr></tbody></table>

each network simulation<fnref n="2" />, it was prohibitive to perform an exhaustive stress-test analysis as in the simple network case. Therefore, we decided to perform stress tests on a limited set of representative links: the three links with the highest flow of vehicles in the base simulation with no disruption, the three ones with the highest UBC and three randomly selected edges with UBC in three classes of values (i.e., high, medium and low). We discuss in the following our simulations related to the realistic available values of travel demand for the DIRIF network. Tab. 4.3 reports the actual values of the considered metrics for the analysed links<fnref n="3" />.

Consistently with our previous analysis on the simple network, Tab. 4.3 shows that rankings of critical links vary significantly when different metrics are used. As an example, *links 95, 93, 94* have a very high value of *STC*, whereas the topological metrics rate them much less critical than *links 802, 803, 607*. As pointed out in the previous section, the indications from the *STC* metric appear to be more reliable, since the metric seems to capture the higher criticality of *links 95, 93, 94* associated with higher demand.

On the small link subset considered in our analysis, taking into account travel times (*TTBC* and *TTBC* entries-exits) does not significantly change rankings with respect to

<footnote n="2">²Stress testing one link from such a large network with 5 capacity drops takes more than 1 hour on an Intel Xeon E3 CPU equipped with 8 GB of RAM running a Matlab implementation of the mesoscopic simulator.</footnote>

<footnote n="3">³Differently from Tab. 4.2, we do not report metric rankings but actual values of the metrics for each analysed link. This was due to the impossibility of performing all the required simulations for computing the *STC* metric for all network links.</footnote>
UBC, since lengths (and therefore the free-flow travel times) happen to be very similar for all considered links. Finally, it is worth noting that in the DIRIF network, UBC values (especially in the entries/exits variations) tend to significantly cluster themselves (i.e., many edges have very similar values of UBC), thus exhibiting a lower discriminant power than in the toy-network case study.

**Takeaways:** In a real-world scenario, stress testing proved to be a realistic and reliable approach to evaluate network resilience. Our evaluation confirms the importance of traffic demand and network dynamics for fine-grained ranking of the most vulnerable road-network links. Stress testing has however the important drawback of requiring very high execution times due to computationally intensive network simulations.

## 4.2 Graph Theory and Big Data Processing towards Online Transport Monitoring

As opposed to purely topological metrics, the system-based stress-testing approach described in Sec. 4.1 confirmed its ability to take into account travel demand and traffic dynamics to reliably identify links that are most critical to day-to-day disruptions and changes of travel demand. In our follow-up work [197], we also extended this approach for handling area-covering disruptive events, similarly to [83].

However, the stress-testing approach has a major drawback: it requires an extremely high number of computationally intensive traffic simulations. Therefore, it is recommended only for vulnerability analysis of small-sized networks, and as a preliminary step in application scenarios that allow for a sufficiently long response time, such as vulnerability-aware network design and maintenance planning (see Sec. 4.3). Another limitation of the stress-testing approach is its limited capability to model complex and hardly predictable disruptions that cannot be solely represented in terms of capacity reduction on single links or areas and/or an increase in travel demand. This is especially true for cascading failures, which combine, usually in sequence, multiple heterogeneous spatio-temporal disturbances causing major anomalies in the mobility system.

In scenarios that require rapid response times, such as recovery of critical infrastructure or emergency evacuations, swiftly detecting complex and unpredictable disturbances through real-time network monitoring is among the few viable strategies. However, effective transport monitoring hinges on access to up-to-date and detailed data on travellers' activity and transportation supply within the network. It also depends on selecting highly responsive metrics capable of quickly identifying or even anticipating the occurrence of disruptive events.

This section describes our research efforts aimed at exploring how topological indicators, specifically, indicators based on degree centrality (Sec. 4.2.1) and betweenness centrality (Sec. 4.2.2), could become suitable metrics for online network monitoring when calculated on a time-varying weighted graph. By this approach, the power of topological analysis combines with the availability of data on a system's dynamics, such as travel demand, speed and traffic flow, the latter being typically associated with system-based approaches. The outcomes of this research have confirmed that continuous data-driven monitoring of specific weighted graph metrics represents a key solution to increasing transport resilience. These studies represent the motivation at the basis of the platform and online tools for real-time monitoring of CIs, which will be described in Chapter 5.
4.2.1 Weighted Degree Centrality Metrics for Transport Network
Resilience

In [198], we proposed a dynamic modelling approach combining complex networks and big
data processing suitable for online continuous monitoring of road traffic networks.

Specifically, we focused on a specific metric, i.e., network heterogeneity, also known as graph irregularity, traditionally used as a dynamic indicator of resilience for gene regulatory, ecological and power supply systems. Heterogeneity provides global information on the diversity of network nodes connectivity [199]. When calculated over a time-varying, weighted graph, the heterogeneity metric can account for traffic dynamics, making it sensitive to perturbations. Furthermore, the repeated observation of specific patterns associated with this metric in the presence of disruptions suggests that it could be valuable for the development of effective, real-time monitoring tools. These tools could enable the automatic and rapid detection of such patterns and, consequently, of the associated disturbances.

The definition of heterogeneity is typically based on degree centrality. The degree central-
ity of a node describes the number of edges (or nodes) that are adjacent to the considered
node. On directed graphs, where each edge has a specific direction, it is possible to dis-
tinguish between in-degree and out-degree centrality by considering the number of edges
entering (respectively, exiting) the analysed node. The definition of this metric can also be
extended to the case of a weighted network (i.e., weighted degree centrality), considering the
sum of all the edge weights connected to the generic node i, i.e., C<sub>i</sub> = &Sigma;<sub>j&isin;N<sub>i</sub></sub> w<sub>i,j</sub>, where
C<sub>i</sub> is the non-normalised weighted degree centrality of node i, N<sub>i</sub> is the set of neighbours
of node i, and w<sub>i,j</sub> represents the weight of the edge connecting nodes i and j. As with
the traditional degree centrality, the weighted degree centrality of the generic node i can be
distinguished into in-degree and out-degree centrality, being equal to the sum of the edge
weights that join the node i to its predecessors (respectively, successors) j.

Different definitions of network heterogeneity have been proposed in complex network
theory. For instance, in [200], the authors use the entropy of nodes degrees, concluding that
such metric captures the disorder of the analysed network and represents an effective measure
of network resilience to random failures. Indeed, variations of this statistical measure can
describe losses of information in a network and thus indicate the impact of disruptions in
terms of reduced information transfer. Similarly, Jacob et al. [199] propose a definition
of heterogeneity based on the degree spectrum rather than degree values. Specifically, a
graph is considered heterogeneous if all the nodes have a different degree. Conversely, a star
network is considered strongly homogeneous because only one node has a different degree
from all others. The study can be extended to the weighted case, by computing a weighted
degree and analysing the in-degree and the out-degree distributions. Gao et al. [201] identify
heterogeneity as a strong predictor of network resilience, with more heterogeneous networks
being more resilient to targeted attacks but less resilient to random failures. Heterogeneity
is computed over a directed weighted graph, according to the density functions of the in-
and out-degree variables (k<sub>in</sub> and k<sub>out</sub>, respectively), i.e., formally:

$$
h = \frac{\sigma_{in} \sigma_{out}}{\langle k \rangle} \tag{4.6}
$$

where σ²ᵢₗ (respectively, σ²ₐₗ) is the variance of the density function of the weighted in-degree P(kᵢₗ) (respectively, out-degree P(kₐₗ)) and ⟨k⟩ is the average degree or network density, defined as follows:

$$
\langle k \rangle = \frac{1}{N} \sum_{i=1}^{N} k_i \tag{4.7}
$$

where N is the number of nodes in the network and kᵢ is the degree centrality of the node i, i.e., kᵢₗ(i) + kₐₗ(i).

To quantify the resilience of different types of multi-dimensional systems, the authors
propose to use a combination of three network metrics derived from the degree centrality,
<figure><img src="image_56.png" /><figcaption>(a) Node's unweighted (top) and weighted (bottom) degree centrality</figcaption></figure>

<caption>(b) Edge weight computation at time slot $t$ based on the (average) observed travel time</caption>

<caption>(c) The average speed profile (m/s) for working days and during the weekend.</caption>

<caption>**Figure 4.4:** The time-varying network modelling approach.</caption>

namely: heterogeneity $h$, density $\langle k \rangle$ and network symmetry $S$. Symmetry is defined as the correlation coefficient of in-degree and out-degree, i.e., formally:

$$S = \frac{\langle k_{in}k_{out} \rangle - \langle k_{in} \rangle \langle k_{out} \rangle}{\sigma_{in}\sigma_{out}} \qquad (4.8)$$

where $\langle k_{in} \rangle$ (respectively, $\langle k_{out} \rangle$) is the average in-degree (respectively, out-degree), $\langle k_{in}k_{out} \rangle$ the scalar product of both vectors (in-degree and out-degree).

We thus explored the possibility to apply the metrics proposed by Gao et al. in the context of road traffic networks, by studying their sensitivity to different disruptions and discovered the existence of relevant patterns, which could be exploited to quickly detect the appearance of disruptions in the network and anticipate recovery actions.

In the following, we provide a brief description of the methodology and the most relevant outcomes of its evaluation. The reader can refer to [198] for more details on the implementation, the dataset and the complete set of performed experiments.

## Network Modelling Approach

We modelled the road network of the city of Lyon, France, as a weighted, time-varying, directed graph, $G(V, E, W^t)$, with nodes $V$ corresponding to road intersections and edges $E$ corresponding to road segments. Edge weights $W^t$ were assumed to indicate dynamic traffic conditions (i.e., average speeds), computed every 30 minutes from available GPS floating car data. The free-flow travel time was used when no observation was available on a given edge at a given time slot in the reference period.

The first step of the proposed methodology consists in building a graph from the road network topology of the analysed city. Edge weights, at time $t$, are computed in the range [0, 1] by converting the average speed observed in that specific time slot into a traffic congestion indicator, with 1 indicating the best traffic conditions and 0 indicating the worst (Fig. 4.4(a)). In free flow conditions, all nodes are assumed to be connected at the best possible level, and all edges receive a weight equal to one. When the observed average travel time is higher compared to the free-flow travel time (fftt), the edge weight becomes closer to zero, according to the ratio of ffftt and the observed average travel time (tt) (i.e., the orange curve in Fig. 4.4(b)), with zero corresponding to the case of a completely congested edge, i.e., a road segment where vehicles are completely stuck. A discretization process was applied to the edge weights in order to model bounded rationality in travellers' route choice behaviour [202]. The idea is that a small travel time increase should produce a negligible impact on the edge weight (i.e., the step-wise function in Fig. 4.4(b)).

The analysis reported below explores the impact of selected disruptions on the road network of Lyon, through the lens of heterogeneity (Eq. 4.6), density (Eq. 4.7) and symmetry
(Eq. 4.8). The metrics are computed at each time step of 30 minutes, via the proposed model, under both normal and disturbed conditions. The analysis aimed to prove that they are indeed very sensitive to modifications of the travel demand, minor perturbations and extreme weather conditions. The results of the analysis revealed clear patterns in how the road traffic network reacts to day-to-day disruptions and more extreme disturbances, as evidenced by variations in the aforementioned metrics with respect to typical values observed in normal conditions. These findings are consistent with previous studies on static networks in other domains, such as Gao et al. [201].

Experimental Results

We performed an extensive evaluation of the previously described degree-centrality metrics in relation to four different disruptions that occurred in the period covered by the available GPS data<fnref n="4" />. The specific days were chosen based on the severity of the event or observed significant differences in speed or speed variance between disrupted and typical days.

On Monday, December 18th, 2017, a heavy snowfall in Lyon caused severe disruptions in the road network, which lasted until the afternoon when traffic returned to normal. On Tuesday, December 19th, 2017, the subway service was shut down for maintenance from 7:30 am to 4:15 pm. We studied the impact of this event on the road network because of the potential modal shift it generated. On Wednesday, April 4th, 2018, several protests took place in Lyon, affecting normal traffic circulation on the road network. Finally, on Saturday, June 2nd, 2018, the Fourvière tunnel, which is crossed by over 100,000 vehicles daily, was closed for renovations in the north-south direction for three days. Although some of these disruptions affected the entire network, others were localised in specific areas of the city. We report here results related to December 18th (heavy snowfall) and June 2nd (closed tunnel) and refer readers to our paper [198] for further details on the other days.

We observed significant variance of the aggregate speed variable computed on our dataset across the different days of the week (transparent margins in Fig. 4.4(c)). Therefore, we compared the speed profiles of disrupted days with the typical speed profile of the same day of the week (without major disruptions). For instance, if a disturbance occurred on a Monday, we compared the traffic conditions of that day with a typical (non-disrupted) Monday. A similar approach was used to study the time profile of the average degree (i.e., network density).

We begin by examining the typical speed profiles for vehicles on the same weekdays when the snowfall and closed tunnel events occurred. These profiles are represented by the blue curves in Fig. 4.5(a) for the snowfall scenario, and Fig. 4.5(b) for the closed tunnel scenario. Upon comparison, we find that the typical speed in the snowfall scenario (Fig. 4.5(a)) is lower than that in the closed tunnel scenario (Fig. 4.5(b)). This is because the first scenario occurs on a weekday (Monday), while the second one occurs on a weekend day (Saturday), with obviously different traffic conditions. The same trend appears in the typical degree profile displayed in blue (Figures 4.5(c) and 4.5(d)). The histograms in Figures 4.5(e) and 4.5(f) report the degree distributions, focusing on the morning peak time, 7:30 am. This is the time slot when we consistently noticed a significant gap between the typical average speed profiles and those related to disrupted events. For all studied cases (Figures 4.5(e) and 4.5(f)), we observe a higher number of nodes possessing a degree value between 6 and 9 in the typical distribution than in the disrupted one. However, for the ranges grouping small degree values below 3, the trend is reversed: there are more nodes for such degree values in the deteriorated situation than in the normal one. This result is not surprising and indicates that, in the presence of disruptions, the degree distributions shift towards lower values (worse traffic conditions). This is especially the case for the first disruption (Fig.

<footnote n="4">⁴The data, acquired by the FCD provider Media Mobile (now Be-Mobile), were available from October 2017 to September 2018 (see also Sec. 3.1, where the same data were used for speed estimation from UAPD data).</footnote>
<figure><img src="image_57.png" /><figcaption>**Figure 4.5:** Average speed profile, average degree centrality profile and degree distribution during two disruptive events: snow on December 18th, 2017 (left) and tunnel closure on June 2nd, 2018 (right).</figcaption></figure>

4.5(e)), which has the most significant impact on the average speed profile (Fig. 4.5(a)).

The comparison of the gap between the average degree centrality curves in the two disrupted scenarios confirms this result. Whereas the area between the curves is equal to 5.24 for the first disruption (Fig. 4.5(c)), it is only equal to 2.87 for the second one (Fig. 4.5(d)). In general, the average degree is always lower in disrupted conditions than in normal conditions on the same day of the week, showing a certain sensitivity of this global metric to events affecting the network. It is also worth highlighting that the gap is larger in the average degree centrality than in the average speed profiles, likely due to the weight discretization process, introduced to model bounded rationality in route choice.

We now focus on analysing the patterns associated with the combination of network heterogeneity, symmetry and density, in both normal and disrupted conditions. Similarly to [201], we visualise the values of these metrics in a three-dimensional plot (Fig. 4.6) in
<figure><img src="image_58.png" /><figcaption>**Figure 4.6:** Analysis of disruptions in the heterogeneity-symmetry-density space.</figcaption></figure>

order to characterise the network under different traffic conditions.

In the presence of a disruption (symbol × in Fig. 4.6), the graph density and its heterogeneity decrease, while the symmetry increases compared to reference situations (symbol ● in Fig. 4.6). The previous observations concerning the degree distribution in the presence of disruptions (i.e., values concentrating in lower ranges than in normal conditions due to reduced speeds in the presence of disruptions) easily explain the reduction in network density. The decrease in heterogeneity (Eq. 4.6) can be instead explained via the lower observed values of the standard deviation of the in-degree and out-degree distributions. Finally, the increase in the symmetry value reflects an increase in the correlation coefficient of in-degree and out-degree distributions. When a disruption occurs, both in-degree and out-degree values become generally lower, thus reducing the corresponding average degrees. However, this reduction is accompanied by a stronger decrease in the standard deviation of the degrees, which explains the increase in symmetry. Remarkably, the impact of disruptions (×) varies in intensity. As Fig. 4.6(a) highlights, the effect of snowfall (red symbols) is stronger than the other considered scenarios. Finally, the tunnel closure scenario (blue points and lines) is rather isolated in the plot space with respect to the other disruptions. The former occurred on a Saturday, while all others happened on weekdays when typical travel times are lower (Fig. 4.4(c)), explaining the difference in network states.

Our dynamic graph modelling approach also allows, differently from [201], to follow the evolution over time of the metrics in the three-dimensional plot. Specifically, we compared heterogeneity, symmetry and density every 30 minutes on a typical Monday and on the day (a Monday) of the heavy snow (Fig. 4.6(b)). The figure shows, with the same colour, but
with a different symbol, the triple (h, 〈k〉, S) at each given time slot, for both the typical and the disrupted situations. It can be noted that, in all cases, the time slots on the disrupted day exhibit lower heterogeneity, lower density and higher or similar symmetry values, confirming the already observed global trend. This aspect denotes that the disruption had negative effects on traffic throughout the whole day but with increased intensity in specific moments of the day. Specifically, on a regular Monday, the worst traffic conditions are observed at 7:00 am, i.e., during the morning peak hour (green ● points in Fig. 4.6(b)), while good traffic conditions can be observed, e.g., at 7:00 pm, in accordance with the speed profile of working days (Fig. 4.4(c)). When considering the day disrupted by the snowfall, the × points always exhibit higher symmetry, lower density and lower heterogeneity with respect to the regular day, independently of the time slot. The most significant changes between the typical state and the disrupted one can be observed during morning, i.e., at 7:00 am and 9:00 am (Fig. 4.6(b)), corresponding to the moments when the snow had actually the most disruptive impact on the transport network.

We also conducted a spatial analysis by separately computing the values of the three metrics in four districts of the city, at 7:30 am. Fig. 4.6(d) presents the per-district network density, heterogeneity, and symmetry on December 18th. It is relevant to underline that all areas exhibit the usual pattern in the presence of disruption: density and heterogeneity decrease, while the symmetry increases. For the typical Monday, the network states are close, even overlapping, for the 6th and the 8th districts under normal conditions. During the snowfall instead, in all cases, an increase of symmetry can be observed, but a significant reduction of heterogeneity is only observed in the 6th district, more urbanised and probably more affected by the snowfall. Similar considerations could be observed also for the other disruptive events, in both the local and global analyses.

In conclusion, the proposed approach utilising a dynamic, graph-based model and incorporating the concepts of heterogeneity, symmetry, and average degree revealed interesting insights into the behaviour of road traffic networks exposed to perturbations. The conducted experiments identified a specific pattern, where lower values of density and heterogeneity and higher values of symmetry appear to be associated with more significant deterioration of traffic conditions. This indicates a disrupted behaviour in contrast to typical baseline scenarios. Monitoring these topological indicators can thus provide valuable real-time information about a potential incoming perturbation, thus helping in preparing and recovering from it. Furthermore, these indicators can be used to compare different subsets of network nodes, such as different neighbourhoods within a city, to assess their robustness to disruptions.

### 4.2.2 Correlation of Betweenness Centrality and Traffic Metrics

In our efforts to explore useful metrics for monitoring transport networks, we investigated in [203] the properties of edge BC (Eq. 4.1) when computed on dynamically weighted graphs. We followed a data-driven modelling approach similar to the one described in Sec. 4.2.1 but with a different weighting procedure. Instead of using a congestion indicator, we assigned the average travel time to each edge, periodically computed every 30 minutes. We chose this approach to be able to determine temporal shortest paths in the definition of BC. In a first case study, we computed average travel times from real GPS FCD, and in the second study, we retrieved them through simulations using realistic data on demand and supply for the city of Lyon, France. In the following, we report the main outcomes related to the first case study. More details can be found in [203].

Our analysis confirmed the existence of a non-negligible correlation between BC and traffic flow, as already observed in related literature [204, 205]. In addition, we explored and clarified the nature of such correlation from a spatio-temporal perspective, as well as a function of structural and behavioural traffic properties. Similarly to the topological static approach discussed in Sec. 4.1.1, we considered different variations of edge BC, namely, (i) UBC, in which the graph is considered static and unweighted; (ii) Free-Flow weighted Be-
<figure><img src="image_59.png" /><figcaption>**Figure 4.7:** Correlation of per-edge average traffic flow and different BC metrics (UBC, FFBC, TTBC).</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>**Figure 4.8:** Per-edge temporal correlation between BC and flow (a), and zoom on a specific region (b,c).</figcaption></figure>

tweenness Centrality (FFBC), computed considering static free-flow travel times as weights on each edge; (iii) TTBC, which uses our dynamic data-driven approach with travel time information computed every 30 minutes and associated to the edges as weights.

In Fig. 4.7(a) and Fig. 4.7(b), we present the values of what we called *spatial correlation*. The latter was obtained by performing a linear least squares regression between the selected edge BC metrics (computed on a static graph) and the flows observed via loop detectors on the corresponding network edges. Each point corresponds to an edge with associated values of hourly flow (on the x-axis) and the edge's value of BC (on the y-axis). The results show a 0.34 value of $r^2$ for UBC and a 0.40 value of $r^2$ for FFBC, pointing to some mild capability of static BC metrics to explain flow distribution on the network.

We also found that TTBC is sensitive to traffic flow dynamics when computed on a time-varying weighted graph with average travel time weights: the correlation of edge BC and traffic flow varies, rather significantly, over time, decreasing during peak hours associated to higher demand, and, consequently, congestion. Fig. 4.7(c) shows the evolution over the different time slots of the $r^2$ coefficient of the linear regression between flow and three metrics: the static (UBC and FFBC) and the dynamic (TTBC) ones, the latter computed using average taxi travel times. TTBC appears to be the most sensitive to traffic dynamics among the three different metrics. Indeed, we observe an important drop of correlation during rush hours 07:30-09:00 for TTBC whereas the correlation for the FFBC metric remains practically constant, with mild variations only due to the dynamics of the traffic flow.

This phenomenon was partially explained by the existence of anti-correlation patterns
(i.e., reduction of BC in the presence of increased traffic flow, and vice-versa) for some edges of the analysed network, which was confirmed both in simulated and real-world experiments. To identify these patterns, we studied what we called *per-edge temporal correlation*. This correlation was computed by considering, for a given edge *l* of the graph, the vector corresponding to the **TTBC** values over the available time slots and the vector of the corresponding observed vehicle flows (measured as vehicles/hour). Hence, we calculated the linear correlation score (which can be either positive or negative) for each edge *l* of the network. The edge temporal correlation values are spatially visualised in Fig. 4.8(a). In the figure, highly anti-correlated edges can be spotted (orange/red links) from the positively-correlated ones (blue/green links), while yellow links correspond to weakly-correlated edges.

We analysed the dynamics related to **TTBC** and flows by focusing on a specific region of the graph, including two roads with a mirror behaviour in terms of temporal correlation, i.e., Quai Dr. Gailleton (QDG) and Quai Claude Bernard (QCB). Since the **FFBC** of QDG is much higher than the **FFBC** of QCB (see values reported in Fig. 4.8(b)), QDG should typically attract a larger flow than QCB, as confirmed by the observed average flow of QDG equals to 2,679 veh/hour, significantly higher than the one observed on QCB (967 vehicles/hour). A high value of **FFBC** combined with a high average flow suggests that the link can be often expected to be close to critical flow. We confirmed from the available speed data that congestion is actually observed on QDG at two different time slots, i.e., 08:00 and 08:30. Our data also indicates that when congestion happens on this link, an increase in travel time is also observed on the link. This is coupled to a decrease of **TTBC**, which thus exhibits an anti-correlation behaviour with respect to traffic flow. Conversely, the **QCB** link, which is characterised by lower demand, appears to be in free-flow conditions during the whole observation period. As travel time is rather high for QDG during congestion, QCB becomes then a viable alternative for drivers to avoid congestion on QDG. Consequently, we observe that **TTBC** increases as well as its flow. This explains why the **QCB** edge exhibits a highly positive temporal correlation. After the congestion phase (08:00 and 08:30), the flow globally decreases in the area. As QDG becomes more fluid, the corresponding travel time decreases to free-flow travel time and **TTBC** increases, i.e., thus maintaining an anti-correlated tendency. Contextually, **QCB** appears to lose its attractiveness compared to QDG as the flow decreases on this edge. **TTBC** decreases as well thus keeping the positively correlated trend. The analysis is confirmed by the evolution of the flow and the **TTBC**, which are known from the available data for both QDG and QCB as reported in Fig. 4.8(c).

To summarise, **FFBC** gives information on critical edges, i.e., road segments on which a high flow and possible congestion should be expected. Then, by studying the dynamic behaviour of **TTBC**, we highlighted different types of edges. On the one hand, anti-correlated edges seem to be critical by nature with higher **FFBC**. Even when congestion is observed on these links, the flow tends to remain high, whereas **TTBC** decreases. On the other hand, positively correlated edges appear to correspond to alternative routes in terms of travel time that become attractive when nearby roads are congested. In an intermediate situation, *neutral* edges (low correlation) are characterised by either low or medium/high demand, being able to dispatch such flow, without becoming congested.

Overall, the results of this study let us speculate that, if (i) people exploit the knowledge about current traffic conditions, (ii) alternative paths are available, and (iii) no specific constraint obliges travellers to follow specific routes, vehicles can re-distribute themselves in the network more efficiently and better exploit it. In this particular case, it can be assumed that edges with higher **BC** in a certain time slot represent critical (or very attractive) components of the network where traffic flow will most likely significantly increase in the near future. These properties could thus justify using this metric for the implementation of next-generation proactive, data-driven urban monitoring systems that can use alternative route recommendations to reduce traffic congestion and vulnerabilities.

These conclusions led us to investigate: (i) the possibility of computing in real-time, on weighted time-varying graphs, **BC** on large networks, which is a very challenging task due
to the high processing time of state-of-the-art BC computation algorithms; (ii) considering traffic control solutions that could take into account real-time information on traffic flow and BC to protect specific areas of the network, expected to attract massive traffic in the future. The research related to these two problems related to real-time monitoring and control is described in the next chapter. The following section focuses instead on managing disruptive situations occurring with a known (relatively low) frequency via multi-modal transport solutions for robust and sustainable network design.

## 4.3 Vulnerability-aware Multi-modal Network Design

The knowledge of a transport network's vulnerabilities, acquired, e.g., with the tools introduced in Sec. 4.1, can be harnessed in decision-making systems to enhance the network's ability to respond to specific disruptions. This corresponds to the *knowing what to do* cornerstone of resilience, as illustrated in Fig. 1.4(b). Within this research scope, we have explored *on-demand and multi-modal mobility services* as a means (i) to ensure greater robustness in the face of atypical events and (ii) to improve adaptability to large variations in travel demand that have already occurred.

Concerning the first point, i.e., increasing robustness to different kinds of disruptive scenarios, we have explored, in collaboration with the University of New South Wales, Sidney, Australia, the integration of on-demand mobility solutions into an existing large-scale urban transport network by taking into account the occurrence probability and expected traffic dynamics of a multitude of both regular and atypical scenarios [206]. In the same scope of research, we have also proposed using metrics for multi-layer networks to schedule the construction of new transit lines a transport operator might want to deploy to enrich its mobility offer [207]. Such methodologies typically fall within the broad field of transit network design. The novelty of these research activities [206, 207] lies in the consideration of multi-modality as a way for improving transport resilience. As an additional contribution, we considered the requirement of quick computation to allow for a rapid re-configuration of the transport supply, towards possible applications even in emergency situations.

Regarding the second point, i.e., enhancing adaptability in the presence of extreme variations of travel demand, we have explored, in collaboration with Ghent University, Belgium, optimisation and machine-learning-based solutions [208] for adapting transit supply to the decreased use of public transportation. The methodology was tested in the COVID-19 pandemic scenario, during which lock-downs and safety-related restrictions led to an over-capacitated and less efficient public transport system.

In the following, we summarise the methodological contribution and the main results concerning the optimal deployment of an on-demand mobility service to increase transport resilience [206]. The other mentioned contributions are left to be explored by the interested reader in papers [207, 208].

### 4.3.1 Park-and-Ride for Resilient On-Demand Urban Mobility

In [206], we focused on designing a specific form of on-demand ride-sharing service, namely, Park-and-Ride (PR), to improve the resilience of an urban transport network with respect to disruptive scenarios.

PR systems are transportation services enabling people to drive their private vehicles to designated *pick-up* locations, usually situated in the outskirts of a city. These pick-ups have parking lots where travellers can commute to public or private ride-sharing modes such as buses, subway, or on-demand shuttles. People can use these modes to move to *drop-off* locations within the city and reach close-by destinations. The idea is to allow people to reduce the costs of their trips and decrease the number of private vehicles circulating in the urban area. In opposition to traditional fixed-schedule/fixed-routes transit systems, ride-
sharing services offer flexibility in terms of routes, departure/arrival time, fleet size, etc., which can be helpful to reduce the total demand for car trips and increase city resilience [209].

Our study considered a **PR** facility location problem that captures users' entire trip chains – including access, transit and egress trips – in an urban mobility network. The **PR** system is assumed to comprise a shuttle service that complements the existing mobility network. User mobility outside of our **PR** system is captured through the *reserve mode*. In designing the **PR** facility, both typical mobility conditions and recurrent disruptive scenarios are taken into account. The objective is to maximise the overall ridership of the proposed **PR** system across all considered scenarios. Disruptions are modelled as scenarios with predetermined probabilities of occurrence, which must be known in advance. The aim was to specifically account for weather-related events (such as heavy rain and snowfall), pollution peaks, and other special events that could adversely affect regular transportation. The frequency of these disruptive events can be estimated in advance for a given urban area, in relation to typical traffic conditions. In our study, we quantified the attractiveness of the **PR** system by incorporating a simple model of user mode-choice behaviours into the formulation. This model is based on a multinomial logit choice model and seeks to maximise the expected ridership within the **PR** system at the expense of users' reserve travel option.

The proposed methodology provides: (i) the optimal location of pick-ups, i.e., the car parks where commuters can transfer from their vehicle to shared mobility solutions (e.g., shuttle buses); (ii) the optimal location of drop-offs, i.e., the transit stops where users exit the PR system and continue to their destination by other means, such as walking, micro-mobility (e.g., bike-sharing, e-scooters) solutions or the regular public transportation; and (iii) the expected flow share associated to the PR pick-up/drop-off pairs from the set of OD pairs, in each considered scenario. We assume a given budget constraint for constructing the car parks and transit stops.

Given the varying and hardly predictable nature of perturbations, as well as the large size and complexity of typical metropolitan mobility networks, guaranteeing a low computation time in solving the optimisation problem was considered a critical requirement to satisfy. By quickly solving the optimisation problem whenever a novel set/probability distribution of disturbances is identified, a different configuration of the park-and-ride system can be dynamically enacted (e.g., changing the location of the drop-offs) by transport authorities.

The proposed solution builds on and significantly extends the approach of Aros-Vera et al. [210], where the authors only focused on the pick-up and transit part of users' trip, thus ignoring the impact of drop-off location with respect to the final destination.

The main contributions of our study are: (i) the mathematical formulation of an integrated PR system, which captures the entire user trip chain from the origin to the destination via pick-up and drop-off nodes in a mobility network and accounting for mode choice; (ii) the incorporation of a stochastic programming approach to take into account recurrent disruptions in an urban mobility network; (iii) the development of a customised Lagrangian Relaxation Algorithm (LRA) able to provide competitive solutions, by out-performing commercial Mixed-Integer Linear Programming (MILP) solver CPLEX, based on a branch-and-cut search, for large-scale mobility networks in a restricted computational time; (iv) the implementation of the proposed PR system on a realistic instance representing the city of Lyon, France, which provides key insights for mobility service providers in urban areas. In the following, we briefly detail points (i) and (iv), leaving the interested reader to paper [206] for more details on the remaining contributions.

## Mathematical formulation of an integrated PR system

We modelled the **PR** problem as a stochastic problem which aims at considering a set $\Omega$ of different scenarios $\omega$, with a probability of occurrence $p^\omega$, representing both normal conditions and recurrent disruptions, to optimise our **PR** system in an existing transport network. The stochastic nature of the problem allows making resilient **PR** design, by tak-
ing into account the probability of occurrence of multiple possible scenarios of disruptions. Scenario-dependent variables and parameters are denoted with a superscript $^ω$ in the following. The transport network is represented by a directed graph with $|E|$ edges representing the roads and $|N|$ nodes corresponding to the road intersections. In each scenario, we determined the set of the OD pairs $(r, s)$ with a non-null total travel demand $d_{rs}^ω$. This set is called $W^ω = \{(r, s) : d_{rs}^ω > 0, r \in N, s \in N\}$.

For each OD pair, it is possible to identify the sets of the accessible pick-ups $i$ and drop-offs $j$. To be attractive, the generic PR pick-up $i$ has to be reached from the origin $r$ in a reasonable time $t_{ri}^ω$ lower than a fixed time $t^{access}$. Similarly, the destination $s$ must be reached from the drop-off $j$ in a given time $t_{js}^ω$ lower than the fixed egress time $t^{egress}$. The travel times $t_{ri}^ω$ (respectively, $t_{js}^ω$) correspond to the shortest travel time between the origin $r$ and the pick-up $i$ of our road network (respectively, the drop-off $j$ and the destination $s$). Ultimately, we introduced the set of the potential pick-ups $P_r^ω = \{i \in P : t_{ri}^ω \le t^{access}\}$ among all the pick-ups $P \subseteq N$ and the set of potential drop-offs $D_s^ω = \{j \in D : t_{js}^ω \le t^{egress}\}$ among all the drop-offs $D \subseteq N$. From these potential pick-ups and drop-offs, we extracted for each OD pair the set $\Sigma_{rs}^ω = \{(i, j) : i \in P_r^ω, j \in D_s^ω, (r, s) \in W^ω\}$ of the PR alternatives which correspond to the possible combinations of pick-ups and drop-offs at each OD pair in each scenario $ω$.

We defined a binary decision variable $y_i \in \{0, 1\}, \forall i \in P \cup D$, indicating the facility locations by being equal to one when a car park (pick-up) or a transit stop (drop-off) is open. Furthermore, we defined the decision variable $x_{rijs}^ω \in \mathbb{R}^+$, $\forall ω \in \Omega, \forall (r, s) \in W^ω, \forall (i, j) \in \Sigma_{rs}^ω \cup \{\mathbb{R}_{rs}\}$, as the flow which determines the part of the scenario-dependent demand $d_{rs}^ω$ associated to our PR system on the OD pair $(r, s)$ using the pick-up $i$ and the drop-off $j$. The decision variable $x_{rs}^{ω, \mathbb{R}_{rs}}$ defines instead the part of the users that choose the reserve mode on the OD pair $(r, s)$.

To determine the mode choice via a market share model based on a logit formulation, we used, for each itinerary, a generalised cost $g_{rijs}^ω$ composed of the travel cost between $r$ and $s$ using the proposed PR system, that we considered equal to the travel time, between the three parts of the path: (i) from the origin $r$ to the pick-up $i$, i.e., $c_{ri}^{access,ω}$; (ii) from the pick-up $i$ to the drop-off $j$, i.e., $c_{ij}^{route,ω}$; (iii) and from the drop-off $j$ to the destination $s$, i.e., $c_{js}^{walk}$. When using the reserve mode to realise the entire itinerary, symbolised as $\mathbb{R}_{rs}$, the travel cost $c_{rs}^{\mathbb{R}_{rs},ω}$, equal to the generalised cost $g_{rs}^{\mathbb{R}_{rs}}$, is modelled as the travel time to reach the destination $s$ from the origin $r$ by car. The probability that users choose the pair $i$ and $j$ of the PR system corresponds to the $x_{rijs}^ω$ decision variable, which can be expressed as follows:

$$x_{rijs}^{\omega} = \frac{y_i y_j e^{-\theta g_{rijs}^{\omega}}}{\sum_{(m,n) \in \Sigma_{rs}^{\omega} \cup \{\mathbb{R}_{rs}\}} y_m y_n e^{-\theta g_{rmns}^{\omega}}} \quad \forall (r,s) \in W, \forall (i,j) \in \Sigma_{rs}^{\omega} \cup \{\mathbb{R}_{rs}\}. \quad (4.9)$$

In other words, variable $x_{rijs}^ω$ represents the probability of using the reserve mode between the OD pair $(r, s)$, when $i$ and $j$ are equal to $\mathbb{R}_{rs}$, as well as the probability of using the proposed PR system by using a shuttle joining the pick up $i$ to the drop off $j$ for the OD pair $(r, s)$. Because both pick up and drop off have to be open to be considered as an option in the mode choice, the exponential of the utility for all the transport mode possibilities are multiplied by the pick up and the drop off location decision variables $y_i$ and $y_j$. The closure of a pick up, a drop off or both will remove the consideration of the associated transport modes in the logit computation. We denoted $θ$ the logit parameter corresponding to the user's sensitivity to the generalised cost. The higher $θ$, the more people are sensitive to the travel cost and will choose the cheaper transport mode. Finally, $-\theta g_{rs}^{\mathbb{R}_{rs},ω}$ (i.e., $g_{rs}^{\mathbb{R}_{rs},ω} = g_{rijs}^ω$ with $i$ and $j$ equal to $\mathbb{R}_{rs}$) corresponds to the utility of the reserve mode and $-\theta g_{rijs}^ω$ represents the utility of the PR system from the pick-up $i$ to the drop-off $j$ for the path going from the origin $r$ to the destination $s$. Although simple, our utility function considers the travel time preference in the mode choice. A more complex function could more accurately
represent this mode choice, for instance, by being trip purpose-dependent [211].

A budget B is allowed to the construction of the pick-ups and the drop-offs whose costs $c_i^{\text{loc}}$ mostly depend on the nature of the facility (Eq. 4.10b). Whereas car parks must be constructed at the pick-up locations for users to park their cars and embark on the shuttles of the PR system, transit stops are sufficient to leave the shuttle and reach the final destination by walking.

The resulting PR facility location formulation summarised in (4.10) is a MILP.

$$ \max \sum_{\omega \in \Omega} p^{\omega} \sum_{(r,s) \in W^{\omega}} d_{rs}^{\omega} \sum_{(i,j) \in \Sigma_{rs}^{\omega}} x_{rijs}^{\omega} \quad (4.10a) $$

$$ \text{s.t.} \quad \sum_{i \in P \cup D} c_i^{\text{loc}} y_i \leq B \quad (4.10b) $$

$$ \sum_{(i,j) \in \Sigma_{rs}^{\omega} \cup \{R_{rs}\}} x_{rijs}^{\omega} = 1 \quad \forall \omega \in \Omega, \forall (r, s) \in W^{\omega} \quad (4.10c) $$

$$ x_{rijs}^{\omega} \leq x_{rmns}^{\omega} \frac{e^{-\theta g_{rijs}^{\omega}}}{e^{-\theta g_{rmns}^{\omega}}} + (2 - y_m - y_n) \quad \forall \omega \in \Omega, \forall (r, s) \in W^{\omega}, \\ \forall (m, n), (i, j) \in \Sigma_{rs}^{\omega} \cup \{R_{rs}\} : \\ (i, j) \neq (m, n) \quad (4.10d) $$

$$ x_{rijs}^{\omega} \leq y_i \quad \forall \omega \in \Omega, \forall (r, s) \in W^{\omega}, \forall (i, j) \in \Sigma_{rs}^{\omega} \quad (4.10e) $$

$$ x_{rijs}^{\omega} \leq y_j \quad \forall \omega \in \Omega, \forall (r, s) \in W^{\omega}, \forall (i, j) \in \Sigma_{rs}^{\omega} \quad (4.10f) $$

$$ y_{R_{rs}} = 1 \quad \forall \omega \in \Omega, \forall (r, s) \in W^{\omega} \quad (4.10g) $$

$$ x_{rijs}^{\omega} \geq 0 \quad \forall \omega \in \Omega, \forall (r, s) \in W^{\omega}, \forall (i, j) \in \Sigma_{rs}^{\omega} \quad (4.10h) $$

$$ y_i \in \{0, 1\} \quad \forall i \in P \cup D \quad (4.10i) $$

By maximising the objective function (4.10a), the optimal solution will maximise the portion $x_{rijs}^{\omega}$ of the demand $d_{rs}^{\omega}$ using our PR system in each scenario $\omega$. The higher the probability $p^{\omega}$ related to the occurrence of a given scenario $\omega$, the higher the weight this scenario has in the objective function. The first constraint (4.10b) determines the ability to build the facilities by satisfying the global cost of the open parks. Constraint (4.10c) ensures that the whole demand is served for each origin-destination pair. Users must perform these trips through the PR system or using the reserve mode. The second constraint (4.10d) requires that the flow share of each mobility alternative obeys a linearised logit model. In [206], we provided the mathematical proofs, based on the ones proposed by [210], that lead to Eq. 4.10d from Eq. 4.9, which controls the proportion of the users for each mobility alternative. For each OD pair, this choice is governed by the utility of mobility alternatives through the PR system, which is a function of opened PR facilities and of the reserve mode. For a given PR path $(r, i, j, s)$, if the pick-up location $y_i$ or the drop-off location $y_j$ are closed, the constraint is inactive due to linking constraints (4.10e) and (4.10f), forcing the flow share to be lower than (or equal to) one when parks are open, and zero otherwise. For a given OD pair $(r, s)$, it is only possible to have a positive share of travellers using path $(r, i, j, s)$ if and only if $y_i = y_j = 1$. The decision variable representing the portion of users going from the origin $r$ to the destination $s$ through the pick-up $i$ and the drop-off $j$, $x_{rijs}^{\omega}$, must be positive (4.10h). Constraints (4.10g) and (4.10i) set the domain of binary variables $y_i$ and fix all such variables to one for the mobility alternative corresponding to the reserve mode.

The complexity of the problem highly increases with the graph size due to the number of variables. For instance, with a small network composed of 59 nodes used in our numerical experiments, the formulation corresponds to 267 000 constraints and 15 000 variables; the problem grows to 15 000 000 constraints and 270 000 variables for a larger one, composed of
<figure><img src="image_61.png" /><figcaption>**Figure 4.9:** Network segmentation (a); and the network used for the real instance $I^R$ (b).</figcaption></figure>

135 nodes. These figures made the problem hard or impossible to solve with commercial solvers and led us to develop a decomposition approach based on Lagrangian Relaxation (LR). The proposed LR, widely discussed in [206], provides upper bounds on the original problem and uses customised heuristic algorithms to iteratively generate lower bounds during the solution of the LR problem.

Experimental Results in a Realistic Case Study

The proposed formulation was solved in different case studies with both a commercial CPLEX solver and our LRA. All of the case studies were based on the road network of Lyon, France, composed of 10905 nodes and 19703 edges. The geography of the city, which is characterised by a peninsula surrounded by two rivers, makes the city centre particularly exposed to traffic disruptions and weather events, due to the limited number of alternative routes. The city centre is also the location of major cultural events, tourist attractions, and commercial establishments that attract large numbers of people. Finally, the pedestrianisation of Lyon's peninsula is currently an active topic of discussion in the city. For these reasons, we investigated the potential of an alternative, on-demand transport mode to improve accessibility and reduce the use of private vehicles. We focused on trips from outside the peninsula to the peninsula. All trip origins (r) and potential locations of pick-ups (i) were placed outside the peninsula, while all potential locations of drop-offs (j) and destinations (s) were considered within the first and second neighbourhoods that make up the peninsula.

Travel costs were determined using the FCD from October 2017 to September 2018, already introduced in Sections 4.2.1 and 3.1. The observed travel times were used to compute travel costs between each origin and potential pick-up locations (i.e., $c_{ri}^{access,\omega}$ costs) and between each pick-up and potential drop-off locations (i.e., $c_{ij}^{route,\omega}$ costs). The travel costs between the drop-offs and the destinations (i.e., $c_{js}^{walk}$) were computed by assuming a walking speed of 1.5m/s over the shortest paths of the road network. The travel costs associated with the reserve mode (i.e., $c_{rs}^{Rs,\omega}$) were calculated as the average observed travel times along the shortest path for each OD pair. For the sake of simplicity, we ignored any additional cost associated with the reserve transport mode (e.g., fuel, parking, congestion, pollution).
<table><caption>**Table 4.4:** Results on I<sup>R</sup> after 1 hour, with a heuristic time constraint of 20 minutes.</caption><thead><tr><th rowspan="2">Instance</th><th colspan="2">Parameters</th><th colspan="2">MILP</th><th colspan="2">LRA</th></tr><tr><th>t<sup>access</sup></th><th>t<sup>egress</sup></th><th>Budget (%)</th><th>Objective</th><th>Gap (%)</th><th>LB</th><th>Gap (%)</th><th>Gain (%)</th></tr></thead><tbody><tr><td rowspan="4">I<sup>R</sup></td><td rowspan="4">2</td><td rowspan="4">8</td><td>5</td><td>287.956</td><td>2.9</td><td>902.346</td><td>51.5</td><td>68.1</td></tr><tr><td>10</td><td>1 156.906</td><td>0.3</td><td>1 352.844</td><td>27.3</td><td>14.4</td></tr><tr><td>15</td><td>1 408.596</td><td>0.2</td><td>1 534.226</td><td>17.6</td><td>8.1</td></tr><tr><td>20</td><td>1 618.571</td><td>0.0</td><td>1 618.593</td><td>13.0</td><td>0.0</td></tr></tbody></table>

In the realistic case study, we considered four scenarios: i) a snowfall, ii) a public transport disruption with two subway lines stopped, iii) a typical weekday, and iv) a typical weekend day (see Sec. 4.2.1). Each scenario presents a different speed profile representative of the traffic dynamics during the corresponding event. Regardless of the studied scenario, we weighted the graph using the travel times recorded at 7:30 am, when the impacts of disruptions were the worst. The occurrence probability of each scenario was computed as follows. The snowfall was considered to have an impact similar to rainfall and we fixed its occurrence to 35%. The subway disruption scenario was considered exemplary of a public transport failure, representing a rather recurrent event. We fixed its probability to 30%. Finally, the probability for the normal weekday (respectively, weekend day) was set to 25% (respectively, 10%). Regarding the demand $d_{rs}^{\omega}$, we used realistic data reconstructed by other members of our laboratory [212] by relying on a methodology combining survey-based information, simulations and measures of flows collected via loop detectors installed in the city of Lyon.

In our experiments, we assumed that the construction cost $c_i^{\text{loc}}$ of pick-up nodes is ten times higher than that of drop-off nodes, as the former typically involves constructing a car park while the latter may only require a transit stop. The budget allocated to facility location construction, denoted as $B$, was defined as a percentage of the total cost for opening all facility locations. In all our analyses, we considered four different budgets corresponding to 20%, 40%, 60%, or 80% of the total construction cost. The trip generalised costs $g_{rijs}^{\omega}$ were determined by considering the travel costs $c_{rs}^{R,\omega}$, $c_{ri}^{\text{access},\omega}$, $c_{ij}^{\text{route},\omega}$, $c_{js}^{\text{walk}}$, as well as the demand $d_{rs}^{\omega}$ and the construction costs of pick-up and drop-off nodes $c_i^{\text{loc}}$. We conducted multiple sensitivity analyses (reported in [206]) to compare the performance of the proposed LRA approach with a direct MILP approach. Finally, only the morning peak hour was considered in this study to configure the PR system, by considering that the role of pick-ups and drop-offs would be reversed during the evening commute.

We assumed that the set of origins and destinations differs from the set of potential pick-up and drop-off nodes (as shown in Fig. 4.9(a)). To distribute the potential pick-up and drop-off nodes uniformly, we employed a segmentation of the city called Iris sectors (as shown in Fig. 4.9(b)). These sectors were developed by the French Institute of Statistics to divide the conurbation of Lyon into small geographical areas, each containing approximately 2,000 inhabitants<fnref n="5" />. We assumed that only one pick-up or drop-off node could be placed in each Iris sector. The location of the potential pick-up or drop-off was chosen as the road network node closest to the Iris area centroid. As for origin and destination nodes, we selected them to correspond to one-tenth of the road intersection nodes in Lyon's road network and uniformly distributed them throughout the city network. The final instance I<sup>R</sup> included 745 origins and 171 destinations with 155 potential pick-ups and 25 potential drop-offs. The access and egress times were respectively fixed at t<sup>access</sup> = 2 minutes and t<sup>egress</sup> = 8 minutes. Such a problem corresponds to about 9 062 000 variables and induced around 77 000 constraints, according to Formulation 4.10⁶.

<footnote n="5">⁵https://www.insee.fr/fr/metadonnees/definition/c1523</footnote>

<footnote n="6">⁶The proposed formulations and algorithms were implemented in Python on a machine with 16 Gb of RAM and a CPU of 4.20 GHz (Intel(R) Core(TM) i7-7700K). All mixed-integer and/or linear programs were solved using CPLEX's Python API.</footnote>
<table><caption>**Table 4.5:** Park-and-ride flow shares with different budgets, using $\theta = 0.1$, $t^{\text{access}} = 2$ minutes, and $t^{\text{egress}} = 8$ minutes.</caption><thead><tr><th rowspan="2">Description</th><th rowspan="2">Prob.</th><th colspan="4">Park-and-ride flow share</th></tr><tr><th>Budget 5%</th><th>Budget 10%</th><th>Budget 15%</th><th>Budget 20%</th></tr></thead><tbody><tr><td>Snowfall/Heavy Rainfall</td><td>0.35</td><td>42.67%</td><td>67.41%</td><td>76.84%</td><td>86.35%</td></tr><tr><td>Public transport disruption</td><td>0.30</td><td>46.95%</td><td>75.78%</td><td>82.70%</td><td>86.90%</td></tr><tr><td>Week day</td><td>0.25</td><td>54.82%</td><td>72.05%</td><td>86.25%</td><td>88.98%</td></tr><tr><td>Week-end</td><td>0.1</td><td>44.33%</td><td>69.56%</td><td>78.65%</td><td>82.79%</td></tr></tbody></table>

Tab. 4.4 summarises the results obtained with the realistic instance. The sensitivity parameter to the generalised cost $\theta$ was chosen equal to 0.1. Whatever the allocated budget (5%, 10%, 15% and 20%), the LRA always outperforms the MILP approach, by always finding a feasible solution with a higher value of the objective function. For the smallest budget, the gain between both solutions even exceeds 50%. Although the gap of the LRA is important due to the small reduction of the upper bound, the lower bound, representing the best found feasible solution, is very interesting to solve the complex problem of allocating PR, respecting a market share model, on a large-scale network.

Table 4.5 summarises the flow shares using the PR system, depending on the allocated budget for the facilities construction and the scenario. The higher the budget, the more people use the PR due to the increase of possible path alternatives induced by the increase in the number of open pick-ups and drop-offs (graphically reported with the related budget in Fig. 4.10). The larger set of alternatives reduces the travel cost of the PR mode by locating pick-ups and drop-offs closer to origins and destinations. It can be noticed how the PR can properly handle recurrent disruptions, by being highly attractive even in the presence of snowfall or public transport disruption. Although the flow share is not the same for all the scenarios, the order of magnitude remains the same. Additional experiments also proved that an increase in access and egress times generally leads to increased attractiveness of the PR option by serving a larger part of the demand for the same allocated budgets. Increasing the access and egress times makes PR alternatives accessible to more origins and destinations, thus augmenting the part of the demand served with the same amount of opened pick-ups and drop-offs. In these configurations with higher values of access and egress time, the MILP approach could not provide a solution contrary to our LRA, proving again the interest of our proposed methodology.

In conclusion, the results obtained on the realistic instance proved the ability of the proposed LRA to provide information about the implementation of a budget-constrained optimised PR system, adapted to recurrent disruptions in a large-scale network. The results also provided insights about the flow share distribution between the reserve mode and the PR for a given allocated budget. While flow shares of this magnitude may be unrealistically high, this can be explained by the fact that only a single concurrent mobility mode was considered – the reserve mode – which may not be sufficient to fully capture competition effects across travel options in an urban multi-modal transportation network. This study has focused on the design of PR systems with the aim to improve the resilience of urban mobility networks in the presence of recurrent disruptions. However, the proposed facility location approach can be generalised beyond PR systems to the design of flexible multi-modal and on-demand mobility networks. Indeed, the proposed framework could be used to model three-link trip chains in a generic context. We are currently working to incorporate fleet sizing and vehicle or line capacity in the proposed framework. As another limitation of the approach, it is worth noting that the addition of the PR system may influence transit flows and subsequently traffic congestion in the network. The incorporation of congestion effects in the proposed framework may be achieved by including a traffic and/or transit assignment module in our formulations. Such extensions of the proposed formulation will be explored in future works.
<figure><img src="image_62.png" /><figcaption>**Figure 4.10:** Solutions found for the real instance: green nodes are origins (potential pick-ups) and red nodes are the destinations (potential drop-offs).</figcaption></figure>
# Chapter 5

## Real-time and Large-scale
Solutions for Resilient Urban
Networks

As discussed in previous chapters, the resilience of transportation networks can significantly benefit from the availability of extensive, multi-source human mobility data. The latter typically arise from interactions between individuals and various communication devices, sensors, and applications. These data are stored and curated for later use in offline settings by traffic operators, telecommunication companies, and other service providers. As highlighted in Chapters 3 and 4, such data not only deepen our understanding of how transport networks function but also enable advanced strategies for vulnerability assessment and resilience-oriented design.

Modern cities have augmented this capability by generating real-time data streams through IoT devices and digital infrastructures. These real-time data can be utilised across various domains for numerous applications. These include, for instance, dynamic route optimisation and immediate intervention in the case of accidents; guidance for emergency vehicles and evacuation procedures during disasters; real-time repositioning of buses, taxis, and other transit options to adapt to dynamic travel demand fluctuations; real-time crowd management during special events; and urgent public health warnings, to name but a few. In this respect, Chapter 4 explored the potential of combining traffic data with complex network theory for resilience purposes. In particular, this synergy enables continuous surveillance of transportation networks, thus facilitating the quick identification of unexpected and difficult-to-predict disruptions [198, 213].

However, to unlock the full potential of both real-time data streams and extensive repositories of historical data, robust and heterogeneous technological support is imperative. From a software architectural point of view, it is necessary to provide flexible, scalable, and reliable solutions to accommodate the multifaceted needs that can arise in urban environments. Proper software infrastructure is required to manage a diverse array of envisioned services across various domains, which are often difficult to anticipate but usually rely on the same sources of data and share the same basic non-functional needs. Moreover, this support framework must provide general-purpose and complex data modelling capabilities to provide the necessary building blocks to implement the most variegated city applications from the available data. Importantly, it should be adept at managing the massive volumes of data generated within urban settings, especially in real-time scenarios. Finally, the infrastructure should be capable of adapting to fluctuating levels of usage load for these applications, ensuring that services remain accessible and responsive at all times.

Such support could thus take the form of a platform designed for the continuous and
seamless life-cycle of services, from design and development to deployment and operation. Over the past decade, several platforms designed for specific problems in various application domains have been introduced. Despite some advancements, we identified a clear gap<fnref n="1" /> for general-purpose city-oriented platforms that could be highly adaptable and easy to customise for supporting city resilience according to a complete data-driven perspective.

Starting from these needs, in this chapter, we resume the research on software architectural, algorithmic and technological solutions carried out within the ANR project PROMENADE, culminated with the development of an open-source prototype platform of the same name<fnref n="2" />. The platform was designed to facilitate the monitoring and management of large-scale urban infrastructures towards data-driven resilience and sustainability, focusing on a transportation-related prototypical implementation. In addition, we developed and integrated into the platform two solutions for real-time monitoring and management of large-scale traffic networks, i.e., rapid computation of betweenness centrality and real-time traffic zonal control. These solutions, based on the research results described in Chapter 4, represented two case-study applications of the proposed platform that confirmed its relevance for resilience enhancement and allowed for its performance evaluation.

The chapter is organised as follows. Sec. 5.1 details the requirements and modelling approach adopted for the architectural design of the PROMENADE platform. The platform global architecture and its main components are detailed in Sec. 5.2. The implementation of the platform core and its specialisation for large-scale road traffic networks is described in Sec. 5.3. This section also describes the realistic case studies that we considered for the evaluation of the platform prototype, by relying on road traffic datasets in the context of real-time road network monitoring. Sec. 5.4 describes our contribution to the exact computation of the betweenness centrality metric on very large-scale networks, a metric that we used for road monitoring purposes. Finally, Sec. 5.5 describes the proposed multi-agent cooperative traffic zonal control solution.

The results presented in this chapter have been the outcome of an intense research collaboration with University of Sannio, Italy, and have been at the core of the work of two PhD students (Cecile Daniel and Lorenzo Goglia) and two postdoctoral fellows (Mohammed Amine Merzoug and Pierre Lemaire) whom I have co-advised.

This chapter includes content from the following papers:

1. C. Colarusso, A. De Iasio, A. Furno, L. Goglia, M. A. Merzoug, and E. Zimeo, "PROMENADE: A big data platform for handling city complex networks with dynamic graphs", Future Generation Computer Systems, vol. 137, pp. 129-145, 2022.

2. C. Daniel, A. Furno, L. Goglia, and E. Zimeo, "Fast cluster-based computation of exact betweenness centrality in large graphs," Journal of Big Data, vol. 8, no. 1, pp. 1-39, 2021.

3. C. Daniel, P. Lemaire, A. Ladino, A. Furno, N.-E. E. Faouzi, and S. Hassas, "COMFORT: Cooperative multi-agent framework for large-scale routing-based traffic control", in 102nd Annual Meeting of the Transportation Research Board (TRB), 2023.

## 5.1 Modelling Approach and Platform Requirements

Cities are highly complex environments with diverse technical infrastructures co-existing and interacting continuously. For instance, city lights and water infrastructures are essential for the proper functioning of the road system. In addition to this structural interdependency, monitoring the vulnerabilities of a given city also means considering its inherent dynamism:

<footnote n="1">¹The interested reader can refer to the related work section from [214].</footnote>

<footnote n="2">²The source code of is available on GitHub, https://github.com/licit-lab/promenade_platform</footnote>