Simulation-Based Optimization for the Operation of Toll Plaza at Car Park Exit with Mixed Types of Tollbooths and Waiting-Time-Dependent Service

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.


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. e 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. e 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. erefore, 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.
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. erefore, 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. e 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. ird, the positions of POF and PAE tollbooths will affect the efficiency if the queues stretch to the bottleneck between transition area and approaching lanes. erefore, 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. us, 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. us, 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.
e traffic efficiency evaluations of POF and PAE tollbooths are derived from a mathematical model that describes the queueing process at the toll plaza. 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. e operation is described with the dynamic traffic states [1][2][3] or queue evolution [4][5][6][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][9][10][11]. e 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][19][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.

Simulation-Based Optimization Approach.
e 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. e simulation model provides the results of the optimization objective function.
is 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. e 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]. e 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]. e differences of these models are mainly in the approximating functions. e 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]. e 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. e offset between the starting times of two consecutive stages is defined as roll period. Roll periods also can be fixed [32][33][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.

Objectives and Contributions.
e 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. e 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. en, 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). e 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). e 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. e problem is decomposed by the rolling horizon approach into subproblems that are solved via KMA. e 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. e 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.

Problem Description and Model Formulation
In this study, the operation scheme of the toll plaza is determined in the time horizon [0, T], where T is the ending time of the cars' arrival at the car park exit. We discretize the time horizon [0, T] into N number of δT time periods } where T � NδT. Period n can also be denoted as [T n−1 , T n ], where T 0 � 0, T n � nδT, and T N � T. δT can be one hour in day-to-day operations. Before we formulate the optimization problem, the following assumption is proposed initially.

Assumption 1.
e 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 [0, T], 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. en, 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.

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. erefore, we need to balance the traffic efficiency and the cost for the toll plaza at the car park exit.
Let B denote the label set of tollbooths and B � 1, . . . , j, . . . , N T where N T denotes the number of tollbooths (toll lanes). e number of tollbooths is described as |B| where | · | is the cardinality of the set. 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. e 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 d [n] t and d t denote the total delays for the cars arriving in period n and entire time horizon [0, T], respectively and we have where d [n] t (B POF n , B PAE n ) represents the total delay under the scheme B POF n , B PAE n in period n and can be obtained from the nth subsimulation. Hereinafter, the superscript [n] on the notations of variables indicates that the variables are in period n (i.e., in the nth subsimulation).
We use the manpower and electricity expense to represent the operational cost. ere is a tollman for each PAE tollbooth in service. e operational cost of the toll plaza under the scheme B POF n , B PAE n in period n, c [n] p , is given by where c m and c e denote the manpower and electricity expense for each tollbooth in the unit time of the studied time horizon, respectively. (T l ) [n] denotes the moment when the last car leaves the toll plaza in the nth subsimulation. e operational cost of the toll plaza for the entire time horizon, c p , is given by Let α and α n denote the set of the indicators for tollbooth types in the entire time horizon and period n, respectively, j � 1 if tollbooth j is a POF one in service for period n. en, the sets B IN n , B POF n , and B PAE n can be derived from the set B and the variable vector α n . Hence, the tollbooth operation problem, i.e., discrete-time DTAP can be described as a simulationbased integer programming (SIMIP) model, presented as follows: [SIMIP-DTAP] where α * and α * n denote the optimal tollbooth allocation schemes for the entire time horizon and the period n, respectively. In the objective function, d t (α 1 , . . . , α N ) and c p (α 1 , . . . , α N ) are derived from equations (1) and (3), respectively. w and β are the weight for operational cost and the value of time (VOT), respectively. Note that the outputs are stochastic from simulation model. erefore, E[·] 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 j at the time δT in the (n − 1)th subsimulation, i.e., q [n−1] j (δT), which is derived from the inputs α n−1 and q [n− 2] . Here, q [n] j (·) denotes the queue length in front of tollbooth j with respect to the time in the nth subsimulation. q [n] represents the matrix of the queue length evolution for all tollbooths in period n. e queue length in our study is defined as the total number of cars staying in an area at a moment. e 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. e queue length in front of a tollbooth is the total number of cars which have selected this tollbooth at a moment. e derivation of all outputs from the simulation model such as d [n] t , (T l ) [n] , and q [n] j (·) is elaborated in the next section.

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. e customers arrive randomly with a probability distribution and can be presented as a random arrival profile. e 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. e above regulations are also captured in the simulation.

