Optimizing Multi-Terminal Customized Bus Service With Mixed Fleet

The customized bus (CB) transit is recognized as an effective transportation mode offering more flexible and demand-responsive service than traditional bus transit with fixed route and schedule, especially during the peak hours. The novelty of this study is the development of a mixed integer non-linear model for optimizing multi-terminal CB service in an urban setting. According to the estimated spatiotemporal passenger demand, the objective total cost, consisting of supplier’s and users’ costs, is minimized subject to capacity and time window constraints. A mixed bus fleet with various bus sizes is employed to accommodate passenger demand, which increases vehicle utilization and reduces supplier’s cost. The inconvenience of passengers caused by early arrival at the destination is treated as penalty and considered in users’ cost. The study optimization problem is combinatorial with many decision variables including trip assignment, bus routing and associated timetables, and fleet size. A hybrid genetic algorithm (HGA) which integrates the features of genetic algorithm (GA) and simulated annealing (SA) is developed to effectively search for the optimal solution. A real-world CB network is employed to demonstrate the applicability of the developed model and explore the relation between the model parameters and optimized results. It was found that the total cost can be reduced by 16.5% after employing multiple terminals and a mixed bus fleet.


I. INTRODUCTION
Rapid urbanization has led transportation demand drastically increases, which deteriorates the service quality and efficiency of transportation systems, especially during the peak hours. Commuters are tired of lengthy travel time in crowded vehicles. Transit agencies thus promote demand responsive transit (DRT) as an alternative mode to elevate service quality by offering greater mobility and accessibility during the peak hours. With the rapid development of mobility as a service (MaaS) and personal mobile devices, an emerging DRT service, called customized bus (CB), has been initiated [1], [2].
The CB concept originates from car-sharing, subscription bus, and dial-a-ride services, which could be generalized as a dial-a-ride problem (DARP). Unlike traditional DARP a door-to-door service primarily offered to elderly and disabled people [3] with smaller vehicles, the CB service is operated on a stop-to-stop basis to serve greater volume of The associate editor coordinating the review of this manuscript and approving it for publication was Yanbo Chen . passengers with larger vehicles considering spatiotemporal demand [4].
The process of acquiring the CB service begins by collecting travel requests from passengers, in which the information related to origination and destination (OD), desired arrival time, and subscription duration must be given. The operators then process the requests to develop a service plan (e.g., stop location, routing, and scheduling) which would yield the objective (e.g., minimize the total cost in this study) subject to a set of practical constraints. The generated bus routing and scheduling information is then fed back to users for reservation. The service plan will be updated on a short-term basis (e.g., one to few weeks) subject to the demand change [5]. It is worth noting that CB does not serve walk-in passengers without reservation.
To optimize the study CB problem, a sound model with an efficient solution algorithm is developed. Several practical issues concerned by the CB suppliers but missed in previous studies are addressed, which include: (1) The desired arrival time at destination requested by passengers is ensured, while the inconvenience caused by early arrival at the destination is treated as penalty in user's cost. (2) Considering spatiotemporal passenger demand, a mixed bus fleet with various sizes is employed to increase the vehicle utilization and reduce supplier cost. (3) To reduce deadheading travel distance and associated cost, the multi-terminal operation is introduced.
The remainder of this paper is organized as follows. In Section II, an overview of related research in the existing literature is summarized. Section III discusses the design process for CB service and model assumptions of the problem. The derivation of the objective function and the related constraints are also presented. The developed solution algorithm is discussed in Section IV, and its effectiveness is assessed in Section V. Section VI presents a case study, in which the optimized solutions and results of sensitivity analysis are discussed. Finally, the research findings are summarized and future research is discussed in Section VII.

II. LITERATURE REVIEW
Most previous studies focused on the CB service considered single terminal operation and the sizes of buses are identical. Few of them jointly optimized trip assignment, bus route and the associated timetable, mixed bus fleet for a multi-terminal CB network. The findings in the literature review are discussed next.

