Romain Guigourès

Dominique Gay

Marc Boullé

Fabrice Clérot

Fabrice Rossi

Country-scale Exploratory Analysis of Call Detail Records through the Lens of Data Grid Models

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Telco operators' activities generate massive volume of data, mainly from three sources: networks, service platforms and customers data bases. Particularly, the use of mobile phones generates the so called Call Detail Records (CDRs), containing information about end-point antenna stations, date, time and duration of the calls (the content of the calls is excluded). While this data is initially stored for billing purpose, useful information and knowledge (related to human mobility [START_REF] Wang | Human mobility, social ties, and link prediction[END_REF][START_REF] Becker | Human mobility characterization from cellular network data[END_REF], social interactions [START_REF] Vieira | Characterizing dense urban areas from mobile phone-call data: Discovery and social dynamics[END_REF] and economic activities) might be derived from the large sets of CDRs collected by the operators.

domain, e.g., [START_REF] Frías-Martínez | An agent-based model of epidemic spread using human mobility and social network information[END_REF] suggest a model for epidemic spread. Mobile phones may also provide other types of data (e.g., the Nokia Mobile Data Challenge [START_REF] Laurila | From big smartphone data to worldwide research: The mobile data challenge[END_REF]), like applications events, WLAN connection data, etc. For instance, [START_REF] Jiang | A review of urban computing for mobile phone traces: current methods, challenges and opportunities[END_REF] pre-processed phone activities of one million users to obtain information about their approximative temporal location, then mined daily motifs from the spatio-temporal data to infer human activities. Finally, smart phones are or will be equipped with accelerometers and/or gyroscopes providing data about physical activities of users: [START_REF] Lockhart | The benefits of personalized smartphone-based activity recognition models[END_REF] suggest a complete system of activity recognition based on smartphone accelerometers with potential application to health monitoring. Research work related to data grid models: We are not coclustering data (objects × attributes) like pioneering work of Hartigan [START_REF] Hartigan | Direct clustering of a data matrix[END_REF]. Data grid models are related to the work of Dhillon et al. [START_REF] Dhillon | Information-theoretic co-clustering[END_REF] who have proposed an informationtheoretic coclustering approach for two discrete random variables Y 1 and Y 2 : the loss in Mutual

Information M I(Y 1 , Y 2 ) -M I(Y M 1 , Y M 2
) is minimized to obtain a locally-optimal grid with a user-defined number of clusters for each dimension. This is limited to two variables and requires to choose the number of clusters per variable. Going beyond 2D matrices, recent significant progress has been done in multi-way tensor analysis [START_REF] Sun | Beyond streams and graphs: dynamic tensor analysis[END_REF][START_REF] Kolda | Scalable tensor decompositions for multi-aspect data mining[END_REF]. Dealing with k-adic data, (also known as co-occurrence data, like contingency table), [START_REF] Peng | Temporal relation co-clustering on directional social network and author-topic evolution[END_REF] suggest a coclustering method for social network and temporal sequence (with pre-discretization of time). The Information Bottleneck (IB) method [START_REF] Tishby | The information bottleneck method[END_REF] stems from another informationtheoretic paradigm: given the joint probability P (X, Y ), IB aims at grouping X into clusters T in order to both compress X and keep as much information as possible about Y . IB also minimizes a difference in Mutual Information: M I(T, X) -βM I(T, Y ), where β is a positive Lagrange multiplier. Wang et al. [START_REF] Wang | Information bottleneck co-clustering[END_REF] build upon IB and suggest a coclustering method for two categorical variables. Extending IB for more than two categorical variables, Slonim et al. [START_REF] Slonim | Agglomerative multivariate information bottleneck[END_REF] have suggested the agglomerative multivariate IB that allows constructing several interacting systems of clusters simultaneously; the interactions among variables are specified using a Bayesian network structure. To the best of our knowledge, our summarization approach is the only one to combine the following advantages: it is parameter-free, scalable and can be applied to mixed-type attributes (categorical, numerical, thus multiple types of time dimensions without pre-processing). Therefore, the same generic method can be used to analyze network graph, temporal sequence and mobility data.

Impacts on economic strategy

Besides the high-level knowledge extracted from country-scale data and confirmed by local sociologists from the University of Bouaké in Ivory Coast, these studies have also a strong impact on future economic development strategy, mainly in two identified branches:

-Network planning strategy: In 2014, there are around 20M inhabitants in Ivory Coast and the mobile service penetration rate is 84% -with a still growing mobile phone market in a context of demographic growth. The analyses of the first two case studies and the resulting map projections (that can be seen as the network of calls available at various granularities, see Sections 5.1 and 5.2) are considered as an additional input for network planning and investment; for instance to help network designer in answering questions about how many and where the next antennas have to be set while preserving the quality of service at a reasonable cost.

Exploratory analysis through data grid models

Data grid models aim at estimating the joint distribution between K variables of mixed-types (categorical as well as numerical). The main principle is to simultaneously partition the values taken by the variables, into groups/clusters of categories for categorical variables and into intervals for numerical variables. The result is a multidimensional (K-d) data grid whose cells are defined by a part of each partitioned variable value set. Notice that in all rigor, we are working only with partitions of variable value sets. However, to simplify the discussion we will sometime use a slightly incorrect formulation by mentioning a "partition of a variable" and a "partitioned variable". In order to choose the "best" data grid model M * (given the data) from the model space M, we use a Bayesian Maximum A Posteriori (MAP) approach. We explore the model space while minimizing a Bayesian criterion, called cost. The cost criterion implements a trade-off between the accuracy and the robustness of the model and is defined as follows:

cost(M ) = -log(p(M | D) posterior ) ∝ -log(p(M ) prior × p(D | M ) likelihood )
Thus, the optimal grid M * is the most probable one (maximum a posteriori) given the data. Due to space limitation, the details about the cost criterion and the optimization algorithm (called khc) are available in appendix of [START_REF] Guigourès | Country-scale exploratory analysis of call detail records through the lens of data grid models[END_REF]. Hereafter, we focus on the tools for exploiting the grid and their applications on largescale CDRs data. The key features to keep in mind are: (i) khc is parameter-free, i.e., there is no need for setting the number of clusters/intervals per dimension; (ii) khc provides an effective locally-optimal solution to the data grid model construction efficiently, in sub-quadratic time complexity (O(N √ N log N ) where N is the number of data points). Figure 1 illustrates the input data and output results of khc on an examplary mobility data stemming from CDRs. Fig. 1: From 3D mobility data, stemming from CDRs data, to data grid.

Data grid exploitation and visualization

Because of the very large number observations in CDRs data, the optimal grid M * computed by khc can be made of hundreds of parts per dimension, i.e., millions of cells, which is difficult to exploit and interpret. To alleviate this issue, we suggest a grid simplification method together with several criteria that allow us to choose the granularity of the grid for further analysis, to rank values in clusters and to gain insights in the data through meaningful visualizations.

Dissimilarity index and grid structure simplification. We suggest a simplification method of the grid structure that iteratively merge clusters or adjacent intervals -choosing the merge generating the least degradation of the grid quality. To this end, we introduce a dissimilarity index between clusters or intervals which characterize the impact of the merge on the cost criterion.

Definition 1 (Dissimilarity index). Let c .1 and c .2 be two parts of a variable partition of a grid model M . Let M c.1∪c.2 be the grid after merging c .1 and c .2 . The dissimilarity ∆(c .1 , c .2 ) between the two parts c .1 and c .2 is defined as the difference of cost before and after the merge:

∆(c .1 , c .2 ) = cost(M c.1∪c.2 ) -cost(M ) (1)
When merging clusters that minimize ∆, we obtain the sub-optimal grid M (with a coarser grain, i.e. simplified) with minimal cost degradation, thus with minimal information loss w.r.t. the grid M before merging. Performing the best merges w.r.t. ∆ iteratively over the K variables without distinction, starting from M * until the null model M ∅ , K agglomerative hierarchies are built and the end-user can stop at the chosen granularity that is necessary for the analysis while controlling either the number of clusters/cells or the information ratio kept in the model. The information ratio of the grid M is defined as follows:

IR(M ) = cost(M ) -cost(M ∅ ) cost(M * ) -cost(M ∅ ) (2)
where M ∅ is the null model (the grid with a single cell).

Typicality for ranking categorical values in a cluster. When the grid is coarsen during the hierarchical agglomerative process, the number of clusters per categorical dimension decreases and the number of values per cluster increases.

It could be useful to focus on the most representative values among thousands of values of a cluster. In order to rank values in a cluster, we define the typicality of a value as follows.

Definition 2 (Typical values in a cluster). For a value v in a cluster c of the partition Y M of dimension Y given the grid model M , the typicality of v is defined as:

τ (v, c) = 1 1 -P Y M (c) × c j ∈Y M c j =c P Y M (c j )(cost(M |c \ v, c j ∪ v) -cost(M )) ( 3
)
where P Y M (c) is the probability of having a point with a value in cluster c, c\v is the cluster c from which we have removed value v, c j ∪v is the cluster c j to which we add value v and M |c \ v, c j ∪ v the grid model M after the aforementioned modifications.

Intuitively, the typicality evaluates the average impact in terms of cost on the grid model quality of removing a value v from its cluster c and reassigning it to another cluster c j = c. Thus, a value v is representative (say typical) of a cluster c if v is "close" to c and "different in average" from other clusters c j = c. Notice that this measure does not introduce any numerical encoding of the categories of the categorical variable under study.

Insightful visualizations with Mutual Information. It is common to visualize 2D coclustering results using 2D frequency matrix or heat map. For KD coclustering, it is useful to visualize the frequency matrix of two variables while selecting a part of interest for each of K -2 other variables. We also suggest an insightful measure for co-clusters to be visualized, namely, the Contribution to Mutual Information (CMI) -providing additional valuable visual information inaccessible with only frequency representation. Notice that such visualizations are also valid whatever the variable of interest.