Assumptions and Simulation Inputs.
To simplify the representation of the queueing system, the following behavioral assumptions are proposed in the simulation. 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. e exit is divided by three cross sections, denoted as CS 1 , CS 2 , and CS 3 , respectively, in Figure 3. CS 1 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. CS 2 is the place where cars start to be served at the toll plaza and CS 3 is the place where cars leave the toll plaza. e studied area is determined as the one between CS 1 and CS 3 .
N A and N T denote the numbers of approaching lanes and toll lanes, respectively. ere are N L more toll lanes at the left-hand side of approaching lanes and N R more toll lanes at the right-hand side of approaching lanes. Hence, Figure 3. 6 Journal of Advanced Transportation We denote K [n] l and λ n as the number and arrival rate of car platoons in period n, respectively, and we have ? ? ?? � ? ? [?] .
e arrival pattern is generated via a given stochastic process with the arrival rate, as the representation of the cumulative profile at CS 1 for toll plaza, denoted by A [n] (t [n] ; λ n δT), t [n] ∈ [T n−1 , T n ] as the cumulative profile at CS 1 for toll plaza, where t [n] denotes the time in the nth subsimulation. Let (K l ) [n] j denote the number of cars choosing tollbooth j in the period n and we have j denote the moment in the nth subsimulation when the last car leaves the tollbooth j. In addition, we define that In conclusion, the main inputs in the simulation are the number of cars in period n, K [n] l , and the label sets of POF and PAE tollbooths in service in the period n, B POF n and B PAE n .

Driver Behaviors and Output
Evaluation. Now we focus on driver behaviors in the simulation and the output evaluation. Consider the Kth car in the period n (denoted by (K) [n] ) which arrives at CS 1 at time t [n] K with the parking duration c [n] K . e parking duration is randomly generated via empirical distribution.
Firstly, the type of the car is determined. If c [n] K < t f , the driver has the chance to be free of charge if one exits the toll plaza immediately, where t f 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 (K) [n] is PAE car. Otherwise, the driver must pay. Whether the car (K) [n] is PAE or POF car is determined using Bernoulli distribution in terms of his or her preference to the POF system. Let P r denote the proportion of drivers who prefer to the POF system and the car type is generated by P r . en, 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 ω e[n] j (K) denote the estimated waiting time for tollbooth j by the car (K) [n] . Here, we assume that drivers only know the average not the variance of a tollbooth's service times. Hence, ω e[n] j (K) is given by where s j denotes the average service time for tollbooth j. q [n] j (·) denotes the queue length in front of tollbooth j with respect to the time in the nth subsimulation. q [n] j (t [n] K ) is the summation of the newly formed queue length in the period n and the queue length in the last subsimulation (n − 1) at the

t [n]
K + δT). e detail for the derivation of q [n] j (·) is presented in Appendix A.
Once the Kth car makes a decision at CS 1 , one is labelled as the kth car choosing tollbooth j and denoted by (k) [n] j . en, 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 j start to serve it. Let τ j , ω B[n] j (k), and ω Q[n] j (k) denote the travel time from CS 1 to tollbooth j when there is no queue in front of tollbooth j, the waiting time for the car (k) [n] j due to the blockage of the queues in front of other tollbooths, and the waiting time for the car (k) [n] j due to the busyness of tollbooth j, respectively. e derivation of ω B[n] j (k) and j (k) as the waiting time for the car (k) [n] j and we have e service time for kth car is determined by the type of tollbooth j, and c [n] K , τ j , and ω [n] j (k) of the car. Let s [n] j (k) and d [n] j (k) denote the service time and delay for the car (k) [n] j , 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 n, d [n] t , is given by Total delay for the cars arriving in the entire time horizon [0, T], d t , is given by In addition, other outputs from the nth subsimulation are obtained that the next subsimulation (n + 1) needs. ey 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 n, T [n] l , is given by where for tollbooth i ∉ B IN n , we have (T l ) [n] i � 0. Matrix of the queue length evolution for all tollbooths in period n, q [n] , is given by where (T l ) [n] denotes the moment when the last car leaves the toll plaza in the nth subsimulation and Matrix of the queue length evolution for all tollbooths in the entire horizon, q, is given by where 8