A. TRANSIT NETWORK DESIGN (TND)
Transit network design (TND) problem can be formulated as a mathematical model, so that bus routing and associated attributes including frequencies, timetables and stop spacing can be optimized. Chien and Spasovic [6] optimized stop spacing, bus headway, and fare considering demand elasticity for a grid transit network, which maximize supplier's profit. Baaj and Mahmassani [7] minimized total cost, which was yielded by optimized bus routes and frequencies.
To yield greater benefit (e.g., reduced cost and increased profit), Chien et al. [8] optimized the integration of service patterns to elevate service quality and to stimulate demand. Chien et al. [9] assessed conventional and subscription bus systems considering stochastic demand and non-additive value of time. Chien and Qin [10] optimized stop locations which minimized total cost considering a realistic road network and spatial demand. Chien [11] optimized bus size, route and headway that minimized total cost for a feeder service. DiJoseph and Chien [12] optimized stop locations, headway and fare that maximized total profit for a feeder bus service.
Ulusoy et al. [13], Ulusoy and Chien [14] optimized bus service patterns and frequencies which minimized total cost considering demand elasticity. Later, Qu et al. [15] extended that work by considering time-dependent demand. Chen et al. [16] jointly optimized the bus routes, frequencies, locations of wireless power devices and battery capacity for an electric feeder bus system, which minimized the total cost.

B. DEMAND-RESPONSIVE TRANSIT (DRT)
Demand-responsive transit (DRT) is a form of shared private transport (or quasi-public) for groups travelling where buses alter their routes based on particular demand rather than using fixed routes or timetables. Navidi et al. [17] indicated that DRT may improve the mobility and quality of service by reducing the users perceived travel time. Wang et al. [18] found that DRT meets the demand of niche markets (e.g., disabled, elderly and live in areas with low population density). Cordeau [19] optimized vehicle assignment and routing, which minimized the operation cost. Later, Parragh [20] extended that work by considering heterogeneous vehicles and demand.
Molenbruch et al. [21] analyzed operational effects of service level variations for a dial-a-ride problem. It was found that the tradeoffs between service quality and costs should encourage service providers to make informed decisions regarding potential changes in the service level they offer. Shen et al. [22] optimized bus routing and scheduling that minimized the total cost of a demand-responsive feeder transit. Beaudry et al. [23] optimized bus routing and scheduling for a paratransit to serve patients with medical appointments at hospitals, which minimized the cost of patients.

C. CUSTOMIZED BUS (CB) TRANSIT
Liu and Ceder [24] presented a systematic and comprehensive analysis on a CB system. Liu et al. [1] assessed the overall performance metrics (e.g., travel costs, travel time, and fuel consumption) of CB, private car (PC), and traditional public transport (PT). It was found that CB is a more attractive commuting alternative than PT and PC.
Wang et al. [2] found that CB is very attractive especially to commuters traveling in the peak hours because of its features in reducing travel time and coverage of spatiotemporal demand. It was also found that CB could serve as a supplement mode in the areas with low PT coverage.
With a density-based clustering algorithm, Qiu et al. [25] estimated the origin and destination demand for a CB network. Ma et al. [26] optimized the CB route considering operation cost and social benefit. Li et al. [27] optimized bus routing which minimized the operating cost of CB service providers. Cao and Wang [28] optimized trip assignment that minimized total cost. Ma et al. [29] optimized locations of stops, routes and timetables that minimized total cost including operation, environmental and congestion cost.
Guo et al. [30] optimized stop locations, bus routes, and trip assignment that minimized total cost, and then extended that work by considering service time window [31]. Huang et al. [34] designed a two-phase model to optimize trip assignment, routes, timetables, which maximized total profit. Tong et al. [5] developed a commodity flow model to optimize trip assignment and bus routing subject to spatiotemporal window constraints, which maximized ridership. Lyu et al. [32] developed a framework called CB-Planner to jointly optimize stop locations, bus routing and timetables, and probabilistic mode choice, which maximized total profit. Han et al. [33] presented an analytical model that minimized total cost and demonstrated that a mixed-size of bus fleet is cost-effective than a single-sized bus fleet. According to the findings discussed above, Table 1 summarizes the features of key studies focusing on the CB service.