M I(Y M 1 ; Y M 2 ) = J1 i1=1 J2 i2=1 M I i1i2 where M I i1i2 = p(c i1i2 ) log p(c i1i2 ) p(c i1. )p(c .i2 ) ( 4
)
where M I i1i2 represent the contribution of cell c i1i2 to the mutual information, p(c i1i2 ) is the observed joint probability of points in cell c i1i2 and p(c i1. )p(c .i2 ) is the expected probability in case of independence, i.e., the product of marginal probabilities.

Thus, if M I i1i2 > 0 then p(c i1i2 ) > p(c i1. )p(c .i2 ) and we observe an excess interaction between c i1. and c .i2 located in cell c i1i2 defined by parts i

1 of Y M 1 and i 2 of Y M 2 . Conversely, if M I i1i2 < 0, then p(c i1i2 ) < p(c i1. )p(c .i2
), and we observe a deficit of interactions in cell c i1i2 . Finally, if M I i1i2 = 0, then either p(c i1i2 ) = 0 in which case the contribution to MI and there is no interaction or p(c i1i2 ) = p(c i1. )p(c .i2 ) and the quantity of interactions in c i1i2 is that expected in case of independence between the partitioned variables.

The visualization of cells' CMI highlight valuable information that is local to the K -2 selected parts and bring complementary insights to exploit the summary provided by the grid.

Exploration results

