Probabilistic Model of Delay Propagation along the Train Flow

Abstract

In this chapter, we propose a probabilistic model for train delay propagation. There are deduced formulas for the probability distributions of arrival headways and knock-on delays depending on distributions of the primary delay duration and the departure headways. We prove some key mathematical statements. The obtained formulas allow to predict the frequency of train arrival delays and to determine the optimal traffic adjustments. Several important special cases of initial probability distributions are considered. Results of the theoretical analysis are verified by comparison with statistical data on the train traffic at the Russian railways.

1. Introduction

The trains’ movement is subject to a variety of random factors which leads to unplanned delays. This causes the scattering of the arrival times, hence, the inconvenience to passengers and consignees. Knowledge of the arrival times’ distribution properties leads to the possibility of predicting the characteristics of the train traffic and making correct decisions on the transportation process management. This makes it possible to improve the punctuality of train traffic and save resources, in particular, electric power.

The properties of the arrival headways distributions allow us to estimate the probability of delays emergence and theirs characteristics, which are important from a practical point of view. Probabilistic modeling of the delay propagation process along the train flow is the main tool for solving this problem.

The models for the distribution of delays in a dense train flow are divided into two classes. These are deterministic and stochastic models. Stochastic models take into account the unpredictable nature of obstacles in the railway. A mathematical model, proposed in the present chapter, make it possible to determine the probability distributions of the arrival headways of two consecutive trains at the station. The distribution properties are analyzed for different scattering of input random variables (the primary delay and the initial headways). Comparison of theoretical distributions with real statistics of train traffic on the Russian railways is performed.

2. Literature review

A substantial volume of literature is devoted to study of the train delays effect on the railway functioning. Deterministic models for primary and knock-on delays description were proposed in [1, 2]. These models based on the application of graph theory allow adjust the train traffic schedule. However, such approach considering the different characteristics of train traffic (e.g., travel and dwell times, headways, etc.) as deterministic values does not take into account the uncertainties that arise in reality.

Stochastic modeling takes the influence of random factors (e.g., see [3, 4, 5, 6, 7, 8]) into account. Authors of [7] determine a probabilistic distribution of the arrival times. The problem of finding a distribution of arrival train delays is examined in [8]. It should be noted that in these papers, special cases of primary delay distribution are considered. It is supposed in [8] that the random duration of the primary delay corresponds to some generalization of the exponential law. The paper [7] employs discretization of the delay distribution.

Some of the researchers have analyzed statistical data on deviations of the train arrival times from the planned ones. In particular, the papers [9, 10, 11] show that scattering of these deviations correspond to the exponential distribution.

3. Description of models and analysis of the arrival headways distribution

3.1. The first model

Trains follow one path one after another in one direction from station A to station B with the same average speed $v_0$ . Let the total number of trains is n. The distance from the train j to the train (j − 1) is denoted by $X_j+s_0$ , where j = 2, 3, …, n, $s_0>0$ is the minimal safe distance between trains, and $X_2$ , $X_3$ , …, $X_n$ are the random variables (without any assumptions about their distributions). All trains have the same destination station.

Let us also introduce the notations: ${\mu }_j=X_j/v_0$ , $t_0=s_0/v_0$ . Suppose that train 1 departs from station A at the time $t=0$ . Then, the moment $T^{\left(m\right)}$ of departure train m can be found as (as shown at Figure 1):

$$T^{\left(m\right)}=\sum _{j=2}^m{\mu }_j+\left(m-1\right)t_0, m=2,3,\dots ,n$$ E1

Figure 1.

Departure times of trains 1, 2, and 3 from station A.

Assume that at some point in time, train 1 makes unplanned stop. The duration of this stop is random value $\tau$ . The subsequent trains suffer knock-on delays, when the value $\tau$ is large enough. Following train stops when the distance to the front train is reduced to $s_0$ . It is assumed that as soon as the front train restore running, then the next one immediately follows it. The following problem is considered: to find out the probability distribution of the random arrival headway between the trains (k − 1) and k at the destination B (denote this headway as ${\nu }_k$ ), assume that only the first train makes an unplanned stop. In other words, we need to find the (cumulative) distribution functions $W_k\left(t\right)=\mathrm{{\rm P}}\left({\nu }_k<t\right)$ , k = 2, 3, …, n. Call this problem by the first problem.

3.2. The second model

Suppose that train 1 was delayed at station A at the moment $t=0$ and waited for a random time $\tau$ . If $\tau <{\mu }_2$ , then trains 2, 3, and so on, depart at the planned times: $T^{\left(2\right)}$ , $T^{\left(3\right)}$ , etc. If $\tau >{\mu }_2$ , then train 2 will be delayed and will depart at the time $\tau +t_0>T^{\left(2\right)}.$ Train 3 departs according to the same rule depending on the delay time of train 2, and so on. In this formulation, ${\nu }_k$ is actual departure headway between the trains with numbers (k − 1) and k. It is required to determine the distribution functions $W_k\left(t\right)$ of random variables ${\nu }_k$ , k = 2, 3, …, n.

Example 1. Let n = 5, ${\mu }_k=2$ , $k=\overline{2,5}$ , $t_0=1$ . The moments of planned departures of trains satisfy the equalities $T^{\left(k\right)}=3\left(k-1\right)$ , $k=\overline{1,5}$ . Figure 2 shows the process of headways ${\nu }_k$ forming, $k=\overline{2,5},$ depending on the six values of the interval $\tau$ . The dots represent real train departure times that result from the primary delay $\tau$ .

Figure 2.

The headways ${\nu }_k$ for some values $\tau$ .

The basic model assumptions are follows: (1) only train 1 is exposed to primary delay $\tau$ . (2) $T^{\left(k\right)}-T^{\left(k-1\right)}>t_0$ , k = 2, 3, …, n.

Denote by $R^{\left(k\right)}$ the real departure time of the train with number k, which depends on $\tau$ and $t_0$ .

We suppose that the departure times of trains satisfy the following two rules. Let k be fixed, $2\le k\le n$ . The first rule: if $R^{\left(k-1\right)}\le T^{\left(k\right)}-t_0$ , then $R^{\left(k\right)}=T^{\left(k\right)}$ . The second rule: if $R^{\left(k-1\right)}\ge T^{\left(k\right)}-t_0$ , then $R^{\left(k\right)}=R^{\left(k-1\right)}+t_0$ . Obviously, $R^{\left(k\right)}\ge T^{\left(k\right)}$ .

