Abstract
This study presents an approach of simulation-based optimization to the operation of the toll plaza at the car park exit. We first propose a simulation model, as the representation of the queueing system for the toll plaza with mixed-type customers and servers where the service time is dependent on the waiting time of customer. Then, a simulation-based integer programming model is developed to design more traffic-efficient yet cost-effective operation schemes. It is decomposed by a rolling horizon approach into subproblems which are all solved via the Kriging metamodel algorithm. A numerical example is presented to illustrate the model and offer insight on how to achieve traffic efficiency and cost-effectiveness.
1. Introduction
Pay-at-exit (PAE) is a manual parking fee payment system. Cars’ plates are identified by the camera when drivers enter a parking area and drivers pay a fee on exiting based on the duration of stay. The manual payment process is relatively long and could result in long queues at the exits when the departing car flows are high. In recent years, some parking operators have installed a new parking fee payment system, which is called pay-on-foot (POF). In the POF system, drivers are required to walk to pay at centralized pay stations before they return to the cars, or pay in their smartphone after scanning the Quick Response (QR) code inside the car park. Upon exiting the car park, the cars’ plates become the verified paid tickets as exit passes. The POF system could not only reduce the manpower expense, but also fasten the exit time and reduce the delays of customers at the exit.
For the convenience of customers, parking operators usually set up both the POF and PAE tollbooths in the toll plazas at the car park exits. Drivers could pay cash at the PAE tollbooths if they have not paid at centralized pay stations or in their smartphone. POF tollbooths serve POF cars exclusively. But PAE tollbooths can serve either the POF or PAE cars. Therefore, the long queues on the toll lanes will form if the numbers of POF and PAE tollbooths are not suited to the proportions of POF and PAE cars, as Figure 1(a) illustrates.
(a)
(b)
(c)
To improve the service at the toll plaza, there are more tollbooths than the approaching lanes. Hence, a transition area is needed where the number of approaching lanes has widened to the number of toll lanes in front of the toll plaza. In this area, the conflicts derived from the lane changing behaviors incur delays for drivers, as Figure 1(b) presents. In addition, when the queues in front of some tollbooths stretch to the bottleneck between transition area and approaching lanes (see Figure 1(c)), drivers have to wait before the transition area even though they plan to head to the toll lanes with a few queues.
Therefore, improper allocation of POF and PAE tollbooths in service will incur serious delays for drivers, causing traffic inefficiency of the exit. However, the efficiency evaluation of POF or PAE tollbooth is not simply multiplying average efficiency by the number of tollbooths. The reasons are threefold. First, the car type is not exclusive to the tollbooth type. PAE tollbooths can serve either the POF or PAE cars. Second, the service time is waiting-time-dependent. In most cases, PAE tollbooths have more service time than POF ones. However, if a POF car unfortunately overstays a specified time due to the long queues after paying inside the car park, the driver has to make an extra payment in the smartphone after scanning the QR code at POF tollbooth. In this case, the service time in POF tollbooth considerably increases. In addition, car parks grant drivers certain time for free parking and a driver needs to pay only if the one overstays the time. If a car stays within a predetermined time, the car does not have to pay and becomes PAE car arriving at the exit, spending only a few seconds in PAE tollbooths. Third, the positions of POF and PAE tollbooths will affect the efficiency if the queues stretch to the bottleneck between transition area and approaching lanes. Therefore, a more sophisticated mathematical model describing the queueing process is needed for the efficiency evaluation of toll plaza.
More tollbooths in service can improve the traffic efficiency, while it will lead to higher operational costs. Thus, the tradeoff between traffic efficiency and operational cost should be taken into account. In addition, since the arrival intensity of the cars fluctuates throughout the day, the allocation scheme of POF and PAE tollbooths in service should be time-varying in day-to-day operation. Car park operators usually estimate the exiting intensity of cars in a period before the beginning of the period, especially for the car park of airport where the number of cars is highly dependent on the flight schedule. Thus, they adjust the tollbooth operation scheme gradually, dividing the whole tollbooth operation problem into several stages, where each stage is only for exiting cars in the coming period. Only the plan for the near future is updated and executed. Based on this practice, the allocation is changed on a period-by-period basis.
Consequently, in this study, the operation of toll plaza is formulated as a discrete-time dynamic tollbooth allocation problem (DTAP) where the traffic efficiency and operational cost for the toll plaza are balanced in the objective function. The traffic efficiency evaluations of POF and PAE tollbooths are derived from a mathematical model that describes the queueing process at the toll plaza.
1.1. Literature Review
1.1.1. Operation Analysis of Toll Plaza
Few researchers have studied the operation of the toll plaza at a car park exit, especially the one with both POF and PAE tollbooths. Nevertheless, various studies have evaluated and proposed strategies to improve the operation of the highway toll collection system, which is similar to the parking fee collection system. The operation is described with the dynamic traffic states [1–3] or queue evolution [4–7]. After the introduction of Electronic Toll Collection (ETC) system, research efforts have been made to appropriately allocate the ETC and manual toll collection (MTC) tollbooths [8–11]. The operation of toll plaza on highway with ETC and MTC tollbooths is similar to that at car park exit with POF and PAE tollbooths. ETC tollbooths serve ETC cars exclusively while MTC tollbooths can serve either the ETC or MTC cars. However, the billing rule is the key distinction between the highway fare collection system and the parking fee collection system. For highway application, the fare is distance-based and is independent of the amount of time when a car waits at the toll plaza while for parking application, the fare is waiting-time-based. Some car parks grant drivers certain time for free parking. A driver needs to pay only if the one overstays the time. In addition, if a POF car unfortunately overstays a specified time after paying inside the car park, the driver has to make an extra payment in the smartphone after scanning the QR code at POF tollbooth.
Hence, it is important to consider the dependence of the service time in tollbooths on the waiting time in front of the toll plaza when we model the operation of the parking fee collection system.
As the aforementioned studies indicate, queueing analysis is the mainstream of evaluating the operation of toll plaza. Approaches to queueing analysis usually fall into two categories: analytical modeling and simulation. Analytical modeling is based on the queueing theory, which uses known probability distributions to describe the cars’ interarrival and tollbooths’ service patterns. A large body of literature studies exists on the queueing system with dependent services, which are related to arrival rate [12], waiting time [13], queue length [14, 15], or workload [16, 17]. Nevertheless, the queueing theory fails when the probability distribution of either cars’ interarrival or tollbooths’ service is not all mathematically explicit or when the car’s arrival rate is greater than the processing capacity in which the queue tends to be infinite. Simulation tools [4, 10, 18–20] could capture the behaviors of individual cars, including deceleration, waiting, acceleration, lane choice, lane changing, and paying. Although emulating the entire queueing process could be complex, simulation provides an access to evaluate the system performance when the interarrival of car platoons and service patterns of servers do not yield to the distributions with explicit mathematical forms, and the servers are of mixed types and are dependent on the waiting times of customers.
1.1.2. Simulation-Based Optimization Approach
The optimization for the operation of toll plaza tends to be developed based on the analytical model for the queueing system [8, 21]. Although simulation has been traditionally used as a tool to understand and experiment with a system, connecting the simulation model to the optimization engine also gives an effective solution to the optimization problem [22, 23]. Simulation-based optimization is an approach whereby an optimization engine provides the decision variables for the simulation model. The simulation model provides the results of the optimization objective function. This process will continue iteratively between the simulation model and the optimization engine until it results in a satisfactory solution or a termination due to prescribed conditions [24, 25].
Hence, the discrete-time DTAP in this paper developed by an integer programming model is formulated based on the simulation results from the queueing system. The simulation-based optimization approach has been widely used to solve the problems of congestion pricing, traffic signal control, transit scheduling, vehicle sharing, supply chain management, liner shipping, etc. [23]. However, the above simulation-based optimization problem is computationally expensive and cannot be solved by derivative-based solvers because the objective function and/or constraint set have to be treated as black boxes for the algebraic description of the simulation is not directly available. To overcome the challenges of simulation-based optimization, significant research efforts have been done to generate surrogate models of the black-box functions [26]. The surrogate method includes the response surface method, multivariate adaptive regression splines, the regression polynomials method, the Kriging method, the radial basis function (RBF) method, and the neural network method [27, 28]. The differences of these models are mainly in the approximating functions. The Kriging model is a widely used surrogate model. We use the Kriging metamodel to solve the discrete-time DTAP in this study because it is more flexible than polynomial regression models in fitting arbitrary smooth response functions and less sensitive than other metamodels (such as RBF) to a small change in the design of the experiment [23].
The complexity of discrete-time DTAP is also given by the number of integer variables, which increases with the number of periods [29]. Furthermore, it is difficult to obtain an optimal or even a near-optimal solution, since a long time horizon needs to be considered [30]. One way to handle this problem is to use a rolling horizon heuristic. Such a procedure only considers a portion of the entire time horizon, solves the reduced problem, and fixes parts of the solution. Rolling horizon schemes have been applied to several problems within operational problems under uncertainty (see Chand et al. [31] for an extensive review). Rolling horizon approaches decompose the time horizon into several stages. The offset between the starting times of two consecutive stages is defined as roll period. Roll periods also can be fixed [32–34] or event-based [30]. Usually, car park operators change the tollbooth plan at fixed intervals for the convenience of tollmen’s scheduling, even though such operation cannot give an immediate response to the intensity of exiting cars. Hence, in this study, we utilize a rolling horizon approach to decompose the discrete-time DTAP based on fixed roll period.
1.2. Objectives and Contributions
The objective of this study is to determine the optimal operation scheme of POF and PAE tollbooths for a toll plaza at the car park exit. The problem is formulated as a discrete-time DTAP that a simulation-based integer programming (SIMIP) model is formulated where the objective function balances the traffic efficiency and operational cost. Then, the discrete-time DTAP is decomposed into period-based subproblems via a rolling horizon approach. Each subproblem is iteratively solved by the Kriging metamodel-based algorithm (KMA). The contributions of this study are twofold. First, the car park exit is simulated as a queueing system where the customers and servers are of mixed types, and the service time is dependent on the waiting time of customer (hereinafter we use “customer,” “car,” and “driver” interchangeably, and “server” and “tollbooth” interchangeably). The simulation also takes into account the time-varying allocations of POF and PAE tollbooths with respect to arrival intensities of car platoons and the blockage from the queueing spillover on the adjacent toll lanes in the transition area. Second, the operation of a toll plaza is modeled as a discrete-time DTAP where a SIMIP is formulated. The problem is decomposed by the rolling horizon approach into subproblems that are solved via KMA.
The remainder of this paper is organized as follows. Section 2 develops a SIMIP model to describe the discrete-time DTAP for the toll plaza where a simulation model is formulated as the representation of the queueing system. The validation for the simulation model is also presented. A rolling horizon approach is proposed in Section 3 where the problem is decomposed into subproblems which are solved via KMA. Section 4 demonstrates the simulation results and recommended operation scheme for the toll plaza in the numerical example. Finally, the conclusions and future work are discussed in Section 5.
2. Problem Description and Model Formulation
In this study, the operation scheme of the toll plaza is determined in the time horizon , where is the ending time of the cars’ arrival at the car park exit. We discretize the time horizon into number of time periods, where . Period can also be denoted as , where , , and . can be one hour in day-to-day operations. Before we formulate the optimization problem, the following assumption is proposed initially.
Assumption 1. The allocation of PAE and POF tollbooths can only be changed at the beginning of each period.
It will heavily take up the memory of the computer if we simulate the queueing system in the entire time horizon , as shown in Figure 2(a). In addition, the allocation of PAE and POF tollbooths changes at the beginning of each period according to the arrival intensity of car platoons. Hence, we discretize the simulation into a number of subsimulations on a period-by-period basis, as Figure 2(b) presents. In simulation for each period, inputs are the car platoons arriving in this period and the queueing evolution vectors from the last period, while outputs are queueing evolution vectors and the performance measures in this period.
In Section 2.1 the SIMIP model for discrete-time DTAP is presented. Then, we describe the queueing system of the toll plaza with mixed-type customers and servers where the service time is dependent on the waiting time of customer in Section 2.2.
(a)
(b)
2.1. SIMIP Model for Discrete-Time DTAP
According to the brief analysis in the Introduction, the number of tollbooths in service, the proportion, and the location of POF and PAE tollbooths in service affect the waiting in queues in front of tollbooths and also the blockage from the queueing spillover on the adjacent tollbooths. Hence, we should find the appropriate tollbooth operation scheme which delivers the highest traffic efficiency. More tollbooths in service can improve the traffic efficiency of the toll plaza. However, it will also lead to higher operational cost. Therefore, we need to balance the traffic efficiency and the cost for the toll plaza at the car park exit.
Let denote the label set of tollbooths and where denotes the number of tollbooths (toll lanes). The number of tollbooths is described as where is the cardinality of the set. denotes the label set of tollbooths in service in the period . The set is divided into two disjoint subsets and , where and denote the label sets of POF and PAE tollbooths in service in the period , respectively.
In this paper, we use the total delay for the cars arriving in the entire time horizon as the measure of traffic efficiency, where the delay for a car is defined as the waiting time before receiving a service plus the service time in the toll plaza. The waiting for a car includes the waiting due to the busyness of tollbooth the car has chosen and the one due to the blockage of queues in front of other tollbooths. Let and denote the total delays for the cars arriving in period and entire time horizon , respectively and we havewhere represents the total delay under the scheme in period and can be obtained from the subsimulation. Hereinafter, the superscript on the notations of variables indicates that the variables are in period (i.e., in the subsimulation).
We use the manpower and electricity expense to represent the operational cost. There is a tollman for each PAE tollbooth in service. The operational cost of the toll plaza under the scheme in period , , is given bywhere and denote the manpower and electricity expense for each tollbooth in the unit time of the studied time horizon, respectively. denotes the moment when the last car leaves the toll plaza in the subsimulation. The operational cost of the toll plaza for the entire time horizon, , is given by
Let and denote the set of the indicators for tollbooth types in the entire time horizon and period , respectively, where and . if tollbooth is out of service for period . if tollbooth is a PAE one in service for period . if tollbooth is a POF one in service for period . Then, the sets , , and can be derived from the set and the variable vector . Hence, the tollbooth operation problem, i.e., discrete-time DTAP can be described as a simulation-based integer programming (SIMIP) model, presented as follows:
[SIMIP-DTAP]where and denote the optimal tollbooth allocation schemes for the entire time horizon and the period , respectively. In the objective function, and are derived from equations (1) and (3), respectively. and are the weight for operational cost and the value of time (VOT), respectively. Note that the outputs are stochastic from simulation model. Therefore, indicates the mean or median values after multiple parallel simulations. Constraint (6) is proposed to guarantee that at least one PAE tollbooth is on service, since the POF tollbooth cannot serve PAE cars while PAE tollbooth can serve both PAE and POF cars. In addition, a server cannot be closed or switch to another type in a period if it is still busy at the end of the last period, as presented in constraint (7) where denotes the queue length in front of tollbooth at the time in the subsimulation, i.e., , which is derived from the inputs and .
Here, denotes the queue length in front of tollbooth with respect to the time in the subsimulation. represents the matrix of the queue length evolution for all tollbooths in period . The queue length in our study is defined as the total number of cars staying in an area at a moment. The queue length for the exit is the total number of cars in the exit at a moment, including the cars moving, queueing, and served at tollbooth and blocked by queues at other tollbooths. The queue length in front of a tollbooth is the total number of cars which have selected this tollbooth at a moment.
The derivation of all outputs from the simulation model such as , , and is elaborated in the next section.
2.2. Simulation Model for Queueing System of Toll Plaza
In this section, we develop the simulation model to represent the queueing system on a period-by-period basis and obtain the performance of the toll plaza. In each period of the queueing system (see Figure 3), the locations of PAE and POF servers in service are fixed. The customers arrive randomly with a probability distribution and can be presented as a random arrival profile. The service times for customers in a server are also random variables that are dependent with the types of customers and the server, and the customers in queue. Customers are served based on First-In-First-Out (FIFO) rule in a server. In addition, the customers’ behaviors before receiving a service are also taken into account such as server selection and queueing due to the blockage from the queueing spillover on the adjacent servers.
In addition, car parks grant drivers certain time for free parking. If a car stays within a predetermined time, the car does not have to pay and becomes PAE car arriving at the exit. In addition, POF cars to exit the car park within a required time after they pay at centralized pay stations or in their smartphone. Otherwise, the POF cars have to pay by cash at PAE tollbooths or in smartphone at POF tollbooths for the extra times. The above regulations are also captured in the simulation.
2.2.1. Assumptions and Simulation Inputs
To simplify the representation of the queueing system, the following behavioral assumptions are proposed in the simulation.
Assumption 2. Travel speeds of cars are homogeneous. Decelerations and accelerations of cars are not considered.
Assumption 3. Before a driver arrives at the toll lanes, he or she selects the tollbooth in service with the shortest expected waiting time to wait. Once choosing a toll lane to queue on, he or she will not switch to another tollbooth.
Assumption 4. Blockage from the queueing spillover on the adjacent toll lanes is evaluated, while the delay due to the conflicts from the lane changing is not taken into account.
The exit is divided by three cross sections, denoted as , , and , respectively, in Figure 3. is the place where drivers decide which tollbooth they will wait for exit. It is located upstream of the exit and is assumed to be rarely reached by queues. is the place where cars start to be served at the toll plaza and is the place where cars leave the toll plaza. The studied area is determined as the one between and .
and denote the numbers of approaching lanes and toll lanes, respectively. There are more toll lanes at the left-hand side of approaching lanes and more toll lanes at the right-hand side of approaching lanes. Hence, , as presented in Figure 3.
We denote and as the number and arrival rate of car platoons in period , respectively, and we have 𝜆𝑛𝛿𝑇 = 𝐾𝑙[𝑛]. The arrival pattern is generated via a given stochastic process with the arrival rate, as the representation of the cumulative profile at for toll plaza, denoted by , as the cumulative profile at for toll plaza, where denotes the time in the th subsimulation. Let denote the number of cars choosing tollbooth in the period and we have . Let denote the moment in the subsimulation when the last car leaves the tollbooth . In addition, we define that and .
In conclusion, the main inputs in the simulation are the number of cars in period , , and the label sets of POF and PAE tollbooths in service in the period , and .
2.2.2. Driver Behaviors and Output Evaluation
Now we focus on driver behaviors in the simulation and the output evaluation. Consider the car in the period (denoted by ) which arrives at at time with the parking duration . The parking duration is randomly generated via empirical distribution.
Firstly, the type of the car is determined. If , the driver has the chance to be free of charge if one exits the toll plaza immediately, where denotes the free-parking time of the car park. He or she does not need to pay at centralized pay stations or pay in the smartphone after scanning the QR code inside the car park and the car is PAE car. Otherwise, the driver must pay. Whether the car is PAE or POF car is determined using Bernoulli distribution in terms of his or her preference to the POF system. Let denote the proportion of drivers who prefer to the POF system and the car type is generated by .
Then, the driver compares the queues of tollbooths in service. Since the service time of a POF server is usually shorter than that of a PAE server, the driver chooses the tollbooth with the shortest waiting time that he or she estimates (see Assumption 3). If more than one server has the shortest waiting time, the driver can select any one of those servers with identical probabilities. Let denote the estimated waiting time for tollbooth by the car . Here, we assume that drivers only know the average not the variance of a tollbooth’s service times. Hence, is given bywhere denotes the average service time for tollbooth . denotes the queue length in front of tollbooth with respect to the time in the subsimulation. is the summation of the newly formed queue length in the period and the queue length in the last subsimulation at the same time, . The detail for the derivation of is presented in Appendix A.
Once the car makes a decision at , one is labelled as the car choosing tollbooth and denoted by . Then, the car may be blocked and wait for the discharge of the queues from other tollbooths in the transition area, and afterwards, it queues on the toll lane until tollbooth start to serve it. Let , , and denote the travel time from to tollbooth when there is no queue in front of tollbooth , the waiting time for the car due to the blockage of the queues in front of other tollbooths, and the waiting time for the car due to the busyness of tollbooth , respectively. The derivation of and is presented in Appendix A. We denote as the waiting time for the car and we have
The service time for car is determined by the type of tollbooth , and , , and of the car. Let and denote the service time and delay for the car , respectively. We have
After capturing the movement of each individual car, the performance measures of the toll plaza can be described. We use the total delay of car platoons in order to reflect the overall time costs of drivers’ waiting. Total delay for the cars arriving in period , , is given by
Total delay for the cars arriving in the entire time horizon , , is given by
In addition, other outputs from the subsimulation are obtained that the next subsimulation needs. They are the vector of the moments when the last car leaves each server, and the matrix of the queue length evolution for each server. Vector of the moments when the last car leaves corresponding servers in period , , is given bywhere for tollbooth , we have . Matrix of the queue length evolution for all tollbooths in period , , is given bywhere denotes the moment when the last car leaves the toll plaza in the subsimulation and . If the tollbooth has , then we have , . In addition, for tollbooth , we have , .
Matrix of the queue length evolution for all tollbooths in the entire horizon, , is given bywhere
In conclusion, the main outputs from the subsimulation can be described as black-box functions, presented as follows:where and are defined.
The total delay for car platoons from the simulation is presented as follows:
2.2.3. Simulation Procedure
The subsimulation () is performed as follows (Figures 4 and 5; Table 1):
(a)
(b)
(a)
(b)
(c)
[SIM-n] Step 0: initialization for environment. Give the sets of POF and PAE tollbooths in service, and , respectively, the free-parking time, , the no-extra-pay time after POF cars paying inside the car park, , and the travel time from to tollbooth when there is no queue in front of tollbooth , . Input from the last subsimulation. Step 1: arrival information generation. Generate the arrival time for each car and the total number of cars arriving, in the period at via a stochastic process with arrival intensity , parking duration via empirical distribution in Figure 4(a), and the travel time from car park to , via empirical distribution in Figure 4(b). Set . Step 2: car type generation. If , the car is labelled as PAE one. If , the type of car is determined using Bernoulli distribution according to the surveyed preference of drivers to the POF system; Step 3: tollbooth selection. PAE car chooses the PAE tollbooth with the shortest estimated waiting time at time . That is, find a PAE tollbooth label with the probability . POF car could use either the PAE or POF tollbooth and will choose the tollbooth with the shortest estimated waiting time. That is, use equation (1) to find a tollbooth label with the probability . Step 4: waiting time computation. Once the car becomes the car chooses the tollbooth, i.e., , we can compute the waiting time using equation (A.8) and the waiting time using equation (A.11). Step 5: service time generation. Generate the time to be served, , at tollbooth . The service time depends on the type of tollbooths, the type of cars, and the computed waiting time in equation (9), as presented in Table 1. Step 6: move to next arriving car. Obtain the delay of the car , , in equation (10). And . Step 7: simulation termination criterion. If , simulation terminates and give the outputs , , and ; return to Step 2 otherwise.
2.2.4. Validation for Simulation Model
(1) Studied Car Park. Terminal 3A of Chongqing Jiangbei International Airport, China came into operation in August 2017. According to the planning document, there will be approximately 8000 and 12300 passengers arriving at the terminal during peak hours in years 2020 and 2040, respectively. In future, a large number of cars will come to pick up passengers in the car park especially when subway is out of service after 11:00 p.m. To accommodate the considerably heavy flows leaving the car park in the future, the toll plaza for the car park exit is designed with 16 tollbooths: 13 for cars, 2 for coaches, and 1 for motorcycles, as Figure 6(a) presents. However, there is no need for operators to keep all 13 tollbooths in service since the flows from the park are not high so far. Usually, 2 PAE and 2 POF tollbooths are in service.
(a)
(b)
Video observation was adopted to capture the arrival and service features of car platoons at the toll plaza for this car park exit. Travel times from car park to are estimated via the distances between all parking spaces and . However, the car park did not provide the data of parking durations for us due to the lack of detectors for parking spaces. Fortunately, we investigated the operation of P7 car park for Shanghai Hongqiao International Airport in 2017, and large quantities of parking duration samples were manually obtained by videos at the parking spaces. The data were used in this study and assumed to be able to reflect the parking condition in the car park for Terminal 3A of Chongqing Jiangbei International Airport.
In addition, the P7 car park for Shanghai Hongqiao International Airport has 4 PAE tollbooths and 1 staff only one in 2015, as Figure 6(b) illustrates. The billing rules between two parks are the same. Hence, the service times for PAE tollbooths in this park can be supplemented to the database of this study.
In sum, to develop and calibrate the simulation model in our study, the arrival profiles were collected from the car park exit for Terminal 3A of Chongqing Jiangbei International Airport. The parking durations were obtained in the P7 car park for Shanghai Hongqiao International Airport. The data for service times in POF and PAE tollbooths were observed in both two car park exits. The data are presented in Appendix B.
We use queue length to validate the simulation model via hypothesis test. The queue length evolutions are recorded at the same time we collected cars’ arrival data by video at 23:00–01:00 on fourteen consecutive days in November 2018, as presented in the previous section. During the observation, 1 PAE and 3 POF tollbooths are in service where the sets for PAE and POF tollbooths in service are and , respectively. The data of arrival and queues are divided by one hour, and then, 28 samples are captured.
(2) Validation. We validate the simulation model via hypothesis test. The queue length evolutions are recorded at the same time we collected cars’ arrival data by video at 23:00–01:00 on fourteen consecutive days in November 2018, as presented in Appendix. During the observation, 1 PAE and 3 POF tollbooths are in service where the sets for PAE and POF tollbooths in service are and , respectively. The data of arrival and queues are divided by one hour, and then, 28 samples are captured.
In each sample, the arrival profile is fixed as presented in Figure 7. We can adopt queue length at per second in the evolution or the average queue length to compare the results from 1000 parallel simulations and the one from observation. However, it is time-consuming and unnecessary for the validation to compare the queue length at every second in the evolution (the case of 432 cars arriving is shown in Figure 8). Hence, we use the average queue length (AQ) as the key measure for validation in 28 samples.
The hypothesis testing for each sample (or experiment) () is presented as , where and denote the AQ from observation and the mean of AQs from 1000 parallel simulations, respectively. The boxplots for simulation results and corresponding observed results are shown in Figure 9. Given the 0.95 confidence level, the significance, lower bound (LB), and upper bound (UB) of confidence interval (CI) for each sample are listed in Table 2.
According to Table 2, only four hypotheses () are rejected while others are accepted. The simulation model is valid for the description of the queueing system at car park exit. Admittedly, the average queue length as the key measure is sensitive since it only can be integer. It is nevertheless easy to observe in experiments and is an acceptable alternative for validation.
3. Solution Algorithm
To solve the simulation-based integer programming model for discrete-time DTAP, SIMIP-DTAP, we propose a solution algorithm with two steps. The first step is to decompose the discrete-time DTAP into subproblems by periods where the arrival intensity and allocation scheme are static. The second one is to solve each period-based subproblem. The two steps are detailed in the next two sections, respectively.
3.1. Rolling Horizon Approach to Decomposition of Discrete-Time DTAP
This section presents a rolling horizon (RH) approach to decompose the discrete-time DTAP, handling scheme dynamics and output stochasticity from simulation. The complexity of discrete-time DTAP is mainly given by the number of integer variables, which increases with the number of periods and subsimulations [29]. In addition, it is difficult to obtain an optimal or even a near-optimal solution, since a long time horizon needs to be considered [30]. Furthermore, car park operators usually estimate the exiting intensity of cars in a period before the beginning of the period. Thus, they divide the whole tollbooth operation problem into several stages, where each stage is only for exiting cars in the coming period. Only the plan for the near future is updated and executed. Based on this practice, an RH approach is utilized to decompose our problem, which is helpful in decreasing the computation time. Note that considering the required computation time, the model in a real-world tollbooth operation needs to be executed a few minutes before the onset of each roll period.
The RH approach is a decomposition method where the allocation problem at a stage is developed conditionally on the optimal allocation scheme from the last stage. Therefore, the discrete-time DTAP is decomposed into static subproblems on the stage-by-stage basis. According to Assumption 1, this decomposition coincides with discretization for the simulation on the period-by-period basis. The RH approach is presented in Figure 10, where only three stages are shown for instance. The parameter is the stage horizon which is stochastic, i.e., the time elapsed from when the period starts to when the last car leaves the toll plaza in the subsimulation. The parameter is the roll period for each stage and ; The parameter denotes the start time for the stage (). Hereinafter, we use and interchangeably.
The computational steps of RH approach to solving discrete-time DTAP are presented as follows.
[RH-DTAP] Step 0: initialization. Give the stage , subsimulation , time , optimal allocation scheme at stage , , and optimal output from stage , . Step 1: solve the static tollbooth allocation problem (STAP-n) based on inputs and , and subsimulation at stage . Obtain the optimal allocation scheme and output . Step 2: update the stage, subsimulation, and time, , , and . Step 3: iteration termination criterion. If , iteration terminates and give the optimal allocation scheme sequence ; return to Step 1 otherwise.
In RH, the subproblem STAP-n () is also described as a SIMIP model, presented as follows:
[SIMIP-STAP-n]where and are the optimal results from [SIMIP-STAP-] as well as the constants in [STAP-]. We define that and . is obtained using the following equation:
Note that here is different from but close to that in equation (2). The solution of the subproblem STAP-n is explicated in the next section.
3.2. Kriging Metamodel Algorithm for Solving Period-Based Subproblem
Due to the stochasticity of the outputs from the simulation, the SIMIP-STAP-n is computationally difficult to yield the exact solution. An approximation to the solution can be derived by existing algorithms. In this paper, the Kriging metamodel is used as a surrogate. The metamodel method is favored in many literature studies due to the determinacy of the response, thus making it computationally efficient. The metamodel is an approximation that is inexpensive to compute to describe the original model with the complicated and computationally expensive response.
3.2.1. Framework
The integer (global) optimization problem can be written in the following form:where denotes the objective function defined in equation (19). and denote the lower and upper bounds on the variable , respectively. The optimization model defined in (24) is replaced by the sum of a constant value and a Gaussian random error term as follows:where is the mean value of and is a stationary Gaussian random process with mean zero, variance , and nonzero covariance. The covariance is expressed aswhere is the covariance between two sample points and . is the Kriging function, or Gaussian exponential correlation function, equivalently. The Kriging function is defined as follows:where denotes the vector of unknown correlation parameters and and denote the variable in and , respectively. denotes the dimension of and here in this paper. Note that the parameter can be regarded as measuring the importance of the variables and . The larger of , the greater the Euclidean distance between and , hence the lower the correlation between them.
Let denote the number of the known sample points for variable . The predictive estimate of the response function at unknown point is given bywhere , is the unit column, and is the estimated correlation matrix and is diagonally symmetric wherewhere
is the correlation vector between the unknown point and sample points where
Therefore, in order to yield the value of , the estimators and , , need to be obtained, given bywhere is the determinant of matrix [35]. According to equation (33), is obtained based on , , , while , , , are estimated via in terms of equation (34). Hence, as claimed in Kleijnen [36], this estimation is a difficult mathematical problem, and the divide-and-conquer methodology [23] is used here for estimation. The algorithm to estimate the above parameters PE is presented as follows:
[PE] Step 0: initialize the parameters of , , , , and obtain the defined in equation (31) and defined in equation (30) Step 1: compute the value of in terms of equation (33) Step 2: substitute obtained in Step 0 and obtained in Step 1 into equation (34), solve equation (34), and update , , , Step 3: update and with the updated , , , , and substitute this updated into equation (33) Step 4: if the convergence criteria have not been satisfied, then return to Step 1; otherwise, stop and output the estimators of and , , ,
If the problem has other constraints in addition to the bound constraint (25), it can be formulated as the following form:subject to
In order to use Kriging metamodel, the nonpositive constraint (36) is added as the penalty term to the objective function value at every candidate point , [37, 38] which is rewritten aswhere denotes the penalty augmented objective function (PAOF). represents the worst feasible objective function value so far. is the constraint violation function, and denotes the penalty factor which is positive and adjusted. This definition guarantees that the penalty augmented function values of the infeasible points are larger than the worst feasible objective function value.
3.2.2. Kriging Metamodel Algorithm
Based on the principle interpreted in Section 3.2.1, SIMIP-STAP-n can be transformed into an integer (global) optimization problem through approximating the outputs of simulation by metamodels.
The Kriging metamodel algorithm (KMA) starts with an initial experimental design. The initial experimental design is to produce points which are repeatedly generated via a symmetric Latin hypercube design [39]. In Latin hypercube design, each variable of is stratified into equal intervals, and in each subinterval, a sample point is randomly generated.
Then, the optimization works iteratively and, in each iteration, new points for calculating the next expensive function are added to the set of already sampled points. A candidate point-based approach is applied to determine the new sample site. Candidates are generated by (a) uniformly selecting points and (b) perturbing the best point found so far. The candidate point in the iteration is denoted by , .
In the given set of randomly generated candidate points, the parameters of the Kriging metamodel are computed by the algorithm PE and the is then used to approximate the at every candidate point . A good candidate point for function evaluation ideally should have a low estimated function value (since the goal is to minimize it) and should be far away from previously evaluated points (since this promotes a more global search). Hence, the selected point has the best weighted score based on two criteria: (1) estimated function value obtained from the response surface model (response surface criterion ) and (2) minimum distance from previously evaluated points (distance criterion ) [40]. The weighted score is computed aswhere , is the weight for the response surface criterion, and is the weight for the distance criterion. The response surface criterion for every candidate point is given bywhere and . The distance criterion for every candidate point is given bywhere , , and .
The candidate point with the lowest weighted score is chosen as the new sample site in the next iteration until the maximum number of allowed function evaluations has been reached.
4. Numerical Example
In this section, a numerical example is presented for the car park exit for Terminal 3A of Chongqing Jiangbei International Airport, China, as Figure 6(a) demonstrates. The toll plaza for the car park exit has 7-lane approach and 13 tollbooths for cars where , , and . We number the tollbooths from No. 1 to No. 13 consecutively (from the left- to right-hand side in the toll plaza). for the toll lanes from No. 7 to No. 12 are 7, 6, 5, 5, 4, and 4 cars, respectively. According to our investigation, nearly 70% of drivers prefer the POF system.
The datasets for the random parameters in the simulation are established according to Appendix B. The algorithm is coded in MATLAB and run on a personal computer with Intel (R) Core (TM) i7-8700 3.20 GHz CPU.
4.1. Convergence Analysis for Simulation Model
More parallel evaluations for the simulation model will lead to more accurate distribution of results, but they are computationally expensive. Therefore, in the second part, we need to perform convergence analysis such that an appropriate number of evaluations are found which balances the accuracy and computation performance.
In the convergence analysis, average and total delays are chosen as the key measures under the arrival intensities of 1000, 2000, and 3000 cars in an hour. To accommodate the arrival intensity, 2 PAE and 5 POF tollbooths are assumed to be in service where the sets for PAE and POF tollbooths in service are and , respectively. The median, mean, and standard deviation of the key measures are obtained with different parallel evaluations for the simulations, as demonstrated in Figures 11–13.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
As Figures 11–14 illustrate, after 150 parallel evaluations of simulations, the convergence patterns for the median, mean, and standard deviation of key measures are smooth. Thus, 150 parallel simulations are adequate for the subsequent performance analysis and optimization for the operation of toll plaza. Since the distributions of average and total delays represent either skewed or normal patterns, the median rather than mean is considered as the statistical quantity in the study when we develop SIMIP-DTAP.
(a)
(b)
(c)
4.2. Appropriate Operation Schemes for Toll Plaza
In this part, we find the solution for SIMIP-DTAP as the recommended tollbooth operation schemes, with different weights of operational costs. The horizon is 8 hours. The sequence of arrival intensities is cars per hour. These intensities are higher than the current situation shown in Section 2.2 and are for future consideration. The manpower and electricity expense for each tollbooth are ￥20.00 and ￥1.00 per hour, respectively, where ￥ denotes the Chinese currency unit. The VOT for all people in one car is ￥50.00 per hour. In KMA for each stage of RH-DTAP, 200 function evaluations are allowed at the maximum, with 20 initial points and 20 new sample points in each iteration.
The recommended tollbooth operation scheme derived from the solution of discrete-time DTAP and its performance are investigated with different weights of operational cost, as Figures 15–17 presented. It is worth noting that(i)PAE tollbooths are usually more welcome than POF ones in spite of the different weights between traffic efficiency and cost-effectiveness(ii)As the arrival intensity increases, the proportion of PAE tollbooths to POF ones grows
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
These two statements are drawn because PAE tollbooths are more flexible that they can serve both PAE and POF cars. When all POF tollbooths have long queues, POF cars can choose PAE tollbooths as alternatives. The average service times for POF and PAE tollbooths are 4 and 17 seconds, respectively. The ratio of POF to PAE cars is 7 : 3 if we do not take into account the cars who supposed to be free-parking and naturally become PAE cars. Therefore, the estimated balanced ratio (EBR) of POF to PAE tollbooths is roughly obtained as 1 : 1.8. When the arrival intensity is low, the ratio of POF to PAE tollbooths is close to EBR. When the arrival intensity is high, however, the ratio of POF to PAE tollbooths is far higher than EBR. The toll plaza still needs more PAE tollbooths than EBR to satisfy the heavy demand even if car operators want to control operational cost.
In addition, the lengths of all periods are identical and the allocation of PAE and POF tollbooths is assumed to be changed only at the beginning of each period. Therefore, the proposed allocation schemes do not have immediate response to the event when a bursty flow arrives in the middle of the period. An event-based dynamic tollbooth allocation problem will be proposed in the further research.
5. Conclusions and Recommendations
This paper investigated traffic-efficient yet cost-effective operation of POF and PAE tollbooths at the toll plaza for the car park, formulated as a discrete-time dynamic tollbooth allocation problem (DTAP) where a simulation-based integer programming model is developed. To describe the traffic dynamic at the car park exit, a simulation model was first proposed to represent the queueing system of the toll plaza with the mixed-type customers (cars) and servers (tollbooths), where the service time is dependent on the waiting times of customers. The simulation also took into account the time-varying tollbooth allocations with respect to arrival intensities of car platoons, and the queueing due to the blockage from the queueing spillover on the adjacent servers in the transition area before the toll plaza. Then, the simulation model is validated at the studied car park exit. Based on the results from the simulation, an integer programming model was developed to minimize a weighted sum of the traffic efficiency measure and operational cost for operation schemes of toll plaza. It was decomposed by rolling horizon approach into period-based subproblems that are iteratively solved via the Kriging metamodel algorithm.
A numerical example from the studied car park exit was presented to illustrate the proposed simulation model and discrete-time DTAP. The recommended operation schemes are derived. In recommended schemes, PAE tollbooths are usually more welcome than POF ones in all recommended schemes with different weights between traffic efficiency and cost-effectiveness, since PAE tollbooths are more flexible that they can serve both PAE and POF cars. In particular, when the arrival intensity is high, the toll plaza needs far more PAE tollbooths to satisfy the heavy demand even if car operator wants to control operational cost.
Further research directions can be derived from the simulation methodology proposed. Either the decelerations or accelerations in the driving behaviors of car following, lane changing, and gap acceptance are not taken into account in the simulation. Neither are the competitions of right-of-way among all cars. The microscopic simulation needs to be more sophisticated. In addition, the lengths of all periods are identical, and the allocation of PAE and POF tollbooths is assumed to be changed only at the beginning of each period. Therefore, the proposed allocation schemes do not have immediate response to the event when a bursty flow arrives in the middle of the period. An event-based dynamic tollbooth allocation problem will be proposed in the further research.
Appendix
A. Evaluation of Queue Length in System and Waiting Time for Individual Car
In order to help understand the relationships of the inputs and outputs among all periods, we first establish the arrival and departure cumulative profiles of the cross sections for the entire time horizon and of each period . Then, the relationship among profiles in periods is described, and the performance measures in each period are quantified, respectively.
The cumulative profiles for tollbooth at , , and are denoted as , , and , respectively. As demonstrated in Figure 18, these cumulative profiles are formulated in the rectangular coordinate where the represents the time and the represents the cumulative number of cars choosing tollbooth in the simulation. These profiles are further divided by the period, as shown in Figure 6. The cumulative profiles for tollbooth at , , and are denoted as , , and , respectively, and are formulated in the rectangular coordinate where the -axis represents the time and the represents the cumulative number of cars choosing tollbooth in the subsimulation, respectively. The relationships between profiles and , and and are described in equations (A.1) and (A.2), respectively:where denotes the moment in the subsimulation when the last car leaves the tollbooth , and where the superscript indicates the inverse function. In addition, we define that and .
A.1. Queue Length with respect to Time
Now, we focus on evaluating the performance measures in each period. Let denote the queue length in front of tollbooth with respect to the time in the subsimulation. Queue length in a subsimulation is the summation of the newly formed queue length and the queue length in the last subsimulation at the same time, which is described as follows:where denotes the newly formed queue length in front of tollbooth by the arriving cars during the period with respect to the time in the subsimulation. can be captured by iteratively obtaining and . In addition, is defined. The newly formed queue length in front of tollbooth by the arriving cars during the subsimulation at any time , , is given by equations (A.4), (A.5), or (A.6) with respect to different conditions as follows.
If , we have
If , we have
If , we have
A.2. Waiting Time due to Blockage of the Queues from Other Servers
Let denote the waiting time for the car due to the blockage of the queues in front of other tollbooths. For the car , the queueing blockage is presented in three scenarios as follows.
If there exists a toll lane where , such that , then the car is blocked by the queue in front of tollbooth , as Figure 19 presents. denotes the minimum queue length in front of tollbooth which blocks the path for other cars reaching their objective toll lanes. The car has to wait until the queue in front of this car for any tollbooth between toll lanes and is shorter than . In other words, the queue behind the car (the car platoon arriving at later than the car ) for any tollbooth between toll lanes and does not affect the movement of this car. Let denote the queue length for tollbooth in front of the car at time and we have
(a)
(b)
(c)
The moment when the car blocked by the queue can restart to move to the objective toll lane, , is given by . Then, is equal to .
If , then the car waits in the queue for the service and the blockage effect from other queues does not exist.
If there exists a toll lane where , such that , then the car is blocked by the queue in front of tollbooth . The car has to wait until the queue for any tollbooth between toll lanes and is shorter than . The moment when the car can restart to move, , is given by .
To sum up, the waiting time for the car due to the blockage of the queues in front of the other tollbooths is given by
A.3. Waiting Time due to Busyness of Selected Server
To help evaluate the waiting time due to the busyness of selected server, we first introduce a virtual profile for at tollbooth in the case that there is no queue in front of tollbooth . The profile is denoted by , presented in Figure 1(a) and obtained using the following equation:where represents the travel time from to tollbooth when there is no queue in front of tollbooth .
Let denote the waiting time for the car due to the busyness of tollbooth , and we have . In the case that the car arrives at the tollbooth by following the tail of the queue, the service starts at the departure time of the car . Otherwise, the service start time of the car is . Therefore, we have
is given by
B. Data for Studied Car Park
B.1. Arrival Profile
Cars’ arrival data were collected by video at 23:00–01:00 on fourteen consecutive days in November 2018. During the period, large numbers of flights arrived at the airport while subway was out of service. Many cars were expected to come and pick up passengers in the car park. The arrival samples cover from 20 to 650 cars per hour. Unfortunately, in these peak hours, no more than 650 cars exit the park per hour and the heavy-flow samples cannot be captured.
The main reason for underexpected car number is that many drivers illegally waited at the roadside lanes of the roadway around the airport rather than in the car park and then they pick up the passengers at the arrival floor in order to avoid the parking fees. When these waiting behaviors were strictly prohibited by the traffic police in future, the cars had to pour into the car park and the toll plaza at the exit would face the massive car flows. The car platoon’s arrivals at the toll plaza are represented in Figure 8. They are approximated as Poisson processes in the case of low flow via Kolmogorov–Smirnov test. Then, we assume that Poisson arrival is still fit for moderate and high flows.
Indeed, flows departing airport’s car park follow no typical stochastic process sometimes due to its high dependence on flight schedule. In addition, it is extremely hard to capture the stochastic processes with explicit mathematical forms at all levels of arrival intensities. Hence, the assumption of Poisson arrival for car platoons does not fully reflect reality; it nevertheless is an available alternative when we simulate the flow with time-varying arrival intensities.
B.2. Parking Duration and Travel Time from Car Park to
Both of the two car parks grant drivers certain time for free parking. If a car stays in the garage for less than 15 minutes, the one does not have to pay. Hence, parking duration and travel time from car park to of a car affect whether the car itself has an opportunity to be free of charge when exiting the toll plaza. The distribution of parking durations is obtained from the parking occupancy data. The data provide the times of entering and leaving the parking space for each car. In the car park, nearly 15% of cars park for less than 10 minutes which have the chance to be free-parking. In addition, under 15% of cars stay for more than 5 hours in which the longest parking car parks for over 10 days. For the majority of cars parks no more than 5 hours and full-length parking duration figure cannot show the distribution features clearly, we demonstrate a part of the parking duration distribution (less than 5 hours), as presented in Figure 4(a). The travel time distribution from car park to is estimated via the distances between all parking spaces and , as presented in Figure 4(b).
B.3. Service Time
At PAE tollbooths, free-parking cars wait for the car plates automatically verified by the cameras. The distribution of the service times for free-parking cars is presented in Figure 5(c). If the cars need to pay (i.e., PAE cars stay more than 15 minutes or POF cars stay more than 20 minutes after smartphone paying in the car park), the service time includes the time waiting for the car plates automatically verified by the cameras and the time paid at the toll plaza (i.e., QR code paid time in POF tollbooths or cash paid time in PAE tollbooths). The distributions of the service times when cars have to pay at POF and PAE tollbooths are presented in Figures 5(a) and 5(b), respectively. Usually, the QR code paid time is shorter than cash paid time, but some drivers spend far more time on QR scanning and smartphone payment than cash payment to tollmen.
Data Availability
Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request. Meanwhile, the code for the proposed model is also available.
Disclosure
The authors take sole responsibility for all views and opinions expressed in this paper.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work has been supported by the National Natural Science Foundation for Young Scholars of China (Grant no. 71901190) and Special Project for Technology Innovation and Application Development of Chongqing, China (Grant no. cstc2019jscx-tjsbX0013).