Each application of khc (available at http://www.khiops.com) for the various case studies data is achieved within a day of computation on a commodity computer -which confirms the efficiency of the method.

Analysis of call network between antennas

The application of data grid models on the CDRs provides a segmentation with 1150 clusters, that corresponds to nearly one antenna per cluster. This is due to the large amount of data -471 millions CDRs. Indeed, the number of calls is so high for each antenna that the distribution of calls originating from (resp. terminating to) each antenna can be distinguished from each other. In order to Fig. 2: Evolution of the information kept in the data grid model w.r.t. the number of clusters using the ascending hierarchical post-processingfrom optimal data grid M * (100%) to the null model M ∅ (0%). obtain a more interpretable segmentation, we apply the post-processing introduced in the Section 4.1. Figure 2 shows the information ratio (see definition 1) versus the number of clusters for all intermediate models obtained during the ascending hierarchical post-processing. Interestingly, the resulting Pareto curve shows that very informative models are obtained with few clusters. In our study, we decrease the number of clusters until keeping 60% of the model informativity -corresponding to 20 clusters, an admissible number for the interpretation. Throughout the simplification process, both partitions of source and target antennas stay identical. Thus we consider only the partition of source antennas for the rest of the study. Those clusters are projected on a map of Ivory Coast in Figure 3. Antennas are identified using dots, which color matches with the cluster they belong to.

The first observation is the strong correlation between the clusters and the geography of the country. Indeed, antennas from a same cluster are close to each other. The size of the clusters is almost the same in terms of area and match with the administrative zones of the country, except for Abidjan, the economic capital, which is split into four clusters. This is due to the high concentration of antennas in the city (32% of the ivorian antennas) and the dense phone traffic (34% of the calls).

We use the typicality (see defini-Fig. 3: Twenty clusters displayed on Ivory Coast map. There is one color per cluster.

tion 2) to rank the antennas of each cluster. The place, where the antenna with the highest typicality is located, is used to label the cluster. On the map in Figure 3, the size of the dots are proportional to the antenna typicality. Most typical antennas are located in the main cities of Ivory Coast. This phenomenon has already been observed in [START_REF] Blondel | Regions and borders of mobile telephony in Belgium and in the Brussels metropolitan zone[END_REF] and [START_REF] Guigourès | Segmentation of towns using call detail records[END_REF]: the discovered clusters match with the area of influence of the main cities of a country. We observe few exceptions: the cluster of the city of Sassandra contains the antennas of the city of Divo, while Divo is almost 4 times bigger than Sassandra (population wise) and is the sixth Ivorian city. Antennas in Divo are 40% less typical than the ones in Sassandra, meaning that allocating them to another cluster would be less costly for the criterion. Actually, calls emitted from Divo are significant in direction to other regions of Ivory Coast whereas calls from Sassandra are more internal to its region. In more formal terms, the calls distributions of the antennas in Divo are closer to the marginal distribution than to its cluster's distribution. This observation is not really surprising because Divo has experienced a recent growth of its population, due to migrations within the country [START_REF] Gnabéli | La production d'une identité autochtone en Côte d'Ivoire[END_REF]. Divo is also located in an area specialized in the intensive farming, that attracts seasonal workers from other parts of Ivory Coast. Now, focusing on the segmentation of Abidjan: the city is divided into four parts with a strong socioeconomic correlation. The first cluster -in red in Figure 3 -covers central Abidjan, including the Central Business District (le Plateau), the transport hub (Adjamé) and the embassies and upper class area (Cocody). The second cluster -in light green -is located in the South of the city. The covered neighborhoods are mainly residential areas and ports. Note that this cluster and the previous one are separated by a strip of sea, except for its North part that is included in the previous cluster. This very localized neighborhood matches with the party area of Abidjan. Finally, the last two clusters group antennas located in two areas with a similar profile: these are lower class neighborhoods. These clusters are separated not only because they are located in different parts of the city but especially because their call distribution differs: Abobo in dark blue and Yopougon in grey in the Figure 3.

Traffic between clusters. Now, we Fig. 4: Analysis of Excess of calls between clusters of antennas analyze the distribution of calls between clusters of antennas using the contribution to the mutual information. We suggest to visualize the lacks and excesses of calls between the clusters, compared to the expected traffic in case of independence. Whatever the granularity level of the clustering, we observe a strong excess of calls from the clusters to themselves and weaker excesses and lacks between clusters. Studying the traffic within the clusters has a limited interest. We only focus on the inter-clusters traffic. To visualize the traffic between clusters, we use a finer clustering than previously. Here, we have 355 clusters for 95% informativity (see Figure 2). Figure 4 depicts the excesses of traffic between clusters -highlighted with red segments. The end points of the segments are drawn at the positions of the most representative antennas of the associated clusters (i.e with the highest typicalities). The opacity of a segment is proportional to the value of the contribution to mutual information and its width is proportional to the number of calls between clusters. The biggest cities -like Bouaké, San Pedro and Man -are clearly marked on the map: they are regional capitals, a fact that is confirmed and highlighted by the call traffic visualization. The case of Bouaké is particularly interesting: although it is not the country capital, its national influence seems bigger than the one of Yamoussoukro, the actual capital. Yamoussoukro is twice smaller than Bouaké (population wise) and is a quite recent city where there is no major economical activity, contrary to Bouaké. This fact can explain our observation. We also observe that excess of traffic between major cities is a rare phenomenon. Cities are more like phone hubs, except in the West of the country around Soubré. This area is not a densely populated area but corresponds to a region with important migration flows. Finally, in Abidjan, we observe important excesses of traffic within neighborhoods, but not between neighborhoods.

Temporal Analysis of the Calls Distribution

From previous section, we learnt that the correlation between source and destination antennas is very high. The evolution of the calls distribution over time might be the same for both sets of antennas. Therefore, to track the evolution of traffic over time, we only study the evolution of the originating calls: one call is described by the emitting antenna and a day count (stemming from the date).

Again, the clustering of antennas resulting from the optimal data grid is also too fine for an easy interpretation (1051 clusters of antennas and 140 intervals for the day count). We coarse the grain of the grid with our hierarchical post-processing so that the informativity of the model is 80%, with ten clusters of antenna and twenty time segments. Since, missing values are abundant in this data, i.e., some antennas emitted no call during some time periods, consequently, we obtain time segments that are strongly correlated with missing data. For the same reason, antennas are grouped together because they experienced an absence of calls during one or several similar periods. In the Figure 5, the colored antennas belong to clusters having experienced simultaneous absences of calls. We observe that the green, orange, light blue and purple clusters are located in localized area. The missing data appear during short periods for these clusters. This grouping might be due to localized technical issues on the network. The antennas of the yellow cluster are spread over the country. These antennas are grouped because they have been activated at the same date. This use case provides a better understanding the dysfunctions in the network over the year.

Analysis of output communications w.r.t. week day and hour

Our objective is to build simultaneously a partition of the antennas, a partition of the week days and a discretization of the hour, i.e., a triclustering. For the same reasons as previously, we only keep the emitting antennas. At the finest level, we obtain a triclustering with 806 clusters of emitting antennas, 7 clusters of days and 22 time segments. Again, these results must be simplified to ease the interpretation. However, we fix the numbers of clusters of days and time segments, since they are acceptable for the analysis and we only reduce the number of clusters of antennas. With four clusters of antennas, we keep 51% of the informativity of the model.

Antennas are displayed on the map Fig. 6: Clusters on the map of Ivory Coast. Dots are antennas. There is one color per cluster.

of Figure 6. We also build a calendar (see Figure 7) for each cluster with days in columns and time segments in lines. The color of the cells indicates the excesses (red) or the lacks (blue) of traffic emitted from the corresponding cluster. The lacks and excesses are measured using the contribution to the mutual information (see definition 3) between the cluster and the cross product of the cluster of weekday and the time segment:

M I(X M 1 ; X M 2 × X M 3 )
, with X M 1 the partition of the antennas, X M 2 the partitions of the weekdays and X M 3 the discretization of the time. Now we focus on the analysis of each cluster of antennas that we can easily label manually: Abidjan -Le Plateau (yellow) This cluster covers exactly the Central Business District of Abidjan. In the calendar of Figure 7, we observe an excess of calls from the Monday to the Friday, between 8-9am and 4-5pm. The rest of the time, there is a low lack of traffic emitted from this area. In other words, during the office hours, the phone traffic is higher than expected and lower the rest of the time. This is expected and representative of this type of area: a non-residential business district.

Economic zones (red) The antennas of this cluster are located either in the commercial areas of the cities or in areas with a strong economic activity, like plantations or mines. In Abidjan, these antennas are located in industrial zones (South and North-West), the shopping districts (North of the business district) and the universities and embassies neighborhood (East). The traffic in these areas is mainly in excess from the Monday to the Saturday between 9 am and 5 pm. The correlation is very strong between the working hours and the calls traffic on these areas.

Urban residential areas (blue)

The antennas belonging to this cluster are mainly located in the cities like Abidjan, Bouaké and Yamoussoukro. If we focus on Abidjan, we observe that the cluster covers the residential neighborhood located in the West and in the North-East of the city. At a finer level of partition of the antennas, this cluster would be split according to the socioeconomic class of the neighborhood: the upper class neighborhood in the East of the city is separated from the lower class neighborhoods, located in the North and the West. The calendar shows lacks of calls during the office hours and excesses the weekend, the night and the early morning during the week. This is correlated with the presence of people in residential areas. Note that the excesses of calls start around 8 pm, while it stops around 5 pm in the Central Business district or in economic areas. This time lag is due to the cheaper price of calls after 8 pm.

The countryside (green) The antennas of this cluster are spread over the country, except in Abidjan and other cities in general. The calendar for this cluster is quite similar to the one of the urban residential areas, except that the excess periods are limited to the early evening and the whole Sunday.

User mobility analysis w.r.t. week day and hour

Among the 50000 anonymized users, we focus on mobile users characterized by a frequent use of a large set of distinct antennas: after filtering, 6894 users are under study. For these 4-d data (user, antenna, week day and hour), khc operates a tetra-clustering: as a result, users with the same mobility profile are grouped together, i.e., users who have connected to similar groups of antennas, on similar days of the weeks at similar time periods.

At the finest grain, we obtain 237 clusters of users, 218 clusters of antennas and three time segments ,while week days remain as singletons. Again, the granularity prevent us from an easy interpretation, and we simplify the model. We keep 50% of informativity, that enables a reduction of the numbers of clusters of users and antennas to 40, and the numbers of groups of week days and hour segments to two. The week is divided in two parts: the working days and the weekend. For the hour dimension, the split occurs around 6 pm. The intervals are 0 am -6 pm and 6 pm -12 am. Note that the bound at midnight is artificial, because the day start as this time. The cut at 6 pm is the last in the hierarchy of the time segmentation. Then it would have been more relevant to consider a day from 6 pm to 6 pm the next day. Nevertheless, it is easier to have an interpretations on a "usual" time period between 0 am and 12 pm. Therefore we keep the following segmentation: 0 am -6 pm, 6 pm -12 pm. To illustrate the characterization of users' behaviors in terms of mobility provided by the grid, we focus on a group of users. The maps of Figure 8 shows the excesses and lacks of traffic in Abidjan during the week, for both periods of the day and for the selected group of users. The colors correspond to the mutual information M I(X M 1 ; X M 2 × X M 3 × X M 4 ) where X M 1 is the partition of antennas;

Definition 3 ( 1 and Y M 2 (

 312 Contribution to mutual information). Given the K -2 selected parts c i3...i K , the mutual information between two partitioned variables Y M from the partition M of Y 1 and Y 2 variables induced by the grid model M ) is defined as:

Fig. 5 :

 5 Fig. 5: Antennas activity clusters projected on Ivory Coast map. Colored clusters show inactivity periods while grey clusters indicate antennas whose traffic is complete over the period.

Fig. 7 :

 7 Fig. 7: Calendars of excesses (red) and lacks of calls emitted from each of the four clusters of antennas, in function of the weekday and the daytime.

  (a) Working days before 6 pm (b) Working days after 6 pm

Fig. 8 :

 8 Fig. 8: For a group of user, excesses and lacks of uses of antennas according to the day of the week and the time of the day. Focus on Abidjan.

-

  Yield management pricing strategy: a part of the pricing policy, called Bonus Zone, established in Ivory Coast offers discount prices (from 10% to 90%) to calling users depending on the location and hour of the emitting call. Maps and calendars resulting from the last two case studies on temporal distribution of output calls (see Section 5.3) and on mobility data (see Section 5.4) that are available at various granularities, provide valuable information to economic analysts in order to design optimized spatio-temporal pricing policy in the context of Bonus Zone.

	i.e. approximatively 55 millions records. The records are described by the four
	following variables: anonymized user ID (50000 categorical values); connexion
	antenna (1214 categorical values); time of call (minute precision); date fo call
	(from 2012/12/12 to 2012/12/24).
	From this data, we consider the user trajectories (identified by user ID) inside
	the network for the following analysis:
	1. Analysis of user mobility w.r.t. week day and hour. We consider the user ID,
	antennas, week day and hour. This data set can be considered as a set of
	spatio-temporal footprints, where each user ID is associated with a sequence
	of antenna usage over the time dimension. Here again, the time dimension
	is represented by two variables and the data of the whole period is folded up
	to week day and hour.
	Case studies on communication data. Communication data consists in 471 mil-
	lions mobile calls and covers a 5-month period (from 2011, December 1st to 2012,
	April 28th). The records are described by the four following variables: emitting
	antenna (1214 categorical values); receiving antenna (1216 categorical values);
	time of call (with hour precision); date of call (from 2011/12/01 to 2012/04/28).
	From this data set, we consider three subsets for:
	1. Analysis of call network between antennas. Considering emitting antennas,
	receiving antennas and the calls made between antennas, the data set can
	be seen as a directed multigraph where nodes are antennas and links are the
	calls between antennas.
	2. Analysis of output traffic w.r.t. date of call. We consider emitting antennas
	and the number of days for each call from referral to first day of recording.
	This data set can be considered as a temporal event sequence spanning over
	the whole observation period, where the time is the number of days passed
	and the events are the emitting antenna IDs.
	3. Analysis of output traffic w.r.t. week day and hour of call. We consider emit-
	ting antennas, the day of the week (stemming from the date and considered
	as a numerical variable) and the hour of the day for each call. Here the time
	dimension is represented by two variables and the data of the whole period
	are folded up to week day and hour.
	Case studies on mobility data. Mobility data consists in mobility traces of 50000
	users over a 2-week period (from 2012 December 12th to 2012 December 24th),

Data description & studies. The CDRs data under study come from the Orange D4D challenge 4 (Data For Development

[START_REF] Blondel | Data for development: the D4D challenge on mobile phone data[END_REF]

). We consider several case studies on two anonymized CDRs data sets from Ivory Coast, namely communication data and mobility data:

Romain Guigourès was with Orange Labs when this work began.

http://d4d.orange.com/en/home

X M

2 , the partition of the weekdays; X M 3 the discretization of the daytime; and X M 4 , the selected partition of the users. The selected group of users mainly connects to the antennas located in the East of Abidjan after 6 pm during the working days, while they rarely connect to the same antennas before 6 pm the same days. Then, it can be assumed that the selected cluster of users is composed people living in the same area. This hypothesis is reinforced by the socioeconomic nature of this part of Abidjan: it is a residential area. The contributions to mutual information of the other clusters of antennas are smaller. Three areas experience excesses of traffic before 6 pm and lacks after 6 pm. They correspond to the business district (Le Plateau), the embassies and universities neighborhood and the industrial zone located in the West of the city. The common feature of all these areas is their economic activity during the day. To sum up, we can assume that the users of the selected cluster are similar in that they live in the same area and work during the week in three localized area of Abidjan. Similar observations stand for several other clusters of users -thus we are able to summarize users' mobility behavior.

Conclusion

Motivated by two key points of economic development strategy of a telco in emerging countries, we have instantiated a generic methodology for exploratory analysis of CDRs data. Our method is based on a joint distribution estimation technique providing the user analyst with a summary of the data in a parameterfree way. We have also suggested several tools for exploring and exploiting the summary at various granularities and highlighting its relevant components. We have demonstrated the applicability of the method on graph data, temporal sequence data as well as user mobility data stemming from country-scale CDRs data. The results of the exploratory analysis are currently considered as valuable additional input to improve network planning strategy and pricing strategy.


Stefan Haar

Simon Theissing

Forecasting Passenger Loads in Transportation Networks

Keywords: Stochastic hybrid automata, Transportation networks, Fokker-Planck Equation

This work is part of an ongoing effort to understand the dynamics of passenger loads in modern, multimodal transportation networks (TNs) and to mitigate the impact of perturbations. The challenge is that the percentage of passengers at any given point of the TN that have a certain destination, i.e. their distribution over different trip profiles, is unknown. We introduce a stochastic hybrid automaton model for multimodal TNs that allows to compute how such probabilistic load vectors are propagated through the TN, and develop a computation strategy for forecasting the network's load a certain time into the future.

Introduction

We continue here the work begun in [START_REF] Haar | A hybrid-dynamical model for passenger-flow in transportation systems[END_REF] for capturing both the discrete vehicle movements and continuous passenger transfers in a multimodal public transportation network (TN). In [START_REF] Haar | A hybrid-dynamical model for passenger-flow in transportation systems[END_REF], a deterministic hybrid automaton (DHA) model was used, so as to overcome via fluidification the state space explosion that makes fully discrete models intractable. For the specification, we used both discrete and continuous Petri nets (PNs) as basic modelling blocks [START_REF] David | Discrete, Continuous, and Hybrid Petri Nets[END_REF], where the marking of the continuous places and the flows between them were vectorial instead of scalar. In fact, we integrated the numbers of passengers belonging to different trip profiles, i.e. having different destinations, as components of vector markings and -flows, with routing matrices relating them. Now -and this is the starting point of the present work -a real TN is everything but deterministic. On the one hand, there are highly unpredictable asynchronous events for which statistical data is hard to obtain. It is thus difficult to include them in the daily network operation [START_REF] Schmöecker | Metro service delay recovery: comparison of strategies and constraints across systems[END_REF], e.g. by means of minute-by-minute or hourly forecasts. Typical examples are passenger incidents. Now, note that -apart from few exceptions -such incidents originate locally, in one mode or line, and then propagate to other modes or lines by passenger transfers. These transfers are predictable, not necesseraly deterministic, if one knows the destination or trip profile of the passengers; but in general, this can only be known through probabilistic estimates. Finally, there are the "continuous" passenger arrival processes for which statistical data is easier to obtain: How many passenger will arrive at a station at what time? According to which route, including which vehicle missions, will they travel?

Here, we will extend the DHA model from [START_REF] Haar | A hybrid-dynamical model for passenger-flow in transportation systems[END_REF] in that we will replace all deterministic passenger arrival processes by their stochastic counterparts, and, in doing so, introduce a stochastic hybrid automaton (SHA) model with jumps between its discrete modes, at a priori equidistantly-spaced discrete points in time, defined a priori. The literature reveals many predecessors of our SHA model, notably in the past two decades, with every approach introducing the uncertainty at a different point in the model dynamics. For instance, the authors of [START_REF] Trivedi | Fspns: Fluid stochastic petri nets[END_REF] extended the dynamics underlying a PN-DHA model in that the jumps between the discrete modes are either exponentially distributed or immediate; with a weighting function as a means to resolve conflicts among simultaneously enabled immediate transitions [START_REF] Balbo | Introduction to generalized stochastic petri nets[END_REF]. However, in this modelling approach the discrete jumps are decoupled from the continuous states, and the latter evolve according (acc.) to deterministic differential balance equations; the authors of [START_REF] Horton | Fluid stochastic petri nets: Theory, applications, and solution techniques[END_REF] bridge the gap between the continuous states and the mode jumps by means of a guard function. Finally, the deterministic balance equations were replaced by normally distributed balance equations in [START_REF] Wolter | Modelling hybrid systems with fluid stochastic petri nets[END_REF]. Outside the framework of PNs, the authors of [START_REF] Hu | Towards a theory of stochastic hybrid systems[END_REF] introduced an SHA model that exhibits state-driven (forced) jumps between the discrete modes subjected to sets of stochastic differential equations (SDEs), one such set per mode. This approach was extended in [START_REF] Bujorianu | Toward a general theory of stochastic hybrid systems[END_REF] in that the mode transitions are no longer limited to forced jumps, but can be initiated by spontaneous jumps with state-dependent transition rates as well. The author of [START_REF] Bect | A unifying formulation of the fokker-planck-kolmogorov equation for general stochastic hybrid systems[END_REF], then showed how the SHA model from [START_REF] Bujorianu | Toward a general theory of stochastic hybrid systems[END_REF] can be formulated in an equivalent system of integro-differential equations together with boundary conditions. Also notice that the authors of [START_REF] Prandini | A numerical approximation scheme for reachability analysis of stochastic hybrid systems with state-dependent switchings[END_REF] presented a grid-based asymptotic approximation method for a backward reachability problem subjected to the dynamics of an SHA model that encounters spontaneous jumps between its discrete modes; with a system of SDEs assigned to every mode. That system is approximated by a Markov chain following a space and subsequent time discretization; whereas in our approach the discretization of the time precedes the discretization of the space, and the latter comes along with a numerical integration of the continuous states in a discrete mode.

In the rest of this paper, we will discuss the specification of TN's infrastructure in our vectorial SHA model, and the vehicle operation as well as the routing of all passenger flows thereon (Sec. 2 on p. 3). We will then elaborate the SHA model's time-continuous dynamics, before pinpointing all mode transition times to an equidistantly-spaced mesh, which can be regarded as a first major step to render forecasts of the model's hybrid state feasible (Sec. 3 on p. 8). Next, we will integrate the SHA model's discrete-time approximation into a computation strategy which can be used to forecast TN's passenger loads (Sec. 4 on p. 16). Finally, we will provide a short summary and an outlook on future work (Sec. 5 on p. 19).

Capturing the Transportation Network at Hand

TN's Infrastructure

We capture TN's infrastructure in our automaton model by a finite set of transportation grids that accommodate TN's different lines and transportation modes (Sec. 2.1.1), another finite set S of stations (Sec. 2.1.2 on p. 4), and an interface that defines all passenger transfer possibilities between the vehicles operated in the transportation grids and the stations (Sec. 2.1.3 on p. 4).

A Transportation Grid

Every transportation grid accommodates the lines of a particular transportation mode such as the rail grid of a metro system: It is made up of a finite set of discrete waypoints that are connected by tracks. In particular, a track that connects the waypoint w 1 to another waypoint w 2 indicates the possibility of a vehicle movement from w 1 to w 2 ; no more (such as a length, a curvature, or a slope) and no less. Thus, every pair of a waypoint is connected by at maximum two tracks, namely one track in each direction. Since by convention every waypoint moreover is either empty or occupied by at maximum one vehicle at a time, then also note that conflicting vehicle movements can occur iff two or more vehicles try to access the same waypoint at the same time.

We specify a particular transportation grid, say g, in form of a Petri net model (with the 'token flow' left out), in which all waypoints (= places) in g are represented as simple circles and all tracks (= transitions) as simple boxes. One edge then connects w 1 to the track t and another edge t to w 2 iff t shall implement the possibility of a vehicle movement from w 1 to w 2 . Moreover, we enumerate every track in order to specify the resolution of all possible conflicting vehicle movements in a deterministic manner.

As an example, consider Fig. It schematically depicts the transportation grid for the orbital line of a people mover: A vehicle at w 1 can go to w 2 via t 1,2 , from w 2 to w 3 via t 2,3 , and from w 3 back to w 1 via t 3,1 , respectively. Alternatively, a vehicle at w 1 can take t 1,3 , and thus directly move from w 1 to w 3 without the necessity to pass by w 2 . Then note that a conflict between two vehicle movements might occur iff one vehicle, say a 1 , wants to access w 3 from w 2 via t 2,3 , and at the same time another vehicle, say a 2 , wants to access w 3 from w 1 via t 1,3 . Assuming this to be true, we privilege a 2 over a 1 due to the fact that a vehicle movement taking t 1,3 has a higher priority (lower integer number inscribed to t 1,3 ) than another vehicle movement taking t 2,3 .

A Station

Every station s ∈ S is made up of a finite set of capacity-limited gathering points (GPs) that accommodate the transferring and waiting passengers, and another finite set of elements, called corridors The station has two GPs and three corridors: Among the two GPs there is one entrance area that can accommodate up to 50 passengers (written in brackets next to it), and one platform that can accommodate up to 200 passengers. Moreover we can say that a passenger at the entrance area can transfer to the platform since a corridor connects both in the right direction. The unambiguous interpretation of the remaining two corridors depends on the context of the station, though; i.e., on the interface specification between this station and all transportation grids of TN; see below. Nevertheless, from the labels written next to all corridors we can anticipate that the entrance area accommodates all passengers who enter TN through this station, and that all passengers who leave the platform do so in order to board a vehicle.

Connecting the Transportation Grids and the Stations

The interface between the transportation grids and the stations defines which corridor in a station is connected to which waypoint in a transportation grid. In doing so, it defines all passenger transfer possibilities between the GPs of the stations and the vehicles that are stopped at particular waypoints in the transportation grids in our SHA model. As an example, look at Fig. 3.

Station Transportation Grid

Platform p It depicts the interface between a simple station, say s, and an extract of a transportation grid: If a vehicle a is stopped at the waypoint w i , then all passengers on-board a can alight from it to the single platform p in s through c a . Similarly, the passengers can board a from p through c b .

In general, the interface between all stations and transportation grids is specified in form of a set of edges that connects the corresponding graphs, in which every edge either connects a corridor c in a station to a waypoint w in a transportation grid indicating the possibility of a passenger transfer from the single parent vertex of c (= GP) to a vehicle that is stopped at w, or vice versa.

We refer to this combined graph (of all stations, transportation grids, and the interface between both) as the infrastructure of our model, and demand that the passengers can access a vehicle a that is stopped at a particular waypoint only from one GP of a particular station, say s, if at all. Accordingly, all passengers on-board can alight from a to one (not necessarily the same) GP of s if at all.

TN's Vehicle Operation

At any time every vehicle (out of a finite set of vehicles) is either in the driving condition parked, stopped, or driving. Now every vehicle in operation, i.e., stopped or driving, executes a mission (Sec. 2.2.3 on p. 6) in a run (Sec. 2.2.2 on p. 6) acc. to a dispatch plan (Sec. 2.2.1).

The Dispatch Plan

The central element of the vehicle operation in our SHA model is a dispatch plan that defines at what time which vehicle is supposed to process which run in form of an ordered list. It can cover part of a day, a complete day, or even more. Moreover, note that in practice different dispatch plans might have been specified for different operational modes (nominal operation, perturbed mode of operation, etc.). However, in our model there is only one dispatch plan, which requires us to merge two or more real-world dispatch plans if necessary; or to compute for every different dispatch plan a separate forecast.

A Vehicle Run

Every vehicle run specifies a sequence of vehicle missions that have to be executed in the given order by a vehicle that processes the run; which should result in a coherent driving cycle of fixed-route dead headings and transportation services.

A Vehicle Mission

The fixed route that underlies every vehicle mission is specified as a path in a transportation grid. An indication of stops at the waypoints along that path, minimum and maximum dwell times for these stops, and deterministic driving times in between the waypoints then complement the specification of a particular vehicle mission.

As an example, look at Fig. 4. It depicts the graphical specification of a sample vehicle mission, say x , for the orbital transportation grid from Fig. 1 on p. 3: A vehicle a that executes x is supposed to pass trhough the orbital grid starting from and ending at w 1 . Moreover, note that a is supposed to stop at w 1 and w 3 since a separate pair of a minimum and a maximum dwell time is assigned to both (comma-separated values written in brackets next to w 1 and w 3 ). In particular, a cannot depart from w 3 before its dwell time at w 3 has exceeded 10 seconds. Assuming that passengers can board a/alight from a at w 3 , then a departs from w 3 before its dwell time at w 3 has exceeded the maximum dwell time of 60 seconds iff no more passengers want to board/alight from it. On the contrary, a departs from w 3 after the maximum dwell time of 60 seconds has elapsed iff some passengers still want to alight from it or its next waypoint is blocked by another vehicle. In this way, the minimum dwell times impose hard constraints on the vehicle operation, whereas the maximum dwell times have to be respected iff possible.

Wayp. w 1 [10, ∞] Track t 1,2 [52] Wayp. w 2 Track t 2,3 [72] Wayp. w 3 [10, 60] Track t 3,1 [41] Wayp. w 1 [10, 0]

The Routing of All Passenger Flows in TN

Look at Fig. 5 which depicts an extract of a sample infrastructure (graph) together with the specification of three different trip profiles 1, 2, and 3: All passengers who travel acc. to the first trip profile enter TN through the station S 1 , and arrive at its platform next. From this platform, they then would like to board a vehicle that executes the mission x 1 , which passes through and stops at all three waypoints w 1 , w 2 , and w 3 . On-board such a vehicle stopped at w 2 , these passengers would like to alight from it to a GP in the station S 2 (not depicted). In contrast to the first trip profile, all passengers who travel acc. to the second or third trip profile would like to alight from a vehicle that is stopped at w 3 to a GP in the station S 3 (not depicted). Then note that, in principle, both vehicle missions x 1 and x 2 provide a transportation service from S 1 to S 3 . However in contrast to x 1 , a vehicle that executes x 2 is supposed to skip w 2 and thus provides an express service from S 1 to S 3 . Now the difference between the passengers who travel acc. to the second trip profile from the passengers who travel acc. to the third trip profile is that the latter show a preference for x 2 ; whereas all passengers of the second trip profile are indifferent to both vehicle missions and thus board whatever vehicle mission is available.

Station S 1

Transportation Grid

Station S 2

Station S 3

Platform Access 1,2,3

x 1 : 1,2 x 2 : 2,3 Board 1,2,3 Wayp. w 1 1,2,3 x 1 , x 2 Track t 1,2 1,2,3
Wayp.

w 2 1,2,3 x 1 , x 2 Track t 3,1 x 1 , x 2 Track t 2,3 2,3
Wayp. w 3 2,3

x 1 Alight 1 1 x 1 , x 2 Alight 2,3 2,3
Fig. 5. Sample infrastructure (graph) of the SHA model together with the specification of three different trip profiles (1, 2, and 3), in which all passengers of the trip profile 3 prefer to board a vehicle that executes the mission x 2 in order to travel from S 1 to S 3 over a vehicle that executes x 1

In general, every passenger travels acc. to a particular trip profile out of a finite set of trip profiles Y. Every trip profile y ∈ Y in turn specifies a path in the SHA model's infrastructure together with labels assigned to every corridor that is dedicated to the boarding of a vehicle so as to account for the passengers' different mission preferences. Then with goal to ease the graphical specification of all trip profiles, we can inscribe the infrastructure graph with the paths and the stops specified for the different vehicle missions. So in the example from above, the fact that x 1 is inscribed to the corridor "Alight" tells us that a vehicle that executes x 1 stops at w 2 ; which in turn implies that all passengers on-board that vehicle can alight from it in the given direction (or the GP this corridor points to; not depicted).

From the graphical specification of all trip profiles, we deduce routing matrices; for every corridor in every station one. Refer to Tab. 1 for an overview of some shorthands pertaining to matrix/vector operations used here and throughout the rest of the paper: Initially, every routing matrix R (•) is a diagonal matrix with zeros and ones only. More precisely, the element R (c) [i, j], with i, j ∈ Y, of the routing matrix R (c) assigned to the corridor c is one iff i = j and the trip profile i ∈ Y covers c; in which we assume that all trip profiles are enumerated from 1 to |Y|. We next might rewrite the columns of some routing matrices to account for the re-routing of some passenger flows. As a result, these (constant) routing matrices are not diagonal any more, in which their interpretation/specification goes as follows: R (c) [i, j] is the relative amount of passengers who enter c acc. to the trip profile i and who leave it acc. to the trip profile j. We then demand that every column R (c) [•, j], with j ∈ {1, 2, . . . , |Y|}, of R (c) either sums up to one or to zero, which will ensure the conservation of all passenger flows in the set up of the balance equations at a later point.

Table 1 An overview of some shorthands pertaining to matrix/vector operations used throughout the paper

Shorthand

Meaning

a[i] Element in the i-th row of a column vector a A T Transpose of a matrix A A[i, j]
Element in the i-th row and the j-th column of a matrix A A[•, j] j-th column of a matrix A

Capturing the Probability Flow

Our automaton model belongs to the class of continuous-time hybrid-state automata with deterministic-time and state-driven probabilistic jumps among a finite set of discrete modes. In general, every mode in such a model refers to a discrete state to which is assigned to a continuous domain with a continuous differential dynamics2 . An edge-labelled graph -here and in the following referred to as a mode graph -captures all jump conditions between the different modes that might have been defined in terms of the continuous states that have to enter particular target sets or time-thresholds that have to elapse. Now in our use case, we limit every mode to a discrete state in TN's vehicle operation that tells us

• the driving condition of every vehicle (parked, stopped, or driving),

• the position of every vehicle, in which, by convention, the position of a driving vehicle is identical to the waypoint it moves towards, and

• the dispatch, the run, and the mission of every stopped and every driving vehicle.

Then note that for every mode we can define a set of balance equations that capture the SHA model's continuous passenger flow dynamics in this particular mode in that they relate the passenger loads of the GPs among themselves and some stopped vehicles (Sec. 3.1 on p. 9). We can next use these balance equations to numerically integrate estimations for all passenger loads so as to evaluate the probability of occurrence of the SHA model's transition to another mode (Sec. 3.2 on p. 10). However at the latest, we then notice that we have to confine all mode transition times to finite sets (Sec. 3.3 on p. 12) for any feasible computation (Sec. 3.4 on p. 14).

The Balance Equations for a Particular Mode

Remember that we impose a very crucial constraint on the specification of the interface for the passenger transfers between the GPs in the stations and the vehicles stopped at the waypoints of the transportation grids (cf. Sec. 2.1.3 on p. 4): If a vehicle a is stopped at a waypoint, then a passenger transfer to/from a is possible from/to one and the same station only; which does not imply that the GP that accommodates all passengers who alight from a must be identical to the GP that accommodates all passengers who board a, though. In this context, we say that a vehicle a is docked to a station s ∈ S iff (i) a is stopped at a waypoint that can be accessed from a GP in s acc. to the specification of the infrastructure's interface, and (ii) some passengers want to alight from or board a acc. to the specification of the passenger (re-)routing. For a particular mode of our automaton, we can thus decompose the set of all balance equations for all passenger loads -which define the automaton's continuous passenger flow dynamics as long as no jump to another mode occurs -w.r.t. the different stations: For every station s ∈ S we set up a (decoupled) system of n := n 1 + n 2 (Itô-) stochastic differential balance equations for the |Y|-dimensional passenger load vectors3 of all n 1 ∈ N >0 GPs in s and all n 2 ∈ N ≥0 vehicles docked to s. For time τ > 0, this system can be written down in the form of dX

s (τ ) = A s (X s (τ )) dτ + B s dW s (τ ), (1)
with the (n × |Y|)-dimensional state vector X s that groups all n 1 +n 2 passenger load vectors; the n-dimensional drift vector A s , and the (n × m)-dimensional diffusion matrix B s iff we assume that m ∈ N >0 trip profiles enter TN through s; in which case W s is an m-dimensional Wiener process4 [START_REF] Szabados | An elementary introduction to the wiener process and stochastic integrals[END_REF]. In other words, we assume that the uncertainty that is inherent to the continuous passenger flow dynamics in form of (1) comes along with the passengers who enter TN through its stations (= probabilistic passenger arrival processes). However, all passenger flows that originate within TN such as a flow between two GPs, a flow between a GP and a stopped vehicle, or a passenger flow that leaves TN from one of its stations, are supposed to be specified in a deterministic manner given the respective passenger loads.

With that said, we now look at the exemplary specification of the balance equations for the station s from Fig. 3 on p. 5 in a particular mode, in which we assume that one vehicle, say a, is stopped at the waypoint w i , and some passengers want to alight from/board a to/from the platform p: First of all, note that every corridor c ∈ C s := {c a , c b , c e , c A } in s must have a maximum throughput φ max (c), with φ max : C s → R ≥0 that limits the number of passengers per second who can pass through it. This maximum throughput thus limits the (scalar) magnitude (

φ m (c, τ ), with φ m : C s × R ≥0 → R ≥0 , of the passenger flow (vector) φ (c, τ ), with φ : C s × R ≥0 → R |Y| ≥0 , through c at time τ ≥ 0 s.t. φ m (c, τ ) = i∈Y φ (c, τ ) [i]
) 2
Note that the passenger load of a must neither be negative nor exceed a's capacity limit κ a > 0. For any τ ≥ 0 and any i ∈ Y, we thus have to require that M (a, τ ) [i] → 0 imply φ (c a , τ ) [i] → 0, and j∈Y M (a, τ ) [j] → κ a implies φ m (c b , τ ) → 0; in which case we say that φ (c a , τ ) is demand-sensitive, and φ (c b , τ ) is capacity-sensitive. Moreover, we need the passenger flow assigned to every corridor to be routing-proper, i.e., for any c ∈ C s , j ∈ Y, and τ ≥ 0, it holds that i∈Y R (c) [i, j] = 0 implies φ (c, τ ) [i] = 0. Now the specification of the balance equation for the passenger load of the platform p is a little bit more complex, though, due to the fact that we have to integrate into it the uncertainty that comes along with the passenger flow φ (c A , τ ) through c A (= arrival process). This is captured by adding, to the impact that φ (c A , τ ) has on M (p, τ ), as expressed by the product R (c A ) φ (c A , τ ), another product term R (c A ) D (c A ). This latter term connects the multidimensional Wiener process W s with M (p, τ ) through the product of R (c A ) with the diffusion matrix

D (c A ) 8 , with D : c A → R |Y|×|Ws| : dM (p, τ ) :=
Drift-and diffusion term for the spec. of the arrival process through

c A R (c A ) φ (c A , τ ) dτ + D (c A ) dW s +R (c a ) φ (c a , τ ) dτ - φ (c e , τ ) dτ -φ (c b , τ ) dτ (3)
In (3), we assume that all passenger flows are routing-proper as well as -iff applicable -capacity-and demand-sensitive, in a way that is similar to the specification of the passenger flows in the balance equation ( 2) for the passenger load of a above. However, note that this requirement taking alone cannot ensure the non-negativity and the capacity limit of p any more. Instead, we have to explicitly enforce both during the numerical integration of (1), by defining appropriate boundary conditions. Similarly, we have to explicitly enforce the maximum throughput of c A .

Numerical Integration of Balance Equations

In practice, there are two dominant approaches to computing the temporal change of an initial distribution 9 subjected to (1). On the one hand, there is the Monte Carlo method [START_REF] Mackay | Introduction to monte carlo methods[END_REF]: From the pool of all possible solution paths, some paths are selected "randomly" and possibly several times. In this approach, one thus would have to ensure that the correlation between the selected solution paths matches the correlation found in the actual stochastic process. Of course, this is a difficult task especially if the statistical data is not at hand, as is often the case in degraded modes of operation. The second approach is the one we will use here: Integrate the balance equations from (1) into systems of linear parabolic partial differential equations (one for every station s ∈ S):

∂ ∂τ pdf (X s (τ )) = - N i=1 ∂ ∂X s,i A s (X s (τ ) , τ ) [i] pdf (X s (τ )) + 1 2 N i=1 N j=1 Ψ s [i, j] ∂ 2 ∂X s,i X s,j pdf (X s (τ )) . (4)
with A s , B s from (1) and the abbreviation Ψ s := B s B T s . This system is also known as the (multidimensional) Fokker-Planck equation (FPE) or the Kolmogorov forward equation; it describes the time evolution of the probability density function (pdf) for X s .

Eqn. ( 4) defines a continuity equation in form of a system of deterministic parabolic partial differential equations: By introducing the probability flux

f (X s (τ )) := A s (X s (τ )) pdf (X s (τ )) - 1 2 Ψ s         ∂ ∂X s,1 pdf (X s (τ )) ∂ ∂X s,2 pdf (X s (τ )) . . . ∂ ∂X s,N pdf (X s (τ ))         , (5)
we can rewrite (4) acc. to

∂ ∂τ pdf (X s (τ )) + div (f (X s (τ ))) = 0, (6)
in which div (•) denotes the divergence operator. Now for this continuity equation we can derive reflecting boundary conditions for the numerical integration of (4) that confine the probability flux to a closed convex polytope K s , which captures the complete passenger load space associated with s in the mode at hand10 . Here, we sketch only the major steps in their derivation in an informal manner, and refer to the literature [5, chapter 5] for more details. The starting point is the insight that the cumulative probability of X s adopting a value in K s at any time instant τ ≥ 0 must be one. Thus, the time derivative of this cumulative probability must be zero. Then, after some transformations employing the divergence theorem, we obtain the reflecting boundary condition

n (X s (τ )) • f (X s (τ )) = 0, ∀ X s (τ ) ∈ ∂K s and τ ≥ 0, (7)
in which n (X s (τ )) denotes a unit vector in the outward orthogonal direction of K s w.r.t. the boundary ∂K s of K s at X s (τ ) ∈ ∂K s . For Cartesian coordinates the scalar product of f (X s (τ )) with n (X s (τ )) simplifies to n (X s (τ )) T f (X s (τ )).

From this, it is possible to include ( 4) and [START_REF] Horton | Fluid stochastic petri nets: Theory, applications, and solution techniques[END_REF] subjected to (1) in form of a boundary value problem in a numerical integration scheme that preserves the conservation of the probability flux. We do not go into further details here, but mention that we go for the finite volume method in combination with the Crank-Nicolson method in an on-going implementation. Finally, note that we can enforce the maximum throughput of every corridor that implements a passenger arrival process, in adjusting the barriers for some reflecting boundary conditions in our time-discrete computation strategy; more on this later.

Discretization of Time

Remember that we resolve conflicting vehicle movements in a deterministic manner if necessary. Thus, two different types of mode transitions can occur in our automaton model in principle: deterministic-timed-and probabilistic mode transitions.

Deterministic-timed Mode Transitions. They involve vehicle dispatches, vehicle arrivals, and vehicle departures assuming that the latter are not tied to any passenger loads. They are thus completely conditioned on (time) constraints that are imposed on the vehicles' dispatch times, as well as their elapsed driving-and dwell times (the driving time of vehicle a 1 must equal τ 1 seconds, the dwell of vehicle a 2 must exceed τ 2 seconds, and so forth).

Probabilistic Mode Transitions. They involve departures of vehicles that are docked to a station because less than one passenger wants to alight from/board a vehicle, or less than one passenger can board a vehicle (because of the vehicle's limited capacity), and so on. They are thus conditioned on constraints that are imposed on some passenger loads (= passenger load constraints) in form of closed convex polytopes; although they might involve time constraints as well. Moreover, note that the departure events of two vehicles, which are docked to two different stations, are mutually exclusive due to the set up of our balance equations (for every station one separate and decoupled system).

Vehicle Load Tree. The possible departure times of vehicles docked to a station form an uncountable set, which renders the computation of our continuous-time automaton model intractable. In order to overcome this burden, we thus discretize the time in the computation of the model's hybrid state. More specifically, we compute the time evolution of the model's vehicle load starting recursively from some initial one from one equidistantly-spaced discrete point in time to the next; in which a particular vehicle load refines a mode of our automaton model in that it also defines elapsed dwell times for all stopped vehicles, and elapsed driving times for all driving vehicles. Thus, several vehicle loads can belong to one and the same mode. We then a) assume that the SHA model can change its vehicle load only at equidistantly-spaced discrete points in time. In addition, we b) replace every equality constraint for a vehicle arrival to be true by an inequality constraint s.t. the elapsed driving time of a vehicle has to exceed a certain time threshold τ > 0; and does not have to equal τ . Finally, we c) ceil every dispatch time (from the dispatch plan) to the next closest discrete point in time. As a result, we obtain a vehicle load tree G; an edge-and vertex-labelled directed rooted tree (cf. Fig. 6): Every vertex captures a particular vehicle load L of TN. This L is stored as a tuple (L , n), in which L is a vehicle load itself, and n ∈ N >0 is a non-negative integer.