Journal of Advanced Transportation
In conclusion, the main outputs from the subsimulation can be described as black-box functions, presented as follows: ; λ n δT , P r , n � 1, . . . , N, where q [0] � 0 and T [0] l � 0 are defined. e total delay for car platoons from the simulation is presented as follows: Step 0: initialization for environment. Give the sets of POF and PAE tollbooths in service, B POF n and B PAE n , respectively, the free-parking time, t f , the no-extra-pay time after POF cars paying inside the car park, t e , and the travel time from CS 1 to tollbooth j when there is no queue in front of tollbooth j, τ j . Input q [n− 1] from the last subsimulation.
Step 1: arrival information generation. Generate the arrival time t [n] K for each car (K) [n] and the total number of cars arriving, K [n] l in the period t [n] ∈ [0, δT] at CS 1 via a stochastic process with arrival intensity λ n , parking duration c [n] K via empirical distribution in Figure 4(a), and the travel time from car park to CS 1 , (t d ) [n] K via empirical distribution in Figure 4(b). Set (K) [n] � 1.
Step 2: car type generation. If c [n] K < t f , the car is labelled as PAE one. If c [n] K ≥ t f , the type of car (K) [ (9), as presented in Table 1.
Step 7: simulation termination criterion. If l , simulation terminates and give the outputs d [n] t , q [n] , and T [n] l ; return to Step 2 otherwise.

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. 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 CS 1 are estimated via the distances between all parking spaces and CS 1 . 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. e 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. e 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. e parking durations were obtained in the P7 car park for Shanghai Hongqiao International Airport. e data for service times in POF and PAE tollbooths were observed in both two car park exits. e data are presented in Appendix B.
We use queue length to validate the simulation model via hypothesis test. e 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 10 { } and 1, 2, 5 { }, respectively. e 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. e 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  Figure 5(b) Figure 5(c) Figure 5(b) POF Figure 5 POF tollbooths in service are 10 { } and 1, 2, 5 { }, respectively. e 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. e hypothesis testing for each sample (or experiment) σ (σ � 1, . . . , 28) is presented as H 0σ : μ 1σ � μ 2σ , where μ 1σ and μ 2σ denote the AQ from observation and the mean of AQs from 1000 parallel simulations, respectively. e 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 H 0σ (σ � 3, 5, 7, 12) are rejected while others are accepted. e 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.

Solution Algorithm
To solve the simulation-based integer programming model for discrete-time DTAP, SIMIP-DTAP, we propose a solution algorithm with two steps. e first step is to decompose the discrete-time DTAP into subproblems by periods where the arrival intensity and allocation scheme are static. e second one is to solve each period-based subproblem. e two steps are detailed in the next two sections, respectively.

Rolling Horizon Approach to Decomposition of Discrete-Time DTAP.
is section presents a rolling horizon (RH) approach to decompose the discrete-time DTAP, handling scheme dynamics and output stochasticity from simulation. e 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. us, 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. e 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. erefore, 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. e RH approach is presented in Figure 10, where only three stages are shown for instance. e parameter (T l ) [n] is the stage horizon which is stochastic, i.e., the time elapsed from when the period n starts to when the last car leaves the toll plaza in the nth subsimulation. e parameter r is the roll period for each stage and r � δT; e parameter T n−1 denotes the start time for the stage s (s � n). Hereinafter, we use s and n interchangeably. e computational steps of RH approach to solving discrete-time DTAP are presented as follows. [RH-DTAP] Step 0: initialization. Give the stage s � 1, subsimulation n � 1, time t � 0, optimal allocation scheme at stage (s − 1), α * n−1 , and optimal output from stage (s − 1), q [n− 1] * .
We define that q [0] * � 0 and α [0] * j � 0. c [n] p (α n ) is obtained using the following equation: Note that c [n] p (α n ) here is different from but close to that in equation (2). e solution of the subproblem STAP-n is explicated in the next section.

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. e metamodel method is favored in many literature studies due to the determinacy of the response, thus making it computationally efficient. e metamodel is an approximation that is inexpensive to compute to describe the original model with the complicated and computationally expensive response.

Framework.
e integer (global) optimization problem can be written in the following form: subject to x l j ≤ x j ≤ x u j , ∀j � 1, . . . , J, where G(x) denotes the objective function defined in equation (19). x l j and x u j denote the lower and upper bounds on the variable x j , respectively. e 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 G(x) andZ(x) is a stationary Gaussian random process with mean zero, variance σ 2 , and nonzero covariance. e covariance is expressed as where cov[Z(x m ), Z(x h )] is the covariance between two sample points x m and x h . ψ(x m , x h ) is the Kriging function, or Gaussian exponential correlation function, equivalently. e Kriging function is defined as follows: where λ � [λ 1 , . . . , λ j , . . . , λ J ] denotes the vector of unknown correlation parameters and x m j and x h j denote the jth variable in x m and x h , respectively. J denotes the dimension of x and here J � |B| in this paper. Note that the parameter λ j can be regarded as measuring the importance of the variables x m j and x h j . e larger of λ j , the greater the Euclidean distance between x m j and x h j , hence the lower the correlation between them.
Let n s denote the number of the known sample points for variable x.
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 x m , [37,38] which is rewritten as where G P (·) denotes the penalty augmented objective function (PAOF). G max represents the worst feasible objective function value so far. z( is the constraint violation function, and w p denotes the penalty factor which is positive and adjusted. is definition guarantees that the penalty augmented function values of the infeasible points are larger than the worst feasible objective function value.

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.
e Kriging metamodel algorithm (KMA) starts with an initial experimental design. e initial experimental design is to produce n s points which are repeatedly generated via a symmetric Latin hypercube design [39]. In Latin hypercube design, each variable of x is stratified into n s equal intervals, and in each subinterval, a sample point is randomly generated.
en, 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. e rth candidate point in the iteration is denoted by χ r , r � 1, . . . , R.
In the given set of randomly generated candidate points, the parameters of the Kriging metamodel are computed by the algorithm PE and the G(χ r ) is then used to approximate the G(χ r ) at every candidate point χ r . 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 V R ) and (2) minimum distance from previously evaluated points (distance criterion V D ) [40]. e weighted score is computed as where w R + w D � 1, w R ≥ 0 is the weight for the response surface criterion, and w D ≥ 0 is the weight for the distance criterion. e response surface criterion for every candidate point χ r is given by where G max � max r�1,...,R G(χ r ) and s min � min r�1,...,R G(χ r ). e distance criterion for every candidate point χ r is given by where Δ(χ r ) � min m�1,...,n s χ r − x m , Δ max � max t�1,...,T Δ(χ r ), and Δ min � min t�1,...,T Δ(χ r ). e 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.

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. e toll plaza for the car park exit has 7-lane approach and 13 tollbooths for cars where N A � 7, N R � 0, and N L � 6. We number the tollbooths from No. 1 to No. 13 consecutively (from the left-to right-hand side in the toll plaza). q B i 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. e datasets for the random parameters in the simulation are established according to Appendix B. e algorithm is coded in MATLAB and run on a personal computer with Intel (R) Core (TM) i7-8700 3.20 GHz CPU.

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. erefore, 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 10, 11 { } and 1, 2, 5, 6, 9 { }, respectively. e median, mean, and standard deviation of the key measures are obtained with different parallel evaluations for the simulations, as demonstrated in Figures 11-13.
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. us, 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.      maximum, with 20 initial points and 20 new sample points in each iteration. e 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 ese 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. e average service times for POF and PAE tollbooths are 4 and 17 seconds, respectively. e 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. erefore, 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. e toll plaza still needs more PAE tollbooths than EBR to satisfy the heavy demand even if car operators want to control operational cost.     Journal of Advanced Transportation

Appropriate Operation
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. erefore, 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.

Conclusions and Recommendations
is 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. e simulation also took into account the timevarying 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.
en, 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. e 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. e 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. erefore, 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.