Dynamic Modeling for Metro Passenger Flows on Congested Transfer Routes

Abstract

With the rapid development of urbanization, the metro becomes more and more important for people’s travel in big cities. To quantitatively describe metro passenger flows on congested transfer routes, this paper introduces a dynamic model based on automated data from the automatic fare collection (AFC) and automatic vehicle location (AVL) systems. An expectation maximization (EM) algorithm is proposed to compute the maximum likelihood estimates of unknown parameters in our model. Our model can yield a systematic analysis of one-transfer passenger flows on both population and individual aspects. Important characteristics, including transfer time, boarding probabilities, walking time, passenger-to-train assignment probabilities, and total travel time, can be inferred using only the AFC and AVL data. We provide a case study on the Beijing metro. Detailed analysis results based on our model are given. We also present a cross-validation method to validate our model with real data.

1. Introduction

The metro is an important component of the urban transit system. With the advantages of fast speed, punctual operation, and large load, it has become the first choice for most commuters in many big cities. There are over ten million people traveling by metro in Beijing every day [1]. In the meantime, with the rapid development of urbanization, limited numbers of metro lines and trains can no longer meet people’s travel needs. Congestion on platforms and in carriages often occurs, especially in the morning and evening peak periods. To alleviate the problem of congestion, metro administrations should develop reasonable management measures based on effective evaluation of passenger flows. Therefore, it is necessary to study and quantify characteristics of passenger flows, especially during peak hours.

The automatic fare collection (AFC) and automatic vehicle location (AVL) systems in the metro provide a great convenience for collecting data to study passenger flows. We can obtain a passenger’s OD (origin-destination) station, tap-in time, and tap-out time from the AFC data. The route and schedule of each train can be obtained from the AVL data. With the help of AFC data, [2] proposed an estimation model based on a Bayesian inference formulation to enable the study of the main stages affecting passenger travel time, i.e., walking, waiting, transfer, and in-vehicle times. Ref. [3] used AFC and AVL data to propose a comprehensive approach to estimate passenger congestion costs in terms of equivalent travel time loss. Other studies on metro passenger flows based on automated AFC and AVL data can be found in [4,5,6,7,8], among others.

Several authors studied passenger-to-train assignment models and related inferential issues concerning passenger flows on no-transfer metro routes based on automated data [9,10,11]. This paper focuses on congested transfer routes that are more challenging than no-transfer cases. There are limited papers on analysis of metro transfer routes in the literature. Refs. [12,13] developed optimization models for minimizing the transfer waiting time through modifying train schedules. They did not provide methods to infer passenger flows based on automated data. Ref. [14] proposed an equilibrium model based on train operating diagrams and used it to estimate spatial and temporal distributions of the passenger flows at transfer stations. However, this paper did not estimate specific characteristics of passenger flows such as the boarding probability, transfer time, and walking time. Ref. [15] proposed a complete method to estimate the distributions of walking time and waiting time of transfer passengers based on automated data. However, their method is valid only during off-peak periods. Actually, the issue of estimating those distributions becomes much more complicated during peak hours because of the left behind phenomenon. Ref. [16] proposed a data-driven passenger itinerary inference model to infer left behind, route choice, and probabilities of boarding feasible trains. Their model requires a certain assumption on transfer time that is difficult to verify with automated data. Some important characteristics including walking time and total travel time are not provided in [16]. Ref. [11] presented a new passenger-to-train assignment model, and extended it to one-transfer routes. However, they did not provide feasible methods to estimate the input parameters of such a model for transfer cases. Furthermore, Ref. [11] method is based on the data of the same tap-in time, which cannot be straightforwardly applied to dynamic passenger flows.

In this paper we propose a dynamic modeling method for passenger flows on congested transfer routes based on automated data. Our model constructs a probabilistic quantification relationship between a passenger’s tap-in and tap-out times. An expectation maximization (EM) algorithm [17] is proposed to compute the maximum likelihood estimates (MLEs) of unknown parameters in our model. Our model can yield a systematic analysis of one-transfer passenger flows on both population and individual aspects. In particular, inference and uncertainty quantification of transfer time, an important feature on transfer routes, are presented via our model. Other characteristics on a specific one-transfer route including boarding probabilities, walking time, passenger-to-train assignment probabilities, and total travel time can also be inferred. Our method requires only AFC and AVL data on the route under minimal assumptions. We provide a case study on the Beijing metro. Detailed analysis results based on our model are given. Our method can be validated based on real data via cross-validation.

The major contributions of this paper can be summarized as follows.

(i). We provide reasonable statistical modeling and inferential techniques for passenger flows on metro transfer routes, which constitute a comprehensive analysis on passenger flow characteristics, at both the aggregate and individual levels.

(ii). The proposed methods are completely data-driven, with only automated data. We present quantification indices as much as possible by making full use of automated data.

(iii). Based on rigorous mathematical and statistical derivation, we provide a general idea to analyze metro passenger flows. This idea can be extended to other situations.

The rest of this paper is organized as follows. Section 2 provides a literature review. Section 3 introduces the dynamic model for one-transfer routes and provides parameter estimation methods based on AFC and AVL data. Section 4 discusses applications of the proposed model to inferences for passenger flow characteristics. Section 5 provides the case study and validation of our model. Section 6 ends this paper with some conclusions.

2. Research Background

In the literature, research on metro systems is focused on two aspects: the train operation level and the passenger travel level. At the train operation level, the time interval between train departures is an important factor to describe the efficiency of passengers’ travel. A well-organized metro system often requires a reasonable train schedule. Therefore, how to alleviate the metro congestion by adjusting the train schedule has become a hot topic. Ref. [18] proposed a train timetable rescheduling method based on deep learning. This method enables real-time rescheduling of train schedules after random disturbances. Ref. [19] propose a mixed-integer optimization method to generate well-synchronized timetables, thereby reducing the transfer waiting time for transfer passengers. In fact, there are many uncontrollable random factors in metro systems. Optimization of metro system parameters should be based on accurate estimation of the uncertainty caused by these factors, which requires reasonable statistical models. However, such models have not been well developed in the literature, especially for the transfer routes.

At the level of passenger travel, quantitative indices of passenger travel behavior reflect the current situation in the metro system. Statistical or machine learning methods are used to estimate such indices. Ref. [20] developed an ensemble deep learning approach to forecast the dynamic metro passenger flow. This approach was designed for the entire metro network, and cannot accurately estimate the passenger flow indices on a certain route. Ref. [5] analyzed distributional characteristics of the travel time with AFC data, and proposed a method for estimating path-selecting proportions. Ref. [21] provided a passenger-to-train assignment method, and estimated the walking time and waiting time of metro passengers. Although the two papers analyzed the distributional characteristics of the travel time, the impact of passengers’ tap-in times was ignored in their studies. In fact, passengers who enter the station at different tap-in times face different congestion conditions. Hence, to estimate real-time characteristics of the metro system, we need the study on dynamic models for passenger flows.

Transfer is often an indispensable part of long-distance metro travel. However, the problem of estimating dynamic metrics on congested transfer routes has not received a lot of attention in the literature. This paper aims to fill such a gap, and completely data-driven approaches are proposed in the following sections.

3. Dynamic Models for Transfer Routes

3.1. Likelihood Model of Tap-Out Times