In what follows, we use the notation $I\left(x\in A\right)=\left\{\begin{array}{l}1,\mathit{\text{if}}x\in A,\\ 0,\mathit{\text{if}}x\in R\backslash A,\end{array}\right.$ where A is an arbitrary set on the real line R.

Suppose that the total number of trains is equal to $n\ge 2$ . Formally, we set ${\nu }_k=0$ if $k>n$ . Let us proceed to the formulation of the obtained results. We note that the proofs of the majority of the assertions are not given here due to the condition on the size. They take up a lot of space and will be published in our other work.

Theorem 1. 1. If $\tau <{\mu }_2$ , then ${\nu }_2={\mu }_2+t_0-\tau$ , ${\nu }_k={\mu }_k+t_0$ , $3\le k\le n$ .

2. Let k be a fixed integer, $2\le k\le n$ . If $\tau \ge {\sum }_{j=2}^k{\mu }_j$ , then ${\nu }_2=\dots ={\nu }_k=t_0$ .

3. If ${\sum }_{j=2}^k{\mu }_j\le \tau <{\sum }_{j=2}^{k+1}{\mu }_j$ , then

$${\nu }_{k+1}=I\left(k+1\le n\right)\left[\sum _{j=2}^{k+1}{\mu }_j+t_0-\tau \right],$$ E2

$${\nu }_m=I\left(m\le n\right)\left({\mu }_m+t_0\right), m=k+2,\dots ,n$$ E3

Theorem 2. Let $n\ge 2$ . For any k, $2\le k\le n$ , the following formula holds

$$W_k\left(t\right)=I\left(t>t_0\right)\left[\mathrm{{\rm P}}\left({\mu }_k<t-t_0,\tau <\sum _{j=2}^{k-1}{\mu }_j\right)+\mathrm{{\rm P}}\left(\tau +t-t_0>\sum _{j=2}^k{\mu }_j,\tau \ge \sum _{j=2}^{k-1}{\mu }_j\right)\right],$$ E4

in particular,

$$W_2\left(t\right)=I\left(t>t_0\right)\mathrm{{\rm P}}\left(\tau +t-t_0>{\mu }_2\right)$$ E5

Let us introduce the notations, $G\left(x\right)=\mathrm{{\rm P}}\left(\tau <x\right)$ , $\overline{G}\left(x\right)=\mathrm{{\rm P}}\left(\tau >x\right)$ . Note that $G\left(x\right)+\overline{G}\left(x\right)+\mathrm{{\rm P}}\left(\tau =x\right)=1$ . We denote by $g\left(x\right)$ the density function of $\tau$ in the case when it is absolutely continuous.

Further, some corollaries of Theorem 2 are formulated.

Corollary 1. Let ${\mu }_j$ , $2\le j\le n$ , be arbitrary positive numbers, then for $2\le k\le n$

$$W_k\left(t\right)=I\left(t_0<t\le {\mu }_k+t_0\right)\overline{G}\left(\sum _{j=2}^k{\mu }_j-t+t_0\right)+I\left(t>{\mu }_k+t_0\right),$$ E6

in particular,

$$W_2\left(t\right)=I\left(t>t_0\right)\overline{G}\left({\mu }_2-t+t_0\right).$$ E7

Example 2. Let the primary delay $\tau$ have exponential distribution, that is,

$$g\left(t\right)=I\left(t\ge 0\right)\lambda e^{-\mathit{\lambda t}}, \lambda >0$$ E8

As initial parameters, we take the following quantities.

$$\lambda =0.4, t_0=3, {\mu }_2=5, {\mu }_3=6, {\mu }_4=10$$ E9

Graphs of the functions $W_2\left(t\right)$ from Eq. (7), $W_3\left(t\right)$ and $W_4\left(t\right)$ from Eq. (6) with the parameters (Eq. (9)) are depicted in Figure 3.

Figure 3.

Behavior of the functions $W_k\left(t\right)$ , k = 2, 3, 4.

It should be noted that in this and the subsequent examples, we use the following measures for the values: ${\mu }_k$ , $T^{\left(k\right)}$ , $t_0$ , $\tau$ , ${\tau }_k$ , ${\nu }_k$ , $T$ , $b$ , $\mathrm{{\rm E}}{\nu }_k$ (minutes, min); $\lambda$ (1/min); $\mathrm{D}{\nu }_k$ (min²). The product $\mathit{\alpha \beta }$ (as mean of ${\mu }_k$ ), where $\alpha$ is a shape parameter, $\beta$ is a scale parameter (in min).

Corollary 2. Let ${\mu }_j=T$ , $2\le j\le n$ , be a positive constant, then for $2\le k\le n$

$$W_k\left(t\right)=I\left(t_0<t\le T+t_0\right)\overline{G}\left(\left(k-1\right)T-t+t_0\right)+I\left(t>T+t_0\right),$$ E10

in particular,

$$W_2\left(t\right)=I\left(t>t_0\right)\overline{G}\left(T-t+t_0\right).$$ E11

Example 3. Let $\tau$ has density (Eq. (8)). As initial parameters, we take the following quantities:

$$\lambda =0.4, t_0=4, T=8.$$ E12

Graphs of the functions $W_2\left(t\right)$ from Eq. (11), $W_3\left(t\right)$ and $W_4\left(t\right)$ from Eq. (10) with the parameters (Eq. (12)) are depicted in Figure 4.

Figure 4.

Behavior of the functions $W_k\left(t\right)$ , k = 2, 3, 4.

Figures 3 and 4 show that in the case of constant ${\mu }_j$ , the primary delay $\tau$ practically does not affect the fourth train and all subsequent ones. This is consistent with the equality $\underset{k\to \infty }{\text{lim}}{W}_k\left(t\right)=I\left(t>t_0+T\right)$ which, as it is not difficult to verify, follows from Eq. (10).

Remark 1. It is known that the distribution of sum of the independent random variables is the convolution of their distributions. The convolution of distribution functions $F_1$ and $F_2$ is determined by the formula $\left(F_1\ast F_2\right)\left(x\right)={\int }_{-\infty }^{\infty }{F}_1\left(x-y\right){\mathit{dF}}_2\left(y\right)$ , where the integral sign means the improper Riemann-Stieltjes integral. We consider exceptionally piecewise-continuous distribution functions, then the indicated integral exists with the exception of the case when $F_1$ and $F_2$ have at least one common discontinuity point (e.g., [12]). The convolution operation is permutable. In the case, when $F_1=F_2=\dots =F_m=F$ , we shall use the following notations: $F^{\ast 2}≔F\ast F$ , $F^{\ast m}≔F\ast F^{\ast \left(m-1\right)}$ , $m\ge 2$ . By definition, we assume that $F^{\ast 1}≔F$ . The convolution $\left(f_1\ast f_2\right)\left(x\right)$ of densities $f_1$ and $f_2$ is defined as the improper Riemann integral ${\int }_{-\infty }^{\infty }{f}_1\left(x-y\right)f_2\left(y\right)\mathit{dy}$ .

Corollary 3. Let ${\mu }_j$ , $2\le j\le n$ , be independent identically distributed random variables with a continuous distribution function $\mathrm{\Psi }\left(x\right)$ . Let $\tau$ be independent of ${\mu }_j$ , $2\le j\le n$ . Then

$$W_2\left(t\right)=I\left(t>t_0\right){\int }_{-\infty }^{\infty }\overline{G}\left(z-t+t_0\right)d\mathrm{\Psi }\left(z\right),$$ E13

$$W_k\left(t\right)=I\left(t>t_0\right)\left\{\mathrm{\Psi }\left(t-t_0\right)+{\int }_{-\infty }^{\infty }\left[{\int }_{t-t_0}^{\infty }\overline{G}\left(z+u-t+t_0\right)d\mathrm{\Psi }\left(z\right)\right]d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right)\right\}, 3\le k\le n.$$ E14

Corollary 4. Let ${\mu }_j$ , $2\le j\le n$ , be independent identically distributed random variables with a density function $\psi \left(x\right)$ . Let $\tau$ be independent of all ${\mu }_j$ and has a density function $g\left(x\right)$ . Then