D. SOLUTION ALGORITHM
Since the study CB optimization problem is combinatorial that is non-deterministic polynomial-time hard (NP-hard) [34]- [36], an efficient algorithm is desired to search for the optimal solution. Genetic algorithm (GA) is a stochastic search algorithm that finds the optimal solution space through drawing on the natural laws (i.e. survival of the fittest) of the biological world. GA does not require to calculate the derivatives of the objective function, and its performance in optimizing complex mathematical models (e.g. vehicle routing problems [36]- [39]) is promising [8].
Many studies have employed GA to search for optimal solution. Uchimura et al. [40] solved a routing problem for an advanced public transit system. Shrivastava et al. [41] optimized transit routes and frequencies. Bagloee and Ceder [42] optimized transit routes and schedules. Zhang et al. [43] developed a GA to solve a multi-trip DARP.
However, traditional GA suffers some issues while searching for the solution (e.g., premature solution or local optimal solution [44]). To cope with this problem, some studies integrated GA with features of other algorithms to enhance the performance. Masmoudi et al. [37] proposed a hybrid GA with a local search strategy to solve a DARP. Chen et al. [16] developed a nested genetic algorithm (NGA) for an electric feeder bus system. Rekiek et al. [45] optimized bus routes for serving disabled passengers with a group genetic algorithm (GGA). Zhao and Zeng [44] developed an algorithm which integrated GA and simulated annealing (SA), called GASA, to optimize vehicle routing and headway for a large-scale transit network. Song et al. [46] improved GASA and applied that to minimize vehicle emission on arterial roads. The results suggested that GASA has better efficiency and robustness than GA.

III. METHODOLOGY
In this section, we formulate a model for optimizing the CB service. The passengers are classified into groups based on the OD information and desired arrival time. To optimize the CB service, the decision variables include trip assignment, routes, timetables, and fleet size. The objective is to minimize total cost, consisting of the supplier and user costs, subject to capacity and time window constraints.

A. PASSENGERS PARTITION
Assuming that the stop and terminal locations, passenger OD demand, and available parking spaces at each terminal are given. For each OD pair, we partition passengers into groups based on desired arrival time at the destination stop. Passengers with the same OD are grouped in a set, denoted as C h , where 1 ≤ h ≤ n, and n is the number of OD pairs. For each C h , passengers are arranged as an ascending order according to the desired arrival time. They are partitioned into α groups based on pre-specified time interval (e.g. 15 minutes). After looping for n sets, the total number of groups m is equal to n multiplied by α.
Based on the earliest desired arrival time in each group, a service time window (e.g., 15 minutes) on the drop-off stop is constructed.

B. ASSUMPTIONS
As shown in Figure 1, the passenger travel time includes wait time and in-vehicle time. The wait time is minor and negligible as the real-time bus arrival information is known via advanced traveler information systems. Therefore, the travel time concerned here is in-vehicle time, which is determined by the drop-off time minus the pick-up time. The penalty incurred by early arrival passengers at destination is determined by the elapse time from the drop-off time to the desired arrival time, multiplied by the value of time.
To formulate the proposed model, some assumptions are described as follows: • The travel distance and associated travel time between any pair of nodes are known, and the dwell time at stops is fixed.
• A bus initially parks at a terminal and may return to a different one after the journey. Parking fees at terminals are considered in the supplier cost.
• The spatiotemporal passenger demand distribution is fixed within the designated service period.
The model is defined on a complete graph G = (N , A).
is a set of nodes including stops and terminals, where N + is a set of pick-up stops; N − is a set of drop-off stops; and W is a set of terminals. A is a set of links that connect every two vertices in N . The notations used to formulate the proposed model, including variables and parameters, are defined and listed in Table 2.