We consider a one-transfer trip in metro systems, which consists of two segments. Let segment 1 and segment 2 denote the segments from the origin to the transfer and from the transfer to the destination, respectively. The time between the tap-in and tap-out of a passenger can be divided into the following five periods:

(i).

Boarding time: the time from entering the tap-in fare gate to boarding a train on segment 1;

(ii).

Travel time on the train of segment 1;

(iii).

Transfer time: the time from getting off the train on segment 1 to boarding the train on segment 2;

(iv).

Travel time on the train of segment 2;

(v).

Egress time: the time to walk to the exit fare gate after alighting from the train of segment 2.

The flow chart of such a trip is shown in Figure 1.

For a passenger on the transfer route, denoted by Route A in the following, there may be more than one possible itinerary due to crowded platforms and carriages, especially during peak hours, and it is not determined which train he/she could take with only his/her tap-in and tap-out times. Similar to [11], we need to infer the passenger-to-train assignment and other quantification indices of passenger flows on Route A using the corresponding AVL and AFC data. We then present the following notation and definitions.

Suppose there are n passengers on Route A. Let ${\mathit{T}}_A^{in}={\left(t_{A1}^{in},t_{A2}^{in},\cdots ,t_{An}^{in}\right)}^{\mathrm{T}}$ and ${\mathit{T}}_A^{out}={\left(t_{A1}^{out},t_{A2}^{out},\cdots ,t_{An}^{out}\right)}^{\mathrm{T}}$ denote the sets of their tap-in and tap-out times, respectively. The egress times of the passengers are independently and identically distributed from the log-normal distribution $LN\left({\mu }_A,{\sigma }_A^2\right)$ with probability density function

$$f\left(x|{\mu }_A,{\sigma }_A\right)=exp\{-{\left(lnx-{\mu }_A\right)}^2/2{\sigma }_A^2\}/{\left(2\pi {\sigma }_A^2{x}^2\right)}^{1/2},x>0,$$

(1)

where $-\infty <{\mu }_A<\infty$ and $0<{\sigma }_A<\infty$. This log-normal assumption is widely used and has been verified in the literature (e.g., [15]). Actually, our method in this paper can also be applied to the cases of other egress time distributions such as in [11].

Considering the periodicity of the train schedule, we define a train departure interval as a unit time interval $\left(0,1\right]$, which implies that the departure times of trains at the origin station are $1,2,3,\dots$. We first handle the passengers whose tap-in times lie in $\left(0,1\right]$, i.e., $t_{Ai}^{in}\in \left(0,1\right],i=1,\dots ,n$. General cases for peak hours will be discussed in Section 3.4. Mark the first train passengers can take after arriving at the origin station as Train 1 of segment 1. If a passenger takes Train 1 on segment 1, then the first train he/she can take on segment 2 is recorded as Train 1 on segment 2. Define:

$$\begin{array}{cc}\hfill & D{T}_j^{\left(1\right)}:\mathrm{departure}\mathrm{time}\mathrm{of}\mathrm{Train}j\mathrm{from}\mathrm{the}\mathrm{origin}\mathrm{station}\mathrm{on}\mathrm{segment}1,\hfill \\ \hfill & A{T}_j^{\left(1\right)}:\mathrm{arrival}\mathrm{time}\mathrm{of}\mathrm{Train}j\mathrm{at}\mathrm{the}\mathrm{transfer}\mathrm{station}\mathrm{on}\mathrm{segment}1,\hfill \\ \hfill & D{T}_k^{\left(2\right)}:\mathrm{departure}\mathrm{time}\mathrm{of}\mathrm{Train}k\mathrm{from}\mathrm{the}\mathrm{transfer}\mathrm{station}\mathrm{on}\mathrm{segment}2,\hfill \\ \hfill & A{T}_k^{\left(2\right)}:\mathrm{arrival}\mathrm{time}\mathrm{of}\mathrm{Train}k\mathrm{at}\mathrm{the}\mathrm{terminal}\mathrm{station}\mathrm{on}\mathrm{segment}2,\hfill \\ \hfill & m_i:\mathrm{the}\mathrm{number}\mathrm{of}\mathrm{feasible}\mathrm{itineraries}\mathrm{for}\mathrm{the}\mathrm{passenger}i\mathrm{of}\mathrm{Route}\mathrm{A},\mathrm{then}{m}_i=\#S_i,\hfill \\ \hfill & \mathrm{where}{S}_i=\{(j,k):t_{Ai}^{in}\le D{T}_j^{\left(1\right)},A{T}_j^{\left(1\right)}\le D{T}_k^{\left(2\right)},A{T}_k^{\left(2\right)}\le t_{Ai}^{out}\}\hfill \\ \hfill & =\{(j,k):A{T}_j^{\left(1\right)}\le D{T}_k^{\left(2\right)},A{T}_k^{\left(2\right)}\le t_{Ai}^{out}\},\hfill \\ \hfill & K_i:\mathrm{the}\mathrm{number}\mathrm{of}\mathrm{feasible}\mathrm{itineraries}\mathrm{for}\mathrm{the}\mathrm{passenger}i\mathrm{of}\mathrm{Route}\mathrm{A}\mathrm{on}\mathrm{segment}2,\hfill \\ \hfill & \mathrm{i}.\mathrm{e}.,K_i=\#\left\{k:\mathrm{there}\mathrm{exists}j\mathrm{s}.\mathrm{t}.(j,k)\in S_i\right\},\mathrm{and}\mathrm{let}K=\underset{1⩽i⩽n}{max}{K}_i,\hfill \\ \hfill & m_i^{\left(k\right)}:\mathrm{the}\mathrm{number}\mathrm{of}\mathrm{feasible}\mathrm{itineraries}\mathrm{for}\mathrm{the}\mathrm{passenger}i\mathrm{of}\mathrm{Route}\mathrm{A}\mathrm{on}\mathrm{segment}1\hfill \\ \hfill & \mathrm{who}\mathrm{takes}\mathrm{Train}k\mathrm{on}\mathrm{segment}2,\mathrm{i}.\mathrm{e}.,m_i^{\left(k\right)}=\#\left\{j:(j,k)\in S_i\right\},k=1,\dots ,K_i,\hfill \\ \hfill & \mathrm{where}\#A\mathrm{denotes}\mathrm{the}\mathrm{number}\mathrm{of}\mathrm{elements}\mathrm{in}\mathrm{the}\mathrm{set}A,\mathrm{and}\mathrm{let}{m}^{\left(k\right)}=\underset{1⩽i⩽n}{max}{m}_i^{\left(k\right)},\hfill \\ & J=\underset{1⩽k⩽K}{max}{m}^{\left(k\right)}.\hfill \end{array}$$

We have $m_i={\displaystyle \sum _{k=1}^{K_i}}{m}_i^{\left(k\right)}$ and $1⩽m_i^{\left(1\right)}⩽\cdots ⩽m_i^{\left(K_i\right)}$. For example, if $S_i=\{(1,1),(1,2),(1,3),(2,2),(2,3),(3,3)\}$, we have $m_i=6,K_i=3,m_i^{\left(1\right)}=1,m_i^{\left(2\right)}=2$, and $m_i^{\left(3\right)}=3$.