$$W_2\left(t\right)=I\left(t>t_0\right){\int }_{-\infty }^{\infty }\left({\int }_{z-t+t_0}^{\infty }g\left(x\right)\mathit{dx}\right)\psi \left(z\right)\mathit{dz},$$ E15

$$W_k\left(t\right)=I\left(t>t_0\right)\left\{{\int }_{-\infty }^{t-t_0}\psi \left(z\right)\mathit{dz}+{\int }_{-\infty }^{\infty }\left[{\int }_{t-t_0}^{\infty }\left({\int }_{z+u-t+t_0}^{\infty }g\left(x\right)\mathit{dx}\right)\psi \left(z\right)\mathit{dz}\right]{\psi }^{\ast \left(k-2\right)}\left(u\right)\mathit{du}\right\},$$ E16

$$3\le k\le n.$$

Remark 2. The integration limit “ $-\infty$ ” can be replaced by 0 in Corollaries 3 and 4 if ${\mu }_j\ge 0$ . On the other hand, we may consider in these corollaries the case when ${\mu }_j$ takes values of different signs. From a practical point of view, such an approach is acceptable if the probability that these random quantities take negative values is small enough. This assumption allows to consider, for example, models in which the random variables ${\mu }_j$ are normally distributed with a variance small enough and to use the property that the class of normal distributions is closed with respect to the convolution operation.

Example 4. Let $\tau$ has the density (Eq. (8)), and all ${\mu }_j$ have the same gamma density

$$\psi \left(t\right)=I\left(t>0\right)\frac{e^{-t/\beta }{t}^{\alpha -1}}{\mathrm{\Gamma }\left(\alpha \right){\beta }^{\alpha }},$$ E17

where $\alpha >0,\beta >0$ , $\mathrm{\Gamma }\left(\alpha \right)={\int }_0^{\infty }{x}^{\alpha -1}{e}^{-x}\mathit{dx}$ is gamma function. Put

$$\lambda =0.3, t_0=5, \alpha =14, \beta =0.5.$$ E18

One can show that in the example under consideration it follows from Eqs. (15) and (16) that

$$W_k\left(t\right)=I\left(t>t_0\right)\left[1-\frac{\mathrm{\Gamma }\left(\alpha ,\left(t-t_0+b\right)/\beta \right)}{\mathrm{\Gamma }\left(\alpha \right)}+{\mathit{ae}}^{\lambda \left(t-t_0+b\right)}{\left(\frac{1}{1+\mathit{\lambda \beta }}\right)}^{\left(k-1\right)\alpha }\frac{\mathrm{\Gamma }\left(\alpha ,\left(1+\mathit{\lambda \beta }\right)\left(t-t_0+b\right)/\beta \right)}{\mathrm{\Gamma }\left(\alpha \right)}\right],$$

where $\mathrm{\Gamma }\left(\alpha ,y\right)={\int }_y^{\infty }{x}^{\alpha -1}{e}^{-x}\mathit{dx}$ is incomplete gamma function. Graphs of the distribution functions $W_k\left(t\right),2\le k\le 5,$ with the parameters (Eq. (18)) are depicted in Figure 5.

Figure 5.

Behavior of the functions $W_k\left(t\right)$ , k = 2, 3, 4, 5.

It is not difficult to verify that for $W_k\left(t\right)$ from Example 4 the following formula holds:

$$\begin{array}{l}\underset{k\to \infty }{\text{lim}}{\left.W_k\left(t\right)\right|}_{a=1,b=0}\\ =\underset{k\to \infty }{\text{lim}}\left[I\left(t>t_0\right)\left(1-\frac{\mathrm{\Gamma }\left(\alpha ,\left(t-t_0\right)/\beta \right)}{\mathrm{\Gamma }\left(\alpha \right)}+e^{\lambda \left(t-t_0\right)}{\left(\frac{1}{1+\mathit{\lambda \beta }}\right)}^{\left(k-1\right)\alpha }\frac{\mathrm{\Gamma }\left(\alpha ,\left(1+\mathit{\lambda \beta }\right)\left(t-t_0\right)/\beta \right)}{\mathrm{\Gamma }\left(\alpha \right)}\right)\right]\\ =W_{\infty }\left(t\right)≔I\left(t>t_0\right)\left[1-\frac{\mathrm{\Gamma }\left(\alpha ,\left(t-t_0\right)/\beta \right)}{\mathrm{\Gamma }\left(\alpha \right)}\right].\end{array}$$

It can be seen from Figure 5, curves $W_4\left(t\right),W_5\left(t\right)$ and so on are practically merged. Hence, in the case under consideration, one can draw the following conclusion: primary delay $\tau$ affects to fifth and all successive trains approximately like on the fourth one.

Remark 3. We define the 0-fold convolution as a generalized function with the following property: the equality ${\int }_{-\infty }^{\infty }v\left(t\right){\psi }^{\ast 0}\left(t\right)\mathit{dt}=v\left(0\right)$ holds for any bounded continuous function $v\left(t\right)$ . Then, Eq. (16) for $k=2$ coincides with Eq. (15).

We do not give proofs for the statements of Section 3 because of limitations on the volume. We will make this in another work.

4. Some results on the knock-on delays

Denote by N the random number of knock-on delays (within the framework of the model under consideration).

Lemma 1. For each fixed integer m, $1\le m\le n-1$ ,

$$\mathrm{{\rm P}}\left(N\ge m\right)=\mathrm{{\rm P}}\left(\tau >\sum _{j=2}^{m+1}{\mu }_j\right).$$ E19

Proof. Easily seen: $\left\{N=0\right\}=\left\{t_0\le \tau +t_0\le T^{\left(2\right)}\right\},\left\{N=m\right\}=\left\{T^{\left(m+1\right)}<\tau +{\mathit{mt}}_0\le T^{\left(m+2\right)}-t_0\right\},$

m = 1, 2, …, n – 2, $\left\{N=n-1\right\}=\left\{\tau +\left(n-1\right)t_0>T^{\left(n\right)}\right\}.$ This implies that

$$\mathrm{{\rm P}}\left(N\ge m\right)=\mathrm{{\rm P}}\left(\tau +{\mathit{mt}}_0>T^{\left(m+1\right)}\right)=\mathrm{{\rm P}}\left(\tau >\sum _{j=2}^{m+1}{\mu }_j\right).\square$$

Here and below, the sign □ denotes the end of the proof.

The corollaries of this lemma are given below. Their proofs are simple and therefore we do not present them.

Corollary 5. If ${\mu }_j=T$ is a constant value, $2\le j\le n$ , then for every fixed integer m, $1\le m\le n-1$ , we have the equality $\mathrm{{\rm P}}\left(N\ge m\right)=\overline{G}\left(\mathit{mT}\right)$ .

Corollary 6. If ${\mu }_j=T$ is a constant value, $2\le j\le n$ , and $\tau$ is exponentially distributed with parameter $\lambda$ , then for every fixed integer m, $1\le m\le n-1$ , the following equality holds,

$$\mathrm{{\rm P}}\left(N\ge m\right)=e^{-\mathit{\lambda mT}}.$$

