A Model for the Stop Planning and Timetables of Customized Buses

Customized buses (CBs) are a new mode of public transportation and an important part of diversified public transportation, providing advanced, attractive and user-led service. The operational activity of a CB is planned by aggregating space–time demand and similar passenger travel demands. Based on an analysis of domestic and international research and the current development of CBs in China and considering passenger travel data, this paper studies the problems associated with the operation of CBs, such as stop selection, line planning and timetables, and establishes a model for the stop planning and timetables of CBs. The improved immune genetic algorithm (IIGA) is used to solve the model with regard to the following: 1) multiple population design and transport operator design, 2) memory library design, 3) mutation probability design and crossover probability design, and 4) the fitness calculation of the gene segment. Finally, a real-world example in Beijing is calculated, and the model and solution results are verified and analyzed. The results illustrate that the IIGA solves the model and is superior to the basic genetic algorithm in terms of the number of passengers, travel time, average passenger travel time, average passenger arrival time ahead of schedule and total line revenue. This study covers the key issues involving operational systems of CBs, combines theoretical research and empirical analysis, and provides a theoretical foundation for the planning and operation of CBs.


Introduction
The conflict between increasing traffic demand and a relatively lagged traffic supply is becoming increasingly prominent with the rapid development of economy and society and the accelerating process of urbanization. Traffic congestion, traffic environment pollution, traffic accidents, energy consumption and societal fairness problems are widespread in all large and medium-sized cities in our country. These difficult problems perplex the majority of urban travelers and government administration. To satisfy the increasing and changing traffic demand for citizens' travel and minimize the network traffic load, Beijing proposed building a a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 "public transport city" in 2009 and stated that managers should promote the construction of a green travel system that coordinates rail transit as the backbone, ground buses as the main body, walking and cycling. In recent years, CBs have come online in many large and mediumsized cities in China due to the rapid development of information and communication technology, particularly the widespread use of smartphone apps. The first CB was implemented in Beijing in October 2013; since December 2015, 287 CB lines have been implemented in Beijing. The average demand for a CB has increased from 200 in 2013 to more than 100,000 at present. With such a large-scale passenger demand for CBs, many large and medium-sized cities in China and other developing countries are now faced with determining how to scientifically plan, design, and operate CBs, improve the level of CB service, reduce the costs of operating CBs, facilitate passenger travel, and alleviate urban road congestion and smog.
The car-sharing problem, which is regarded as the earliest concept of CBs, was first conceived in Zurich in 1948 [1]. Kirby and Bhatt [2] analyzed ten cases of a subscription bus, which represents the implementation of the car-sharing idea in public transport. The authors provided guidelines on the planning, organization and operation of subscription bus services. Kirby and Bhatt [3] then noted seven main features of subscription bus services, such as a relatively large concentration of at least 50 fairly long trips with proximate origins or destinations, specialized organization operation management, constant adjustment of lines and schedules to meet demand, and guaranteed seats for personalized service. Based on the work of Kirby and Bhatt, a cost comparison between several types of subscription bus services was made by Bautz [4]. Bautz indicated that the least cost services were provided when the private carriers operated a subscription bus with the support of the government. McKnight and Paaswell [5] designed a subscription bus service in Chicago and indicated that it could help reduce the peak demand and deficit on certain commuter railroad lines. The authors discussed key elements of the subscription bus service that was successfully implemented in the public transport sector. Based on the travel demand of diversification and humanization, CBs came into being as a new travel mode. Shaheen et al. reviewed the development of CBs from the concept of "car sharing" and summarized the advantages of CBs [6], but CB networks and timetable modeling was not examined. Eiro et al. [7], Martinez et al. [8], and Lopes et al. [9] proposed methods based on clustering analysis and a multi-agent model to solve the network planning and timetable problems of a CB in Lisbon, Portugal. However, the dynamic and real-time of passenger demand were not considered in these articles. De Lorimier et al. [10] used the multi-hierarchical regression analysis method to determine the decisive effect of a CB system in Montreal on the effectiveness of vehicle use, providing a reference for the establishment or expansion of CB networks. Tao Liu [11] systematically reviewed the background and implementation of CBs, analyzed various stages of the CB operation planning process and summarized the shortcomings of CBs. The abovementioned studies regard CBs as a new mode of transportation, operation pattern and planning; comparative analyses between CBs and conventional bus are provided, but the stop planning and timetables of CBs are not involved. Based on passenger demand, operator characteristics, and social benefits, this article further examines the model and solution method for CB stop planning and timetables.
Bus stops design affects the alighting and boarding time of passengers and the dwell times of buses; thus, it is necessary to examine the stop design for CBs. In this section, we briefly review the relevant literature on bus stop design problems. Rodrigo Fernández [12] presented a microscopic model for the operation of public transport stops, mainly bus stops and light rail transit stations, from the perspective of traffic analysis. The authors showed that stop design should incorporate the possibility of allowing exact arrival and departure patterns for both vehicles and passengers. The proposed model provides more detailed information about stop operations. Gu W, Cassidy M J, and Li Y [13] considered curbside bus stops of the type that serve multiple bus routes and that are isolated from the effects of traffic signals. A Markov chain embedded in the bus queuing process was used to develop steady-state queuing models for two special cases of this stop type. The methods for the two special cases were used to derive a closed-form or parsimonious approximation model for general cases. To estimate the service time at a curbside bus stop, a compound Poisson service time estimation model (CPSTM) was proposed by Bian B, Zhu N, Ling S, et al. [14]. The model considered the interactions among arriving buses and the number of boarding and alighting passengers. Four different scenarios that occur in curbside bus stops were established, and the corresponding probability models were formulated based on the Poisson process. The abovementioned studies demonstrate that bus stop planning should consider the capacity, congestion, queuing, and serve time of bus stops. These modes and methods have provided great help in solving the stop planning and timetables of CBs.
Although there is little research on CB network design and timetables, it is worthwhile to provide a rather extensive review of conventional bus design models. Lampkin [15] and De Hsu [16] proposed methods for determining public transport networks and lines by classifying public transport lines and stops based on demand. Ceder, A. and N. H. M. Wilson [17] summarized research on public transport network design and proposed and solved a new model incorporating the benefits of passengers and operators. The proposed approach was easier to implement and required a smaller data set. Patanik [18] divided the route network design problem into two stages. First, a set of candidate routes competing for the optimal solution is generated. Second, the optimal set is selected using a genetic algorithm. Bielli et al. [19] described a method for transit line design using a genetic algorithm. In the proposed method, each iteration is based on the computational performance index of the distribution result, including the distribution demand in the existing network. To calculate the fitness function value, these indexes serve as input for each network in the multi-criteria analysis. The transit route network design problem was formulated as an optimization problem of minimizing the sum of the operating cost and the generalized travel cost by Agrawai [20]. In this study, two parallel genetic algorithm models are proposed. Through a case study, these models are tested with respect to computation time, speedup, and efficiency. Xiaolei M A, WU, YaoJan, et al. [21] proposed a data mining method capable of identifying travel patterns for individual transit riders using a large smart card dataset. The travel patterns could be beneficial to understanding the variability of urban travel behavior and facilitating network design. For more information regarding other models of network design and timetables, the reader is referred to Reference [22], Reference [23] and Reference [24].
The bus network problem always consists of a tremendous number of stops and stations, which makes the scale of the problem too large and intractable to obtain an optimal solution. GA is a flexible meta-heuristic algorithm that has shown outstanding performance in solving large-scale public transit problems. Mazloumi [25] solved the minimum cost solution using genetic algorithm and ant colony algorithm. The experiments indicated that optimal solution were obtained. A large city public transport network design problem consisting of a complex network topology, multiple modes and many-to-many travel demand was solved with a parallel GA by Ernesto Cipriani et al. [26]. Arbex [27] proposed an alternating objective genetic algorithm to efficiently solve the multi-objective public transport network design and frequency setting problem. The algorithm could overcome the problems of a large search space and multiple constraints. Ngamchai and Lovell [28] proposed a new model demonstrating how genetic algorithms can be manipulated to help optimize bus transit routing design, incorporating unique service frequency settings for each route. The model was applied on a benchmark network to test its efficiency, and performance results were presented. The results showed that the proposed model is more efficient than the binary-coded genetic algorithm benchmark. Based on solution methodologies for transit network planning and scheduling in the abovementioned studies, this paper further examines solution methodologies for the stop planning and timetables of CBs. In addition, the abovementioned studies have shown that genetic algorithms are highly efficient, accurate and scalable in models of CB stop planning and timetables. Therefore, the model proposed in this paper uses an improved immune genetic algorithm (IIGA).
The reviewed literature reveals that the existing network planning and scheduling approaches focus on conventional buses. Although a CB is part of an urban public transit system, it is different from a conventional bus. Thus, the abovementioned models and methods are not completely suitable for CBs. In modern large cities, the commuting distance is always long. The commuter hopes to enjoy the features of a high-quality bus travel service, such as on-time departure, no transfers, and on-time arrival. However, the complicated road conditions that occur during peak hours makes the traditional bus system consistently unreliable and unpunctual. CBs provide personalized, advanced and flexible demand-responsive transit service to commuters or other specific clientele. CBs constitute a demand-based transit system that aggregates the travel demand of individual passengers. CBs have fixed stops, lines, vehicles, timetables and prices; furthermore, CBs have the characteristics of serving one person, one direct stop and lane flexibility. CB stops include boarding stops in the boarding zone and alighting stops in the alighting zone. There are no intermediate stops. In addition, timetables of every stop in CBs are open to passengers. Based on large-scale travel demand, this paper proposes a model for the stop planning and timetables of CB. The model considers the comprehensive benefit of passengers, operators, and society. An IIGA is then designed to solve the model, and a real-world CB example based in Beijing is provided. The main improvements of the IIGA are multiple population design and transport operator design, memory library design, mutation probability design and crossover probability design, and the fitness calculation of the gene segment. Due to the characteristics of CBs, vehicle can use public transport lanes at peak times. The line in the middle area is not fixed, except for boarding stops and alighting stops, which must search for the shortest path based on traffic flow conditions. Therefore, this paper studies only CB stops in the boarding zone and alighting zone and does not use the driving line in the middle of the zone. This research provides an approach for public transport companies to establish CBs in metropolitan areas.
The remainder of this paper is organized as follows. The CB stop model and the formulation for stops and timetables of CBs are described in Section 2. Section 3 builds the model for the stop planning and timetables of CBs, and Section 4 solves the mode by improving the GA. The effectiveness of the IIGA is demonstrated by applying the algorithm to a real-world scenario in Section 5. Section 6 presents the conclusions of this research.