We then obtain L from (L , n) in that we increment every elapsed driving-and dwell time in L by n. The vertex label assigned to L gives the discrete time step that elapses before our automaton model -if at all -is in L starting from the vehicle load in the root of G at τ = 0; it equals the height of (L , n) in G. On the other hand, every edge label, say from the vehicle load L 1 to L 2 , gives a convex polytope for TN's passenger load that equals the complete passenger load space 11 iff the transition of the automaton's vehicle load from L 1 to L 2 is not conditioned on a probabilistic mode transition, and a subset of the passenger load space otherwise; in which (latter) case this subset is adopted from the specification of the passenger load constraint for the respective mode transition. It then follows that branches in the vehicle load tree arise from probabilistic mode transitions only.

L 0 L 1 τ > 2.9 ∆ τ , x ∈ X L 2 τ ≥ 0 ∆ τ , x ∈ X 2 (a) Mode graph L 0 , 0 0 L 0 , 1 X 1 L 2 , 0 X 2 1 L 0 , 2 X 2 L 2 , 0 X 2 2 L 1 , 0 X 3 (b)
Vehicle load tree Fig. 6. Comparison of a mode graph and a vehicle load tree, in which X denotes the complete passenger load space, X 2 ⊂ X, and ∆ τ > 0 is the time step for the equidistant discretization