Let ${\rho }_{jk}\left(t_i^{in}\right)$ denote the probability that a passenger whose tap-in time is $t_i^{in}$ takes Train j on segment 1 and Train k on segment 2 for $j=1,\dots ,J$, $k=1,\dots ,K$. If the passenger takes Train k on segment 2 and $J>m_i^{\left(k\right)}$, then we let ${\rho }_{jk}\left(t_i^{in}\right)=0,j=m_i^{\left(k\right)}+1,\cdots ,J,k=1,\dots ,K_i$. Moreover, if $K>K_i$, then we let ${\rho }_{jk}\left(t_i^{in}\right)=0,j=1,\dots ,J,k=K_i+1,\dots ,K$. Therefore, for each passenger on Route A entering the station at time $t^{in}\in (0,1]$, he/she has at most $JK$ feasible itineraries, and the corresponding probabilities are ${\rho }_{11}\left(t^{in}\right),\cdots ,{\rho }_{1K}\left(t^{in}\right),\dots ,{\rho }_{J1}\left(t^{in}\right),\dots ,{\rho }_{JK}\left(t^{in}\right)$.

Let

$${\eta }_j\left(t^{in}\right)={\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t^{in}\right),j=1,\cdots ,J,$$

(2)

denote the probability that a passenger whose tap-in time is $t^{in}$ takes Train j on segment 1. Let

$${\gamma }_k\left(t^{in}\right)={\displaystyle \sum _{j=1}^J}{\rho }_{jk}\left(t^{in}\right),k=1,\cdots ,K,$$

(3)

denote the probability that a passenger whose tap-in time is $t^{in}$ takes Train k on segment 2. Furthermore, let ${\rho }_{k|j}\left(t^{in}\right)$ denote the conditional probability that a passenger whose tap-in time is $t^{in}$ takes Train k on segment 2 conditionally on the event that he/she takes Train j on segment 1, and let ${\tilde{\rho }}_{j|k}\left(t^{in}\right)$ denote the conditional probability that a passenger whose tap-in time is $t^{in}$ takes Train j on segment 1 conditionally on the event that he/she takes Train k on segment 2. By the conditional probability formula, we have

$$\begin{array}{c}\hfill {\rho }_{k|j}\left(t^{in}\right)={\displaystyle \frac{{\rho }_{jk}\left(t^{in}\right)}{{\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t^{in}\right)}}={\displaystyle \frac{{\rho }_{jk}\left(t^{in}\right)}{{\eta }_j\left(t^{in}\right)}},j=1,\cdots ,J,k=1,\cdots ,K,\end{array}$$

(4)

$$\begin{array}{c}\hfill {\tilde{\rho }}_{j|k}\left(t^{in}\right)={\displaystyle \frac{{\rho }_{jk}\left(t^{in}\right)}{{\displaystyle \sum _{j=1}^J}{\rho }_{jk}\left(t^{in}\right)}}={\displaystyle \frac{{\rho }_{jk}\left(t^{in}\right)}{{\gamma }_k\left(t^{in}\right)}},j=1,\cdots ,J,k=1,\cdots ,K.\end{array}$$

(5)

Let $E_{ik}$ denote the event that passenger i takes Train k on segment 2. Given the tap-in time of passenger i, the distribution function of $t_{Ai}^{out}$ can be written as

$$\begin{array}{cc}\hfill P\left(t_{Ai}^{out}\le t\right)& ={\displaystyle \sum _{k=1}^K}P\left(E_{ik}\right)P\left(t_{Ai}^{out}\le t|E_{ik}\right)\hfill \\ & ={\displaystyle \sum _{k=1}^K}\left({\displaystyle \sum _{j=1}^J}{\rho }_{jk}\left(t_{Ai}^{in}\right)\right)P\left(t_{Ai}^{out}\le t|E_{ik}\right)\hfill \\ & ={\displaystyle \sum _{k=1}^K}\left({\displaystyle \sum _{j=1}^J}{\rho }_{jk}\left(t_{Ai}^{in}\right)\right)P\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}\le t-A{T}_k^{\left(2\right)}|E_{ik}\right)\hfill \\ & ={\displaystyle \sum _{k=1}^K}\left({\displaystyle \sum _{j=1}^J}{\rho }_{jk}\left(t_{Ai}^{in}\right)\right)F_A\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}\right),\hfill \end{array}$$

where $F_A$ is the cumulative distribution of $LN\left({\mu }_A,{\sigma }_A^2\right)$. Taking the derivative of the above equation with respect to $t_{Ai}^{out}$, we have the tap-out time $t_{Ai}^{out}$ which follows a mixture probability distribution with density

$$p_A\left(t_{Ai}^{out}\right)={\displaystyle \sum _{j=1}^J}{\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t_{Ai}^{in}\right)f\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}|{\mathit{\theta }}_A\right),$$

(6)

where ${\mathit{\theta }}_A={\left({\mu }_A,{\sigma }_A\right)}^{\mathrm{T}}$ and f is defined in (1). According to (6), the joint density of ${\mathit{T}}_A^{out}={\left(t_{A1}^{out},t_{A2}^{out},\dots ,t_{An}^{out}\right)}^{\mathrm{T}}$ can be given. Here, we write it as the following likelihood form,

$$\begin{array}{ccc}\hfill L_A\left({\mathit{\phi }}_A\mid {\mathit{T}}_A^{out}\right)& =& {\displaystyle \prod _{i=1}^n}\left\{{\displaystyle \sum _{j=1}^J}{\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t_{Ai}^{in}\right)f\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}|{\mathit{\theta }}_A\right)\right\}\hfill \\ & =& {\displaystyle \prod _{i=1}^n}\left\{{\displaystyle \sum _{k=1}^K}f\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}|{\mathit{\theta }}_A\right)\right\}\left({\displaystyle \sum _{j=1}^J}{\rho }_{jk}\left(t_{Ai}^{in}\right)\right),\hfill \end{array}$$

(7)

where ${\mathit{\phi }}_A={\left({\mathit{\theta }}_A^{\mathrm{T}},{\mathit{\rho }}^{\mathrm{T}}\right)}^{\mathrm{T}},{\mathit{\theta }}_A={\left({\mu }_A,{\sigma }_A\right)}^{\mathrm{T}},\mathrm{and}\mathit{\rho }(·)={\left({\rho }_{11}\left(·\right),\dots ,{\rho }_{JK}\left(·\right)\right)}^{\mathrm{T}}$.

It can be seen from (7) that, by maximizing it over the parameter space, we can only obtain the estimators of ${\displaystyle {\displaystyle \sum _{j=1}^J}}{\rho }_{jk}\left(t\right),k=1,\dots ,K$, i.e., the probabilities of taking the trains on segment 2, while the value of each ${\rho }_{jk}\left(t\right),j=1,\dots ,J,k=1,\dots ,K$ is inestimable. For example, assuming $J=K=2$, we have $\mathit{\rho }\left(t\right)={\left({\rho }_{11}\left(t\right),{\rho }_{12}\left(t\right),{\rho }_{21}\left(t\right),{\rho }_{22}\left(t\right)\right)}^{\mathrm{T}}$. By maximizing (7), we can obtain estimators of ${\rho }_{11}\left(t\right)+{\rho }_{21}\left(t\right)$ and ${\rho }_{12}\left(t\right)+{\rho }_{22}\left(t\right)$, and cannot estimate ${\rho }_{11}\left(t\right),{\rho }_{12}\left(t\right),{\rho }_{21}\left(t\right),{\rho }_{22}\left(t\right)$. To overcome such an identifiability issue, we introduce a method based on auxiliary data, and discuss it in the next subsection.