Formulation for the Stops and Timetables of CBs
CBs have fixed stops, lines, vehicles, timetable and prices and have the advantages of serving one person, one direct stop and lane flexibility. Because passenger travel demand is diverse, it is difficult to meet different passenger travel demands simultaneously. Hence, many factors must be considered in planning CB lines. There are three modes of CB stops: (1) One bus on one line This model can meet passenger demand when passenger travel demand is small, as shown in Fig 1. (2) Multiple buses on one line The departure times of CB should be broader when the passenger travel demand is greater and the passenger travel time is more scattered. This model is shown in Fig 2. (3) Multiple buses on multiple lines  Additional lines should be opened when the CB network in a city is more mature, passenger demand is greater, and boarding and alighting stops are distributed in multiple zones. This model is shown in Fig 3. CB lines can be represented by the following mathematical models. This paper assumes that a CB line consists of λ 1 boarding stops and λ 2 alighting stops. Boarding stops are represented by s 1 ; s 2 ; . . . ; s l 1 , and alighting stops are represented by r 1 ; r 2 ; . . . ; r l 2 . The vector lx represents all stops passed by each bus: lx ¼ fs 1 ; s 2 ; . . . s l 2 ; r 1 ; r 2 ; . . . r l 2 g. The lines of customized buses between O and D are represented by the matrix LX, where m is the number of vehicles.
. . . If two boarding (alighting) stops are the same for the vehicle-namely,s h j ¼ s h jþ1 (r h j ¼ r h jþ1 )the vehicle parks at this stop once. Accordingly, the number of boarding (alighting) stops is ultimately less than λ 1 (λ 2 ). If all stops and ordering are the same for two vehicles, the two vehicles are in one line; only the departure times are different, as shown in Fig 2. If the solution is m = 1, only one vehicle participates in the operation, and LX is a vector this time, as shown in Fig 1. According to the mathematical model of CB lines, the mathematical model of the relative stop timetable can also be represented by a matrix. The arrival time of all stops in the CB line can be expressed as The arrival time of all stops can also be represented by the following matrix for these vehicles: Among them, vehicle h between stops j+1 and j satisfies the following formulas: where v avg is the CB travel speed and T(x) is the CB dwell time at stop x. All stop dwell times T(x) of vehicle h are related to the number of boarding and alighting stops as follows:

Model for the Stop Planning and Timetables of CBs
This paper builds a model for the stop planning and timetables of CBs considering the passenger, operator, and societal perspectives based on passenger travel space and time demand.

Objective function
CB operation must consider the passenger, operator, and societal perspectives. Passengers seek to minimize their commuting time and maximize their scheduling time, operators seek to maximum operational income while minimizing operational costs, and societal benefits should be considered in selecting the number of operating CBs.
The objective function is calculated according to the three aspects of passengers, operators, and societal benefits. (1) Whether vehicle h meets the travel demand of passenger i There are three main standards used to determine whether vehicle h meets the travel demand of passenger i: (1) the boarding and alighting stops in the line for passenger i; (2) whether the actual vehicle arrival time is earlier than the passenger's expected time of arrival; and (3) passengers who screen stops will prefer to choose vehicle h when some vehicles pass a stop. The condition for meeting the first standard is x s i ; y r i 2 lx h . To meet the second standard, one must first determine whether the passenger's alighting stop is one of the stops. If r h j À y r i ¼ 0 ðj ¼ 1; 2; Á Á Á ; l 2 Þ, passenger i alights at stop j of vehicle h. If t sd i ! t rh j -namely, the passenger's expected latest time of arrival is not earlier than the passenger alighting time at stop j of vehicle h-passenger i can take vehicle h; otherwise, passenger i cannot take vehicle h. If t sd i ) t rh j -namely, the difference between t sd i and t rh j is greater than critical value T-the passenger does not take vehicle h because the wasted time is excessively high for the passenger. To meet the last standard, passengers are considered to give preference to the vehicle whose actual starting time is latest. If two or more vehicles can meet the passenger travel demand, passengers choose the vehicle for which t rh j is maximum. Whether passenger i takes the CB h is described as follows: indicates whether the actual arrival time vehicle h is earlier than the expected time of The time cost of passengers is divided into three parts. The time cost of passengers is represented by C t . First, this paper assumes that the passenger waiting time is zero. Second, the time wasted when a passenger's actual arrival time is earlier than the passenger's expected arrival time is represented by C d (i). Finally, travel cost is represented by C l (i). The travel time cost of passenger i in travel is calculated as follows: The wasted time is calculated as follows: The traveling cost is calculated as follows: In this paper, the time value is calculated by the income approach, which reflects that the wasted time of passengers during travel leads to reduced income because the passengers are unable to work during that time. The calculation method is as follows: where NC(i) is the personal annual income of passenger i and NT(i) is the annual working time of passenger i.
The income approach reflects the difference in time valuation between different passengers. For the income inequality of passengers in the case of the same price, travel time savings is more important for a higher passenger revenue. Furthermore, the travel time value is proportional to passenger income as follows: where α g (i) is the loss factor of the travel time of passenger i. For any type of traffic, the CB can arrive at the destination in advance. The wasted time value when the CB arrives at the destination in advance is proportional to passenger income as follows: where α f (i) is the loss factor of wasted time for passenger i when the CB arrives at the destination ahead of schedule.
For the loss factor of time, it is higher when passenger i considers that the period is no value. Thus, α g (i) is typically greater than α f (i).
The passenger's total time cost is described as follows: (3) Passenger travel value Commuter travel is aimed at going to work, creating societal value, and obtaining a certain salary. Thus, passenger travel also has a certain value for passengers. Passenger travel value is represented by R x .
This paper assumes that passenger travel value and travel time are related to wages as follows: where α x (i) is the proportion of passenger travel value and the wages per unit time and Δt(i) is the travel time of passenger i. CB passenger travel value is described as follows: Thus, all CB passenger travel values are described as follows:

Operators' perspective.
From the operators' perspective, this paper considers that operating a CB consists mainly of operational income and operational cost.
(1) Operational income The operational income of a CB is derived from the fare income of the CB (other income can be neglected). The CB fare is based on the operating distance, which is an integer. The fare income is calculated as follows: where R s is societal benefits, ω p is the fare per unit travel in kilometers and P g is the starting fare for the CB.
(2) Operational cost The operational cost of a CB refers to all expenses in the form of money associated with the consumption of passenger travel, which consists of driver wages and welfare, vehicle fuel and depreciation fees in this paper. The operational cost is calculated as follows: where C y is the operational cost, C pd is the driver cost, C yr is the fuel cost, and C cz is the depreciation cost.
The driver cost consists of fixed wages and variable wages. Variable wages are related to travel kilometers. The driver cost is calculated as follows: where C pdg is the driver's fixed cost, t h cc is the time required to drive vehicle h from the parking lot to the first stop, t h cd is the time required to drive vehicle h from the end of the transport mission back to the parking lot, and ω d is the unit time cost of driving the vehicle.
Fuel cost is proportional to driving distance as follows: where ω r is the fuel cost per unit distance, L h cc is the driving distance of vehicle h from the parking lot to the first stop, and L h cd is the driving distance of vehicle h from the end of the transport mission back to the parking lot.
The value of the vehicle will gradually decrease because of wear and tear. This decrease in value is the depreciation cost, which is considered the average vehicle value of fixed assets according to the service life of the vehicle in this paper. The depreciation cost for each operation of the vehicle is calculated as follows: where C h g is the purchase cost of vehicle h, ε h is the service life of vehicle h, and n h c c is the number of times that vehicle h is used in one day.

Societal perspective.
This paper considers that the societal benefits mainly include road congestion and pollutant emissions. Because the CB service object contains private commuters, if the passengers do not take the CB, then the commuter who drives a private car will increase road congestion and increase pollution gas emissions.
(1) Reduction in traffic congestion cost The congestion cost is proportional to the road area and the driving distance as follows: where ω yj is the congestion cost per unit area, S zd is the occupied road area of vehicles, and L zx is the vehicle driving distance. Thus, for a certain number of passengers, choosing to take a CB will reduce congestion costs relative to choosing private cars as follows: where S car zd is the occupied road area of cars and S bus zd is the occupied road area of a CB. (2) Reduction in environmental pollution cost For a certain number of passengers, choosing to take a CB will reduce the environmental pollutant cost relative to choosing private cars as follows: where o car wb is the external unit cost of the discharged pollutant per kilometer traveled by car and o bus wb is the external unit cost of the discharged pollutant per kilometer traveled by CB. In summary, societal benefits include the reduction in traffic congestion and environmental pollution costs as follows: where R w is the operational income.