Computation of the Vehicle Load Tree. The iterative computation of a vehicle load tree from one time layer to the next, with a fixed time step of ∆ τ > 0 seconds, is straightforward; and we thus describe the major steps of every (iterative) run only in an informal manner here: Starting point is the tree G 0 = (V 0 , E 0 , f e,0 , f v,0 ) with a single vertex v 0 s.t. V 0 = {v 0 } and E 0 = ∅. This vertex captures the automaton's vehicle load L 0 at time τ = 0 s.t. v 0 = (L 0 , 0) and f v (v 0 ) = 0. For a summary of the shorthands V i , E i , f v,i , and f e,i pertaining to the specification of the vehicle load tree G i , with i ∈ N ≥0 , refer to Tab. 2 on p. 14. The run i > 0 then computes G i from G i-1 in that it first adopts the specification of G i-1 s.t. G i := G i-1 , and then processes for every vertex v ∈ andM is the mode that is contained in L the following 5 steps: (i) Compute the set F of all enabled mode transitions from M to any other mode M , in which a mode transition from M to M is enabled iff all time constraints (independently of the passenger load constraints) are met until time step i (ii) Define F d ⊆ F as the subset of all deterministic-timed mode transitions, and F p ⊆ F as the subset of all probabilistic mode transitions s.t.