C. OBJECTIVE FUNCTION
The objective total cost formulated as Eq. (1) is a mixed integer non-linear programming model (MINLP), consisting of supplier's cost denoted as Z 1 and user's cost denoted as Z 2 . Thus, where x wjk and x ijk are binary variables, which are equal to 1 if bus k travels from nodes w and i to j. Z 2 consists of travel time cost and early arrival penalty incurred by passengers. The in-vehicle time from stop i to stop j by bus k is T jk less T ik . Thus, travel time cost is the product of in-vehicle time and value of passengers' time λ v . The early arrival time for passenger p at drop-off stop j by bus k is T pj less T jk . Thus, early arrival penalty is the product of (T pj -T jk ) and value of early arrival time denoted as λ e . Thus, (3) where y ik and y jk are binary variables, which are equal to 1 if stops i and j are served by bus k.
To minimize the total cost, a set of practical constraints concerning trip assignment, bus capacity, routing and timetabling are considered. Eq. (4) restricts that each stop could be served by more than one bus, while Eqs. (5-6) state that if a bus k is used, the bus starts and ends at one terminal at most.
Eqs. (7)(8) restrict that each bus route starts and ends at terminals. Eq. (9) ensures the continuity of a bus route.
The subtour elimination constraint is expressed by Eq. (10), which ensures that every bus route must begin and end at a terminal. u ik and u jk are arbitrary non-negative integers, which are employed to eliminate the subtour. Thus, Eq. (11) indicates that there is no passenger in the bus when it departs from and arrives at a terminal. Eq. (12) calculates the load of bus k departing from stop j, which is denoted as L jk . As formulated as Eq. (13), H jk is the number of passengers boardings/alightings served by bus k at stop j, which is positive if stop j is a pick-up stop; otherwise, it is negative if stop j is a drop-off stop. Eq. (14) ensures that the number of passengers in a bus is less than or equal to the capacity.
156460 VOLUME 8, 2020 The arrival time of bus k at node j denoted as T jk can be determined by Eq. (15). The service time window constraint is formulated as Eq. (16), while Eq. (17) ensures that the pair and precedence constraints sustain while searching for the optimal solution.

IV. SOLUTION ALGORITHM
The optimization of the proposed model discussed in Section III is a combinatorial problem. The solution is difficult to be optimized using the exact algorithm, especial for large-scale networks [31]. In order to efficiently search for the solution, a hybrid genetic algorithm (HGA) that integrated the features of GA and SA is proposed.

A. FITNESS FUNCTION
For maximization problems, the fitness function is same as the objective function. However, for minimization problems, the fitness function is the reciprocal of the objective function.
To ensure that the solutions always satisfy bus capacity and time window constraints, a large penalty Z p is introduced. Thus, the fitness function can be formulated as follows: where F and Z are fitness value and objective value, respectively. The minimum objective value (e.g., total cost) represented by the maximum fitness value would suggest the best solution found by HGA.

B. CHROMOSOME ENCODING
Integer coding is used for encoding bus routes. The chromosome is encoded as an integer string, whose length depends on fleet size and number of passenger groups. Each string consists of several substrings, and each substring is a sequence of genes representing stops of a bus route that begins and ends at a terminal. The chromosome is randomly generated and an example is shown in Figure 2 in which 5 passenger groups and 2 buses are presented, while 1, 2, 3, 4 and 5 represent the pick-up stops, and 6, 7, 8, 9 and 10 represent the drop-off stops, respectively. Note that the gene with 0 represents a bus terminal.

C. GA OPERATION
There are three GA operators discussed next, including selection, crossover, and mutation. Selection -The tournament selection [47] is applied here, which randomly select a number of chromosomes from the population and then choose the one yielding the best fitness value. This process is repeated until the number of chromosomes required (e.g., population size) is reached. Crossover -A random two-point crossover with a minor alteration to respect the pair and precedence constraints is applied, and the process is shown in Figure 3. In Phase 1, a gene in the chromosome is randomly selected (e.g., the fifth gene marked with an arrow). Then, the terminals (represented by 0s) associated with a bus route, marked in green, containing the gene are chosen as crossover points. In Phase 2, the genes between two crossover points in parent chromosomes 1 and 2 are swapped. And the remaining genes are inherited in child chromosomes. In Phase 3, the replicated genes in child chromosome 1 (i.e., 4 and 9 marked in orange) are deleted, and the missing genes (i.e., 4 and 9 marked in blue) are inserted in child chromosome 2.
Mutation -A random two-point mutation is applied and show in Figure 4. The genes of two randomly selected pick-up stops assigned to different bus routes in parent chromosome (e.g., 4 and 3 marked in orange) are exchanged, as well as the genes of corresponding drop-off stops (e.g., 9 and 8 marked VOLUME 8, 2020 in green). The rest of genes in the child chromosome are inherited by the parent chromosome.