Objective function, constraints and assumptions
CB stop design must satisfy constraints at the service level, such as the number of stops, as well as constraints on the lines and stop mileage and the requirements for the vehicle load factor. In addition, the stop design should ensure that passengers take the same CB when they are at the same starting stop. Based on the foregoing considerations, the objective function of the model is as follows: Constraint (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22): The CB will operate when the number of passengers is at least λ 3 on a line, and the rated passenger capacity of the CB is λ 4 .
Constraint : The upper bound λ 6 indicates that the CB line should not be excessively long, which avoids the problem of long travel time. The lower bound λ 5 indicates that the CB line should not be extremely short. When the travel path is sufficiently short, the original public transport system can also meet passenger demands, and the high CB fare reduces the competitiveness of this mode of transportation.

Solution Method
The genetic algorithm (GA) is based on the concept of natural selection. A GA searches for the optimal solution by simulating the natural process of evolution and is one of the most commonly used artificial intelligence algorithms. GAs have been successfully applied in many fields due to their simplicity, robustness and adaptability to parallel distributed processing. However, Belew R K [29] believed GAs still have certain shortcomings, such as premature convergence and lack of local search ability. The immune genetic algorithm (IGA) is a new type of combined optimization method that incorporates the characteristics of the biological immune system based on a GA. An IGA has high antigen recognition, self-regulation and other functions based on the robustness of the GA, which performs well in terms of search speed, global search capability and local search capability. Jiao L and Wang L [30] verified that the IGA is feasible and effective and is conducive to alleviating the degeneration phenomenon in the original GA based on examples of the traveling salesman problem (TSP). The variables are too numerous and difficult to obtain for the model proposed in Section 3; thus, this paper solves the model for the stop planning and timetables of CBs using a heuristic algorithm with an IGA, which searches the feasible solution space rapidly and exhibits an excellent global searching ability. Because the IGA is a single-population design, its global search ability is poor, and it easily falls into a local optimal solution. The fixed crossover operator and mutation operator can prevent the algorithm from being adjusted according to the actual situation of the group and from reflecting the different requirements of the population's antibodies in different evolutionary stages, thus resulting in a poor optimized effect. This paper proposes an IIGA according to the characteristics of the formulated model. The main improvements are as follows: 1. Multiple-population design. The introduction of multiple-population design can enrich the evolutionary direction of different populations and introduce the superior antibody of each subpopulation using the migration operator, which can effectively prevent the premature convergence phenomenon.
2. Memory library design. The introduction of memory library design can ensure that better antibody can be inherited by the next generation, preventing the antibody degradation phenomenon.
3. Designing adaptive crossover probability and mutation probability. The evolution of the crossover probability and that of the mutation probability are determined according to the diversity and excellence of the previous population generation, which prevents the premature convergence phenomenon.
4. Fitness calculation of the gene segment. Because the model solution is composed of a number of stops and the first stop time, it should be as far as possible to reduce the crossover and mutation operation for a vehicle of excellent performance; thus, the convergence to the optimal solution will be accelerated.