V i-1 s.t. f v (L, n) = i -1, v := (L, n),
F = F d ∪ F p and F d ∩ F p = ∅ (iii) If F is empty: (a) Define v := (L, n + 1)
, and X as the complete passenger load space

(b) Append v as a child of v to G i (c) Define f e,i (v, v ) := X, and f v,i (v) := i (iv) If F d = ∅:
(a) Define X as the complete passenger load space (b) Identify the mode transition f ∈ F d with the biggest impact on L in terms of the number of vehicles that change their operational states 12 (position, mission, etc.) (c) W.r.t. M , let M be the mode following f , and A be the set of all vehicles that change their operational states due to f (d) Compute the vehicle load L :

• Adopt M for the mode that is contained in L • Set the elapsed driving/dwell time of every driving/stopped vehicle a ∈ A in L to zero • For the elapsed driving/dwell time of every other driving/stopped vehicle a ∈ A in L , take the corresponding time from L and increment it by one (e) Define v := (L , 0) (f) Continue with step iii.b (v) If F = F p , then do for every f ∈ F :

• Define X as the convex polytope that captures the passenger load constraint imposed on f • Continue with step iv.c Table 2 Shorthands pertaining to the computation of a vehicle load tree G assuming that V is the vertex set of G, and E is its edge set Shorthand Meaning

