Abstract
Robust optimal design of circulation systems (e.g., roads for vehicles or corridors for pedestrians) relies on an accurate steady-state traffic flow model that considers the effect of randomly changing environmental factors (e.g., daily periodicity and weather). Most analytical models assume that the customer interarrival time and service time of circulation facilities follow the exponential distribution with fixed rate parameters, which is unrealistic in most cases. In this paper, we develop a stationary PH /PH state-dependent queuing model in a randomly changing environment (RE), which is represented by a Markov chain. The model simultaneously considers the general randomness of arrival and service, the randomly varying rate parameters, and the state-dependent service (the travel time increases with the number of customers). The existing matrix analytic scheme (MAS) algorithm is extended to solve the proposed model because it avoids the explicit calculation of probability distributions. The space complexity of the algorithm is only linear in the number of RE states and is independent of the enormous (four-dimensional) state space of the Markov process. Its time complexity is a linear function of the product of the queue capacity and the number of RE states. Our model is validated versus simulation estimates. The obtained conditional performance measures can accurately capture the queue accumulation and dissipation and reveal the effect of randomly changing environments. Numerical experiments provide some interesting findings. (1) The proposed stationary model coincides with the transient M()/G// fluid queuing model under special conditions. (2) Under high traffic intensities, increasing the randomness in the duration time of the RE state leads to an obvious growth in the conditional queue length. (3) An increase in the facility length leads to an increase or a decrease in the average output rate, depending on whether the congestion dissipates effectively in one cycle. (4) A larger width is required to obtain the maximum average output rate for traffic demand with a greater nonuniformity.
1. Introduction
Congestions in traffic circulation systems (e.g., roads for vehicles and corridors for pedestrians) have become a global issue that has raised extensive concern. For road traffic, congestion is getting worse across the world and congestion costs are rapidly increasing in recent years [1]. As for pedestrian traffic, several large-scale stampede accidents occurred in the crowd during the past decade and caused severe casualties [2]. Efficiency and safety are the primary goals when designing traffic circulation systems because of the increased congestion. A robust design of traffic circulation facilities (such as road, corridors, and passageway) that meets the demand of customers (vehicles or pedestrians) under various situations is the key to improving efficiency and guaranteeing safety. Optimal design decisions rely on an accurate steady-state traffic flow model comprehensively considering factors that affect the performance of traffic circulation systems. Therefore, accurately modeling traffic circulation systems becomes the primary task.
However, modeling traffic circulation systems in the real world is a great challenge as they possess several complex characteristics. First, both the traffic demand and the service ability of traffic circulation facilities are random. Second, the service ability of circulation facilities is state-dependent. State dependence describes the phenomenon where an increasing number of customers leads to slower velocity due to limited land space. Third, a congestion propagation effect exists in a circulation network. Congestions in one facility will propagate in an upstream direction and affect the adjacent facilities. In addition, circulation systems are affected by randomly changing environmental factors such as daily periodicity, large events, accidents, and extreme weather conditions [3]. As a result, system parameters, such as the traffic demand and the service ability, vary randomly because of these environmental factors.
Scholars have proposed a lot of models for traffic circulation systems. The macroscopic traffic flow model focuses on the aggregate features of traffic flows. Microscopic models, such as the car-following models and lane-changing models for vehicles [4] and the microsimulation tools for pedestrians [5, 6, 7] can present more detailed characteristics but need more parameters and are less computationally efficient. To solve the accuracy-efficiency trade-off, mesoscopic models, such as the queuing model and the gas-kinetic model, are proposed. The queuing model is the most popular mesoscopic model applied in this area, and its application in modeling pedestrian or vehicle circulation systems has been demonstrated by Lovas [8], Parlar and Sharafali [9], Arita and Schadschneider [10], Tréca et al. [11], Hasani Goodarzi et al. [12], and so on. Queuing models can further be categorized into two groups: transient models and stationary models. Transient models, such as the M/G// model [13], aim to capture the dynamic performance and are often used in real-time management and control. In contrast with transient queuing models, stationary queuing models focus on obtaining the steady performance and are more useful in planning and design.
Most of the stationary queuing models in this field have considered randomness and state-dependence. For example, Yuhaski et al. [14] applied the exponential distribution to describe the random interarrival times and the service time depends on the number of customers (also called system state, denoted by ). Hu et al. [15] established a PH/PH// state-dependent queuing model where the interarrival time and service time with general randomness follow the phase-type (PH) distribution. PH distribution theoretically can approach to any positive random variable infinitely. Readers can also refer to the study of Neuts [16] and Buchholz et al. [17] for more details on PH distribution and PH-based queues. Although PH distribution can address the general randomness in arrival and service, the change in arrival or service rate is ignored in these research studies. Nowadays, more and more attention is paid to the varying parameters in traffic circulation systems, such as the study of [3], Gerum and Baykal-Gürsoy [18], and Yang et al. [19]. The robust design of transportation facilities relies on more realistic stationary models. The goal of this paper is to propose a stationary queuing model for traffic circulation systems by considering the general randomness, state dependence, and the effect of randomly changing environmental factors. We apply the PH distribution to address the general randomness in traffic demand and service ability. Meanwhile, the state-dependent service time is used to depict customers’ velocity that depends on the number of customers in the system. To further consider the effect of the randomly changing environmental factors, we introduce a random environment (RE) represented by a Markov chain. The queuing system is modulated by the RE; hence, the interarrival time and the service time change accordingly. In this study, the congestion feedback effect is ignored. Moreover, we focus on circulation facilities that are roughly rectangular and have a traffic flow in one direction. Complex situations, such as circulation networks with a feedback effect or bidirectional traffic flow, will be discussed in our future work.
The contribution of the paper is twofold. The first contribution lies in proposing a PH /PH// state-dependent queuing model in an RE. The proposed model comprehensively considers the general randomness, state dependence, and the effect of randomly changing environmental factors. Compared with existing microscopic simulation models, the proposed queuing model is time-efficient, easy to calibrate, and possesses a stable numerical solution. Compared with the stationary PH/PH()// model in the findings Hu et al. [15] which applies fixed arrival and service rates, the proposed model applies arrival and service rates varying with RE states to address the effect of randomly changing environmental factors. Unlike the transient M()/G// model based on the exponential distribution in the study of Hu et al. [13], our stationary model captures the steady performance and is more suitable for facility planning and design. Moreover, the proposed model is more practical owing to the superior characteristics of the PH distribution. Meanwhile, the duration time of each RE state in the proposed model is random as the RE is represented by a Markov chain, while the length of different phases is deterministic in the M()/G// model which cannot reflect the randomness of environment states (e.g., incident durations and peak hours). Compared with the BMAP/PH/ and MAP/PH/ models in an RE in the studies of Kim et al. [20] and Wu et al. [21], the proposed model applies a state-dependent service time and a finite space capacity to account for the congestion in traffic circulation systems. The difference between this paper and some of the most related references is shown in Table 1.
The second contribution of the paper is that we extend the matrix analytic scheme (MAS) algorithm in the study of Baumann and Sandmann [32]. Due to the introduction of the state-dependence and the RE, the proposed model has a four-dimensional state space which is considerably large. Conventional methods based on the explicit computation of the stationary distribution, such as LU-decomposition, Gaussian elimination, and Gauss–Seidel method [33, 34, 35] and the extended reversed compound agent theorem method proposed by Balsamo and Marin [36] require massive computer storage. Compared with these methods, the efficiency of the MAS algorithm is improved remarkably as the explicit computation of the stationary distribution is not necessary. Therefore, we extend it to solve the proposed model in an RE. The space complexity of the extended MAS algorithm is only linear in the number of RE states and is independent of the enormous state space of the Markov process. Its time complexity is a linear function of the product of the queue capacity and the number of RE states. We also validate the proposed model and algorithm versus simulation estimates.
The conditional performance measures of the proposed model provide insights into the varying performance of traffic circulation systems. Critical information can be revealed, such as the effect of randomly changing environmental factors and the formation, duration, and dissipation of congestions. When applying fixed arrival and service parameters, such information is greatly smoothed out. The practical value of the proposed model is embodies in offering a solid base for robust optimal design. As the proposed queuing model comprehensively considers several vital factors that influence the system performance, the design method based on it is expected to be superior in accommodating traffic demands under various situations.
The rest of this paper is organized as follows. In Section 2, we review the related works. Then, preparations for the model are presented in Section 3. The mathematical model is formulated in Section 4. Section 5 contains numerical experiments and Section 6 concludes the paper.
2. Literature Review
Queuing models with changing arrival or service rate are applied in various domains, for example, the time-dependent Markov queuing models in the study by Nelson and Taaffe [37]; the time-dependent fluid queues in the study by Liu and Whitt [22]; and queues modulated by an RE in the study by Baykal-Gürsoy et al. [3] and Kim et al. [20].
Queuing models in an RE provide a good choice for capturing the steady performance of queuing systems whose parameters change with environmental factors. An RE that represents the varying environmental factors has several states. For example, periods of heavy, medium, and light traffic can be regarded as an RE with three states (the number of RE states is 3). Every RE state lasts a random period of time before transferring to another state. At a given RE state , the system operates as a classic queue. When the RE transfers to another state, the queue will immediately change its parameters, such as arrival rate or service rate.
Queuing models in an RE have received considerable attention in the literature. Eisen and Tainiter [38] first considered a queuing process with two mean arrival and service rates. Queuing models in an RE where the interarrival time and service time are subject to exponential distribution M, such as M/M/1, M/M/, and M/M/ queuing models were later explored by Yechiali and Naor [27]; Neuts [39]; O’cinneide and Purdue [40]; Krenzler and Daduna [30]; and Yu and Liu [29]. More recent works can be found in the study by Pang et al. [41] and Naumov and Samouylov [42]. Using the M/M/ queue in an RE in modeling traffic systems was conducted by Baykal-Gürsoy et al. [3]. Generalization work was carried out by Neuts [28], who discussed the M/G/1 model under an RE. Sztrik [43] studied the M/G/ blocking system in an RE. To account for the random arrival process with generality and diversity, queues with PH distribution, MAP (Markovian arrival process), and BMAP (batch Markovian arrival process) were further explored (e.g., Krieger et al. [44]; Kim et al. [20]; Wu et al. [21]; and Kim et al. [31].
As mentioned in Section 1, state dependence exists in most traffic circulation systems. Due to the limited land space, the velocity of vehicles or pedestrians is greatly affected by the number of vehicles or pedestrians in the circulation facility. State-dependent queuing models based on exponential distribution and their applications can be found in the study by Conway and Maxwell [45]; Yuhaski et al. [14]; Cheah and Smith [46]; Smith [47]; Mitchell and Smith [48]; Weiss et al. [49]; Jain and Smith [24]; and Smith and Cruz [25]. However, the exponential arrival interval is based on the hypothesis that the arrival process is a Poisson process. This means the variation coefficient of arrival interval equals 1, which is inconsistent with the reality in most cases according to Jiang et al. [50] and Chen et al. [51]. Given the strong hypothesis on the exponentially distributed interarrival time in the M/G// state-dependent queuing model, Jiang et al. [52] developed a discrete event simulation model based on a G/G()// state-dependent queuing system. Hu et al. [15] established a PH/PH()// state-dependent queuing model for a pedestrian corridor and Zhu et al. [26] proposed the PH/PH()// state-dependent queuing network model, which can handle interarrival times and service times with various squares of the coefficient of variation (abbreviated as SCV).
We aim to propose a stationary PH/PH// state-dependent queuing model in an RE which comprehensively considers the general randomness in arrival and service, the state-dependent service time and the effect of randomly changing environmental factors. It is very difficult to solve the proposed PH()/PH// state-dependent queuing model in an RE. However, some efficient solutions provide a solid foundation. For example, the matrix analytic method is used in the study by Neuts [53] to solve the single M/G/1 queue in an RE. The extended reversed compound agent theorem based on the product-form theory is proposed in the study by Balsamo and Marin [36] for a continuous time Markov chain in an RE. Meanwhile, quasi-birth-and-death (QBD) processes have been developed to deal with a broader range of PH-based queuing models. The study by Baumann and Sandmann [32] applied the MAS algorithm to compute the stationary expectations of a level-dependent QBD process. The advantage of the MAS algorithm is that it does not explicitly compute the stationary distribution . Therefore, we extend the MAS algorithm to solve the proposed model in this paper.
3. Preparations
In this section, we first present the problem. Then, some main assumptions concerning the queuing system and the calibration of parameters are discussed.
3.1. Problem Statement
The traffic circulation facilities analyzed in this paper are presented in Figure 1. and are two dimension parameters. When the analyzed facility is for pedestrians (such as a corridor or hallway), means length and is the width (in meters). For road traffic, is the length of the road segment in miles and is the number of lanes in one direction. The capacity of the facility can be expressed as follows:where is the largest integer that is less than or equal to ; is a constant determined by the density of the appropriate context. is typically 5 peds/m2 for pedestrian traffic [14] and ranges from 185 -265 veh/mile-lane for road traffic, according to different values of jam density [24]. Note that the customary unit in each field is used for convenience.
The customers arrive at the facility with interarrival times A and the service time is denoted as B. When A and B are described by a PH distribution, then the circulation system can be described as a PH/PH// state-dependent queuing system with parallel-serial servers (details can be found in the study by Hu et al. [15]).
As mentioned before, the traffic queuing system is always affected by randomly changing environmental factors. The RE has several states and transfers randomly or circularly among these RE states. When the RE transfers from one state to another, the parameters of the queuing system, such as arrival rates and service rates, will change. For example, the arrival rate will vary with the daily period of traffic flow, and the service rate will be affected by bad weather conditions or accidents. When the RE is introduced, the interarrival time distribution at RE state is expressed as and the service time distribution at RE state is . A traffic circulation system affected by the environmental factors can then be described as a PH/PH// state-dependent queuing system in an RE. Our main objective is to establish the analytical queuing model for the PH()/PH// state-dependent queuing system in an RE and design an algorithm to solve it.
For convenience, the main notation is listed in Table 2. Note that there are two terms of state in this paper, system state and RE state . The former is the number of customers in a queuing system, while the latter is defined as the status of the RE. ”State-dependent/dependence” always refers to the system state.
3.2. Main Assumptions
The circulation facility we analyze is a straight horizontal facility. The pedestrians and vehicles are assumed to be evenly distributed in the facility. Although the model can address bidirectional or multidirectional passenger flows, we focus on a unidirectional traffic flow.
The circulation system is described as a ”first in first out” queue, which means overtaking is not allowed. This assumption is appropriate, except for the situation where the traffic density is very small. Another assumption associated with the PH/PH// queuing model is that an arriving customer will be rejected and lost when the system state is equal to the capacity . This assumption is only accurate when the customer has the choice to leave the facility.
Additionally, the interarrival time A and the total service time B are assumed to follow the simplified PH distribution. The general PH distribution has a large number of parameters. Determining all these parameters (the number can be very large) is difficult in practice. However, the mean value and the SCV (denoted as ) of a random variable can be obtained in most cases.
SCV, or the square of the coefficient of variation, is calculated by the following equation:where and are the mean and standard deviation of the random variable. When , the variables for the interarrival time and service time are constant; when , the interarrival times and service time can be described by an exponential distribution. The larger is, the greater randomness exists.
As the mean value and SCV for a random variable are available in most cases, [35] and Weerstra [54] used a simplified PH, which is determined by and . Although not the emphasis of this paper, we will introduce the simplified PH briefly for the sake of understanding. Readers can also refer to the related papers and books.
A random variable that follows the simplified PH has an -order representation as , where is a vector and is an matrix. and can be calculated according to the mean value and SCV of the random variable by the following method:(1)When , , matrix can be expressed as follows: where . For , is set to 30 to simplify the queuing analysis.(2)When , , and the matrix can be expressed as follows:
For interarrival time A and service time B that follow the simplified PH distribution, we need their mean value ( and ) and SCV ( and ) to determine their PH representation. Given an RE state , we assume is the arrival rate of pedestrians or vehicles and is the state-dependent service rate when the system state is , . Since the arrival rate is the inverse of the mean interarrival time and the service rate is the inverse of the mean service time, we have and .
Assume and represent the SCV of interarrival times and the SCV of service times, respectively. Then, the interarrival time and the service time can be expressed as PH and PH . Their PH representations PH and PH can thus be calculated according to (3) or (4). The calibration of and will be presented in the next subsection.
Here, the SCV of service time is assumed to be independent of system state . Nevertheless, Hu et al. [15] indicated that the SCV of service time is also state-dependent, as the fluctuation range of the velocity also varies with . This relation can be considered by applying the state-dependent SCV of service time in Hu et al. [15].
Next, we present assumptions concerning the RE. There are two kinds of RE in the literature, asynchronous RE and synchronous RE. A synchronous RE assumes that the change in RE state is synchronized with the customers’ departure or (and) arrival in the queuing system. More popular is the asynchronous RE, which transitions independently of the behavior in the queuing system, as considered in the study of, e.g., Krieger et al. [44] and Yang et al. [55]. We assume the RE is asynchronous in this paper, as the change in RE state, such as a change in weather conditions, is uncorrelated with the behavior of the traffic queue.
A change in the RE can change the interarrival time, service time, capacity of the queuing system, and even the service discipline. However, we consider only the first two cases in this research. The capacity and the service discipline of the queuing system always remain unchanged. In addition, the interarrival times distribution and service time distribution of a customer are assumed to depend on only the RE state at the time the arrival or service process begins [28]. In other words, once an arrival or service for a customer has been initiated, the parameters will remain the same until completion.
3.3. Parameter Calibrations
We present the calibration method briefly in this subsection. A similar but more detailed calibration process for parameters not affected by an RE can be found in the study by Hu et al. [15].
For the arrival process, the arrival rate and SCV of interarrival times at RE state are easy to obtain based on field data. In some cases where the hourly traffic volume and the peak hour factor at RE state are given, we can also calculate and according to Equations (9)-(12) in Hu et al. [15]. Then, based on (3) or (4), we can acquire and for the PH interarrival times.
For the service process, the service rate can be obtained as follows:where and are, respectively, the mean service time and the mean velocity of a customer.
By introducing the representative points that vary with the random environment to the mean velocity models proposed by Yuhaski et al. [14], we obtain the mean velocity models that consider the state dependence as well as the effect of the random environment. The linear model is as follows:
The exponential model is as follows:where
The necessary representative points are (1, ), (, ), and (, ). According to the research on pedestrian velocity by Tregenza [56]; = 1.5 m/s for = 1, = 0.64 m/s for and = 0.25 m/s for . For road traffic, Smith and Cruz [25] proposed the following coordinates: = 62.5 mph for = 1, = 48.0 mph for and = 20.0 mph for . Note: (1) The original models proposed by Yuhaski et al. [14] are based on the velocity-density curves presented in Tregenza [56], where the effect of an RE is not considered (the subscript is absent in the original models). (2) The above data for representative points are for a unidirectional traffic flow; data for bidirectional and multidirectional passengers flow can be found in the study by Mitchell and Smith [48] and Hu et al. [15]. Details on the relationship between speed, travel time, flow rate vary and density when the above state-dependent velocity models are applied can found in Hu et al. [13].
Equations (6) and (7) indicate that we can obtain the mean velocity models based on the values for , , and at each RE state . These data can be obtained by field or simulation experiments. Thus, we can obtain the service rate .
As for the SCV for service time, Hu et al. [15] proposed a similar exponential model based on the standard deviation of velocity. However, the process is complicated, and a considerable amount of data is needed for the calibration. For simplicity, is assumed to be independent of system state and to depend only on RE state in this paper. Then, by applying (3) or (4), we can calculate and for the PH service time.
4. Mathematical Model
After the preparation work, we propose the PH /PH // state-dependent queuing model in an RE in this section. Then, the main conditional performance measures of concern are derived. Furthermore, the MAS algorithm is extended to solve the proposed model and obtain the performance measures.
4.1. Model Development
Let , be an irreducible continuous time Markov chain that represents the RE. It has a finite state space . The infinitesimal generator for is . is the invariant probability vector of and satisfies , , where denotes a column vector of appropriate size consisting of all ones. As the RE is represented by an irreducible continuous time Markov chain, the duration of each RE state can follow a random distribution, such as an exponential, Erlang or PH distribution.
The interarrival time distribution PH or PH can be represented by a random variable describing the time until absorption of a Markov process , . has one absorbing state and transient states under RE state . The infinitesimal generator for can be expressed as follows:where is obtained by equations (3) or (4) and , and for , , is a nonnegative vector and satisfies .
Similarly, the service time distribution PH , or PH can be described by a continuous time Markov process , . has one absorbing state and transient states under RE state . Its infinitesimal generator is as follows:where (obtained by equations (3) or (4)) satisfies , and for , , is a nonnegative vector and .
Property 2 in the study by Hu et al. [15] indicates that if B PH or PH , then B/ PH or PH . Therefore, the state-dependent PH /PH // queuing system with parallel-serial servers can be approximately transformed into a PH /PH /1/ state-dependent queuing system with a single server, see Figure 2.
Notably, although the model is proposed for a circulation facility without obvious service desks, it is also applicable for queues with obvious parallel service desks, such as ticket vending machines and ticket windows. This is because parallel queue models can also be similarly converted into state-dependent queue models as discussed in Hu [57].
The process of the PH /PH /1/ state-dependent queuing system can be described by a four-dimensional continuous time Markov chain , where(i) is the system state, or number of customers, in the queuing system, ;(ii) is the RE state of the RE, ;(iii) is the phase of the PH arrival process, ;(iv) is the phase of the PH service process, .
The state space for the Markov chain is as follows:where stands for the system state at time , is the RE state, is the arrival phase and is the service phase. This state space can also be divided into different levels: 0, 1, 2, …, C. Enumerating the state of in lexicographic order, we obtain the following:where
Level consists of states, and level consists of states. Note that there is no service phase for level as the service process has not yet started. Level represents the states when there are customers in the system. A step from level to level represents an arrival, and a step from level to level represents a departure.
The queuing system can be described by a finite level-dependent QBD process, which leads to the following generator matrix of the Markov chain:where diag diag .
Note: is the Kronecker product and is defined as a special Kronecker product. For example, means that each element in is multiplied by an identity matrix , whose dimension depends on and , the coordinate of .
is a block matrix. The nonzero (tri-)diagonal blocks of , such as , are also block matrices with small blocks, and these small blocks represent the transition rate matrix among different RE states of the RE.
For the finite states, the stationary distribution of the Markov chain always exists [55]. The stationary probability vector of the Markov chain with generator can be obtained by solving the following global balance equation:where . represents the steady-state probability vector for level . , where represents the steady-state probability vector for system state 0 and environment state , and it can be further expanded into . represents the steady-state probability vector for level . is the steady-state probability vector for system state and environment state . The length of is .
In addition, the stationary probability satisfies , where is a level-dependent rate matrix with nonnegative elements. Meanwhile, (a by matrix with all zeros).
4.2. Performance Measures
By solving the global balance equation (15), we can obtain the stationary probability vector . The performance measures of interest can then be calculated. Generally, stationary average performance measures, such as average queue length , average blocking probability and average output rate (the average here refers to the mean value of the results for different RE states), convey all the information of interest about the queuing system. However, for queuing models in RE, these average results provide little information for analyzing circulation systems with varying parameters. The essential characteristic that system performance changes with the RE state is concealed by the average value. In practice, we are more interested in the conditional results that reflect the changing performance measures at different RE states. Therefore, this section focuses on deducing the conditional performance measures and designing a solution algorithm.
To explore the conditional performance measures, the conditional probability is necessary. Let represent the conditional probability that the system state is equal to on the condition that the RE state is ,
Based on this conditional probability, we can deduce the main conditional performance measures, such as the conditional means and variances of the queue lengths, the conditional blocking probability , and the conditional output rate at RE state . Note that these conditional performance measures represent the average performance for a given RE state.
To calculate these performance measures, one method is to first solve the global balance (15). As a limited linear equation set, (15) can be directly solved by LU-decomposition, Gaussian elimination, Gauss–Seidel method, matrix continued fractions, or matrix analytic method. Based on the stationary distribution , we can obtain the above performance measures. However, huge computer storage is required due to the explicit calculation of . To avoid excessive storage demand, we extend the MAS algorithm applied by Baumann and Sandmann [32] to compute the stationary results.
The algorithm works in two steps: First, according to , the expectation for the stationary performance measures can be expressed as follows:where is a vector function of the system state and RE state . Second, the expectation factor without the storage of is computed using the Horner rule (refer to Baumann and Sandmann [32] for details).
To satisfy the required form of the expectation factor, in equation (16) is transformed into the following form:where is denoted as a special column vector that satisfies . The length of is equal to that of . The elements in , which corresponds to , are ones, and the total number of ones is , while the others are all zeros. For example, . Now, we can express the vector function for different performance measures as follows:(1)For in equation (12), ,(2)For in equation (13), ,(3)For in equation (14), ,(4)For in equation (15),
The extended MAS algorithm is presented in Figure 3. As the proposed model is described as a four-dimensional continuous time Markov chain, the transform of the state of the RE is incorporated in the extended MAS algorithm. The time complexity of the extended MAS algorithm is a linear function of (the product of queue capacity and the number of RE states ), i.e., O . As the length of the level-dependent matrix is , where and are not larger than 30 (as stated in Section 3.2), the space complexity of the algorithm is a linear function of the number of RE states , i.e., O .
5. Numerical Experiments
In this section, we first verify the accuracy of the proposed model and then illustrate its applications in modeling road and pedestrian traffic circulation systems. Meanwhile, the sensitivity of the parameters that might influence the performance measures of the circulation systems is analyzed.
In all the experiments, for convenience, we assume the RE has an exponential duration at each RE state unless otherwise specified. The exponential duration of each RE state has a mean value equal to −1/ (the opposite of the inverse of ), where represents the diagonal entry of . For example, the morning peak period lasting 90 minutes can be represented by an RE state whose duration follows an exponential distribution with a mean value of 90 minutes ( = −1/90).
For convenience, we use = to represent the arrival rate vector, the entries of which are the arrival rates at different RE states. = is the velocity vector of a single customer under the RE. Similarly, and denote the vector for the SCV of the interarrival times and service time, respectively. In the control group of the following experiments, we use the average arrival rate or average velocity. Here, the average refers to the weighted average, where the invariant probability vector of the RE is the weighting.
5.1. Model Verification
In this subsection, we verify the accuracy of the proposed queuing model by comparing it with three other models, the M/M/1/ queuing model in an RE in Neuts [39]; the M/G// fluid queuing model developed by Hu et al. [13] and the discrete-event simulation model. In addition, the comparing between the proposed model and the PH/PH()// queuing model in Hu et al. [15] is presented in Section Sec:exper2.
First, we compare the proposed PH/PH// state-dependent queuing model in a RE with the M/M/1/ queuing model in an RE in the study by Neuts [39]. To carry out the comparison on the same benchmark, the entries of and of the proposed PH()/PH// model are set to 1. In addition, the first is set equal to 1, the second is set to be a sufficiently large number and the effect of state-dependence on service time is eliminated. The arrival rate and service rate are the same as those used in the study by Neuts [39]. The conditional means and variances of the queue length are obtained by the two methods, as shown in Table 3.
Table 3 shows that the results of the two methods are almost the same, which indicates that the proposed model can accurately reproduce the performance of queuing systems with changing arrival and service rates. In addition, the results imply that the M/M/1/ queuing model in an RE can be approximated by the PH/PH// state-dependent queuing model in an RE. This supposition is coincident with the work of Hu et al. [15]; which proved that the M/G/1/ queuing model can be approximated by the PH/PH()// state-dependent queuing model.
In the following experiment, we compare the proposed model with the transient M/G// fluid queuing model developed by Hu et al. [13]. The two models are used to analyze a road circulation system with a varying traffic demand. A road section of 850 meters is studied. The arrival rate of the M()/G// fluid queuing model changes with time determinately, while the RE state of the proposed model changes randomly. Therefore, we also explore the impact of the random length of each RE state by changing the value of the SCV of the duration time of the RE states (SCVRE) in this experiment. Four values for the SCVRE, 0.5, 1, 2 and 4, are tested. When the SCVRE equals 1, the duration time of the RE state follows an exponential distribution. When the SCVRE is 0, the duration time of each RE state is a constant value, similar to the situation in the M()/G// fluid queuing model. The entries of and of the proposed PH()/PH// model are set to 1. Two situations with low and high traffic demand are analyzed. The conditional queue length obtained by the proposed model is compared with the dynamic queue length calculated by the transient M()/G// fluid queuing model, as shown in Figure 4.
(a)
(b)
For the low traffic demand in Figure 4(a), changing the SCVRE has little effect on the queue length and the result of the two models are almost the same. In Figure 4(b), the conditional queue length increases apparently with the SCVRE during the morning peak and the result of the proposed model are all larger than that of the M/G// fluid queuing model. This indicates that increasing the randomness in the duration time of the RE state has an obvious influence on the performance measures of a circulation system under high traffic intensities.
This experiment indicates that the proposed stationary model coincides with the transient M/G// fluid queuing model under the following conditions: (1) the SCV of the interarrival times of the proposed model is equal to 1; (2) the SCVRE of the proposed model is equal to 0; and (3) the conditional queue length of the proposed model is close to 0 at the first and the last RE state. This implies that the proposed stationary PH/PH// state-dependent queuing model in a RE may be useful in conducting transient analysis after certain transformations, which will be explored in our future work. Although the two methods obtain similar results under the same condition, the proposed model is superior to the M()/G// fluid queuing model in two aspects. First, the result of the proposed model is an accurate solution while that of the M()/G// fluid queuing model is an approximation. Second, the M()/G// fluid queuing model can only deal with the situation where the SCV of the interarrival times equals 1 and the varying arrival rate changes with time deterministically. However, the proposed model can capture the performance in different situations where the SCV of the interarrival times are much larger than 1 and the arrival rate changes randomly. Both the SCV of the interarrival times and the SCVRE have a remarkable influence on the performance measures under high traffic demand (The influence of the SCV of the interarrival times will be examined in the following experiment). It also should be mentioned that the calculation time of the proposed model is 8% less than that of the M()/G// fluid queuing model.
Next, we will validate the accuracy of the proposed model by comparing it with the discrete-event simulation model. A pedestrian corridor with a length of 5 meters and a width of 3 meters is analyzed. Assume the arrival and service process of the corridor is affected by an RE with four states ( = 4) that repeat cyclically. Therefore, the invariant probability vector of the RE is = 1/4(1, 1, 1, 1), which means the probability that the RE is at state , , is 1/4. Each RE state has an exponential duration with an average value of 40 s (, ). Therefore, the infinitesimal generator can be expressed as follows:
The parameters for arrival and service under the RE are as follows: =(5, 3, 3, 2), =(1.5, 1.5, 1.5, 1.5), and =(1.5, 1.5, 1.4, 1.4), =(1.5, 1.5, 1.5, 1.5). The exponential model for pedestrian velocity in (7) is applied.
The discrete-event simulation model for such a corridor with changing arrival and service rates is set up in Simulink in MATLAB. A total of 5000 pedestrians are simulated, and the first 1000 pedestrians are excluded from the data collection. The simulation is repeated 10 times to calculate the average results. The conditional performance measures, , and , are obtained from the simulation model and the proposed model. The results, as well as the average relative error, are listed in Table 4 (the simulation result is taken as the benchmark).
The data in Table 4 indicate that the relative error is quite small, with an average value of 3.71%. This result further verifies the accuracy of the proposed model. Although the analytical model and the simulation model produce similar results, there is a big difference between the time efficiency of the two methods. On an Intel(R) Core(TM) i7-2600 CPU @ 3.40 GHz, the running time for the analytical model is 0.13 s while that for simulation is 313.57 s.
5.2. Vehicle Queues in an RE
In this subsection, we apply the proposed model to analyze vehicle queues with different arrival rates and SCVs of interarrival times. Consider a road segment of 1 mile consisting of one lane (, ) and has a unidirectional traffic flow. The constant is 200 veh/mile-lane. Therefore, the capacity of the road segment is = 200. The 24 hours of a day are treated as a circular RE with 8 states. The infinitesimal generator of the RE is as follows:and = 1/8(1, 1, 1, 1, 1, 1, 1, 1). Thus, the duration of each RE state follows an exponential distribution with an expected duration time of 3 hours. The first RE state represents the time of 4:00–7:00 am, the second, 7:00–10:00 am, and so on.
Three arrival rate vectors are applied, =(1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500), =(1200, 2200, 2100, 1700, 1400, 1300, 1100, 1000), and =(900, 2400, 2300, 2200, 1500, 1000, 900, 800). The three arrival rate vectors have the same total traffic volume in a day while the nonuniformity increases sequentially. According to our survey at Commonwealth Avenue in Singapore and the research of Zhang [58], the SCV of the interarrival times for vehicles varies from 0.5 to 15. Therefore, we set three vectors for the SCV of the interarrival times, =(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5), =(5, 5, 5, 5, 5, 5, 5, 5), and =(10, 10, 10, 10, 10, 10, 10, 10). We consecutively number the experiments with parameter pairs (, ), (, ), (, ), (, ), …, (, ), that is, Group 1, 2, …, 9. Note that the first three groups, which apply the average arrival rate, are essentially equivalent to the PH/PH // queuing model in Hu et al. [15].
All 9 groups of experiments have the same service rate vector and SCV of service time. In each experiment, the service rate and the SCV of service time are assumed to be independent of the RE state. The SCV of service time is 1. The service rate is obtained based on the exponential velocity model in (7), where the velocities for the three representative points are = 55 mph, = 48 mph, and = 20 mph (similar to the values in Jain and Smith [24]). The conditional queue length and the conditional coefficient of variation of the queue length CV are presented in Figure 5. The average time cost for one experiment is 0.8 s.
(a)
(b)
Figure 5 indicates that the nonuniformity in traffic demand causes remarkable fluctuation in the conditional performance measures. The greater the nonuniformity is, the greater the fluctuation in conditional performance will be. By contrast, the conditional mean and SCV of queue length for the first three groups are the same at all RE states. However, Figures 5(a) and 5(b) show different tendencies. In Figure 5(a), the conditional queue length of Groups 4 to 9 at RE states 2, 3, and 4 (the peak period) is much larger than that of Groups 1 to 3. The maximum conditional queue length reaches 180 during the peak period for Group 7 while that for Groups 1 to 3 is always 38. This indicates that the conditional queue length during the peak period can be tremendously underestimated by applying an average arrival rate since the fluctuating traffic intensity levels off. According to the three vectors of arrival rate and the data in Figure 5, we can see that a one percent increase in conditional arrival rate gives rise to a more than one percent increase in the conditional queue length during the peak period. The conditional CV of queue length in Figure 5(b) generally shows an opposite tendency. The conditional CV of queue length tends to be smaller at RE states with larger conditional queue length. When congestion occurs (the number of vehicles waiting in the queue is too large), the freedom of motion is almost completely limited and less fluctuation exists in queue length. When the congestion begins to dissipate (the accumulated queue length begins to decline) after the peak period, the conditional CV of queue length increases.
Figure 5 also shows the effect of the SCV of interarrival times. In Figure 5(a), the effect of the SCV of interarrival times varies. When the traffic intensity is low, the effect of the SCV of interarrival times is quite small. As the traffic intensity increases, the influence of the SCV of interarrival times tends to increase. Figure 5(a) shows that an increasing SCV of interarrival times leads to a decrease in at RE states 2, 3, and 4 in Groups 7 to 9, which appears to be unreasonable. The reason is as follows. A larger SCV of interarrival times means vehicles arrive more randomly. When the conditional queue length is close to the capacity of the road, most vehicles are rejected due to the capacity limitation. Consequently, the number of vehicles that get into the system decreases, resulting in a slight decrease in the queue length. Notably, this unexpected phenomenon occurs only when the expected queuing length is close to the capacity. Otherwise, the system has sufficient available capacity to accommodate the additional vehicles. An additional experiment shows that the PH/PH // state-dependent queuing model in the study by Hu et al. [15] presents a similar phenomenon. These results suggest a precondition shared by both models that the studied queuing system is a loss queue with limited capacity. The newly arriving vehicles will leave the system when they find that the system is full. In Figure 5(b), generally, a larger SCV of interarrival times results in a larger conditional CV of queue length.
A close look at Figure 5(a) provides additional interesting information. Although RE states 3 and 4 have smaller arrival rates than that of RE state 2, their conditional mean and CV of queue length are sometimes larger than that of RE state 2. For example, is larger than in Group 7. This result implies that the performance of an RE state is not determined solely by its arrival rate.
To verify this conclusion, we design a specific scenario that focuses on 3-mile road segment with one lane. We focus on the 4 hours around the peak period. Assume that the 4 hours are treated as an RE with 8 states (each RE state lasts 0.5 hours). Two cases are analyzed, where the arrival rate vectors are (2150, 2400, 2400, 2100, 2000, 1800, 1800, 1800) and (2150, 2600, 2600, 2100, 2000, 1800, 1800, 1800). The two arrival rate vectors differ only in the arrival rate for the peak period (RE states 2 and 3). The SCV of the interarrival times is 5 in both cases. The other settings remain the same as in the previous experiment. The conditional queue length is obtained by the proposed method. We also use the existing PH/PH// queuing model in the study by Hu et al. [15] for comparison. In the control group, the queue length for each RE state is calculated by the PH/PH()// queuing model with the arrival rate for each RE state in isolation. For example, 2150 is used to calculate the queue length for RE state 1 and 2400 is used for RE state 2 in case 1. This result is referred to as the piecewise result. The results of the two methods are presented in Figure 6.
(a)
(b)
Figure 6 shows that the conditional queue length is very different from the piecewise queue length. At the peak period of RE states 2 and 3, the piecewise queue length is much larger than the conditional queue length. After the peak period, the piecewise queue length is smaller because the piecewise queue length is the result of isolated RE states and is not affected by other RE states. Therefore, the queue length at RE states 2 and 3 is close to the capacity. The queue length is not influenced by RE state 1, where the traffic intensity is low. RE states 2 and 3 do not have a prolonged impact, and the queue length drops sharply after the peak period.
By contrast, the conditional queue length obtained by the proposed model reflects the traffic intensity of the current and previous RE states. The conditional queue length for RE states 2 and 3 is influenced by the unsaturated RE state 1. Therefore, it is much smaller than the piecewise result. Meanwhile, the traffic intensity at RE states 2 and 3 has a significant after-effect on the queue length of the subsequent RE states. Because of the aftereffect, the conditional queue length for RE states 4 to 8 is much larger than that of the piecewise result.
The difference between the two lines during RE states 4 to 8 quantifies the aftereffect of the peak period. The queue accumulated during the peak period (RE states 2 and 3) is dissipated in the subsequent RE states. The difference between Figures 6(a) and 6(b) reveals that different traffic intensities during the peak period have different aftereffects. In case 1, the difference between the two lines (the conditional queue length and the piecewise queue length) is small and is close to zero at RE state 8, which means the effect of the peak period, i.e., congestion, has almost completely dissipated. In case 2, the difference between the two lines is much larger due to the larger traffic intensity during the peak period. The difference between the two lines is still very large at RE state 8, which means the effect of the peak period is not fully dissipated and will continue to affect the performance of the subsequent RE states. In circulation systems with periodic RE, this type of information is particularly important. To avoid deterioration of the system performance, we should ensure that the accumulated queue or the after-effect of one peak period is completely dissipated before the next peak period arrives.
5.3. Varying Corridor Parameters for Pedestrian Queues in an RE
In this subsection, we examine the effect of varying width and length of a circulation facility. Here, we analyze passenger queues where the traffic demand is nonuniform, such as that in metro stations, stadiums, and comprehensive transportation hubs. According to a rough survey of Chengdu metro stations, the average headway of trains is 240 s. For outbound or transfer facilities in metro stations, most passenger arrive in the beginning 40 s of one headway and after that the arrival rate will be very small. We apply an RE with three states to model the headway of metro trains. The duration of the three RE states is, respectively, 40 s, 40 s, and 160 s. The invariant probability vector of this RE is = 1/6(1, 1, 4), and the infinitesimal generator of the RE is as follows:
First, we analyze the effect of corridor width. A corridor of 15 meters in length is analyzed. The width varies from 0.5 to 7.5 meters. We apply three arrival rate vectors with increasing nonuniformity, namely, Case 1 (2.8, 2.4, 1.1), Case 2 (3.4, 3, 0.8), and Case 3 (4.5, 3.5, 0.4). RE state 1 is the peak period and RE state 3 is the off-peak period. The three cases have the same average arrival rates as 1.6 ped/s. The PH/PH// queuing model in the study by Hu et al. [15] applying the average arrival rare is used as a control. The average output rate for varying widths is obtained and presented in Figure 7. Here, the result of our model is the average output rate for three RE states.
Figure 7 shows that the average output rate always increases with the width until it equals the average arrival rate. For the arrival rate vector with greater nonuniformity, a larger width is required to obtain the maximum average output rate. The curves also indicate that the output rate of the PH/PH// queuing model is more sensitive to width compared with our model. When the width increases, the output rate of the PH/PH()// queuing model grows faster. For Cases 1 to 3 where the arrival nonuniformity gets greater, the average output rate grows slower when the width increases. This is because the conditional output rate at RE state 3 becomes smaller sequentially for Cases 1 to 3.
Next, we study the effect of corridor length. A corridor with a constant width of 1.5 meters is analyzed. The length varies from 3 meters to 90 meters. We apply two groups of experiments, one with low traffic demand and one with high traffic demand. In Group (a), the average arrival rate is 1.6 ped/s. Three arrival rate vectors are Case 1 (2.8, 2.4, 1.1), Case 2 (3.4, 3, 0.8), and Case 3 (4.5, 3.5, 0.4). In Group (b), the average arrival rate is 1.75 ped/s. The three arrival rate vectors with increasing nonuniformity are Case 1 (2.1, 2, 1.6), Case 2 (3.1, 2.6, 1.2), and Case 3 (3.6, 3.3, 0.9). In both groups, RE state 1 is the peak period. The PH/PH// queuing model in Hu et al. [15] applying the average arrival rare is used as a control. The average output rate is calculated and presented in Figure 8.
(a)
(b)
Figures 8(a) and 8(b) show different tendencies as the length grows. This is because the growth in length has two opposite effects on the output rate. On the one hand, there is more land space and more pedestrians are allowed to get into the corridor when the length gets larger. This is the positive effect that improves the output rate. On the other hand, more pedestrians getting into the corridor leads to a longer accumulated queue and worsens the congestion. This is the negative force that brings down the output rate. In Figure 8(a) where the traffic demand is low, most of the accumulated queue dissipates at RE state 3 and the positive effect is larger. Therefore, the average output rate for both methods increases with the corridor length. In Figure 8(b), for the PH/PH// queuing model and Case 1 of our model, the traffic demand is at a high level for the three RE states and the accumulated queue cannot dissipate effectively. Therefore, the negative effect is larger and the output rate decreases when the length grows. For Case 2 and Case 3 of our model where the arrival nonuniformity is large, the positive effect is dominant at the first stage. As the arrival rate is relatively small at RE state 3 in these two cases, the congestion greatly dissipates. However, when the corridor length exceeds a certain value, the negative effect plays a leading role as the accumulated queue becomes too long and cannot dissipate timely during the off-peak period.
Meanwhile, the average output rate gets smaller from Case 1 to Case 3 as the nonuniformity in traffic demand increases (Figure 8(a)). This is because increasing nonuniformity in traffic demand leads to more pedestrians failing to get in the corridor during the peak period due to congestion. In Figure 8(b), when the length is smaller than 10, the situation is the same as that in Figure 8(a). However, when the length is large enough, the average output rate becomes larger from Case 1 to Case 3 due to the above reason.
The experiment on varying lengths suggests that nonuniformity in traffic demand affects the optimal length of circulation facilities. In the situation where the congestion accumulated at the peak period can dissipate greatly during the off-peak period, the average output rate increases with the length. However, if the congestion accumulated at the peak period cannot dissipate effectively within one cycle, the average output rate declines as the length grows.
6. Conclusion
This paper proposes a stationary PH /PH// state-dependent queuing model operating in an RE to model traffic circulation systems under the influence of randomly changing environmental factors, such as the daily traffic period, accidents, weather conditions and so on. The RE is represented by a Markov chain and the duration time of each RE state is random. At the same time, the proposed model considers the general randomness in the arrival and service process by using the PH distribution. It also explicitly models congestion by applying a state-dependent service time. The MAS algorithm is extended to solve the proposed model and obtain both the average and the conditional performance measures. The time complexity of the extended MAS algorithm is a linear function of the product of the queue capacity and the number of RE states. The space complexity of the extended algorithm is only a linear function of the number of RE states and is independent of the enormous (four-dimensional) state space of the Markov process for the queuing system. The accuracy of the proposed model is verified through comparisons with the M/M/1/ queuing model in an RE in the study by Neuts [39]; the transient M()/G// fluid queuing model in the study by Hu et al. [13] and the discrete-event simulation model.
The obtained conditional performance measures can accurately capture the queue accumulation and dissipation and reveal the effect of randomly changing environments. Such information is particularly important when analyzing circulation systems with nonuniform traffic demand. Numerical experiments on vehicle and pedestrian queues are conducted. Some interesting results are revealed.(1)The proposed stationary PH/PH// state-dependent queuing model coincides with the transient M()/G// fluid queuing model under the following conditions: the SCV of the interarrival times of the proposed model equals 1, the SCVRE of the proposed model is close to 0, and the conditional queue length of the proposed model is close to 0 at the first and the last RE state. This result implies that the proposed stationary model may be useful in conducting transient analysis after certain transformations, which will be explored in our future work.(2)Under high traffic intensities, increasing the randomness in the duration time of the RE state leads to an obvious growth in the conditional queue length. This effect is quite small under low traffic intensities.(3)Nonuniformity in traffic demand plays an important role in determining the optimal length of circulation facilities. When the congestion accumulated at the peak period can dissipate greatly during the off-peak period, the average output rate increases with the length. However, if the congestion accumulated at the peak period cannot dissipate effectively within one cycle, the average output rate declines as the length grows.(4)The average output rate always increases with the width until it equals the average arrival rate. For traffic demand with greater nonuniformity, a larger width is required to obtain the maximum average output rate.(5)The SCV of interarrival times has a remarkable effect on the conditional mean queue length under high traffic intensities. Moreover, a larger SCV of interarrival times results in a larger conditional CV of queue length. The conditional CV of queue length tends to be smaller at RE states with a larger conditional mean queue length. When the congestion begins to dissipate after the peak period, the conditional CV of queue length increases.
The proposed model serves as a solid base for robust optimal design for traffic circulation systems owing to its comprehensively considerations. In addition to the circulation facilities discussed in this paper, the proposed model can also be applied to analyze facilities with parallel servers, such as ticket vending machines and toll stations, due to the state-dependent service time. Applications in other fields, such as logistics and manufacturing, can also be explored, taking advantage of the generality of the PH distribution and the Markov RE.
Nevertheless, some issues remain to be solved in the future. We assume that only the arrival rate and service rate change with the RE. In cases such as accidents or facility disruptions, the capacity of the circulation system or even the service discipline varies. Queuing models that can account for these complex situations will be universal but also more complex. The circulation facilities are described as finite capacity queuing systems without feedback. However, in most actual situations, the newly arrived customers have to visit the system repeatedly when the system is full and they have no other option. This feedback effect will considerably influence the performance of the congested facility, as well as that of upstream facilities. The queuing network model with the feedback effect will be discussed in our following research.
Data Availability
The data that support the findings of this study are available from the corresponding author.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was funded by the National Natural Science Foundation of China (No. 71901183), the Applied Basic Research Programs of Science and Technology Department of Sichuan Province (No. 2021YJ0066), and the Talent Project of Xihua University (No. RZ2100000800). The authors are particularly grateful to Professor Der-Horng LEE from National University of Singapore and Professor Jiuh-Biing Sheu from National Taiwan University for their valuable comments.