Designing the IIGA
(1) Antibody coding The model decision variables are LX and all vehicles departure times t sh 1 . Because the number of decision variables is greater, a binary code cannot be adopted; thus, this paper adopts real code, and the length of the antibody is the number of decision variables. For each vehicle, all stops lx h and starting time t sh 1 are taken as a unit to form an antibody according to the order of the vehicles, as shown in Fig 4. (2) Initial antibody population generation The initial antibody population size M and the length of the antibody in population (λ 1 + λ 2 + 1)m, m 2 [1, m max ] are set. Each antibody randomly selects a different length structure-namely, it randomly selects the value of m; then, the population contains a number of different vehicles.
The specific methods for determining the internal variables of the antibody are as follows: First, all passengers with a given boarding stop demand and alighting stop demand are placed  is randomly generated from a range of values. The initial population attempts to set the stops nearby because the stop number rule is designed according to the ascending sequential order of the X-and Y-coordinates of the stops.
In this paper, multiple-population design divides the antibody population into three subgroup antibody populations-namely, an advanced antibody population, a behindhand antibody population and an intermediate antibody population. Under the condition that the crossover probability and mutation probability have been obtained, the different subgroup antibody populations must be multiplied by different correction values, which increases the crossover probability and mutation probability and reduces the crossover probability and mutation probability of the behindhand antibody population. The correction values of the advanced, behindhand, and intermediate antibody populations are γ 5 , γ 6 and γ 7 , respectively. The process involves the independent calculation of each antibody affinity and immune genetic operation. The subgroup antibody population will exchange several optimal antibodies every few evolutionary algebraic operations; in this manner, a more excellent antibody is introduced at the same time and enriches the diversity of the antibody population, which prevents immature and early convergence. (

3) Affinity between antibodies and antigens
The affinity between antibodies and antigens can be obtained by changing the objective function and constraints: A more optimized antibody is associated with a greater objective function value of the corresponding objective function and greater affinity. The calculated formula is as follows: where pc is a penalty constant, η n (v) = {0,1} indicates whether antibody V violates the nth constraint, and I is the number of constraints. Thus, this paper proposes the concept of gene segments. Because the solution of the model is composed of the stops of all vehicles and the first stop time, the solution should minimize the crossover and mutation operations for a vehicle of the excellent performance and thereby speed up the convergence of the optimal solution. A vehicle in an antibody can be viewed as a gene segment. The fitness calculation of the gene segment is described as follows: where v h is the gene values, which are fixed from (λ 1 +λ 2 +1)(h−1)+1 to (λ 1 +λ 2 +1)h and are 0 otherwise. (4) Transport operator In this paper, the design concept of the transport operator is as follows: The algorithm defines transport interval γ 8 and transport quantity γ 9 , where γ 8 is the frequency of transport behavior between subgroup antibody populations and γ 9 is the number of each transport antibody. When the IGA evolves to γ 8 over multiple generations, the best γ 9 antibodies for each subgroup antibody population are transferred to other antibody groups.
(5) Memory library design The memory library design can prevent the destruction of better antibodies by a crossover operator or mutation operator and retain better antibodies in the subsequent generation to accelerate the algorithm convergence. The memory library design is performed as follows: 1. Setting the initial memory library. The memory library size is J; the previous J antibodies are placed into the initial memory library of each subspecies population in descending order of the subpopulation antibody affinity.
2. Updating the memory library. After genetic manipulation, the affinity of the subpopulation antibody is sorted. If the affinity of the newly generated antibody is higher than that of antibodies in the original memory library, these new antibodies will be placed into the memory library; then, the same number of antibodies for which the affinity is the worst in the memory bank will be replaced. Otherwise, the memory library remains the same. The same antibody is not selected to ensure the diversity of the memory library.
3. The antibodies of the memory library in the next generation. After each memory library is updated, the antibodies in the memory library are directly inherited by the next generation.

(6) Antibody of promotion and inhibition
In the iterative process of the IGA, the antibody concentration is increased to a certain value that is inhibited, improving the production and selection probability of the antibody for which the concentration is lower. The affinity between antibodies is calculated to determine the concentration of the antibody and thus determine whether the antibody is inhibited or promoted.
1) Calculation of the affinity between antibodies Lee W, Lee S, Lee B H, et al. [31] proposed that the affinity between antibodies is calculated based on Shannon information entropy theory. First, the degree of similarity between the antibodies is calculated. Each antibody corresponds to a solution X i . The affinity degree A ij between antibody X i and antibody X j indicates the degree of similarity between the antibodies. The Shannon information entropy is shown in Fig 5. The information entropy of the gene at location j is as follows: where p ij is the appearance probability of the jth position in the antibody for allelic gene i. The diversity of the average information entropy is as follows: According to the definition of entropy, the affinity between the antibodies is as follows: According to the affinity between antibodies, the antibody concentration can be calculated as follows: where and T acl is the set threshold for antibody affinity (0.9 T acl 1).