D. SA OPERATION
The purpose of SA operation in HGA is to prevent GA trapped in a local optimum, which permits the GA process accepting a slightly worse solution with a probability based on the Metropolis criterion [46] as formulated in Eq. (23).
where Pro is the probability of accepting a new solution, and F new and F old represent the fitness values of the new and the old solutions, respectively. The term γ T is the temperature of the current iteration where γ represents changing rate of temperature, and T is the temperature of the previous iteration. If F new ≥ F old , the new solution is accepted; otherwise, it can be accepted subject to a certain probability. The higher the current temperature, the greater the probability is. However, the temperature decreases at a rate γ , and the probability of accepting a worse solution reduces. The integration of the SA procedure into GA is discussed next. First, for each solution (chromosome), two random genes of the pick-up stops are chosen, and a new solution is discovered by exchanging the two genes, as well as the genes of the corresponding drop-off stops. Second, the fitness value of the new solution is compared to the old one. The selection of a solution is determined by Eq. (23).
A step-by-step procedure to search for the optimal solution using the proposed HGA is shown in Figure 5 and discussed below.
Step 1 -Parameter setting. Specify the parameters for GA and SA, including initial temperature, the rate of temperature change, population size, maximum iterations, and crossover and mutation probability.
Step 2 -Initialization. Create the initial population as the first generation.
Step 3 -Evaluation. Determine the origin and destination terminals and bus type for each route and then calculate the value of fitness function with Eq. (22).
Step 4 -SA Operation.Implement the SA process to generate a new population.
Step 5 -GA Operation. Implement the GA process to reproduce new solutions.
Step 6 -Termination. Check if the specified termination criteria (i.e., maximum number of iterations) is attained. If not, go to Step 4; otherwise, terminate the algorithm and report the best solution.

V. ALGORITHM PERFORMANCE ANALYSIS
In this section, we assess the performance of HGA against LINGO and GA using a set of CB networks with different number of stops. The networks are classified into three scales: small (e.g., 8 and 10 nodes), medium (e.g., 16, 24 and 32 nodes), and large (e.g., 40 and 48 nodes). The minimized total cost with the optimized decision variables were found by LINGO (Global optimal solver, 18.0 version) as well as by GA and HGA (MATLAB R2018b on a laptop computer, Intel Core i5, 8G, 1.8GHz).

A. PARAMETER SETTING
The HGA parameters consist of population size, maximum number of iterations, crossover and mutation probability, initial annealing temperature and the rate of temperature change, which are listed in Table 3. The parameter values were determined by the effectiveness to search the optimal solution and stability of the yielded results in the simulation analysis. Table 3, population size and maximum number of iterations increase as the network size increases. In general, GA needs more iteration to converge to an acceptable result than HGA, and HGA requires less population size than GA to yield the acceptable result. Crossover and mutation probability are 0.9 and 0.1 respectively with both GA and HGA. The initial annealing temperature and the rate of temperature change with HGA are 1000 and 0.95, respectively.

B. RESULTS ANALYSIS
Ten simulation runs were executed for each experiment so that the randomness of GA and HGA results may be assessed. The result statistics are summarized in Table 4, while the average computation times with LINGO, GA, and HGA are illustrated in Figure 6. It is worth noting that the upper bound of computation time is set as 10,800 seconds. Figure 6 suggests that LINGO is able to optimize small-scale CB networks (e.g., 8 and 10 nodes). As the number of nodes increases, the computation time exponentially escalated and failed to find the solution within the time limit. GA and HGA are able to find solutions within the time limit, and HGA outperformed GA in terms of shorter computation time as well as less minimized total cost and randomness of the results. As the number of nodes increases, the result randomness increases (e.g., see SD) with GA and HGA. The differences between the minimized costs denoted as Z G and Z H with GA and HGA, respectively, against that with LINGO denoted as Z L , are represented by '' GA '' and '' HGA '', formulated as Eqs. (24) and (25). Thus,