3.2. Likelihood Model Combined with Auxiliary Data

Let Route B represent the route from the origin to the transfer. Passengers on Route A share segment 1 with those on Route B. Therefore, it can be assumed that they have the same boarding probabilities on segment 1. In this subsection, we introduce the AFC data on Route B as auxiliary data that help to estimate each component of $\mathit{\rho }$ in (7).

Let ${\mathit{T}}_B^{in}={\left(t_{B1}^{in},t_{B2}^{in},\dots ,t_{B{n}_0}^{in}\right)}^{\mathrm{T}}$ and ${\mathit{T}}_B^{out}={\left(t_{B1}^{out},t_{B2}^{out},\dots ,t_{B{n}_0}^{out}\right)}^{\mathrm{T}}$ denote the vectors of passengers’ tap-in and tap-out times on Route B, respectively. Like in Section 3.1, we assume $t_{Bi}^{in}\in (0,1]$ for $i=1,\dots ,n_0$. Assume further that the distribution of egress time on Route B is $LN({\mu }_B,{\sigma }_B^2)$. Given $t_{Bi}^{in}$, the tap-out time $t_{Bi}^{out}$ follows the following mixture probability distribution with density

$$p_B\left(t_{Bi}^{out}\right)={\displaystyle \sum _{j=1}^J}f\left(t_{Bi}^{out}-A{T}_j^{\left(1\right)}|{\mathit{\theta }}_B\right)\left({\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t_{Bi}^{in}\right)\right),$$

where ${\mathit{\theta }}_B={\left({\mu }_B,{\sigma }_B\right)}^{\mathrm{T}}$, and f is defined in (1). The likelihood function based on ${\mathit{T}}^{out}={\left({{\mathit{T}}_A^{out}}^{\mathrm{T}},{{\mathit{T}}_B^{out}}^{\prime }\right)}^{\mathrm{T}}$ is

$$\begin{array}{ccc}\hfill L\left(\mathit{\phi }\mid {\mathit{T}}^{out}\right)& =& {\displaystyle \prod _{i=1}^n}\left\{{\displaystyle \sum _{k=1}^K}f\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}|{\mathit{\theta }}_A\right)\left({\displaystyle \sum _{j=1}^J}{\rho }_{jk}\left(t_{Ai}^{in}\right)\right)\right\}\hfill \\ & & ·{\displaystyle \prod _{i=1}^{n_0}}\left\{{\displaystyle \sum _{j=1}^J}f\left(t_{Bi}^{out}-A{T}_j^{\left(1\right)}|{\mathit{\theta }}_B\right)\left({\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t_{Bi}^{in}\right)\right)\right\},\hfill \end{array}$$

(8)

where $\mathit{\phi }={\left({\mathit{\theta }}_A^{\mathrm{T}},{\mathit{\theta }}_B^{\mathrm{T}},{\mathit{\rho }}^{\prime }\right)}^{\mathrm{T}}$. Maximum likelihood estimators (MLEs) of unknown parameters $\mathit{\phi }$, including each ${\rho }_{jk}\left(t\right),j=1,\dots ,J,k=1,\dots ,K$, can be computed by maximizing (8). The use of ${\mathit{T}}_B^{in}$ and ${\mathit{T}}_B^{out}$ overcomes the identifiability issue with only the data on Route A.

Note that $\mathit{\rho }(·)$ are unknowns functions, and that it is hard to find their MLEs. We use the quadratic logistic function to parameterize them, and assume

$$\begin{array}{cc}& {\rho }_{jk}\left(t|\mathit{\beta }\right)={\displaystyle \frac{exp\left({\mathit{\beta }}_{jk}^{\prime }\tilde{\mathit{t}}\right)}{1+{\displaystyle {\displaystyle \sum _{j=1}^{J-1}}{\displaystyle \sum _{k=1}^K}}exp\left({\mathit{\beta }}_{jk}^{\prime }\tilde{\mathit{t}}\right)+{\displaystyle {\displaystyle \sum _{k=1}^{K-1}}}exp\left({\mathit{\beta }}_{Jk}^{\prime }\tilde{\mathit{t}}\right)}},\hfill \\ & j=1,\dots ,J,k=1,\dots ,K,(j,k)\ne (J,K),\hfill \\ & {\rho }_{JK}\left(t|\mathit{\beta }\right)={\displaystyle \frac{1}{1+{\displaystyle {\displaystyle \sum _{j=1}^{J-1}}{\displaystyle \sum _{k=1}^K}}exp\left({\mathit{\beta }}_{jk}^{\prime }\tilde{\mathit{t}}\right)+{\displaystyle {\displaystyle \sum _{k=1}^{K-1}}}exp\left({\mathit{\beta }}_{Jk}^{\prime }\tilde{\mathit{t}}\right)}},\hfill \end{array}$$

(9)

where $\tilde{\mathit{t}}={\left(1,t,t^2\right)}^{\mathrm{T}},\mathit{\beta }={\left({\mathit{\beta }}_{11}^{\mathrm{T}},\dots ,{\mathit{\beta }}_{J(K-1)}^{\mathrm{T}}\right)}^{\mathrm{T}}$, and ${\mathit{\beta }}_{jk}={\left({\beta }_{0,jk},{\beta }_{1,jk},{\beta }_{2,jk}\right)}^{\mathrm{T}}$ for $j=1,\dots ,J,k=1,\dots ,K,(j,k)\ne (J,K)$. Plugging (9) into (8), we have

$$\begin{array}{ccc}\hfill L\left({\mathit{\phi }}_t\mid {\mathit{T}}^{out}\right)& =& {\displaystyle \prod _{i=1}^n}\left\{{\displaystyle \sum _{k=1}^K}f\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}|{\mathit{\theta }}_A\right)\left({\displaystyle \sum _{j=1}^J}{\rho }_{jk}\left(t_{Ai}^{in}|\mathit{\beta }\right)\right)\right\}\hfill \\ & & ·{\displaystyle \prod _{i=1}^{n_0}}\left\{{\displaystyle \sum _{j=1}^J}f\left(t_{Bi}^{out}-A{T}_j^{\left(1\right)}|{\mathit{\theta }}_B\right)\left({\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t_{Bi}^{in}|\mathit{\beta }\right)\right)\right\},\hfill \end{array}$$

(10)

where ${\mathit{\phi }}_t={\left({\mathit{\theta }}_A^{\mathrm{T}},{\mathit{\theta }}_B^{\mathrm{T}},{\mathit{\beta }}^{\mathrm{T}}\right)}^{\mathrm{T}}$ is a $(3JK+1)$-dimensional vector of unknown parameters.

3.3. The EM Algorithm

This section proposes the EM algorithm [17] to compute MLEs of ${\mathit{\phi }}_t$ in (10). We first introduce two types of hidden variables as follows.