2) Antibody viability
Viability is the ability to survive to the next iteration of the antibody and is calculated as follows: where ρ v is the antibody concentration between the antibody and antibody v. (7) Selection operator In this paper, the operation selection method is a roulette wheel. The probability of antibodies being selected is based on the antibody viability e v . The roulette wheel method operates as follows: A random number r is randomly generated between 0 and 1 and compared with the cumulative probability P v of each antibody viability. If P v−1 <r P v , antibody v will be inherited by the next generation.
The cumulative probability of antibody v is as follows: According to the roulette wheel selection, the probability p v of the antibody is selected as follows: (8) Adaptive crossover operator based on the gene segment In this paper, an improved single-point crossover operator can increase the crossover probability of the diversity of antibodies and reduce the crossover operation of the gene segments with excellent performance. The steps of this operator are shown in Fig 6. Step 1. Select two antibodies to cross: Determine whether the two antibodies cross based on the crossover probability. The probability of antibody crossover is calculated by an adaptive method and is determined by the antibody population diversity of the current generation. According to the expression of the logistic curve equation, the crossover probability is defined as follows: Step 2. Select the gene segment in the crossover location: When selecting the shorter antibody from two crossover antibodies, the probability that each segment is selected depends on the fitness of the gene segment. The probability of segment h is selected as follows: Step 3. Randomly select the crossover position in the gene segment.
Step 4. Determine the antibody sequence after exchanging each crossover point. (9) Adaptive mutation operator based on the gene segment The mutation operator design resembles the crossover operator, reducing the crossover probability of the antibody diversity and reducing the crossover operation of the excellent gene segment. The corresponding steps are shown in Fig 7. Step 1. Select two antibodies to implement the mutation: The antibody mutation probability, which is calculated by an adaptive approach and is determined by the antibody population current diversity, is used to determine whether mutation is needed.
The mutation probability is defined as follows: Step 2. Select a gene segment in the mutation location: The probability that the segment is selected depends on the fitness of the gene segment. The probability that segment h is selected is as follows: Step 3. Randomly select the mutation position in the gene segment.
Step 4. Change the gene value of the mutation point.
(10) Immune operator The immune operator is divided into three main parts: the extraction vaccine, vaccination and immune selection. In this paper, the extraction vaccine is based on the passenger travel demand, with γ 1 typical boarding stops, γ 2 typical alighting stops, and γ 3 dense departure times. The method for selecting typical boarding stops and typical alighting stops is as follows: γ 1 typical boarding stops and γ 2 typical alighting stops for which the demand is maximum are selected after summarizing and analyzing the passenger travel demand of boarding stops and alighting stops. The dense departure times are selected as follows: The travel time is estimated according to the distance between the boarding zone and alighting zone, then a detour coefficient α rx is set and detour travel time is α rx times as long as the distance of the travel time.γ 3 expected arrival times are then chosen, which causes all passengers to require as much as possible in the range of γ 4 minutes before and after expected arrival times. γ 3 dense departure times are calculated based on the detour travel time.
The vaccination method is performed as follows: The gene locus of the required vaccination antibody is divided into three parts. The first part is the boarding stop coding region-namely, the range of [(λ 1 + λ 2 + 1)h + 1, (λ 1 + λ 2 + 1)h + λ 1 ] (h = 0,1,..,m max − 1) for the antibody. The second part is the alighting stop coding region-namely, the range of [(λ 1 + λ 2 + 1)h + λ 1 + 1, (λ 1 + λ 2 + 1)h + λ 1 +λ 2 ] (h = 0,1,..,m max − 1) for the antibody. The third part is the departure time coding region-namely, {(λ 1 + λ 2 + 1)h + λ 1 + λ 2 + 1} (h = 0,1,..,m max − 1) for the antibody. One mutation point is randomly selected from these three regions, and the gene values of these genes are changed. A value is selected accordingly and randomly from the values of typical boarding stop numbers, typical alighting stop numbers and dense departure times in the vaccine, and these values must be replaced. The immunity process is shown in Fig 8. The probability of immunity selection is set to a certain value p b ; thus, antibodies that must implement the immunity selection operation based on p b are randomly selected. (11) Termination criterion of the IIGA The termination criteria used in this paper is a certain number of iterations that determines a sufficiently large positive number T, and the total number of iterations for the algorithm does not exceed T. Due to the limited computation time and capacity, the number of iterations is not infinite. The advantages of this criterion are that the computation time can be easily controlled.

Flow of the IIGA
The steps for designing the IIGA in this paper are as follows. The flowchart is shown in Fig 9. Step 1. Antigen recognition. The problem of stop planning and timetable creation for CBs must be abstracted to the form of an antigen in accord with the behavior of the immune system, i.e., the problem is solved.
Step 2. Generate the initial antibody population. Three initial antibody populations are generated based on the principle of generating the initial antibody population and randomly generating the initial time in a reasonable range.
Step 3. Calculate the affinity. The affinity between the antibody and antigen in each subpopulation is calculated.
Step 4. Determine whether the antibody is consistent with the transport condition. If the antibody is consistent with the transport condition, the transport operator must implement the transport operation according to the algorithm. Otherwise, advance to Step 5.
Step 5. The antibody population with better affinity is stored in the memory library of the subpopulations, and the antibody genes in the memory library are directly inherited by the next generation population.
Step 6. Promote and inhibit the antibody. The antibody concentration is calculated by calculating the affinity between antibodies; the antibody viability is then calculated.
Step 7. According to the design of the algorithm, the antibody performs the selection, crossover, mutation and immune operations, and the antibody is updated.
Step 8. Determine whether the algorithm is consistent with the termination condition. If the algorithm is consistent with the termination condition, the optimal solution is the output. Otherwise, return to Step 2.