VI. CASE STUDY
To assess the effectiveness and applicability of the developed model and algorithm, a real-world CB network in Xi'an City is applied. The study network consists of 19 stops and 5 terminals as shown in Figure 7. The locations of pick-up and drop-off stops cover large communities and three business areas. Five major stations in Xi'an are chosen as the candidate terminals. Table 5 lists a sample of 326 requests (from home to work), which include the pick-up and drop-off stops and the desired arrival time in the AM peak from 8:00 AM to 9:30 AM. The value of passengers' time is 6 $/pass-hr, and the value of early arrival time is 9 $/pass-hr. Bus types and associated fixed and variable costs are listed in Table 6. The parameters of GA and HGA have been determined earlier and listed in Table 3.

A. RESULTS ANALYSIS
Based on the collected data of passenger reservation, we found that the desired arrival time at destination varies within the study period. Therefore, the time interval applied for aggregating demand will affect the supplier and user costs.     Table 7, the optimal results with varying time interval suggest that the demand estimated on a 10-min basis would yield the least total cost. The optimal CB routes and schedules are summarized in the Appendix.

As shown in
It was indicated that as the time interval decreases, the supplier cost increases, but the user cost reduces. The increase of the supplier cost results from the increase of fleet size. The decrease of users' cost mainly results from reduced travel time because the number of stops per route decreases. The trade-off between the supplier and user costs can be observed in Table 7. The total cost found by GA and HGA over iteration is illustrated in Figure 8. The results suggest that HGA is more efficient because the total cost found by HGA converges significantly quicker than that with GA. To demonstrate the effect of multi-terminal, mixed bus fleet, and early arrival penalty to the minimized total cost, ten simulation runs for each of five scenarios are summarized in Table 8.
Scenario 1 -Single terminal and single bus size without early arrival penalty.
Scenario 2 -Multi-terminal and single bus size without early arrival penalty. Scenario 3 -Single terminal and multi-bus-size without early arrival penalty.
Scenario 4 -Multi-terminal and multi-bus-size without early arrival penalty.
Scenario 5 -Multi-terminal and multi-bus-size with early arrival penalty.
In Scenario 1, the minimized total cost is 4,210 $, consisting of supplier cost of 2,335 $ and user cost of 1,875 $. The average travel time and early arrival time per passenger are 45.8 min and 7.8 min, while the average travel and deadheading distances per bus are 28.9 and 8.1 km, respectively. The optimal bus capacity is 15 seat/veh, and the fleet size is 28. The suggested terminal location is node 20.
Under Scenario 2 with multiple terminals, the minimized total cost can be reduced by 5.9% (from 4,210 $ to 3,961 $). The supplier cost significantly reduced by 13% (from 2,335 $ to 2,030 $), albeit the user cost slightly increased. The supplier cost reduction mainly resulted from reduced variable cost since the average bus travel distance decreased from 28.9 km to 24.8 km (14%), and the average deadheading distance decreased from 8.1 km to 3.6 km (56%).
Under Scenario 3 with mixed bus fleet, the minimized total cost can be reduced by 11.3% (from 4,210 $ to 3,734 $). The supplier cost significantly reduced by 20.3% (from 2,335 $ to 1,861 $) while the user cost remained the same. The supplier cost reduction was mainly from the reduced fleet size (from 28 buses to 22 buses).
Comparing the results between Scenario 4 and Scenario 2 (and 3), we found that the minimized total cost can be further reduced by 11.2% (and 5.8%). The supplier cost significantly reduced by 18.8% (and 11.4%), while the user cost slightly reduced by 3.3% (and 0.3%). The results suggest that the integration of mixed bus fleet and multi-terminal can further decrease the total cost, mainly from supplier's side.
Finally, comparing the results between Scenario 5 and Scenario 4, we found that with early arrival penalty, the total cost can be reduced by 1.1% (from 3,516 $ to 3,478 $). The user cost reduced by 2.6% (from 1,868 $ to 1,820 $), while the supplier cost slightly increased. The user cost reduced mainly from the reduction in early arrival penalty (from 382 $ to 317 $). The supplier cost increased mainly caused by the increase of fleet size (from 22 buses to 24 buses). It indicated that the joint effect of multi-terminal, mixed bus fleet, and early arrival penalty shall be jointly considered to yield the least total cost. To summarize the results, the minimized total costs under five scenarios are shown in Figure 9.