$$\begin{array}{cc}\hfill & U_{i,k}:\mathrm{the}\mathrm{indicator}\mathrm{of}\mathrm{the}\mathrm{event}\mathrm{that}\mathrm{passenger}i\mathrm{on}\mathrm{Route}\mathrm{A}\mathrm{took}\mathrm{Train}k\mathrm{on}\mathrm{segment}2,\hfill \\ \hfill & \mathrm{i}.\mathrm{e}.,U_{i,k}=1\mathrm{if}\mathrm{this}\mathrm{passenger}\mathrm{took}\mathrm{Train}k\mathrm{on}\mathrm{segment}2;U_{i,k}=0\mathrm{otherwise},\hfill \\ \hfill & V_{i,j}:\mathrm{the}\mathrm{indicator}\mathrm{of}\mathrm{the}\mathrm{event}\mathrm{that}\mathrm{passenger}i\mathrm{on}\mathrm{Route}\mathrm{B}\mathrm{took}\mathrm{Train}j,\hfill \\ \hfill & \mathrm{i}.\mathrm{e}.,V_{i,j}=1\mathrm{if}\mathrm{this}\mathrm{passenger}\mathrm{took}\mathrm{Train}j;V_{i,j}=0\mathrm{otherwise}.\hfill \end{array}$$

According to (10), the complete likelihood function can be written as

$$\begin{array}{cc}\hfill L_c\left({\mathit{\phi }}_t\mid {\mathit{T}}^{out},\mathit{U},\mathit{V}\right)& ={\displaystyle \prod _{i=1}^n}{\displaystyle \prod _{k=1}^K}{\left\{f\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}|{\mathit{\theta }}_A\right)\left({\displaystyle \sum _{j=1}^J}{\rho }_{jk}\left(t_{Ai}^{in}|\mathit{\beta }\right)\right)\right\}}^{U_{i,k}}·\hfill \\ \hfill & {\displaystyle \prod _{i=1}^{n_0}}{\displaystyle \prod _{j=1}^J}{\left\{f\left(t_{Bi}^{out}-A{T}_j^{\left(1\right)}|{\mathit{\theta }}_B\right)\left({\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t_{Bi}^{in}|\mathit{\beta }\right)\right)\right\}}^{V_{i,j}},\hfill \end{array}$$

where $\mathit{U}={\left\{U_{i,k}\right\}}_{i=1,\dots ,n,k=1,\dots ,K}$ and $\mathit{V}={\left\{V_{i,j}\right\}}_{i=1,\dots ,n_0,j=1,\dots ,J}$. The corresponding complete log-likelihood function is

$$\begin{array}{cc}\hfill l_c\left({\mathit{\phi }}_t\mid {\mathit{T}}^{out},\mathit{U},\mathit{V}\right)& ={\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{k=1}^K}{U}_{i,k}\left\{log\left(f\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}|{\mathit{\theta }}_A\right)\right)+log\left({\displaystyle \sum _{j=1}^J}{\rho }_{jk}\left(t_{Ai}^{in}|\mathit{\beta }\right)\right)\right\}+\hfill \\ \hfill & {\displaystyle \sum _{i=1}^{n_0}}{\displaystyle \sum _{j=1}^J}{V}_{i,j}\left\{log\left(f\left(t_{Bi}^{out}-A{T}_j^{\left(1\right)}|{\mathit{\theta }}_B\right)\right)+log\left({\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t_{Bi}^{in}|\mathit{\beta }\right)\right)\right\}.\hfill \end{array}$$

The steps of the EM algorithm are as follows.

(I) Given the initial values ${\mathit{\phi }}_t^{\left(0\right)}={\left({{\mathit{\theta }}_A^{\left(0\right)}}^{\mathrm{T}},{{\mathit{\theta }}_B^{\left(0\right)}}^{\mathrm{T}},{{\mathit{\beta }}^{\left(0\right)}}^{\mathrm{T}}\right)}^{\mathrm{T}}$ of the parameters.

(II) Expectation step: after the qth iteration, compute the expectation for the hidden variables $U_{i,k}$ and $V_{i,j}$ in the log-likelihood function, we have

$$\begin{array}{cc}\hfill Q\left({\mathit{\phi }}_t|{\mathit{\phi }}_t^{\left(q\right)}\right)& =E\left(l_c\left({\mathit{\phi }}_t|{\mathit{T}}^{out}\right)|{\mathit{T}}^{out},{\mathit{\phi }}_t^{\left(q\right)}\right)\hfill \\ \hfill & =l_1\left({\mathit{\theta }}_A\right)+l_2\left({\mathit{\theta }}_B\right)+l_3\left(\mathit{\beta }\right),\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill & l_1\left({\mathit{\theta }}_A\right)={\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{k=1}^K}{C}_{ik}^{\left(q\right)}logf\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}|{\mathit{\theta }}_A\right),\hfill \\ \hfill & l_2\left({\mathit{\theta }}_B\right)={\displaystyle \sum _{i=1}^{n_0}}{\displaystyle \sum _{j=1}^J}{D}_{ij}^{\left(q\right)}logf\left(t_{Bi}^{out}-A{T}_j^{\left(1\right)}|{\mathit{\theta }}_B\right),\hfill \\ \hfill & l_3\left(\mathit{\beta }\right)={\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{k=1}^K}{C}_{ik}^{\left(q\right)}log\left({\displaystyle \sum _{j=1}^J}{\rho }_{jk}\left(t_{Ai}^{in}|\mathit{\beta }\right)\right)+{\displaystyle \sum _{i=1}^{n_0}}{\displaystyle \sum _{j=1}^J}{D}_{ij}^{\left(q\right)}log\left({\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t_{Bi}^{in}|\mathit{\beta }\right)\right),\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill & C_{ik}^{\left(q\right)}=E\left(U_{i,k}|{\mathit{\phi }}_t^{\left(q\right)}\right)={\displaystyle \frac{f\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}|{\mathit{\theta }}_A^{\left(q\right)}\right)\left({\displaystyle {\displaystyle \sum _{j=1}^J}}{\rho }_{jk}\left(t_{Ai}^{in}|{\mathit{\beta }}^{\left(q\right)}\right)\right)}{{\displaystyle {\displaystyle \sum _{k=1}^K}}f\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}|{\mathit{\theta }}_A^{\left(q\right)}\right)\left({\displaystyle {\displaystyle \sum _{j=1}^J}}{\rho }_{jk}\left(t_{Ai}^{in}|{\mathit{\beta }}^{\left(q\right)}\right)\right)}},\hfill \\ \hfill & D_{ik}^{\left(q\right)}=E\left(V_{i,j}|{\mathit{\phi }}_t^{\left(q\right)}\right)={\displaystyle \frac{f\left(t_{Bi}^{out}-A{T}_j^{\left(1\right)}|{\mathit{\theta }}_B^{\left(q\right)}\right)\left({\displaystyle {\displaystyle \sum _{k=1}^K}}{\rho }_{jk}\left(t_{Bi}^{in}|{\mathit{\beta }}^{\left(q\right)}\right)\right)}{{\displaystyle {\displaystyle \sum _{j=1}^J}}f\left(t_{Bi}^{out}-A{T}_j^{\left(1\right)}|{\mathit{\theta }}_B^{\left(q\right)}\right)\left({\displaystyle {\displaystyle \sum _{k=1}^K}}{\rho }_{jk}\left(t_{Bi}^{in}|{\mathit{\beta }}^{\left(q\right)}\right)\right)}}.\hfill \end{array}$$

(III) Maximization step: maximizing $Q\left({\mathit{\phi }}_t|{\mathit{\phi }}_t^{\left(q\right)}\right)$ yields the $(q+1)$th iteration of parameter ${\mathit{\phi }}_t$, where ${\mathit{\theta }}_A^{(q+1)}=argmax{l}_1\left({\mathit{\theta }}_A\right)$, ${\mathit{\theta }}_B^{(q+1)}=argmax{l}_2\left({\mathit{\theta }}_B\right)$, and ${\mathit{\beta }}^{(q+1)}=argmax{l}_3\left(\mathit{\beta }\right)$.