f v (v) Height of v ∈ V in G f e (v, v )
Subset of the complete passenger load space that has to be entered by the automaton's passenger load so as to induce a jump from v to v , with (v, v ) ∈ E

• v Set • v := {v ∈ V : (v , v) ∈ E} of all parents of v ∈ V v • Set v • := {v ∈ V : (v, v ) ∈ E} of all children of v ∈ V

Computing Transition Probabilities

So far we have described how estimations for the passenger loads of all GPs in a station, and of all vehicles docked to that station, can be propagated in time in a particular mode of our automaton model by the numerical integration of deterministic partial differential equations. We then introduced the vehicle load of TN, and described how -starting from some initial vehicle load L 0 at time τ = 0 -all possible solution paths of that vehicle load can be computed in a time-discrete manner, in which all mode transition times are confined to an equidistantly-spaced time-layered mesh. The next step is thus to connect both; the continuous-time propagations of all passenger load estimations in the different modes on the one hand, and all possible solution paths for the vehicle load on the other hand. In doing so, we describe how a discrete-time approximation for the hybrid state (= particular vehicle load of TN + all passenger load vectors) of our (continuous-time) automaton can be computed 13 . Note first that a missing edge from v ∈ V to v ∈ V in the automaton's vehicle load tree G = (V , E, f e , f v ) implies that the corresponding transition of the automaton's vehicle load cannot occur, given the choice of the discrete time step ∆ τ > 0. On the other hand, the existence of this edge individually only implies that the corresponding transition may occur, but does not necessarily have to do so. In this context focussing on the SHA model's discrete state first, note that

P v := v ∈ • v P v |v P (v) (8)
must hold for every v ∈ V , in which P (v ) is the marginal probability that the vehicle load of the automaton at time step f v (v ) is v , and P (v |v) is the (marginal) conditional probability that the automaton's vehicle load at

f v (v ) is v , given that it is v at f v (v) = f v (v ) -1.
This conditional probability of course depends on the automaton's passenger load; see below. Next note that at any time step, our automaton model must be in some vehicle load, which is identical to saying that v ∈ fv -1 (i)

P (v) = 1, ∀i ∈ N ≥0 . (9)
Now coming to the integration of the SHA's continuous state into the discrete framework for its vehicle load defined by ( 8) and ( 9) and thus to the computation of the SHA's complete hybrid state, note that at any time step i ∈ N ≥0 , the automaton's passenger load must be in a non-empty subset from the complete passenger load space with a positive probability that we store in the density pdf (i). Starting from the automaton's initial hybrid state in form of pdf (0) for its passenger load and v 0 := (L 0 , 0) for its vehicle load, with P (v 0 ) = 1, v 0 := (L 0 , 0), and f v (v 0 ) = 0; the automaton's hybrid state at i ∈ N >0 computes from the preceding one at i -1 as follows:

• The automaton's passenger load at i computes acc. to

pdf (i) = v ∈ f -1 v (i-1) pdf o (v) P (v) . (10)
Therein, pdf o (v), with v ∈ V , is the result of the numerical integration of (4) on p. 11 from time step f v (v) to f v (v) + 1; with pdf (f v (v)) as initial density; (4) subjected to (7) on p. 11; and (4) set up for all stochastic balance equations (for every station s ∈ S one) acc. (1) on p. 9 for the mode that is contained in v.

• The probability that the automaton is in the vehicle load v ∈ V , with f v (v ) = i, can be computed from the probability P (v) that the automaton is in v ∈ • v at f v (v) = i -1 and the conditional probability P (v |v) acc. to [START_REF] Hu | Towards a theory of stochastic hybrid systems[END_REF]. Therein, P (v |v) captures the extend acc. to which the automaton's passenger load enters the target set f e (v, v ) at i assuming that the automaton's vehicle load is v. It computes acc. to

P v |v := fe(v,v ) pdf o (v) dV. ( 11
)
4 Forecasting Passenger Loads

The Problem Statement

We start from an observation of TN's (network) state at time τ = 0 and a question about its passenger load (= objective of the forecast). The former is captured by an exact knowledge of its vehicle load L 0 at τ = 0 and an estimation of its passenger load (= passenger load of every GP in every station, and of every vehicle operated in the infrastructure at hand) for τ = 0 in form of the density pdf (0). From the latter, we can deduce a prediction horizon τ * , a target set for TN's passenger load X * , and a threshold α * for the probability that TN's passenger load adopts a value from X * . For instance, a particular question might read as follows:

"Will the passenger load of the platform p in the station s ∈ S exceed 200 passengers (→ X * ) with a probability greater than 0.7 (→ α * ) within the next 20 minutes (→ τ * )?"

The General Solution

We compute α(0) := X * pdf (0) dV , i.e., the probability of TN's passenger loads to take on a value from X * at τ = 0. Then, proceed as follows: If α(0) > α * , we have already obtained an affirmative answer to our question and we can stop here. Otherwise, we make use of our time-discrete automaton model in that we provide α(0), the time index i = 0, and the vehicle load tree G 0 := (V 0 , E 0 , f e,0 , f v,0 ) as input to the following iteration that terminates iff α(j) > α * or (j + 1) ∆ τ > τ * , with j ∈ N ≥0 ; here V 0 := {v 0 }, with v 0 := (L 0 , 0), E 0 := ∅, and f v,0 (v 0 ) := 0: (i) Compute the vehicle load tree G i+1 from G i as described in Sec. 3.3 (ii) Compute pdf (i + 1) acc. to [START_REF] Prandini | A numerical approximation scheme for reachability analysis of stochastic hybrid systems with state-dependent switchings[END_REF]

on p. 15 (iii) Compute α(i + 1) := X * pdf (i + 1) dV (iv) If α(i + 1) ≤ α * , then compute P (v) acc. to (8) on p. 15 for every v ∈ V i+1 s.t. f v,i+1 = i + 1 (v) Set i := i + 1
Thus by construction, we know the binary true/false answer to our initial question upon termination of this iterative computation. So the next question then is whether or not the iterative computation itself is feasible in a reasonable amount of time (number of trivial computation steps) with reasonable memory constraints, which is probably not the case unless some further simplifications are made. This will lead us to a more efficient implementation of the iterative computation of the forecast from above (Sec. 4.3) together with a re-definition of all balance equations for the numerical integration of the passenger flows (Sec. 4.4 on p. 17).

The Efficient Implementation of the Iterative Forecasting Algorithm

First of all, note that in order to compute pdf (i + 1), with i ∈ H * and H * := 1, 2, . . . , τ * ∆ τ

, acc. to Eqn. 10 on p. 15, we only need to know the set of all leaf nodes of the vehicle load tree G i , say V i , together with the density pdf (i -1) and P (v) for every v ∈ V i . So, instead of computing a vehicle load tree in the iterative computation of the forecast, as described above, anew every time, it makes sense to propagate only this reduced piece of information from one iteration to the next. However, taking this measure alone, we still end up with a potentially huge set V i given even only very few passenger ride possibilities on-board concurrently moving vehicles; and this is a big problem since we do not only have to compute the set V i itself, but for every v ∈ V i we also have to numerically solve partial differential equations in the form of Eqn. 1 on p. 9 (so far, one for every station). What we have to do thus is to identify all vehicle movements that do not affect the forecast and to ignore them. In this context, note that -starting from the time step (i -1) ∆ τ until the end of the prediction horizon at time τ * -we can compute how far the operational state (including the position) of every vehicle can evolve from the constant driving times for the individual tracks and the minimum dwell times for the waypoints in the specification of the different vehicle missions. Combining this knowledge with the specification of the complete passenger re-routing and the definition of the target set X * -from which we can deduce a target zone in form of a set of GPs in the stations and vehicles, -we can thus compute whether or not a passenger outside this target zone can penetrate it within the prediction horizon, given the passenger's current position. In this way, we can compute at every iterative step i ∈ H * not only the set of all vehicle movements that we have to take into consideration in the computation of V j , but also for every v ∈ V j the set of all stations that we have to include in the numerical integration of the passenger loads; that is, excluding those stations that are outside the target zone and from which passengers cannot penetrate the target zone within the remaining prediction horizon.

Redefining all Balance Equations for the Numerical Integration