B. SENSITIVITY ANALYSIS
A sensitivity analysis based on the condition with Scenario 5 is performed to explore the relation between model parameters and optimized results. Figure 10 suggests that as the value of early arrival time λ e increases, the supplier cost would increase because more buses are needed to reduce passenger travel time and early arrival time.   Figure 11 shows that λ e significantly affects the fleet composition. As λ e increases from 3 to 39 $/pass-hr, the number of minibuses (10-seat) increases from 9 buses to 18 buses, while the number of larger buses (15-seat and 20-seat) VOLUME 8, 2020 decreases from 14 buses to 10 buses. The increase of fleet size as discussed earlier would reduce passenger travel time and early arrival time, while using more smaller vehicles may reduce the supplier cost and increase vehicle utilization. In addition, we alter the variable cost via a cost multiplier shown in Figure 12 and investigate the variation in costs, travel time and fleet size. The results indicate that as the cost multiplier increases, the total cost increases. fleet size is expected to decrease, and thus larger buses are preferable, albeit the user cost would increase because of escalated travel time and early arrival time. Figure 13 shows that the cost multiplier would significantly influence on the fleet composition. For example, as the multiplier increases from 0.2 to 2, the number of small buses (10-seat and 15-seat) decreases from 27 buses to 14 buses, while the number of larger vehicles (20-seat) increases from 2 buses to 8 buses. As discussed earlier, the decrease of fleet size using larger vehicles would increase passenger travel time and early arrival time. Figure 14 suggests that as the value of passengers' time λ v increases, the supplier cost would increase because more buses are needed to reduce travel time and early arrival time.    Figure 15 shows that the increase of λ v significantly affect fleet composition. As λ v increases from 2 to 30 $/pass-hr, the number of minibuses (10-seat) increases from 10 buses to 21 buses, while the number of larger vehicles (15-seat and 20-seat) decreases from 13 buses to 7 buses. The increase of fleet size as discussed earlier would reduce passenger travel time and early arrival time, while using smaller vehicles would reduce the supplier cost and increase vehicle utilization.

VII. CONCLUSION
Considering the spatiotemporal passenger demand, this study developed a model to optimize the CB service, which minimized the total cost subject to capacity and time window constraints. The decision variables included bus routing and the associated timetables, and passenger trip assignment, while multi-terminal, mixed bus fleet and early arrival penalty were considered. A hybrid genetic algorithm (HGA) was developed, which can effectively search for the acceptable solution with minor randomness.
The optimal solution was found through optimizing a real-world CB service in Xi'an China. A sensitivity analysis was conducted, and the impacts of model parameters to the optimized results were explored.
The results suggest that employing multiple terminals and mixed bus fleet individually, the minimized total cost can be reduced by 5.9% and 11.3%, respectively. However, employing multiple terminals with mixed bus fleet may reduce the cost by 16.5%, and the major cost reduction is on the supplier's side. Moreover, as early arrival penalty is also considered, the cost can be further reduced by 17.4%.
The results of sensitivity analysis suggest that as the value of time (e.g. early arrival time and travel time) increases, an increase in supplier cost resulting from the increased fleet size shall be expected, and the early arrival time and average travel time may significantly decrease. As fleet size increases, smaller buses are preferable to reduce cost and improve vehicle utilization. On the other hand, as the bus operating cost increases, fleet size decreases, and larger vehicle is preferable. This situation would increase the user cost because of increased travel time, resulting from increased route length and number of stops per route.
As an immediate extension of this study, the developed model may be enhanced by considering more realistic conditions, such as stochastic travel time between stops, the impact of demand elasticity to mode/route choice, and dynamic dispatching strategies. This extension would permit the model to optimize CB operation and service planning subject to more practical constraints. APPENDIX See Table 9.
QIAN SUN received the B.S. degree in transportation engineering from the School of Automobile, Chang'an University, Xi'an, China, in 2013, where she is currently pursuing the Ph.D. degree in transportation planning and management with the School of Transportation Engineering. Her research interests include transit route planning, demand-responsive transit network optimization, intelligent transportation modeling, and transportation system analysis.