(IV) Iteration: repeat steps (II) and (III) until convergence.

It can be seen that the EM algorithm simplified the original optimization problem through introducing hidden variables. Note that its performance relies on the choice of the initial values. Multiple initial values can be used to better approximate the global solution. Our experience based on numerical simulations shows that the proposed EM algorithm is efficient to compute the MLEs.

3.4. Periodic Modeling

In this paper, we focus on modeling for the transfer routes during peak hours. With train departure intervals of equal lengths, the passenger flow has stable performance. Therefore, the distribution f of egress times remains, and the boarding probabilities $\mathit{\rho }\left(t\right)$ behave periodically. Specifically, consider train intervals $(D{T}_j^{\left(1\right)},D{T}_{j+1}^{\left(1\right)}],\dots ,(D{T}_{j+l-1}^{\left(1\right)},D{T}_{j+l}^{\left(1\right)}]$ in such a period. The functions $\mathit{\rho }\left(t\right)$ in any two of these intervals are identical. We construct a periodic model for $\mathit{\rho }\left(t\right)$ to analyze the AFC data of passengers who enter the metro system during peak hours.

Without loss of generality, let $j=1$, $D{T}_k^{\left(1\right)}=k-1$ for $k=1,\dots ,l+1$. For $t\in (D{T}_1^{\left(1\right)},D{T}_{l+1}^{\left(1\right)}]=(0,l]$, assume that

$$\begin{array}{cc}& {\rho }_{jk}\left(t|\mathit{\beta }\right)={\displaystyle \frac{exp\left({\mathit{\beta }}_{jk}^{\prime }\left\{\dot{\mathit{t}}\right\}\right)}{1+{\displaystyle \sum _{j=1}^{J-1}}{\displaystyle \sum _{k=1}^K}exp\left({\mathit{\beta }}_{jk}^{\prime }\left\{\dot{\mathit{t}}\right\}\right)+{\displaystyle \sum _{k=1}^{K-1}}exp\left({\mathit{\beta }}_{Jk}^{\prime }\left\{\dot{\mathit{t}}\right\}\right)}},\hfill \\ & j=1,\dots ,J,k=1,\dots ,K,(j,k)\ne (J,K),\hfill \\ & {\rho }_{JK}\left(t|\mathit{\beta }\right)={\displaystyle \frac{1}{1+{\displaystyle \sum _{j=1}^{J-1}}{\displaystyle \sum _{k=1}^K}exp\left({\mathit{\beta }}_{jk}^{\prime }\left\{\dot{\mathit{t}}\right\}\right)+{\displaystyle \sum _{k=1}^{K-1}}exp\left({\mathit{\beta }}_{Jk}^{\prime }\left\{\dot{\mathit{t}}\right\}\right)}},\hfill \end{array}$$

(11)

where $\left\{\dot{\mathit{t}}\right\}={\left(1,\dot{t},{\dot{t}}^2\right)}^{\mathrm{T}}$ and $\dot{t}=t-\lfloor t\rfloor$ if $t-\lfloor t\rfloor >0$, $\dot{t}=1$ otherwise. $\lfloor ·\rfloor$ denotes the floor function. $\mathit{\beta }={\left({\mathit{\beta }}_{11}^{\mathrm{T}},\dots ,{\mathit{\beta }}_{J(K-1)}^{\mathrm{T}}\right)}^{\mathrm{T}}$, and ${\mathit{\beta }}_{jk}={\left({\beta }_{0,jk},{\beta }_{1,jk},{\beta }_{2,jk}\right)}^{\mathrm{T}}$ for $j=1,\dots ,J,k=1,\dots ,K,(j,k)\ne (J,K)$ are unknown parameters.

Clearly, (9) is a special case of (11) with $l=1$. Replacing (9) with (11) in (10), we obtain the periodicity model of $({\mathit{T}}_A^{in},{\mathit{T}}_A^{out},{\mathit{T}}_B^{in},{\mathit{T}}_B^{out})$, where $t_{Ai}^{in}\in \left(0,l\right],i=1,\cdots ,n$ and $t_{Bi}^{in}\in \left(0,l\right],i=1,\cdots ,n_0$. The method to estimate the unknown parameters for $l=1$ in previous subsections can also be used here.

Figure 2 is a process flow chart of the proposed methodology. The first part presents the model and the method for estimating an unknown parameter in the model. The second part provides applications of the model and estimated parameters, which will be discussed in the next section.

4. Applications of Analysis of Transfer Routes

In this section, we present some applications of our dynamic model in the previous sections for the analysis of transfer routes during peak hours. First of all, we should point out that the parameters in (8) can be viewed as measures that describe some aspects of characteristics of transfer passenger flows. Therefore, our estimators can directly be used to infer these measures. Furthermore, our method can be used to infer the transfer time and to predict the travel time on the transfer route.

4.1. Inference for the Transfer Time

Recall that a passenger’s transfer time is the time from getting off the train on segment 1 to boarding the train on segment 2. It contains transfer walking time and waiting time at the platform of the transfer station, and thus is an important index to quantify the congestion condition of the transfer channel and segment 2. Based on AFC and AVL data, this subsection studies inferences for transfer times via the proposed dynamic model.

First, we consider the population aspect, i.e., inference for the transfer time when only the tap-in time is given. By definition, the transfer time $T^{transfer}$ of a passenger takes its value on a finite set $\{D{T}_k^{\left(2\right)}-A{T}_j^{\left(1\right)},j=1,\dots ,J,k=1,\dots ,K\}$. Given the passenger’s tap-in time $t_A^{in}$, it is not hard to derive the distribution of the transfer time

$$P\left(T^{transfer}=D{T}_k^{\left(2\right)}-A{T}_j^{\left(1\right)}|t_A^{in}\right)={\rho }_{jk}\left(t_A^{in}\right)j=1,\dots ,J,k=1,\dots ,K.$$

Therefore, the conditional mean and conditional standard deviation of $T^{transfer}$ are given by

$$\begin{array}{cc}\hfill & E\left(T^{transfer}|t_A^{in}\right)={\displaystyle \sum _{j=1}^J}{\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t_A^{in}\right)\left(D{T}_k^{\left(2\right)}-A{T}_j^{\left(1\right)}\right),\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & sd\left(T^{transfer}|t_A^{in}\right)={\left\{{\displaystyle \sum _{j=1}^J}{\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t_A^{in}\right){\left(D{T}_k^{\left(2\right)}-A{T}_j^{\left(1\right)}\right)}^2-{\left(E\left(T^{transfer}|t_A^{in}\right)\right)}^2\right\}}^{1/2}.\hfill \end{array}$$

(13)

The MLEs of ${\rho }_{jk}\left(t_A^{in}\right)$ obtained in Section 3 can be used to compute these quantities.

Next, we consider the individual aspect, i.e., inference for the transfer time when the tap-in and tap-out times are given. To study this problem, we need the modified passenger-to-train assignment model [11],

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p_{ijk}& =P\left(\mathrm{passenger}i\mathrm{takes}\mathrm{Train}j\mathrm{on}\mathrm{segment}1\mathrm{and}\mathrm{Train}k\mathrm{on}\mathrm{segment}2|t_{Ai}^{in},t_{Ai}^{out}\right)\hfill \\ \hfill & ={\displaystyle \frac{f\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}\right){\rho }_{jk}\left(t_{Ai}^{in}\right)}{{\displaystyle \sum _{j=1}^J}{\displaystyle \sum _{k=1}^K}f\left(t_{Ai}^{out}-A{T}_k^{\left(2\right)}\right){\rho }_{jk}\left(t_{Ai}^{in}\right)}}.\hfill \end{array}\end{array}$$