Case Study
Passenger travel demand is concentrated for Liyuan and Guomao. Liyuan covers an area of approximately 18 km 2 , and Guomao covers an area of approximately 3 km 2 . After demand acquisition, the demand is 568 in the morning peak period from Liyuan to Guomao; the travel densities in Liyuan and Guomao are 47 and 189 people/km 2 , respectively; and the travel demand is dense. After screening, these requirements are distributed to 57 residential districts in the Liyuan community and 57 office buildings in Guomao, which are scattered in the two According to the coding mode implemented in this paper, t sd i 2 f135; 150; 165; 180; 195g. Every passenger boarding time t s is 5 s, the alighting time t x is 5 s, and the necessary dwell time is 30 s. Passengers' personal annual income NC(i) is 120,000 yuan, the annual number of working days is 250, the annual working time NT(i) is an average of 120,000 min according to 8 working hours per day, the loss factor of travel time α g (i) is 100%, the loss factor of arriving Stop Planning and Timetables of Customized Buses at the destination ahead of schedule α f (i) is 50%, and the proportion of passenger travel value and wages per unit time α x (i) is 80%.
According to Beijing CB fees, the travel fee consists of a seat fee and credit card fee for each take. The standard seat fee within 20 km (inclusive) is 8 yuan and increases by 3 yuan with each additional 5 km (inclusive). The average CB speed is 40 km, the driver fixed cost C pdg is 50 yuan, the unit time cost of driving a vehicle is ω d 1 yuan/min, the average passenger weight is 60 kg, the vehicle weight is 11,500 kg, the fuel cost of unit distance ω r is 2.5×10 −4 yuan/ (km•kg), the vehicle purchase fee C h g is 700,000 yuan, the vehicle service life ε h is 8 years, the number of vehicle uses per day n h c c is 2, and the congestion cost of unit area ω yj is 0.02. Passenger boarding stop and alighting stop demand data are calculated by the K-means algorithm in MATLAB according to the distance between each boarding (alighting) stop. The clustering results meet the demand for a distance within 500 m because passengers can walk to alternate stops within this distance. After K-means clustering, 37 and 14 alternative stops are identified for Liyuan and Guomao, respectively. The location of the center of mass at which the clustering results are generated is an alternative stop. For the simple calculation, the alternative site number selects the original stop number that is closest to the center of mass, which is no longer renumbered.

Solution to the example using the IIGA
The parameters of the model and algorithm are set as shown in Table 1. The range of the crossover probability is [0.48, 0.9), and the range of the mutation probability is (0, 0.1].
The optimal solution of the model is obtained by the IIGA as follows:  Table 2.
Similarly, the boarding stops, alighting stops and relevant information can be drawn for the second, third and fourth CBs. Passenger information and vehicle information are collected for four lines, as shown in Table 3.

Analysis and comparison of results
The results of the IIGA and IGA are compared in terms of the number of stops, number of passengers, vehicle driving distance, travel time, average passenger travel time, average passenger arrival time ahead of schedule, and total line revenue in Table 4.  The results of the comparison in Table 4 illustrate that for the IIGA relative to the IGA, the numbers of boarding stops and alighting stops increased by 7 and 1, respectively; the number of passengers increased by 31; the vehicle driving distance increased by 1.2 km; the travel time increased by 16 min; the average passenger travel time increased by 1.8 min; the average passenger arrival time ahead of schedule decreased by 0.9 min; and the total line revenue increased by 220.3 yuan. The increase in the numbers of boarding stops and alighting stops indicates that a CB line covers a larger area; therefore, the distance, travel time and average travel time of the vehicle increase. The increase in the number of passengers causes the total line revenue to increase. Moreover, the average passenger arrival time ahead of schedule meets the passenger demand. In terms of the speed of the algorithm, the IGA required 532 s to complete, whereas the IIGA required 268 s. An iteration diagram of the total line revenue and vehicle travel distance is shown in Fig 12. The IIGA also significantly outperforms the IGA in terms of the CB total line revenue. The IIGA shows advantages with respect to the number of stops, number of passengers, average passenger arrival time ahead of schedule and total line revenue, but the vehicle driving distance, travel time and average passenger travel time are inferior. In addition, the IIGA is significantly faster than the IGA, with a 49% faster running speed. Thus, the proposed model has strong potential for use in solving CB problems.

Conclusions
Based on an analysis of domestic and international research and the current development of CBs in China and considering passenger travel data, this paper studies the problems associated with the operation of CBs, such as stop selection, line planning and timetable setting, and establishes a model for the stop planning and timetable of CBs. The improved immune genetic algorithm (IIGA) is used to solve the model, which includes the following features: 1) multiple population design and transport operator design, 2) memory library design, 3) mutation probability design and crossover probability design, and 4) the fitness calculation of the gene segment. Finally, a real-world example set in Beijing is calculated, and the model and solution results are verified and analyzed. The results illustrate that the IIGA solves the model and is superior to the basic genetic algorithm in terms of the number of passengers, travel time, average passenger travel time, average passenger arrival time ahead of schedule and total line revenue. This study covers the key issues involving the operational system of CBs, combines theoretical research and empirical analysis, and provides a theoretical foundation for the planning and operation of CBs.
Future research will focus on the following aspects: 1. This paper does not consider the time delay caused by road congestion, and the arrival time of each stop is ideal. Future studies can combine road traffic flow data and arrival prediction theory to determine the exact arrival time of each stop, thus improving the punctuality rate and service level of CBs.
2. Because CBs are still in the initial stage of development for most cities and are simply attached to existing transit fleets, in this paper, the vehicle demands between two traffic zones optimize computation without considering fleet size restrictions. Therefore, in the future, the problem should be linked to fleet size, and relationships among travel demands, number of lines, boarding and alighting stops, departure times and fleet size should be considered.
Supporting Information S1 Dataset. The passenger flow survey data. (XLSX)