Corollary 7. If ${\mu }_2$ , …, ${\mu }_n$ are independent identically distributed random variables with a density function $\psi$ , then for every fixed integer m, $1\le m\le n-1$ , we have the equality

$$\mathrm{{\rm P}}\left(N\ge m\right)={\int }_{-\infty }^{\infty }\overline{G}\left(u\right){\psi }^{\ast m}\left(u\right)\mathit{du}.$$

In what follows, ${\tau }_1=\tau$ is the delay duration of the first train, ${\tau }_k$ , k = 2, …, n, is the knock-on delay of the k-th train. The problem is to find the distribution functions $G_k\left(t\right)=\mathrm{{\rm P}}\left({\tau }_k<t\right)$ , k = 2, 3, …, n. Note that the solution of this problem, which we call by the second problem, allows us to find the distribution of the deviations of the real arrival times from the planned ones.

In what follows, we will use the notation $a\vee b$ instead of $\text{max}\left(a,b\right)$ .

Theorem 3. The following formula holds:

$${\tau }_k=\left({\tau }_{k-1}-{\mu }_k\right)\vee 0, 2\le k\le n.$$ E20

Corollary 8. The following formula holds:

$${\tau }_k=\left(\tau -\sum _{j=2}^k{\mu }_j\right)\vee 0, 2\le k\le n.$$ E21

It should be noted that within the framework of our model the deviation of the real arrival time from the planned one for k-th train coincides with ${\tau }_k$ , $1\le k\le n$ . Figure 6 illustrates this statement.

Figure 6.

Deviation arrival times from the schedule: delays ${\tau }_1$ and ${\tau }_2$ .

The dotted lines (lines $\overline{1}$ and $\overline{2}$ ) represent the scheduled trajectories of trains 1 and 2, solid lines (1 and 2) depict the real trajectories taking into account the delays. It can be seen that the arrival time of the train 1 differs from the schedule at $\tau$ and the train 2 on the ${\tau }_2$ .

Denote ${\overline{\mu }}_k={\sum }_{j=2}^k{\mu }_j$ , $2\le k\le n$ . As it follows from the assumption that the random variables ${\mu }_2,\dots ,{\mu }_k$ have the same distribution function $\mathrm{\Psi }\left(t\right)$ . They are mutually independent. The random variable ${\overline{\mu }}_k$ has the distribution function ${\mathrm{\Psi }}^{\ast \left(k-1\right)}\left(t\right)$ .

Corollary 9. The distribution function of ${\tau }_k$ has the following form:

$$G_k\left(t\right)=I\left(t>0\right)\mathrm{{\rm P}}\left(\tau -{\overline{\mu }}_k<t\right), 2\le k\le n.$$ E22

The next Corollaries 10 and 11 follow from Corollary 9 in an obvious way.

Corollary 10. Let ${\mu }_j>0$ , $2\le j\le n$ be some constant values. Then

$$G_k\left(t\right)=I\left(t>0\right)G\left(t+{\overline{\mu }}_k\right).$$ E23

Corollary 11. Let ${\mu }_j=T>0$ , $2\le j\le n$ be a constant value. Then

$$G_k\left(t\right)=I\left(t>0\right)G\left(t+\left(k-1\right)T\right).$$ E24

Corollary 12. Let ${\mu }_j$ , $2\le j\le n$ be independent identically distributed random variables with a continuous distribution function $\mathrm{\Psi }\left(t\right)$ . Let $\tau$ be independent of ${\mu }_2,\dots ,{\mu }_n$ . Then

$$G_k\left(t\right)=I\left(t>0\right){\int }_{-\infty }^{\infty }G\left(t+y\right)d{\mathrm{\Psi }}^{\ast \left(k-1\right)}\left(y\right).$$ E25

Corollary 13. Let ${\mu }_j$ , $2\le j\le n$ be independent identically distributed random variables with a density function $\psi \left(t\right)$ . Let $\tau$ be independent of ${\mu }_j$ , $2\le j\le n$ and has a density function $g\left(t\right)$ . Then $G_k\left(t\right)=I\left(t>0\right){\int }_{-\infty }^{\infty }\left({\int }_{-\infty }^{t+y}g\left(z\right)\mathit{dz}\right){\psi }^{\ast \left(k-1\right)}\left(y\right)\mathit{dy}$ .

5. Proof of Theorem 3 and its corollaries

Lemma 2. The following formula is valid:

$${\tau }_2=\left(\tau -{\mu }_2\right)\vee 0.$$ E26

Proof. Let $t>0$ be the time spent by the train on the path length (distance to the place, where an unplanned stop of the train 1 occurred). We show the equality ${\tau }_2=0$ holds under the condition $\tau \le {\mu }_2$ . The departure time of the train 1 after stopping is $t+\tau$ . The time point when train 2 reaches $s$ can be written as ${\mu }_2+t_0+t$ . The knock-on delay of train 2 will not occur, i.e., ${\tau }_2=0$ , in the case, when the indicated time points are separated by the value $r\ge t_0$ , i.e., $r={\mu }_2+t_0+t-\left(t+\tau \right)\ge t_0$ , or, which is the same thing, $\tau \le {\mu }_2$ . The considered case is illustrated in Figure 7a.

Figure 7.

Two traffic scenarios: (a) lack of the knock-on delay; (b) knock-on delay occurs.

The knock-on delay of the duration ${\tau }_2=\tau -{\mu }_2$ will occur when $\tau >{\mu }_2$ . Indeed, since trains after a random stop depart simultaneously, then the equality $t+\tau ={\mu }_2+t_0+\left(t-t_0\right)+{\tau }_2$ holds, i.e., $\tau ={\mu }_2+{\tau }_2$ . The case under consideration is illustrated in Figure 7b. Thus, the validity of Eq. (26) is shown. □

Proof of Theorem 3. We shall use the method of mathematical induction. The equality (Eq. (20)) for k = 2 is established by Lemma 2. Let Eq. (20) be satisfied. We show that:

$${\tau }_{k+1}=\left({\tau }_k-{\mu }_{k+1}\right)\vee 0,2\le k+1\le n.$$ E27

It follows from the inductive hypothesis that ${\tau }_k=0$ under the condition ${\tau }_{k-1}\le {\mu }_k$ . But if the delay of the k-th train is 0, then the next train does not undergo any delay, that is, ${\tau }_{k+1}=0$ . The present case is illustrated in Figure 8.

Figure 8.

The case, when ${\tau }_{k-1}\le {\mu }_k$ : knock-on delays of k-th and consecutive trains are not observed.

In the case, when ${\tau }_{k-1}>{\mu }_k$ , a knock-on delay of the k-th train occurs and equals to ${\tau }_k={\tau }_{k-1}-{\mu }_k$ (according to the inductive hypothesis). Further, two cases are possible: either (1) a delay ${\tau }_k$ entails a delay ${\tau }_{k+1}$ , or (2) ${\tau }_{k+1}=0$ .

Case 1. If the k-th train is delayed, then (k + 1)-th one will be delayed only if ${\tau }_k>{\mu }_{k+1}$ , and its delay duration is ${\tau }_{k+1}={\tau }_k-{\mu }_{k+1}$ (this fact follows from the equality of the moments of departure of the k-th and (k + 1)-th trains after an unscheduled stop: $T^{\left(k\right)}+t-\left(k-1\right)t_0+{\tau }_k=T^{\left(k+1\right)}+t-{\mathit{kt}}_0+{\tau }_{k+1}$ ). Case 1 is illustrated in Figure 9a.