(14)

Consequently, the conditional probability distribution, conditional mean, and conditional standard deviation of $T^{transfer}$ are given by

$$\begin{array}{ccc}& & P_i\left(T^{transfer}=D{T}_k^{\left(2\right)}-A{T}_j^{\left(1\right)}|t_{Ai}^{in},t_{Ai}^{out}\right)=p_{ijk}j=1,\dots ,J,k=1,\dots ,K,\hfill \\ & & E_i\left(T^{transfer}|t_{Ai}^{in},t_{Ai}^{out}\right)={\displaystyle \sum _{j=1}^J}{\displaystyle \sum _{k=1}^K}{p}_{ijk}\left(D{T}_k^{\left(2\right)}-A{T}_j^{\left(1\right)}\right),\hfill \\ & & s{d}_i\left(T^{transfer}|t_{Ai}^{in},t_{Ai}^{out}\right)={\left({\displaystyle \sum _{j=1}^J}{\displaystyle \sum _{k=1}^K}{p}_{ijk}{\left(D{T}_k^{\left(2\right)}-A{T}_j^{\left(1\right)}\right)}^2-{\left(E_i\left(T^{transfer}|t_{Ai}^{in},t_{Ai}^{out}\right)\right)}^2\right)}^{1/2}.\hfill \end{array}$$

(15)

These quantities can be calculated by (14) and the MLEs of f and ${\rho }_{jk}$. The passenger-to-train assignment model and individual inference for $T^{transfer}$ can be used to infer the movement of a specified passenger in the metro system according to his/her tap-in and tap-out times.

4.2. Prediction and Validation

This subsection discusses how to predict the tap-out time of a transfer passenger when his/her tap-in time is given. Hence, the total travel time can be predicted given one’s tap-in time. The prediction method can also be used to validate our method via cross-validation (CV).

By (6), given $t_A^{in}$, the conditional density function of the tap-out time $t_A^{out}$ of the passenger is

$$p_A\left(t_A^{out}\right)={\displaystyle \sum _{j=1}^J}{\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t_A^{in}\right)f\left(t_A^{out}-A{T}_k^{\left(2\right)}|{\mathit{\theta }}_A\right).$$

(16)

When the parameters are known, the predictor of $t_A^{out}$ that minimizes the mean squared prediction error is the following conditional mean:

$$E\left(t_A^{out}\right|t_A^{in})={\displaystyle \sum _{j=1}^J}{\displaystyle \sum _{k=1}^K}{\rho }_{jk}\left(t_A^{in}\right)\int tf\left(t-A{T}_k^{\left(2\right)}|{\mathit{\theta }}_A\right)dt.$$

(17)

Plugging the MLEs of the parameters into (17), we get the empirical version of the best predictor. With the MLEs, the conditional distribution of $t_A^{out}$ can be estimated via (16), and this can yield prediction intervals of $t_A^{out}$. Specifically, given a confidence level $\alpha \in (0,1)$, the $100\left(1-\alpha \right)\%$ prediction interval of the tap-out time of a passenger with the tap-in time $t_A^{in}$ can be given by

$$[c\left(\alpha /2\right),c\left(1-\alpha /2\right)],$$

(18)

where $c\left(\alpha \right)$ represents the lower $\alpha$-quantile of the estimator of the distribution (16).

In addition, the accuracy of the above prediction method on a test sample can be used to validate the proposed model. When there does not exist a test sample, the CV method [22] can be applied. We randomly divide the current AFC data of Route A into two parts. One is used as the training sample to train the model parameters by the method in Section 3, and the other part acts as the test sample that is used to evaluate the prediction accuracy.

5. Case Study

5.1. Data Description

There are 24 metro lines of the total mileage 727 km in Beijing at the end of 2020. With an increase of 499 km between 2010 and 2020, the Beijing metro has entered an era of rapid growth [23]. The congestion problem broadly exists in the Beijing metro during peak hours. Analysis of metro data is necessary to study passenger flows on congested routes, which is helpful for metro management.

In this section, the proposed dynamic model is applied to analyze the passenger flow on the following route (Route A) in the Beijing metro during peak hours:

• Segment 1: from Wuzixueyuanlu Station to Hujialou Station (Line 6, 8 stops),
• Segment 2: from Hujialou Station to Liangmaqiao Station (Line 10, 3 stops),

as shown in Figure 3.

The data used to analyze Route A are the AFC and AVL data from 7:00 a.m. to 9:00 a.m. in ten weekdays in 2021. The auxiliary data are from:

• Route B. from Wuzixueyuanlu Station to Hujialou Station (Line 6, 8 stops),

during the same period. The sample sizes of Route A and Route B are 2494 and 4084, respectively. Figure 4 shows the frequency histogram of Route A passengers’ tap-in and tap-out times during peak hours.

According to the data, the departure interval between two trains at Wuzixueyuanlu is about 3 min. The train travel times of segments 1 and 2 are about 21 and 6 min, respectively. The train whose departure time is 7:03 is set as Train 1 of segment 1. The first train departure interval is $\left(7\text{:}00,7\text{:}03\right]$. After removing several large outliers, we find that each passenger who enters in the interval taps out before 7:41. Therefore, all possible trains for these passengers are those shown in

Table 1. Train schedule.

	Wuzixueyuanlu Station		Hujialou Station
Train	Departure Time	Arrival Time	Departure Time	Arrival Time
1	7:03:00	7:24:00	7:26:24	7:32:24
2	7:06:00	7:27:00	7:29:24	7:35:24
3	7:09:00	7:30:00	7:32:24	7:38:24

, which implies $J=K=3$. We can further infer that the set of all feasible itineraries is

$$S=\left\{\right(1,1),(1,2),(1,3),(2,2),(2,3),(3,3\left)\right\}$$

(19)

from this table. By periodicity, this holds for all passengers whose tap-in times are from 7:00 to 9:00.

5.2. Parameter Estimation

Based on the AFC and AVL data on Route A and Route B, we build the likelihood model and use the proposed EM algorithm to compute MLEs of unknown parameters. The MLEs are ${({\widehat{\mu }}_A,{\widehat{\sigma }}_A)}^{\mathrm{T}}={(0.2329,0.5620)}^{\mathrm{T}}$ and ${({\widehat{\mu }}_B,{\widehat{\sigma }}_B)}^{\mathrm{T}}={(0.3003,0.5063)}^{\mathrm{T}}$. From this, we can get the mean and standard deviation of the distribution of egress time on Route A are $4.4346$ and $0.9008$, respectively; those for Route B are $4.6047$ and $0.8296$, respectively. Figure 5 shows the corresponding probability density functions. There is not much difference between the egress times on the two routes.