In the worst case, the system of balance equations derived for a station s ∈ S in a particular mode of the automaton acc. to (1) on p. 9 has the dimension (n 1 + n 2 ) |Y|, in which n 1 is the number of GPs in s, n 2 is the number of vehicles docked to s, and Y is the finite set of different trip profiles introduced during the specification of the passenger routing. Then note that in practice (1) with even very few dimensions renders the numerical integration of the boundary value problem in form of (4) on p. 11 subjected to (7) on p. 11 intractable given reasonable computation constraints. Now one approach to overcome this computational burden is to decouple all passenger flows during the numerical integration of all passenger loads associated with s, and to merge trip profiles with common last miles: First of all, we assume that all n := n 1 + n 2 passenger load vectors associated with s were estimated for the forecast at hand in isolation in form of separate and uncorrelated probability densities. Moreover, we assume that these estimations do not violate any additional constraints (non-existent in the original setting of the balance equations) w.r.t. the maximum number of passengers who can be located at a particular GP or on-board a vehicle as a function of their trip profile; see below. Assuming this to be true, we then integrate these estimations into a finite set of decoupled one-and twodimensional systems of balance equations for the passenger loads associated with s, which we classify either as arrival processes, passenger transfers, or outflows: Arrival Process. For every corridor in s that implements a passenger arrival process, say c, to one of its GPs, say p, set up a one-dimensional balance equation acc. to (1) for the free capacity of p. This free capacity can be computed from p's (constant) total capacity minus its cumulative passenger load. Then note that we can ensure the maximum throughput φ max (c) assigned to c during the numerical integration of this one-dimensional balance equation by reducing the total capacity of p s.t. it equals the product of φ max (c) with ∆ τ > 0 if necessary; in which ∆ τ is the discrete time step that underlies the computation of the vehicle load tree, and thus defines the horizon of the numerical integration (from one vehicle load to the next) 14 .

Passenger Transfer. For every passenger transfer between two discrete points associated with s (a GP in s or a vehicle docked to s), say from p 1 to p 2 , set up two one-dimensional balance equations acc. to (1); one for the cumulative passenger load of p 1 and the other for the free capacity of p 2 . Ensure that the corresponding passenger flow assigned to the corridor which connects p 1 to p 2 is capacity-and demand-sensitive (cf. Sec. 3.1 on p. 9). Recall that, by assumption, neither of the two balance equations can have any diffusion terms since every modelled uncertainty is captured by a passenger arrival process or the estimations for all initial passenger loads.

Outflow. For every corridor in s that implements a flow of passengers leaving s from one of its GPs, say p, set up a one-dimensional balance equation acc. to (1) (no diffusion term) for the cumulative passenger load of p. Ensure that the specification of the passenger flow assigned to c is demand-sensitive. Now that we have replaced the system of SDEs for the station s in a particular mode by a finite set of decoupled systems of one-and two dimensional SDEs, we have to specify their proper initialization including the passenger loads and the capacity constraints: Imagine that two corridors, say c 1 and c 2 , discharge into the same GP in s, say p; further, that two corridors, say c 3 and c 4 , originate in p; and that all four corridors implement a passenger transfer to/from p from/to two other GPs in s. If we set up four decoupled balance equations for all four passenger transfer possibilities as described above, we then have to i) distribute the initial estimation for p's passenger load among the four balance equations, ii) numerically integrate the four balance equations, and then iii) join the integrated estimations into one single density that captures the passenger load of p at the end of the integration horizon. We achieve this by breaking down p's total capacity into distinct parts; one part, say κ 1 , that defines the total capacity of p in the balance equation set

Fig. 2 .

 2 Fig. 2. Sample station with two capacity-limited gathering points (double circles) and three corridors (double boxes): the entrance area can hold up to 50 passengers and the platform 200 passengers

Fig. 3 .

 3 Fig.3. Every dashed edge specifies the possibility of a passenger transfer from a GP in a station to a vehicle stopped at a waypoint in a transportation grid (through a corridor), or vice versa

Fig. 4 .

 4 Fig. 4. Vehicle mission for the orbital transportation grid from Fig. 1 on p. 3: the single number written in brackets next to every track defines a constant driving time, and the pair of two comma-separated numbers and written in brackets next to every waypoint define minimum and maximum dwell times at that waypoint in the respective order

  and φ m (c, τ ) ≤ φ max (c) for every c ∈ C s , and every τ ≥ 0; in which φ (c, τ ) [i] gives the flow of all passengers who enter c at τ acc. to the trip profile i ∈ Y 5 . It then follows that 6 the product R (c) [j, •] φ (c, τ ) of the j-th row of the routing matrix R (c) assigned to c with φ (c, τ ), gives the flow of all passengers who leave c at τ acc. to the trip profile j ∈ Y. As balance equation for the passenger load M (a, τ ), with M : {a, p} × R ≥0 → R |Y| ≥0 , of a at τ ≥ 0 we thus write down 7 ) dM (a, τ ) := R (c b ) φ (c b , τ ) dτ -φ (c a , τ ) dτ.

  1 , with each corridor in s implementing the possibility of a passenger movement away from/to a GP in s w.r.t. another GP or the outdoor area of s. Once again, we use a Petri net model (with the 'token flow' left out) to capture the connection between all GPs and all corridors in s.

			In
	particular, every vertex in either captures a particular GP (= place) in s (double
	circle) or a particular corridor (= transition; double box); and not a waypoint or a
	track.
	As an example, look at the simple station from Fig. 2.
		Cross
	Step down	turnstiles	Board
	stairs		vehicle
	Entrance area	Platform
	[50]		[200]

We can associate with a corridor a limited throughput for a directed flow of passengers that is conserved between two discrete points; and this metaphor perfectly fits into our gas dynamic approach for all passenger flows

System of differential equations

Y is the finite set of different trip profiles (cf. Sec. 2.3 on p. 6), and |Y| is the cardinality thereof

A continuous-time stochastic process with independent and stationary increments Wt -Ws whose law is Gaussian with parameter t -s

Refer to Tab. 1 on p. 8 for an overview of some matrix/vector operations employed here

From the specification of the passenger routing, which is not specified for this example here

In integral form since this balance equation is supposed to be integrated into (1

Tuning parameter from estimation; in general not time-invariant

Estimation for the passenger loads of all GPs in a particular station and of all vehicles which are docked to this station in our use case

Cartesian product of |Y|-dimensional real vector spaces: for the passenger load vector of every GP in s and of every vehicle docked to s one; in which Y is the set of different trip profiles and |Y| the cardinality thereof; lower-and upper bounded by all non-negativity-and capacity constraints

Cartesian product of a finite set of closed intervals on the real number axis; for every GP in every station and for every vehicle one interval that is lower bounded by zero and upper bounded by the capacity limit of the respective GP or vehicle

This mode transition is unique

Approximation that should asymptotically converge to the exact (analytical) solution by reducing the discretization time step for the discrete mode graph and by increasing the granularity of all numerical integrations

The reflecting boundary condition from[START_REF] Horton | Fluid stochastic petri nets: Theory, applications, and solution techniques[END_REF] still ensures the non-negativity of p's (cumulative) passenger load and p's total capacity

Acknowledgement

This research work has been carried out under the leadership of the Technological Research Institute SystemX, and therefore granted with public funds within the scope of the French Program "Investissements d'Avenir".

up for the passenger transfer associated with c 1 and the remaining part κ 2 that gives the total capacity of p in the balance equation set up for c 2 , respectively. In particular, set κ i := φmax(c i ) φmax(c 1 )+φmax(c 2 ) , where i ∈ {1, 2} and φ max (c i ) is the maximum throughput of c i . Thus, the part of p's total capacity that is attributed to the balance equation set up for c i is proportional to the maximum throughput of c i and inversely proportional to the sum of the throughputs of all corridors that discharge into p. This means that we have to extend the denominator in the definition of κ i if further corridors (implementing passenger transfers or arrival processes) discharge into p. In this context assuming that one such additional corridor implements a passenger arrival process, we might have to further reduce p's total capacity in the set up of the respective balance equation as so as to ensure the maximum throughput as described above. We next use the specification of the passenger (re-)routing to distribute the estimation of p's initial passenger load among the two (decoupled) balance equations set up for the passenger transfers associated with c 3 and c 4 . In fact, we know from the specification of the passenger (re-)routing which group of passengers (which trip profile) wants to take which of the two corridors c 3 or c 4 if at all. We thus extract from the initial estimation the number of all those passengers who want to take c 3 and assign it as the initial (cumulative) passenger load to p in the balance equation set up for c 3 . We do likewise for c 4 . Upon completion of the numerical integration of the balance equations for c 3 and c 4 , we then recover the individual trip profile from the computed density of the cumulative passenger loads by extrapolation; more on this will be described in a paper in preparation.

Summary & Outlook

In this paper we have introduced a continuous-time stochastic hybrid automaton model and showed how its discrete time approximation can be used to forecast passenger loads in multimodal transportation networks: All passengers are grouped into a finite set of trip profiles that route the respective (continuous) passenger flows. All vehicles on the other hand are either parked, or executing missions with deterministic driving times between every pair of two consecutive discrete points (whereas docking times at stations vary according to the passenger load). Uncertainty explicitly enters the model in form of the passengers who start their trips in the network and the estimations for the initial passenger loads of all stations and vehicles.

Two bottlenecks have to be tackled in the iterative computation (from one timelayer to the next) of every forecast: First, there is the vehicle load tree that captures all possible solution paths of the automaton's discrete state that in turn captures the operational state of all vehicles. Every concurrent service run (= vehicle movement with the possibility of passengers on-board) introduces some extra branches which tremendously increase the number of numerical integrations that have to be performed in every time-layer. Now for both bottlenecks we have used two workarounds that we plan to improve on in the future. Moreover, we intend to develop algorithms to efficiently compute forecasts not only for particular passenger loads but also more complex circumstances, such as the travel times between origin-and destination pairs.