Figure 9.

Two traffic scenarios: (a) the case, when ${\tau }_{k-1}>{\mu }_k$ , ${\tau }_k>{\mu }_{k+1}$ : all trains up to (k + 1)-th are detained; (b) the case, when ${\tau }_{k-1}>{\mu }_k$ , ${\tau }_k\le {\mu }_{k+1}$ : all trains from (k + 1)-th up to n-th are not delayed.

Case 2. If the k-th train is delayed, then (k + 1)-th one will not be delayed ( ${\tau }_{k+1}=0$ ) only if ${\tau }_k\le {\mu }_{k+1}$ . Case 2 is illustrated in Figure 9b. Note that if the knock-on delay of the k-th train occurs, a conflict of the k-th train with (k + 1)-th is described similar to the interaction of trains 1 and 2 (see Lemma 2). All described cases lead to Eq. (20).□

Proof of Corollary 8. We indicate that Eq. (21) is similar to Eq. (20). According to the statement of Theorem 3, we have:

$${\tau }_2=\left(\tau -{\mu }_2\right)\vee 0, {\tau }_3=\left({\tau }_2-{\mu }_3\right)\vee 0, {\tau }_k=\left({\tau }_{k-1}-{\mu }_k\right)\vee 0.$$ E28

Using the method of mathematical induction and taking into account that ${\overline{\mu }}_{k-1}+{\mu }_k={\overline{\mu }}_k$ , we obtain Eq. (21) from Eq. (28).□

Proof of Corollary 9. It follows from Corollary 8 that ${\tau }_k=0$ if $\tau \le {\overline{\mu }}_k$ (see, e.g., Figure 8), and ${\tau }_k=\tau -{\overline{\mu }}_k$ if $\tau >{\overline{\mu }}_k$ (see, e.g., Figure 9a). Using the law of total probability, we obtain the following chain of equalities:

$$\begin{array}{l}{G}_k\left(t\right)=\mathrm{{\rm P}}\left({\tau }_k<t\right)=I\left(t>0\right)\mathrm{{\rm P}}\left({\tau }_k<t\right)\\ =I\left(t>0\right)\left(\mathrm{{\rm P}}\left(\left.{\tau }_k<t\right|\tau \le {\overline{\mu }}_k\right)\mathrm{{\rm P}}\left(\tau \le {\overline{\mu }}_k\right)+\mathrm{{\rm P}}\left(\left.{\tau }_k<t\right|\tau >{\overline{\mu }}_k\right)\mathrm{{\rm P}}\left(\tau >{\overline{\mu }}_k\right)\right)\\ =I\left(t>0\right)\mathrm{{\rm P}}\left(\tau -{\overline{\mu }}_k\le 0\right)+I\left(t>0\right)\mathrm{{\rm P}}\left(0<\tau -{\overline{\mu }}_k<t\right)=I\left(t>0\right)\mathrm{{\rm P}}\left(\tau -{\overline{\mu }}_k<t\right).\square \end{array}$$

Proof of Corollary 12. Apply the well-known assertion to Eq. (22): if $Y_1$ and $Y_2$ are independent random variables, then for any function of two variables $f\left(\cdot ,\cdot \right)$ and any $c\in R$ , the following equality holds: $\mathrm{{\rm P}}\left(f\left(Y_1,Y_2\right)<c\right)={\int }_{-\infty }^{\infty }\mathrm{{\rm P}}\left(f\left(y,Y_2\right)<c\right){\mathit{dF}}_1\left(y\right)$ , where $F_1$ is the distribution function of $Y_1$ . Consequently, $G_k\left(t\right)=I\left(t>0\right){\int }_{-\infty }^{\infty }\mathrm{{\rm P}}\left(\tau -y<t\right)d{\mathrm{\Psi }}^{\ast \left(k-1\right)}\left(y\right)$ . This implies Eq. (25).□

Proof of Corollary 13. The assertion follows from Eq. (25).□

Note that the function $G_k\left(t\right)$ has a jump at zero which is equal to:

$G_k\left(0+\right)={\int }_{-\infty }^{\infty }{G}^{+}\left(y\right)d{\mathrm{\Psi }}^{\ast \left(k-1\right)}\left(y\right)$ , where $G^{+}\left(y\right)=\underset{t\to 0+}{\text{lim}}G\left(t+y\right)$ .

In the case, when $\tau$ and ${\mu }_j$ are absolutely continuous, it follows from Eq. (25) that

$$G_k\left(t\right)=I\left(t>0\right){\int }_{-\infty }^{\infty }\left({\int }_{-\infty }^{t+y}g\left(z\right)\mathit{dz}\right){\psi }^{\ast \left(k-1\right)}\left(y\right)\mathit{dy},$$ E29

where $g\left(\cdot \right)$ and $\psi \left(\cdot \right)$ are the density functions of $\tau$ and ${\mu }_1$ , respectively, ${\psi }^{\ast j}\left(y\right)$ is the j-fold convolution of the density $\psi \left(\cdot \right)$ . In this case, we also have

$$g_k\left(t\right)≔I\left(t>0\right)G_k^{\prime }\left(t\right)=I\left(t>0\right){\int }_{-\infty }^{\infty }g\left(t+y\right){\psi }^{\ast \left(k-1\right)}\left(y\right)\mathit{dy}.$$ E30

If we assume that $\tau \ge 0$ , then we deduce from Eq. (29) that

$$G_k\left(t\right)=I\left(t>0\right){\int }_{-t}^{\infty }\left({\int }_0^{t+y}g\left(z\right)\mathit{dz}\right){\psi }^{\ast \left(k-1\right)}\left(y\right)\mathit{dy},$$ E31

and we get from Eq. (30),

$$g_k\left(t\right)=I\left(t>0\right){\int }_{-t}^{\infty }g\left(t+y\right){\psi }^{\ast \left(k-1\right)}\left(y\right)\mathit{dy}.$$ E32

6. Corollary of Theorem 2 when the distribution of primary delay is a mixture of exponential and one-point distributions

Consider the cumulative distribution function of the following type:

$$G\left(x\right)\equiv \mathrm{{\rm P}}\left(\tau <x\right)=I\left(x\ge b\right)\left(1-{\mathit{ae}}^{-\lambda \left(x-b\right)}\right),$$ E33

where $0\le a\le 1$ , $b\ge 0$ , and $\lambda >0$ are some parameters. Such distribution function is considered, for example, in [8]. It is easy to see that $G\left(x\right)=\left(1-a\right)G_0\left(x-b\right)+\mathit{aG}\left(x-b;\lambda \right)$ , where $G_0\left(x\right)$ is the distribution function of the degenerate distribution concentrated at the point $x=0$ , $G\left(x;\lambda \right)=I\left(x\ge 0\right)\left[1-e^{-\mathit{\lambda x}}\right]$ .

Let us find out the form of the distribution functions (Eqs. (13) and (14)) in the case of Eq. (33), when the function $\mathrm{\Psi }$ is continuous. In what follows, we mean that $n\ge 3$ .

