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

Authors: Shanchuan Yu; Yuchuan Du; Jindong Wang; Yishun Li; Yong Zhu
Journal: Journal of Advanced Transportation
 (2021)

---

## 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.


---

## Body

## 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 Figure1(a) illustrates.Figure 1
Causes for delays at the car park exit. (a) Queueing on the toll lanes. (b) Lane changing in the transition area. (c) Bottleneck at the approaching lanes.
(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 Figure1(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. Section2 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.


## 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.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. Section2 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 horizon0,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, 1,2,…,n,…,N where T=NδT. Period n can also be denoted as Tn−1,Tn, where T0=0, Tn=nδT, and TN=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.
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 horizon0,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.InSection 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.Figure 2
Simulation for the entire time horizon and subsimulations based on periods. (a) Simulation for the entire time horizon. (b) Simulations on a period-by-period basis.
(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.LetB denote the label set of tollbooths and B=1,…,j,…,NT where NT denotes the number of tollbooths (toll lanes). The number of tollbooths is described as B where ⋅ is the cardinality of the set. BnIN denotes the label set of tollbooths in service in the period n. The set BnIN is divided into two disjoint subsets BnPOF and BnPAE, where BnPOF and BnPAE denote the label sets of POF and PAE tollbooths in service in the period n, 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. Letdtn and dt denote the total delays for the cars arriving in period n and entire time horizon 0,T, respectively and we have(1)dtB1POF,B1PAE,…,BNPOF,BNPAE=∑n=1NdtnBnPOF,BnPAE,where dtnBnPOF,BnPAE represents the total delay under the scheme BnPOF,BnPAE 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. There is a tollman for each PAE tollbooth in service. The operational cost of the toll plaza under the schemeBnPOF,BnPAE in period n, cpn, is given by(2)cpnBnPOF,BnPAE=cmBnPAEδT+ceBnPAE+BnPOFδT,if1≤n≤N−1,cmBnPAEmaxδT,Tln+ceBnPAE+BnPOFmaxδT,Tln,ifn=N,where cm and ce denote the manpower and electricity expense for each tollbooth in the unit time of the studied time horizon, respectively. Tln denotes the moment when the last car leaves the toll plaza in the nth subsimulation. The operational cost of the toll plaza for the entire time horizon, cp, is given by(3)cpB1POF,B1PAE,…,BNPOF,BNPAE=∑n=1NcpnBnPOF,BnPAE.Letα and αn denote the set of the indicators for tollbooth types in the entire time horizon and period n, respectively, where α=α1,…,αN and αn=…,αjn,…,j∈B. αjn=−1 if tollbooth j is out of service for period n. αjn=0 if tollbooth j is a PAE one in service for period n. αjn=1 if tollbooth j is a POF one in service for period n. Then, the sets BnIN, BnPOF, and BnPAE 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 simulation-based integer programming (SIMIP) model, presented as follows:[SIMIP-DTAP](4)α∗=α1∗,…,αN∗∈argminα1,…,αNpcfα=E1−wβdtα1,…,αN+wcpα1,…,αN,(5)subject toαjn=−1,0,1,j∈B,n=1,…,N,(6)∏j∈Bαjn=0,n=1,…,N,(7)Eqjn−1δT;αn−1,qn−2αjn−αjn−1=0,j∈B,n=2,…,N,where α∗ and αn∗ denote the optimal tollbooth allocation schemes for the entire time horizon and the period n, respectively. In the objective function, dtα1,…,αN and cpα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. Therefore, 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 qjn−1δT;αn−1,qn−2 denotes the queue length in front of tollbooth j at the time δT in the n−1th subsimulation, i.e., qjn−1δT, which is derived from the inputs αn−1 and qn−2.Here,qjn⋅ denotes the queue length in front of tollbooth j with respect to the time in the nth subsimulation. qn represents the matrix of the queue length evolution for all tollbooths in period n. 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 asdtn, Tln, and qjn⋅ 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 Figure3), 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.Figure 3
Queueing system for a toll plaza.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 asCS1, CS2, and CS3, respectively, in Figure 3. CS1 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. CS2 is the place where cars start to be served at the toll plaza and CS3 is the place where cars leave the toll plaza. The studied area is determined as the one between CS1 and CS3.NA and NT denote the numbers of approaching lanes and toll lanes, respectively. There are NL more toll lanes at the left-hand side of approaching lanes and NR more toll lanes at the right-hand side of approaching lanes. Hence, NT=NA+NR+NL, as presented in Figure 3.We denoteKln and λn as the number and arrival rate of car platoons in period n, 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 CS1 for toll plaza, denoted by Antn;λnδT, tn∈Tn−1,Tn as the cumulative profile at CS1 for toll plaza, where tn denotes the time in the nth subsimulation. Let Kljn denote the number of cars choosing tollbooth j in the period n and we have Kln=∑j∈BnINKljn. Let Tljn denote the moment in the nth subsimulation when the last car leaves the tollbooth j. In addition, we define that Tlj0=0 and Klj0=0.In conclusion, the main inputs in the simulation are the number of cars in periodn, Kln, and the label sets of POF and PAE tollbooths in service in the period n, BnPOF and BnPAE.


#### 2.2.2. Driver Behaviors and Output Evaluation

Now we focus on driver behaviors in the simulation and the output evaluation. Consider theKth car in the period n (denoted by Kn) which arrives at CS1 at time tKn with the parking duration cKn. The parking duration is randomly generated via empirical distribution.Firstly, the type of the car is determined. IfcKn<tf, the driver has the chance to be free of charge if one exits the toll plaza immediately, where tf 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 Kn is PAE car. Otherwise, the driver must pay. Whether the car Kn is PAE or POF car is determined using Bernoulli distribution in terms of his or her preference to the POF system. Let Pr denote the proportion of drivers who prefer to the POF system and the car type is generated by Pr.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 Assumption3). If more than one server has the shortest waiting time, the driver can select any one of those servers with identical probabilities. Let ωjenK denote the estimated waiting time for tollbooth j by the car Kn. Here, we assume that drivers only know the average not the variance of a tollbooth’s service times. Hence, ωjenK is given by(8)ωjenK=qjntKn⋅s¯j,where s¯j denotes the average service time for tollbooth j. qjn⋅ denotes the queue length in front of tollbooth j with respect to the time in the nth subsimulation. qjntKn 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 same time, qjn−1tKn+δT. The detail for the derivation of qjn⋅ is presented in Appendix A.Once theKth car makes a decision at CS1, one is labelled as the kth car choosing tollbooth j and denoted by kjn. 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 j start to serve it. Let τj, ωjBnk, and ωjQnk denote the travel time from CS1 to tollbooth j when there is no queue in front of tollbooth j, the waiting time for the car kjn due to the blockage of the queues in front of other tollbooths, and the waiting time for the car kjn due to the busyness of tollbooth j, respectively. The derivation of ωjBnk and ωjQnk is presented in Appendix A. We denote ωjnk as the waiting time for the car kjn and we have(9)ωjnk=ωjBnk+ωjQnk.The service time forkth car is determined by the type of tollbooth j, and cKn, τj, and ωjnk of the car. Let sjnk and djnk denote the service time and delay for the car kjn, respectively. We have(10)djnk=ωjnk+sjnk.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 periodn, dtn, is given by(11)dtn=∑j∈B∑k=1Kljndjnk,n=1,…,N.Total delay for the cars arriving in the entire time horizon0,T, dt, is given by(12)dt=∑n=1Ndtn.In addition, other outputs from thenth subsimulation are obtained that the next subsimulation n+1 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 n, Tln, is given by(13)Tln=TljnNT×1=Tl1n,Tl2n,…,TlNTnT,n=1,…,N,where for tollbooth i∉BnIN, we have Tlin=0. Matrix of the queue length evolution for all tollbooths in period n, qn, is given by(14)qn=qjntnNT×Tln=q1n1q1n2⋯q1nTlnq2n1q2n2⋯q2nTln⋮⋮⋱⋮qNTn1qNTn2…qNTnTln,n=1,…,N,where Tln denotes the moment when the last car leaves the toll plaza in the nth subsimulation and Tln=maxj∈BTljn. If the tollbooth i∈BnIN has Tlin<Tln, then we have qintn=0, tn∈Tlin,Tln. In addition, for tollbooth i∉BnIN, we have qintn=0, tn∈0,Tln.Matrix of the queue length evolution for all tollbooths in the entire horizon,q, is given by(15)q=q^1,…,q^n,…,q^N−1,qN,where(16)q^n=qjntnNT×δT=q1n1q1n2⋯q1nδTq2n1q2n2⋯q2nδT⋮⋮⋱⋮qNTn1qNTn2…qNTnδT,n=1,…,N−1.In conclusion, the main outputs from the subsimulation can be described as black-box functions, presented as follows:(17)dtn=dtnBnPOF,BnPAE,qn−1;Antn;λnδT,Pr,n=1,…,N,qn=qnBnPOF,BnPAE,qn−1;Antn;λnδT,Pr,n=1,…,N,Tln=TlnBnPOF,BnPAE,qn−1;Antn;λnδT,Pr,n=1,…,N,where q0=0 and Tl0=0 are defined.The total delay for car platoons from the simulation is presented as follows:(18)dt=dtB1POF,B1PAE,…,BNPOF,BNPAE;A1t1;λ1δT,…,ANtN;λNδT,Pr.


#### 2.2.3. Simulation Procedure

Thenth subsimulation (n=1,…,N) is performed as follows (Figures 4 and 5; Table 1):Figure 4
Probability densities of (a) parking durations and (b) travel times from car park toCS1.
(a)(b)Figure 5
Probability densities of service times at toll plaza. (a) Probability density of car plate automatically verified plus QR code paid time. (b) Probability density of car plate automatically verified plus cash paid time. (c) Probability density of car plate automatically time.
(a)(b)(c)Table 1
Service times in all possible scenarios.
Parking durationCar typeCondition of free-parking (for PAE car) or no-extra-pay (for POF car)Tollbooth typeService timecKn<tfPAEcKn+tdKn+τj+wjnk>tfPAEFigure5(b)cKn+tdKn+τj+wjnk≤tfFigure5(c)cKn≥tfPAE—PAEFigure5(b)POFtdKn+τj+wjnk>tePAEFigure5(b)POFFigure5(a)tdKn+τj+wjnk≤tePAE or POFFigure5(c)[SIM-n]Step 0: initialization for environment. Give the sets of POF and PAE tollbooths in service,BnPOF and BnPAE, respectively, the free-parking time, tf, the no-extra-pay time after POF cars paying inside the car park, te, and the travel time from CS1 to tollbooth j when there is no queue in front of tollbooth j, τj. Input qn−1 from the last subsimulation.Step 1: arrival information generation. Generate the arrival timetKn for each car Kn and the total number of cars arriving, Kln in the period tn∈0,δT at CS1 via a stochastic process with arrival intensity λn, parking duration cKn via empirical distribution in Figure 4(a), and the travel time from car park to CS1, tdKn via empirical distribution in Figure 4(b). Set Kn=1.Step 2: car type generation. IfcKn<tf, the car is labelled as PAE one. If cKn≥tf, the type of car Kn 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 timetKn. That is, find a PAE tollbooth label j∈I1=i∈BnPAE|ωienK≤ωlenK,∀l∈BnPAE with the probability 1/I1. 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 j∈I2=i∈BnIN|ωienK≤ωlenK,∀l∈BnIN with the probability 1/I2.Step 4: waiting time computation. Once the carKn becomes the car k chooses the jth tollbooth, i.e., kjn, we can compute the waiting time ωjBnk using equation (A.8) and the waiting time ωjQnk using equation (A.11).Step 5: service time generation. Generate the time to be served,sjnk, at tollbooth j. The service time depends on the type of tollbooths, the type of cars, and the computed waiting time ωjnk in equation (9), as presented in Table 1.Step 6: move to next arriving car. Obtain the delay of the carkjn, djnk, in equation (10). And Kn=Kn+1.Step 7: simulation termination criterion. IfKn>Kln, simulation terminates and give the outputs dtn, qn, and Tln; 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.Figure 6
Two car parks for data sources. (a) Car park for Terminal 3A of Chongqing Jiangbei International Airport. (b) P7 car park for Shanghai Hongqiao International Airport.
(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 toCS1 are estimated via the distances between all parking spaces and CS1. 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 Figure6(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 are10 and 1,2,5, 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 10 and 1,2,5, 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 Figure7. 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.Figure 7
Comparison of queue lengths at every second in the evolution between 1000 simulations and observation in one case.Figure 8
Observed arrival profiles for low-flow case (the lines with different colors represent the different samples of arrivals we observed).The hypothesis testing for each sample (or experiment)σ (σ=1,…,28) is presented as H0σ:μ1σ=μ2σ, where μ1σ and μ2σ 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.Figure 9
Comparison of queue lengths at every second in the evolution between simulations and observation in one case.Table 2
Significances and 0.95 confidence intervals for validation samples.
Experiment no.1234567891011121314H0σ00101010000100Significance0.670.430.000.900.000.220.020.280.220.690.420.000.160.93LB of CI2.942.973.843.955.785.936.906.936.987.968.978.7110.9211.96UB of CI3.043.083.934.045.856.026.997.027.078.059.068.7811.0112.05Experiment no.1516171819202122232425262728H0σ00000000000000Significance0.970.240.470.540.940.620.200.630.360.510.060.090.790.30LB of CI11.9512.9814.9416.9417.9519.9720.9121.9722.9423.9425.9125.9131.9331.90UB of CI12.0513.0715.0317.0318.0620.0521.0222.0623.0224.0326.0026.0132.0532.03According to Table2, only four hypotheses H0σ (σ=3,5,7,12) 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.


## 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.LetB denote the label set of tollbooths and B=1,…,j,…,NT where NT denotes the number of tollbooths (toll lanes). The number of tollbooths is described as B where ⋅ is the cardinality of the set. BnIN denotes the label set of tollbooths in service in the period n. The set BnIN is divided into two disjoint subsets BnPOF and BnPAE, where BnPOF and BnPAE denote the label sets of POF and PAE tollbooths in service in the period n, 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. Letdtn and dt denote the total delays for the cars arriving in period n and entire time horizon 0,T, respectively and we have(1)dtB1POF,B1PAE,…,BNPOF,BNPAE=∑n=1NdtnBnPOF,BnPAE,where dtnBnPOF,BnPAE represents the total delay under the scheme BnPOF,BnPAE 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. There is a tollman for each PAE tollbooth in service. The operational cost of the toll plaza under the schemeBnPOF,BnPAE in period n, cpn, is given by(2)cpnBnPOF,BnPAE=cmBnPAEδT+ceBnPAE+BnPOFδT,if1≤n≤N−1,cmBnPAEmaxδT,Tln+ceBnPAE+BnPOFmaxδT,Tln,ifn=N,where cm and ce denote the manpower and electricity expense for each tollbooth in the unit time of the studied time horizon, respectively. Tln denotes the moment when the last car leaves the toll plaza in the nth subsimulation. The operational cost of the toll plaza for the entire time horizon, cp, is given by(3)cpB1POF,B1PAE,…,BNPOF,BNPAE=∑n=1NcpnBnPOF,BnPAE.Letα and αn denote the set of the indicators for tollbooth types in the entire time horizon and period n, respectively, where α=α1,…,αN and αn=…,αjn,…,j∈B. αjn=−1 if tollbooth j is out of service for period n. αjn=0 if tollbooth j is a PAE one in service for period n. αjn=1 if tollbooth j is a POF one in service for period n. Then, the sets BnIN, BnPOF, and BnPAE 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 simulation-based integer programming (SIMIP) model, presented as follows:[SIMIP-DTAP](4)α∗=α1∗,…,αN∗∈argminα1,…,αNpcfα=E1−wβdtα1,…,αN+wcpα1,…,αN,(5)subject toαjn=−1,0,1,j∈B,n=1,…,N,(6)∏j∈Bαjn=0,n=1,…,N,(7)Eqjn−1δT;αn−1,qn−2αjn−αjn−1=0,j∈B,n=2,…,N,where α∗ and αn∗ denote the optimal tollbooth allocation schemes for the entire time horizon and the period n, respectively. In the objective function, dtα1,…,αN and cpα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. Therefore, 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 qjn−1δT;αn−1,qn−2 denotes the queue length in front of tollbooth j at the time δT in the n−1th subsimulation, i.e., qjn−1δT, which is derived from the inputs αn−1 and qn−2.Here,qjn⋅ denotes the queue length in front of tollbooth j with respect to the time in the nth subsimulation. qn represents the matrix of the queue length evolution for all tollbooths in period n. 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 asdtn, Tln, and qjn⋅ 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 Figure3), 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.Figure 3
Queueing system for a toll plaza.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 asCS1, CS2, and CS3, respectively, in Figure 3. CS1 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. CS2 is the place where cars start to be served at the toll plaza and CS3 is the place where cars leave the toll plaza. The studied area is determined as the one between CS1 and CS3.NA and NT denote the numbers of approaching lanes and toll lanes, respectively. There are NL more toll lanes at the left-hand side of approaching lanes and NR more toll lanes at the right-hand side of approaching lanes. Hence, NT=NA+NR+NL, as presented in Figure 3.We denoteKln and λn as the number and arrival rate of car platoons in period n, 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 CS1 for toll plaza, denoted by Antn;λnδT, tn∈Tn−1,Tn as the cumulative profile at CS1 for toll plaza, where tn denotes the time in the nth subsimulation. Let Kljn denote the number of cars choosing tollbooth j in the period n and we have Kln=∑j∈BnINKljn. Let Tljn denote the moment in the nth subsimulation when the last car leaves the tollbooth j. In addition, we define that Tlj0=0 and Klj0=0.In conclusion, the main inputs in the simulation are the number of cars in periodn, Kln, and the label sets of POF and PAE tollbooths in service in the period n, BnPOF and BnPAE.


### 2.2.2. Driver Behaviors and Output Evaluation

Now we focus on driver behaviors in the simulation and the output evaluation. Consider theKth car in the period n (denoted by Kn) which arrives at CS1 at time tKn with the parking duration cKn. The parking duration is randomly generated via empirical distribution.Firstly, the type of the car is determined. IfcKn<tf, the driver has the chance to be free of charge if one exits the toll plaza immediately, where tf 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 Kn is PAE car. Otherwise, the driver must pay. Whether the car Kn is PAE or POF car is determined using Bernoulli distribution in terms of his or her preference to the POF system. Let Pr denote the proportion of drivers who prefer to the POF system and the car type is generated by Pr.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 Assumption3). If more than one server has the shortest waiting time, the driver can select any one of those servers with identical probabilities. Let ωjenK denote the estimated waiting time for tollbooth j by the car Kn. Here, we assume that drivers only know the average not the variance of a tollbooth’s service times. Hence, ωjenK is given by(8)ωjenK=qjntKn⋅s¯j,where s¯j denotes the average service time for tollbooth j. qjn⋅ denotes the queue length in front of tollbooth j with respect to the time in the nth subsimulation. qjntKn 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 same time, qjn−1tKn+δT. The detail for the derivation of qjn⋅ is presented in Appendix A.Once theKth car makes a decision at CS1, one is labelled as the kth car choosing tollbooth j and denoted by kjn. 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 j start to serve it. Let τj, ωjBnk, and ωjQnk denote the travel time from CS1 to tollbooth j when there is no queue in front of tollbooth j, the waiting time for the car kjn due to the blockage of the queues in front of other tollbooths, and the waiting time for the car kjn due to the busyness of tollbooth j, respectively. The derivation of ωjBnk and ωjQnk is presented in Appendix A. We denote ωjnk as the waiting time for the car kjn and we have(9)ωjnk=ωjBnk+ωjQnk.The service time forkth car is determined by the type of tollbooth j, and cKn, τj, and ωjnk of the car. Let sjnk and djnk denote the service time and delay for the car kjn, respectively. We have(10)djnk=ωjnk+sjnk.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 periodn, dtn, is given by(11)dtn=∑j∈B∑k=1Kljndjnk,n=1,…,N.Total delay for the cars arriving in the entire time horizon0,T, dt, is given by(12)dt=∑n=1Ndtn.In addition, other outputs from thenth subsimulation are obtained that the next subsimulation n+1 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 n, Tln, is given by(13)Tln=TljnNT×1=Tl1n,Tl2n,…,TlNTnT,n=1,…,N,where for tollbooth i∉BnIN, we have Tlin=0. Matrix of the queue length evolution for all tollbooths in period n, qn, is given by(14)qn=qjntnNT×Tln=q1n1q1n2⋯q1nTlnq2n1q2n2⋯q2nTln⋮⋮⋱⋮qNTn1qNTn2…qNTnTln,n=1,…,N,where Tln denotes the moment when the last car leaves the toll plaza in the nth subsimulation and Tln=maxj∈BTljn. If the tollbooth i∈BnIN has Tlin<Tln, then we have qintn=0, tn∈Tlin,Tln. In addition, for tollbooth i∉BnIN, we have qintn=0, tn∈0,Tln.Matrix of the queue length evolution for all tollbooths in the entire horizon,q, is given by(15)q=q^1,…,q^n,…,q^N−1,qN,where(16)q^n=qjntnNT×δT=q1n1q1n2⋯q1nδTq2n1q2n2⋯q2nδT⋮⋮⋱⋮qNTn1qNTn2…qNTnδT,n=1,…,N−1.In conclusion, the main outputs from the subsimulation can be described as black-box functions, presented as follows:(17)dtn=dtnBnPOF,BnPAE,qn−1;Antn;λnδT,Pr,n=1,…,N,qn=qnBnPOF,BnPAE,qn−1;Antn;λnδT,Pr,n=1,…,N,Tln=TlnBnPOF,BnPAE,qn−1;Antn;λnδT,Pr,n=1,…,N,where q0=0 and Tl0=0 are defined.The total delay for car platoons from the simulation is presented as follows:(18)dt=dtB1POF,B1PAE,…,BNPOF,BNPAE;A1t1;λ1δT,…,ANtN;λNδT,Pr.


### 2.2.3. Simulation Procedure

Thenth subsimulation (n=1,…,N) is performed as follows (Figures 4 and 5; Table 1):Figure 4
Probability densities of (a) parking durations and (b) travel times from car park toCS1.
(a)(b)Figure 5
Probability densities of service times at toll plaza. (a) Probability density of car plate automatically verified plus QR code paid time. (b) Probability density of car plate automatically verified plus cash paid time. (c) Probability density of car plate automatically time.
(a)(b)(c)Table 1
Service times in all possible scenarios.
Parking durationCar typeCondition of free-parking (for PAE car) or no-extra-pay (for POF car)Tollbooth typeService timecKn<tfPAEcKn+tdKn+τj+wjnk>tfPAEFigure5(b)cKn+tdKn+τj+wjnk≤tfFigure5(c)cKn≥tfPAE—PAEFigure5(b)POFtdKn+τj+wjnk>tePAEFigure5(b)POFFigure5(a)tdKn+τj+wjnk≤tePAE or POFFigure5(c)[SIM-n]Step 0: initialization for environment. Give the sets of POF and PAE tollbooths in service,BnPOF and BnPAE, respectively, the free-parking time, tf, the no-extra-pay time after POF cars paying inside the car park, te, and the travel time from CS1 to tollbooth j when there is no queue in front of tollbooth j, τj. Input qn−1 from the last subsimulation.Step 1: arrival information generation. Generate the arrival timetKn for each car Kn and the total number of cars arriving, Kln in the period tn∈0,δT at CS1 via a stochastic process with arrival intensity λn, parking duration cKn via empirical distribution in Figure 4(a), and the travel time from car park to CS1, tdKn via empirical distribution in Figure 4(b). Set Kn=1.Step 2: car type generation. IfcKn<tf, the car is labelled as PAE one. If cKn≥tf, the type of car Kn 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 timetKn. That is, find a PAE tollbooth label j∈I1=i∈BnPAE|ωienK≤ωlenK,∀l∈BnPAE with the probability 1/I1. 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 j∈I2=i∈BnIN|ωienK≤ωlenK,∀l∈BnIN with the probability 1/I2.Step 4: waiting time computation. Once the carKn becomes the car k chooses the jth tollbooth, i.e., kjn, we can compute the waiting time ωjBnk using equation (A.8) and the waiting time ωjQnk using equation (A.11).Step 5: service time generation. Generate the time to be served,sjnk, at tollbooth j. The service time depends on the type of tollbooths, the type of cars, and the computed waiting time ωjnk in equation (9), as presented in Table 1.Step 6: move to next arriving car. Obtain the delay of the carkjn, djnk, in equation (10). And Kn=Kn+1.Step 7: simulation termination criterion. IfKn>Kln, simulation terminates and give the outputs dtn, qn, and Tln; 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.Figure 6
Two car parks for data sources. (a) Car park for Terminal 3A of Chongqing Jiangbei International Airport. (b) P7 car park for Shanghai Hongqiao International Airport.
(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 toCS1 are estimated via the distances between all parking spaces and CS1. 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 Figure6(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 are10 and 1,2,5, 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 10 and 1,2,5, 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 Figure7. 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.Figure 7
Comparison of queue lengths at every second in the evolution between 1000 simulations and observation in one case.Figure 8
Observed arrival profiles for low-flow case (the lines with different colors represent the different samples of arrivals we observed).The hypothesis testing for each sample (or experiment)σ (σ=1,…,28) is presented as H0σ:μ1σ=μ2σ, where μ1σ and μ2σ 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.Figure 9
Comparison of queue lengths at every second in the evolution between simulations and observation in one case.Table 2
Significances and 0.95 confidence intervals for validation samples.
Experiment no.1234567891011121314H0σ00101010000100Significance0.670.430.000.900.000.220.020.280.220.690.420.000.160.93LB of CI2.942.973.843.955.785.936.906.936.987.968.978.7110.9211.96UB of CI3.043.083.934.045.856.026.997.027.078.059.068.7811.0112.05Experiment no.1516171819202122232425262728H0σ00000000000000Significance0.970.240.470.540.940.620.200.630.360.510.060.090.790.30LB of CI11.9512.9814.9416.9417.9519.9720.9121.9722.9423.9425.9125.9131.9331.90UB of CI12.0513.0715.0317.0318.0620.0521.0222.0623.0224.0326.0026.0132.0532.03According to Table2, only four hypotheses H0σ (σ=3,5,7,12) 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.


## 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 asCS1, CS2, and CS3, respectively, in Figure 3. CS1 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. CS2 is the place where cars start to be served at the toll plaza and CS3 is the place where cars leave the toll plaza. The studied area is determined as the one between CS1 and CS3.NA and NT denote the numbers of approaching lanes and toll lanes, respectively. There are NL more toll lanes at the left-hand side of approaching lanes and NR more toll lanes at the right-hand side of approaching lanes. Hence, NT=NA+NR+NL, as presented in Figure 3.We denoteKln and λn as the number and arrival rate of car platoons in period n, 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 CS1 for toll plaza, denoted by Antn;λnδT, tn∈Tn−1,Tn as the cumulative profile at CS1 for toll plaza, where tn denotes the time in the nth subsimulation. Let Kljn denote the number of cars choosing tollbooth j in the period n and we have Kln=∑j∈BnINKljn. Let Tljn denote the moment in the nth subsimulation when the last car leaves the tollbooth j. In addition, we define that Tlj0=0 and Klj0=0.In conclusion, the main inputs in the simulation are the number of cars in periodn, Kln, and the label sets of POF and PAE tollbooths in service in the period n, BnPOF and BnPAE.


## 2.2.2. Driver Behaviors and Output Evaluation

Now we focus on driver behaviors in the simulation and the output evaluation. Consider theKth car in the period n (denoted by Kn) which arrives at CS1 at time tKn with the parking duration cKn. The parking duration is randomly generated via empirical distribution.Firstly, the type of the car is determined. IfcKn<tf, the driver has the chance to be free of charge if one exits the toll plaza immediately, where tf 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 Kn is PAE car. Otherwise, the driver must pay. Whether the car Kn is PAE or POF car is determined using Bernoulli distribution in terms of his or her preference to the POF system. Let Pr denote the proportion of drivers who prefer to the POF system and the car type is generated by Pr.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 Assumption3). If more than one server has the shortest waiting time, the driver can select any one of those servers with identical probabilities. Let ωjenK denote the estimated waiting time for tollbooth j by the car Kn. Here, we assume that drivers only know the average not the variance of a tollbooth’s service times. Hence, ωjenK is given by(8)ωjenK=qjntKn⋅s¯j,where s¯j denotes the average service time for tollbooth j. qjn⋅ denotes the queue length in front of tollbooth j with respect to the time in the nth subsimulation. qjntKn 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 same time, qjn−1tKn+δT. The detail for the derivation of qjn⋅ is presented in Appendix A.Once theKth car makes a decision at CS1, one is labelled as the kth car choosing tollbooth j and denoted by kjn. 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 j start to serve it. Let τj, ωjBnk, and ωjQnk denote the travel time from CS1 to tollbooth j when there is no queue in front of tollbooth j, the waiting time for the car kjn due to the blockage of the queues in front of other tollbooths, and the waiting time for the car kjn due to the busyness of tollbooth j, respectively. The derivation of ωjBnk and ωjQnk is presented in Appendix A. We denote ωjnk as the waiting time for the car kjn and we have(9)ωjnk=ωjBnk+ωjQnk.The service time forkth car is determined by the type of tollbooth j, and cKn, τj, and ωjnk of the car. Let sjnk and djnk denote the service time and delay for the car kjn, respectively. We have(10)djnk=ωjnk+sjnk.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 periodn, dtn, is given by(11)dtn=∑j∈B∑k=1Kljndjnk,n=1,…,N.Total delay for the cars arriving in the entire time horizon0,T, dt, is given by(12)dt=∑n=1Ndtn.In addition, other outputs from thenth subsimulation are obtained that the next subsimulation n+1 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 n, Tln, is given by(13)Tln=TljnNT×1=Tl1n,Tl2n,…,TlNTnT,n=1,…,N,where for tollbooth i∉BnIN, we have Tlin=0. Matrix of the queue length evolution for all tollbooths in period n, qn, is given by(14)qn=qjntnNT×Tln=q1n1q1n2⋯q1nTlnq2n1q2n2⋯q2nTln⋮⋮⋱⋮qNTn1qNTn2…qNTnTln,n=1,…,N,where Tln denotes the moment when the last car leaves the toll plaza in the nth subsimulation and Tln=maxj∈BTljn. If the tollbooth i∈BnIN has Tlin<Tln, then we have qintn=0, tn∈Tlin,Tln. In addition, for tollbooth i∉BnIN, we have qintn=0, tn∈0,Tln.Matrix of the queue length evolution for all tollbooths in the entire horizon,q, is given by(15)q=q^1,…,q^n,…,q^N−1,qN,where(16)q^n=qjntnNT×δT=q1n1q1n2⋯q1nδTq2n1q2n2⋯q2nδT⋮⋮⋱⋮qNTn1qNTn2…qNTnδT,n=1,…,N−1.In conclusion, the main outputs from the subsimulation can be described as black-box functions, presented as follows:(17)dtn=dtnBnPOF,BnPAE,qn−1;Antn;λnδT,Pr,n=1,…,N,qn=qnBnPOF,BnPAE,qn−1;Antn;λnδT,Pr,n=1,…,N,Tln=TlnBnPOF,BnPAE,qn−1;Antn;λnδT,Pr,n=1,…,N,where q0=0 and Tl0=0 are defined.The total delay for car platoons from the simulation is presented as follows:(18)dt=dtB1POF,B1PAE,…,BNPOF,BNPAE;A1t1;λ1δT,…,ANtN;λNδT,Pr.


## 2.2.3. Simulation Procedure

Thenth subsimulation (n=1,…,N) is performed as follows (Figures 4 and 5; Table 1):Figure 4
Probability densities of (a) parking durations and (b) travel times from car park toCS1.
(a)(b)Figure 5
Probability densities of service times at toll plaza. (a) Probability density of car plate automatically verified plus QR code paid time. (b) Probability density of car plate automatically verified plus cash paid time. (c) Probability density of car plate automatically time.
(a)(b)(c)Table 1
Service times in all possible scenarios.
Parking durationCar typeCondition of free-parking (for PAE car) or no-extra-pay (for POF car)Tollbooth typeService timecKn<tfPAEcKn+tdKn+τj+wjnk>tfPAEFigure5(b)cKn+tdKn+τj+wjnk≤tfFigure5(c)cKn≥tfPAE—PAEFigure5(b)POFtdKn+τj+wjnk>tePAEFigure5(b)POFFigure5(a)tdKn+τj+wjnk≤tePAE or POFFigure5(c)[SIM-n]Step 0: initialization for environment. Give the sets of POF and PAE tollbooths in service,BnPOF and BnPAE, respectively, the free-parking time, tf, the no-extra-pay time after POF cars paying inside the car park, te, and the travel time from CS1 to tollbooth j when there is no queue in front of tollbooth j, τj. Input qn−1 from the last subsimulation.Step 1: arrival information generation. Generate the arrival timetKn for each car Kn and the total number of cars arriving, Kln in the period tn∈0,δT at CS1 via a stochastic process with arrival intensity λn, parking duration cKn via empirical distribution in Figure 4(a), and the travel time from car park to CS1, tdKn via empirical distribution in Figure 4(b). Set Kn=1.Step 2: car type generation. IfcKn<tf, the car is labelled as PAE one. If cKn≥tf, the type of car Kn 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 timetKn. That is, find a PAE tollbooth label j∈I1=i∈BnPAE|ωienK≤ωlenK,∀l∈BnPAE with the probability 1/I1. 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 j∈I2=i∈BnIN|ωienK≤ωlenK,∀l∈BnIN with the probability 1/I2.Step 4: waiting time computation. Once the carKn becomes the car k chooses the jth tollbooth, i.e., kjn, we can compute the waiting time ωjBnk using equation (A.8) and the waiting time ωjQnk using equation (A.11).Step 5: service time generation. Generate the time to be served,sjnk, at tollbooth j. The service time depends on the type of tollbooths, the type of cars, and the computed waiting time ωjnk in equation (9), as presented in Table 1.Step 6: move to next arriving car. Obtain the delay of the carkjn, djnk, in equation (10). And Kn=Kn+1.Step 7: simulation termination criterion. IfKn>Kln, simulation terminates and give the outputs dtn, qn, and Tln; 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.Figure 6
Two car parks for data sources. (a) Car park for Terminal 3A of Chongqing Jiangbei International Airport. (b) P7 car park for Shanghai Hongqiao International Airport.
(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 toCS1 are estimated via the distances between all parking spaces and CS1. 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 Figure6(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 are10 and 1,2,5, 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 10 and 1,2,5, 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 Figure7. 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.Figure 7
Comparison of queue lengths at every second in the evolution between 1000 simulations and observation in one case.Figure 8
Observed arrival profiles for low-flow case (the lines with different colors represent the different samples of arrivals we observed).The hypothesis testing for each sample (or experiment)σ (σ=1,…,28) is presented as H0σ:μ1σ=μ2σ, where μ1σ and μ2σ 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.Figure 9
Comparison of queue lengths at every second in the evolution between simulations and observation in one case.Table 2
Significances and 0.95 confidence intervals for validation samples.
Experiment no.1234567891011121314H0σ00101010000100Significance0.670.430.000.900.000.220.020.280.220.690.420.000.160.93LB of CI2.942.973.843.955.785.936.906.936.987.968.978.7110.9211.96UB of CI3.043.083.934.045.856.026.997.027.078.059.068.7811.0112.05Experiment no.1516171819202122232425262728H0σ00000000000000Significance0.970.240.470.540.940.620.200.630.360.510.060.090.790.30LB of CI11.9512.9814.9416.9417.9519.9720.9121.9722.9423.9425.9125.9131.9331.90UB of CI12.0513.0715.0317.0318.0620.0521.0222.0623.0224.0326.0026.0132.0532.03According to Table2, only four hypotheses H0σ (σ=3,5,7,12) 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 Assumption1, 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 Tln 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. The parameter r is the roll period for each stage and r=δT; The parameter Tn−1 denotes the start time for the stage s (s=n). Hereinafter, we use s and n interchangeably.Figure 10
The rolling horizon approach.The computational steps of RH approach to solving discrete-time DTAP are presented as follows.[RH-DTAP]Step 0: initialization. Give the stages=1, subsimulation n=1, time t=0, optimal allocation scheme at stage s−1, αn−1∗, and optimal output from stage s−1, qn−1∗.Step 1: solve thenth static tollbooth allocation problem (STAP-n) based on inputs αn−1∗ and qn−1∗, and nth subsimulation at stage s. Obtain the optimal allocation scheme αn∗ and output qn∗.Step 2: update the stage, subsimulation, and time,s=s+1, n=n+1, and t=t+δT.Step 3: iteration termination criterion. Ifs>N, iteration terminates and give the optimal allocation scheme sequence α∗=α1∗,…,αN∗; return to Step 1 otherwise.In RH, the subproblem STAP-n (n=1,…,N) is also described as a SIMIP model, presented as follows:[SIMIP-STAP-n](19)αn∗∈argminαnpcfnαn=E1−wβdtnαn;qn−1∗+wcpnαn,(20)subject toαjn=−1,0,1,j∈B,(21)∏j∈Bαjn=0,(22)Eqjn−1∗δTαjn−αjn−1∗=0,j∈B,where qn−1∗ and αjn−1∗ are the optimal results from [SIMIP-STAP-n−1] as well as the constants in [STAP-n]. We define that q0∗=0 and αj0∗=0. cpnαn is obtained using the following equation:(23)cpnαn=cmBnPAEmaxδT,Tln+ceBnPAE+BnPOFmaxδT,Tln,n=1,…,N.Note thatcpnαn 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:(24)minimizeGx,(25)subject toxjl≤xj≤xju,∀j=1,…,J,where Gx denotes the objective function defined in equation (19). xjl and xju denote the lower and upper bounds on the variable xj, respectively. The optimization model defined in (24) is replaced by the sum of a constant value and a Gaussian random error term as follows:(26)Gx=ρ+Zx,where ρ is the mean value of Gx andZx is a stationary Gaussian random process with mean zero, variance σ2, and nonzero covariance. The covariance is expressed as(27)covZxm,Zxh=σ2ψxm,xh,where covZxm,Zxh is the covariance between two sample points xm and xh. ψxm,xh is the Kriging function, or Gaussian exponential correlation function, equivalently. The Kriging function is defined as follows:(28)ψxm,xh=exp−∑j=1Jλjxjm−xjh2,where λ=λ1,…,λj,…,λJ denotes the vector of unknown correlation parameters and xjm and xjh denote the jth variable in xm and xh, 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 xjm and xjh. The larger of λj, the greater the Euclidean distance between xjm and xjh, hence the lower the correlation between them.Letns denote the number of the known sample points for variable x. The predictive estimate G˜x of the response function Gx at unknown point x is given by(29)G˜x=ρ˜+R˜TxΨ˜−1G−1ρ˜,where G=Gx1,Gx2,…,GxnsT, 1 is the ns×1 unit column, and Ψ˜ is the estimated correlation matrix and is diagonally symmetric where(30)Ψ˜=ψ˜x1,x1ψ˜x1,x2⋯ψ˜x1,xnsψ˜x2,x1ψ˜x2,x2⋯ψ˜x2,xns⋮⋮⋱⋮ψ˜xns,x1ψ˜xns,x2⋯ψ˜xns,xns∈ℛns×ns,where(31)ψ˜xm,xh=exp−∑j=1Jλ˜jxjm−xjh2.R˜x is the correlation vector between the unknown point x and sample points where(32)R˜x=ψ˜x,x1,ψ˜x,x2,…,ψ˜x,xnsT.Therefore, in order to yield the value ofG˜x, the estimators ρ˜ and λ˜1,λ˜2,…,λ˜K, λ˜2, λ˜K need to be obtained, given by(33)ρ˜=1TΨ˜−11−11TΨ˜−1G,(34)λ˜1,λ˜2,…,λ˜J=argmaxλ˜1,λ˜2,…,λ˜J>0−12nsln1nsG−1ρ˜TΨ˜−1G−1ρ˜+lndetΨ˜,where det is the determinant of matrix [35]. According to equation (33), ρ˜ is obtained based on λ˜1,λ˜2,…,λ˜J, λ˜1, λ˜2, λ˜J while λ˜1,λ˜2,…,λ˜J, λ˜1, λ˜2, λ˜J 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λ˜1,λ˜2,…,λ˜K, λ˜1, λ˜2, λ˜K, and obtain the ψ˜xm,xh 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 λ˜1,λ˜2,…,λ˜K, λ˜1, λ˜2, λ˜KStep 3: updateψ˜xm,xh and Ψ˜ with the updated λ˜1,λ˜2,…,λ˜J, λ˜1, λ˜2, λ˜J, 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 λ˜1,λ˜2,…,λ˜J, λ˜1, λ˜2, λ˜JIf the problem has other constraints in addition to the bound constraint (25), it can be formulated as the following form:(35)minimizeGx,subject to(36)Hjx≤0,∀j=1,…,J,(37)xjl≤xj≤xju,∀j=1,…,J.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 xm, [37, 38] which is rewritten as(38)GPxm=Gmax+wpzxm,ifzxm>0,Gxm,otherwise,where GP⋅ denotes the penalty augmented objective function (PAOF). Gmax represents the worst feasible objective function value so far. zxk=∑j=1mmax0,Hjxk2 is the constraint violation function, and wp 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 Section3.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 producens points which are repeatedly generated via a symmetric Latin hypercube design [39]. In Latin hypercube design, each variable of x is stratified into ns 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. Therth 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 theG˜χ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 VR) and (2) minimum distance from previously evaluated points (distance criterion VD) [40]. The weighted score is computed as(39)Wr=wRVRχr+wDVDχr,where wR+wD=1, wR≥0 is the weight for the response surface criterion, and wD≥0 is the weight for the distance criterion. The response surface criterion for every candidate point χr is given by(40)VRχr=G˜χr−GminGmax−Gmin,ifGmax≠Gmin,1,otherwise,where Gmax=maxr=1,…,RG˜χr and smin=minr=1,…,RG˜χr. The distance criterion for every candidate point χr is given by(41)VDχr=Δmax−ΔχrΔmax−Δmin,ifΔmax≠Δmin,1,otherwise,where Δχr=minm=1,…,nsχr−xm, Δmax=maxt=1,…,TΔχr, and Δmin=mint=1,…,TΔχr.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.


## 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 Assumption1, 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 Tln 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. The parameter r is the roll period for each stage and r=δT; The parameter Tn−1 denotes the start time for the stage s (s=n). Hereinafter, we use s and n interchangeably.Figure 10
The rolling horizon approach.The computational steps of RH approach to solving discrete-time DTAP are presented as follows.[RH-DTAP]Step 0: initialization. Give the stages=1, subsimulation n=1, time t=0, optimal allocation scheme at stage s−1, αn−1∗, and optimal output from stage s−1, qn−1∗.Step 1: solve thenth static tollbooth allocation problem (STAP-n) based on inputs αn−1∗ and qn−1∗, and nth subsimulation at stage s. Obtain the optimal allocation scheme αn∗ and output qn∗.Step 2: update the stage, subsimulation, and time,s=s+1, n=n+1, and t=t+δT.Step 3: iteration termination criterion. Ifs>N, iteration terminates and give the optimal allocation scheme sequence α∗=α1∗,…,αN∗; return to Step 1 otherwise.In RH, the subproblem STAP-n (n=1,…,N) is also described as a SIMIP model, presented as follows:[SIMIP-STAP-n](19)αn∗∈argminαnpcfnαn=E1−wβdtnαn;qn−1∗+wcpnαn,(20)subject toαjn=−1,0,1,j∈B,(21)∏j∈Bαjn=0,(22)Eqjn−1∗δTαjn−αjn−1∗=0,j∈B,where qn−1∗ and αjn−1∗ are the optimal results from [SIMIP-STAP-n−1] as well as the constants in [STAP-n]. We define that q0∗=0 and αj0∗=0. cpnαn is obtained using the following equation:(23)cpnαn=cmBnPAEmaxδT,Tln+ceBnPAE+BnPOFmaxδT,Tln,n=1,…,N.Note thatcpnαn 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:(24)minimizeGx,(25)subject toxjl≤xj≤xju,∀j=1,…,J,where Gx denotes the objective function defined in equation (19). xjl and xju denote the lower and upper bounds on the variable xj, respectively. The optimization model defined in (24) is replaced by the sum of a constant value and a Gaussian random error term as follows:(26)Gx=ρ+Zx,where ρ is the mean value of Gx andZx is a stationary Gaussian random process with mean zero, variance σ2, and nonzero covariance. The covariance is expressed as(27)covZxm,Zxh=σ2ψxm,xh,where covZxm,Zxh is the covariance between two sample points xm and xh. ψxm,xh is the Kriging function, or Gaussian exponential correlation function, equivalently. The Kriging function is defined as follows:(28)ψxm,xh=exp−∑j=1Jλjxjm−xjh2,where λ=λ1,…,λj,…,λJ denotes the vector of unknown correlation parameters and xjm and xjh denote the jth variable in xm and xh, 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 xjm and xjh. The larger of λj, the greater the Euclidean distance between xjm and xjh, hence the lower the correlation between them.Letns denote the number of the known sample points for variable x. The predictive estimate G˜x of the response function Gx at unknown point x is given by(29)G˜x=ρ˜+R˜TxΨ˜−1G−1ρ˜,where G=Gx1,Gx2,…,GxnsT, 1 is the ns×1 unit column, and Ψ˜ is the estimated correlation matrix and is diagonally symmetric where(30)Ψ˜=ψ˜x1,x1ψ˜x1,x2⋯ψ˜x1,xnsψ˜x2,x1ψ˜x2,x2⋯ψ˜x2,xns⋮⋮⋱⋮ψ˜xns,x1ψ˜xns,x2⋯ψ˜xns,xns∈ℛns×ns,where(31)ψ˜xm,xh=exp−∑j=1Jλ˜jxjm−xjh2.R˜x is the correlation vector between the unknown point x and sample points where(32)R˜x=ψ˜x,x1,ψ˜x,x2,…,ψ˜x,xnsT.Therefore, in order to yield the value ofG˜x, the estimators ρ˜ and λ˜1,λ˜2,…,λ˜K, λ˜2, λ˜K need to be obtained, given by(33)ρ˜=1TΨ˜−11−11TΨ˜−1G,(34)λ˜1,λ˜2,…,λ˜J=argmaxλ˜1,λ˜2,…,λ˜J>0−12nsln1nsG−1ρ˜TΨ˜−1G−1ρ˜+lndetΨ˜,where det is the determinant of matrix [35]. According to equation (33), ρ˜ is obtained based on λ˜1,λ˜2,…,λ˜J, λ˜1, λ˜2, λ˜J while λ˜1,λ˜2,…,λ˜J, λ˜1, λ˜2, λ˜J 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λ˜1,λ˜2,…,λ˜K, λ˜1, λ˜2, λ˜K, and obtain the ψ˜xm,xh 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 λ˜1,λ˜2,…,λ˜K, λ˜1, λ˜2, λ˜KStep 3: updateψ˜xm,xh and Ψ˜ with the updated λ˜1,λ˜2,…,λ˜J, λ˜1, λ˜2, λ˜J, 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 λ˜1,λ˜2,…,λ˜J, λ˜1, λ˜2, λ˜JIf the problem has other constraints in addition to the bound constraint (25), it can be formulated as the following form:(35)minimizeGx,subject to(36)Hjx≤0,∀j=1,…,J,(37)xjl≤xj≤xju,∀j=1,…,J.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 xm, [37, 38] which is rewritten as(38)GPxm=Gmax+wpzxm,ifzxm>0,Gxm,otherwise,where GP⋅ denotes the penalty augmented objective function (PAOF). Gmax represents the worst feasible objective function value so far. zxk=∑j=1mmax0,Hjxk2 is the constraint violation function, and wp 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 Section3.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 producens points which are repeatedly generated via a symmetric Latin hypercube design [39]. In Latin hypercube design, each variable of x is stratified into ns 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. Therth 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 theG˜χ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 VR) and (2) minimum distance from previously evaluated points (distance criterion VD) [40]. The weighted score is computed as(39)Wr=wRVRχr+wDVDχr,where wR+wD=1, wR≥0 is the weight for the response surface criterion, and wD≥0 is the weight for the distance criterion. The response surface criterion for every candidate point χr is given by(40)VRχr=G˜χr−GminGmax−Gmin,ifGmax≠Gmin,1,otherwise,where Gmax=maxr=1,…,RG˜χr and smin=minr=1,…,RG˜χr. The distance criterion for every candidate point χr is given by(41)VDχr=Δmax−ΔχrΔmax−Δmin,ifΔmax≠Δmin,1,otherwise,where Δχr=minm=1,…,nsχr−xm, Δmax=maxt=1,…,TΔχr, and Δmin=mint=1,…,TΔχr.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.


## 3.2.1. Framework

The integer (global) optimization problem can be written in the following form:(24)minimizeGx,(25)subject toxjl≤xj≤xju,∀j=1,…,J,where Gx denotes the objective function defined in equation (19). xjl and xju denote the lower and upper bounds on the variable xj, respectively. The optimization model defined in (24) is replaced by the sum of a constant value and a Gaussian random error term as follows:(26)Gx=ρ+Zx,where ρ is the mean value of Gx andZx is a stationary Gaussian random process with mean zero, variance σ2, and nonzero covariance. The covariance is expressed as(27)covZxm,Zxh=σ2ψxm,xh,where covZxm,Zxh is the covariance between two sample points xm and xh. ψxm,xh is the Kriging function, or Gaussian exponential correlation function, equivalently. The Kriging function is defined as follows:(28)ψxm,xh=exp−∑j=1Jλjxjm−xjh2,where λ=λ1,…,λj,…,λJ denotes the vector of unknown correlation parameters and xjm and xjh denote the jth variable in xm and xh, 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 xjm and xjh. The larger of λj, the greater the Euclidean distance between xjm and xjh, hence the lower the correlation between them.Letns denote the number of the known sample points for variable x. The predictive estimate G˜x of the response function Gx at unknown point x is given by(29)G˜x=ρ˜+R˜TxΨ˜−1G−1ρ˜,where G=Gx1,Gx2,…,GxnsT, 1 is the ns×1 unit column, and Ψ˜ is the estimated correlation matrix and is diagonally symmetric where(30)Ψ˜=ψ˜x1,x1ψ˜x1,x2⋯ψ˜x1,xnsψ˜x2,x1ψ˜x2,x2⋯ψ˜x2,xns⋮⋮⋱⋮ψ˜xns,x1ψ˜xns,x2⋯ψ˜xns,xns∈ℛns×ns,where(31)ψ˜xm,xh=exp−∑j=1Jλ˜jxjm−xjh2.R˜x is the correlation vector between the unknown point x and sample points where(32)R˜x=ψ˜x,x1,ψ˜x,x2,…,ψ˜x,xnsT.Therefore, in order to yield the value ofG˜x, the estimators ρ˜ and λ˜1,λ˜2,…,λ˜K, λ˜2, λ˜K need to be obtained, given by(33)ρ˜=1TΨ˜−11−11TΨ˜−1G,(34)λ˜1,λ˜2,…,λ˜J=argmaxλ˜1,λ˜2,…,λ˜J>0−12nsln1nsG−1ρ˜TΨ˜−1G−1ρ˜+lndetΨ˜,where det is the determinant of matrix [35]. According to equation (33), ρ˜ is obtained based on λ˜1,λ˜2,…,λ˜J, λ˜1, λ˜2, λ˜J while λ˜1,λ˜2,…,λ˜J, λ˜1, λ˜2, λ˜J 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λ˜1,λ˜2,…,λ˜K, λ˜1, λ˜2, λ˜K, and obtain the ψ˜xm,xh 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 λ˜1,λ˜2,…,λ˜K, λ˜1, λ˜2, λ˜KStep 3: updateψ˜xm,xh and Ψ˜ with the updated λ˜1,λ˜2,…,λ˜J, λ˜1, λ˜2, λ˜J, 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 λ˜1,λ˜2,…,λ˜J, λ˜1, λ˜2, λ˜JIf the problem has other constraints in addition to the bound constraint (25), it can be formulated as the following form:(35)minimizeGx,subject to(36)Hjx≤0,∀j=1,…,J,(37)xjl≤xj≤xju,∀j=1,…,J.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 xm, [37, 38] which is rewritten as(38)GPxm=Gmax+wpzxm,ifzxm>0,Gxm,otherwise,where GP⋅ denotes the penalty augmented objective function (PAOF). Gmax represents the worst feasible objective function value so far. zxk=∑j=1mmax0,Hjxk2 is the constraint violation function, and wp 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 Section3.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 producens points which are repeatedly generated via a symmetric Latin hypercube design [39]. In Latin hypercube design, each variable of x is stratified into ns 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. Therth 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 theG˜χ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 VR) and (2) minimum distance from previously evaluated points (distance criterion VD) [40]. The weighted score is computed as(39)Wr=wRVRχr+wDVDχr,where wR+wD=1, wR≥0 is the weight for the response surface criterion, and wD≥0 is the weight for the distance criterion. The response surface criterion for every candidate point χr is given by(40)VRχr=G˜χr−GminGmax−Gmin,ifGmax≠Gmin,1,otherwise,where Gmax=maxr=1,…,RG˜χr and smin=minr=1,…,RG˜χr. The distance criterion for every candidate point χr is given by(41)VDχr=Δmax−ΔχrΔmax−Δmin,ifΔmax≠Δmin,1,otherwise,where Δχr=minm=1,…,nsχr−xm, Δmax=maxt=1,…,TΔχr, and Δmin=mint=1,…,TΔχr.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 Figure6(a) demonstrates. The toll plaza for the car park exit has 7-lane approach and 13 tollbooths for cars where NA=7, NR=0, and NL=6. We number the tollbooths from No. 1 to No. 13 consecutively (from the left- to right-hand side in the toll plaza). qiB 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 are10,11 and 1,2,5,6,9, 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.Figure 11
Median and mean of the average delays with different simulation times (a) when 1000 cars arrive, (b) when 2000 cars arrive, and (c) when 3000 cars arrive.
(a)(b)(c)Figure 12
Standard deviation of the average delays with different simulation times (a) when 1000 cars arrive, (b) when 2000 cars arrive, and (c) when 3000 cars arrive.
(a)(b)(c)Figure 13
Median and mean of the total delays with different simulation times (a) when 1000 cars arrive, (b) when 2000 cars arrive, and (c) when 3000 cars arrive.
(a)(b)(c)As Figures11–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.Figure 14
Standard deviation of the total delays with different simulation times (a) when 1000 cars arrive, (b) when 2000 cars arrive, and (c) when 3000 cars arrive.
(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 is1500,2500,3500,4000,4000,3500,3000,2000 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 Figures15–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 growsFigure 15
Recommended tollbooth operation scheme and its performance with weight of operational cost = 0.2. (a) Queue length evolution with time. (b) Best tollbooth operation scheme. (c) Costs in each period.
(a)(b)(c)Figure 16
Recommended tollbooth operation scheme and its performance with weight of operational cost = 0.5. (a) Queue length evolution with time. (b) Best tollbooth operation scheme. (c) Costs in each period.
(a)(b)(c)Figure 17
Recommended tollbooth operation scheme and its performance with weight of operational cost = 0.8. (a) Queue length evolution with time. (b) Best tollbooth operation scheme. (c) Costs in each period.
(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.


## 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 are10,11 and 1,2,5,6,9, 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.Figure 11
Median and mean of the average delays with different simulation times (a) when 1000 cars arrive, (b) when 2000 cars arrive, and (c) when 3000 cars arrive.
(a)(b)(c)Figure 12
Standard deviation of the average delays with different simulation times (a) when 1000 cars arrive, (b) when 2000 cars arrive, and (c) when 3000 cars arrive.
(a)(b)(c)Figure 13
Median and mean of the total delays with different simulation times (a) when 1000 cars arrive, (b) when 2000 cars arrive, and (c) when 3000 cars arrive.
(a)(b)(c)As Figures11–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.Figure 14
Standard deviation of the total delays with different simulation times (a) when 1000 cars arrive, (b) when 2000 cars arrive, and (c) when 3000 cars arrive.
(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 is1500,2500,3500,4000,4000,3500,3000,2000 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 Figures15–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 growsFigure 15
Recommended tollbooth operation scheme and its performance with weight of operational cost = 0.2. (a) Queue length evolution with time. (b) Best tollbooth operation scheme. (c) Costs in each period.
(a)(b)(c)Figure 16
Recommended tollbooth operation scheme and its performance with weight of operational cost = 0.5. (a) Queue length evolution with time. (b) Best tollbooth operation scheme. (c) Costs in each period.
(a)(b)(c)Figure 17
Recommended tollbooth operation scheme and its performance with weight of operational cost = 0.8. (a) Queue length evolution with time. (b) Best tollbooth operation scheme. (c) Costs in each period.
(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.



---

*Source: 6674037-2021-02-23.xml*