In the meantime the EM algorithm yields the MLEs of $\mathit{\beta }$ in the logistic model (11), which further give the probabilities ${\rho }_{jk}\left(t^{in}\right),(j,k)\in S$ in (19). Their curves with $t^{in}$ varying in the first train departure interval are shown in Figure 6. It can be seen from Figure 6 that passengers are more likely to take train combinations $(1,1)$ and $(1,2)$ when they tap in the station earlier, and that they are more likely to take train combinations $(2,2)$ and $(3,3)$ when their tap-in times are closer to the train departure time. Overall, the later a passenger’s tap-in time is, the more likely he/she takes the later trains on segment 1 and segment 2. For an in-depth study, we also need to discuss the boarding probabilities and conditional boarding probabilities on segment 1 and segment 2. Curves of related boarding probabilities in (2)–(5) are presented in Figure 7. Some interesting phenomena can be found from the two figures. It can be seen from (a) and (b) in Figure 7 that, as expected, the probabilities of boarding the first (third) train on both segments 1 and 2 become smaller (larger) as $t^{in}$ increases. From Figure 6, or (a) in Figure 7, we can find that passengers who tap in after 7:02 have very small probabilities of boarding the first train on segment 1. This may make the values of the conditional probabilities ${\rho }_{k|1}\left(t^{in}\right)$ shown by (c) in Figure 7 inaccurate. From (d) in Figure 7, the conditional probabilities ${\rho }_{2|2}\left(t^{in}\right)$ have an increasing trend as $t^{in}$ goes to 7:03. A possible reason is that a passenger who enters late but boards the second train has a fast walking speed, and thus he/she has larger probability of boarding his/her first feasible train (Train 2) on segment 2. From (e) and (f) in Figure 7, we can see that the trajectories of the conditional probabilities ${\tilde{\rho }}_{j|k}\left(t^{in}\right)$ are consistent with our expectation. Specifically, given the train which passengers take on segment 2, we can deduce that passengers with an earlier tap-in time have a higher probability of taking Train 1 on segment 1; passengers with a later tap-in time have a higher probability of taking Train 2 or Train 3 on segment 1.

Plugging the MLEs into (6), we can get the conditional density functions of $t^{out}$ conditionally on $t^{in}$. Figure 8 presents these function curves given $t^{in}=$ 7:00:02, 7:01:00, 7:02:00, and 7:03:00, respectively. They are all mixture density functions, and have more than one peak. According to (17), we can also derive the conditional mean of $t^{out}$, shown in Figure 9. As expected, it is an increasing function with a non-linear pattern. By periodicity, it can be extended to the situation when $t^{in}$ lies between 7:00 and 9:00. This curve can be used to predict the tap-out time of a passenger given his/her tap-in time.

5.3. Inference for $T^{transfer}$

In this subsection, we infer the transfer times of passengers on Route A. First, consider the population aspect. According to Table 1, all possible values of $T^{transfer}$ are in $\{2.4min,5.4min,8.4min\}$. Let $T_1^{*}=2.4min$, $T_2^{*}=5.4min$, and $T_3^{*}=8.4min$. We have $P(T^{transfer}=T_1^{*}|t^{in})={\rho }_{11}\left(t^{in}\right)+{\rho }_{22}\left(t^{in}\right)+{\rho }_{33}\left(t^{in}\right)$, $P(T^{transfer}=T_2^{*}|t^{in})={\rho }_{12}\left(t^{in}\right)+{\rho }_{23}\left(t^{in}\right)$, and $P(T^{transfer}=T_3^{*}|t^{in})={\rho }_{13}\left(t^{in}\right)$. Figure 10 shows the curves of $P(T^{transfer}=T_i^{*}|t^{in}),i=1,2,3$ with respect to the tap-in time $t^{in}$. It can be seen that, when $t^{in}$ is close to 7:03, $P(T^{transfer}=T_1^{*}|t^{in})$ has an increasing trend. A possible reason is that, some passengers in a hurry picked up speed when seeing the train coming, and when they missed the train, they were more likely to walk fast through the transfer channel. By (12) and (13), we show the conditional mean and standard deviation of $T^{transfer}$ in Figure 11. Figure 11 also indicates that the passengers whose tap-in times are close to 7:03 have a shorter $T^{transfer}$ on average, which is consistent with the above inferring result. Specifically, passengers who tap in the origin station earlier are more likely to be on Train 1 of segment 1. They have less time to wait and are therefore less likely to be in a hurry to transfer, resulting in a longer transfer time. Conversely, passengers who tap in the origin station later will usually spend some time waiting for a train at the origin station, so they will have to walk faster if they want to catch Train 2 in segment 2, thus reducing their transfer times.

Next, consider the individual aspect. By (14), we compute the six passenger-to-train assignment probabilities $p_{ijk},(j,k)\in S$ for given $t_i^{in}$ and $t_i^{out}$, and show their contour maps in Figure 12. Moreover, by (15), we present the contour map of the average transfer time for given $t_i^{in}$ and $t_i^{out}$ in Figure 13. From these figures, we can infer the behavior of an individual passenger whose tap-in and tap-out times are known on the transfer route. For example, assuming that passenger i has a tap-in time of 7:02:30 and a tap-out time of 7:42:00, we can draw the conclusion of $p_{i11}<p_{i12}<p_{i22}<p_{i13}<p_{i33}<p_{i23}$ from Figure 12.

5.4. Validation via CV

We use the CV method to validate the proposed model. The AFC data of passengers on Route A are randomly divided into a training set and a test set with proportions of 80% and 20%, respectively. By (18), we use the training data to construct the 95% prediction interval of the tap-out time for each tap-in time during the peak hours; see Figure 14. It can be seen that about 95% of the test data lie within the corresponding prediction intervals. This provides a validation of our model. Furthermore, an interesting phenomenon can be found from Figure 14: when a passenger’s tap-in time is close to 7:00 or 9:00, the time when the peak hours begin or end, his/her tap-out time becomes closer to the lower limit of the prediction interval. The reason is that, such a passenger has a relatively short travel time since the congestion status is not so serious at that time. We can also see that the passengers whose tap-in times lie between 7:30 and 8:30 have more even tap-out times between the upper and lower limits of the prediction interval.

6. Conclusions

This paper proposes a dynamic logistic model to study passenger flow characteristics on one-transfer metro routes. Estimates of parameters in our model can be computed by the EM algorithm based on automated AFC and AVL data. A number of passenger flow characteristics including transfer time, boarding probabilities, and travel time can be inferred and quantified with the proposed model. We also present a CV method to validate our model. Real applications to the Beijing metro show the effectiveness of the proposed model.

The proposed model and its derived inferential methods possess the following advantages. First, they are completely data-driven, with only automated data. We do not require additional information including survey data, speed information of passengers, and layouts of stations. Second, our methods provide a comprehensive analysis on passenger flow characteristics. Using reasonable statistical modeling and inferential techniques, we present quantification indices as much as possible by making full use of automated data. The obtained results can be used for the optimal management of passenger flows. For instance, we can estimate the transfer times on many transfer routes, and consider adjusting the train schedule or implementing passenger limiting at the transfer stations with long transfer times. Third, the basic idea of the dynamic logistic model is to derive the conditional distribution of a passenger’s tap-out time conditionally on his/her tap-in time. This is a very general idea, and can be followed in more complex situations.

There are several future directions:

(i)

Extensions of the proposed method to the situations where there are multiple transfers and/or multiple route selections between the origin and destination;

(ii)

Generalizations of the proposed method, possibly combined with deep learning or other machine learning methods; for analyzing the whole metro network;

(iii)

Combinations with optimization methods for improving metro management.