Lemma 3. Let the function G be defined by Eq. (33), and $\mathrm{\Psi }$ be continuous. Then

$$W_2\left(t\right)=I\left(t>t_0\right)\left(\mathrm{\Psi }\left(t-t_0+b\right)+{\mathit{ae}}^{\lambda \left(t-t_0+b\right)}{\int }_{t-t_0+b}^{\infty }{e}^{-\mathit{\lambda z}}d\mathrm{\Psi }\left(z\right)\right),$$ E34

$$\begin{array}{l}{W}_k\left(t\right)=I\left(t>t_0\right)\left\{\mathrm{\Psi }\left(t-t_0\right)+{\mathit{ae}}^{\lambda \left(t-t_0+b\right)}\left[{\int }_b^{\infty }{e}^{-\mathit{\lambda u}}d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right){\int }_{t-t_0}^{\infty }{e}^{-\mathit{\lambda z}}d\mathrm{\Psi }\left(z\right)\right.\right.\\ \left.+{\int }_{-\infty }^b{e}^{-\mathit{\lambda u}}\left({\int }_{t-t_0-u+b}^{\infty }{e}^{-\mathit{\lambda z}}d\mathrm{\Psi }\left(z\right)\right)d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right)\right]\\ +\left.{\int }_{-\infty }^b\left(\mathrm{\Psi }\left(t-t_0-u+b\right)-\mathrm{\Psi }\left(t-t_0\right)\right)d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right)\right\},k\ge 3.\end{array}$$ E35

Proof. According to Eq. (33), one may conclude that function $\overline{G}\left(x\right)$ has a unique discontinuity point $x=b$ . Hence, the integral ${\int }_{-\infty }^{\infty }\overline{G}\left(z-t+t_0\right)d\mathrm{\Psi }\left(z\right)$ exists for any continuous distribution function $\mathrm{\Psi }$ . Note that if $\mathrm{\Psi }\left(z\right)$ had a discontinuity point $z=t_1$ , then the function $\overline{G}\left(z-t+t_0\right)$ would also be discontinuous at the point $z=t_1$ for $t=t_0+t_1-b$ , and then the considered integral would not exist (see Remark 1). Since

$$\overline{G}\left(x\right)=I\left(x<b\right)+I\left(x\ge b\right){\mathit{ae}}^{-\lambda \left(x-b\right)},$$ E36

then

$${\int }_{-\infty }^{\infty }\overline{G}\left(z-t+t_0\right)d\mathrm{\Psi }\left(z\right)=\mathrm{\Psi }\left(t-t_0+b\right)+{\mathit{ae}}^{\lambda \left(t-t_0+b\right)}{\int }_{t-t_0+b}^{\infty }{e}^{-\mathit{\lambda z}}d\mathrm{\Psi }\left(z\right).$$ E37

In accordance with Eq. (13), the relation (Eq. (34)) is proved.

Let $k\ge 3$ . It follows from Eq. (14) that

$$W_k\left(t\right)=I\left(t>t_0\right)\left[\mathrm{\Psi }\left(t-t_0\right)+{\int }_{-\infty }^{\infty }V\left(u\right)d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right)\right],$$ E38

where $V\left(u\right)={\int }_{t-t_0}^{\infty }\overline{G}\left(z+u-t+t_0\right)d\mathrm{\Psi }\left(z\right)$ . Given Eq. (36), it is easy to see that

$$V\left(u\right)=V_1\left(u\right)+V_2\left(u\right),$$ E39

$V_1\left(u\right)={\mathit{ae}}^{-\lambda \left(u-t+t_0-b\right)}{\int }_{t-t_0}^{\infty }I\left(z+u-t+t_0\ge b\right)e^{-\mathit{\lambda z}}d\mathrm{\Psi }\left(z\right)$ , $V_2\left(u\right)={\int }_{t-t_0}^{\infty }I\left(z+u-t+t_0<b\right)d\mathrm{\Psi }\left(z\right)$ .

By using equalities

$$\begin{array}{l}\left\{\left(u,z\right):u\ge b,z>t-t_0,z\ge t-t_0-u+b\right\}=\left\{\left(u,z\right):u\ge b,z>t-t_0\right\},\\ \left\{\left(u,z\right):u<b,z>t-t_0,z\ge t-t_0-u+b\right\}=\left\{\left(u,z\right):u<b,z\ge t-t_0-u+b\right\},\end{array}$$

we receive

$$\begin{array}{l}{\int }_{-\infty }^{\infty }{V}_1\left(u\right)d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right)={\mathit{ae}}^{\lambda \left(t-t_0+b\right)}\left[{\int }_b^{\infty }\left({\int }_{t-t_0}^{\infty }{e}^{-\lambda \left(z+u\right)}d\mathrm{\Psi }\left(z\right)\right)d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right)\right.\\ \left.+{\int }_{-\infty }^b\left({\int }_{t-t_0-u+b}^{\infty }{e}^{-\lambda \left(z+u\right)}d\mathrm{\Psi }\left(z\right)\right)d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right)\right].\end{array}$$ E40

Since $\left\{\left(u,z\right):u\ge b,z>t-t_0,z<t-t_0-u+b\right\}=\varnothing$ ,

$$\left\{\left(u,z\right):u<b,z>t-t_0,z<t-t_0-u+b\right\}=\left\{\left(u,z\right):u<b,t-t_0<z<t-t_0-u+b\right\},$$

then

$${\int }_{-\infty }^{\infty }{V}_2\left(u\right)d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right)={\int }_{-\infty }^b\left({\int }_{t-t_0}^{t-t_0-u+b}d\mathrm{\Psi }\left(z\right)\right)d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right).$$ E41

It follows from Eqs. (39)–(41) that

$$\begin{array}{l}{\int }_{-\infty }^{\infty }V\left(u\right)d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right)={\mathit{ae}}^{\lambda \left(t-t_0+b\right)}\left[{\int }_b^{\infty }{e}^{-\mathit{\lambda u}}d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right){\int }_{t-t_0}^{\infty }{e}^{-\mathit{\lambda z}}d\mathrm{\Psi }\left(z\right)\right.\\ \left.+{\int }_{-\infty }^b\left({\int }_{t-t_0-u+b}^{\infty }{e}^{-\mathit{\lambda z}}d\mathrm{\Psi }\left(z\right)\right)e^{-\mathit{\lambda u}}d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right)\right]+{\int }_{-\infty }^b\left(\mathrm{\Psi }\left(t-t_0-u+b\right)-\mathrm{\Psi }\left(t-t_0\right)\right)d{\mathrm{\Psi }}^{\ast \left(k-2\right)}\left(u\right).\end{array}$$ E42

The equalities Eq. (38) and Eq. (42) entail Eq. (35).□

Below we give without a proof a corollary of Lemma 3 in the case when ${\mu }_j$ are not random variables, and they are equal to the same constant.

Corollary 14. Let ${\mu }_j=T>0$ , $2\le j\le n$ , be a constant. Let function $G$ be defined by Eq. (33). Then, for $2\le k\le n$ , the following formula holds:

$$\begin{array}{l}{W}_k\left(t\right)=I\left(0\le b\le \left(k-2\right)T\right)\left[I\left(0<t-t_0\le T\right){\mathit{ae}}^{-\lambda \left(\left(k-1\right)T-t+t_0-b\right)}+I\left(t-t_0>T\right)\right]\\ +I\left(\left(k-2\right)T<b<\left(k-1\right)T\right)\left[I\left(0<t-t_0\le \left(k-1\right)T-b\right){\mathit{ae}}^{-\lambda \left(\left(k-1\right)T-t+t_0-b\right)}\right.\\ \left.+I\left(t-t_0>\left(k-1\right)T-b\right)\right]+I\left(b\ge \left(k-1\right)T\right)I\left(t>t_0\right).\end{array}$$ E43

Furthermore,

$$\begin{array}{l}\mathrm{{\rm E}}{\nu }_k=I\left(0\le b\le \left(k-2\right)T\right)\left[t_0+T-\frac{a}{\lambda }{e}^{-\lambda \left(\left(k-2\right)T-b\right)}\left(1-e^{-\mathit{\lambda T}}\right)\right]+I\left(b\ge \left(k-1\right)T\right)t_0\\ +I\left(\left(k-2\right)T<b<\left(k-1\right)T\right)\left[t_0+\left(k-1\right)T-b+\frac{a}{\lambda }\left(e^{-\lambda \left(\left(k-1\right)T-b\right)}-1\right)\right].\end{array}$$ E44

$$\begin{array}{l}\mathrm{D}{\nu }_k=I\left(0\le b\le \left(k-2\right)T\right)\frac{a}{{\lambda }^2}{e}^{-\lambda \left(\left(k-2\right)T-b\right)}\left[2\left(1-e^{-\mathit{\lambda T}}\right)-{\mathit{ae}}^{-\lambda \left(\left(k-2\right)T-b\right)}{\left(1-e^{-\mathit{\lambda T}}\right)}^2-2\lambda {\mathit{Te}}^{-\mathit{\lambda T}}\right]\\ +I\left(\left(k-2\right)T<b<\left(k-1\right)T\right)\frac{a}{{\lambda }^2}{e}^{-\lambda \left(\left(k-2\right)T-b\right)}\left[2\left(e^{\lambda \left(\left(k-2\right)T-b\right)}-e^{-\mathit{\lambda T}}\right)\right.\\ \left.-{\mathit{ae}}^{-\lambda \left(\left(k-2\right)T-b\right)}{\left(e^{\lambda \left(\left(k-2\right)T-b\right)}-e^{-\mathit{\lambda T}}\right)}^2-2\lambda \left(\left(k-1\right)T-b\right)e^{-\mathit{\lambda T}}\right].\end{array}$$ E45

Example 5. Figure 10 depicts the graphs of the functions $W_k\left(t\right)$ defined by Eq. (43) with k = 2, 3 for the following parameters:

$$a=1, b=0, \lambda =0.26, t_0=4, T=7.$$ E46

Figure 10.

Behavior of the functions $W_2\left(t\right)$ and $W_3\left(t\right)$ .

We calculated the values of $\mathrm{{\rm E}}{\nu }_k$ and $\mathrm{D}{\nu }_k$ using the formulas (44) and (45) (see Table 1).

	$k=2$	$k=3$	$k=5$	$k=8$	$k=10$
$\mathrm{{\rm E}}{\nu }_k$	7.77702	10.47779	10.98629	10.99994	10.99999
$\mathrm{D}{\nu }_k$	5.68009	2.33067	0.068156	0.00029	7.63176 × 10^−6

Table 1.

The behavior $\mathrm{{\rm E}}{\nu }_k$ and $\mathrm{D}{\nu }_k$ with growth of the parameter k.

Remark 4. It can be easily seen that the larger k, $W_k\left(t\right)$ from Eq. (43) is closer to $W\left(t\right)≔I\left(b\ge 0\right)I\left(t>t_0+T\right)$ . This agrees with Figure 10 and the formulas (44) and (45) due to which we have $\mathrm{{\rm E}}{\nu }_k\to t_0+T$ , $\mathrm{D}{\nu }_k\to 0$ as $k\to \infty$ , and also with the results of calculations in Table 1.

Let the random variable $\tau$ be distributed with the density (Eq. (33)) with parameters $a=1$ , $b=0$ . Now, we find the condition on the parameter T, under which the probability that at least m of knock-on delays will occur would not exceed a given probability p. Note that the departure headway is equal to $T+t_0$ .

According to Corollary 6, it is necessary to solve the inequality $\text{exp}\left(-\mathit{\lambda mT}\right)\le p$ . As a result, we obtain the desired condition:

$$T\ge \left(1/\left(\mathit{m\lambda }\right)\right)\text{ln}\left(1/p\right)$$ E47

(see also [13]). Denote by $T\left(m,p,\lambda \right)$ the minimal T satisfying the inequality (Eq. (47)).

Example 6. Let us fix $\lambda =0.26$ . The behavior of $T\left(m,p,\lambda \right)$ as a function of the continuous parameter m with $p=0.1$ and $p=0.05$ is shown in Figure 11a. Obviously, $T\left(m,p,\lambda \right)$ is the decreasing function with respect to the argument p. Exact calculations can be made using the formula:

$$T\left(m,p,\lambda \right)=\left(1/\left(\mathit{m\lambda }\right)\right)\text{ln}\left(1/p\right).$$ E48

Figure 11.

The behavior of the function $T\left(m,p,\lambda \right)$ : (a) $\lambda =0.26$ , $p=0.1$ , or $p=0.05$ ; (b) $p=0.1$ , $\lambda =0.26$ , or $\lambda =0.15$ .

Let $p=0.1$ . The behavior of $T\left(m,p,\lambda \right)$ as a function of the continuous parameter m with $\lambda =0.26$ and $\lambda =0.15$ is shown in Figure 11b. In accordance with Eq. (48), $T\left(m,p,\lambda \right)$ is the decreasing function with respect to the argument $\lambda$ . In the case of exponential density $g\left(t\right)$ , we have $\mathrm{{\rm E}}\tau =1/\lambda$ . Therefore, the decrease of $\lambda$ leads to increase in the average of primary delay and the departure headways (if we want to reduce the number of knock-on delays).

We also obtain the corollaries of Lemma 3 in the case when ${\mu }_j$ are distributed according to the gamma-law with the density (Eq. (17)).

Corollary 15. If primary delay $\tau$ has an exponential distribution $g\left(t\right)=I\left(t>0\right)\lambda e^{-\mathit{\lambda t}}$ and ${\mu }_k$ , $2\le k\le n$ , has the density (Eq. (17)), then the following formulas are true:

$$G_k\left(t\right)=I\left(t>0\right)\left(1-e^{-\mathit{\lambda t}}{\left(\mathit{\lambda \beta }+1\right)}^{-\left(k-1\right)\alpha }\right),$$ E49

$$g_k\left(t\right)=I\left(t>0\right){\left(\mathit{\lambda \beta }+1\right)}^{-\left(k-1\right)\alpha }\lambda e^{-\mathit{\lambda t}}.$$ E50

Remark 5. The function $g_k\left(t\right)=I\left(t>0\right)G_k^{\prime }\left(t\right)$ is not a density, in particular, because of ${\int }_{-\infty }^{\infty }{g}_k\left(t\right)\mathit{dt}={\int }_0^{\infty }{g}_k\left(t\right)\mathit{dt}=G_k\left(\infty \right)-G_k\left(0+\right)=1-h_k\ne 1$ , where $h_k$ is the jump of the function $G_k\left(t\right)$ at the origin. At the same time, the function ${\tilde{g}}_k\left(t\right)≔\frac{1}{1-h_k}{g}_k\left(t\right)$ is a density.

Corollary 15 can be reformulated as follows.

Corollary 15*. Let primary delay $\tau$ is exponentially distributed with a parameter $\lambda$ , and ${\mu }_k$ , $2\le k\le n$ , have the same gamma distribution with the density (Eq. (17)). Then, ${\tau }_k$ has the distribution function of the form Eq. (33) with $a={\left(\mathit{\lambda \beta }+1\right)}^{-\left(k-1\right)\alpha }$ , b = 0, and, consequently,

$$\mathrm{{\rm P}}\left({\tau }_k=0\right)=G_k\left(0+\right)=1-{\left(\mathit{\lambda \beta }+1\right)}^{-\left(k-1\right)\alpha }.$$

Remark 6. Let $\mathrm{{\rm P}}\left({\tau }_2=0\right)=p$ , $0<p<1.$ Then by Corollary 15*, $\mathrm{{\rm P}}\left({\tau }_k=0\right)=1-{\left(1-p\right)}^{k-1},$ $3\le k\le n$ . Hence, $\mathrm{{\rm P}}\left({\tau }_k=0\right)\to 1$ as $k\to \infty$ .

Example 7. Let ${\mu }_2,{\mu }_3,\dots$ be independent random variables having the same density function (Eq. (17)). We perform three series of experiments and investigate a behavior of distribution of the arrival time deviations ${\tau }_k$ with various combinations of parameters: $\alpha$ , $\beta$ , k. The results are presented in graphical form in Figures 12–15 The functions $G_k\left(t\right)$ are calculated by formula (49), and the functions $g_k\left(t\right)$ by formula (50). Note that product $\mathit{\alpha \beta }$ is the mean of ${\mu }_k$ . Parameter $\lambda$ is equal to 0.25 and $\mathit{\alpha \beta }=7$ as it observes in reality.

Figure 12.

Behavior of distribution $G_2\left(t\right)$ when (a) $\alpha =0.5,3,8$ and (b) $\beta =0.1,0.5,1$ .

Figure 13.

Behavior of distribution $G_3\left(t\right)$ when (a) $\alpha =0.5,3,8$ and (b) $\beta =0.1,0.5,1$ .

Figure 14.

Behavior of distribution $G_4\left(t\right)$ when (a) $\alpha =0.5,3,8$ and (b) $\beta =0.1,0.5,1$ .

Figure 15.

Behavior of distributions $G_k\left(t\right)$ (a) and densities $g_k\left(t\right)$ (b), $k=2,3,4$ and $\alpha =3$ .

7. Comparison with statistics of real train traffic

Let us consider the following random variable: the deviation of the real moment of arrival at a certain station from the scheduled one. Denote it by $\xi$ . Statistical analysis of data on this random variable, received from the Russian railways, has led to the conclusion that in many cases, they obey the modified exponential law with the distribution function of the form Eq. (33) with $b=0.$ Using data on the suburban trains of the direction “Moscow-Tver” for the period: January, 11–15, February, 1–6, 2016, we obtained a sample from the distribution of $\xi$ of the size $n=50$ with the sample mean 1.44 and sample variance 2.7. We tested the hypothesis that $\xi$ obeys distribution (Eq. (33)) with $\lambda =0.35$ and $a=0.64$ . To this end, we applied the Kolmogorov goodness-of-fit test with the significance level $\alpha =0.05$ and obtained the fit between the hypothesis and the sample data (see Figure 16).

Figure 16.

The empirical distribution function and the calculated function of the form Eq. (33).

Remark 7. It should be noted that in considered example the deviation $\xi$ is nonnegative. But in reality, it can frequently be both positive and negative. Positive values are due to arisen delay. Negative values occur due to the fact that sometimes early arrivals take place.

Remark 8. Although the hypothetical distribution function from Figure 16 is constructed for deviations without any details about the train number k, it is well correlated with the graph of the function $G_2\left(t\right)$ with $\alpha =0.5$ from Figure 12.

This allows us to assume that the distribution of the deviation $\xi$ is mainly determined by the distribution of the delay ${\tau }_2$ .

Remark 9. It was verified that if the length of the random variables ${\mu }_j$ have the same gamma distribution, any variation of the parameters of this distribution ( $\alpha$ and $\beta$ ) has a rather small influence on behavior of output distribution (see Figures 12–15).

Remark 10. Since the primary delay has a great influence on formation of the output distribution of deviations from the schedule ( ${\tau }_k$ ), then a knowledge of the primary delay distribution in each particular situation allows to predict the distribution of knock-on delays.

One important practical effect of the considered model is that it enables us to estimate the standard deviation (SD) of the actual arrival delays at the destination station. As an example, we calculated this parameter for the suburban railway line. The data analyzed were collected at the Tver station in the period of January 2016 and February 2016.

Example 8. Due to statistical data, we can consider that $\tau$ has the exponential distribution with the parameter $\lambda =0.25$ (i.e., $\tau$ has the distribution function (Eq. (33)) with $\lambda =0.25$ , $a=1$ , $b=0$ ), and ${\mu }_2$ has gamma distribution with the density function (Eq. (17)), where $\alpha =0.6$ , $\beta =11.7$ . Using formulas (49) and (50) with $k=2$ , we have:

$${\mathit{SD}}^2={\int }_{-\infty }^{\infty }{\left(t-a_2\right)}^{2}d{G}_2\left(t\right)={\int }_{-\infty }^{\infty }{t}^{2}d{G}_2\left(t\right)-a_2^2={\int }_0^{\infty }{t}^2{g}_2\left(t\right)\mathit{dt}-a_2^2\approx 10.987.$$

Here $a_2={\int }_{-\infty }^{\infty }td{G}_2\left(t\right)=\frac{1}{\lambda }{\left(\lambda \beta +1\right)}^{-0.6}\approx 1.763,{\int }_{-\infty }^{\infty }{t}^{2}d{G}_2\left(t\right)=\frac{2}{{\lambda }^2}{\left(\lambda \beta +1\right)}^{-0.6}\approx 14.088$ ,

$${\int }_0^{\infty }{t}^2{g}_2\left(t\right)\mathit{dt}-a_2^2=\frac{2}{{\lambda }^2}{\left(\lambda \beta +1\right)}^{-0.6}-\frac{1}{{\lambda }^2}{\left(\lambda \beta +1\right)}^{-1.2}\approx 10.987.$$

Thus, theoretical $\mathit{SD}\approx \sqrt{10.987}\approx 3.315\text{min}$ . This corresponds with the real statistics which shows the SD amount is 3.32 min for the mentioned station.

8. Conclusions

The mathematical model of train traffic proposed in the chapter allows us to find conditions on initial headways, which provide a smallness of frequency of a large number of delays. In other words, the formulas for the distributions of arrival headways obtained in the chapter enable to optimize the frequency of arriving train delays.

Acknowledgments

The research is funded by the JSC Russian Railways (grant 2016 for development of the scientific school).